Repository: kaspersgit/ml_2_sql (language: python)

Issue:
Bump pillow from 9.1.1 to 9.3.0
Bumps [pillow](https://github.com/python-pillow/Pillow) from 9.1.1 to 9.3.0.
- [Release notes](https://github.com/python-pillow/Pillow/releases)
- [Changelog](https://github.com/python-pillow/Pillow/blob/main/CHANGES.rst)
- [Commits](https://github.com/python-pillow/Pillow/compare/9.1.1...9.3.0)

---
updated-dependencies:
- dependency-name: pillow dependency-type: direct:production ...

Signed-off-by: dependabot[bot] <support@github.com>


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 0707758..2999d25 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,10 +10,15 @@ input/configuration/*.json
 trained_models/*
 !trained_models/example_titanic/
 
+manual_testing.py
+
 # Editors
 .vscode/
 .idea/
 
+# Ruff
+.ruff*
+
 # Vagrant
 .vagrant/
 
diff --git a/Readme.md b/Readme.md
index 4b06d0a..0a0579a 100644
--- a/Readme.md
+++ b/Readme.md
@@ -13,20 +13,23 @@
 5. [Troubleshooting](#troubleshooting)
 
 
-# What is it?
-## Background
-Due to SQL being the main language being used for data manipulation and thus has extensive support in terms of compute and scheduling,
-why not perform inference with a machine learning model written in SQL code? The big limitation here is SQL itself, that's why we attempt to use
-machine learning models which have a simple structure it is writable in SQL. One additional benefit of this is that the model is interpretable,
-if you can write down the model in a basic logical language (SQL) you should be able to understand it (with limitation ofcourse).
-
+# In short
 This project tries to make the process simple enough for any SQL user to train a model, check the performance and deploy that model in SQL.
 
+## Philosophy:
+- Automated and easy to use 
+- Model in SQL (avoid Data Science debt for a complex pipeline)
+- Explainable Boosting Machine on par with other boosted methods
+
+## Background
+An automated machine learning tool which trains, graphs performance and saves the model in SQL. Using interpretable ML models (from interpretML) to train models which are explainable and transparent in how they come to their prediction. SQL infrastructure is the only requirement to put a model into production.
+This tool can be used by anybody, but is aimed for people who want to easily train a model and want to use it in their (company) SQL system. 
+
 ## Current state
 - Only EBM is implemented (decision tree, logistic regression and rule set not yet)
 - Automated model training is working for binary classification and regression
 - SQL creation of model is working fully for binary clasification
-- SQL for regression and the whole process for multiclass classification is wip
+- SQL for regression and the whole process for multiclass classification is work in progress
 - Up/down sampling not fully implemented yet 
 
 </br>
@@ -56,21 +59,33 @@ This project tries to make the process simple enough for any SQL user to train a
    
     Mac/Linux:
    `source .ml2sql/bin/activate` or `deactivate`
+4. For usage of this tool the virtual environment can be deactivated.
 
 ## Try it out demo
-1. Make sure virtual environment is active (step 3 of Pre requisites)
-2. In the terminal in the root of this folder run: `python run.py`
-3. Follow the instructions on screen by selecting the demo data and (similarly named) config file
-4. Check the output in the newly created folder
+1. In the terminal in the root of this folder run: 
+  - `python3 run.py` (Mac/Linux)
+  - `python run.py` (Windows)
+2. Follow the instructions on screen by selecting the demo data and (similarly named) config file
+3. Check the output in the newly created folder
 
 ## Try it out using own data
-1. Make sure virtual environment is active (step 3 of Pre requisites)
-2. Save csv file containing target and all features in the `input/data/` folder (more on input data at [input data](#data))
-3. Save a settings json file in the `input/configuration/` (explained below at [configuration json](#configuration-json))
-4. In the terminal run: `bash run.sh`
-5. Follow the instruction on screen
-6. The output will be saved in the folder `trained_models/<current_date>_<your_model_name>/`
-7. The `.sql` file will contain a SQL Case When statement imitating the decision tree/EBM  
+1. Save csv file containing target and all features in the `input/data/` folder (more on input data at [input data](#data))
+2. Save a settings json file in the `input/configuration/` (explained below at [configuration json](#configuration-json))
+3. In the terminal in the root of this folder run: 
+  - `python3 run.py` (Mac/Linux)
+  - `python run.py` (Windows)
+4. Follow the instruction on screen
+5. The output will be saved in the folder `trained_models/<current_date>_<your_model_name>/`
+6. The `.sql` file in the `model` folder will contain a SQL written model  
+
+## Testing already trained model on a new dataset
+1. Make sure new dataset is exactly the same format as the dataset the model was trained on
+2. Save dataset in the `input/data/` folder (more on input data at [input data](#data))
+3. In the terminal in the root of this folder run: 
+  - `python3 test_model.py` (Mac/Linux)
+  - `python test_model.py` (Windows)
+4. Follow the instructions on screen
+5. The output will be saved in the folder `trained_models/<selected_model>/tested_datasets/<selected_dataset>/`
 
 </br>
 
@@ -98,8 +113,8 @@ Dictionary of parameters that can be used with model of choice (optional). Check
 - any other value, no calibration applied
 
 `sql_split` options:
-- `false`, outputs the SQL model as one column by adding all separate scores up directly
-- `true`, outputs the SQL model as one column for each feature and a total score columns afterwards. This might be needed to avoid some memory related (stackoverflow) error.
+- `false`, outputs the SQL model as one SELECT statement, using column aliases within the same select statement
+- `true`, outputs the SQL model as several CTEs, this can be used if column aliases can't be referenced within the same SELECT statement
 
 `file_type` options (optional):
 - `png`, output of features importance graphs will be static .png (smaller file).
@@ -116,7 +131,6 @@ Dictionary of parameters that can be used with model of choice (optional). Check
 `time_sensitive_column` options (optional):
 - Name of date column
   - used when `cv_type = timeseriesplit`  
-  - used when out-of-time dataset is created (not implemented yet)
 
 `upsampling` options (optional, should not be used without calibration):
 - `true`, applying the SMOTE(NC) algorithm on the minority class to balance the data
@@ -132,29 +146,32 @@ Name of target column (required)
 ## Notes
 - Any NULL values should be imputed before using this script
 - Data imbalance treatments (e.g. oversampling + model calibration) not fully implemented
-- Resampling (almost) always makes the trained model ill calibrated
 - Multiclass and regression are experimental
 
 ## TODO list
-- Add decision tree
-- Add logistic regression
-- Add Skope rules
-- Add calibration (platt scaling/isotonic regression)
-- Add changelog and versioning
-- Implement null handling (there is an implementation mentioned [here](https://github.com/interpretml/interpret/issues/18))
-- Make multi class classification EBM work fully
-- Make regression EBM work fully
-- Removing outliers by using quantiles (e.g. only keeping 1 - 99 % quantiles)
-- Spatial Cross-validation discovery
-- Extend logging granularity (add model parameters)
-- Use menu function bash for model type choosing
-- Add target single unique value check
-- Replace/improve `classification_report` and `confusion_matrix` due to dependance on threshold
-- Add MCC, cohen kappa and other metrics plotted with threshold
-- Add SQL translation for decision rule
-- Add random seed to config file
-- Make a proper testing procedure
-- Makefile (or bash script) to setup virt env and install packages (instead of running 3 setup lines)
-- Improve ReadMe
-- Add more example files (binary/regression/multiclass)
-- Add tool to make configuration json
+- EBM
+  - SQL interaction terms first term with least values
+  - SQL interaction terms group by to merge similar score bounds
+  - Make multi class classification EBM work fully
+  - Make regression EBM work fully
+
+- Checks and config
+  - Add check if variables and target are finite 
+  - Add variables being NULL checked
+  - Add random seed to config file
+  - Implement null handling (there is an implementation mentioned [here](https://github.com/interpretml/interpret/issues/18))
+
+- Performance monitoring
+  - Add performance summary for easy and quick comparison (including label count, auc pr & roc, best f1-score, etc)
+  - Add feature over/under fitting plot (https://towardsdatascience.com/which-of-your-features-are-overfitting-c46d0762e769)
+  - Make distribution plot grouped instead of overlaid or stacked (maybe switch to plotly histogram)
+
+- Other 
+  - Add calibration (platt scaling/isotonic regression)
+  - Add changelog and versioning
+  - Extend logging granularity 
+  - Add SQL translation for decision rule
+  - Make a proper testing procedure
+  - Add more example files (binary/regression/multiclass)
+
+# Troubleshooting
\ No newline at end of file
diff --git a/docs/index.md b/docs/index.md
new file mode 100644
index 0000000..000ea34
--- /dev/null
+++ b/docs/index.md
@@ -0,0 +1,17 @@
+# Welcome to MkDocs
+
+For full documentation visit [mkdocs.org](https://www.mkdocs.org).
+
+## Commands
+
+* `mkdocs new [dir-name]` - Create a new project.
+* `mkdocs serve` - Start the live-reloading docs server.
+* `mkdocs build` - Build the documentation site.
+* `mkdocs -h` - Print help message and exit.
+
+## Project layout
+
+    mkdocs.yml    # The configuration file.
+    docs/
+        index.md  # The documentation homepage.
+        ...       # Other markdown pages, images and other files.
diff --git a/input/configuration/example_binary_titanic.json b/input/configuration/example_binary_titanic.json
index e73b173..960f236 100644
--- a/input/configuration/example_binary_titanic.json
+++ b/input/configuration/example_binary_titanic.json
@@ -8,6 +8,7 @@
     "Fare"
   ],
   "model_params": {
+      "interactions": 5
   },
   "post_params": {
     "calibration": "false",
diff --git a/input/configuration/example_multiclass_faults.json b/input/configuration/example_multiclass_faults.json
index e979be3..1dd0b63 100644
--- a/input/configuration/example_multiclass_faults.json
+++ b/input/configuration/example_multiclass_faults.json
@@ -29,6 +29,7 @@
     "SigmoidOfAreas"
   ],
   "model_params": {
+	"max_bins": 5
   },
   "post_params": {
     "calibration": "false",
diff --git a/input/data/example_multiclass_faults.csv b/input/data/example_multiclass_faults.csv
index eb72caa..92cd27d 100644
--- a/input/data/example_multiclass_faults.csv
+++ b/input/data/example_multiclass_faults.csv
@@ -1,4 +1,4 @@
-,X_Minimum,X_Maximum,Y_Minimum,Y_Maximum,Pixels_Areas,X_Perimeter,Y_Perimeter,Sum_of_Luminosity,Minimum_of_Luminosity,Maximum_of_Luminosity,Length_of_Conveyer,TypeOfSteel_A300,TypeOfSteel_A400,Steel_Plate_Thickness,Edges_Index,Empty_Index,Square_Index,Outside_X_Index,Edges_X_Index,Edges_Y_Index,Outside_Global_Index,LogOfAreas,Log_X_Index,Log_Y_Index,Orientation_Index,Luminosity_Index,SigmoidOfAreas,target
+X_Minimum,X_Maximum,Y_Minimum,Y_Maximum,Pixels_Areas,X_Perimeter,Y_Perimeter,Sum_of_Luminosity,Minimum_of_Luminosity,Maximum_of_Luminosity,Length_of_Conveyer,TypeOfSteel_A300,TypeOfSteel_A400,Steel_Plate_Thickness,Edges_Index,Empty_Index,Square_Index,Outside_X_Index,Edges_X_Index,Edges_Y_Index,Outside_Global_Index,LogOfAreas,Log_X_Index,Log_Y_Index,Orientation_Index,Luminosity_Index,SigmoidOfAreas,target
 0,42,50.0,270900,270944,267,17,44,24220,76,108,1687,1,0,80.0,0.0498,0.2415,0.1818,0.0047,0.4706,1.0,1.0,2.4265,0.9031,1.6435,0.8182,-0.2913,0.5822,Pastry
 1,645,651.0,2538079,2538108,108,10,30,11397,84,123,1687,1,0,80.0,0.7647,0.3793,0.2069,0.0036,0.6,0.9667,1.0,2.0334,0.7782,1.4624,0.7931,-0.1756,0.2984,Pastry
 2,829,835.0,1553913,1553931,71,8,19,7972,99,125,1623,1,0,100.0,0.971,0.3426,0.3333,0.0037,0.75,0.9474,1.0,1.8513,0.7782,1.2553,0.6667,-0.1228,0.215,Pastry
diff --git a/requirements.txt b/requirements.txt
index c322586..37868cc 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,9 +2,11 @@ appnope==0.1.3
 backcall==0.2.0
 Brotli==1.0.9
 certifi==2022.12.7
+cffi==1.15.1
 charset-normalizer==2.0.12
 click==8.1.2
 cloudpickle==2.0.0
+colorama==0.4.6
 cycler==0.11.0
 dash==2.3.1
 dash-core-components==2.0.0
@@ -19,14 +21,15 @@ Flask==2.1.1
 Flask-Compress==1.11
 fonttools==4.33.1
 gevent==21.12.0
+ghp-import==2.1.0
 greenlet==1.1.2
 idna==3.3
 imageio==2.17.0
 imbalanced-learn==0.9.0
 imblearn==0.0
 importlib-metadata==4.11.3
-interpret==0.2.7
-interpret-core==0.2.7
+interpret==0.4.2
+interpret-core==0.4.2
 ipykernel==6.13.0
 ipython==7.32.0
 itsdangerous==2.1.2
@@ -37,45 +40,60 @@ jupyter-client==7.2.2
 jupyter-core==4.10.0
 kaleido==0.2.1
 kiwisolver==1.4.2
+lazy_loader==0.2
 lime==0.2.0.1
-llvmlite
+littleutils==0.2.2
+llvmlite==0.40.0
+Markdown==3.3.7
 MarkupSafe==2.1.1
 matplotlib==3.5.1
 matplotlib-inline==0.1.3
+mergedeep==1.3.4
+mkdocs==1.4.2
 multiprocess==0.70.12.2
 nest-asyncio==1.5.5
 networkx==2.6.3
-numba==0.55.1
-numpy==1.21
+numba==0.57.0
+numpy==1.21.0
+outdated==0.2.2
 packaging==21.3
-pandas==1.3.5
+pandas==1.5.3
+pandas-flavor==0.5.0
 parso==0.8.3
 pathos==0.2.8
+patsy==0.5.3
 pexpect==4.8.0
 pickleshare==0.7.5
 Pillow==9.3.0
+pingouin==0.5.3
 plotly==5.7.0
 pox==0.3.0
 ppft==1.6.6.4
 prompt-toolkit==3.0.29
 psutil==5.9.0
 ptyprocess==0.7.0
+pycparser==2.21
 Pygments==2.11.2
 pyparsing==3.0.8
 python-dateutil==2.8.2
 pytz==2022.1
 PyWavelets==1.3.0
+PyYAML==6.0
+pyyaml_env_tag==0.1
 pyzmq==22.3.0
 requests==2.27.1
 SALib==1.4.5
 scikit-image==0.19.2
 scikit-learn==1.0.2
 scipy==1.7.3
+seaborn==0.12.2
 shap==0.40.0
 six==1.16.0
 sklearn==0.0
 skope-rules==1.0.1
 slicer==0.0.7
+statsmodels==0.14.0
+tabulate==0.9.0
 tenacity==8.0.1
 threadpoolctl==3.1.0
 tifffile==2021.11.2
@@ -83,10 +101,13 @@ tornado==6.1
 tqdm==4.64.0
 traitlets==5.1.1
 treeinterpreter==0.2.3
-typing-extensions==4.2.0
+typing_extensions==4.2.0
 urllib3==1.26.9
+watchdog==2.2.1
 wcwidth==0.2.5
 Werkzeug==2.1.1
+xarray==2023.5.0
 zipp==3.8.0
 zope.event==4.5.0
 zope.interface==5.4.0
+pywin32==304; sys_platform == 'win32'
\ No newline at end of file
diff --git a/run.py b/run.py
index a9b5ee2..8f4d75d 100644
--- a/run.py
+++ b/run.py
@@ -21,6 +21,7 @@
 # List files in input/data/ directory
 data_dir = "input/data/"
 files = os.listdir(data_dir)
+files.sort()
 
 print("Files in input/data/:")
 for i, file in enumerate(files, 1):
@@ -29,7 +30,7 @@
 # Ask for CSVPATH
 csv_path = None
 while csv_path is None:
-    csv_file_index = input("\nSelect path to csv file for training the model: ")
+    csv_file_index = input("\nSelect CSV file for training the model: ")
     try:
         csv_file_index = int(csv_file_index) - 1
         csv_path = os.path.join(data_dir, files[csv_file_index])
@@ -42,6 +43,7 @@
 # List files in input/configuration/ directory
 configuration_dir = "input/configuration/"
 files = os.listdir(configuration_dir)
+files.sort()
 
 print(f"\n\nFiles in {configuration_dir}:")
 for i, file in enumerate(files, 1):
@@ -50,7 +52,7 @@
 # Ask for JSONPATH
 json_path = None
 while json_path is None:
-    json_file_index = input("\nSelect path to json file for model configuration: ")
+    json_file_index = input("\nSelect JSON file for training configuration: ")
     try:
         json_file_index = int(json_file_index) - 1
         json_path = os.path.join(configuration_dir, files[json_file_index])
@@ -60,7 +62,7 @@
 print(f"Configuration file {json_path} will be used for modelling")
 
 # Ask for MODEL_TYPE
-model_types = ["Explainable Boosting Machine", "Decision Tree", "Decision Rule"]
+model_types = ["Explainable Boosting Machine", "Decision Tree", "Decision Rule", "Logistic/Linear regression"]
 model_type = None
 while model_type is None:
     print("\nWhat type of model do you want?")
@@ -76,17 +78,28 @@
         print("Invalid option")
 
 print(f"Algorithm chosen for modelling: {model_type}")
+
+# Rename model choice for easier reference in code later on
 model_type = 'ebm' if model_type == 'Explainable Boosting Machine' else model_type
+model_type = 'l_regression' if model_type == 'Logistic/Linear regression' else model_type
 model_type = model_type.lower().replace(" ", "_")
-print(model_type)
 
 # Model name
-model_name = input("\nGive it a name: ")
+unique_name = False
+while not unique_name: 
+    model_name = input("\nGive it a name: ")
+
+    # Current date
+    current_date = datetime.today().strftime('%Y%m%d')
 
-# Current date
-current_date = datetime.today().strftime('%Y%m%d')
+    full_model_name = f"{current_date}_{model_name}"
 
-full_model_name = f"{current_date}_{model_name}"
+    # Check if folder already exists with this name
+    if os.path.isdir(f"trained_models/{full_model_name}"):
+        print('Folder with this name already exists please try another')
+    else:
+        unique_name = True
+        
 
 # Make directory with current data and model name
 try:
@@ -100,9 +113,13 @@
 
 print("Starting script to create model")
 
-# Run python script with given input
+# Run script with given input using python in the venv (so venv does not need to be activated beforehand)
 print("\nRunning command:")
-command = f"python main.py --name trained_models/{full_model_name} --data_path {csv_path} --configuration {json_path} --model {model_type}"
+if sys.platform == 'win32':
+    command = f".ml2sql\Scripts\python.exe scripts/main.py --name trained_models/{full_model_name} --data_path {csv_path} --configuration {json_path} --model {model_type}"
+else:
+    command = f".ml2sql/bin/python scripts/main.py --name trained_models/{full_model_name} --data_path {csv_path} --configuration {json_path} --model {model_type}"
+
 print(command)
 
 os.system(command)
diff --git a/create_config.py b/scripts/create_config.py
similarity index 53%
rename from create_config.py
rename to scripts/create_config.py
index 75669c5..88eed5a 100644
--- a/create_config.py
+++ b/scripts/create_config.py
@@ -1,8 +1,42 @@
 import csv
 import json
+import os
+
+def get_input(options, message):
+    print(f"\n{message}")
+    choice = None
+    while choice is None:
+        for i, option in enumerate(options):
+            print(f"{i + 1}. {option}")
+        response = input("Enter the number of your choice: ")
+        try:
+            choice = options[int(response) - 1]
+        except (ValueError, IndexError):
+            print("Invalid option, try again.")
+    return choice
+
+# List files in input/data/ directory
+data_dir = "input/data/"
+files = os.listdir(data_dir)
+
+print("Files in input/data/:")
+for i, file in enumerate(files, 1):
+    print(f"{i}. {file}")
+
+# Ask for CSVPATH
+csv_path = None
+while csv_path is None:
+    csv_file_index = input("\nSelect path to csv file to create config for: ")
+    try:
+        csv_file_index = int(csv_file_index) - 1
+        csv_path = os.path.join(data_dir, files[csv_file_index])
+    except (ValueError, IndexError):
+        print("Invalid option, try again.")
+
+print(f"CSV file {csv_path} will be used to create a config file for")
 
 # Load the example CSV file
-with open("input/data/example_binary_titanic.csv") as f:
+with open(csv_path) as f:
     reader = csv.reader(f)
     headers = next(reader)
 
@@ -60,6 +94,9 @@
 # Create the features list from the selected features indices
 features = [features_dict[i] for i in selected_features]
 
+# Non used columns
+[features_dict.pop(i) for i in selected_features]
+
 # Fill in the JSON template with the selected features
 input_params = {
   "features": features,
@@ -99,15 +136,46 @@
     else:
         if len(input_params[key]) > 1:
             for key2 in input_params[key]:
-                input_value = input(f"Enter the value for {key2}: ")
+                if key2 == 'calibration':
+                    # input_value = get_input(["true","false"], "Calibration (true/false): ")
+                    input_value = "false"
+                elif key2 == 'sql_split':
+                    input_value = get_input(["true","false"], "SQL split (true/false): ")
+                elif key2 == 'file_type':
+                    input_value = get_input(["html","png"], "File type (html/png): ")
+                elif key2 == 'cv_type':
+                    input_value = get_input(["timeseriessplit","kfold"], "CV type (timeseriessplit/kfold): ")
+                elif key2 == 'max_rows':
+                    input_value = int(input("\nMax rows (number): "))
+                elif key2 == 'time_sensitive_column':
+                    if input_params['pre_params']['cv_type'] == 'timeseriessplit':
+                        if len(features_dict) == 0:
+                            print('No columns to select as time sensitive, switching to kfold')
+                            input_params['pre_params']['cv_type'] = 'kfold'
+                        print("\nWhich column would you like to use as the time sensitive column?")
+                        for i, header in features_dict.items():
+                            print(f"{i}. {header}")
+                        date_col_index = int(input("Enter the index of the time sensitive column: "))
+                        input_value = features_dict[date_col_index]
+                    else:
+                        input_value = '_'
+                elif key2 == 'upsampling':
+                    # input_value = get_input(["true","false"], "Upsampling (true/false): ")
+                    input_value = "false"
+                
+                # Fill in value in dict
                 input_params[key][key2] = input_value
 
 # convert the dictionary to a JSON string
 json_string = json.dumps(input_params, indent=4)
 
+# name config file
+config_name = input("\nName config file:")
+config_name = config_name.split(".")[0]
+
 # write the JSON string to a file
-with open("input_params.json", "w") as file:
+with open(f"input/configuration/{config_name}.json", "w") as file:
     file.write(json_string)
 
 # Print the resulting JSON
-print(json_string)
+print(json_string)
\ No newline at end of file
diff --git a/main.py b/scripts/main.py
similarity index 75%
rename from main.py
rename to scripts/main.py
index fe6c99a..3534b18 100644
--- a/main.py
+++ b/scripts/main.py
@@ -1,17 +1,21 @@
 # Load packages
 import logging
-import sys
+import pandas as pd
+import json
 
-from utils.modelling.main_modeler import *
+# Main modelling function
+from utils.modelling.main_modeler import make_model
 
 # The translations to SQL (grey as we refer to them dynamically)
 from utils.output_scripts import decision_tree_as_code
 from utils.output_scripts import decision_rule_as_code
 from utils.output_scripts import ebm_as_code
 
-from utils.helper_functions.config_handling import *
-from utils.helper_functions.parsing_arguments import *
-from utils.pre_process import *
+from utils.helper_functions.checks import checkInputDataHard
+
+from utils.helper_functions.config_handling import config_handling
+from utils.helper_functions.parsing_arguments import GetArgs
+from utils.pre_processing.pre_process import pre_process_kfold
 
 def main(args):
     # get given name from the first given argument
@@ -20,11 +24,12 @@ def main(args):
     # set logger
     logging.basicConfig(format='%(asctime)s %(message)s', filename=given_name+'/logging.log', level=logging.DEBUG)
     logging.getLogger('matplotlib.font_manager').disabled = True # ignore matplotlibs font warnings
+    logging.info(f'Script input arguments: \n{args}')
 
-    # path to data
-    data_ = pd.read_csv(args.data_path)
+    # Load in data
+    data = pd.read_csv(args.data_path)
 
-    # get target and features columns
+    # Get configuration file
     with open(args.configuration) as json_file:
         configuration = json.load(json_file)
 
@@ -34,17 +39,14 @@ def main(args):
     # Handle the configuration file
     target_col, feature_cols, model_params, pre_params, post_params = config_handling(configuration, logging)
 
-    # pre processing
-    # make copy (doesn't change anything but for future use)
-    data = data_.copy()
-
     # Log parameters
     logging.info(f'Configuration file content: \n{configuration}')
 
+    # Perform some basic checks
+    checkInputDataHard(data, configuration)
+
     # set model type based on target value
-    if data[target_col].nunique() == 1:
-        raise Exception("Target column needs more than 1 unique value")
-    elif (data[target_col].dtype == 'float') | ((data[target_col].dtype == 'int') & (data[target_col].nunique() > 10)):
+    if (data[target_col].dtype == 'float') | ((data[target_col].dtype == 'int') & (data[target_col].nunique() > 10)):
         model_type = 'regression'
     else:
         model_type = 'classification'
@@ -82,10 +84,12 @@ def main(args):
         argvals = None
     else:
         # (Dev) script is not being run through the terminal
+        # make sure pwd is the root folder (not in scripts)
+
         # Command line arguments used for testing
-        argvals = '--name trained_models/test ' \
-                  '--data_path input/data/example_multiclass_faults.csv ' \
-                  '--configuration input/configuration/example_multiclass_faults.json ' \
+        argvals = '--name ../trained_models/test ' \
+                  '--data_path ../input/data/example_binary_titanic.csv ' \
+                  '--configuration ../input/configuration/example_binary_titanic.json ' \
                   '--model ebm'.split() # example of passing test params to parser
 
         # settings
@@ -96,7 +100,7 @@ def main(args):
     random_seed = 42
 
     # Get arguments from the CLI
-    args = GetArgs(argvals)
+    args = GetArgs('main', argvals)
 
     # Run main with given arguments
     main(args)
\ No newline at end of file
diff --git a/scripts/setup.py b/scripts/setup.py
new file mode 100644
index 0000000..d938a0b
--- /dev/null
+++ b/scripts/setup.py
@@ -0,0 +1,24 @@
+import os
+import subprocess
+import sys
+
+VENV_NAME = ".ml2sql"
+
+# Create virtual environment
+if sys.platform.startswith("win"):
+    print(f"Creating virtual environment '{VENV_NAME}'...")
+    subprocess.run(["python", "-m", "venv", VENV_NAME], shell=True)
+else:
+    print(f"Creating virtual environment '{VENV_NAME}'...")
+    subprocess.run(["python3", "-m", "venv", VENV_NAME])
+
+# Install packages in virtual environment
+print(f"Installing packages in virtual environment...")
+if sys.platform.startswith("win"):
+    install_packages = os.path.join(VENV_NAME, "Scripts", "pip")
+else:
+    install_packages = os.path.join(VENV_NAME, "bin", "pip")
+
+subprocess.run([f"{install_packages}", "install", "-r", "requirements.txt"])
+
+print("Virtual environment created and packages installed successfully.")
\ No newline at end of file
diff --git a/utils/__init__.py b/scripts/utils/__init__.py
similarity index 100%
rename from utils/__init__.py
rename to scripts/utils/__init__.py
diff --git a/utils/feature_selection/correlations.py b/scripts/utils/feature_selection/correlations.py
similarity index 60%
rename from utils/feature_selection/correlations.py
rename to scripts/utils/feature_selection/correlations.py
index cbf5746..24746b6 100644
--- a/utils/feature_selection/correlations.py
+++ b/scripts/utils/feature_selection/correlations.py
@@ -7,6 +7,19 @@
 
 
 def plotClustermap(dfc, matrix_type, given_name, file_type, logging):
+    """
+    Plot a clustermap of a given correlation matrix/dataframe and save it as either a png or an html file.
+
+    Args:
+        dfc (pandas DataFrame): A correlation matrix/dataframe.
+        matrix_type (str): A string representing the type of matrix being plotted (e.g., 'correlation', 'covariance').
+        given_name (str): A string representing the name of the project being worked on.
+        file_type (str): A string representing the type of file to save the plot as (either 'png' or 'html').
+        logging (logging.Logger): A logging object used for logging events.
+
+    Returns:
+        None
+    """
 
     # using correlation matrix/dataframe as input
     labels = dfc.columns
@@ -105,15 +118,26 @@ def plotClustermap(dfc, matrix_type, given_name, file_type, logging):
         # or as html file
         fig.write_html(f'{given_name}/feature_info/{matrix_type}_clustermap.html')
 
-    print(f'Created and saved {matrix_type} clustermap')
-    logging.info(f'Created and saved {matrix_type} clustermap')
-
 def plotPearsonCorrelation(df, given_name, file_type, logging):
+    """
+    Creates a Pearson correlation matrix using numerical columns of a Pandas DataFrame,
+    and generates a clustermap plot of the matrix using the plotClustermap function.
+
+    Args:
+        df (Pandas DataFrame): The DataFrame containing the data to analyze.
+        given_name (str): The name of the project, used in the output file name.
+        file_type (str): The file type to save the output as ('png' or 'html').
+        logging (logging.Logger): The logger object to log progress and errors.
+
+    Returns:
+        None
+
+    Raises:
+        None
+
+    """
     # Numerical values
     # Creating pearson correlation matrix
-    print(f'Creating Pearson correlation matrix')
-    logging.info(f'Creating Pearson correlation matrix')
-
     data_ = df.select_dtypes(exclude=['object', 'category', 'bool']).copy()
 
     # if above dataframe is empty then skip function
@@ -131,11 +155,29 @@ def plotPearsonCorrelation(df, given_name, file_type, logging):
     # Plot matrix
     plotClustermap(data_corr, matrix_type, given_name, file_type, logging)
 
+    # success message
+    print(f'Created Pearson correlation matrix (for numerical features)')
+    logging.info(f'Created Pearson correlation matrix (for numerical features)')
+
 
 def cramers_corrected_stat(confusion_matrix):
-    """ calculate Cramers V statistic for categorial-categorial association.
-        uses correction from Bergsma and Wicher,
-        Journal of the Korean Statistical Society 42 (2013): 323-328
+    """
+    Calculates the corrected Cramer's V statistic for categorical-categorical association.
+
+    Parameters:
+    -----------
+    confusion_matrix: array-like, shape (n_categories, n_categories)
+        A confusion matrix representing the association between two categorical variables.
+
+    Returns:
+    --------
+    cramers_v: float
+        The corrected Cramer's V statistic for the given confusion matrix.
+
+    References:
+    -----------
+    1. Bergsma, Wicher, and Marcel A. Croon. "Marginal models for dependent, clustered, and longitudinal categorical data." Springer Science & Business Media, 2009.
+    2. Bergsma, Wicher, and Marcel A. Croon. "Reply to 'A note on the gamma statistic for measuring nominal association'." Journal of the Korean Statistical Society 42.3 (2013): 323-328.
     """
     chi2 = ss.chi2_contingency(confusion_matrix)[0]
     n = confusion_matrix.to_numpy().sum()
@@ -147,17 +189,48 @@ def cramers_corrected_stat(confusion_matrix):
     return np.sqrt(phi2corr / min( (kcorr-1), (rcorr-1)))
 
 def plotCramervCorrelation(df, given_name, file_type, logging):
+    """
+    Plot CramerV correlation matrix for nominal categorical values.
+
+    Parameters:
+    -----------
+    df : pandas.DataFrame
+        Dataframe containing categorical variables for which correlation is to be computed
+    given_name : str
+        Given name for the output file containing the correlation matrix plot
+    file_type : str
+        Type of the output file, e.g., png, pdf, etc.
+    logging : logging.Logger
+        Logging object for recording the progress and errors during runtime.
+
+    Returns:
+    --------
+    None
+
+    Raises:
+    -------
+    None
+
+    Notes:
+    ------
+    This function computes the Cramer's V correlation matrix for the categorical variables in the input dataframe. 
+    If the dataframe does not have at least two categorical variables, the function will skip the computation and return.
+    The correlation matrix is plotted using the seaborn clustermap function, which creates a heatmap of the correlation 
+    matrix, and performs hierarchical clustering of the variables.
+
+    References:
+    -----------
+    - Bergsma, W., & Wicher, M. (2013). 
+      Journal of the Korean Statistical Society, 42(3), 323-328.
+    """
     # Nominal Categorical values
     # Creating CramerV correlation matrix
-    print(f'Creating CramerV correlation matrix')
-    logging.info(f'Creating CramerV correlation matrix')
-
     data_ = df.select_dtypes(include=['category', 'object', 'bool']).copy()
 
     # if above dataframe is empty then skip function
     if data_.shape[1] < 2:
-        print(f'Less than 2 categorical variables found')
-        logging.info(f'Less than 2 categorical variables found')
+        print(f'Skip CramerV correlation matrix due to fewer than 2 categorical variables found')
+        logging.info(f'Skip CramerV correlation matrix due to fewer than 2 categorical variables found')
         return
 
     cols = data_.columns
@@ -182,3 +255,6 @@ def plotCramervCorrelation(df, given_name, file_type, logging):
     # Plot matrix
     plotClustermap(df_cramerv, matrix_type, given_name, file_type, logging)
 
+    # Success message
+    print(f'Created CramerV correlation matrix (for categorical values)')
+    logging.info(f'Created CramerV correlation matrix (for categorical values)')
diff --git a/utils/helper_functions/checks.py b/scripts/utils/helper_functions/checks.py
similarity index 52%
rename from utils/helper_functions/checks.py
rename to scripts/utils/helper_functions/checks.py
index c5b2134..2fd4e2e 100644
--- a/utils/helper_functions/checks.py
+++ b/scripts/utils/helper_functions/checks.py
@@ -1,8 +1,38 @@
-from sklearn.tree import _tree
 import pandas as pd
 import numpy as np
 
+
+def checkTargetHard(target):
+    if target.nunique() == 1:
+        raise Exception("Target column needs more than 1 unique value")
+
+def checkFeaturesHard(features):
+    featNullCount = features.isnull().sum()
+    nullf = featNullCount[featNullCount > 0]
+    if len(nullf) > 0:
+        raise Exception(f"NULL values not allowed, found the following: \n{nullf}")
+
+def checkInputDataHard(data, config):
+    """
+    Checks at start ensuring target and feature columns are good to go
+    """
+    checkTargetHard(data[config['target']])
+    checkFeaturesHard(data[config['features']])
+
 def checkAllClassesHaveLeafNode(clf):
+    """
+    Check if all classes are represented by a leaf node in a given decision tree classifier.
+
+    Parameters:
+    -----------
+    clf : sklearn.tree.DecisionTreeClassifier object
+        The decision tree classifier to be checked.
+
+    Returns:
+    --------
+    bool
+        True if all classes are represented by a leaf node, False otherwise.
+    """
     n_nodes = clf.tree_.node_count
     children_left = clf.tree_.children_left
     children_right = clf.tree_.children_right
diff --git a/utils/helper_functions/config_handling.py b/scripts/utils/helper_functions/config_handling.py
similarity index 73%
rename from utils/helper_functions/config_handling.py
rename to scripts/utils/helper_functions/config_handling.py
index fc99715..a00797e 100644
--- a/utils/helper_functions/config_handling.py
+++ b/scripts/utils/helper_functions/config_handling.py
@@ -1,7 +1,33 @@
 # Handle the configuration file
-import json
 
 def config_handling(configuration, logging):
+    """
+    Handles the configuration file and extracts necessary information.
+
+    Parameters:
+    -----------
+    configuration : dict
+        A dictionary containing configuration information such as target column, features columns, 
+        model related parameters, pre processing related parameters, and post modeling related parameters.
+
+    logging : logger object
+        A logger object to log the progress and debugging information.
+
+    Returns:
+    --------
+    tuple
+        A tuple containing the following:
+        - target_col : str
+            The name of the target column.
+        - feature_cols : list
+            A list containing the names of the feature columns.
+        - model_params : dict
+            A dictionary containing model related parameters.
+        - pre_params : dict
+            A dictionary containing pre processing related parameters.
+        - post_params : dict
+            A dictionary containing post modeling related parameters.
+    """
     # target column
     target_col = configuration['target']
 
diff --git a/scripts/utils/helper_functions/parsing_arguments.py b/scripts/utils/helper_functions/parsing_arguments.py
new file mode 100644
index 0000000..ac524ec
--- /dev/null
+++ b/scripts/utils/helper_functions/parsing_arguments.py
@@ -0,0 +1,64 @@
+import argparse
+
+def SetArgParser():
+    """
+    Create and configure an argument parser for command-line interface.
+
+    Returns:
+    parser: An argparse.ArgumentParser object.
+
+    Usage example:
+    parser = SetArgParser()
+    args = parser.parse_args()
+    """
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--name", type=str, help="Enter project name",
+                        nargs='?', default='no_name', const='no_name')
+    parser.add_argument("--data_path", type=str, help="Enter path to csv file",
+                        nargs='?', default='no_data', const='no_data')
+    parser.add_argument("--configuration", type=str, help="Enter path to json file",
+                        nargs='?')
+    parser.add_argument("--model_name", type=str, help="Enter model type",
+                        nargs='?', default='decision_tree', const='full')
+
+    return parser
+
+def GetArgs(script, argv):
+    """
+    Parse command-line arguments using argparse and return the parsed arguments.
+
+    Args:
+    - script (str): 'main' or 'test_model' dependent on what will run.
+    - argv (list of str): List of command-line arguments to parse.
+
+    Returns:
+    - args (argparse.Namespace): An object containing the parsed arguments.
+      The object has the following attributes:
+      - name (str): The name of the project. Default: 'no_name'.
+      - data_path (str): The path to the CSV file containing the data.
+      - configuration (str): The path to the JSON file containing the model configuration.
+      - model_name (str): The name of the model to use.
+    """
+    parser = argparse.ArgumentParser()
+
+    if script == 'main':
+      parser.add_argument("--name", type=str, help="Enter project name",
+                          nargs='?', default='no_name', const='no_name')
+      parser.add_argument("--data_path", type=str, help="Enter path to csv file",
+                          nargs='?', const='no_data')
+      parser.add_argument("--configuration", type=str, help="Enter path to json file",
+                          nargs='?', const='no_config')
+      parser.add_argument("--model_name", type=str, help="Enter model type",
+                          nargs='?', const='full')
+    
+    if script == 'test_model':
+      parser.add_argument("--model_path", type=str, help="Enter path to model",
+                          nargs='?', const='no_name')
+      parser.add_argument("--data_path", type=str, help="Enter path to csv file",
+                          nargs='?', const='no_data')
+      parser.add_argument("--destination_path", type=str, help="Enter path to folder for saving files",
+                          nargs='?', const='no_data')
+
+    args = parser.parse_args(argv)
+
+    return args
\ No newline at end of file
diff --git a/utils/modelling/__init__.py b/scripts/utils/modelling/__init__.py
similarity index 100%
rename from utils/modelling/__init__.py
rename to scripts/utils/modelling/__init__.py
diff --git a/utils/modelling/calibration.py b/scripts/utils/modelling/calibration.py
similarity index 63%
rename from utils/modelling/calibration.py
rename to scripts/utils/modelling/calibration.py
index 143b1c3..da66814 100644
--- a/utils/modelling/calibration.py
+++ b/scripts/utils/modelling/calibration.py
@@ -4,6 +4,31 @@
 from sklearn.linear_model import LogisticRegression
 
 def calibrateModel(model, X, y, logging, method='auto', final_model=False):
+    """
+    Calibrate a classification model.
+
+    Parameters
+    ----------
+    model : object
+        A trained classification model.
+    X : array-like or sparse matrix of shape (n_samples, n_features)
+        The input data.
+    y : array-like of shape (n_samples,)
+        The target values.
+    logging : object
+        A logging object for logging messages.
+    method : str, optional
+        The calibration method to use. If 'auto', the method is chosen automatically based on the size of the calibration dataset. If 'sigmoid', the logistic regression method is used. If 'isotonic', the isotonic regression method is used. Default is 'auto'.
+    final_model : bool, optional
+        If True, a calibrated model and calibration function are returned. Default is False.
+
+    Returns
+    -------
+    model_cal : object
+        A calibrated classifier.
+    cal_reg : object, optional
+        A calibration function. Returned only if `final_model` is True.
+    """
     # set calibration function depending on calibration data size
     if method == 'auto':
         if len(X) > 1000:
diff --git a/scripts/utils/modelling/main_modeler.py b/scripts/utils/modelling/main_modeler.py
new file mode 100644
index 0000000..2cd19ec
--- /dev/null
+++ b/scripts/utils/modelling/main_modeler.py
@@ -0,0 +1,187 @@
+import joblib
+import random
+import numpy as np 
+
+from sklearn.model_selection import train_test_split
+from utils.modelling.performance import *
+from utils.modelling.calibration import *
+
+# The actual algorithms (grey as we refer to them dynamically)
+from utils.modelling.models import ebm
+from utils.modelling.models import decision_rule
+from utils.modelling.models import decision_tree
+from utils.modelling.models import l_regression
+
+def make_model(given_name, datasets, model_name, model_type, model_params, post_params, logging):
+    """
+    Train and save a model, and generate performance plots.
+
+    Parameters
+    ----------
+    given_name : str
+        The name of the model.
+    datasets : dict
+        A dictionary containing the training and test datasets.
+    model_name : str
+        The name of the model to be used.
+    model_type : str
+        The type of the model (classification or regression).
+    model_params : dict
+        A dictionary containing the hyperparameters of the model.
+    post_params : dict
+        A dictionary containing the postprocessing parameters.
+    logging : logger
+        A logger object used for logging.
+
+    Returns
+    -------
+    A trained machine learning model that can be used to make predictions on new data.
+
+    """
+    # unpack datasets
+    X_train = datasets['cv_train']['X']
+    y_train = datasets['cv_train']['y']
+    X_test = datasets['cv_test']['X']
+    y_test = datasets['cv_test']['y']
+    X_all = datasets['final_train']['X']
+    y_all = datasets['final_train']['y']
+
+
+    # check if X is a list (CV should be applied in that case)
+    if isinstance(X_train, list):
+
+        y_test_pred_list = list()
+        y_test_prob_list = list()
+
+        # Save all trained models in a dictionary
+        clf_dict = {'test': {}}
+
+        for fold_id in range(len(X_train)):
+            print(f'Fold {fold_id} \nTrain model on test data')
+            logging.info(f'Fold {fold_id} \nTrain model on test data')
+
+            # Check if model needs to be calibrated
+            if post_params['calibration'] != 'false':
+                X_slice_train, X_cal, y_slice_train, y_cal = train_test_split(X_train[fold_id],y_train[fold_id], test_size=0.2,random_state=random.randint(0,100))
+            else:
+                X_slice_train, y_slice_train = X_train[fold_id], y_train[fold_id]
+
+            # Train the model
+            clf = globals()[model_name].trainModel(X_slice_train, y_slice_train, model_params, model_type, logging)
+
+            # Save model of this fold in a dict
+            clf_dict['test'][fold_id] = clf
+
+            if post_params['calibration'] != 'false':
+                clf = calibrateModel(clf, X_cal, y_cal, logging, method=post_params['calibration'], final_model=False)
+
+            # discrete predictions
+            y_test_pred_list.append(clf.predict(X_test[fold_id]))
+
+            if model_type == 'classification':
+            # probability predictions
+                # Binary classification
+                if len(clf.classes_) == 2:
+                    y_test_prob_list.append(clf.predict_proba(X_test[fold_id])[:,1])
+                elif len(clf.classes_) > 2:
+                    y_test_prob_list.append(clf.predict_proba(X_test[fold_id]))
+
+        # Merge list of prediction lists into one list
+        y_test_pred = np.concatenate(y_test_pred_list, axis=0)
+
+        if model_type == 'classification':
+            # Merge list of prediction probabilities lists into one list
+            y_test_prob = np.concatenate(y_test_prob_list, axis=0)
+
+    # If just regular train/test split has been applied
+    else:
+        clf = globals()[model_name].trainModel(X_train, y_train, model_params, model_type, logging)
+
+        # discrete prediction
+        y_test_pred = clf.predict(X_test)
+
+        if model_type == 'classification':
+            # probability prediction
+            y_test_prob = clf.predict_proba(X_test)[:,1]
+
+    # train model one last time on all samples (upsampled)
+    print('Train final model on all data')
+    logging.info('Train final model on all data')
+
+    # Check if model needs to be calibrated
+    if post_params['calibration'] != 'false':
+        X_all, X_cal, y_all, y_cal = train_test_split(X_all, y_all, test_size=0.2, random_state=123)
+
+    # Train final model
+    clf = globals()[model_name].trainModel(X_all, y_all, model_params, model_type, logging)
+
+    # Save in dict
+    clf_dict['final'] = clf
+
+    # Save target column name as part of model
+    clf.target = y_all.name
+
+    # Save model in pickled format
+    filename = f'{given_name}/model/{model_name}_{model_type}.sav'
+    joblib.dump(clf, open(filename, 'wb'))
+
+    # train set prediction of final model
+    y_all_pred = clf.predict(X_all)
+
+    if model_type == 'classification':
+    # probability predictions
+        # Binary classification
+        if len(clf.classes_) == 2:
+            y_all_prob = clf.predict_proba(X_all)[:,1]
+        elif len(clf.classes_) > 2:
+            y_all_prob = clf.predict_proba(X_all)
+
+
+    if post_params['calibration'] != 'false':
+        cal_clf, cal_reg = calibrateModel(clf, X_cal, y_cal, logging, method=post_params['calibration'], final_model=True)
+
+        # Save model in pickled format
+        filename = given_name + '/model/ebm_calibrated_{model_type}.sav'.format(model_type=model_type)
+        joblib.dump(cal_clf, open(filename, 'wb'))
+
+        # train set prediction of final model
+        y_all_pred = cal_clf.predict(X_all)
+
+        if model_type == 'classification':
+            # probability predictions
+            # Binary classification
+            if len(clf.classes_) == 2:
+                y_all_prob = clf.predict_proba(X_all)[:, 1]
+            elif len(clf.classes_) > 2:
+                y_all_prob = clf.predict_proba(X_all)
+
+    # concat y train, X test and X train to single lists
+    y_test_concat = np.concatenate(y_test, axis=0)
+    y_train_concat = np.concatenate(y_train, axis=0)
+    X_test_concat = np.concatenate(X_test, axis=0)
+    X_train_concat = np.concatenate(X_train, axis=0)
+
+    post_datasets = {
+        'X_test_list': X_test,
+        'X_train_list': X_train,
+        'X_test_concat': X_test_concat,
+        'X_train_concat': X_train_concat,
+        'X_all': X_all,
+        'y_train_list': y_train,
+        'y_train_concat': y_train_concat,
+        'y_all': y_all,
+        'y_all_pred': y_all_pred,
+        'y_test_concat': y_test_concat,
+        'y_test_pred': y_test_pred,
+        'y_test_list': y_test,
+    }
+
+    if model_type == 'classification':
+        post_datasets['y_all_prob'] = y_all_prob
+        post_datasets['y_test_prob'] = y_test_prob
+        post_datasets['y_test_prob_list'] = y_test_prob_list
+
+    # Performance and other post modeling plots
+    postModellingPlots(clf_dict, model_name, model_type, given_name, post_datasets, post_params, logging)
+
+    return clf
diff --git a/utils/modelling/models/__init__.py b/scripts/utils/modelling/models/__init__.py
similarity index 100%
rename from utils/modelling/models/__init__.py
rename to scripts/utils/modelling/models/__init__.py
diff --git a/utils/modelling/models/decision_rule.py b/scripts/utils/modelling/models/decision_rule.py
similarity index 53%
rename from utils/modelling/models/decision_rule.py
rename to scripts/utils/modelling/models/decision_rule.py
index 548fac4..b030bab 100644
--- a/utils/modelling/models/decision_rule.py
+++ b/scripts/utils/modelling/models/decision_rule.py
@@ -4,6 +4,33 @@
 from utils.modelling.performance import *
 
 def trainModel(X_train, y_train, params, model_type, logging):
+    """
+    Train a decision list classifier on the given training data.
+
+    Parameters
+    ----------
+    X_train : pandas.DataFrame
+        The input features for training the model.
+    y_train : pandas.Series
+        The target variable for training the model.
+    params : dict
+        A dictionary containing the parameters for the decision list classifier.
+    model_type : str
+        The type of model to train. Only 'classification' is currently supported.
+    logging : logging.Logger
+        A logger object to record messages during the training process.
+
+    Returns
+    -------
+    clf : DecisionListClassifier
+        The trained decision list classifier.
+
+    Raises
+    ------
+    ValueError
+        If an unsupported model type is provided.
+
+    """
     if 'feature_names' not in params.keys():
         params['feature_names'] = X_train.columns
     if model_type == 'classification':
@@ -19,6 +46,33 @@ def trainModel(X_train, y_train, params, model_type, logging):
     return clf
 
 def featureImportanceSave(clf, given_name, logging):
+    """
+    Save the feature importance graphs for the trained model.
+
+    Parameters
+    ----------
+    clf : object
+        The trained decision list classifier object.
+
+    given_name : str
+        The name of the file to which the feature importance graphs will be saved.
+
+    logging : logging.Logger
+        The logging object.
+
+    Returns
+    -------
+    None
+
+    Notes
+    -----
+    The feature importance graphs include the overall feature importance graph and the feature specific explanation graphs.
+
+    The feature importance graphs are saved in the following files:
+        - overall feature importance graph - {given_name}/decisions_rule.html
+        - feature specific explanation graphs - {given_name}/explain_{feature_name}.html
+
+    """
 
     clf_global = clf.explain_global()
 
diff --git a/scripts/utils/modelling/models/decision_tree.py b/scripts/utils/modelling/models/decision_tree.py
new file mode 100644
index 0000000..6a7ef31
--- /dev/null
+++ b/scripts/utils/modelling/models/decision_tree.py
@@ -0,0 +1,123 @@
+from sklearn import tree
+import pickle
+from utils.modelling.performance import *
+
+def trainModel(X_train, y_train, params, model_type, logging):
+    """
+    Trains a decision tree model using the given training data and parameters.
+
+    Parameters
+    ----------
+    X_train : pandas.DataFrame
+        Training features.
+    y_train : pandas.Series
+        Target variable.
+    params : dict
+        Parameters to configure the decision tree.
+    model_type : str
+        Type of the model. Can be 'regression' or 'classification'.
+    logging : logging.Logger
+        Logger object to record messages.
+
+    Returns
+    -------
+    clf : DecisionTreeRegressor or DecisionTreeClassifier
+        Trained decision tree model.
+
+    Raises
+    ------
+    ValueError
+        If an unsupported model type is provided.
+
+    Notes
+    -----
+    If the model type is 'regression', an instance of `DecisionTreeRegressor` is created
+    and trained. If the model type is 'classification', an instance of `DecisionTreeClassifier`
+    is created and trained. If an unsupported model type is provided, a `ValueError` is raised.
+
+    The trained decision tree model is returned.
+
+    Example:
+    --------
+    >>> from sklearn.datasets import load_breast_cancer
+    >>> from sklearn.model_selection import train_test_split
+    >>> from sklearn.tree import DecisionTreeClassifier
+    >>> X, y = load_breast_cancer(return_X_y=True)
+    >>> X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
+    >>> params = {'max_depth': 3, 'min_samples_split': 2}
+    >>> clf = trainModel(X_train, y_train, params, 'classification', logging)
+    """
+    if model_type == 'regression':
+        clf = tree.DecisionTreeRegressor(**params)
+    elif model_type == 'classification':
+        clf = tree.DecisionTreeClassifier(**params)
+    else:
+        print('Only regression or classification available')
+        logging.warning('Only regression or classification available')
+
+    clf.fit(X_train, y_train)
+    logging.info(f'Model params:\n {clf.get_params}')
+
+    print('Trained decision tree \n')
+    logging.info('Trained decision tree')
+
+    return clf
+
+def plotTreeStructureSave(clf, given_name):
+
+    plt.figure(figsize=(30,30))
+
+    tree.plot_tree(clf, fontsize=10, feature_names=clf.feature_names_in_, class_names=clf.classes_)
+    plt.savefig(f'{given_name}/tree_plot.png')
+
+    print('Tree structure plot saved')
+
+
+def featureImportanceSave(clf, given_name, file_type, logging):
+    """
+    Generates and saves a bar plot of feature importance using Plotly.
+    
+    Parameters:
+    -----------
+    clf: DecisionTreeClassifier or DecisionTreeRegressor object
+        The trained decision tree model.
+    given_name: str
+        The directory name where the plot should be saved.
+    file_type: str {'png', 'html'}
+        The type of file in which the plot should be saved.
+    logging: logging.Logger object
+        A logging object to log information or errors.
+        
+    Returns:
+    --------
+    None
+    
+    Raises:
+    -------
+    None
+    
+    Notes:
+    ------
+    - This function requires the Plotly package.
+    - The plot is saved in the given directory with the name 'gini_feature_importance.png' if 'file_type' is 'png', 
+      or 'gini_feature_importance.html' if 'file_type' is 'html'.
+    """
+    importance_df = pd.DataFrame({'importance':clf.feature_importances_, 'feature':clf.feature_names_in_}).sort_values('importance', ascending=True).reset_index(drop=True)
+    importance_non_zero = importance_df[importance_df['importance'] > 0]
+    plotly_fig = px.bar(importance_non_zero, x='importance', y='feature')
+
+    # Update size of figure
+    plotly_fig.update_layout(xaxis_title='Importance', yaxis_title='Feature',
+                      title=f'Feature importance',
+                      width=1000,
+                      height=800)
+
+    if file_type == 'png':
+        plotly_fig.write_image(f'{given_name}/gini_feature_importance.png')
+    elif file_type == 'html':
+        plotly_fig.write_html(f'{given_name}/gini_feature_importance.html')
+
+    print('Gini feature importance plot saved')
+
+def postModelPlots(clf, given_name, file_type, logging):
+    featureImportanceSave(clf, given_name, file_type, logging)
\ No newline at end of file
diff --git a/scripts/utils/modelling/models/ebm.py b/scripts/utils/modelling/models/ebm.py
new file mode 100644
index 0000000..3de1f59
--- /dev/null
+++ b/scripts/utils/modelling/models/ebm.py
@@ -0,0 +1,141 @@
+from interpret.glassbox import ExplainableBoostingClassifier, ExplainableBoostingRegressor
+import pickle
+import random
+from sklearn.model_selection import train_test_split
+from utils.modelling.performance import *
+from utils.modelling.calibration import *
+
+def trainModel(X_train, y_train, params, model_type, logging):
+    """
+    Trains an Explainable Boosting Machine (EBM) model on the given training data.
+
+    Parameters
+    ----------
+    X_train : pandas.DataFrame
+        Training data features.
+    y_train : pandas.Series
+        Training data target.
+    params : dict
+        Parameters for the EBM model.
+    model_type : str
+        Type of model to train, either 'regression' or 'classification'.
+    logging : logging.Logger
+        Logger instance for logging information.
+
+    Returns
+    -------
+    clf : EBM model
+        Trained EBM model.
+
+    Raises
+    ------
+    None
+
+    Notes
+    -----
+    This function trains an Explainable Boosting Machine (EBM) model on the given training data
+    using the provided parameters and model type. The trained model is returned and can be used
+    for prediction and feature importance analysis.
+
+    If 'feature_names' is not included in the provided params dictionary, the function will use
+    the column names of the X_train DataFrame as feature names.
+
+    If the model_type parameter is not 'regression' or 'classification', an error message is printed
+    and a warning is logged.
+
+    The function logs the model parameters and a message indicating that the EBM model has been
+    trained successfully.
+    """
+    if 'feature_names' not in params.keys():
+        params['feature_names'] = X_train.columns
+    if model_type == 'regression':
+        clf = ExplainableBoostingRegressor(**params)
+    elif model_type == 'classification':
+        clf = ExplainableBoostingClassifier(**params)
+    else:
+        print('Only regression or classification available')
+        logging.warning('Only regression or classification available')
+
+    clf.fit(X_train, y_train)
+    logging.info(f'Model params:\n {clf.get_params}')
+
+    print('Trained explainable boosting machine \n')
+    logging.info('Trained explainable boosting machine')
+
+    return clf
+
+
+def featureExplanationSave(clf, given_name, file_type, logging):
+    """
+    Saves feature-specific and overall feature importance graphs for a given Explainable Boosting Machine (EBM) model.
+
+    Parameters
+    ----------
+    clf : EBM model
+        Trained EBM model for which to generate feature importance graphs.
+    given_name : str
+        Name for the output files.
+    file_type : str
+        Type of file to save the output graphs, either 'png' or 'html'.
+    logging : logging.Logger
+        Logger instance for logging information.
+
+    Returns
+    -------
+    None
+
+    Raises
+    ------
+    None
+
+    Notes
+    -----
+    This function generates and saves feature-specific and overall feature importance graphs for a given
+    Explainable Boosting Machine (EBM) model. The overall feature importance graph is saved as a png or
+    html file, depending on the file_type parameter. The file is saved in a directory with the given_name
+    parameter as its name.
+
+    The function also generates and saves feature-specific explanation graphs for each feature group in the
+    EBM model. The file is saved in a directory with the given_name parameter as its name. Feature names
+    are reformatted by replacing certain characters with underscores.
+
+    The function logs a message indicating how many features have been saved.
+    """
+    
+    clf_global = clf.explain_global()
+
+    # Save overall feature importance graph
+    plotly_fig = clf_global.visualize()
+
+    if file_type == 'png':
+        plotly_fig.write_image('{given_name}/1_overall_feature_importance.png'.format(given_name=given_name))
+    elif file_type == 'html':
+        plotly_fig.write_html('{given_name}/1_overall_feature_importance.html'.format(given_name=given_name))
+
+    # Save feature specific explanation graphs
+    for index, value in enumerate(clf.term_features_):
+        plotly_fig = clf_global.visualize(index)
+
+        # Formatting feature name
+        # if combined feature create combined feature name
+        if len(value) == 2:
+            feature_name = f'{clf.feature_names[value[0]]}_x_{clf.feature_names[value[1]]}'
+        else:
+            feature_name = clf.feature_names[index]
+
+        chars = "\\`./ "
+        for c in chars:
+            if c in feature_name:
+                feature_name = feature_name.replace(c, "_")
+
+        if file_type == 'png':
+            plotly_fig.write_image(f'{given_name}/explain_{feature_name}.png')
+        elif file_type == 'html':
+            # or as html file
+            plotly_fig.write_html(f'{given_name}/explain_{feature_name}.html')
+
+    print('Explanation plots of {n_features} features saved'.format(n_features=index+1))
+    logging.info('Explanation plots of {n_features} features saved'.format(n_features=index+1))
+
+def postModelPlots(clf, given_name, file_type, logging):
+    featureExplanationSave(clf, given_name, file_type, logging)
diff --git a/utils/modelling/models/ebm.py b/scripts/utils/modelling/models/l_regression.py
similarity index 73%
rename from utils/modelling/models/ebm.py
rename to scripts/utils/modelling/models/l_regression.py
index c0db680..c510b76 100644
--- a/utils/modelling/models/ebm.py
+++ b/scripts/utils/modelling/models/l_regression.py
@@ -1,30 +1,30 @@
-from interpret.glassbox import ExplainableBoostingClassifier, ExplainableBoostingRegressor
 import pickle
-import random
-from sklearn.model_selection import train_test_split
+from interpret.glassbox import LinearRegression, LogisticRegression
 from utils.modelling.performance import *
-from utils.modelling.calibration import *
 
 def trainModel(X_train, y_train, params, model_type, logging):
-    if 'feature_names' not in params.keys():
-        params['feature_names'] = X_train.columns
+    # if 'feature_names' not in params.keys():
+    #     params['feature_names'] = X_train.columns
     if model_type == 'regression':
-        clf = ExplainableBoostingRegressor(**params)
+        clf = LinearRegression(**params)
+        clf_name = 'Linear regression'
     elif model_type == 'classification':
-        clf = ExplainableBoostingClassifier(**params)
+        clf = LogisticRegression(**params)
+        # Hard code classes_
+        clf.classes_ = list(set(y_train))
+        clf_name = 'Logistic regression'
     else:
         print('Only regression or classification available')
         logging.warning('Only regression or classification available')
 
     clf.fit(X_train, y_train)
-    logging.info(f'Model params:\n {clf.get_params}')
+    logging.info(f'Model non default params:\n {clf.kwargs}')
 
-    print('Trained explainable boosting machine \n')
-    logging.info('Trained explainable boosting machine')
+    print(f'Trained {clf_name.lower()} \n')
+    logging.info(f'Trained {clf_name.lower()}')
 
     return clf
 
-
 def featureExplanationSave(clf, given_name, file_type, logging):
 
     clf_global = clf.explain_global()
@@ -58,4 +58,4 @@ def featureExplanationSave(clf, given_name, file_type, logging):
     logging.info('Explanation plots of {n_features} features saved'.format(n_features=index+1))
 
 def postModelPlots(clf, given_name, file_type, logging):
-    featureExplanationSave(clf, given_name, file_type, logging)
+    featureExplanationSave(clf, given_name, file_type, logging)
\ No newline at end of file
diff --git a/scripts/utils/modelling/performance.py b/scripts/utils/modelling/performance.py
new file mode 100644
index 0000000..d031709
--- /dev/null
+++ b/scripts/utils/modelling/performance.py
@@ -0,0 +1,981 @@
+import pandas as pd
+import numpy as np
+from matplotlib import pyplot as plt
+import plotly.express as px
+import plotly.graph_objects as go
+from plotly.subplots import make_subplots
+import pingouin
+from sklearn.calibration import calibration_curve
+from sklearn.metrics import brier_score_loss
+from sklearn.metrics import classification_report
+from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score, mean_absolute_percentage_error, explained_variance_score, max_error, median_absolute_error, mean_squared_log_error, mean_poisson_deviance, mean_gamma_deviance
+
+
+# The actual algorithms (grey as we refer to them dynamically)
+from utils.modelling.models import ebm
+from utils.modelling.models import decision_rule
+from utils.modelling.models import decision_tree
+from utils.modelling.models import l_regression
+
+def plotConfusionMatrixStatic(given_name, y_true, y_pred, data_type, logging):
+    """
+    Plot and save a confusion matrix for given predictions.
+
+    Parameters
+    ----------
+    given_name : str
+        The name of the experiment or model.
+    y_true : array-like of shape (n_samples,)
+        True labels of the samples.
+    y_pred : array-like of shape (n_samples,)
+        Predicted labels of the samples.
+    data_type : str
+        The type of the data (train, val, test, etc.)
+    logging : logging object
+        The object to log the status and messages.
+
+    Returns
+    -------
+    None
+
+    Notes
+    -----
+    Saves the confusion matrix in {given_name}/performance/{data_type}_confusion_matrix.png
+
+    """
+
+    from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
+    fig, ax = plt.subplots(figsize=(10, 8))
+
+    disp = ConfusionMatrixDisplay.from_predictions(y_true, y_pred, ax=ax)
+    plt.xticks(rotation=90)
+    plt.savefig('{given_name}/performance/{data_type}_confusion_matrix.png'.format(given_name=given_name, data_type=data_type), bbox_inches='tight')
+
+    print(f'Created and saved confusion matrix for {data_type} data')
+    logging.info(f'Created and saved confusion matrix for {data_type} data')
+
+# Mainly meant for binary classification
+def plotConfusionMatrixSlider(given_name, y_true, y_prob, data_type, logging):
+    conf_matrices = []
+    steps=[]
+
+    threshold_list = np.arange(0.0, 1.05, 0.05)
+
+    for threshold in threshold_list:
+        y_pred = [1 if x > threshold else 0 for x in y_prob]
+        
+        df = pd.DataFrame({'y_true': y_true, 'y_pred': y_pred})
+
+        # confusion matrix
+        cm = []
+        for p in [0,1]:
+            new_row = []
+            for a in [1,0]:
+                new_row.append(len(df[(df['y_pred'] == p) & (df['y_true'] == a)])/len(y_pred))
+            cm.append(new_row)
+
+        # Add to list of matrices
+        conf_matrices.append(cm)
+
+        # Create sklearn summary report 
+        report = classification_report(y_true, y_pred, output_dict=True, zero_division=0)
+
+        # Create a list of lists for the matrix
+        metrics_matrix = []
+
+        # Iterate over the keys and values in the summary report dictionary
+        for key, value in report.items():
+            if key == 'accuracy':
+                metrics_matrix.append(['', '', '', '', ''])
+                metrics_matrix.append(['Accuracy', '', '', value, ''])
+            elif key.startswith(('macro', 'weighted')):
+                metrics_matrix.append([key.capitalize(), value['precision'], value['recall'], value['f1-score'], value['support']])
+            else:
+                row = [key, value['precision'], value['recall'], value['f1-score'], value['support']]
+                metrics_matrix.append(row)
+
+        # Formatting options
+        precision_format = '{:.3f}'
+        integer_format = '{:d}'
+        empty_string_format = ''
+
+        # Format the matrix elements
+        formatted_matrix = []
+        for row in metrics_matrix:
+            formatted_row = [precision_format.format(cell) if isinstance(cell, float)
+                            else integer_format.format(cell) if isinstance(cell, int)
+                            else empty_string_format if cell == ''
+                            else cell
+                            for cell in row]
+            formatted_matrix.append(formatted_row)
+
+        # transpose matrix
+        metrics = list(zip(*formatted_matrix))
+
+        # Add values to the different steps
+        steps.append(dict(method = "restyle",
+                    args = [{'z': [ cm ], #in the initial fig update z and text
+                        'text': [cm],
+                        'cells.values':[metrics]}],
+                        label=round(threshold,2),
+                        ))
+        
+        if threshold == 0.5:
+            conf_matrix = cm
+            table_metrics = metrics
+
+    # Manually add labels
+    labels = ["1", "0"]
+    labels_r = labels.copy()
+    labels_r.reverse()
+
+    # Make subplots top table bottom heatmap
+    fig = make_subplots(rows=1, cols=2,
+                        specs=[[{'type': 'table'},
+                            {'type': 'heatmap'}]])
+
+    # Create heatmap 
+    heatmap = go.Heatmap(x=labels, y=labels_r,
+                        z=conf_matrix,
+                        text=conf_matrix,
+                        texttemplate="%{text:.3f}",
+                        hovertemplate="<br>".join([
+                            "Predicted: %{y}",
+                            "Actual: %{x}",
+                            "Share of total cases: %{z}",
+                            "<extra></extra>"
+                        ])
+    )
+
+    # Create table
+    metrics_table = go.Table(header=dict(values=['          ', 'Precision', 'Recall', 'F1-Score', 'Support']),
+                    cells=dict(values=table_metrics) 
+                    )
+
+    # Add table and heatmap to figure 
+    fig.add_trace(metrics_table, row=1, col=1)
+    fig.add_trace(heatmap, row=1, col=2)
+
+
+
+    fig.update_layout(width=1000, height=600,
+                    xaxis_title="Actual",
+                    yaxis_title="Predicted",
+                    xaxis_type="category", 
+                    yaxis_type="category", 
+                    xaxis_side="top",
+                    title_text= "Normalized confusion matrix", 
+                    title_x=0.5,
+                    )
+
+    sliders = [dict(
+        active=10,
+        currentvalue={"prefix": "Threshold: "},
+        pad={"t": 50}, #sort it out later
+        steps=steps
+    )]
+
+    fig.update_layout(sliders=sliders)
+
+    fig.write_html(f'{given_name}/performance/{data_type}_confusion_matrix.html', auto_play=False)
+
+    print(f'Created and saved confusion matrix for {data_type} data')
+    logging.info(f'Created and saved confusion matrix for {data_type} data')
+
+
+def plotConfusionMatrix(given_name, y_true, y_prob, y_pred, file_type, data_type, logging):
+    # If html is wanted and binary classification
+    # Make confusion matrix plot with slider
+    if (file_type=='html') & (len(set(y_true))==2):
+        plotConfusionMatrixSlider(given_name, y_true, y_prob, data_type, logging)
+    
+    # Otherwise make 'simple' static confusion matrix plot
+    else:
+        plotConfusionMatrixStatic(given_name, y_true, y_pred, data_type, logging)
+
+def classificationReportSave(given_name, y_true, y_pred, data_type, logging):
+    """
+    Save the classification report as a CSV file.
+
+    Parameters
+    ----------
+    given_name : str
+        The name of the project or experiment.
+    y_true : array-like of shape (n_samples,)
+        Ground truth (correct) target values.
+    y_pred : array-like of shape (n_samples,)
+        Estimated targets as returned by a classifier.
+    data_type : str
+        The type of data being evaluated (e.g. 'train', 'test').
+    logging : logging.Logger
+        An instance of the logging.Logger class to record the event.
+
+    Returns
+    -------
+    None
+        This function does not return any value but saves the classification report as a CSV file.
+
+    """
+
+    from sklearn.metrics import classification_report
+
+    report = classification_report(y_true, y_pred, output_dict=True)
+    report_df = pd.DataFrame(report).transpose()
+    report_df.to_csv('{given_name}/performance/{data_type}_classification_report.csv'.format(given_name=given_name, data_type=data_type))
+
+    print('Classification report saved')
+    logging.info('Classification report saved')
+
+def plotClassificationCurve(given_name, y_true, y_prob, curve_type, data_type, logging):
+    """
+    Plots the ROC or Precision-Recall curve for a binary classification model, and saves the plot image.
+
+    Parameters:
+    -----------
+    given_name : str
+        Path where the plot image should be saved.
+    y_true : array-like of shape (n_samples,)
+        True binary labels for the samples.
+    y_prob : array-like of shape (n_samples,)
+        Estimated probabilities of the positive class for the samples.
+    curve_type : str
+        Type of curve to plot. Either 'ROC' for Receiver Operating Characteristic curve, or 'PR' for Precision-Recall curve.
+    data_type : str
+        Type of data being plotted. This string is included in the plot title.
+    logging : logging object
+        Logging object for recording progress and error messages.
+
+    Returns:
+    --------
+    mean_auc : float
+        Mean area under the curve (AUC) value for all folds or for the fitted model.
+    """
+    from sklearn.metrics import roc_curve, precision_recall_curve, auc
+    curve_type = curve_type.lower()
+
+    if curve_type == 'roc':
+        title = 'ROC curve - {data_type} data'.format(data_type=data_type)
+        xlabel = 'False Positive Rate'
+        ylabel = 'True Positive Rate (Recall)'
+        diagCor1 = [0, 1]
+        diagCor2 = [0, 1]
+
+        curve_ = 'roc'
+
+    elif ('recall' in curve_type and 'precision' in curve_type) or curve_type == 'pr':
+        title = 'Precision-Recall curve - {data_type}'.format(data_type=data_type)
+        xlabel = 'True Positive Rate (Recall)'
+        ylabel = 'Precision'
+
+        # add the random line (= share of positives in sample)
+        if isinstance(y_true, list):
+            pos_share = sum([sum(el) for el in y_true]) / sum([len(el) for el in y_true])
+        else:
+            pos_share = sum(y_true)/len(y_true)
+
+        diagCor1 = [0, 1]
+        diagCor2 = [pos_share, pos_share]
+
+        curve_ = 'pr'
+
+    # Create plot
+    fig = go.Figure()
+
+    # set axis range to 0 - 1
+    fig.update_layout(xaxis_range=[0, 1], yaxis_range=[0, 1],
+                      xaxis_title=xlabel, yaxis_title=ylabel,
+                      title=title,
+                      width=1000,
+                      height=800)
+
+    # Add diagonal reference line
+    fig.add_shape(type="line",
+                  xref="paper", yref="paper",
+                  x0=diagCor1[0], x1=diagCor1[1], y0=diagCor2[0], y1=diagCor2[1],
+                  line=dict(
+                      color="black",
+                      width=2,
+                      dash="dot",
+                  )
+    )
+
+    # Save the auc(s)
+    auc_list = list()
+
+    if isinstance(y_prob, list):
+        for fold_id in range(len(y_prob)):
+
+            if curve_ == 'roc':
+                xVar, yVar, thresholds = roc_curve(y_true[fold_id], y_prob[fold_id])
+            elif curve_ == 'pr':
+                yVar, xVar, thresholds = precision_recall_curve(y_true[fold_id], y_prob[fold_id])
+
+            # Calculate area under curve (AUC)
+            auc_list.append(auc(xVar, yVar))
+            fig.add_trace(go.Scatter(x=xVar, y=yVar, mode='lines', name=f'Fold {fold_id}'))
+    else:
+
+        if curve_ == 'roc':
+            xVar, yVar, thresholds = roc_curve(y_true, y_prob)
+        elif curve_ == 'pr':
+            yVar, xVar, thresholds = precision_recall_curve(y_true, y_prob)
+
+        # Calculate area under curve (AUC)
+        auc_list.append(auc(xVar, yVar))
+        fig.add_trace(go.Scatter(x=xVar, y=yVar, mode='lines', name=f'Fitted model'))
+
+    # add (average) auc in image
+    fig.add_annotation(x=0.5, y=0,
+                       text=f"Mean AUC: {np.mean(auc_list)}",
+                       showarrow=False,
+                       yshift=10)
+
+    fig.write_image(f'{given_name}/performance/{data_type}_{curve_type}_plot.png')
+
+    print(f'Created and saved {curve_type} plot for {data_type} data')
+    logging.info(f'Created and saved {curve_type} plot for {data_type} data')
+
+    return np.mean(auc_list)
+
+def plotCalibrationCurve(given_name, y_true, y_prob, data_type, logging):
+    """
+    Plot the calibration curve for a set of true and predicted values.
+
+    Parameters
+    ----------
+    given_name : str
+        The name of the given data.
+    y_true : array-like of shape (n_samples,)
+        The true target values.
+    y_prob : array-like of shape (n_samples,)
+        The predicted probabilities of target.
+    data_type : str
+        The type of the given data (train or test).
+    logging : logging.Logger
+        An instance of the `logging.Logger` class.
+
+    Returns
+    -------
+    None
+
+    Notes
+    -----
+    This function creates and saves a calibration plot that shows the relationship
+    between the true target values and the predicted probabilities of target. The plot
+    is saved in a subdirectory named `performance` inside the directory specified by
+    `given_name`.
+
+    The function uses the `calibration_curve` function from scikit-learn to calculate
+    the fraction of positives and the mean predicted value for a range of predicted
+    probabilities. The number of points in the plot is determined by the `n_bins`
+    parameter of `calibration_curve`.
+
+    The function also calculates the Brier score loss for each fold of cross-validation
+    and adds the mean Brier score loss to the plot.
+
+    Examples
+    --------
+    >>> y_true = [0, 1, 1, 0, 1, 0]
+    >>> y_prob = [0.2, 0.6, 0.7, 0.4, 0.8, 0.1]
+    >>> plotCalibrationCurve("data", y_true, y_prob, "test", logging.getLogger())
+    """
+
+    # Create plot
+    fig = go.Figure()
+
+    # set axis range to 0 - 1
+    fig.update_layout(xaxis_range=[0, 1], yaxis_range=[0, 1],
+                      xaxis_title='Predicted probability', yaxis_title='Share being positive',
+                      title=f'Calibration plot (reliability curve) - {data_type} data',
+                      width=1000,
+                      height=800)
+
+    # Add diagonal reference line
+    fig.add_shape(type="line",
+                  xref="paper", yref="paper",
+                  x0=0, x1=1, y0=0, y1=1,
+                  line=dict(
+                      color="black",
+                      width=2,
+                      dash="dot",
+                  )
+    )
+
+    # Save the brier score loss
+    bsl_list = list()
+
+    if isinstance(y_prob, list):
+        for fold_id in range(len(y_prob)):
+            # summaries actuals and predicted probs to (bins) number of points
+            fraction_of_positives, mean_predicted_value = \
+                calibration_curve(y_true[fold_id], y_prob[fold_id], n_bins=10, strategy='quantile')
+
+            # Calculate area under curve (AUC)
+            bsl_list.append(brier_score_loss(y_true[fold_id], y_prob[fold_id]))
+            fig.add_trace(go.Scatter(x=mean_predicted_value, y=fraction_of_positives, mode='markers+lines', name=f'Fold {fold_id}'))
+    else:
+        # summaries actuals and predicted probs to (bins) number of points
+        fraction_of_positives, mean_predicted_value = \
+            calibration_curve(y_true, y_prob, n_bins=10, strategy='quantile')
+
+        # Calculate area under curve (AUC)
+        bsl_list.append(brier_score_loss(y_true, y_prob))
+        fig.add_trace(go.Scatter(x=mean_predicted_value, y=fraction_of_positives, mode='markers+lines', name=f'Fitted model'))
+
+    # add (average) auc in image
+    fig.add_annotation(x=0.5, y=0,
+                       text=f"Mean Brier Score Loss: {np.mean(bsl_list)}",
+                       showarrow=False,
+                       yshift=10)
+
+    fig.write_image(f'{given_name}/performance/{data_type}_calibration_plot.png')
+
+    print(f'Created and saved calibration plot for {data_type} data')
+    logging.info(f'Created and saved calibration plot for {data_type} data')
+
+    return
+
+# for multiclass classification WIP
+def multiClassPlotCalibrationCurvePlotly(given_name, actuals, probs, title, bins=10):
+    """
+    Plot the calibration curve for a set of true and predicted values
+
+    :param actuals: true target value
+    :param probs: predicted probabilities per class as pandas dataframe
+    :param bins: how many bins to divide data in for plotting
+    :param savefig: boolean if plot should be saved
+    :param saveFileName: str path to which to save the plot
+    :return: calibration plot
+    """
+    calibration_df = pd.DataFrame(columns=['classification','fraction_of_positives','mean_predicted_value'])
+
+    # loop through different classes
+    for cl in probs.columns:
+        y_true = np.where(actuals == cl, 1, 0)
+        y_prob = probs[cl]
+
+        # summarise actuals and predicted probs to (bins) number of points
+        fraction_of_positives, mean_predicted_value = \
+            calibration_curve(y_true, y_prob, n_bins=bins, strategy='quantile')
+
+        calibration_df = calibration_df.append(pd.DataFrame({'classification':cl, 'fraction_of_positives':fraction_of_positives, 'mean_predicted_value':mean_predicted_value}), ignore_index=True)
+
+    # Create scatter plot
+    fig = px.scatter(calibration_df, x='mean_predicted_value', y='fraction_of_positives', title='Calibration plot (reliability curve)', color='classification',width=1200, height=800)
+
+
+    # Make trace be line plus dots
+    for l in range(len(fig.data)):
+        fig.data[l].update(mode='markers+lines')
+
+    # set axis range to 0 - 1
+    # fig.update_layout(xaxis_range=[-0.1,1.1], yaxis_range=[-0.1,1.1], xaxis_title='Predicted probability', yaxis_title='Fraction of positives')
+
+    # Add diagonal reference line
+    fig.add_shape(type="line",
+                  xref="paper", yref="paper",
+                  x0=0, x1=1, y0=0, y1=1,
+                  line=dict(
+                      color="black",
+                      width=2,
+                      dash="dot",
+                  )
+    )
+
+    fig.show(renderer='browser')
+
+def plotProbabilityDistribution(given_name, y_true, y_prob, data_type, logging):
+    """
+    Plot the probability distribution of true and predicted values.
+
+    Parameters
+    ----------
+    given_name : str
+        Given name to create and save plot.
+    y_true : array-like of shape (n_samples,)
+        True target value.
+    y_prob : array-like of shape (n_samples,)
+        Predicted probabilities.
+    data_type : str
+        Type of data, such as "train", "test", or "validation".
+    logging : logger object
+        Logging object for printing errors and warnings.
+
+    Returns
+    -------
+    None
+    """
+    import plotly.figure_factory as ff
+    df = pd.DataFrame({'actuals': y_true, 'prob': y_prob})
+    do = df[df['actuals']==1]['prob']
+    dont = df[df['actuals']==0]['prob']
+
+    # Catching any kind of exception
+    try:
+        # Create distplot with custom bin_size
+        fig = ff.create_distplot([do,dont], ['1','0'], colors=['green','red'], bin_size=.01)
+    except Exception as e:
+        print(f'Could not create distribution plot because of \n{e}')
+        logging.info(f'Could not create distribution plot because of \n{e}')
+        return
+
+    # Update size of figure
+    fig.update_layout(xaxis_title='Predicted probability', yaxis_title='Frequency',
+                      title=f'Distribution plot - {data_type} data',
+                      width=1000,
+                      height=800)
+
+    # fig.show(renderer='browser')
+
+    fig.write_image(f'{given_name}/performance/{data_type}_distribution_plot.png')
+
+    print(f'Created and saved probability distribution plot')
+    logging.info(f'Created and saved probability distribution plot')
+
+    return
+
+def plotDistribution(given_name, groups, values, data_type, logging):
+    """
+    Create a plotly histogram figure displaying the distribution of values for each group.
+
+    Parameters
+    ----------
+    given_name: str
+    A string representing the name of the plot to be created.
+    groups: array-like of shape (n_samples,)
+    An array-like object containing the group assignments of each sample.
+    values: array-like of shape (n_samples,)
+    An array-like object containing the values to be plotted.
+    data_type: str
+    A string representing the type of data being plotted (e.g., training, validation, test).
+    logging: logger object
+    A logger object that can be used to record messages.
+
+    Returns
+    -------
+    None
+
+    Raises
+    ------
+    Exception: If the histogram cannot be created.
+
+    Example
+    -------
+    >>> groups = [0, 0, 1, 1, 2, 2]
+    >>> values = [1, 2, 2, 3, 3, 3]
+    >>> plotDistribution("myplot", groups, values, "training", logging.getLogger())
+    """
+    df = pd.DataFrame({'groups': groups, 'values': values})
+    
+    min_value = min(values)
+    max_value = max(values)
+    bin_size = (max_value - min_value) / 50
+    
+    colors = ['red','green','blue','purple','orange']
+
+    # Catching any kind of exception
+    try:
+        # Create distplot 
+        fig = go.Figure()
+        for g in range(len(set(groups))):
+            X = df[df['groups']==g]['values']
+            fig.add_trace(go.Histogram(
+                x=X,
+                histnorm='probability density',
+                name=str(g), # name used in legend and hover labels
+                xbins=dict( # bins used for histogram
+                    start=min_value,
+                    end=max_value,
+                    size=bin_size
+                ),
+                marker_color=colors[g],
+                opacity=0.75
+            ))
+    except Exception as e:
+        print(f'Could not create distribution plot because of \n{e}')
+        logging.info(f'Could not create distribution plot because of \n{e}')
+        return
+
+    fig.update_layout(
+        xaxis_title='Predicted probability',
+        yaxis_title='Probability density',
+        title=f'Distribution plot - {data_type} data',
+        width=1000,
+        height=800,
+        bargap=0.2, # gap between bars of adjacent location coordinates
+        bargroupgap=0.1 # gap between bars of the same location coordinates
+    )
+
+    fig.write_image(f'{given_name}/performance/{data_type}_distribution_plot.png')
+
+    print(f'Created and saved probability distribution plot')
+    logging.info(f'Created and saved probability distribution plot')
+
+    return
+
+"""
+Metrics and plot which are related to regression
+"""
+
+def plotYhatVsYSave(given_name, y_true, y_pred, data_type, logging):
+    """
+    Plots a scatter plot of predicted values (y_pred) against true values (y_true)
+    and saves the plot in a specified directory.
+
+    Parameters:
+    -----------
+    given_name : str
+        Name of the directory where the plot will be saved.
+    y_true : array-like of shape (n_samples,)
+        Ground truth (correct) target values.
+    y_pred : array-like of shape (n_samples,)
+        Estimated targets as returned by a classifier.
+    data_type : str
+        Type of data, e.g. train, validation, or test.
+    logging : logging.Logger
+        Logging instance to log information about the execution.
+
+    Returns:
+    --------
+    None
+    """
+
+    plot_df = pd.DataFrame({'y_true':y_true, 'y_pred':y_pred})
+    fig = px.scatter(plot_df, 'y_true', 'y_pred')
+    fig.write_image('{given_name}/performance/{data_type}_scatter_yhat_vs_y.png'.format(given_name=given_name, data_type=data_type))
+
+    print(f'Scatter plot of yhat vs y saved for {data_type}')
+    logging.info(f'Scatter plot of yhat vs y saved for {data_type}')
+
+def plotQuantileError(given_name, y_true, y_pred, data_type, logging):
+
+    # Add prediction error
+    plot_df = pd.DataFrame({'y_true':y_true, 'y_pred':y_pred, 'y_diff':y_pred - y_true})
+
+    # Make 20 quantiles
+    plot_df['quantile'] = pd.qcut(plot_df['y_true'], 20, duplicates='drop')
+
+    # Sort by quantile and make quantile colunn a string for plotting
+    plot_df = plot_df.iloc[plot_df['quantile'].cat.codes.argsort()]
+    plot_df['quantile'] = plot_df['quantile'].astype(str)
+
+    fig = px.box(plot_df, 'quantile', 'y_diff')
+    fig.update_layout(
+        xaxis_title='Actual',
+        yaxis_title='Prediction - actual',
+        title=f'Quantile error plot - {data_type} data',
+        width=1000,
+        height=800
+    )
+
+    fig.write_image('{given_name}/performance/{data_type}_quantile_error_plot.png'.format(given_name=given_name, data_type=data_type))
+
+    print(f'Quantile error plot saved for {data_type}')
+    logging.info(f'Quantile error plot saved for {data_type}')
+
+def regressionMetricsTable(given_name, y_true, y_pred, X_all, data_type, logging):
+    y_true_pos = np.clip(y_true, 0, None)
+    y_pred_pos = np.clip(y_pred, 0, None)
+
+    # Calculate metrics
+    mae = mean_absolute_error(y_true, y_pred)
+    mse = mean_squared_error(y_true, y_pred)
+    rmse = mse ** 0.5
+    r2 = r2_score(y_true, y_pred)
+    adj_r2 = 1 - (1 - r2) * (len(y_true) - 1) / (len(y_true) - X_all.shape[1] - 1)
+    mape = mean_absolute_percentage_error(y_true, y_pred)
+    evs = explained_variance_score(y_true, y_pred)
+    me = max_error(y_true, y_pred)
+    medae = median_absolute_error(y_true, y_pred)
+    msle = mean_squared_log_error(y_true_pos, y_pred_pos)
+    rmsle_value = msle ** 0.5
+
+    header = ["Metric"]
+
+    # Define the rows of the table
+    rows = [
+        ["Mean Absolute Error (MAE)"],
+        ["Mean Squared Error (MSE)"],
+        ["Root Mean Squared Error (RMSE)"],
+        ["R-squared (R2)"],
+        ["Adjusted R-squared (Adj R2)"],
+        ["Mean Absolute Percentage Error (MAPE)"],
+        ["Explained Variance Score (EVS)"],
+        ["Max Error"],
+        ["Median Absolute Error (MedAE)"],
+        ["Mean Squared Log Error (MSLE)"],
+        ["Root Mean Squared Log Error (RMSLE)"]
+    ]
+
+    # Create the table
+    fig = go.Figure(data=[go.Table(header=dict(values=header), cells=dict(values=rows))])
+
+    # Add metrics data
+    metric_values = [
+        [mae],
+        [mse],
+        [rmse],
+        [r2],
+        [adj_r2],
+        [mape],
+        [evs],
+        [me],
+        [medae],
+        [msle],
+        [rmsle_value]
+    ]
+
+    metric_names = [row for row in rows]
+
+    # Add metrics data as a new column
+    cells_values = [metric_names, metric_values]
+    column_names = ["Metric", "Value"]
+
+    fig.add_trace(go.Table(
+        header=dict(values=column_names),
+        cells=dict(values=cells_values, format=["",".3"])
+    ))
+
+    # Update table layout
+    fig.update_layout(
+        title="Regression Performance Metrics",
+        margin=dict(l=10, r=10, t=50, b=10),
+    )
+
+    # Save table
+    fig.write_image(f'{given_name}/performance/{data_type}_regression_metrics.png')
+
+    print(f'Created and saved regression metrics table')
+    logging.info(f'Created and saved regression metrics table')
+
+    return
+    
+
+# To be used for feature exploration (TODO add more colors)
+def plotDistributionViolin(given_name, groups, values, data_type, logging):
+    """
+    Plot violin plot of distribution of values across different groups.
+
+    Parameters
+    ----------
+    given_name : str
+    Name of the figure to be saved.
+
+    groups : array_like
+    Group labels of the data.
+
+    values : array_like
+    Values of the data.
+
+    data_type : str
+    Type of the data.
+
+    logging : logging.Logger
+    Logger object for logging messages.
+
+    Returns
+    -------
+    None
+    Returns nothing, but saves the violin plot figure as a png file.
+
+    Notes
+    -----
+    This function uses the Plotly library to create the violin plot. It first creates a pandas DataFrame from the input data, and then uses a for loop to create a violin plot for each group. The resulting figure is saved as a png file in the directory specified by given_name.
+
+    Examples
+    --------
+    >>> plotDistributionViolin('my_figure', [0, 1, 0, 1], [2.5, 3.7, 4.2, 5.3], 'continuous', logging)
+    """
+    df = pd.DataFrame({'groups': groups, 'values': values})
+    
+    colors = ['red','green','blue','purple','orange']
+
+    # Catching any kind of exception
+    try:
+        # Create Violin plot 
+        fig = go.Figure()
+        for g in range(len(set(groups))):
+            X = df[df['groups']==g]['values']
+            fig.add_trace(go.Violin(
+                x=X,
+                name=str(g), # name used in legend and hover labels
+                marker_color=colors[g],
+                opacity=0.75
+            ))
+    except Exception as e:
+        print(f'Could not create distribution plot because of \n{e}')
+        logging.info(f'Could not create distribution plot because of \n{e}')
+        return
+
+    fig.update_layout(
+        xaxis_title='Value',
+        yaxis_title='Group',
+        title=f'Distribution plot - {data_type} data',
+        width=1000,
+        height=800
+    )
+
+    fig.write_image(f'{given_name}/feature_info/{feature_name}_class_distributions.png')
+
+    print(f'Created and saved feature distribution plot')
+    logging.info(f'Created and saved feature distribution plot')
+
+    return
+
+# TODO below function would need to have the test models available and make the calucations using those
+# So that the X_test is actually not seen yet by the model 
+# Currently the clf['final'] is the model trained over all data
+def scorePartialCorrelation(given_name, clf, post_datasets):
+    for fold_id in range(len(post_datasets['X_test_list'])):
+        clf_fold = clf['test'][fold_id]
+
+        X_test = post_datasets['X_test_list'][fold_id]
+        X_train = post_datasets['X_train_list'][fold_id]
+        y_test = post_datasets['y_test_list'][fold_id]
+        y_train = post_datasets['y_train_list'][fold_id]
+
+        sl_test = clf_fold.predict_and_contrib(X_test, output='labels')
+        X_test_scores = pd.DataFrame(sl_test[1], columns=clf_fold.feature_names)
+        y_test_label = pd.Series(y_test, name='target')
+
+        sl_train = clf_fold.predict_and_contrib(X_train, output='labels')
+        X_train_scores = pd.DataFrame(sl_train[1], columns=clf_fold.feature_names)
+        y_train_label = pd.Series(y_train, name='target')
+
+    # Define function for partial correlation
+    def partial_correlation(X, y):
+        out = pd.Series(index=X.columns, dtype=float)
+        for feature_name in X.columns:
+            out[feature_name] = pingouin.partial_corr(
+                data=pd.concat([X, y], axis=1).astype(float),
+                x=feature_name,
+                y=y.name,
+                x_covar=[f for f in X.columns if f != feature_name]
+            ).loc['pearson', 'r']
+        return out
+
+    parscore_test = partial_correlation(X_test_scores, y_test_label)
+    parscore_train = partial_correlation(X_train_scores, y_train_label)
+    parscore_diff = pd.Series(parscore_test - parscore_train, name = 'parscore_diff')
+
+    # Plot parshap
+    plotmin, plotmax = min(parscore_train.min(), parscore_test.min()), max(parscore_train.max(), parscore_test.max())
+    plotbuffer = .05 * (plotmax - plotmin)
+    fig, ax = plt.subplots()
+    if plotmin < 0:
+        ax.vlines(0, plotmin - plotbuffer, plotmax + plotbuffer, color='darkgrey', zorder=0)
+        ax.hlines(0, plotmin - plotbuffer, plotmax + plotbuffer, color='darkgrey', zorder=0)
+    ax.plot(
+        [plotmin - plotbuffer, plotmax + plotbuffer], [plotmin - plotbuffer, plotmax + plotbuffer],
+        color='darkgrey', zorder=0
+    )
+    sc = ax.scatter(
+        parscore_train, parscore_test,
+        edgecolor='grey',
+        #   c=[5]*16, 
+          s=50, 
+          cmap=plt.cm.get_cmap('Reds'))
+    ax.set(title='Partial correlation bw Score and target...', xlabel='... on Train data', ylabel='... on Test data')
+    cbar = fig.colorbar(sc)
+    cbar.set_ticks([])
+    for txt in parscore_train.index:
+        ax.annotate(txt, (parscore_train[txt], parscore_test[txt] + plotbuffer / 2), ha='center', va='bottom')
+    fig.savefig(f'{given_name}/performance/PartialCorrelation.png', dpi=300, bbox_inches="tight")
+
+
+def postModellingPlots(clf, model_name, model_type, given_name, post_datasets, post_params, logging):
+    # Performance and other post modeling plots
+    # unpack dict
+    X_all = post_datasets['X_all']
+    y_all, y_all_pred = post_datasets['y_all'], post_datasets['y_all_pred']
+    y_test_concat, y_test_pred = post_datasets['y_test_concat'], post_datasets['y_test_pred']
+    y_test_list = post_datasets['y_test_list']
+
+    if model_type == 'classification':
+        # probabilities only for classifications
+        y_test_prob, y_all_prob, y_test_prob_list = post_datasets['y_test_prob'], post_datasets['y_all_prob'], post_datasets['y_test_prob_list']
+
+        # Threshold dependant
+        plotConfusionMatrix(given_name, y_all, y_all_prob, y_all_pred, post_params['file_type'], data_type='final_train', logging=logging)
+        plotConfusionMatrix(given_name, y_test_concat, y_test_prob, y_test_pred, post_params['file_type'], data_type='test', logging=logging)
+
+        if len(clf['final'].classes_) == 2:
+            # Also create pr curve for class 0
+            y_all_neg = np.array([1 - j for j in list(y_all)])
+            y_all_prob_neg = np.array([1 - j for j in list(y_all_prob)])
+
+            y_test_list_neg = [[1 - j for j in i] for i in y_test_list]
+            y_test_prob_list_neg = [[1 - j for j in i] for i in y_test_prob_list]
+
+            # Threshold independant
+            plotClassificationCurve(given_name, y_all, y_all_prob, curve_type='roc', data_type='final_train', logging=logging)
+            plotClassificationCurve(given_name, y_test_list, y_test_prob_list, curve_type='roc', data_type='test', logging=logging)
+
+            plotClassificationCurve(given_name, y_all, y_all_prob, curve_type='pr', data_type='final_train_class1', logging=logging)
+            plotClassificationCurve(given_name, y_all_neg, y_all_prob_neg, curve_type='pr', data_type='final_train_class0', logging=logging)
+
+            plotClassificationCurve(given_name, y_test_list, y_test_prob_list, curve_type='pr', data_type='test_data_class1', logging=logging)
+            plotClassificationCurve(given_name, y_test_list_neg, y_test_prob_list_neg, curve_type='pr', data_type='test_data_class0', logging=logging)
+
+            plotCalibrationCurve(given_name, y_all, y_all_prob, data_type='final_train', logging=logging)
+            plotCalibrationCurve(given_name, y_test_list, y_test_prob_list, data_type='test', logging=logging)
+
+            plotProbabilityDistribution(given_name, y_all, y_all_prob, data_type='final_train', logging=logging)
+            plotProbabilityDistribution(given_name, y_test_concat, y_test_prob, data_type='test', logging=logging)
+
+            # TODO fix this
+            # scorePartialCorrelation(given_name, clf, post_datasets)
+
+        # If multiclass classification
+        elif len(clf['final'].classes_) > 2:
+            # loop through classes
+            for c in clf['final'].classes_:
+
+                # stating the class
+                print(f'\nFor class {c}:')
+                logging.info(f'\nFor class {c}:')
+
+                # creating a list of all the classes except the current class
+                other_class = [x for x in clf['final'].classes_ if x != c]
+
+                # Get index of selected class in clf['final'].classes_
+                class_index = list(clf['final'].classes_).index(c)
+
+                # marking the current class as 1 and all other classes as 0
+                y_test_list_ova = [[0 if x in other_class else 1 for x in fold_] for fold_ in y_test_list]
+                y_test_prob_list_ova = [[x[class_index] for x in fold_] for fold_ in y_test_prob_list]
+
+                # concatonate probs together to one list for distribution plot
+                y_test_ova = np.concatenate(y_test_list_ova, axis=0)
+                y_test_prob_ova = np.concatenate(y_test_prob_list_ova, axis=0)
+
+                y_all_ova = [0 if x in other_class else 1 for x in y_all]
+                y_all_prob_ova = [x[class_index] for x in y_all_prob]
+
+
+                # Threshold independant
+                # plotClassificationCurve(given_name, y_all_ova, y_all_prob_ova, curve_type='roc', data_type=f'final_train_class_{c}', logging=logging)
+                plotClassificationCurve(given_name, y_test_list_ova, y_test_prob_list_ova, curve_type='roc', data_type=f'test_class_{c}', logging=logging)
+
+                # plotClassificationCurve(given_name, y_all_ova, y_all_prob_ova, curve_type='pr', data_type='final_train_class1', logging=logging)
+                plotClassificationCurve(given_name, y_test_list_ova, y_test_prob_list_ova, curve_type='pr', data_type=f'test_class_{c}', logging=logging)
+
+                # multiClassPlotCalibrationCurvePlotly(given_name, y_all, pd.DataFrame(y_all_prob, columns=clf['final'].classes_), title='fun')
+                plotCalibrationCurve(given_name, y_test_list_ova, y_test_prob_list_ova, data_type=f'test_class_{c}', logging=logging)
+
+                # plotProbabilityDistribution(given_name, y_all_ova, y_all_prob_ova, data_type='final_train', logging=logging)
+                plotProbabilityDistribution(given_name, y_test_ova, y_test_prob_ova, data_type=f'test_class_{c}', logging=logging)
+
+    # if regression
+    elif model_type == 'regression':
+        plotYhatVsYSave(given_name, y_test_concat, y_test_pred, data_type='test', logging=logging)
+        plotYhatVsYSave(given_name, y_all, y_all_pred, data_type='final_train', logging=logging)
+
+        plotQuantileError(given_name, y_test_concat, y_test_pred, data_type='test', logging=logging)
+        plotQuantileError(given_name, y_all, y_all_pred, data_type='final_train', logging=logging)
+
+        regressionMetricsTable(given_name, y_test_concat, y_test_pred, X_all, data_type='test', logging=logging)
+        regressionMetricsTable(given_name, y_all, y_all_pred, X_all, data_type='final_train', logging=logging)
+
+    # Post modeling plots, specific per model but includes feature importance among others
+    globals()[model_name].postModelPlots(clf['final'], given_name + '/feature_importance', post_params['file_type'], logging)
\ No newline at end of file
diff --git a/utils/output_scripts/decision_rule_as_code.py b/scripts/utils/output_scripts/decision_rule_as_code.py
similarity index 100%
rename from utils/output_scripts/decision_rule_as_code.py
rename to scripts/utils/output_scripts/decision_rule_as_code.py
diff --git a/utils/output_scripts/decision_tree_as_code.py b/scripts/utils/output_scripts/decision_tree_as_code.py
similarity index 100%
rename from utils/output_scripts/decision_tree_as_code.py
rename to scripts/utils/output_scripts/decision_tree_as_code.py
diff --git a/utils/output_scripts/ebm_as_code.py b/scripts/utils/output_scripts/ebm_as_code.py
similarity index 51%
rename from utils/output_scripts/ebm_as_code.py
rename to scripts/utils/output_scripts/ebm_as_code.py
index 0d26946..b6d5341 100644
--- a/utils/output_scripts/ebm_as_code.py
+++ b/scripts/utils/output_scripts/ebm_as_code.py
@@ -3,6 +3,136 @@
 from interpret.glassbox import ExplainableBoostingClassifier
 from contextlib import redirect_stdout
 
+def ReduceSingleFeature(df_):
+    # For safety
+    df = df_.copy()
+
+    # drop any duplicates (which should not exist)
+    df_feat_types = df[['feature','feat_type']].drop_duplicates().reset_index(drop=True).copy()
+
+    # Set index
+    df = df.set_index(['feature'])
+
+    # Explode per column 
+    lookup_df = df[['feat_bound', 'score',
+    'score_lower_bound', 'score_upper_bound']].apply(lambda x: x.apply(pd.Series).stack()).reset_index().drop('level_1', axis=1)
+
+    # add feature type back in
+    lookup_df = lookup_df.merge(df_feat_types, on='feature', how='left')
+
+    # sort_values but None will be in the bottom
+    lookup_df.sort_values(['feature','feat_bound'], inplace=True, ascending=False)
+
+    # So reverse the ordering
+    lookup_df = lookup_df.iloc[::-1]
+
+    lookup_df['nr_features'] = 1
+
+    # being able for score arrays to be grouped
+    # for categorical features do not group by score hence we include value in score_hash
+    lookup_df['pre_hash_column'] = np.where(lookup_df['feat_type'] == 'categorical'
+                                            , lookup_df['score'].astype(str) + lookup_df['feat_bound'].astype(str)
+                                            , lookup_df['score'].astype(str))
+    lookup_df['score_hash'] = pd.util.hash_pandas_object(lookup_df['pre_hash_column'], index=False)
+    hash_lookup_df = lookup_df[['score_hash','score']].drop_duplicates('score_hash').reset_index(drop=True).copy()
+
+    # to group by score but keep them in order (not grouping similar scores for different feature values)
+    adj_check = (lookup_df.score_hash != lookup_df.score_hash.shift()).cumsum()
+
+    lookup_df['adj_score'] = adj_check
+    lookup_df_simple = lookup_df.groupby(['feature', 'feat_type', 'score_hash', 'adj_score']).agg({'feat_bound':'first'
+                                                ,'score_lower_bound':'first'
+                                                ,'score_upper_bound':'last'
+                                                }, sort=False).reset_index(drop=False)
+
+    lookup_df_grouped = lookup_df_simple.merge(hash_lookup_df, on='score_hash', how='left')
+    # ordering sets None in last place
+    lookup_df_grouped.sort_values(['feature', 'feat_bound'], inplace=True, ascending=False)
+    # So reverse the ordering
+    lookup_df_grouped = lookup_df_grouped.iloc[::-1].reset_index(drop=True)
+
+    return lookup_df_grouped
+
+
+def RestructureReduceInteractions(df_):
+
+    # for safety
+    df = df_.copy()
+
+    # Restructure
+    # loop over rows
+    for index, row in df.iterrows():
+        # Swap features position if first feature has more bound values then second feature
+        if row['score'].shape[0] > row['score'].shape[1]:
+
+            # loop over column needed to be transposed
+            for c in ['feature', 'feat_bound', 'score', 'feat_type']:
+                if isinstance(df.loc[index,c], list):
+                    df.at[index,c] = [df.at[index,c][1], df.at[index,c][0]]
+                elif isinstance(df.loc[index,c], np.ndarray):
+                    df.at[index,c] = np.transpose(df.at[index,c])
+
+    # Reduce
+
+    # Split list valued columns into two columns
+    df[['feat_1', 'feat_2']] = df['feature'].apply(pd.Series)
+    df[['feat_bound_1', 'feat_bound_2']] = df['feat_bound'].apply(pd.Series)
+    df[['feat_type_1', 'feat_type_2']] = df['feat_type'].apply(pd.Series)
+
+    # Drop unused column 
+    df = df.drop(['feat_bound'], axis=1)
+
+    # Set index
+    df = df.set_index(['feature'])
+
+    # first explode on score and feat_bound_1 (same dimensions)
+    df_explode_1 = df.explode(['score','feat_bound_1'])
+    lookup_df = df_explode_1.explode(['score','feat_bound_2'])
+
+    # sort_values but None will be in the bottom
+    lookup_df = lookup_df.sort_values(['feat_1', 'feat_2','feat_bound_1','feat_bound_2'], ascending=False)
+
+    # So reverse the ordering
+    lookup_df = lookup_df.iloc[::-1]
+
+    lookup_df['nr_features'] = 2
+
+    # being able for score arrays to be grouped
+    # for categorical features do not group by score hence we include value in score_hash
+    lookup_df['pre_hash_column'] = np.where(('categorical' in lookup_df['feat_type'])
+                                            , lookup_df['score'].astype(str) + lookup_df['feat_bound_1'].astype(str) + lookup_df['feat_bound_2'].astype(str)
+                                            , lookup_df['score'].astype(str) + lookup_df['feat_bound_1'].astype(str))
+    lookup_df['score_hash'] = pd.util.hash_pandas_object(lookup_df['pre_hash_column'], index=False)
+    hash_lookup_df = lookup_df[['score_hash','score']].drop_duplicates('score_hash').reset_index(drop=True).copy()
+
+
+    # to group by score but keep them in order (not grouping similar scores for different feature values)
+    adj_check = (lookup_df.score_hash != lookup_df.score_hash.shift()).cumsum()
+
+    lookup_df['adj_score'] = adj_check
+    lookup_df_simple = lookup_df.groupby(['feat_1', 'feat_2', 'feat_type_1', 'feat_type_2', 'score_hash', 'adj_score']).agg({'feat_bound_1':'first'
+                                                , 'feat_bound_2':'first'
+                                                ,'score_lower_bound':'first'
+                                                ,'score_upper_bound':'last'
+                                                }, sort=False).reset_index(drop=False)
+
+    lookup_df_grouped = lookup_df_simple.merge(hash_lookup_df, on='score_hash', how='left')
+
+    # ordering sets None in last place
+    lookup_df_grouped.sort_values(['feat_1', 'feat_2','feat_bound_1','feat_bound_2'], inplace=True, ascending=False)
+
+    # So reverse the ordering
+    lookup_df_grouped = lookup_df_grouped.iloc[::-1].reset_index(drop=True)
+
+    # Group second feature to list again for easier SQL writing
+    lookup_df_group_feat2 = lookup_df_grouped.groupby(['feat_1', 'feat_2', 'feat_type_1', 'feat_type_2','feat_bound_1']).agg({
+        'feat_bound_2': list
+        ,'score': list
+    }).reset_index(drop=False)
+
+    return lookup_df_group_feat2
+    
+ 
 def extractLookupTable(ebm):
     """ Provides reformatted structure of EBM
 
@@ -18,33 +148,36 @@ def extractLookupTable(ebm):
 
 
     for feature_group_index, feature_indexes in enumerate(
-            ebm.feature_groups_
+            ebm.term_features_
     ):
+
         # below includes treatment for missing values
-        model_graph = ebm.additive_terms_[feature_group_index]
+        model_graph = ebm.term_scores_[feature_group_index]
 
         # NOTE: This uses stddev. for bounds, consider issue warnings.
-        errors = ebm.term_standard_deviations_[feature_group_index]
+        errors = ebm.standard_deviations_[feature_group_index]
 
         if len(feature_indexes) == 1:
 
             # hack. remove the 0th index which is for missing values
-            model_graph = model_graph[1:]
-            errors = errors[1:]
+            model_graph = model_graph[1:-1]
+            errors = errors[1:-1]
 
-            bin_labels = ebm.preprocessor_._get_bin_labels(feature_indexes[0])
+            feature_bins = ebm.bins_[feature_group_index][0]
 
             scores = list(model_graph)
             upper_bounds = list(model_graph + errors)
             lower_bounds = list(model_graph - errors)
 
-            # Check for non categorical variable (binary or strings)
-            if (len(bin_labels) > 2) and not (all(isinstance(n, str) for n in bin_labels)):
-                feat_bound = bin_labels[1:-1] + [np.PINF]
-                feat_type = 'numeric'
-            else:
-                feat_bound = bin_labels
+            if isinstance(feature_bins, dict):
+                # categorical
+                feat_bound = list(feature_bins.keys())
                 feat_type = 'categorical'
+        
+            else:
+                feat_bound = np.append(feature_bins, np.PINF)
+                feat_type = 'numeric'
+
 
             score = scores
             score_lower_bound = lower_bounds
@@ -62,30 +195,27 @@ def extractLookupTable(ebm):
             lookup_dicts.append(lookup_dict)
 
         elif len(feature_indexes) == 2:
-            # hack. remove the 0th index which is for missing values
-            model_graph = model_graph[1:, 1:]
-            # errors = errors[1:, 1:]  # NOTE: This is commented as it's not used in this branch.
+            # hack. remove the 0th index which is for missing values and remove last value which is just zeros
+            model_graph = model_graph[1:-1, 1:-1]
 
             score = model_graph
 
-            bin_labels = [ebm.pair_preprocessor_._get_bin_labels(feature_indexes[0]),(ebm.pair_preprocessor_._get_bin_labels(feature_indexes[1]))]
-
-            feat_bound = [[]] * len(bin_labels)
-            feat_type = [[]] * len(bin_labels)
-
-            # not implemented for interaction terms
-            # score_lower_bound = [[]] * len(bin_labels)
-            # score_upper_bound = [[]] * len(bin_labels)
+            feat_bound = [[]] * len(feature_indexes)
+            feat_type = [[]] * len(feature_indexes)
 
             # Make loop over labels to update the score and features bounds
             # Check for non binary variable
-            for b in range(len(bin_labels)):
-                if (len(bin_labels[b]) > 2) and not (all(isinstance(n, str) for n in bin_labels[b])):
-                    feat_bound[b] = bin_labels[b][1:-1] + [np.PINF]
-                    feat_type[b] = 'numeric'
-                else:
-                    feat_bound[b] = bin_labels[b]
+            for b in range(len(feature_indexes)):
+
+                bin_levels = ebm.bins_[feature_indexes[b]]
+                feature_bins = bin_levels[min(len(feature_indexes), len(bin_levels)) - 1]
+            
+                if isinstance(feature_bins, dict):
+                    feat_bound[b] = list(feature_bins.keys())
                     feat_type[b] = 'categorical'
+                else:
+                    feat_bound[b] = np.append(feature_bins, np.PINF) 
+                    feat_type[b] = 'numeric'
 
             lookup_dict = {'nr_features': 2
                            , 'feature': [ebm.feature_names[feature_indexes[0]]
@@ -110,98 +240,88 @@ def extractLookupTable(ebm):
     df_of_lists_single = df_of_lists[df_of_lists['nr_features']==1].reset_index(drop=True)
     df_of_lists_double = df_of_lists[df_of_lists['nr_features']==2].reset_index(drop=True)
 
-    # All transformations below only applied on the single feature
-    # we leave the interaction (double feature) df as it is and has to be treated
-    # differently in the 2sql function
-    df_feat_types = df_of_lists_single[['feature','feat_type']].drop_duplicates().reset_index(drop=True).copy()
-
-    df_of_lists_single = df_of_lists_single.set_index(['feature'])
-    lookup_df = df_of_lists_single[['feat_bound', 'score',
-       'score_lower_bound', 'score_upper_bound']].apply(pd.Series.explode).reset_index()
-
-    # add feature type back in
-    lookup_df = lookup_df.merge(df_feat_types, on='feature', how='left')
-
-    # sort_values but None will be in the bottom
-    lookup_df.sort_values(['feature','feat_bound'], inplace=True, ascending=False)
-    # So reverse the ordering
-    lookup_df = lookup_df.iloc[::-1]
-
-    lookup_df['nr_features'] = 1
+    # Reduce Single featuers lookup
+    lookup_df_single = ReduceSingleFeature(df_of_lists_single)
 
-    # being able for score arrays to be grouped
-    # for categorical features do not group by score hence we include value in score_hash
-    lookup_df['pre_hash_column'] = np.where(lookup_df['feat_type'] == 'categorical'
-                                            , lookup_df['score'].astype(str) + lookup_df['feat_bound'].astype(str)
-                                            , lookup_df['score'].astype(str))
-    lookup_df['score_hash'] = pd.util.hash_pandas_object(lookup_df['pre_hash_column'], index=False)
-    hash_lookup_df = lookup_df[['score_hash','score']].drop_duplicates('score_hash').reset_index(drop=True).copy()
-
-    # to group by score but keep them in order (not grouping similar scores for different feature values)
-    adj_check = (lookup_df.score_hash != lookup_df.score_hash.shift()).cumsum()
-
-    lookup_df['adj_score'] = adj_check
-    lookup_df_simple = lookup_df.groupby(['feature', 'feat_type', 'score_hash', 'adj_score']).agg({'feat_bound':'first'
-                                                ,'score_lower_bound':'first'
-                                                ,'score_upper_bound':'last'
-                                                }, sort=False).reset_index(drop=False)
-
-    lookup_df_grouped = lookup_df_simple.merge(hash_lookup_df, on='score_hash', how='left')
-    # ordering sets None in last place
-    lookup_df_grouped.sort_values(['feature', 'feat_bound'], inplace=True, ascending=False)
-    # So reverse the ordering
-    lookup_df_grouped = lookup_df_grouped.iloc[::-1].reset_index(drop=True)
+    # Check if interactions exist
+    if len(df_of_lists_double) > 0:
+        # Restructure and reduce double feature df 
+        lookup_df_double = RestructureReduceInteractions(df_of_lists_double)
+    else:
+        lookup_df_double = None
 
-    return {'intercept': ebm.intercept_, 'feature_single': lookup_df_grouped, 'feature_double': df_of_lists_double}
+    return {'intercept': ebm.intercept_, 'feature_single': lookup_df_single, 'feature_double': lookup_df_double}
 
-def lookup_df_to_sql(df_dict, split=True):
+def lookup_df_to_sql(model_name, df_dict, split=True):
     # Create list for all feature score names
     feature_list = []
 
-    # start with intercept term
+    # Start with intercept term
+    intercept = df_dict['intercept'][0]
     if not split:
-        print('{intercept} -- intercept'.format(intercept=df_dict['intercept'][0]))
+        print(f'"{model_name}" AS model_name \n, ')
+        print(f'{intercept} AS intercept')
     elif split:
-        print('{intercept} AS intercept'.format(intercept=df_dict['intercept'][0]))
+        # Creating CTE to create table aliases
+        print('WITH feature_scores AS (\nSELECT')
+        print(f'"{model_name}" AS model_name \n, ')
+        print(f'{intercept} AS intercept')
 
     feature_list.append('intercept')
 
     # for single feature
-    single_features = single_features_2_sql(df_dict['feature_single'], split)
+    single_features = single_feature_handling(df_dict['feature_single'])
     feature_list = feature_list + single_features
 
     if len(df_dict['feature_double']) > 0:
-        double_features = double_features_2_sql(df_dict['feature_double'], split)
+        double_features = double_feature_sql_handling(df_dict['feature_double'])
         feature_list = feature_list + double_features
 
     if not split:
-        print('AS score')
+        # Sum up all separate scores
+        print(', ', end='')
+        print(*feature_list, sep=' + ')
+        print(' AS score')
+        
+        # Applying softmax
+        print(', EXP(score)/(EXP(score) + 1) AS probability')
+
     elif split:
-        # sum up all separate scores
+        # Add placeholder for source table
+        print('FROM <source_table> -- TODO replace with correct table')
+
+        # Close CTE and create next SELECT statement
+        print('), add_sum_scores AS (')
+        print('SELECT *')
+
+        # Sum up all separate scores
         print(', ', end='')
         print(*feature_list, sep=' + ')
         print(' AS score')
+        print('FROM feature_scores')
 
-    # Applying softmax
-    print(', EXP(score)/(EXP(score) + 1) AS probability')
+        # Close CTE and make final Select statement
+        print(')')
+        print('SELECT *')
+        # Applying softmax
+        print(', EXP(score)/(EXP(score) + 1) AS probability')
+        print('FROM add_sum_scores')
 
-def single_features_2_sql(df, split=True):
+def single_feature_handling(df):
     feature_nr = 0
     feature_list = []
     for f in df['feature'].unique():
         feature_df = df[df['feature'] == f].reset_index(drop=True)
 
-        if not split:
-            print('+\nCASE')
-        elif split:
-            print(',\nCASE')
+        # Each feature score as seperate column
+        print(',\nCASE')
 
         single_feature_2_sql(feature_df, f)
         feature_nr += 1
         feature_list.append(f'{f}_score')
 
-        if split:
-            print(f'AS {f}_score')
+        # feature score alias
+        print(f'AS {f}_score')
 
     return feature_list
 
@@ -249,109 +369,90 @@ def single_feature_2_sql(df, feature):
 
     print('END')
 
-def double_features_2_sql(df, split=True):
+def double_feature_sql_handling(df):
 
+    feature_nr = 0
     feature_list = []
+    df['double_feature'] = df['feat_1'] + '_x_' + df['feat_2']
+    for f in df['double_feature'].unique():
+        feature_df = df[df['double_feature'] == f].reset_index(drop=True)
 
-    for df_index, df_row in df.iterrows():
+        # Each double feature score as case statement
+        print(',\nCASE')
 
-        feature_name = '_x_'.join(df_row['feature'])
-        if not split:
-            print('+\nCASE')
-        elif split:
-            print(',\nCASE')
-
-        double_feature_2_sql(df_index, df_row)
-
-        feature_list.append(f'{feature_name}_score')
+        double_feature_2_sql(feature_df, f)
+        feature_nr += 1
+        feature_list.append(f'{f}_score')
 
-        if split:
-            print(f'AS {feature_name}_score')
+        # feature score alias
+        print(f'AS {f}_score')
 
     return feature_list
 
-def double_feature_2_sql(df_index, df_row):
-    first_feature = df_row['feature'][0]
-    second_feature = df_row['feature'][1]
-
-    first_feat_bound = df_row['feat_bound'][0]
-    second_feat_bound = df_row['feat_bound'][1]
-
-    first_feat_type = df_row['feat_type'][0]
-    second_feat_type = df_row['feat_type'][1]
-
-    first_feature_nbounds = len(df_row['feat_bound'][0])
-    second_feature_nbounds = len(df_row['feat_bound'][1])
-
-    scores = pd.DataFrame(df_row['score'])
-
-    # first feature
-    for f_ind in range(first_feature_nbounds):
+def double_feature_2_sql(df, double_feature):
+    for index, row in df.iterrows():
         # check if string/category
-        if first_feat_type == 'categorical':
-            print(" WHEN {feature} = '{b}' THEN \n  CASE".format(feature=first_feature
-                                                                 , b=first_feat_bound[f_ind]))
-        elif first_feat_bound[f_ind] == np.PINF:
-            print(' WHEN {feature} > {b} THEN \n   CASE'.format(feature=first_feature,
-                                                                b=first_feat_bound[f_ind - 1]
-                                                                ))
+        if row['feat_type_1'] == 'categorical':
+            print(" WHEN {feature} = '{lb}' THEN \n     CASE".format(feature=row['feat_1'], lb=row['feat_bound_1']))
+
+        # check for last numeric bound
+        elif row['feat_bound_1'] == np.PINF:
+            print(' WHEN {feature} >= {lb} THEN \n      CASE'.format(feature=row['feat_1'],
+                                                                lb=df.loc[index - 1, 'feat_bound_1']))
 
         # otherwise it should be a float
-        elif isinstance(first_feat_bound[f_ind], float):
-            # Check if first bound
-            if f_ind == 0:
-                print(' WHEN {feature} <= {b}'.format(feature=first_feature,
-                                                      b=first_feat_bound[f_ind]), end='')
-                print(' THEN \n  CASE')
+        elif isinstance(row['feat_bound_1'], float):
+            # First bound
+            if index == 0:
+                print(' WHEN {feature} <= {ub} THEN \n      CASE'.format(feature=row['feat_1'],
+                                                                    ub=row['feat_bound_1']))
+
             else:
-                # After first bound set interval
-                print(' WHEN {feature} > {lb} AND {feature} <= {ub} THEN \n     CASE'.format(feature=first_feature,
-                                                                                             lb=first_feat_bound[
-                                                                                                 f_ind - 1],
-                                                                                             ub=first_feat_bound[
-                                                                                                 f_ind]))
+                print(' WHEN {feature} > {lb} AND {feature} <= {ub} THEN \n     CASE'.format(feature=row['feat_1'],
+                                                                                            lb=
+                                                                                            df.loc[index - 1,'feat_bound_1'],
+                                                                                            ub=row['feat_bound_1']))
+
         else:
             raise "not sure what to do"
+        
+        # Nr of bound values 
+        nr_bounds = len(row['feat_bound_2'])
+
+        # Looping over the bound values for the second feature
+        for sf_index in range(nr_bounds):
 
-        # second feature
-        for s_ind in range(second_feature_nbounds):
             # check if string/category
-            if second_feat_type == 'categorical':
-                print("         WHEN {feature} = '{b}' THEN {score}".format(feature=second_feature,
-                                                                            b=second_feat_bound[s_ind],
-                                                                            score=scores.iloc[f_ind, s_ind]
-                                                                            ))
-
-            elif second_feat_bound[s_ind] == np.PINF:
-                print('         WHEN {feature} > {b} THEN {score}'.format(feature=second_feature,
-                                                                          b=second_feat_bound[s_ind - 1],
-                                                                          score=scores.iloc[f_ind, s_ind]
-                                                                          ))
+            if row['feat_type_2'] == 'categorical':
+                # check for manual imputed null values
+                print("         WHEN {feature} = '{lb}' THEN {score}".format(feature=row['feat_2'], lb=row['feat_bound_2'],
+                                                                    score=row['score'][sf_index]))
+
+            # check for last numeric bound
+            elif row['feat_bound_2'][sf_index] == np.PINF:
+                print('         WHEN {feature} >= {lb} THEN {score}'.format(feature=row['feat_2'],
+                                                                    lb=row['feat_bound_2'][sf_index - 1],
+                                                                    score=row['score'][sf_index]))
 
             # otherwise it should be a float
-            elif isinstance(second_feat_bound[s_ind], float):
-                # still use either first or second bound (based on if they are for made for null handling or not
-                if s_ind == 0:
-                    print('         WHEN {feature} <= {b} THEN {score}'.format(feature=second_feature,
-                                                                               b=second_feat_bound[s_ind],
-                                                                               score=scores.iloc[
-                                                                                   f_ind, s_ind]))
-                else:
-                    print('         WHEN {feature} > {lb} AND {feature} <= {ub} THEN {score}'.format(
-                        feature=second_feature,
-                        lb=second_feat_bound[
-                            s_ind - 1],
-                        ub=second_feat_bound[
-                            s_ind],
-                        score=scores.iloc[
-                            f_ind, s_ind]))
+            elif isinstance(row['feat_bound_2'][sf_index], float):
+                # First bound
+                if sf_index == 0:
+                    print('         WHEN {feature} <= {ub} THEN {score}'.format(feature=row['feat_2'],
+                                                                        ub=row['feat_bound_2'][sf_index],
+                                                                        score=row['score'][sf_index]))
 
+                else:
+                    print('         WHEN {feature} > {lb} AND {feature} <= {ub} THEN {score}'.format(feature=row['feat_2'],
+                                                                                            lb=row['feat_bound_2'][sf_index - 1],
+                                                                                            ub=row['feat_bound_2'][sf_index],
+                                                                                            score=row['score'][sf_index]))
 
             else:
                 raise "not sure what to do"
-
+            
             # Add ELSE 0.0 as last entry
-            if s_ind == (second_feature_nbounds - 1):
+            if sf_index + 1 == nr_bounds:
                 print("         ELSE 0.0")
 
         print('     END')
@@ -360,10 +461,15 @@ def double_feature_2_sql(df_index, df_row):
     print('END')
 
 
-def lookup_df_to_sql_multiclass(df, classes):
+def lookup_df_to_sql_multiclass(model_name, df, classes, split):
+    
+    print(f'"{model_name}" AS model_name \n, ')
+    
     class_nr = 0
+    feature_list = {}
     for c in classes:
         feature_nr = 0
+        feature_list[c] = []
         for f in df['feature'].unique():
             feature_df = df[df['feature'] == f].reset_index(drop=True)
 
@@ -372,106 +478,76 @@ def lookup_df_to_sql_multiclass(df, classes):
             elif (feature_nr == 0) & (class_nr > 0):
                 print(', CASE')
             else:
-                print('+\nCASE')
-
-            for index, row in feature_df.iterrows():
-                # check for last numeric bound
-                if row['feat_bound'] == np.PINF:
-                    print(' WHEN {feature} >= {lb} THEN {score}'.format(feature=f,
-                                                                            lb=feature_df.iloc[index - 1, :]['feat_bound'],
-                                                                            score=row['score'][class_nr]))
-                # check if string/category
-                elif isinstance(row['feat_bound'], str):
-                    print(" WHEN {feature} = '{lb}' THEN {score}".format(feature=f, lb=row['feat_bound'],
-                                                                             score=row['score'][class_nr]))
-                # check if boolean casted to int
-                elif row['feat_bound'] == int(row['feat_bound']):
-                    print(' WHEN {feature}::INT = {lb} THEN {score}'.format(feature=f, lb=row['feat_bound'],
-                                                                                score=row['score'][class_nr]))
-                # otherwise it should be a float
-                elif isinstance(row['feat_bound'], float):
-                    if index == 0:
-                        # still use either first or second bound (based on if they are for made for null handling or not
-                        if row['feat_bound'] > 0:
-                            print(' WHEN {feature} <= {ub} THEN {score}'.format(feature=f,
-                                                                                         ub=row['feat_bound'],
-                                                                                         score=row['score'][class_nr]))
-                    else:
-                        # still use either first or second bound (based on if they are for made for null handling or not
-                        if (index == 1) & (row['feat_bound'] > 0) & (feature_df.iloc[index-1, :]['feat_bound'] < 0):
-                            print(' WHEN {feature} <= {ub} THEN {score}'.format(feature=f,
-                                                                                         ub=row['feat_bound'],
-                                                                                         score=row['score'][class_nr]))
+                print(',\nCASE')
 
-                        print(' WHEN {feature} > {lb} AND {feature} <= {ub} THEN {score}'.format(feature=f,
-                                                                                                     lb=feature_df.iloc[index-1, :]['feat_bound'],
-                                                                                                     ub=row['feat_bound'],
-                                                                                                     score=row['score'][class_nr]))
-                else:
-                    raise "not sure what to do"
+            single_feature_2_sql_multiclass(feature_df, f, class_nr)
 
-            print('END')
+            # Feature score as alias
+            print(f'AS {f}_score_{c}')
+            feature_list[c].append(f'{f}_score_{c}')
 
             feature_nr += 1
-
-        print('AS score_{classification}'.format(classification=c))
-
         class_nr += 1
 
-    for c in classes:
-        if c == classes[0]:
-            print(', EXP(score_{classification})'.format(classification=c), end='')
-        else:
-            print(' + EXP(score_{classification})'.format(classification=c), end='')
-    print(' AS total_score')
-
-    for c in classes:
-        print(', EXP(score_{classification}) / (total_score) AS probability_{classification}'.format(classification=c), end='\n')
-
-def lookup_df_to_sql_multiclass_new(df, classes, split=True):
-    class_nr = 0
-    for c in classes:
-        feature_nr = 0
-        feature_list = []
-        for f in df['feature'].unique():
-            feature_df = df[df['feature'] == f].reset_index(drop=True)
+    if split:
+        # Add placeholder for source table
+        print('FROM <source_table> -- TODO replace with correct table')
+        # Close CTE and create next SELECT statement
+        print('), add_sum_scores AS (')
+        print('SELECT *')
 
-            if (feature_nr == 0) & (class_nr == 0):
-                print('CASE')
-            elif (feature_nr == 0) & (class_nr > 0):
-                print(', CASE')
-            elif not split:
-                print('+\nCASE')
-            elif split:
-                print(',\nCASE')
-
-            single_feature_2_sql_multiclass(feature_df, f, class_nr)
-
-            if split:
-                print(f'AS {f}_score_{c}')
-                feature_list.append(f'{f}_score_{c}')
-
-            feature_nr += 1
+    
+    if not split:
+        for c in classes:
+            # Sum up all separate scores
+            print(', ', end='')
+            print(*feature_list[c], sep=' + ')
+            print(f' AS score_{c}')
 
-        if not split:
-            print(f'AS score_{c}')
-        elif split:
-            # sum up all separate scores
+    elif split:
+        for c in classes:
+            # Create CTE with sums of all features scores per class
             print(', ', end='')
-            print(*feature_list, sep=' + ')
+            print(*feature_list[c], sep=' + ')
             print(f' AS score_{c}')
+        # From CTE
+        print('FROM feature_scores')
 
-        class_nr += 1
+    # Summing feature scores to total score per class
+    if not split:
+        for c in classes:
+            if c == classes[0]:
+                print(f', EXP(score_{c})', end='')
+            else:
+                print(f' + EXP(score_{c})', end='')
+        print(' AS total_score')
 
-    for c in classes:
-        if c == classes[0]:
-            print(f', EXP(score_{c})', end='')
-        else:
-            print(f' + EXP(score_{c})', end='')
-    print(' AS total_score')
+    elif split:
+        # Close CTE and create next SELECT statement
+        print('), add_sum_all_scores AS (')
+        print('SELECT *')
 
-    for c in classes:
-        print(f', EXP(score_{c}) / (total_score) AS probability_{c}', end='\n')
+        for c in classes:
+            if c == classes[0]:
+                print(f', EXP(score_{c})', end='')
+            else:
+                print(f' + EXP(score_{c})', end='')
+        print(' AS total_score')
+        
+        # Close CTE and create final SELECT statement
+        print('FROM add_sum_scores')
+        print(')')
+        print('SELECT *')
+
+        
+    # Applying softmax
+    if not split:
+        for c in classes:
+            print(f', EXP(score_{c}) / (total_score) AS probability_{c}', end='\n')
+    elif split: 
+        for c in classes:
+            print(f', EXP(score_{c}) / (total_score) AS probability_{c}', end='\n')
+        print('FROM add_sum_all_scores')
 
 def single_feature_2_sql_multiclass(df, feature, class_nr):
     for index, row in df.iterrows():
@@ -517,11 +593,11 @@ def single_feature_2_sql_multiclass(df, feature, class_nr):
 
     print('END')
 
-def ebm_to_sql(df, classes, split=True):
+def ebm_to_sql(model_name, df, classes, split=True):
     if len(classes) > 2:
-        lookup_df_to_sql_multiclass_new(df['feature_single'], classes, split=True)
+        lookup_df_to_sql_multiclass(model_name, df['feature_single'], classes, split)
     else:
-        lookup_df_to_sql(df, split)
+        lookup_df_to_sql(model_name, df, split)
 
 
 def save_model_and_extras(ebm, model_name, split, logging):
@@ -536,7 +612,7 @@ def save_model_and_extras(ebm, model_name, split, logging):
     # Write printed output to file
     with open('{model_name}/model/ebm_in_sql.sql'.format(model_name=model_name), 'w') as f:
         with redirect_stdout(f):
-            print(f"'{model_name.split('/')[1]}' AS model_name \n, ", end='')
-            ebm_to_sql(lookup_df, ebm.classes_, split)
+            model_name = model_name.split('/')[1]
+            ebm_to_sql(model_name, lookup_df, ebm.classes_, split)
     print('SQL version of EBM saved')
-    logging.info('SQL version of EBM saved')
+    logging.info('SQL version of EBM saved')
\ No newline at end of file
diff --git a/utils/pre_process.py b/scripts/utils/pre_processing/pre_process.py
similarity index 100%
rename from utils/pre_process.py
rename to scripts/utils/pre_processing/pre_process.py
diff --git a/utils/helper_functions/parsing_arguments.py b/utils/helper_functions/parsing_arguments.py
deleted file mode 100644
index d66f162..0000000
--- a/utils/helper_functions/parsing_arguments.py
+++ /dev/null
@@ -1,29 +0,0 @@
-import argparse
-
-def SetArgParser():
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--name", type=str, help="Enter project name",
-                        nargs='?', default='no_name', const='no_name')
-    parser.add_argument("--data_path", type=str, help="Enter path to csv file",
-                        nargs='?', default='no_data', const='no_data')
-    parser.add_argument("--configuration", type=str, help="Enter path to json file",
-                        nargs='?')
-    parser.add_argument("--model_name", type=str, help="Enter model type",
-                        nargs='?', default='decision_tree', const='full')
-
-    return parser
-
-def GetArgs(argv):
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--name", type=str, help="Enter project name",
-                        nargs='?', default='no_name', const='no_name')
-    parser.add_argument("--data_path", type=str, help="Enter path to csv file",
-                        nargs='?', default='no_data', const='no_data')
-    parser.add_argument("--configuration", type=str, help="Enter path to json file",
-                        nargs='?')
-    parser.add_argument("--model_name", type=str, help="Enter model type",
-                        nargs='?', default='decision_tree', const='full')
-
-    args = parser.parse_args(argv)
-
-    return args
\ No newline at end of file
diff --git a/utils/modelling/main_modeler.py b/utils/modelling/main_modeler.py
deleted file mode 100644
index 8eeb943..0000000
--- a/utils/modelling/main_modeler.py
+++ /dev/null
@@ -1,198 +0,0 @@
-import pickle
-import random
-from sklearn.model_selection import train_test_split
-from utils.modelling.performance import *
-from utils.modelling.calibration import *
-
-# The actual algorithms (grey as we refer to them dynamically)
-from utils.modelling.models import ebm
-from utils.modelling.models import decision_rule
-from utils.modelling.models import decision_tree
-
-def make_model(given_name, datasets, model_name, model_type, model_params, post_params, logging):
-    # unpack datasets
-    X_train = datasets['cv_train']['X']
-    y_train = datasets['cv_train']['y']
-    X_test = datasets['cv_test']['X']
-    y_test = datasets['cv_test']['y']
-    X_all = datasets['final_train']['X']
-    y_all = datasets['final_train']['y']
-
-
-    # check if X is a list (CV should be applied in that case)
-    if isinstance(X_train, list):
-        y_test_pred_list = list()
-        y_test_prob_list = list()
-        y_test_list = y_test
-
-        for fold_id in range(len(X_train)):
-            print('Train model on test data')
-            logging.info('Train model on test data')
-
-            # Check if model needs to be calibrated
-            if post_params['calibration'] != 'false':
-                X_slice_train, X_cal, y_slice_train, y_cal = train_test_split(X_train[fold_id],y_train[fold_id], test_size=0.2,random_state=random.randint(0,100))
-            else:
-                X_slice_train, y_slice_train = X_train[fold_id], y_train[fold_id]
-
-            clf = globals()[model_name].trainModel(X_slice_train, y_slice_train, model_params, model_type, logging)
-
-            if post_params['calibration'] != 'false':
-                clf = calibrateModel(clf, X_cal, y_cal, logging, method=post_params['calibration'], final_model=False)
-
-            # discrete predictions
-            y_test_pred_list.append(clf.predict(X_test[fold_id]))
-
-            if model_type == 'classification':
-            # probability predictions
-                # Binary classification
-                if len(clf.classes_) == 2:
-                    y_test_prob_list.append(clf.predict_proba(X_test[fold_id])[:,1])
-                elif len(clf.classes_) > 2:
-                    y_test_prob_list.append(clf.predict_proba(X_test[fold_id]))
-
-        # Merge list of prediction lists into one list
-        y_test_pred = np.concatenate(y_test_pred_list, axis=0)
-
-        if model_type == 'classification':
-            # Merge list of prediction probabilities lists into one list
-            y_test_prob = np.concatenate(y_test_prob_list, axis=0)
-
-        y_test = np.concatenate(y_test, axis=0)
-
-    # If just regular train/test split has been applied
-    else:
-        clf = globals()[model_name].trainModel(X_train, y_train, model_params, model_type, logging)
-
-        # discrete prediction
-        y_test_pred = clf.predict(X_test)
-
-        if model_type == 'classification':
-            # probability prediction
-            y_test_prob = clf.predict_proba(X_test)[:,1]
-
-    # train model one last time on all samples (upsampled)
-    print('Train final model on all data')
-    logging.info('Train final model on all data')
-
-    # Check if model needs to be calibrated
-    if post_params['calibration'] != 'false':
-        X_all, X_cal, y_all, y_cal = train_test_split(X_all, y_all, test_size=0.2, random_state=123)
-
-    clf = globals()[model_name].trainModel(X_all, y_all, model_params, model_type, logging)
-
-    # Save model in pickled format
-    filename = given_name + '/model/ebm_{model_type}.sav'.format(model_type=model_type)
-    pickle.dump(clf, open(filename, 'wb'))
-
-    # train set prediction of final model
-    y_all_pred = clf.predict(X_all)
-
-    if model_type == 'classification':
-    # probability predictions
-        # Binary classification
-        if len(clf.classes_) == 2:
-            y_all_prob = clf.predict_proba(X_all)[:,1]
-        elif len(clf.classes_) > 2:
-            y_all_prob = clf.predict_proba(X_all)
-
-
-    if post_params['calibration'] != 'false':
-        cal_clf, cal_reg = calibrateModel(clf, X_cal, y_cal, logging, method=post_params['calibration'], final_model=True)
-
-        # Save model in pickled format
-        filename = given_name + '/model/ebm_calibrated_{model_type}.sav'.format(model_type=model_type)
-        pickle.dump(cal_clf, open(filename, 'wb'))
-
-        # train set prediction of final model
-        y_all_pred = cal_clf.predict(X_all)
-
-        if model_type == 'classification':
-            # probability predictions
-            # Binary classification
-            if len(clf.classes_) == 2:
-                y_all_prob = clf.predict_proba(X_all)[:, 1]
-            elif len(clf.classes_) > 2:
-                y_all_prob = clf.predict_proba(X_all)
-
-
-    # Performance and other post modeling plots
-    if model_type == 'classification':
-        # Threshold dependant
-        plotConfusionMatrixSave(given_name, y_all, y_all_pred, data_type='final_train', logging=logging)
-        plotConfusionMatrixSave(given_name, y_test, y_test_pred, data_type='test', logging=logging)
-        classificationReportSave(given_name, y_all, y_all_pred, data_type='final_train', logging=logging)
-        classificationReportSave(given_name, y_test, y_test_pred, data_type='test', logging=logging)
-
-        if len(clf.classes_) == 2:
-            # Also create pr curve for class 0
-            y_all_neg = np.array([1 - j for j in list(y_all)])
-            y_all_prob_neg = np.array([1 - j for j in list(y_all_prob)])
-
-            y_test_list_neg = [[1 - j for j in i] for i in y_test_list]
-            y_test_prob_list_neg = [[1 - j for j in i] for i in y_test_prob_list]
-
-            # Threshold independant
-            plotClassificationCurve(given_name, y_all, y_all_prob, curve_type='roc', data_type='final_train', logging=logging)
-            plotClassificationCurve(given_name, y_test_list, y_test_prob_list, curve_type='roc', data_type='test', logging=logging)
-
-            plotClassificationCurve(given_name, y_all, y_all_prob, curve_type='pr', data_type='final_train_class1', logging=logging)
-            plotClassificationCurve(given_name, y_all_neg, y_all_prob_neg, curve_type='pr', data_type='final_train_class0', logging=logging)
-
-            plotClassificationCurve(given_name, y_test_list, y_test_prob_list, curve_type='pr', data_type='test_data_class1', logging=logging)
-            plotClassificationCurve(given_name, y_test_list_neg, y_test_prob_list_neg, curve_type='pr', data_type='test_data_class0', logging=logging)
-
-            plotCalibrationCurve(given_name, y_all, y_all_prob, data_type='final_train', logging=logging)
-            plotCalibrationCurve(given_name, y_test_list, y_test_prob_list, data_type='test', logging=logging)
-
-            plotProbabilityDistribution(given_name, y_all, y_all_prob, data_type='final_train', logging=logging)
-            plotProbabilityDistribution(given_name, y_test, y_test_prob, data_type='test', logging=logging)
-
-        # If multiclass classification
-        elif len(clf.classes_) > 2:
-            # loop through classes
-            for c in clf.classes_:
-                # creating a list of all the classes except the current class
-                other_class = [x for x in clf.classes_ if x != c]
-
-                # Get index of selected class in clf.classes_
-                class_index = list(clf.classes_).index(c)
-
-                # marking the current class as 1 and all other classes as 0
-                y_test_list_ova = [[0 if x in other_class else 1 for x in fold_] for fold_ in y_test_list]
-                y_test_prob_list_ova = [[x[class_index] for x in fold_] for fold_ in y_test_prob_list]
-
-                # concatonate probs together to one list for distribution plot
-                y_test_ova = np.concatenate(y_test_list_ova, axis=0)
-                y_test_prob_ova = np.concatenate(y_test_prob_list_ova, axis=0)
-
-                y_all_ova = [0 if x in other_class else 1 for x in y_all]
-                y_all_prob_ova = [x[class_index] for x in y_all_prob]
-
-
-                # Threshold independant
-                # plotClassificationCurve(given_name, y_all_ova, y_all_prob_ova, curve_type='roc', data_type=f'final_train_class_'{c}, logging=logging)
-                plotClassificationCurve(given_name, y_test_list_ova, y_test_prob_list_ova, curve_type='roc', data_type=f'test_class_{c}', logging=logging)
-
-                # plotClassificationCurve(given_name, y_all_ova, y_all_prob_ova, curve_type='pr', data_type='final_train_class1', logging=logging)
-                plotClassificationCurve(given_name, y_test_list_ova, y_test_prob_list_ova, curve_type='pr', data_type=f'test_class_{c}', logging=logging)
-
-                # multiClassPlotCalibrationCurvePlotly(given_name, y_all, pd.DataFrame(y_all_prob, columns=clf.classes_), title='fun')
-                plotCalibrationCurve(given_name, y_test_list_ova, y_test_prob_list_ova, data_type=f'test_class_{c}', logging=logging)
-
-                # plotProbabilityDistribution(given_name, y_all_ova, y_all_prob_ova, data_type='final_train', logging=logging)
-                plotProbabilityDistribution(given_name, y_test_ova, y_test_prob_ova, data_type=f'test_class_{c}', logging=logging)
-
-    # if regression
-    elif model_type == 'regression':
-        plotYhatVsYSave(given_name, y_test, y_test_pred, data_type='test')
-        plotYhatVsYSave(given_name, y_all, y_all_pred, data_type='final_train')
-
-        adjustedR2 = 1 - (1 - clf.score(X_all, y_all)) * (len(y_all) - 1) / (len(y_all) - X_all.shape[1] - 1)
-        print('Adjusted R2: {adjustedR2}'.format(adjustedR2=adjustedR2))
-        logging.info('Adjusted R2: {adjustedR2}'.format(adjustedR2=adjustedR2))
-
-    # Post modeling plots, specific per model but includes feature importance among others
-    globals()[model_name].postModelPlots(clf, given_name + '/feature_importance', post_params['file_type'], logging)
-
-    return clf
diff --git a/utils/modelling/models/decision_tree.py b/utils/modelling/models/decision_tree.py
deleted file mode 100644
index c96fe5a..0000000
--- a/utils/modelling/models/decision_tree.py
+++ /dev/null
@@ -1,53 +0,0 @@
-from sklearn import tree
-import pickle
-from utils.modelling.performance import *
-
-def trainModel(X_train, y_train, params, model_type, logging):
-    # if 'feature_names' not in params.keys():
-    #     params['feature_names'] = X_train.columns
-    if model_type == 'regression':
-        clf = tree.DecisionTreeRegressor(**params)
-    elif model_type == 'classification':
-        clf = tree.DecisionTreeClassifier(**params)
-    else:
-        print('Only regression or classification available')
-        logging.warning('Only regression or classification available')
-
-    clf.fit(X_train, y_train)
-    logging.info(f'Model params:\n {clf.get_params}')
-
-    print('Trained decision tree \n')
-    logging.info('Trained decision tree')
-
-    return clf
-
-def plotTreeStructureSave(clf, given_name):
-
-    plt.figure(figsize=(30,30))
-
-    tree.plot_tree(clf, fontsize=10, feature_names=clf.feature_names_in_, class_names=clf.classes_)
-    plt.savefig(f'{given_name}/tree_plot.png')
-
-    print('Tree structure plot saved')
-
-
-def featureImportanceSave(clf, given_name, file_type, logging):
-    importance_df = pd.DataFrame({'importance':clf.feature_importances_, 'feature':clf.feature_names_in_}).sort_values('importance', ascending=True).reset_index(drop=True)
-    importance_non_zero = importance_df[importance_df['importance'] > 0]
-    plotly_fig = px.bar(importance_non_zero, x='importance', y='feature')
-
-    # Update size of figure
-    plotly_fig.update_layout(xaxis_title='Importance', yaxis_title='Feature',
-                      title=f'Feature importance',
-                      width=1000,
-                      height=800)
-
-    if file_type == 'png':
-        plotly_fig.write_image(f'{given_name}/gini_feature_importance.png')
-    elif file_type == 'html':
-        plotly_fig.write_html(f'{given_name}/gini_feature_importance.html')
-
-    print('Gini feature importance plot saved')
-
-def postModelPlots(clf, given_name, file_type, logging):
-    featureImportanceSave(clf, given_name, file_type, logging)
\ No newline at end of file
diff --git a/utils/modelling/performance.py b/utils/modelling/performance.py
deleted file mode 100644
index fdb6328..0000000
--- a/utils/modelling/performance.py
+++ /dev/null
@@ -1,281 +0,0 @@
-import pandas as pd
-import numpy as np
-from matplotlib import pyplot as plt
-import plotly.express as px
-import plotly.graph_objects as go
-from sklearn.calibration import calibration_curve
-from sklearn.metrics import brier_score_loss
-
-def plotConfusionMatrixSave(given_name, y_true, y_pred, data_type, logging):
-
-    from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
-    fig, ax = plt.subplots(figsize=(10, 8))
-
-    disp = ConfusionMatrixDisplay.from_predictions(y_true, y_pred, ax=ax)
-    plt.xticks(rotation=90)
-    plt.savefig('{given_name}/performance/{data_type}_confusion_matrix.png'.format(given_name=given_name, data_type=data_type), bbox_inches='tight')
-
-    print('Confusion matrix saved')
-    logging.info('Confusion matrix saved')
-
-def classificationReportSave(given_name, y_true, y_pred, data_type, logging):
-
-    from sklearn.metrics import classification_report
-
-    report = classification_report(y_true, y_pred, output_dict=True)
-    report_df = pd.DataFrame(report).transpose()
-    report_df.to_csv('{given_name}/performance/{data_type}_classification_report.csv'.format(given_name=given_name, data_type=data_type))
-
-    print('Classification report saved')
-    logging.info('Classification report saved')
-
-def plotYhatVsYSave(given_name, y_true, y_pred, data_type, logging):
-
-    plot_df = pd.DataFrame({'y_true':y_true, 'y_pred':y_pred})
-    fig = px.scatter(plot_df, 'y_true', 'y_pred')
-    fig.write_image('{given_name}/performance/{data_type}_scatter_yhat_vs_y.png'.format(given_name=given_name, data_type=data_type))
-
-    print(f'Scatter plot of yhat vs y saved for {data_type}')
-    logging.info(f'Scatter plot of yhat vs y saved for {data_type}')
-
-def plotClassificationCurve(given_name, y_true, y_prob, curve_type, data_type, logging):
-    """
-    :param given_name: string path where to save
-    :param probs: float models predicted probability
-    :param y_true: int 1 or 0 for the y_true outcome
-    :param title: str title of the plot
-    :param curve_type: str either ROC or PR
-    :param savefig: boolean if plot should be saved
-    :param saveFileName: str where to save the plot to
-    :return:
-    """
-    from sklearn.metrics import roc_curve, precision_recall_curve, auc
-    curve_type = curve_type.lower()
-
-    if curve_type == 'roc':
-        title = 'ROC curve - {data_type} data'.format(data_type=data_type)
-        xlabel = 'False Positive Rate'
-        ylabel = 'True Positive Rate (Recall)'
-        diagCor1 = [0, 1]
-        diagCor2 = [0, 1]
-
-        curve_ = 'roc'
-
-    elif ('recall' in curve_type and 'precision' in curve_type) or curve_type == 'pr':
-        title = 'Precision-Recall curve - {data_type}'.format(data_type=data_type)
-        xlabel = 'True Positive Rate (Recall)'
-        ylabel = 'Precision'
-
-        # add the random line (= share of positives in sample)
-        if isinstance(y_true, list):
-            pos_share = sum([sum(el) for el in y_true]) / sum([len(el) for el in y_true])
-        else:
-            pos_share = sum(y_true)/len(y_true)
-
-        diagCor1 = [0, 1]
-        diagCor2 = [pos_share, pos_share]
-
-        curve_ = 'pr'
-
-    # Create plot
-    fig = go.Figure()
-
-    # set axis range to 0 - 1
-    fig.update_layout(xaxis_range=[0, 1], yaxis_range=[0, 1],
-                      xaxis_title=xlabel, yaxis_title=ylabel,
-                      title=title,
-                      width=1000,
-                      height=800)
-
-    # Add diagonal reference line
-    fig.add_shape(type="line",
-                  xref="paper", yref="paper",
-                  x0=diagCor1[0], x1=diagCor1[1], y0=diagCor2[0], y1=diagCor2[1],
-                  line=dict(
-                      color="black",
-                      width=2,
-                      dash="dot",
-                  )
-    )
-
-    # Save the auc(s)
-    auc_list = list()
-
-    if isinstance(y_prob, list):
-        for fold_id in range(len(y_prob)):
-
-            if curve_ == 'roc':
-                xVar, yVar, thresholds = roc_curve(y_true[fold_id], y_prob[fold_id])
-            elif curve_ == 'pr':
-                yVar, xVar, thresholds = precision_recall_curve(y_true[fold_id], y_prob[fold_id])
-
-            # Calculate area under curve (AUC)
-            auc_list.append(auc(xVar, yVar))
-            fig.add_trace(go.Scatter(x=xVar, y=yVar, mode='lines', name=f'Fold {fold_id}'))
-    else:
-
-        if curve_ == 'roc':
-            xVar, yVar, thresholds = roc_curve(y_true, y_prob)
-        elif curve_ == 'pr':
-            yVar, xVar, thresholds = precision_recall_curve(y_true, y_prob)
-
-        # Calculate area under curve (AUC)
-        auc_list.append(auc(xVar, yVar))
-        fig.add_trace(go.Scatter(x=xVar, y=yVar, mode='lines', name=f'Fitted model'))
-
-    # add (average) auc in image
-    fig.add_annotation(x=0.5, y=0,
-                       text=f"Mean AUC: {np.mean(auc_list)}",
-                       showarrow=False,
-                       yshift=10)
-
-    fig.write_image(f'{given_name}/performance/{data_type}_{curve_type}_plot.png')
-
-    print(f'Created and saved {curve_type} plot for {data_type} data')
-    logging.info(f'Created and saved {curve_type} plot for {data_type} data')
-
-    return np.mean(auc_list)
-
-def plotCalibrationCurve(given_name, y_true, y_prob, data_type, logging):
-    """
-    Plot the calibration curve for a set of true and predicted values
-
-    :param actuals: true target value
-    :param probs: predicted probability of target
-    :param bins: how many bins to divide data in for plotting
-    :param savefig: boolean if plot should be saved
-    :param saveFileName: str path to which to save the plot
-    :return: calibration plot
-    """
-
-    # Create plot
-    fig = go.Figure()
-
-    # set axis range to 0 - 1
-    fig.update_layout(xaxis_range=[0, 1], yaxis_range=[0, 1],
-                      xaxis_title='Predicted probability', yaxis_title='Share being positive',
-                      title=f'Calibration plot (reliability curve) - {data_type} data',
-                      width=1000,
-                      height=800)
-
-    # Add diagonal reference line
-    fig.add_shape(type="line",
-                  xref="paper", yref="paper",
-                  x0=0, x1=1, y0=0, y1=1,
-                  line=dict(
-                      color="black",
-                      width=2,
-                      dash="dot",
-                  )
-    )
-
-    # Save the brier score loss
-    bsl_list = list()
-
-    if isinstance(y_prob, list):
-        for fold_id in range(len(y_prob)):
-            # summaries actuals and predicted probs to (bins) number of points
-            fraction_of_positives, mean_predicted_value = \
-                calibration_curve(y_true[fold_id], y_prob[fold_id], n_bins=10, strategy='quantile')
-
-            # Calculate area under curve (AUC)
-            bsl_list.append(brier_score_loss(y_true[fold_id], y_prob[fold_id]))
-            fig.add_trace(go.Scatter(x=mean_predicted_value, y=fraction_of_positives, mode='markers+lines', name=f'Fold {fold_id}'))
-    else:
-        # summaries actuals and predicted probs to (bins) number of points
-        fraction_of_positives, mean_predicted_value = \
-            calibration_curve(y_true, y_prob, n_bins=10, strategy='quantile')
-
-        # Calculate area under curve (AUC)
-        bsl_list.append(brier_score_loss(y_true, y_prob))
-        fig.add_trace(go.Scatter(x=mean_predicted_value, y=fraction_of_positives, mode='markers+lines', name=f'Fitted model'))
-
-    # add (average) auc in image
-    fig.add_annotation(x=0.5, y=0,
-                       text=f"Mean Brier Score Loss: {np.mean(bsl_list)}",
-                       showarrow=False,
-                       yshift=10)
-
-    fig.write_image(f'{given_name}/performance/{data_type}_calibration_plot.png')
-
-    print(f'Created and saved calibration plot')
-    logging.info(f'Created and saved calibration plot')
-
-# for multiclass classification WIP
-def multiClassPlotCalibrationCurvePlotly(given_name, actuals, probs, title, bins=10):
-    """
-    Plot the calibration curve for a set of true and predicted values
-
-    :param actuals: true target value
-    :param probs: predicted probabilities per class as pandas dataframe
-    :param bins: how many bins to divide data in for plotting
-    :param savefig: boolean if plot should be saved
-    :param saveFileName: str path to which to save the plot
-    :return: calibration plot
-    """
-    calibration_df = pd.DataFrame(columns=['classification','fraction_of_positives','mean_predicted_value'])
-
-    # loop through different classes
-    for cl in probs.columns:
-        y_true = np.where(actuals == cl, 1, 0)
-        y_prob = probs[cl]
-
-        # summarise actuals and predicted probs to (bins) number of points
-        fraction_of_positives, mean_predicted_value = \
-            calibration_curve(y_true, y_prob, n_bins=bins, strategy='quantile')
-
-        calibration_df = calibration_df.append(pd.DataFrame({'classification':cl, 'fraction_of_positives':fraction_of_positives, 'mean_predicted_value':mean_predicted_value}), ignore_index=True)
-
-    # Create scatter plot
-    fig = px.scatter(calibration_df, x='mean_predicted_value', y='fraction_of_positives', title='Calibration plot (reliability curve)', color='classification',width=1200, height=800)
-
-
-    # Make trace be line plus dots
-    for l in range(len(fig.data)):
-        fig.data[l].update(mode='markers+lines')
-
-    # set axis range to 0 - 1
-    # fig.update_layout(xaxis_range=[-0.1,1.1], yaxis_range=[-0.1,1.1], xaxis_title='Predicted probability', yaxis_title='Fraction of positives')
-
-    # Add diagonal reference line
-    fig.add_shape(type="line",
-                  xref="paper", yref="paper",
-                  x0=0, x1=1, y0=0, y1=1,
-                  line=dict(
-                      color="black",
-                      width=2,
-                      dash="dot",
-                  )
-    )
-
-    fig.show(renderer='browser')
-
-def plotProbabilityDistribution(given_name, y_true, y_prob, data_type, logging):
-    import plotly.figure_factory as ff
-    df = pd.DataFrame({'actuals': y_true, 'prob': y_prob})
-    do = df[df['actuals']==1]['prob']
-    dont = df[df['actuals']==0]['prob']
-
-    # Catching any kind of exception
-    try:
-        # Create distplot with custom bin_size
-        fig = ff.create_distplot([do,dont], ['1','0'], colors=['green','red'], bin_size=.01)
-    except Exception as e:
-        print(f'Could not create distribution plot because of \n{e}')
-        logging.info(f'Could not create distribution plot because of \n{e}')
-        return
-
-    # Update size of figure
-    fig.update_layout(xaxis_title='Predicted probability', yaxis_title='Frequency',
-                      title=f'Distribution plot - {data_type} data',
-                      width=1000,
-                      height=800)
-
-    # fig.show(renderer='browser')
-
-    fig.write_image(f'{given_name}/performance/{data_type}_distribution_plot.png')
-
-    print(f'Created and saved probability distribution plot')
-    logging.info(f'Created and saved probability distribution plot')
-
-    return


Test patch:
diff --git a/scripts/modeltester.py b/scripts/modeltester.py
new file mode 100644
index 0000000..44c2635
--- /dev/null
+++ b/scripts/modeltester.py
@@ -0,0 +1,128 @@
+import os 
+import joblib 
+import pandas as pd
+import logging
+
+from interpret.glassbox import ExplainableBoostingClassifier, ExplainableBoostingRegressor
+from interpret.glassbox import DecisionListClassifier
+from interpret.glassbox import LinearRegression, LogisticRegression
+
+from utils.helper_functions.parsing_arguments import GetArgs
+from utils.modelling.performance import *
+from sklearn import tree
+
+def apply_model(args):
+
+    # Set destination
+    destination = args.destination_path
+
+    # Set Logger
+    logging.basicConfig(format='%(asctime)s %(message)s', filename=destination+'/logging.log', level=logging.DEBUG)
+    logging.getLogger('matplotlib.font_manager').disabled = True # ignore matplotlibs font warnings
+
+    # Load in model
+    model = joblib.load(open(args.model_path, 'rb'))
+
+    # Load in data
+    df = pd.read_csv(args.data_path)
+    
+    target_col = model.target
+    feature_cols = [ x for x in model.feature_names if " x " not in x ]
+
+    # Perform inference
+    y_true = df[target_col]
+    X = df[feature_cols]
+
+    y_pred = model.predict(X)
+
+    # Create performance metrics graphs
+    if hasattr(model, 'classes_'):
+
+        # If classification make probability prediction
+        y_prob = model.predict_proba(X)
+
+        if len(model.classes_) == 2:
+
+            # for binary only take class one probabilities
+            y_prob = y_prob[:,1]
+
+            # Threshold dependant
+            plotConfusionMatrix(destination, y_true, y_prob, y_pred, file_type='html', data_type='test', logging=logging)
+
+            # Also create pr curve for class 0
+            y_neg = np.array([1 - j for j in list(y_true)])
+            y_prob_neg = np.array([1 - j for j in list(y_prob)])
+
+            # Threshold independant
+            plotClassificationCurve(destination, y_true, y_prob, curve_type='roc', data_type='test', logging=logging)
+
+            plotClassificationCurve(destination, y_true, y_prob, curve_type='pr', data_type='test_class1', logging=logging)
+            plotClassificationCurve(destination, y_neg, y_prob_neg, curve_type='pr', data_type='test_class0', logging=logging)
+
+            plotClassificationCurve(destination, y_true, y_prob, curve_type='pr', data_type='test_data_class1', logging=logging)
+            plotClassificationCurve(destination, y_neg, y_prob_neg, curve_type='pr', data_type='test_data_class0', logging=logging)
+
+            plotCalibrationCurve(destination, y_true, y_prob, data_type='test', logging=logging)
+
+            plotProbabilityDistribution(destination, y_true, y_prob, data_type='test', logging=logging)
+
+
+        # If multiclass classification
+        elif len(model.classes_) > 2:
+
+            # loop through classes
+            for c in model.classes_:
+                # creating a list of all the classes except the current class
+                other_class = [x for x in model.classes_ if x != c]
+
+                # Get index of selected class in model.classes_
+                class_index = list(model.classes_).index(c)
+
+                # marking the current class as 1 and all other classes as 0
+                y_ova = [0 if x in other_class else 1 for x in y]
+                y_prob_ova = [x[class_index] for x in y_prob]  
+
+                # Threshold independant
+                plotClassificationCurve(destination, y_ova, y_prob_ova, curve_type='roc', data_type=f'test_class_{c}', logging=logging)
+
+                plotClassificationCurve(destination, y_ova, y_prob_ova, curve_type='pr', data_type=f'test_class_{c}', logging=logging)
+
+                plotCalibrationCurve(destination, y_ova, y_prob_ova, data_type=f'test_class_{c}', logging=logging)
+
+                plotProbabilityDistribution(destination, y_ova, y_prob_ova, data_type=f'test_class_{c}', logging=logging)
+
+    # if regression
+    else:
+        plotYhatVsYSave(destination, y_true, y_pred, data_type='test', logging=logging)
+
+        plotQuantileError(destination, y_true, y_pred, data_type='test', logging=logging)
+
+        regressionMetricsTable(destination, y_true, y_pred, X, data_type='test', logging=logging)
+
+
+if __name__ == '__main__':
+
+    set_env = 'prod' # either prod or dev
+
+    # Check if this script is run from terminal
+    if set_env == 'prod':
+        # (Prod) script is being run through the terminal
+        argvals = None
+    else:
+        # (Dev) script is not being run through the terminal
+
+        # Command line arguments used for testing
+        argvals = '--model_path ../trained_models/20230622_datathon_ebm/model/ebm_regression.sav ' \
+                  '--data_path ../input/data/20230525_cq_datathon.csv ' \
+                  '--destination_path ../trained_models/20230622_datathon_ebm/tested_datasets/20230525_cq_datathon.csv'.split() # example of passing test params to parser
+
+        # settings
+        pd.set_option('display.max_rows', 500)
+        pd.set_option('display.max_columns', 10)
+
+
+    # Get arguments from the CLI
+    args = GetArgs('test_model', argvals)
+
+    # Run main with given arguments
+    apply_model(args)
\ No newline at end of file
diff --git a/test_model.py b/test_model.py
new file mode 100644
index 0000000..c268538
--- /dev/null
+++ b/test_model.py
@@ -0,0 +1,88 @@
+import os 
+import re
+import sys
+
+# Select data
+# List files in input/data/ directory
+data_dir = "input/data/"
+files = os.listdir(data_dir)
+files.sort()
+
+print("Files in input/data/:")
+for i, file in enumerate(files, 1):
+    print(f"{i}. {file}")
+
+# Ask for CSVPATH
+csv_path = None
+while csv_path is None:
+    csv_file_index = input("\nSelect CSV file for testing the model: ")
+    try:
+        csv_file_index = int(csv_file_index) - 1
+        csv_path = os.path.join(data_dir, files[csv_file_index])
+    except (ValueError, IndexError):
+        print("Invalid option, try again.")
+
+print(f"CSV file {csv_path} will be used for testing model")
+
+# Select model
+# List files in trained_models/ directory
+model_dir = "trained_models/"
+files = []
+
+for root, dirs, filenames in os.walk(model_dir):
+    for dir in dirs:
+        subdir = os.path.join(root, dir)
+        for filename in os.listdir(subdir):
+            if filename.endswith(".sav"):
+                files.append(os.path.join(subdir, filename))
+
+files.sort()
+
+print(f"\nModels in {model_dir}:")
+for i, file in enumerate(files, 1):
+    print(f"{i}. {re.sub(model_dir, '', file)}")
+
+# Ask for ModelPath
+model_path = None
+while model_path is None:
+    model_file_index = input("\nSelect model to apply on test data: ")
+    try:
+        model_file_index = int(model_file_index) - 1
+        model_path = os.path.join(model_dir, files[model_file_index])
+        model_path = files[model_file_index]
+    except (ValueError, IndexError):
+        print("Invalid option, try again.")
+
+print(f"Model {model_path} will be used for testing")
+
+# Create extra folder in trained model folder
+# Make directory with current data and model name
+try:
+    model_folder = model_path.split('trained_models/')[1].split('/')[0]
+    csv_name = csv_path.split('/')[-1]
+    destination_path = f"trained_models/{model_folder}/tested_datasets/{csv_name}"
+
+    if not os.path.exists(destination_path.split(csv_name)[0]):
+        os.makedirs(destination_path.split(csv_name)[0])
+        
+    os.makedirs(destination_path)
+    os.makedirs(f"{destination_path}/performance")
+
+except FileExistsError:
+    sys.exit("Error: tested_datasets directory already exists within trained_models/{model_folder}")
+
+# Trigger python script to do the actual work
+print("\nRunning command:")
+if sys.platform == 'win32':
+    command = f".ml2sql\Scripts\python.exe scripts/modeltester.py --model_path {model_path} --data_path {csv_path} --destination_path {destination_path}"
+else:
+    command = f".ml2sql/bin/python scripts/modeltester.py --model_path {model_path} --data_path {csv_path} --destination_path {destination_path}"
+
+print(command)
+
+os.system(command)
+
+print("\nModel performance outputs can be found in folder:")
+print(f"{os.getcwd()}/trained_models/{model_folder}/tested_datasets/{csv_name}/")
+
+