diff --git a/dlomix/constants.py b/dlomix/constants.py index deac0759..558b547b 100644 --- a/dlomix/constants.py +++ b/dlomix/constants.py @@ -40,3 +40,49 @@ "W": 19, "Y": 20, } + +ALPHABET_MOD = { + "A": 1, + "C": 2, + "D": 3, + "E": 4, + "F": 5, + "G": 6, + "H": 7, + "I": 8, + "K": 9, + "L": 10, + "M": 11, + "N": 12, + "P": 13, + "Q": 14, + "R": 15, + "S": 16, + "T": 17, + "V": 18, + "W": 19, + "Y": 20, + "^": 21, + "}": 22, +} + +MCD_pipeline_parameters = { + # params for full dataset with modified alphabet + #"model_params": {"seq_length": 30, "vocab_dict": ALPHABET_MOD.copy()}, + #"data_params": {"seq_length": 30, 'sequence_col': 'modified_sequence_single_letter', 'target_col': 'median'}, + "model_params": {"seq_length": 30, "vocab_dict": ALPHABET_UNMOD.copy()}, + "data_params": {"seq_length": 30,}, + #TODO upload model in dlomix and change path + "base_model_path": "../data/models/rtmodel_prosit_epoch20", + "test_set_path": "https://raw.githubusercontent.com/wilhelm-lab/dlomix/develop/example_dataset/proteomTools_test.csv" +} + + +QR_pipeline_parameters = { + "model_params": {"seq_length": 30, "vocab_dict": ALPHABET_UNMOD.copy()}, + "data_params": {"seq_length": 30,}, + #specify path of desired model here + "base_model_path": "../../data/models/rtmodel_prosit_epoch20", + "train_val_path": 'https://raw.githubusercontent.com/wilhelm-lab/dlomix/develop/example_dataset/proteomTools_train_val.csv', + "test_path": 'https://raw.githubusercontent.com/wilhelm-lab/dlomix/develop/example_dataset/proteomTools_test.csv', +} \ No newline at end of file diff --git a/dlomix/eval/__init__.py b/dlomix/eval/__init__.py index 4df20838..b5f08758 100644 --- a/dlomix/eval/__init__.py +++ b/dlomix/eval/__init__.py @@ -1,3 +1,5 @@ from .rt_eval import TimeDeltaMetric +from .interval_conformal import IntervalSize, AbsoluteIntervalSize, RelativeCentralDistance, IntervalConformalScore, IntervalConformalQuantile +from .scalar_conformal import ScalarConformalScore, ScalarConformalQuantile -__all__ = [TimeDeltaMetric] +__all__ = [TimeDeltaMetric, IntervalSize, AbsoluteIntervalSize, RelativeCentralDistance, IntervalConformalScore, IntervalConformalQuantile, ScalarConformalScore, ScalarConformalQuantile] diff --git a/dlomix/eval/interval_conformal.py b/dlomix/eval/interval_conformal.py new file mode 100644 index 00000000..a33fcca9 --- /dev/null +++ b/dlomix/eval/interval_conformal.py @@ -0,0 +1,104 @@ +import tensorflow as tf + +class IntervalSize(tf.keras.losses.Loss): + ''' + Size of the prediction interval. + ''' + def __init__(self, name="interval_size", **kwargs): + super().__init__(name=name, **kwargs) + + def call(self, y_true, y_pred): + interval_size = tf.subtract(y_pred[:,1], y_pred[:,0]) + return interval_size + + def get_config(self): + config = {} + base_config = super().get_config() + return {**base_config, **config} + + +class AbsoluteIntervalSize(tf.keras.losses.Loss): + ''' + Absolute size of the prediction interval. + ''' + def __init__(self, name="abs_interval_size", **kwargs): + super().__init__(name=name, **kwargs) + + def call(self, y_true, y_pred): + abs_interval_size = tf.math.abs(tf.subtract(y_pred[:,1], y_pred[:,0])) + return abs_interval_size + + def get_config(self): + config = {} + base_config = super().get_config() + return {**base_config, **config} + + +class RelativeCentralDistance(tf.keras.losses.Loss): + ''' + Distance of the true value from the center of the prediction intverval, + divided by half of the the prediction interval size. + The result is 0 if true value is at center of the prediction interval, + 1 if at (symm.) interval boundary and >1 else. + ''' + def __init__(self, name="relative_central_distance", **kwargs): + super().__init__(name=name, **kwargs) + + def call(self, y_true, y_pred): + # absolute distance between the two interval boundaries + # (if upper boudary is below lower, this is nevertheless positive) + interval_size = tf.math.abs(tf.subtract(y_pred[:,1], y_pred[:,0])) + # absolute distance of the true value from the inverval center/mean + central_dist = tf.math.abs(tf.subtract(y_true[:,0], tf.math.reduce_mean(y_pred, axis=-1))) + # divide by half of interval size + res = tf.divide(central_dist, tf.divide(interval_size, 2.)) + indices = tf.where(tf.math.is_inf(res)) + return tf.tensor_scatter_nd_update( + res, + indices, + tf.ones((tf.shape(indices)[0])) * 0. + ) + + def get_config(self): + config = {} + base_config = super().get_config() + return {**base_config, **config} + + +class IntervalConformalScore(tf.keras.losses.Loss): + ''' + Computes conformal scores for prediction intervals pred_intervals and true values y_true. + ''' + def __init__(self, name="interval_conformal_score", **kwargs): + super().__init__(name=name, **kwargs) + + def call(self, y_true, pred_intervals): + return tf.reduce_max(tf.stack([tf.subtract(pred_intervals, y_true)[:,0], + -tf.subtract(pred_intervals, y_true)[:,1]], 1), -1) + + def get_config(self): + config = {} + base_config = super().get_config() + return {**base_config, **config} + + +class IntervalConformalQuantile(tf.keras.losses.Loss): + ''' + Computes the conformal quantile based on the distribution of conformal scores + for prediction intervals pred_intervals and true values y_true. + ''' + def __init__(self, alpha=0.1, name="interval_conformal_quantile", **kwargs): + super().__init__(name=name, **kwargs) + self.alpha = alpha + + def call(self, y_true, pred_intervals): + scores = tf.reduce_max(tf.stack([tf.subtract(pred_intervals, y_true)[:,0], -tf.subtract(pred_intervals, y_true)[:,1]], 1), -1) + n = tf.cast(tf.shape(y_true)[0], tf.float32) #without casting to float, next line throws an error + q = tf.math.ceil((n + 1.) * (1. - self.alpha)) / n + tfp_quantile = tf.sort(scores, axis=-1, direction='ASCENDING', name=None)[tf.math.minimum(int(q * n), int(n-1))] + return tfp_quantile + + def get_config(self): + config = {'alpha' : self.alpha} + base_config = super().get_config() + return {**base_config, **config} \ No newline at end of file diff --git a/dlomix/eval/scalar_conformal.py b/dlomix/eval/scalar_conformal.py new file mode 100644 index 00000000..72de4168 --- /dev/null +++ b/dlomix/eval/scalar_conformal.py @@ -0,0 +1,38 @@ +import tensorflow as tf + +class ScalarConformalScore(tf.keras.losses.Loss): + ''' + Computes conformal scores for predicted scalar error estimates pred_err and true values y_true. + ''' + def __init__(self, name="scalar_conformal_scores", **kwargs): + super().__init__(name=name, **kwargs) + + def call(self, y_true, pred_err): + scores = tf.divide(tf.math.abs(y_true - pred_err[:,0]), pred_err[:,1]) + return scores + + def get_config(self): + config = {} + base_config = super().get_config() + return {**base_config, **config} + +class ScalarConformalQuantile(tf.keras.losses.Loss): + ''' + Computes the conformal quantile based on the distribution of conformal scores + for scalar error estimates pred_err and true values y_true. + ''' + def __init__(self, alpha=0.1, name="scalar_conformal_quantile", **kwargs): + super().__init__(name=name, **kwargs) + self.alpha = alpha + + def call(self, y_true, pred_err): + scores = tf.divide(tf.math.abs(y_true - pred_err[:,0]), pred_err[:,1]) + n = tf.cast(tf.shape(y_true)[0], tf.float32) #without casting to float, next line throws an error + q = tf.math.ceil((n + 1.)*(1. - self.alpha)) / n + tfp_quantile = tf.sort(scores, axis=-1, direction='ASCENDING', name=None)[int(q * n)] + return tfp_quantile + + def get_config(self): + config = {'alpha' : self.alpha} + base_config = super().get_config() + return {**base_config, **config} \ No newline at end of file diff --git a/dlomix/losses/__init__.py b/dlomix/losses/__init__.py index 81bafa1b..38b84009 100644 --- a/dlomix/losses/__init__.py +++ b/dlomix/losses/__init__.py @@ -1,3 +1,4 @@ from .intensity import masked_spectral_distance, masked_pearson_correlation_distance +from .quantile import QuantileLoss -__all__ = [masked_spectral_distance, masked_pearson_correlation_distance] +__all__ = [masked_spectral_distance, masked_pearson_correlation_distance, QuantileLoss] diff --git a/dlomix/losses/quantile.py b/dlomix/losses/quantile.py new file mode 100644 index 00000000..f3ed451f --- /dev/null +++ b/dlomix/losses/quantile.py @@ -0,0 +1,20 @@ +import tensorflow as tf + +class QuantileLoss(tf.keras.losses.Loss): + ''' + Quantile loss (pinball loss) + ''' + def __init__(self, quantile=tf.constant([[0.1, 0.9]]), name="quantile_loss", **kwargs): + super().__init__(name=name, **kwargs) + self.quantile = quantile + + def call(self, y_true, y_pred): + err = tf.subtract(y_true, y_pred) + return tf.reduce_sum(tf.maximum(self.quantile * err, (self.quantile - 1) * err), axis=-1) + + def get_config(self): + config = { + 'quantile': self.quantile + } + base_config = super().get_config() + return {**base_config, **config} \ No newline at end of file diff --git a/dlomix/models/__init__.py b/dlomix/models/__init__.py index 5f6c4df6..aef129f3 100644 --- a/dlomix/models/__init__.py +++ b/dlomix/models/__init__.py @@ -7,4 +7,4 @@ "PrositRetentionTimePredictor", "DeepLCRetentionTimePredictor", "PrositIntensityPredictor", -] +] \ No newline at end of file diff --git a/dlomix/pipelines/MCDpipeline.py b/dlomix/pipelines/MCDpipeline.py new file mode 100644 index 00000000..283ac623 --- /dev/null +++ b/dlomix/pipelines/MCDpipeline.py @@ -0,0 +1,101 @@ +import pickle +import numpy as np +import tensorflow as tf +from scipy.stats import ks_2samp + +from dlomix.constants import MCD_pipeline_parameters +from dlomix.data.RetentionTimeDataset import RetentionTimeDataset +from dlomix.models import PrositRetentionTimePredictor +from dlomix.eval.scalar_conformal import ScalarConformalScore, ScalarConformalQuantile +from dlomix.losses import QuantileLoss +from dlomix.reports import UncertaintyReport + + +class MCDPipeline: + def __init__(self, alpha=0.1): + self.base_model = None + self.test_dataset = None + self.alpha = alpha + self.res = None + self.label = None + + self._build_base_model() + + def _build_base_model(self): + + self.base_model = PrositRetentionTimePredictor(**MCD_pipeline_parameters["model_params"]) + self.base_model.load_weights(MCD_pipeline_parameters["base_model_path"]).expect_partial() + + def _predict_with_dropout(self, n): + + predictions = [] + for i in range(n): + res = np.concatenate([self.base_model(batch[0], training=True).numpy() for batch in list(self.test_dataset.test_data)]) + predictions.append(res) + return np.column_stack(predictions) + + def load_data(self, data=MCD_pipeline_parameters["test_set_path"], batchsize=32): + + if not (isinstance(data, str) or isinstance(data, np.ndarray)): + raise ValueError( + "Dataset should be provided either as a numpy array or a string pointing to a file." + ) + + self.test_dataset = RetentionTimeDataset( + data_source=data, + **MCD_pipeline_parameters["data_params"], + batch_size=batchsize, + test=True) + + def predict(self, n=10): + + res = {} + self.label = f"PROSIT_MCD_n={n}" + res[self.label] = {} + print("model :", self.label, ", n =", n) + pred = self._predict_with_dropout(n=n) + res[self.label]['data'] = np.array((pred.mean(axis=1), pred.std(axis=1))) + self.res = res + + #return res + + def report(self): + + test_targets = self.test_dataset.get_split_targets(split="test") + avgs, stds = self.res[self.label]['data'][0], self.res[self.label]['data'][1] + print(f'#### {self.label} ####') + + conf_scores = ScalarConformalScore(reduction='none')(test_targets, self.res[self.label]['data'].T).numpy() + conf_quantile = ScalarConformalQuantile()(test_targets, self.res[self.label]['data'].T).numpy() + + print(f"alpha = {self.alpha}, conformal quantile: {conf_quantile:.2f}") + + intervals = np.array([avgs - stds * conf_quantile, avgs + stds * conf_quantile]).T + interval_sizes = intervals[:,1] - intervals[:,0] + within = (test_targets >= intervals[:,0]) & (test_targets <= intervals[:,1]) + + UncertaintyReport.plot_conformal_scores(conf_scores, quantile=conf_quantile) + UncertaintyReport.plot_predictions_with_intervals(test_targets, avgs, intervals) + UncertaintyReport.plot_conformalized_interval_size(interval_sizes) + + pvalue = ks_2samp(interval_sizes[within], interval_sizes[~within]).pvalue # prob. for Null: distr are identical + print(f"p = {pvalue:.5f} : {'Reject' if pvalue < 0.01 else 'Accept'} Null Hypothesis (Distr. identical)") + + UncertaintyReport.plot_conformalized_interval_size_PDFs(interval_sizes, within, pvalue) + + def save_results(self, path="MCD_results.pkl"): + + with open(path, 'wb') as f: + pickle.dump(self.res, f) + + def run_full_pipeline(self, n=10, save=True, save_path='MCD_results.pkl'): + + self.load_data() + self.predict(n=n) + self.report() + if save: + self.save_results(path=save_path) + + + + diff --git a/dlomix/pipelines/QRpipeline.py b/dlomix/pipelines/QRpipeline.py new file mode 100644 index 00000000..1daf1e16 --- /dev/null +++ b/dlomix/pipelines/QRpipeline.py @@ -0,0 +1,149 @@ +import pickle +import os +import numpy as np +import tensorflow as tf +from scipy.stats import ks_2samp + + +from dlomix.constants import QR_pipeline_parameters +from dlomix.data.RetentionTimeDataset import RetentionTimeDataset +from dlomix.models import PrositRetentionTimePredictor +from dlomix.losses import QuantileLoss +from dlomix.reports import UncertaintyReport + +from dlomix.eval.interval_conformal import IntervalSize, AbsoluteIntervalSize, RelativeCentralDistance, \ + IntervalConformalScore, IntervalConformalQuantile + + +class QRpipeline(): + + def __init__(self, frozen, alpha=0.1): + self.model = None + self.train_val_dataset = None + self.test_dataset = None + self.init_predictions = None + self.qr_predictions = None + + self.frozen = frozen + self.alpha = alpha + self.label = "-" + str(frozen) + " layer(s)" if frozen else "complete" + + self._load_model() + + def _load_model(self): + self.model = PrositRetentionTimePredictor(**QR_pipeline_parameters["model_params"]) + self.model.load_weights(QR_pipeline_parameters["base_model_path"]).expect_partial() + + def load_data(self, + train_val_source=QR_pipeline_parameters["train_val_path"], + test_source=QR_pipeline_parameters["test_path"], + val_ratio=0.2, batchsize=32): + + if not (isinstance(train_val_source, str) or isinstance(train_val_source, np.ndarray)): + raise ValueError( + "train_val dataset should be provided either as a numpy array or a string pointing to a file." + ) + if not (isinstance(test_source, str) or isinstance(test_source, np.ndarray)): + raise ValueError( + "train_val dataset should be provided either as a numpy array or a string pointing to a file." + ) + + self.train_val_dataset = RetentionTimeDataset( + data_source=train_val_source, + **QR_pipeline_parameters["data_params"], + batch_size=batchsize, + val_ratio=val_ratio, + test=False) + + self.test_dataset = RetentionTimeDataset( + data_source=test_source, + **QR_pipeline_parameters["data_params"], + batch_size=batchsize, + test=True) + + def predict(self, qr_prediction): + + predictions = self.model.predict(self.test_dataset.test_data) + #predictions = self.test_dataset.denormalize_targets(predictions) + + if qr_prediction: + self.qr_predictions = predictions + else: + predictions = predictions.ravel() + self.init_predictions = predictions + + def adapt_model_for_qr(self): + + if not self.frozen: + #change weights to random since whole model is unfrozen + self.model = PrositRetentionTimePredictor(**QR_pipeline_parameters["model_params"]) + else: + #"freeze" all layers up until the specified val + for layer in self.model.layers[:-self.frozen]: + layer.trainable = False + self.model.output_layer = tf.keras.Sequential([tf.keras.layers.Dense(2)]) + + def train_model_on_ql(self, epochs=10, batch_size=32, checkpoint_dir='checkpoints'): + + train_steps = batch_size * len(self.train_val_dataset.train_data) * epochs + lr_fn = tf.optimizers.schedules.PolynomialDecay(1e-3, train_steps, 1e-6, 2) + opt = tf.optimizers.Adam(lr_fn) + + self.model.compile(optimizer=opt, + loss=QuantileLoss(tf.constant([[self.alpha, 1-self.alpha]])), + metrics=['mean_absolute_error', + RelativeCentralDistance(), + IntervalConformalQuantile(alpha=self.alpha), + AbsoluteIntervalSize()]) + #self.model.build(input_shape=(None, 30)) + + model_checkpoint_callback = tf.keras.callbacks.ModelCheckpoint(filepath=os.path.join(checkpoint_dir, self.label, 'best.hdf5'), + save_weights_only=True, + monitor='val_loss', + mode='min', + save_best_only=True) + + history = self.model.fit(self.train_val_dataset.train_data, + validation_data=self.train_val_dataset.val_data, + callbacks=[model_checkpoint_callback], + epochs=epochs) + #return history + + + def report(self): + + test_targets = self.test_dataset.get_split_targets(split="test") + + conf_scores = IntervalConformalScore(reduction='none')(np.expand_dims(test_targets,-1), self.qr_predictions).numpy() + conf_quantile = IntervalConformalQuantile(alpha=0.1, reduction='none')(np.expand_dims(test_targets,-1), self.qr_predictions).numpy() + + intervals = self.qr_predictions.copy() + intervals[:,0] -= conf_quantile + intervals[:,1] += conf_quantile + interval_sizes = intervals[:,1] - intervals[:,0] + within = (test_targets >= intervals[:,0]) & (test_targets <= intervals[:,1]) + + UncertaintyReport.plot_conformal_scores(conf_scores, quantile=conf_quantile) + UncertaintyReport.plot_predictions_with_intervals(test_targets, self.init_predictions, intervals) + UncertaintyReport.plot_conformalized_interval_size(interval_sizes) + + pvalue = ks_2samp(interval_sizes[within], interval_sizes[~within]).pvalue # prob. for Null: distr are identical + print('####', self.label, '####\nconformal quantile:', conf_quantile, '\n') + + UncertaintyReport.plot_conformalized_interval_size_PDFs(interval_sizes, within, pvalue) + + + def save_data(self, path): + with open(path, 'wb') as f: + pickle.dump(self.qr_predictions, f) + + def run_full_pipeline(self, epochs=10, save=True, save_path='QR_preds.pkl'): + + self.load_data() + self.predict(qr_prediction=False) + self.adapt_model_for_qr() + self.train_model_on_ql(epochs=epochs) + self.predict(qr_prediction=True) + self.report() + if save: + self.save_data(path=save_path) \ No newline at end of file diff --git a/dlomix/pipelines/__init__.py b/dlomix/pipelines/__init__.py index 8d4a3a85..94244f17 100644 --- a/dlomix/pipelines/__init__.py +++ b/dlomix/pipelines/__init__.py @@ -1,3 +1,8 @@ from .pipeline import RetentionTimePipeline +from .MCDpipeline import MCDPipeline +from .QRpipeline import QRpipeline -__all__ = [RetentionTimePipeline] +__all__ = [RetentionTimePipeline, + MCDPipeline, + QRpipeline + ] diff --git a/dlomix/reports/UncertaintyReport.py b/dlomix/reports/UncertaintyReport.py new file mode 100644 index 00000000..7cd8b1d4 --- /dev/null +++ b/dlomix/reports/UncertaintyReport.py @@ -0,0 +1,76 @@ +import numpy as np +import pandas as pd +import seaborn as sns +import matplotlib.pyplot as plt + +class UncertaintyReport(): + + @staticmethod + def plot_conformal_scores(conf_scores, quantile=None, xlim=6, bins=100, range=None, figsize=(6,4)): + fig,ax = plt.subplots(figsize=figsize) + xmax = conf_scores.std() * xlim + if range is not None: + ax.hist(conf_scores[conf_scores < xmax], bins=bins, range=range) + else: + ax.hist(conf_scores[conf_scores < xmax], bins=bins) + ax.set_xlim(range) + if quantile: + ax.axvline(quantile, color='red', alpha=0.5) + ax.set(xlabel='conformal score', ylabel='# per bin', title=f'conformal scores (<{xlim} std.dev. from mean)') + plt.show() + + @staticmethod + def plot_predictions_with_intervals(test_targets, test_estimates, intervals, label=None, figsize=(6,4)): + fig,ax = plt.subplots(figsize=figsize) + if label: + ax.set_title(label) + p = test_targets.argsort() + ax.plot(test_targets[p], intervals[p, 0], alpha=0.1) + ax.plot(test_targets[p], intervals[p, 1], alpha=0.1) + ax.scatter(test_targets[p], test_estimates[p], s=1, alpha=0.1) + ax.plot((test_targets.min(),test_targets.max()), (test_targets.min(),test_targets.max()), alpha=0.3, color='black', linestyle='--') + ax.set(title=f'{label+" " if label else ""}predictions with error intervals, sorted by RT', + xlabel='true retention time', + ylabel='predicted retention time') + plt.show() + + @staticmethod + def plot_interval_size_dist(intervals, xlim=6., bins=100, figsize=(6,4)): + fig,ax = plt.subplots(figsize=figsize) + sizes = intervals[:,1] - intervals[:,0] + xmin, xmax = sizes.mean() - sizes.std() * xlim, sizes.mean() + sizes.std() * xlim + ax.hist(sizes[(sizes >= xmin) & (sizes <= xmax)], bins=bins) + ax.set(xlabel='interval size', ylabel='# per bin') + plt.show() + + @staticmethod + def plot_conformalized_interval_size(interval_sizes, bins=200, figsize=(6,4), range=None, ylim=None): + # plot histogram of conformalized interval size + fig,ax = plt.subplots(figsize=figsize) + if range is not None: + ax.hist(interval_sizes, bins=bins, range=range) + else: + ax.hist(interval_sizes, bins=bins, range=range) + ax.set(xlim=range) + if ylim is not None: + ax.set(ylim=ylim) + ax.set(xlabel='interval size', ylabel='# per bin', title=f'conformalized interval sizes') + plt.show() + + @staticmethod + def plot_conformalized_interval_size_PDFs(interval_sizes, within, pvalue, bins=100, figsize=(6,4), ylim=None, range=None): + # plot PDFs of conformalized interval size depending on target inside/outside conf. intervals + fig,ax = plt.subplots(figsize=figsize) + if range is not None: + ax.hist(interval_sizes[within], bins=bins, range=range, histtype='step', density=True, color='C0', label='inside interval') + ax.hist(interval_sizes[~within], bins=bins, range=range, histtype='step', density=True, color='C1', label='outside interval') + ax.set(xlim=range) + else: + ax.hist(interval_sizes[within], bins=bins, histtype='step', density=True, color='C0', label='inside interval') + ax.hist(interval_sizes[~within], bins=bins, histtype='step', density=True, color='C1', label='outside interval') + if ylim is not None: + ax.set(ylim=ylim) + ax.set(xlabel='interval size', ylabel='fraction per bin', title=f'PDF of conformalized interval sizes') + ax.text(0.98, 0.8, f"identical: p = {pvalue:.5f}", transform=ax.transAxes, ha='right', va='top') + ax.legend() + plt.show() \ No newline at end of file diff --git a/dlomix/reports/__init__.py b/dlomix/reports/__init__.py index ab165911..b966e33d 100644 --- a/dlomix/reports/__init__.py +++ b/dlomix/reports/__init__.py @@ -1,6 +1,8 @@ from .IntensityReport import IntensityReport from .RetentionTimeReport import RetentionTimeReport +from .UncertaintyReport import UncertaintyReport __all__ = ["RetentionTimeReport", "IntensityReport", + "UncertaintyReport", ] diff --git a/dlomix/utils.py b/dlomix/utils.py index 76e183a3..9ae9e728 100644 --- a/dlomix/utils.py +++ b/dlomix/utils.py @@ -1,5 +1,8 @@ import pickle import numpy as np +import tensorflow as tf +import os +import random def save_obj(obj, name): @@ -14,3 +17,18 @@ def load_obj(name): def convert_nested_list_to_numpy_array(nested_list, dtype=np.float32): return np.array([np.array(x, dtype=dtype) for x in nested_list]) + + +def set_global_seed(seed=42): + # https://stackoverflow.com/questions/36288235/how-to-get-stable-results-with-tensorflow-setting-random-seed + os.environ['PYTHONHASHSEED'] = str(seed) + random.seed(seed) + tf.random.set_seed(seed) + np.random.seed(seed) + + os.environ['TF_DETERMINISTIC_OPS'] = '1' + os.environ['TF_CUDNN_DETERMINISTIC'] = '1' + + tf.config.threading.set_inter_op_parallelism_threads(1) + tf.config.threading.set_intra_op_parallelism_threads(1) + \ No newline at end of file diff --git a/notebooks/Example_Pipelines.ipynb b/notebooks/Example_Pipelines.ipynb new file mode 100644 index 00000000..8171f986 --- /dev/null +++ b/notebooks/Example_Pipelines.ipynb @@ -0,0 +1,137 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2023-09-21 15:46:19.002443: I tensorflow/tsl/cuda/cudart_stub.cc:28] Could not find cuda drivers on your machine, GPU will not be used.\n", + "2023-09-21 15:46:19.051031: I tensorflow/tsl/cuda/cudart_stub.cc:28] Could not find cuda drivers on your machine, GPU will not be used.\n", + "2023-09-21 15:46:19.051807: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.\n", + "To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.\n", + "2023-09-21 15:46:19.926196: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT\n" + ] + } + ], + "source": [ + "import os\n", + "import sys\n", + "sys.path.insert(0, '../../dlomix')\n", + "\n", + "from dlomix.pipelines import QRpipeline\n", + "from dlomix.pipelines import MCDPipeline" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "QRP = QRpipeline(frozen=1, alpha=0.1)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "187/187 [==============================] - 17s 81ms/step\n", + "Epoch 1/10\n", + "848/848 [==============================] - 85s 97ms/step - loss: 4.1107 - mean_absolute_error: 11.1265 - relative_central_distance: 6.6045 - interval_conformal_quantile: 11.5184 - abs_interval_size: 17.6196 - val_loss: 3.2001 - val_mean_absolute_error: 10.8963 - val_relative_central_distance: 0.6504 - val_interval_conformal_quantile: 8.9008 - val_abs_interval_size: 19.2391\n", + "Epoch 2/10\n", + "848/848 [==============================] - 87s 103ms/step - loss: 3.2098 - mean_absolute_error: 9.8797 - relative_central_distance: 0.8183 - interval_conformal_quantile: 9.7080 - abs_interval_size: 16.6748 - val_loss: 3.0620 - val_mean_absolute_error: 10.1642 - val_relative_central_distance: 0.6482 - val_interval_conformal_quantile: 8.8622 - val_abs_interval_size: 17.7547\n", + "Epoch 3/10\n", + "848/848 [==============================] - 90s 106ms/step - loss: 3.1757 - mean_absolute_error: 9.8762 - relative_central_distance: 0.7666 - interval_conformal_quantile: 9.6763 - abs_interval_size: 16.7516 - val_loss: 3.0207 - val_mean_absolute_error: 9.9785 - val_relative_central_distance: 0.6498 - val_interval_conformal_quantile: 8.7714 - val_abs_interval_size: 17.3937\n", + "Epoch 4/10\n", + "848/848 [==============================] - 78s 92ms/step - loss: 3.1472 - mean_absolute_error: 9.8832 - relative_central_distance: 0.7538 - interval_conformal_quantile: 9.3518 - abs_interval_size: 16.8408 - val_loss: 3.0158 - val_mean_absolute_error: 10.1754 - val_relative_central_distance: 0.6301 - val_interval_conformal_quantile: 8.4537 - val_abs_interval_size: 17.8983\n", + "Epoch 5/10\n", + "848/848 [==============================] - 82s 97ms/step - loss: 3.1347 - mean_absolute_error: 9.8094 - relative_central_distance: 0.7478 - interval_conformal_quantile: 9.4358 - abs_interval_size: 16.6877 - val_loss: 3.0500 - val_mean_absolute_error: 10.6772 - val_relative_central_distance: 0.5832 - val_interval_conformal_quantile: 7.9791 - val_abs_interval_size: 19.0669\n", + "Epoch 6/10\n", + "848/848 [==============================] - 78s 92ms/step - loss: 3.0917 - mean_absolute_error: 9.7033 - relative_central_distance: 0.7424 - interval_conformal_quantile: 9.2989 - abs_interval_size: 16.5298 - val_loss: 2.9642 - val_mean_absolute_error: 9.9040 - val_relative_central_distance: 0.6335 - val_interval_conformal_quantile: 8.4517 - val_abs_interval_size: 17.3490\n", + "Epoch 7/10\n", + "848/848 [==============================] - 78s 92ms/step - loss: 3.0974 - mean_absolute_error: 9.7681 - relative_central_distance: 0.7313 - interval_conformal_quantile: 9.1211 - abs_interval_size: 16.6768 - val_loss: 2.9377 - val_mean_absolute_error: 9.6069 - val_relative_central_distance: 0.6360 - val_interval_conformal_quantile: 8.7172 - val_abs_interval_size: 16.6724\n", + "Epoch 8/10\n", + "848/848 [==============================] - 82s 97ms/step - loss: 3.1089 - mean_absolute_error: 9.7790 - relative_central_distance: 0.7381 - interval_conformal_quantile: 9.4725 - abs_interval_size: 16.6761 - val_loss: 2.9556 - val_mean_absolute_error: 10.0426 - val_relative_central_distance: 0.6027 - val_interval_conformal_quantile: 8.2148 - val_abs_interval_size: 17.7168\n", + "Epoch 9/10\n", + "848/848 [==============================] - 84s 99ms/step - loss: 3.0770 - mean_absolute_error: 9.7329 - relative_central_distance: 0.7237 - interval_conformal_quantile: 9.1562 - abs_interval_size: 16.6410 - val_loss: 2.9336 - val_mean_absolute_error: 9.8513 - val_relative_central_distance: 0.6180 - val_interval_conformal_quantile: 8.3425 - val_abs_interval_size: 17.2935\n", + "Epoch 10/10\n", + "848/848 [==============================] - 89s 105ms/step - loss: 3.0455 - mean_absolute_error: 9.7010 - relative_central_distance: 0.7216 - interval_conformal_quantile: 9.0561 - abs_interval_size: 16.6396 - val_loss: 2.9649 - val_mean_absolute_error: 10.2913 - val_relative_central_distance: 0.5870 - val_interval_conformal_quantile: 7.9370 - val_abs_interval_size: 18.3152\n", + "187/187 [==============================] - 17s 86ms/step\n" + ] + } + ], + "source": [ + "QRP.run_full_pipeline()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#initialize Quantile Regression pipeline object\n", + "QRP = QRpipeline(frozen=1, alpha=0.1)\n", + "#call the pipeline function to run QR with CP\n", + "QRP.run_full_pipeline()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#initialize Monte Carlo Dropout pipeline object\n", + "MCDP = MCDPipeline(alpha=0.1)\n", + "#call the pipeline function to run MCD with CP\n", + "MCDP.run_full_pipeline()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "dlomix", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.17" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/QR_preds.pkl b/notebooks/QR_preds.pkl new file mode 100644 index 00000000..a780bddd Binary files /dev/null and b/notebooks/QR_preds.pkl differ diff --git a/notebooks/checkpoints/-1 layer(s)/best.hdf5 b/notebooks/checkpoints/-1 layer(s)/best.hdf5 new file mode 100644 index 00000000..8e24b2c3 Binary files /dev/null and b/notebooks/checkpoints/-1 layer(s)/best.hdf5 differ