Repository: AI-SDC/SACRO-ML (language: python)

Issue:
10 add worst case attack code
Closes #10 

Adds code to perform a worst case attack, as well as an example.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..21104add
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,142 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+
+# Saved pickle files
+*.pkl
+
+# Saved joblib files
+*.sav
+
+# Data folder
+/data
+
+# data files
+*.csv
\ No newline at end of file
diff --git a/Contributing.md b/CONTRIBUTING.md
similarity index 100%
rename from Contributing.md
rename to CONTRIBUTING.md
diff --git a/attacks/__init__.py b/attacks/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/attacks/metrics.py b/attacks/metrics.py
new file mode 100644
index 00000000..0cf19535
--- /dev/null
+++ b/attacks/metrics.py
@@ -0,0 +1,220 @@
+'''
+Calculate metrics.
+'''
+
+from typing import Iterable#, Optional, Any
+import numpy as np
+from sklearn.metrics import confusion_matrix
+from sklearn.metrics import roc_auc_score
+from sklearn.metrics import roc_curve
+from scipy.stats import norm
+from scipy import interpolate
+from .mia_extremecase import min_max_disc
+
+# pylint: disable = invalid-name
+
+VAR_THRESH = 1e-2
+
+def div(x: float, y: float, default: float) -> float:
+    """Solve the problem of division by 0 and round up.
+    If y is non-zero, perform x/y and round to 8dp. If it is zero, return the default
+
+    Parameters
+    ----------
+    x: float
+        numerator
+    y: float
+        denominator
+    default: float
+        return value if y == 0
+
+    Returns
+    -------
+        division: float
+            x / y, or default if y == 0
+    """
+    if y != 0:
+        division = round(float(x / y), 8)
+    else:
+        division = float(default)
+    return division
+
+def tpr_at_fpr(
+    y_true: Iterable[float],
+    y_score: Iterable[float],
+    fpr: float=0.001,
+    fpr_perc: bool=False) -> float:
+    """Compute the TPR at a fixed FPR.
+    In particular, returns the TPR value corresponding to a particular FPR. Uses interpolation
+    to fill in gaps.
+
+    Parameters
+    ----------
+    y_true: Iterable[float]
+        actual class labels
+    y_score: Iterable[float]
+        predicted score
+    fpr: float
+        false positive rate at which to compute true positive rate
+    fpr_perc: bool
+        if the fpr is defined as a percentage
+
+    Returns
+    -------
+    tpr: float
+        true positive rate at fpr
+    """
+
+    if fpr_perc:
+        fpr /= 100.
+
+
+    fpr_vals, tpr_vals, thresh_vals = roc_curve(y_true, y_score)
+    thresh_from_fpr = interpolate.interp1d(fpr_vals, thresh_vals)
+    tpr_from_thresh = interpolate.interp1d(thresh_vals, tpr_vals)
+
+    thresh = thresh_from_fpr(fpr)
+    tpr = tpr_from_thresh(thresh)
+
+    return tpr
+
+def expected_auc_var(auc: float, num_pos: int, num_neg: int) -> float:
+    """"Compute variance of AUC under assumption of uniform probabilities
+    uses the expression given as eqn (2) in  https://cs.nyu.edu/~mohri/pub/area.pdf
+
+    Parameters
+    ----------
+
+    auc: float
+        auc value at which to compute the variance
+    num_pos: int
+        number of positive examples
+    num_neg: int
+        number of negative examples
+
+    Returns
+    -------
+    var: float
+        null variance of AUC
+    """
+    p_xxy = p_xyy = 1/3
+    var = (auc * (1 - auc) + (num_pos - 1) * (p_xxy - auc**2) + (num_neg - 1) * (p_xyy - auc**2)) /\
+        (num_pos * num_neg)
+    return var
+
+def get_metrics(clf, # pylint: disable = too-many-locals
+                X_test:np.ndarray,
+                y_test:np.ndarray,
+                permute_rows:bool=True):
+    """
+    Calculate metrics, including attacker advantage for MIA binary.
+    Implemented as Definition 4 on https://arxiv.org/pdf/1709.01604.pdf
+    which is also implemented in tensorFlow-privacy https://github.com/tensorflow/privacy
+
+    Parameters
+    ----------
+    clf: sklearn.Model
+        trained model
+    X_test: np.ndarray
+        test data matrix
+    y_test: np.ndarray
+        test data labels
+
+    Returns
+    -------
+    metrics: dict
+        dictionary of metric values
+
+    Notes
+    -----
+    Includes the following metrics
+
+    True positive rate or recall (TPR)
+    False positive rate (FPR), proportion of negative examples incorrectly classified as positives
+    False alarm rate (FAR), proportion of objects classified as positives that are incorrect,
+        also known as false discovery rate
+    True neagative rate (TNR)
+    Positive predictive value or precision (PPV)
+    Negative predictive value (NPV)
+    False neagative rate (FNR)
+    Accuracy (ACC)
+    F1 Score - harmonic mean of precision and recall.
+    Advantage
+    """
+    metrics = {}
+    if permute_rows:
+        N, _ = X_test.shape
+        order = np.random.RandomState(seed=10).permutation(N) # pylint: disable = no-member
+        X_test = X_test[order, :]
+        y_test = y_test[order]
+    y_pred = clf.predict(X_test)
+    tn, fp, fn, tp = confusion_matrix(y_test, y_pred).ravel()
+    #print('tn', tn, 'fp',fp,'fn', fn,'tp', tp)
+
+    # true positive rate or recall
+    metrics['TPR'] = round(float(tp/(tp + fn)), 8)
+    # false positive rate, proportion of negative examples incorrectly classified as positives
+    metrics['FPR'] = round(float(fp / (fp + tn)), 8)
+    # False alarm rate, proportion of things classified as positives that are incorrect,
+    # also known as false discovery rate
+    metrics['FAR'] = div(fp, (fp + tp), 0)
+    # true negative rate or specificity
+    metrics['TNR'] = round(float(tn / (tn + fp)), 8)
+    # precision or positive predictive value
+    metrics['PPV'] = div(tp, (tp + fp), 0)
+    # negative predictive value
+    metrics['NPV'] = div(tn, (tn + fn), 0)
+    # false negative rate
+    metrics['FNR'] = round(float(fn / (tp + fn)), 8)
+    # overall accuracy
+    metrics['ACC'] = round(float((tp + tn) / (tp + fp + fn + tn)), 8)
+    # harmonic mean of precision and sensitivity
+    metrics['F1score'] = div(2*metrics['PPV']*metrics['TPR'], metrics['PPV']+metrics['TPR'], 0)
+    # Advantage: TPR - FPR
+    metrics['Advantage'] = float(abs(metrics['TPR']-metrics['FPR']))
+
+    #calculate AUC of model
+    y_pred_proba = clf.predict_proba(X_test)[:, 1]
+    metrics['AUC'] = round(roc_auc_score(y_test, y_pred_proba),8)
+
+    # Calculate AUC p-val
+    auc_null_var = expected_auc_var(0.5, y_test.sum(), len(y_test) - y_test.sum())
+    # Probability of getting an AUC higher than that observed
+    prob_higher_auc = 1 - norm.cdf(metrics['AUC'], loc=0.5, scale=np.sqrt(auc_null_var))
+    metrics['P_HIGHER_AUC'] = prob_higher_auc
+
+    fmax, fmin, fdif, pdif = min_max_disc(y_test, y_pred_proba)
+    metrics['FMAX01'] = fmax
+    metrics['FMIN01'] = fmin
+    metrics['FDIF01'] = fdif
+    metrics['PDIF01'] = -pdif # use -log(p) so answer is positive
+
+    fmax, fmin, fdif, pdif = min_max_disc(y_test, y_pred_proba, x_prop=0.2)
+    metrics['FMAX02'] = fmax
+    metrics['FMIN02'] = fmin
+    metrics['FDIF02'] = fdif
+    metrics['PDIF02'] = -pdif # use -log(p) so answer is positive
+
+    fmax, fmin, fdif, pdif = min_max_disc(y_test, y_pred_proba, x_prop=0.01)
+    metrics['FMAX001'] = fmax
+    metrics['FMIN001'] = fmin
+    metrics['FDIF001'] = fdif
+    metrics['PDIF001'] = -pdif # use -log(p) so answer is positive
+
+
+    # Add some things useful for debugging / filtering
+    metrics['pred_prob_var'] = y_pred_proba.var()
+
+    # TPR at various FPR
+    fpr_vals = [0.5, 0.2, 0.1, 0.01, 0.001, 0.00001]
+    for fpr in fpr_vals:
+        tpr = tpr_at_fpr(y_test, y_pred_proba, fpr=fpr)
+        name = f'TPR@{fpr}'
+        metrics[name] = tpr
+
+    fpr, tpr, roc_thresh = roc_curve(y_test, y_pred_proba)
+    metrics['fpr'] = fpr
+    metrics['tpr'] = tpr
+    metrics['roc_thresh'] = roc_thresh
+
+    return metrics
diff --git a/attacks/mia_extremecase.py b/attacks/mia_extremecase.py
new file mode 100644
index 00000000..bcbae0c4
--- /dev/null
+++ b/attacks/mia_extremecase.py
@@ -0,0 +1,84 @@
+"""Method for computing extreme case metrics"""
+from typing import Tuple
+import numpy as np
+from scipy.stats import norm
+
+def min_max_disc(
+    y_true: np.ndarray,
+    pred_probs: np.ndarray,
+    x_prop: float=0.1,
+    log_p: bool=True) -> Tuple[float, float, float, float]: #pylint: disable = line-too-long
+    """
+    Non-average-case methods for MIA attacks. Considers actual frequency of membership
+    amongst samples with highest- and lowest- assessed probability of membership. If an
+    MIA method confidently asserts that 5% of samples are members and 5% of samples are
+    not, but cannot tell for the remaining 90% of samples, then these metrics will flag
+    this behaviour, but AUC/advantage may not. Since the difference may be noisy, a
+    p-value against a null of independence of true membership and assessed membership
+    probability (that is, membership probabilities are essentially random) is also used
+    as a metric (using a usual Gaussian approximation to binomial). If the p-value is
+    low and the frequency difference is high (>0.5) then the MIA attack is successful
+    for some samples.
+
+    Parameters
+    ----------
+        y: np.ndarray
+            true labels
+        yp: np.ndarray
+            probabilities of labels, or monotonic transformation of probabilties
+        xprop: float
+            proportion of samples with highest- and lowest- probability of membership to be
+            considered
+        logp: bool
+            convert p-values to log(p).
+
+    Returns
+    -------
+        maxd: float
+            frequency of y=1 amongst proportion xprop of individuals with highest assessed
+            membership probability
+        mind: float
+            frequency of y=1 amongst proportion xprop of individuals with lowest assessed
+            membership probability
+        mmd: float
+            difference between maxd and mind
+        pval: float
+            p-value or log-p value corresponding to mmd against null hypothesis that random
+            variables corresponding to y and yp are independent.
+
+    Notes
+    -----
+
+    Examples
+    --------
+    >>> y = np.random.choice(2, 100)
+    >>> yp = np.random.rand(100)
+    >>> maxd, mind, mmd, pval = min_max_desc(y, yp, xprop=0.2, logp=True)
+
+    """
+
+    n_examples = int(np.ceil(len(y_true) * x_prop))
+    pos_frequency = np.mean(y_true) # average frequency
+    y_order = np.argsort(pred_probs) # ordering permutation
+
+    # Frequencies
+    # y values corresponding to lowest k values of yp
+    y_lowest_n = y_true[y_order[:n_examples]]
+    # y values corresponding to highest k values of yp
+    y_highest_n = y_true[y_order[-(n_examples):]]
+    maxd = np.mean(y_highest_n)
+    mind = np.mean(y_lowest_n)
+    mmd = maxd - mind
+
+    # P-value
+    # mmd is asymptotically distributed as N(0,sdm^2) under null.
+    sdm = np.sqrt(2 * pos_frequency * (1 - pos_frequency) / n_examples)
+    pval = 1 - norm.cdf(mmd, loc=0, scale=sdm) # normal CDF
+    if log_p:
+        if pval< 1e-50:
+            pval=-115.13
+        else:
+            pval=np.log(pval)
+
+    # Return
+    return maxd, mind, mmd, pval
diff --git a/attacks/report.py b/attacks/report.py
new file mode 100644
index 00000000..4e7fa089
--- /dev/null
+++ b/attacks/report.py
@@ -0,0 +1,304 @@
+"""Code for automatic report generation"""
+import json
+from typing import List, Dict
+import pylab as plt
+import numpy as np
+from fpdf import FPDF
+
+# Adds a border to all pdf cells of set to 1 -- useful for debugging
+BORDER = 0
+
+DISPLAY_METRICS = [
+    "AUC",
+    "ACC",
+    "Advantage",
+    "FDIF01",
+    "PDIF01",
+    "TPR@0.1",
+    "TPR@0.01",
+    "TPR@0.001",
+    "TPR@1e-05"
+]
+
+MAPPINGS = {
+    'PDIF': lambda x: np.exp(-x)
+}
+
+INTRODUCTION = (
+    'This report provides a summary of a series of simulated attack experiments performed '
+    'on the model outputs provided. An attack model is trained to attempt to distinguish '
+    'between outputs from training (in-sample) and testing (out-of-sample) data. The metrics '
+    'below describe the success of this classifier. A successful classifier indicates that the '
+    'original model is unsafe and should not be allowed to be relased from the TRE.\n'
+    'In particular, the simulation splits the data provided into test and train sets (each will '
+    'in- and out-of-sample examples). The classifier is trained on the train set and evaluated '
+    'on the test set. This is repeated with different train/test splits a user-specified number '
+    'of times.\n'
+    'To help place the results in context, the code may also have run a series of baseline '
+    'experiments. In these, random model outputs for hypothetical in- and out-of-sample data are '
+    'generated with identical statistical properties. In these baseline cases, there is no signal '
+    'that an attacker could leverage and therefore these values provide a baseline against '
+    'which the actual values can be compared.\n'
+    'For some metrics (FDIF and AUC), we are able to compute p-values. In each case, shown below '
+    '(in the Global metrics sections) is the number of repetitions that exceeded the p-value '
+    'threshold both without, and with correction for multiple testing (Benjamini-Hochberg '
+    'procedure).\n'
+    'ROC curves for all real (red) and dummy (blue) repetitions are provided. These are shown in '
+    'log space (as reommended here [ADD URL]) to emphasise the region in which risk is highest -- '
+    'the bottom left (are high true positive rates possible with low false positive rates).\n'
+    'A description of the metrics and how to interpret them within the context of an attack is '
+    'given below.'
+)
+
+LOGROC_CAPTION = (
+    'This plot shows the False Positive Rate (x) versus the True Positive Rate (y). '
+    'The axes are in log space enabling us to focus on areas where the False Positive Rate is low '
+    '(left hand area). Curves above the y = x line (black dashes) in this region represent a '
+    'disclosure risk as an attacker can obtain many more true than false positives. '
+    'The solid coloured lines show the curves for the attack simulations with the true model '
+    'outputs. The lighter grey lines show the curves for randomly generated outputs with '
+    'no structure (i.e. in- and out-of- sample predictions are generated from the same '
+    'distributions. Solid curves consistently higher than the grey curves in the left hand '
+    'part of the plot are a sign of concern.'
+)
+
+GLOSSARY = {
+    'AUC': 'Area Under the ROC curve',
+    'True Positive Rate (TPR)': (
+        'The true positive rate is the number of True Positives that are predicted as positive as '
+        'a proportion of the total number of positives. If an attacker has N examples that were '
+        'actually in the training set, the TPR is the proportion of these that they predict as '
+        'being in the training set.'
+    ),
+    'ACC': 'The proportion of predictions that the attacker makes that are correct.'
+}
+
+class NumpyArrayEncoder(json.JSONEncoder):
+    """Json encoder that can cope with numpy arrays"""
+    def default(self, o):
+        """If an object is an np.ndarray, convert to list"""
+        if isinstance(o, np.ndarray):
+            return o.tolist()
+        return json.JSONEncoder.default(self, o)
+
+def _write_dict(pdf, input_dict, indent=0, border=BORDER):
+    """Write a dictionary to the pdf"""
+    for key, value in input_dict.items():
+        pdf.set_font('arial', 'B', 14)
+        pdf.cell(75, 5, key, border, 1, 'L')
+        pdf.cell(indent, 0)
+        pdf.set_font('arial', '', 12)
+        pdf.multi_cell(150, 5, value, border, 'L')
+        pdf.ln(h=5)
+
+def _title(pdf, text, border=BORDER, font_size=24, font_style='B'):
+    """Write a title block"""
+    pdf.set_font('arial', font_style, font_size)
+    pdf.ln(h=5)
+    pdf.cell(0, 0, text, border, 1, 'C')
+    pdf.ln(h=5)
+
+def _subtitle(pdf, text, indent=10, border=BORDER, font_size=12, font_style='B'): # pylint: disable = too-many-arguments
+    """Write a subtitle block"""
+    pdf.cell(indent, border=border)
+    pdf.set_font('arial', font_style, font_size)
+    pdf.cell(75, 10, text, border, 1)
+
+def _line(pdf, text, indent=0, border=BORDER, font_size=11, font_style='', font='arial'): # pylint: disable = too-many-arguments
+    """Write a standard block"""
+    if indent > 0:
+        pdf.cell(indent, border=border)
+    pdf.set_font(font, font_style, font_size)
+    pdf.multi_cell(0, 5, text, border, 1)
+
+def _roc_plot_single(metrics, save_name):
+    """Create a roc_plot for a single experiment"""
+    plt.figure()
+    plt.plot([0, 1], [0, 1], 'k--')
+    plt.plot(metrics['fpr'], metrics['tpr'], 'r', linewidth=2)
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.grid()
+    plt.xlabel("True Positive Rate")
+    plt.ylabel("False Positive Rate")
+    plt.tight_layout()
+    plt.savefig(save_name)
+
+def _roc_plot(metrics, dummy_metrics, save_name):
+    """Create a roc plot for multiple repetitions"""
+    plt.figure()
+    plt.plot([0, 1], [0, 1], 'k--')
+    if dummy_metrics is None or len(dummy_metrics) == 0:
+        do_dummy = False
+    else:
+        do_dummy = True
+
+    # Compute average ROC
+    base_fpr = np.linspace(0, 1, 1000)
+    # base_fpr = np.logspace(-4, 0, 1000)
+    all_tpr = np.zeros((len(metrics), len(base_fpr)), float)
+    for i, metric_set in enumerate(metrics):
+        all_tpr[i, :] = np.interp(base_fpr, metric_set['fpr'], metric_set['tpr'])
+
+    if do_dummy:
+        all_tpr_dummy = np.zeros((len(dummy_metrics), len(base_fpr)), float)
+        for i, metric_set in enumerate(dummy_metrics):
+            all_tpr_dummy[i, :] = np.interp(base_fpr, metric_set['fpr'], metric_set['tpr'])
+
+        for _, metric_set in enumerate(dummy_metrics):
+            plt.plot(
+                metric_set['fpr'],
+                metric_set['tpr'],
+                color='lightsteelblue',
+                linewidth=0.5,
+                alpha=0.5
+            )
+
+    for _, metric_set in enumerate(metrics):
+        plt.plot(metric_set['fpr'], metric_set['tpr'], color='lightsalmon', linewidth=0.5)
+
+
+    tpr_mu = all_tpr.mean(axis=0)
+    plt.plot(base_fpr, tpr_mu, 'r')
+
+    if do_dummy:
+        dummy_mu = all_tpr_dummy.mean(axis=0)
+        plt.plot(base_fpr, dummy_mu, 'b')
+
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.xlabel('False Positive Rate')
+    plt.ylabel('True Positive Rate')
+    plt.tight_layout()
+    plt.grid()
+    plt.savefig(save_name)
+
+def create_mia_report(
+    metadata: Dict,
+    mia_metrics: List,
+    dummy_metrics: List,
+    dummy_metadata: Dict) -> FPDF:
+    '''make a worst case membership inference report
+
+    Parameters
+    ----------
+
+    metadata: dict
+        dictionary of metadata
+
+    mia_metrics: list
+        list of metrics dictionaries
+
+    dummy_metrics: list
+        list of metrics dictionaries
+
+    dummy_metadata: dict
+        metadata for dummy experiments
+
+    Returns
+    -------
+
+    pdf: fpdf.FPDF
+        fpdf document object
+
+    '''
+
+    if dummy_metrics is None or len(dummy_metrics) == 0:
+        do_dummy = False
+    else:
+        do_dummy = True
+
+    _roc_plot(mia_metrics, dummy_metrics, "log_roc.png")
+
+
+    pdf = FPDF()
+    pdf.add_page()
+    pdf.set_xy(0, 0)
+    _title(pdf, "WorstCase MIA attack result report")
+    _subtitle(pdf, "Introduction")
+    _line(pdf, INTRODUCTION)
+    _subtitle(pdf, "Experiment summary")
+    for key, value in metadata['experiment_details'].items():
+        _line(pdf, f"{key:>30s}: {str(value):30s}", font="courier")
+    _subtitle(pdf, "Global metrics")
+    for key, value in metadata['global_metrics'].items():
+        _line(pdf, f"{key:>30s}: {str(value):30s}", font="courier")
+    if do_dummy:
+        _subtitle(pdf, "Baseline global metrics")
+        for key, value in dummy_metadata['global_metrics'].items():
+            _line(pdf, f"{key:>30s}: {str(value):30s}", font="courier")
+
+    _subtitle(pdf, "Metrics")
+    _line(
+        pdf,
+        "The following show summaries of the attack metrics over the repetitions",
+        font="arial"
+    )
+    for metric in DISPLAY_METRICS:
+        vals = np.array([m[metric] for m in mia_metrics])
+        if metric in MAPPINGS:
+            vals = np.array([MAPPINGS[metric](v) for v in vals])
+        text = (
+            f"{metric:>12} mean = {vals.mean():.2f}, var = {vals.var():.4f}, "
+            f"min = {vals.min():.2f}, max = {vals.max():.2f}"
+        )
+        _line(pdf, text, font='courier')
+
+    if do_dummy:
+        _subtitle(pdf, "Baseline metrics")
+        _line(pdf, (
+            'The following show summaries of the attack metrics over the repetitions where there '
+            'is no statistical difference between predictions in the training and test sets. '
+            'Simulation was done with training and test set sizes equal to the real ones'
+        ), font="arial")
+        for metric in DISPLAY_METRICS:
+            vals = np.array([m[metric] for m in dummy_metrics])
+            if metric in MAPPINGS:
+                vals = np.array([MAPPINGS[metric](v) for v in vals])
+            text = (
+                f"{metric:>12} mean = {vals.mean():.2f}, var = {vals.var():.4f}, "
+                f"min = {vals.min():.2f}, max = {vals.max():.2f}"
+            )
+            _line(pdf, text, font='courier')
+
+    pdf.add_page()
+    _subtitle(pdf, "Log ROC")
+    pdf.image('log_roc.png', x = None, y = None, w = 0, h = 140, type = '', link = '')
+    pdf.set_font('arial', '', 12)
+    _line(pdf, LOGROC_CAPTION)
+
+    pdf.add_page()
+    _title(pdf, "Glossary")
+    _write_dict(pdf, GLOSSARY)
+
+
+    return pdf
+
+def create_json_report(mia_metrics, dummy_metrics):
+    '''Create a report in json format for injestion by other tools
+    '''
+    # Initial work, just dump mia_metrics and dummy_metrics into a json structure
+    json_payload = {'metrics': mia_metrics, 'baseline_metrics': dummy_metrics}
+    return json.dumps(json_payload, cls=NumpyArrayEncoder)
+
+def create_lr_report(metadata, mia_metrics):
+    """TODO"""
+    _roc_plot_single(mia_metrics, "log_roc.png")
+    pdf = FPDF()
+    pdf.add_page()
+    pdf.set_xy(0, 0)
+    _title(pdf, "Likelihood Ratio Attack Report")
+    _subtitle(pdf, "Introduction")
+    _subtitle(pdf, "Metadata")
+    for key, value in metadata.items():
+        _line(pdf, f"{key:>30s}: {str(value):30s}", font="courier")
+    _subtitle(pdf, "Metrics")
+    sub_metrics_dict = {key: val for key, val in mia_metrics.items() if isinstance(val, float)}
+    for key, value in sub_metrics_dict.items():
+        if key in MAPPINGS:
+            value = MAPPINGS[key](value)
+        _line(pdf, f"{key:>30s}: {value:.4f}", font="courier")
+    pdf.add_page()
+    _subtitle(pdf, "ROC Curve")
+    pdf.image('log_roc.png', x = None, y = None, w = 0, h = 140, type = '', link = '')
+    return pdf
diff --git a/attacks/worst_case_attack.py b/attacks/worst_case_attack.py
new file mode 100644
index 00000000..8558e6db
--- /dev/null
+++ b/attacks/worst_case_attack.py
@@ -0,0 +1,423 @@
+'''
+worst_case_attack.py
+
+Runs a worst case attack based upon predictive probabilities stored in two .csv files
+'''
+
+import argparse
+import logging
+
+from typing import Tuple, Dict
+
+import numpy as np
+from sklearn.ensemble import RandomForestClassifier
+from sklearn.model_selection import train_test_split
+from scipy.stats import norm
+
+from attacks import metrics
+from attacks import report
+
+logging.basicConfig(level=logging.INFO)
+
+P_THRESH = 0.05
+
+def _generate_array(n_rows: int, beta: float) -> np.ndarray:
+    """Generate a single array of predictions, used when doing baseline experiments
+
+    Parameters
+    ----------
+    n_rows: int
+        the number of rows worth of data to generate
+    beta: float
+        the beta parameter for sampling probabilities
+
+    Returns
+    -------
+    preds: np.ndarray
+        Array of predictions. Two columns, n_rows rows
+
+    Notes
+    -----
+
+    Examples
+    --------
+
+    """
+
+    preds = np.zeros((n_rows, 2), float)
+    for row_idx in range(n_rows):
+        train_class = np.random.choice(2)
+        train_prob = np.random.beta(1, beta)
+        preds[row_idx, train_class] = train_prob
+        preds[row_idx, 1 - train_class] = 1 - train_prob
+    return preds
+
+def generate_arrays(
+    n_rows_in: int,
+    n_rows_out: int,
+    train_beta: float=2,
+    test_beta: float=2) -> Tuple[np.ndarray, np.ndarray]:
+    """Generate train and test prediction arrays, used when computing baseline
+
+    Parameters
+    ----------
+    n_rows_in: int
+        number of rows of in-sample (training) probabilities
+    n_rows_out: int
+        number of rows of out-of-sample (testing) probabilities
+    train_beta: float
+        beta value for generating train probabilities
+    test_beta: float:
+        beta_value for generating test probabilities
+
+    Returns
+    -------
+    train_preds: np.ndarray
+        Array of train predictions (n_rows x 2 columns)
+    test_preds: np.ndarray
+        Array of test predictions (n_rows x 2 columns)
+    """
+    train_preds = _generate_array(n_rows_in, train_beta)
+    test_preds = _generate_array(n_rows_out, test_beta)
+    return train_preds, test_preds
+
+def _make_dummy_data(args: Dict) -> None:
+    """Makes dummy data for testing functionality
+
+    Parameters
+    ----------
+    args: dict
+        Command line arguments
+
+    Returns
+    -------
+
+    Notes
+    -----
+
+    Returns nothing but saves two .csv files
+    """
+    logger = logging.getLogger("dummy-data")
+    logger.info("Making dummy data with %d rows in and %d out", args.n_rows_in, args.n_rows_out)
+    logger.info("Generating rows")
+    train_preds, test_preds = generate_arrays(
+        args.n_rows_in,
+        args.n_rows_out,
+        train_beta=args.train_beta,
+        test_beta=args.test_beta
+    )
+    logger.info("Saving files")
+    np.savetxt("train_preds.csv", train_preds, delimiter=",")
+    np.savetxt("test_preds.csv", test_preds, delimiter=",")
+
+def _run_attack(args, make_report=True):
+    """Method that runs a worst case attack
+    """
+    logger = logging.getLogger("run-attack")
+    logger.info("Loading train predictions form %s", args.in_sample_filename)
+    train_preds = np.loadtxt(args.in_sample_filename, delimiter=",")
+    logger.info("Loaded %d rows", len(train_preds))
+    logger.info("Loading test predictions form %s", args.in_sample_filename)
+    test_preds = np.loadtxt(args.out_sample_filename, delimiter=",")
+    logger.info("Loaded %d rows", len(test_preds))
+    mia_metrics, metadata = attack(args, train_preds, test_preds)
+
+    # do some baseline attacks
+    logger.info("Performing baseline attacks")
+
+    dummy_metrics = []
+    dummy_metadata = None
+    for _ in range(args.dummy_reps):
+        dummy_train, dummy_test = generate_arrays(len(train_preds), len(test_preds))
+        temp_dummy_metrics, dummy_metadata = attack(args, dummy_train, dummy_test)
+        dummy_metrics += temp_dummy_metrics
+
+    if args.dummy_reps > 0:
+        metadata['experiment_details']['n_baseline_experiments_done'] = (
+            f'Done {args.n_reps} reps for {args.dummy_reps} splits of synthetic '
+            f'data (total = {len(dummy_metrics)})'
+        )
+    else:
+        metadata['experiment_details']['n_baseline_experiments_done'] = "None"
+
+    if make_report:
+        pdf = report.create_mia_report(metadata, mia_metrics, dummy_metrics, dummy_metadata)
+        pdf.output(f"{args.report_name}.pdf", 'F')
+        json_details = report.create_json_report(mia_metrics, dummy_metrics)
+        with open(f"{args.report_name}.json", 'w', encoding="utf-8") as json_file:
+            json_file.write(json_details)
+
+def _get_n_significant(p_val_list, p_thresh, bh_fdr_correction=False):
+    '''
+    Helper method to determine if values within a list of p-values are significant at
+    p_thresh. Can peform multiple testing correction.
+    '''
+    if not bh_fdr_correction:
+        return sum([1 for p in p_val_list if p <= p_thresh]) # pylint: disable = consider-using-generator
+    p_val_list = np.asarray(sorted(p_val_list))
+    n_vals = len(p_val_list)
+    hoch_vals = np.array(
+        [(k / n_vals) * P_THRESH for k in range(1, n_vals + 1)]
+    )
+    bh_sig_list = p_val_list <= hoch_vals
+    if any(bh_sig_list):
+        n_sig_bh = (np.where(bh_sig_list)[0]).max() + 1
+    else:
+        n_sig_bh = 0
+    return n_sig_bh
+
+
+def attack( # pylint: disable=too-many-locals
+    args: Dict,
+    train_preds: np.ndarray,
+    test_preds: np.ndarray) -> Tuple[dict, dict]:
+    """Runs the attack based upon the predictions in train_preds and test_preds, and params
+    in args
+
+    Parameters
+    ----------
+
+    args: dict
+        Command line arguments
+    train_preds: np.ndarray
+        Array of train predictions. One row per example, one column per class (i.e. 2)
+    test_preds:  np.ndarray
+        Array of test predictions. One row per example, one column per class (i.e. 2)
+
+    Returns
+    -------
+    mia_metrics: dict
+        dictionary of metrics
+    metadata: dict
+        dictionary of metadata
+
+    """
+    logger = logging.getLogger("attack")
+    logger.info("Sorting probabilities to leave highest value in first column")
+    train_preds = -np.sort(-train_preds, axis=1)
+    test_preds = -np.sort(-test_preds, axis=1)
+
+    logger.info("Creating MIA data")
+    mi_x = np.vstack((train_preds, test_preds))
+    mi_y = np.hstack(
+        (
+            np.ones(len(train_preds)),
+            np.zeros(len(test_preds))
+        )
+    )
+
+    mia_metrics = []
+    for rep in range(args.n_reps):
+        logger.info("Rep %d of %d", rep, args.n_reps)
+        mi_train_x, mi_test_x, mi_train_y, mi_test_y = train_test_split(
+            mi_x, mi_y, test_size=args.test_prop, stratify=mi_y
+        )
+        attack_classifier = RandomForestClassifier()
+        attack_classifier.fit(mi_train_x, mi_train_y)
+
+        mia_metrics.append(
+            metrics.get_metrics(
+                attack_classifier,
+                mi_test_x,
+                mi_test_y
+            )
+        )
+
+    logger.info("Finished simulating attacks")
+    metadata = {
+        'global_metrics': {},
+        'experiment_details': {
+            'n_reps': args.n_reps,
+            'n_in_sample': len(train_preds),
+            'n_out_sample': len(test_preds),
+            'attack model': "Random Forest",
+            'in_sample_filename': args.in_sample_filename,
+            'out_sample_filename': args.out_sample_filename,
+            'p_val_thresh': args.p_thresh
+        }
+    }
+
+    # Compute NULL AUC
+    auc_std = np.sqrt(
+        metrics.expected_auc_var(0.5, mi_test_y.sum(), len(mi_test_y) - mi_test_y.sum())
+    )
+    # Assuming auc is normal, compute the probability of the NULL generating an AUC higher
+    # than the NULL
+    auc_p_vals = [1 - norm.cdf(m['AUC'], loc=0.5, scale=auc_std) for m in mia_metrics]
+
+    metadata['global_metrics']['null_auc_3sd_range'] = \
+        f"{0.5 - 3*auc_std:.4f} -> {0.5 + 3*auc_std:.4f}"
+    metadata['global_metrics']['n_sig_auc_p_vals'] = _get_n_significant(auc_p_vals, args.p_thresh)
+    metadata['global_metrics']['n_sig_auc_p_vals_corrected'] = _get_n_significant(
+        auc_p_vals,
+        args.p_thresh,
+        bh_fdr_correction=True
+    )
+
+    pdif_vals = [np.exp(-m['PDIF01']) for m in mia_metrics]
+    metadata['global_metrics']['n_sig_pdif_vals'] = _get_n_significant(pdif_vals, args.p_thresh)
+    metadata['global_metrics']['n_sig_pdif_vals_corrected'] = _get_n_significant(
+        pdif_vals,
+        args.p_thresh,
+        bh_fdr_correction=True
+    )
+
+    return mia_metrics, metadata
+
+
+def main():
+    '''main method to parse arguments and invoke relevant method'''
+
+    parser = argparse.ArgumentParser(description=(
+        'Perform a worst case attack from saved model predictions'
+        )
+    )
+
+    subparsers = parser.add_subparsers()
+    dummy_parser = subparsers.add_parser('make-dummy-data')
+    dummy_parser.add_argument('--num-rows-in',
+        action='store',
+        dest='n_rows_in',
+        type=int,
+        required=False,
+        default=100,
+        help=(
+            'How many rows to generate in the in-sample file. Default = %(default)d'
+        )
+    )
+
+    dummy_parser.add_argument('--num-rows-out',
+        action='store',
+        dest='n_rows_out',
+        type=int,
+        required=False,
+        default=100,
+        help=(
+            'How many rows to generate in the out-of-sample file. Default = %(default)d'
+        )
+    )
+
+    dummy_parser.add_argument('--train-beta',
+        action='store',
+        type=float,
+        required=False,
+        default=5,
+        dest="train_beta",
+        help=(
+            'Value of b parameter for beta distribution used to sample the in-sample '
+            'probabilities. '
+            'High values will give more extreme probabilities. Set this value higher than '
+            '--test-beta to see successful attacks. Default = %(default)f'
+        )
+    )
+
+    dummy_parser.add_argument('--test-beta',
+        action='store',
+        type=float,
+        required=False,
+        default=2,
+        dest="test_beta",
+        help=(
+            'Value of b parameter for beta distribution used to sample the out-of-sample '
+            'probabilities. High values will give more extreme probabilities. Set this value '
+            'lower than --train-beta to see successful attacks. Default = %(default)f'
+        )
+    )
+
+    dummy_parser.set_defaults(func=_make_dummy_data)
+
+    attack_parser = subparsers.add_parser('run-attack')
+    attack_parser.add_argument('-i', '--in-sample-preds',
+        action='store',
+        dest='in_sample_filename',
+        required=False,
+        type=str,
+        default="train_preds.csv",
+        help=(
+            'csv file containing the predictive probabilities (one column per class) for the '
+            'training data (one row per training example). Default = %(default)s'
+        )
+    )
+
+    attack_parser.add_argument('-o', '--out-of-sample-preds',
+        action='store',
+        dest='out_sample_filename',
+        required=False,
+        type=str,
+        default="test_preds.csv",
+        help=(
+            'csv file containing the predictive probabilities (one column per class) for the '
+            'non-training data (one row per training example). Default = %(default)s'
+        )
+    )
+
+    attack_parser.add_argument('-r', '--n-reps',
+        type=int,
+        required=False,
+        default=5,
+        action="store",
+        dest="n_reps",
+        help=(
+            'Number of repetitions (splitting data into attack model training and testing '
+            'paritions to perform. Default = %(default)d'
+        )
+    )
+
+    attack_parser.add_argument('-t', '--test-prop',
+        type=float,
+        required=False,
+        default=0.3,
+        action="store",
+        dest="test_prop",
+        help=(
+            'Proportion of examples to be used for testing when fiting the attack model. '
+            'Default = %(default)f'
+        )
+    )
+
+    attack_parser.add_argument('--report-name',
+        type=str,
+        action='store',
+        dest='report_name',
+        default='mia_report',
+        required=False,
+        help=(
+            'Filename for the report output. Default = %(default)s. Code will append .pdf and .json'
+        )
+    )
+
+    attack_parser.add_argument('--dummy-reps',
+        type=int,
+        action='store',
+        dest='dummy_reps',
+        default=1,
+        required=False,
+        help=(
+            'Number of dummy datasets to sample. Each will be assessed with --n-reps train and '
+            'test splits. Set to 0 to do no baseline calculations. Default = %(default)d'
+        )
+    )
+
+    attack_parser.add_argument('--p-thresh',
+        action='store',
+        type=float,
+        default=P_THRESH,
+        required=False,
+        dest='p_thresh',
+        help=(
+            'P-value threshold for significance testing. Default = %(default)f'
+        )
+    )
+
+    attack_parser.set_defaults(func=_run_attack)
+
+
+    args = parser.parse_args()
+    print(args)
+    try:
+        args.func(args)
+    except AttributeError:
+        print("Invalid command. Try --help to get more details")
+
+if __name__ == '__main__':
+    main()
diff --git a/data_preprocessing/data_interface.py b/data_preprocessing/data_interface.py
new file mode 100644
index 00000000..36976c24
--- /dev/null
+++ b/data_preprocessing/data_interface.py
@@ -0,0 +1,616 @@
+'''
+data_interface.py
+A set of useful handlers to pull in datasets common to the project and perform the appropriate
+pre-processing
+'''
+
+# pylint: disable=import-error, invalid-name, consider-using-with, too-many-return-statements
+
+import os
+import logging
+from typing import Tuple, List
+from zipfile import ZipFile, BadZipFile
+from collections import Counter
+import pandas as pd
+import numpy as np
+import pylab as plt
+from sklearn.preprocessing import OneHotEncoder, LabelEncoder
+from sklearn.datasets import fetch_openml
+from sklearn.datasets import load_iris
+
+
+# Following is to stop pylint always sqwarking at pandas things
+# pylint: disable = no-member, unsubscriptable-object
+
+
+logging.basicConfig(
+    level = "DEBUG"
+)
+
+logger = logging.getLogger(__file__)
+
+
+PROJECT_ROOT_FOLDER = os.path.dirname(os.path.dirname(__file__))
+logger.info("ROOT PROJECT FOLDER = %s", PROJECT_ROOT_FOLDER)
+
+
+class UnknownDataset(Exception):
+    '''Exception raised if the user passes a name that we don't recognise
+    '''
+
+class DataNotAvailable(Exception):
+    '''Exception raised if the user asks for a dataset that they do not have the data for. I.e.
+    some datasets require a .csv file to have been downloaded.
+    '''
+
+def get_data_sklearn( # pylint: disable = too-many-branches
+    dataset_name: str,
+    data_folder: str = os.path.join(PROJECT_ROOT_FOLDER, "data")
+    ) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    """Main entry method to return data in format sensible for sklearn. User passes a name and that
+    dataset is returned as a tuple of pandas DataFrames (data, labels).
+
+    Parameters
+    ----------
+    dataset_name: str
+        The name of the dataset to load
+    data_folder: str
+        The name of the local folder in which data is stored.
+
+    Returns
+    -------
+
+    X: pd.DataFrame
+        The input dataframe -- rows are examples, columns variables
+    y: pd.DataFrame
+        the target dataframe -- has a single column containing the target values
+
+    Notes
+    -----
+
+    The following datasets are available:
+    mimic2-iaccd (requires data download)
+    in-hospital-mortality (requires data download)
+    medical-mnist-ab-v-br-100 (requires data download)
+    medical-mnist-ab-v-br-500 (requires data download)
+    medical-mnist-all-100 (requires data download)
+    indian liver (requires data download)
+    texas hospitals 10 (requires data download)
+    synth-ae (requires data download)
+    synth-ae-small (requires data download)
+    nursery (downloads automatically)
+    iris (available out of the box via sklearn)
+
+    Datasets can be normalised by adding the following prefixes:
+    standard: standardises all columns to have zero mean and unit variance.
+    minmax: standardises all columns to have values between 0 and 1.
+    round: rounds continues features to have 3dp
+
+    These can be nested.
+
+    Examples
+    --------
+    >>> X, y = get_data_sklearn("mimic2-iaccd") # pull the mimic2-iaccd data
+    >>> X, y = get_data_sklearn("minmax iris") # pull the iris data and round continuous features
+
+
+    """
+    logger.info("DATASET FOLDER = %s", data_folder)
+
+    if dataset_name.startswith("standard"):
+        sub_name = dataset_name.split("standard")[1].strip()
+        feature_df, target_df = get_data_sklearn(sub_name)
+        for column in feature_df.columns:
+            col_mean = feature_df[column].mean()
+            col_std = np.sqrt(feature_df[column].var())
+            feature_df[column] = feature_df[column] - col_mean
+            feature_df[column] = feature_df[column] / col_std
+        return feature_df, target_df
+
+    if dataset_name.startswith("minmax"):
+        sub_name = dataset_name.split("minmax")[1].strip()
+        feature_df, target_df = get_data_sklearn(sub_name)
+        for column in feature_df.columns:
+            col_min = feature_df[column].min()
+            col_range = feature_df[column].max() - col_min
+            feature_df[column] = feature_df[column] - col_min
+            feature_df[column] = feature_df[column] / col_range
+        return feature_df, target_df
+
+    if dataset_name.startswith("round"):
+        sub_name = dataset_name.split("round")[1].strip()
+        logger.debug(sub_name)
+        feature_df, target_df = get_data_sklearn(sub_name)
+        column_dtype = feature_df.dtypes # pylint: disable = no-member
+
+        for i, column in enumerate(feature_df.columns):
+            if column_dtype[i] == "float64":
+                feature_df[column] = feature_df[column].round(decimals=3)
+        return feature_df, target_df
+
+
+
+    if dataset_name == 'mimic2-iaccd':
+        return _mimic_iaccd(data_folder)
+    if dataset_name == 'in-hospital-mortality':
+        return _in_hospital_mortality(data_folder)
+    if dataset_name == 'medical-mnist-ab-v-br-100':
+        return _medical_mnist_loader(data_folder, 100, ['AbdomenCT', 'BreastMRI'])
+    if dataset_name == 'medical-mnist-ab-v-br-500':
+        return _medical_mnist_loader(data_folder, 500, ['AbdomenCT', 'BreastMRI'])
+    if dataset_name == 'medical-mnist-all-100':
+        return _medical_mnist_loader(
+            data_folder,
+            100,
+            ['AbdomenCT', 'BreastMRI', 'CXR', 'ChestCT', 'Hand', 'HeadCT']
+        )
+    if dataset_name == 'indian liver':
+        return _indian_liver(data_folder)
+    if dataset_name == 'texas hospitals 10':
+        return _texas_hospitals(data_folder)
+    if dataset_name == 'synth-ae':
+        return _synth_ae(data_folder)
+    if dataset_name == 'synth-ae-small':
+        x, y = _synth_ae(data_folder, 200)
+        return x, y
+    if dataset_name == 'nursery':
+        return _nursery()
+    if dataset_name == 'iris':
+        return _iris()
+    raise UnknownDataset(dataset_name)
+
+
+def _iris() -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    sklearn iris data - just first two classes
+    '''
+    X, y = load_iris(return_X_y=True, as_frame=True)
+    X = X[y < 2]
+    y = y[y < 2]
+    return X, pd.DataFrame(y)
+
+def _nursery() -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    The sklearn nursery dataset
+    '''
+
+    data = fetch_openml(data_id=26, as_frame=True)
+
+    target_encoder = LabelEncoder()
+    target_vals = target_encoder.fit_transform(data['target'].values)
+    target_dataframe = pd.DataFrame({'target': target_vals})
+
+    feature_encoder = OneHotEncoder()
+    x_encoded = feature_encoder.fit_transform(data['data']).toarray()
+    feature_dataframe = pd.DataFrame(x_encoded, columns=feature_encoder.get_feature_names_out())
+
+
+    return feature_dataframe, target_dataframe
+
+
+# Patched to support non-flattened images. Same behaviour as before except if called with
+# flatten=False explicitly.
+def _images_to_ndarray(
+    images_dir: str,
+    number_to_load: int,
+    label: int,
+    flatten: bool = True) -> Tuple[np.array, np.array]:
+    '''
+    Grab number_to_load images from the images_dir and create a np array and label array
+    '''
+    folder_path = images_dir + os.sep
+    images_names = sorted(os.listdir(folder_path))
+    images_names = images_names[:number_to_load]
+    if flatten:
+        np_images = np.array([plt.imread(folder_path + img).flatten() for img in images_names])
+    else:
+        np_images = np.array([plt.imread(folder_path + img) for img in images_names])
+    labels = np.ones((len(np_images), 1), int) * label
+    return(np_images, labels)
+
+def _medical_mnist_loader( #pylint: disable = too-many-locals
+    data_folder: str,
+    n_per_class: int,
+    classes: List[str]) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    Load Medical MNIST into pandas format
+    borrows heavily from: https://www.kaggle.com/harelshattenstein/medical-mnist-knn
+    Creates a binary classification
+    '''
+
+    base_folder = os.path.join(
+        data_folder,
+        'kaggle-medical-mnist',
+        'archive',
+    )
+
+    zip_file = os.path.join(
+        data_folder,
+        'kaggle-medical-mnist',
+        "archive.zip"
+    )
+
+    logger.debug(base_folder, data_folder)
+    if not any([os.path.exists(base_folder), os.path.exists(zip_file)]):
+        help_message = f"""
+Data file {base_folder} does not exist. Please download fhe file from:
+https://www.kaggle.com/andrewmvd/medical-mnist 
+and place it in the correct folder. It unzips the file first.
+        """
+        raise DataNotAvailable(help_message)
+
+    if os.path.exists(base_folder):
+        pass
+    elif os.path.exists(zip_file):
+        try:
+            with ZipFile(zip_file) as zip_handle:
+                zip_handle.extractall()
+                logger.debug("Extracted all")
+        except BadZipFile:
+            logger.error("Encountered bad zip file")
+            raise
+
+    labels_dict = {
+        0 : 'AbdomenCT',
+        1 : 'BreastMRI',
+        2 : 'CXR',
+        3 : 'ChestCT',
+        4 : 'Hand',
+        5 : 'HeadCT'
+    }
+
+    reverse_labels_dict = {v: k for k, v in labels_dict.items()}
+
+    for i, class_name in enumerate(classes):
+        label = reverse_labels_dict[class_name]
+        x_images, y_images = _images_to_ndarray(
+            os.path.join(base_folder, class_name),
+            n_per_class,
+            label
+        )
+
+        if i == 0:
+            all_x = x_images
+            all_y = y_images
+        else:
+            all_x = np.vstack((all_x, x_images))
+            all_y = np.vstack((all_y, y_images))
+
+    return (pd.DataFrame(all_x), pd.DataFrame(all_y))
+
+
+def _synth_ae(data_folder: str, n_rows: int=5000) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    First norws (default 5000) rows from the Synthetic A&E data from NHS England
+    https://data.england.nhs.uk/dataset/a-e-synthetic-data/resource/81b068e5-6501-4840-a880-a8e7aa56890e # pylint: disable=line-too-long
+    '''
+
+    file_path = os.path.join(
+        data_folder,
+        'AE_England_synthetic.csv' #'A&E Synthetic Data.csv'
+    )
+
+    if not os.path.exists(file_path):
+        help_message = f"""
+Data file {file_path} does not exist. Please download the file from:
+https://data.england.nhs.uk/dataset/a-e-synthetic-data/resource/81b068e5-6501-4840-a880-a8e7aa56890e
+unzip it (7z) and then copy the .csv file into your data folder.
+    """
+        raise DataNotAvailable(help_message)
+
+    input_data = pd.read_csv(file_path, nrows = n_rows)
+    columns_to_drop = [
+        'AE_Arrive_Date', 'AE_Arrive_HourOfDay', 'Admission_Method',
+        'ICD10_Chapter_Code', 'Treatment_Function_Code', 'Length_Of_Stay_Days',
+        'ProvID'
+    ]
+    input_data.drop(columns_to_drop, axis=1, inplace=True)
+
+    # Remove any rows with NAs in the remaining columns
+    input_data.dropna(axis=0, inplace=True)
+
+    # One-hot encode some columns
+    encode_columns = ['Age_Band', 'AE_HRG']
+    encode_data = input_data[encode_columns]
+    input_data.drop(encode_columns, axis=1, inplace=True)
+
+    oh = OneHotEncoder()
+    oh.fit(encode_data)
+    onehot_df = pd.DataFrame(
+        oh.transform(encode_data).toarray(),
+        columns=oh.get_feature_names_out(),
+        index=input_data.index
+    )
+
+    input_data = pd.concat([input_data, onehot_df], axis=1)
+
+    X = input_data.drop(['Admitted_Flag'], axis=1)
+    y = input_data[['Admitted_Flag']]
+
+    return (X, y)
+
+
+def _indian_liver(data_folder: str) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    Indian Liver Patient Dataset
+    https://archive.ics.uci.edu/ml/machine-learning-databases/00225/Indian%20Liver%20Patient%20Dataset%20(ILPD).csv # pylint: disable=line-too-long
+    '''
+    #(https://archive.ics.uci.edu/ml/datasets/ILPD+(Indian+Liver+Patient+Dataset)
+    file_path = os.path.join(
+        data_folder,
+        "Indian Liver Patient Dataset (ILPD).csv"
+    )
+    if not os.path.exists(file_path):
+        help_message = f"""
+Data file {file_path} does not exist. Please download fhe file from:
+https://archive.ics.uci.edu/ml/datasets/ILPD+(Indian+Liver+Patient+Dataset
+and place it in the correct folder.
+        """
+        raise DataNotAvailable(help_message)
+
+    column_names= [
+        'age', 'gender', 'total Bilirubin', 'direct Bilirubin', 'Alkphos',
+        'SGPT', 'SGOT', 'total proteins', 'albumin', 'A/G ratio',  'class'
+    ]
+
+    liver_data = pd.read_csv(file_path, names=column_names, index_col=False)
+
+    liver_data.gender.replace('Male', 0, inplace=True)
+    liver_data.gender.replace('Female', 1, inplace=True)
+
+    liver_data.dropna(axis=0, inplace=True)
+
+    liver_labels = liver_data['class']
+    liver_data.drop(['class'], axis=1, inplace=True)
+
+    label_encoder = LabelEncoder()
+    encoded_labels = label_encoder.fit_transform(liver_labels.values)
+    liver_labels = pd.DataFrame({'class': encoded_labels})
+    return (liver_data, liver_labels)
+
+
+def _in_hospital_mortality(data_folder: str) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    In-hospital mortality data from this study:
+        https://datadryad.org/stash/dataset/doi:10.5061/dryad.0p2ngf1zd
+    '''
+    # Check the data has been downloaded. If not throw an exception with instructions on how to
+    # download, and where to store
+    files = ["data01.csv", "doi_10.5061_dryad.0p2ngf1zd__v5.zip"]
+    file_path = [os.path.join(data_folder, f) for f in files]
+
+    if not any([os.path.exists(fp) for fp in file_path]): # pylint: disable=use-a-generator
+        help_message = f"""
+Data file {file_path[0]} or {file_path[1]} does not exist. Please download the file from:
+https://datadryad.org/stash/dataset/doi:10.5061/dryad.0p2ngf1zd
+and place it in the correct folder. It works with either the zip file or uncompressed.
+        """
+        raise DataNotAvailable(help_message)
+
+    if os.path.exists(file_path[1]):
+        input_data = pd.read_csv(ZipFile(file_path[1]).open("data01.csv"))
+    else:
+        input_data = pd.read_csv(file_path[0])
+    clean_data = input_data.dropna(axis=0, how='any').drop(columns=["group", "ID"])
+    target = 'outcome'
+    labels = clean_data[target]
+    features = clean_data.drop([target], axis=1)
+
+    label_encoder = LabelEncoder()
+    encoded_labels = label_encoder.fit_transform(labels.values)
+    labels = pd.DataFrame({'outcome': encoded_labels})
+
+    return (features, labels)
+
+
+
+
+def _mimic_iaccd(data_folder: str) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    '''
+    Loads the mimic_iaccd data and performs Alba's pre-processing
+    '''
+
+    # Check the data has been downloaded. If not throw an exception with instructions on how to
+    # download, and where to store
+    file_path = os.path.join(
+        data_folder,
+        "mimic2-iaccd",
+        "1.0",
+        "full_cohort_data.csv"
+    )
+
+    if not os.path.exists(file_path):
+        help_message = f"""
+The MIMIC2-iaccd data is not available in {data_folder}.  
+The following file should exist: {file_path}.
+Please download from https://physionet.org/content/mimic2-iaccd/1.0/full_cohort_data.csv
+        """
+        raise DataNotAvailable(help_message)
+
+    # File exists, load and preprocess#
+    logger.info("Loading mimic2-iaccd")
+    input_data = pd.read_csv(file_path)
+
+    logger.info("Preprocessing")
+    # remove columns non-numerical and repetitive or uninformative data for the analysis
+    col = ['service_unit', 'day_icu_intime', 'hosp_exp_flg','icu_exp_flg', 'day_28_flg']
+    # service_num is the numerical version of service_unit
+    # day_icu_intime_num is the numerical version of day_icu_intime
+    # the other columns are to do with death and are somewhat repetitive with censor_flg
+    input_data.drop(col, axis = 1, inplace=True)
+    # drop columns with only 1 value
+    input_data.drop('sepsis_flg', axis=1, inplace=True)
+    # drop NA by row
+    input_data.dropna(axis=0, inplace=True)
+
+    # extract target
+    target = 'censor_flg'
+    y = input_data[target]
+    X = input_data.drop([target, 'mort_day_censored'], axis=1)
+
+    label_encoder = LabelEncoder()
+    encoded_labels = label_encoder.fit_transform(y.values)
+    y = pd.DataFrame({'censor_flag': encoded_labels})
+
+
+    return (X, y)
+
+def _texas_hospitals(data_folder: str) -> Tuple[pd.DataFrame, pd.DataFrame]: # pylint: disable=too-many-statements, too-many-locals
+    '''
+    Texas Hospitals Dataset
+    (https://www.dshs.texas.gov/THCIC/Hospitals/Download.shtm)
+    '''
+    file_list = ["PUDF-1Q2006-tab-delimited.zip",
+                 "PUDF-1Q2007-tab-delimited.zip",
+                 "PUDF-1Q2008-tab-delimited.zip",
+                 "PUDF-1Q2009-tab-delimited.zip",
+                 "PUDF-2Q2006-tab-delimited.zip",
+                 "PUDF-2Q2007-tab-delimited.zip",
+                 "PUDF-2Q2008-tab-delimited.zip",
+                 "PUDF-2Q2009-tab-delimited.zip",
+                 "PUDF-3Q2006-tab-delimited.zip",
+                 "PUDF-3Q2007-tab-delimited.zip",
+                 "PUDF-3Q2008-tab-delimited.zip",
+                 "PUDF-3Q2009-tab-delimited.zip",
+                 "PUDF-4Q2006-tab-delimited.zip",
+                 "PUDF-4Q2007-tab-delimited.zip",
+                 "PUDF-4Q2008-tab-delimited.zip",
+                 "PUDF-4Q2009-tab-delimited.zip"]
+
+    files_path = [os.path.join(
+        data_folder,"TexasHospitals",f) for f in file_list]
+
+    found = [os.path.exists(file_path) for file_path in files_path]
+    not_found = [file_path for file_path in files_path if not os.path.exists(file_path)]
+
+    processed_data_file = "texas_data10_rm_binary.csv"
+    if not all(found):
+        help_message = f"""
+    Some or all data files do not exist. Please accept their terms & conditions,then download the
+    tab delimited files from each quarter during 2006-2009 from:
+    https://www.dshs.texas.gov/THCIC/Hospitals/Download.shtm
+and place it in the correct folder.
+
+    Missing files are:
+    {not_found}
+        """
+        raise DataNotAvailable(help_message)
+
+    if not os.path.exists(os.path.join(data_folder,"TexasHospitals", processed_data_file)):
+
+        logger.info("Processing Texas Hospitals data (2006-2009)")
+
+        #Load data
+        columns_names = [
+            'THCIC_ID',# Provider ID. Unique identifier assigned to the provider by DSHS.
+            # Hospitals with fewer than 50 discharges have been aggregated into the
+            # Provider ID '999999'
+            'DISCHARGE_QTR', #yyyyQm
+            'TYPE_OF_ADMISSION',
+            'SOURCE_OF_ADMISSION',
+            'PAT_ZIP',#Patient’s five-digit ZIP code
+            'PUBLIC_HEALTH_REGION', #Public Health Region of patient’s address
+            'PAT_STATUS', #Code indicating patient status as of the ending date of service for
+            # the period of care reported
+            'SEX_CODE',
+            'RACE',
+            'ETHNICITY',
+            'LENGTH_OF_STAY',
+            'PAT_AGE', #Code indicating age of patient in days or years on date of discharge.
+            'PRINC_DIAG_CODE', #diagnosis code for the principal diagnosis
+            'E_CODE_1', #external cause of injury
+            'PRINC_SURG_PROC_CODE', #Code for the principal surgical or other procedure performed
+            # during the period covered by the bill
+            'RISK_MORTALITY', #Assignment of a risk of mortality score from the All Patient
+            # Refined (APR) Diagnosis Related Group (DRG)
+            'ILLNESS_SEVERITY',#Assignment of a severity of illness score from the All Patient
+            # Refined (APR) Diagnosis RelatedGroup (DRG
+            'RECORD_ID'
+        ]
+        # obtain the 100 most frequent procedures
+        tmp = []
+        for f in files_path:
+            df = [
+                pd.read_csv(ZipFile(f).open(i), sep="\t", usecols=["PRINC_SURG_PROC_CODE"])\
+                     for i in ZipFile(f).namelist() if 'base' in i
+            ][0]
+            df.dropna(inplace=True)
+            tmp.extend(list(df.PRINC_SURG_PROC_CODE))
+        princ_surg_proc_keep = [k for k,v in Counter(tmp).most_common(10)]
+        #remove unecessary variables
+        del tmp
+
+        # Load the data
+        tx_data = pd.DataFrame()
+        for f in files_path:
+            df = [pd.read_csv(ZipFile(f).open(i), sep="\t", usecols=columns_names)
+                  for i in ZipFile(f).namelist() if 'base' in i][0]
+            #keep only those rows with one of the 10 most common principal surgical procedure
+            df = df[df["PRINC_SURG_PROC_CODE"].isin(princ_surg_proc_keep)]
+            #clean up data
+            df.dropna(inplace=True)
+            df.replace('`', pd.NA, inplace=True)
+            df.replace('*', pd.NA, inplace=True)
+            #replace sex to numeric
+            df.SEX_CODE.replace('M',0,inplace=True)
+            df.SEX_CODE.replace('F',1,inplace=True)
+            df.SEX_CODE.replace('U',2,inplace=True)
+            #set to numerical variable
+            for d_code in set(list(df.DISCHARGE_QTR)):
+                df.DISCHARGE_QTR.replace(d_code, ''.join(d_code.split('Q')), inplace=True)
+            df.dropna(inplace=True)
+            #merge data
+            tx_data = pd.concat([tx_data, df])
+        #remove uncessary variables
+        del df
+
+        #Risk moratality, make it binary
+        #1 Minor
+        #2 Moderate
+        #3 Major
+        #4 Extreme
+        tx_data.RISK_MORTALITY.astype(int)
+        tx_data.RISK_MORTALITY.replace(1,0,inplace=True)
+        tx_data.RISK_MORTALITY.replace(2,0,inplace=True)
+        tx_data.RISK_MORTALITY.replace(3,1,inplace=True)
+        tx_data.RISK_MORTALITY.replace(4,1,inplace=True)
+
+        #renumber non-numerical codes for cols
+        cols=['PRINC_DIAG_CODE', 'SOURCE_OF_ADMISSION', 'E_CODE_1']
+        for col in cols:
+            tmp = list(
+                set([x for x in tx_data[col] if not str(x).isdigit() and not isinstance(x, float)]) # pylint: disable=consider-using-set-comprehension
+            )
+            n = max(
+                list(
+                    set([int(x) for x in tx_data[col] if str(x).isdigit() or isinstance(x, float)]) # pylint: disable=consider-using-set-comprehension
+                )
+            )
+            for i, x in enumerate(tmp):
+                tx_data[col].replace(x, n + i, inplace=True)
+        del tmp, n
+        #set index
+        tx_data.set_index('RECORD_ID', inplace=True)
+        #final check and drop of NAs
+        tx_data.dropna(inplace=True)
+        #convert all data to numerical
+        tx_data = tx_data.astype(int)
+        #save csv file
+        tx_data.to_csv(os.path.join(data_folder, "TexasHospitals", processed_data_file))
+    else:
+        logger.info("Loading processed Texas Hospitals data (2006-2009) csv file.")
+        #load texas data processed csv file
+        tx_data = pd.read_csv(os.path.join(data_folder, "TexasHospitals", processed_data_file))
+
+    # extract target
+    var = 'RISK_MORTALITY'
+    labels = tx_data[var]
+    # Drop the column that contains the labels
+    tx_data.drop([var], axis=1, inplace=True)
+
+    label_encoder = LabelEncoder()
+    encoded_labels = label_encoder.fit_transform(labels.values)
+    labels = pd.DataFrame({var: encoded_labels})
+
+
+    return(tx_data, labels)
diff --git a/docs/source/attacks.rst b/docs/source/attacks.rst
index 39d1007a..65ade9ef 100644
--- a/docs/source/attacks.rst
+++ b/docs/source/attacks.rst
@@ -1,9 +1,7 @@
 Attacks
 =======
 
-.. automodule:: experiments
 .. automodule:: metrics.py
 .. automodule:: mia_extremecase.py
-.. automodule:: scenarios.py
-.. automodule:: tfwrapper.py
-
+.. automodule:: worst_case_attack.py
+.. automodule:: report.py
diff --git a/docs/source/datainterface.rst b/docs/source/datainterface.rst
new file mode 100644
index 00000000..f0781fbb
--- /dev/null
+++ b/docs/source/datainterface.rst
@@ -0,0 +1,8 @@
+
+Data Interface
+================
+
+      
+.. automodule:: data_preprocessing.data_interface
+   :members:
+
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 9771b5e5..e2d5ed78 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -16,13 +16,13 @@ This is the class and function reference for safemodel.
    :caption: Contents:
    
    introduction
+   datainterface
    attacks
    safemodel
    safedecisiontree
    saferandomforest
    safekeras
-   
-
+=======
 
 Indices and tables
 ==================
diff --git a/examples/worst_case_attack_example.py b/examples/worst_case_attack_example.py
new file mode 100644
index 00000000..6dfc43c7
--- /dev/null
+++ b/examples/worst_case_attack_example.py
@@ -0,0 +1,164 @@
+"""Examples for using the 'worst case' attack code.
+
+This code simulates a MIA attack providing the attacker with as much information as possible.
+i.e. they have a subset of rows that they _know_ were used for training. And a subset that they
+know were not. They also have query access to the target model.
+
+They pass the training and non-training rows through the target model to get the predictive
+probabilities. These are then used to train an _attack model_. And the attack model is evaluated
+to see how well it can predict whether or not other examples were in the training set or not.
+
+The code can be called from the command line, or accessed programmatically. Examples of both
+are shown below.
+
+In the code, [Researcher] and [TRE] are used in comments to denote which bit is done by whom
+
+Running
+-------
+
+Invoke this code from the root AI-SDC folder with
+python -m examples.worst_case_attack_example
+
+"""
+import os
+from types import SimpleNamespace
+
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.svm import SVC
+from sklearn.datasets import load_breast_cancer
+
+from attacks import worst_case_attack # pylint: disable = import-error
+
+# [Researcher] Access a dataset
+X, y = load_breast_cancer(return_X_y=True, as_frame=False)
+
+# [Researcher] Split into training and test sets
+train_X, test_X, train_y, test_y = train_test_split(X, y, test_size=0.3)
+
+# [Researcher] Define the classifier
+target_model = SVC(gamma=0.1, probability=True)
+
+# [Researcher] Train the classifier
+target_model.fit(train_X, train_y)
+
+# [Researcher] Provide the model and the train and test data to the TRE
+
+# [TRE] Compute the predictions on the training and test sets
+train_preds = target_model.predict_proba(train_X)
+test_preds = target_model.predict_proba(test_X)
+
+# [TRE] Define some attack parameters
+args_dict = {
+    # How many attacks to run -- in each the attack model is trained on a different
+    # subset of the data
+    'n_reps': 2,
+    # Threshold to determine significance of things
+    'p_thresh': 0.05,
+    # Filename arguments needed by the code, meaningless if run programatically
+    'in_sample_filename': None,
+    'out_sample_filename': None,
+    # Proportion of data to use as a test set for the attack model;
+    'test_prop': 0.5
+}
+
+# Convert into a namespace so that the code can access attributes with dot ('.') syntax
+args = SimpleNamespace(**args_dict)
+
+# [TRE] Call attack code
+metrics, metadata = worst_case_attack.attack(args, train_preds, test_preds)
+
+# [TRE] explore the metrics
+# For how many of the reps is the AUC p-value significant, with and without FDR correction. A
+# significant P-value means that the attack was statistically successful at predicting rows at
+# belonging in the training set
+
+print(
+    "Number of significant AUC values (raw):",
+    f"{metadata['global_metrics']['n_sig_auc_p_vals']}/{args.n_reps}"
+)
+
+print(
+    "Number of significant AUC values (FDR corrected):",
+    f"{metadata['global_metrics']['n_sig_auc_p_vals_corrected']}/{args.n_reps}"
+)
+
+# Or the number of repetitions in which the PDIF (0.1) was significant
+print(
+    "Number of significant PDIF values (proportion of 0.1), raw:",
+    f"{metadata['global_metrics']['n_sig_pdif_vals']}/{args.n_reps}"
+)
+
+print(
+    "Number of significant PDIF values (proportion of 0.1), FDR corrected:",
+    f"{metadata['global_metrics']['n_sig_pdif_vals_corrected']}/{args.n_reps}"
+)
+
+# [TRE] to compare the results obtained with those expected by chance, TRE staff can
+# also run dummy attacks
+# [TRE] Generate some fake model predictions in which there is no difference between the
+# distribution of probabilities over the training and testing data. The number of rows to
+# generate are set to match the real example
+dummy_train, dummy_test = worst_case_attack.generate_arrays(
+    len(train_preds),
+    len(test_preds)
+)
+
+# [TRE] runs the attacks on the dummy predictions
+dummy_metrics, dummy_metadata = worst_case_attack.attack(args, dummy_train, dummy_test)
+
+# [TRE] looks at the metric values to compare with those for the model
+print(
+    "(dummy) Number of significant AUC values (raw):",
+    f"{dummy_metadata['global_metrics']['n_sig_auc_p_vals']}/{args.n_reps}"
+)
+
+print(
+    "(dummy) Number of significant AUC values (FDR corrected):",
+    f"{dummy_metadata['global_metrics']['n_sig_auc_p_vals_corrected']}/{args.n_reps}"
+)
+
+# Or the number of repetitions in which the PDIF (0.1) was significant
+print(
+    "(dummy) Number of significant PDIF values (proportion of 0.1), raw:",
+    f"{dummy_metadata['global_metrics']['n_sig_pdif_vals']}/{args.n_reps}"
+)
+
+print(
+    "(dummy) Number of significant PDIF values (proportion of 0.1), FDR corrected:",
+    f"{dummy_metadata['global_metrics']['n_sig_pdif_vals_corrected']}/{args.n_reps}"
+)
+
+print("Programatic example finished")
+print("****************************")
+
+print()
+print()
+print("Command line example starting")
+print("*****************************")
+# Command line version. The same functionality as above, but the attack is run from
+# the command line rather than programatically
+
+# [Researcher] Dump the training and test predictions to .csv files
+np.savetxt("train_preds.csv", train_preds, delimiter=",")
+np.savetxt("test_preds.csv", test_preds, delimiter=",")
+
+# [TRE] Runs the attack. This would be done on the command line, here we do that with os.system
+# [TRE] First they access the help to work out which parameters they need to set
+os.system("python -m attacks.worst_case_attack run-attack --help")
+
+# [TRE] Then they run the attack
+os.system(
+    (
+        "python -m attacks.worst_case_attack run-attack "
+        "--in-sample-preds train_preds.csv "
+        "--out-of-sample-preds test_preds.csv "
+        "--n-reps 10 "
+        "--report-name example_report_risky "
+        "--dummy-reps 1 "
+        "--test-prop 0.1"
+    )
+)
+
+# [TRE] The code produces a .pdf report (example_report_risky.pdf) and a .json file that can be
+# injesetd by the shiny app
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 00000000..3cc9abbf
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,9 @@
+dictdiffer
+fpdf
+joblib
+numpy
+pandas
+pylint
+scikit_learn
+tensorflow
+tensorflow_privacy
diff --git a/safemodel/__pycache__/__init__.cpython-39.pyc b/safemodel/__pycache__/__init__.cpython-39.pyc
deleted file mode 100644
index 2b9812e7..00000000
Binary files a/safemodel/__pycache__/__init__.cpython-39.pyc and /dev/null differ
diff --git a/safemodel/__pycache__/safemodel.cpython-39.pyc b/safemodel/__pycache__/safemodel.cpython-39.pyc
deleted file mode 100644
index 21d3686d..00000000
Binary files a/safemodel/__pycache__/safemodel.cpython-39.pyc and /dev/null differ


Test patch:
diff --git a/.github/workflows/lint-test-on-push.yml b/.github/workflows/lint-test-on-push.yml
new file mode 100644
index 00000000..fe749a67
--- /dev/null
+++ b/.github/workflows/lint-test-on-push.yml
@@ -0,0 +1,28 @@
+name: pylint and unit test
+
+on: [push]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.8"]
+
+    steps:
+      - uses: actions/checkout@v3
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v4
+        with:
+          python-version: ${{ matrix.python-version }}
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
+      - name: Lint with pylint
+        run: |
+          pylint **/*.py --fail-under 7
+      - name: Test with unittest
+        run: |
+          python -m unittest discover tests
\ No newline at end of file
diff --git a/attacks/likelihood_test.py b/attacks/likelihood_test.py
new file mode 100644
index 00000000..65dd93af
--- /dev/null
+++ b/attacks/likelihood_test.py
@@ -0,0 +1,424 @@
+'''
+Likelihood testing scenario from https://arxiv.org/pdf/2112.03570.pdf
+'''
+# pylint: disable = invalid-name
+import argparse
+import json
+import logging
+import importlib
+from typing import Iterable, Any, Dict
+import numpy as np
+from scipy.stats import norm
+from sklearn.ensemble import RandomForestClassifier
+from sklearn.model_selection import train_test_split
+from  sklearn.datasets import load_breast_cancer
+
+from attacks import report
+from attacks import metrics
+
+logging.basicConfig(level=logging.INFO)
+
+N_SHADOW_MODELS = 256 # Number of shadow models that should be trained
+EPS = 1e-16 # Used to avoid numerical issues in logit function
+P_THRESH = 0.05 # default significance threshold
+
+class DummyClassifier:
+    """A Dummy Classifier to allow this code to work with get_metrics"""
+    def predict(self, test_X):
+        """Return an array of 1/0 depending on value in second column"""
+        return 1 * (test_X[:, 1] > 0.5)
+
+    def predict_proba(self, test_X):
+        """Simply return the test_X"""
+        return test_X
+
+def _logit(p: float) -> float:
+    """Standard logit function
+
+    Parameters
+    ----------
+    p: float
+        value to evaluate logit at
+
+    Returns
+    -------
+    li: float
+        logit(p)
+
+    Notes
+    -----
+    If p is close to 0 or 1, evaluating the log will result in numerical instabilities.
+    This code thresholds p at EPS and 1 - EPS where EPS defaults at 1e-16.
+    """
+    if p > 1 - EPS:
+        p = 1 - EPS
+    p = max(p, EPS)
+    li = np.log(p / (1 - p))
+    return li
+
+
+def likelihood_scenario( # pylint: disable = too-many-locals, too-many-arguments
+    shadow_clf,
+    X_target_train: Iterable[float],
+    y_target_train: Iterable[float],
+    target_train_preds: Iterable[float],
+    X_shadow_train: Iterable[float],
+    y_shadow_train: Iterable[float],
+    shadow_train_preds: Iterable[float],
+    n_shadow_models: int=N_SHADOW_MODELS
+):
+    """Implements the likelihood test, using the "offline" version
+    See p.6 (top of second column) for details
+
+    Parameters
+    ----------
+    shadow_clf: sklearn.Model
+        An sklearn classifier that will be trained to form the shadow model. All hyper-parameters
+        should have been set.
+    X_target_train: np.ndarray
+        Data that was used to train the target model
+    y_target_train: np.ndarray
+        Labels that were used to train the target model
+    target_train_preds: np.ndarray
+        Array of predictions produced by the target model on the training data
+    X_shadow_train: np.ndarray
+        Data that will be used to train the shadow models
+    y_shadow_train: np.ndarray
+        Labels that will be used to train the shadow model
+    shadow_train_preds: np.ndarray
+        Array of predictions produced by the target model on the shadow data
+    n_shadow_models: int
+        The number of shadow models to build
+
+    Returns
+    -------
+    mia_scores: np.ndarray
+        Attack probabilities of belonging to the training set or not
+    mia_labels: np.ndarray
+        True labels of belonging to the training set or not
+    mia_cls: DummyClassifier
+        A DummyClassifier that directly returns the scores for compatibility with code
+        in metrics.py
+
+    Examples
+    --------
+
+    >>> X, y = load_breast_cancer(return_X_y=True, as_frame=False)
+    >>> train_X, test_X, train_y, test_y = train_test_split(
+    >>>   X, y, test_size=0.5, stratify=y
+    >>> )
+    >>> rf = RandomForestClassifier(min_samples_leaf=1, min_samples_split=2)
+    >>> rf.fit(train_X, train_y)
+    >>> mia_test_probs, mia_test_labels, mia_clf = likelihood_scenario(
+    >>>     RandomForestClassifier(min_samples_leaf=1, min_samples_split=2, max_depth=10),
+    >>>     train_X,
+    >>>     train_y,
+    >>>     rf.predict_proba(train_X),
+    >>>     test_X,
+    >>>     test_y,
+    >>>     rf.predict_proba(test_X),
+    >>>     n_shadow_models=100
+    >>> )
+    """
+    logger = logging.getLogger("lr-scenario")
+    n_train_rows, _ = X_target_train.shape
+    n_shadow_rows, _ = X_shadow_train.shape
+    indices = np.arange(0, n_train_rows + n_shadow_rows, 1)
+
+    # Combine taregt and shadow train, from which to sample datasets
+    combined_X_train = np.vstack((
+        X_target_train,
+        X_shadow_train
+    ))
+    combined_y_train = np.hstack((
+        y_target_train,
+        y_shadow_train
+    ))
+
+    train_row_to_confidence = {i: [] for i in range(n_train_rows)}
+    shadow_row_to_confidence = {i: [] for i in range(n_shadow_rows)}
+
+    # Train N_SHADOW_MODELS shadow models
+    logger.info("Training shadow models")
+    for model_idx in range(n_shadow_models):
+        if model_idx % 10 == 0:
+            logger.info("Trained %d models", model_idx)
+        # Pick the indices to use for training this one
+        np.random.seed(model_idx) # Reproducibility
+        these_idx = np.random.choice(indices, n_train_rows, replace=False)
+        temp_X_train = combined_X_train[these_idx, :]
+        temp_y_train = combined_y_train[these_idx]
+
+        # Fit the shadow model
+        shadow_clf.set_params(random_state=model_idx)
+        shadow_clf.fit(temp_X_train, temp_y_train)
+
+        # Get the predicted probabilities on the training data
+        confidences = shadow_clf.predict_proba(X_target_train)
+        these_idx = set(these_idx)
+        for i in range(n_train_rows):
+            if i not in these_idx:
+                # If i was _not_ used for training, incorporate the logit of its confidence of
+                # being correct - TODO: should we just be taking max??
+                train_row_to_confidence[i].append(
+                    _logit(
+                        confidences[i, y_target_train[i]]
+                    )
+                )
+        # Same process for shadow data
+        shadow_confidences = shadow_clf.predict_proba(X_shadow_train)
+        for i in range(n_shadow_rows):
+            if i + n_train_rows not in these_idx:
+                shadow_row_to_confidence[i].append(
+                    _logit(
+                        shadow_confidences[i, y_shadow_train[i]]
+                    )
+                )
+
+    # Do the test described in the paper in each case
+    mia_scores = []
+    mia_labels = []
+    logger.info("Computing scores for train rows")
+    for i in range(n_train_rows):
+        true_score = _logit(target_train_preds[i, y_target_train[i]])
+        null_scores = np.array(train_row_to_confidence[i])
+        mean_null = null_scores.mean()
+        var_null = max(null_scores.var(), EPS) # var can be zero in some cases
+        prob = norm.cdf(true_score, loc=mean_null, scale=np.sqrt(var_null))
+        mia_scores.append([1 - prob, prob])
+        mia_labels.append(1)
+
+    logger.info("Computing scores for shadow rows")
+    for i in range(n_shadow_rows):
+        true_score = _logit(shadow_train_preds[i, y_shadow_train[i]])
+        null_scores = np.array(shadow_row_to_confidence[i])
+        mean_null = null_scores.mean()
+        var_null = max(null_scores.var(), EPS) # var can be zeros in some cases
+        prob = norm.cdf(true_score, loc=mean_null, scale=np.sqrt(var_null))
+        mia_scores.append([1 - prob, prob])
+        mia_labels.append(0)
+
+    mia_clf = DummyClassifier()
+    logger.info("Finished scenario")
+    return np.array(mia_scores), np.array(mia_labels), mia_clf
+
+def _run_attack(args: dict) -> None: # pylint: disable = too-many-locals
+    """Runs an attack based on the args parsed from the command line
+    """
+    logger = logging.getLogger("run-attack")
+    logger.info("Reading config from %s", args.json_file)
+    with open(args.json_file, 'r', encoding='utf-8') as f:
+        config = json.loads(f.read())
+
+    logger.info("Loading training data csv from %s", config['training_data_file'])
+    training_data = np.loadtxt(config['training_data_file'], delimiter=",")
+    train_X = training_data[:, :-1]
+    train_y = training_data[:, -1].flatten().astype(int)
+    logger.info("Loaded %d rows", len(train_X))
+
+    logger.info("Loading test data csv from %s", config['testing_data_file'])
+    test_data = np.loadtxt(config['testing_data_file'], delimiter=",")
+    test_X = test_data[:, :-1]
+    test_y = test_data[:, -1].flatten().astype(int)
+    logger.info("Loaded %d rows", len(test_X))
+
+    logger.info("Loading train predictions form %s", config['training_preds_file'])
+    train_preds = np.loadtxt(config['training_preds_file'], delimiter=",")
+    assert len(train_preds) == len(train_X)
+
+    logger.info("Loading test predictions form %s", config['testing_preds_file'])
+    test_preds = np.loadtxt(config['testing_preds_file'], delimiter=",")
+    assert len(test_preds) == len(test_X)
+
+    clf_module_name, clf_class_name = config['target_model']
+    module = importlib.import_module(clf_module_name)
+    clf_class = getattr(module, clf_class_name)
+    clf_params = config['target_hyppars']
+    clf = clf_class(**clf_params)
+    logger.info("Created model: %s", str(clf))
+    mia_test_probs, mia_test_labels, mia_clf = likelihood_scenario(
+        RandomForestClassifier(min_samples_leaf=1, min_samples_split=2, max_depth=10),
+        train_X,
+        train_y,
+        train_preds,
+        test_X,
+        test_y,
+        test_preds,
+        n_shadow_models=args.n_shadow_models
+    )
+    metadata = {
+        'input-config': args.json_file,
+        'n_shadow_models': args.n_shadow_models
+    }
+    logger.info("Computing metrics")
+    mia_metrics = metrics.get_metrics(
+        mia_clf,
+        mia_test_probs,
+        mia_test_labels
+    )
+    # Check for significance of AUC and PDIF
+    logger.info("Creating report")
+    pdif = np.exp(-mia_metrics['PDIF'])
+    metadata['PDIF_sig'] = f"Significant at p={args.p_thresh}" if pdif <= args.p_thresh \
+        else f"Not significant at p={args.p_thresh}"
+    auc_std = np.sqrt(
+        metrics.expected_auc_var(
+            0.5,
+            mia_test_labels.sum(),
+            len(mia_test_labels) - mia_test_labels.sum()
+        )
+    )
+    auc_p = 1 - norm.cdf(mia_metrics['AUC'], loc=0.5, scale=auc_std)
+    metadata['AUC_sig'] = f"Significant at p={args.p_thresh}" if auc_p <= args.p_thresh \
+        else f"Not significant at p={args.p_thresh}"
+    metadata['AUC NULL 3sd range'] = f"{0.5 - 3 * auc_std} -> {0.5 + 3 * auc_std}"
+    pdf = report.create_lr_report(metadata, mia_metrics)
+    pdf.output(args.pdf_name, 'F')
+
+
+def _example(args: Dict) -> None: # pylint: disable = too-many-locals
+    """Runs an example attack using data from sklearn
+    """
+    X, y = load_breast_cancer(return_X_y=True, as_frame=False)
+    train_X, test_X, train_y, test_y = train_test_split(
+        X, y, test_size=0.5, stratify=y
+    )
+    rf = RandomForestClassifier(min_samples_leaf=1, min_samples_split=2)
+    rf.fit(train_X, train_y)
+    mia_test_probs, mia_test_labels, mia_clf = likelihood_scenario(
+        RandomForestClassifier(min_samples_leaf=1, min_samples_split=2, max_depth=10),
+        train_X,
+        train_y,
+        rf.predict_proba(train_X),
+        test_X,
+        test_y,
+        rf.predict_proba(test_X),
+        n_shadow_models=args.n_shadow_models
+    )
+    metadata = {
+        'input-data': "Breast Cancer",
+        'n_shadow_models': args.n_shadow_models
+    }
+    mia_metrics = metrics.get_metrics(
+        mia_clf,
+        mia_test_probs,
+        mia_test_labels
+    )
+    # Check for significance of AUC and PDIF
+    pdif = np.exp(-mia_metrics['PDIF'])
+    metadata['PDIF_sig'] = f"Significant at p={args.p_thresh}" if pdif <= args.p_thresh \
+        else f"Not significant at p={args.p_thresh}"
+    auc_std = np.sqrt(
+        metrics.expected_auc_var(
+            0.5,
+            mia_test_labels.sum(),
+            len(mia_test_labels) - mia_test_labels.sum()
+        )
+    )
+    auc_p = 1 - norm.cdf(mia_metrics['AUC'], loc=0.5, scale=auc_std)
+    metadata['AUC_sig'] = f"Significant at p={args.p_thresh}" if auc_p <= args.p_thresh \
+        else f"Not significant at p={args.p_thresh}"
+    metadata['AUC NULL 3sd range'] = f"{0.5 - 3 * auc_std} -> {0.5 + 3 * auc_std}"
+    pdf = report.create_lr_report(metadata, mia_metrics)
+    pdf.output(args.pdf_name, 'F')
+
+
+def _setup_example_data(_: Any):
+    """Method to create example data and save (including config). Intended to allow users
+    to see how they would need to setup their own data.
+    """
+    X, y = load_breast_cancer(return_X_y=True)
+    train_X, test_X, train_y, test_y = train_test_split(
+        X, y, test_size=0.5, stratify=y
+    )
+    rf = RandomForestClassifier(min_samples_split=2, min_samples_leaf=1)
+    rf.fit(train_X, train_y)
+    train_data = np.hstack((train_X, train_y[:, None]))
+    np.savetxt("train_data.csv", train_data, delimiter=",")
+
+    test_data = np.hstack((test_X, test_y[:, None]))
+    np.savetxt("test_data.csv", test_data, delimiter=",")
+
+    train_preds = rf.predict_proba(train_X)
+    test_preds = rf.predict_proba(test_X)
+    np.savetxt("train_preds.csv", train_preds, delimiter=",")
+    np.savetxt("test_preds.csv", test_preds, delimiter=",")
+
+
+    config = {
+        "training_data_file": "train_data.csv",
+        "testing_data_file": "test_data.csv",
+        "training_preds_file": "train_preds.csv",
+        "testing_preds_file": "test_preds.csv",
+        "target_model": ["sklearn.ensemble", "RandomForestClassifier"],
+        "target_hyppars": {
+            "min_samples_split": 2,
+            "min_samples_leaf": 1
+        }
+    }
+
+    with open('config.json', 'w', encoding='utf-8') as f:
+        f.write(json.dumps(config))
+
+def main():
+    """Main method to parse args and invoke relevant code
+    """
+    parser = argparse.ArgumentParser(add_help=False)
+    parser.add_argument('-s', '--n-shadow-models',
+        type=int,
+        required=False,
+        default=N_SHADOW_MODELS,
+        action='store',
+        dest='n_shadow_models',
+        help=(
+            'The number of shadow models to train (default = %(default)d)'
+        )
+    )
+    parser.add_argument('--report-name',
+        type=str,
+        action="store",
+        dest="pdf_name",
+        required=False,
+        default="lr_report.pdf",
+        help=(
+            'Output name for the .pdf report. Default = %(default)s'
+        )
+    )
+    parser.add_argument('-p', '--p-thresh',
+        type=float,
+        action='store',
+        dest='p_thresh',
+        required=False,
+        default=P_THRESH,
+        help=(
+            'Significance threshold for p-value comparisons. Default = %(default)f'
+        )
+    )
+    subparsers = parser.add_subparsers()
+    example_parser = subparsers.add_parser('run-example', parents=[parser])
+    example_parser.set_defaults(func=_example)
+
+    attack_parser = subparsers.add_parser('run-attack', parents=[parser])
+    attack_parser.add_argument('-j', '--json-file',
+        action="store",
+        required=True,
+        dest="json_file",
+        type=str,
+        help=(
+            'Name of the .json file containing details for the run. Default = %(default)s'
+        )
+    )
+    attack_parser.set_defaults(func=_run_attack)
+
+    example_data_parser = subparsers.add_parser('setup-example-data')
+    example_data_parser.set_defaults(func=_setup_example_data)
+
+    args = parser.parse_args()
+    try:
+        args.func(args)
+    except AttributeError:
+        print("Invalid command. Try --help to get more details")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_data_interface.py b/tests/test_data_interface.py
new file mode 100644
index 00000000..33f2626e
--- /dev/null
+++ b/tests/test_data_interface.py
@@ -0,0 +1,47 @@
+'''
+Test the data interface code
+'''
+
+import unittest
+import pandas as pd
+from data_preprocessing.data_interface import get_data_sklearn, UnknownDataset  # pylint: disable = import-error
+
+
+class TestLoaders(unittest.TestCase):
+    '''
+    Test the data loaders
+    '''
+    def test_iris(self):
+        '''Nursery data'''
+        feature_df, target_df = get_data_sklearn("iris")
+        self.assertIsInstance(feature_df, pd.DataFrame)
+        self.assertIsInstance(target_df, pd.DataFrame)
+
+    def test_unknown(self):
+        '''
+        Test that a nonsense string raises the correct exception
+        '''
+        with self.assertRaises(UnknownDataset):
+            _, _ = get_data_sklearn("NONSENSE")
+
+    def test_standard(self):
+        '''
+        Test that standardisation creates standard features
+        '''
+        feature_df, _ = get_data_sklearn("standard iris")
+        for column in feature_df.columns:
+            temp = feature_df[column].mean()
+            self.assertAlmostEqual(temp, 0.0)
+            temp = feature_df[column].std()
+            self.assertAlmostEqual(temp, 1.0)
+
+    def test_minmax(self):
+        '''
+        Test the minmax scaling
+        '''
+        feature_df, _ = get_data_sklearn("minmax iris")
+        for column in feature_df.columns:
+            temp = feature_df[column].min()
+            self.assertAlmostEqual(temp, 0.0)
+            temp = feature_df[column].max()
+            self.assertAlmostEqual(temp, 1.0)
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
new file mode 100644
index 00000000..e9b8caf6
--- /dev/null
+++ b/tests/test_metrics.py
@@ -0,0 +1,113 @@
+'''
+Test the metrics
+'''
+
+import unittest
+import numpy as np
+from attacks.metrics import get_metrics, tpr_at_fpr
+from attacks.mia_extremecase import min_max_disc
+
+# pylint: disable = invalid-name
+
+PREDICTED_CLASS = np.array([0, 1, 0, 0, 1, 1])
+TRUE_CLASS = np.array([0, 0, 0, 1, 1, 1])
+PREDICTED_PROBS = np.array([
+    [0.9, 0.1],
+    [0.4, 0.6],
+    [0.8, 0.2],
+    [0.55, 0.45],
+    [0.1, 0.9],
+    [0.01, 0.99]
+])
+
+class DummyClassifier:
+    '''
+    Mocks the predict and predict_proba methods
+    '''
+    def predict(self, _):
+        '''
+        return dummy predictions
+        '''
+        return PREDICTED_CLASS
+    def predict_proba(self, _):
+        '''
+        return dummy predicted probabilities
+        '''
+        return PREDICTED_PROBS
+
+class TestMetrics(unittest.TestCase):
+    '''
+    Test the metrics with some dummy predictions
+    '''
+    def test_metrics(self):
+        '''
+        Test each individual metric with dummy data
+        '''
+        clf = DummyClassifier()
+        testX = []
+        testy = TRUE_CLASS
+        metrics = get_metrics(clf, testX, testy, permute_rows=False)
+        self.assertAlmostEqual(metrics['TPR'], 2 / 3)
+        self.assertAlmostEqual(metrics['FPR'], 1 / 3)
+        self.assertAlmostEqual(metrics['FAR'], 1 / 3)
+        self.assertAlmostEqual(metrics['TNR'], 2 / 3)
+        self.assertAlmostEqual(metrics['PPV'], 2 / 3)
+        self.assertAlmostEqual(metrics['NPV'], 2 / 3)
+        self.assertAlmostEqual(metrics['FNR'], 1 / 3)
+        self.assertAlmostEqual(metrics['ACC'], 4 / 6)
+        self.assertAlmostEqual(metrics['F1score'], (8 / 9) / (2 / 3 + 2 / 3))
+        self.assertAlmostEqual(metrics['Advantage'], 1 / 3)
+        self.assertAlmostEqual(metrics['AUC'], 8 / 9)
+
+class TestFPRatTPR(unittest.TestCase):
+    '''
+    Test code that computes TPR at fixed FPR
+    '''
+    def test_tpr(self):
+        '''Test tpr at fpr'''
+        y_true = TRUE_CLASS
+        y_score = PREDICTED_PROBS[:, 1]
+
+        tpr = tpr_at_fpr(y_true, y_score, fpr=0)
+        self.assertAlmostEqual(tpr, 2/3)
+        tpr = tpr_at_fpr(y_true, y_score, fpr=0.001)
+        self.assertAlmostEqual(tpr, 2/3)
+        tpr = tpr_at_fpr(y_true, y_score, fpr=0.1)
+        self.assertAlmostEqual(tpr, 2/3)
+        tpr = tpr_at_fpr(y_true, y_score, fpr=0.4)
+        self.assertAlmostEqual(tpr, 1)
+        tpr = tpr_at_fpr(y_true, y_score, fpr=1.0)
+        self.assertAlmostEqual(tpr, 1)
+
+
+class TestExtreme(unittest.TestCase):
+    '''
+    Test the extreme metrics
+    '''
+    def test_extreme_default(self):
+        '''
+        Tets with the dummy data
+        '''
+        pred_probs = DummyClassifier().predict_proba(None)[:, 1]
+        maxd, mind, mmd, _ = min_max_disc(TRUE_CLASS, pred_probs)
+
+        # 10% of 6 is 1 so:
+        # maxd should be 1 (the highest one is predicted as1)
+        # mind should be 0 (the lowest one is not predicted as1)
+        self.assertAlmostEqual(maxd, 1.0)
+        self.assertAlmostEqual(mind, 0.0)
+        self.assertAlmostEqual(mmd, 1.0)
+
+    def test_extreme_higer_prop(self):
+        '''
+        Tets with the dummy data but increase proportion to 0.5
+        '''
+        pred_probs = DummyClassifier().predict_proba(None)[:, 1]
+        maxd, mind, mmd, _ = min_max_disc(TRUE_CLASS, pred_probs, x_prop=0.5)
+
+        # 10% of 6 is 1 so:
+        # maxd should be 1 (the highest one is predicted as1)
+        # mind should be 0 (the lowest one is not predicted as1)
+        self.assertAlmostEqual(maxd, 2 / 3)
+        self.assertAlmostEqual(mind, 1 / 3)
+        self.assertAlmostEqual(mmd, 1 / 3)