diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index d88f9c7..973f71c 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -41,7 +41,7 @@ jobs:
     - name: Test with pytest
       run: |
         python3 -m pytest
-        python3 -m coverage run --source=./RadClass/ -m pytest
+        python3 -m coverage run --source=./RadClass/,./models/,./scripts/ -m pytest
         python3 -m coverage report
         python3 -m coverage html
         COVERALLS_REPO_TOKEN=${{ secrets.COVERALLS_REPO_TOKEN }} python3 -m coveralls --service=github
diff --git a/README.md b/README.md
index b08bd07..42245fa 100644
--- a/README.md
+++ b/README.md
@@ -25,7 +25,13 @@ Versions 3.6-3.9 are currently supported by tests. The following Python packages
 * h5py
 * numpy
 * progressbar2
+* matplotlib
+* seaborn
 * scipy
+* sklearn
+* hyperopt
+* torch
+* shadow-ssml
 
 Modules can be imported from the repository directory (e.g. `from RadClass.H0 import H0`) or `RadClass` can be installed using pip:
 
diff --git a/requirements.txt b/requirements.txt
index 06d1c3a..8b22315 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,3 +2,10 @@ numpy
 h5py
 progressbar2
 scipy>=1.7.0
+scikit-learn
+hyperopt
+matplotlib
+seaborn
+joblib
+torch
+shadow-ssml
diff --git a/scripts/__init__.py b/scripts/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/scripts/utils.py b/scripts/utils.py
new file mode 100644
index 0000000..d91c826
--- /dev/null
+++ b/scripts/utils.py
@@ -0,0 +1,328 @@
+import numpy as np
+import seaborn as sns
+import matplotlib.pyplot as plt
+# For hyperopt (parameter optimization)
+from hyperopt import Trials, tpe, fmin
+from functools import partial
+# diagnostics
+from sklearn.metrics import confusion_matrix
+# pca
+from sklearn.preprocessing import StandardScaler
+from sklearn.decomposition import PCA
+# Cross Validation
+from sklearn.model_selection import KFold, StratifiedKFold
+
+
+class EarlyStopper:
+    '''
+    Early stopping mechanism for neural networks.
+    Code adapted from user "isle_of_gods" from StackOverflow:
+    https://stackoverflow.com/questions/71998978/early-stopping-in-pytorch
+    Use this class to break a training loop if the validation loss is low.
+    Inputs:
+    patience: integer; forces stop if validation loss has not improved
+        for some time
+    min_delta: "fudge value" for how much loss to tolerate before stopping
+    '''
+
+    def __init__(self, patience=1, min_delta=0):
+        self.patience = patience
+        self.min_delta = min_delta
+        self.counter = 0
+        self.min_validation_loss = np.inf
+
+    def early_stop(self, validation_loss):
+        '''
+        Tests for the early stopping condition if the validation loss
+        has not improved for a certain period of time (patience).
+        Inputs:
+        validation_loss: typically a float value for the loss function of
+            a neural network training loop
+        '''
+
+        if validation_loss < self.min_validation_loss:
+            # keep track of the smallest validation loss
+            # if it has been beaten, restart patience
+            self.min_validation_loss = validation_loss
+            self.counter = 0
+        elif validation_loss > (self.min_validation_loss + self.min_delta):
+            # keep track of whether validation loss has been decreasing
+            # by a tolerable amount
+            self.counter += 1
+            if self.counter >= self.patience:
+                return True
+        return False
+
+
+def run_hyperopt(space, model, data_dict, max_evals=50, verbose=True):
+    '''
+    Runs hyperparameter optimization on a model given a parameter space.
+    Inputs:
+    space: dictionary with each hyperparameter as keys and values being the
+        range of parameter space (see hyperopt docs for defining a space)
+    mode: function that takes params dictionary, trains a specified ML model
+        and returns the optimization loss function, model, and other
+        attributes (e.g. accuracy on evaluation set)
+    max_eval: (int) run hyperparameter optimization for max_val iterations
+    verbose: report best and worse loss/accuracy
+
+    Returns:
+    best: dictionary with returns from model function, including best loss,
+        best trained model, best parameters, etc.
+    worst: dictionary with returns from model function, including worst loss,
+        worst trained model, worst parameters, etc.
+    '''
+
+    trials = Trials()
+
+    # wrap data into objective function
+    fmin_objective = partial(model, data_dict=data_dict)
+
+    # run hyperopt
+    fmin(fmin_objective,
+         space,
+         algo=tpe.suggest,
+         max_evals=max_evals,
+         trials=trials)
+
+    # of all trials, find best and worst loss/accuracy from optimization
+    best = trials.results[np.argmin([r['loss'] for r in trials.results])]
+    worst = trials.results[np.argmax([r['loss'] for r in trials.results])]
+
+    if verbose:
+        print('best accuracy:', 1-best['loss'])
+        print('best params:', best['params'])
+        print('worst accuracy:', 1-worst['loss'])
+        print('worst params:', worst['params'])
+
+    return best, worst
+
+
+def cross_validation(model, X, y, params, n_splits=3,
+                     stratified=False, random_state=None):
+    '''
+    Perform K-Fold cross validation using sklearn and a given model.
+    The model *must* have a fresh_start method (see models in RadClass/models).
+    fresh_start() is used instead of train() to be agnostic to the data needed
+        for training (fresh_start requires a data_dict whereas each model's
+        train could take different combinations of labeled & unlabeled data).
+        This also avoids the need to do hyperparameter optimization (and
+        therefore many training epochs) for every K-Fold.
+    NOTE: fresh_start returns the model and results in a dictionary but
+        does not overwrite/save the model to the respective class.
+        You can manually overwrite using model.model = return.model
+    Hyperparameter optimization (model.optimize) can be done before or after
+        cross validation to specify the (optimal) parameters used by the model
+        since they are required here.
+    NOTE: Fixed default to shuffle data during cross validation splits.
+        (See sklearn cross validation docs for more info.)
+    NOTE: Unlabeled data, if provided, will always be included in the training
+        dataset. This means that this cross validation implementation is
+        susceptible to bias in the unlabeled data distribution. To test for
+        this bias, a user can manually run cross validation as a parent to
+        calling this function, splitting the unlabeled data and adding
+        different folds into X.
+    Inputs:
+    model: ML model class object (e.g. RadClass/models).
+        Must have a fresh_start() method.
+        NOTE: If the model expects unlabeled data but unlabed data is not
+        provided in X/y, an error will likely be thrown when training the model
+        through fresh_start.
+    X: array of feature vectors (rows of individual instances, cols of vectors)
+        This should include all data for training and testing (since the
+        testing subset will be split by cross validation), including unlabeled
+        data if needed/used.
+    y: array/vector of labels for X. If including unlabeled data, use -1.
+        This should have the same order as X. That is, each row index in X
+        has an associated label with the same index in y.
+    params: dictionary of hyperparameters. Will depend on model used.
+        Alternatively, use model.params for models in RadClass/models
+    n_splits: int number of splits for K-Fold cross validation
+    stratified: bool; if True, balance the K-Folds to have roughly the same
+        proportion of samples from each class.
+    random_state: seed for reproducility.
+    '''
+
+    # return lists
+    accs = []
+    reports = []
+
+    if stratified:
+        cv = StratifiedKFold(n_splits=n_splits, random_state=random_state,
+                             shuffle=True)
+    else:
+        cv = KFold(n_splits=n_splits, random_state=random_state,
+                   shuffle=True)
+
+    # separate unlabeled data if included
+    Ux = None
+    Uy = None
+    if -1 in y:
+        U_idx = np.where(y == -1)[0]
+        L_idx = np.where(y != -1)[0]
+        Ux = X[U_idx]
+        Uy = y[U_idx]
+        Lx = X[L_idx]
+        Ly = y[L_idx]
+    else:
+        Lx = X
+        Ly = y
+    # conduct K-Fold cross validation
+    cv.get_n_splits(Lx, Ly)
+    for train_idx, test_idx in cv.split(Lx, Ly):
+        trainx, testx = Lx[train_idx], Lx[test_idx]
+        trainy, testy = Ly[train_idx], Ly[test_idx]
+
+        # construct data dictionary for training in fresh_start
+        data_dict = {'trainx': trainx, 'trainy': trainy,
+                     'testx': testx, 'testy': testy}
+        if Ux is not None:
+            data_dict['Ux'] = Ux
+            data_dict['Uy'] = Uy
+        results = model.fresh_start(params, data_dict)
+        accs = np.append(accs, results['accuracy'])
+        reports = np.append(reports, results)
+
+    # report cross validation results
+    print('Average accuracy:', np.mean(accs))
+    print('Max accuracy:', np.max(accs))
+    print('All accuracy:', accs)
+    # return the results of fresh_start for the max accuracy model
+    return reports[np.argmax(accs)]
+
+
+def pca(Lx, Ly, Ux, Uy, filename):
+    '''
+    A function for computing and plotting 2D PCA.
+    Inputs:
+    Lx: labeled feature data.
+    Ly: class labels for labeled data.
+    Ux: unlabeled feature data.
+    Uy: labels for unlabeled data (all labels should be -1).
+    filename: filename for saved plot.
+        The file must be saved with extension .joblib.
+        Added to filename if not included as input.
+    '''
+
+    plt.rcParams.update({'font.size': 20})
+    # only saving colors for binary classification with unlabeled instances
+    col_dict = {-1: 'tab:gray', 0: 'tab:orange', 1: 'tab:blue'}
+
+    pcadata = np.append(Lx, Ux, axis=0)
+    normalizer = StandardScaler()
+    x = normalizer.fit_transform(pcadata)
+    print(np.mean(pcadata), np.std(pcadata))
+    print(np.mean(x), np.std(x))
+
+    pca = PCA(n_components=2)
+    pca.fit_transform(x)
+    print(pca.explained_variance_ratio_)
+    print(pca.singular_values_)
+    print(pca.components_)
+
+    principalComponents = pca.fit_transform(x)
+
+    fig, ax = plt.subplots(figsize=(10, 8))
+    ax.set_xlabel('Principal Component 1', fontsize=15)
+    ax.set_ylabel('Principal Component 2', fontsize=15)
+    for idx, color in col_dict.items():
+        indices = np.where(np.append(Ly, Uy, axis=0) == idx)[0]
+        ax.scatter(principalComponents[indices, 0],
+                   principalComponents[indices, 1],
+                   c=color,
+                   label='class '+str(idx))
+    ax.grid()
+    ax.legend()
+
+    if filename[-4:] != '.png':
+        filename += '.png'
+    fig.tight_layout()
+    fig.savefig(filename)
+
+
+def multiD_pca(Lx, Ly, Ux, Uy, filename, n=2):
+    '''
+    A function for computing and plotting n-dimensional PCA.
+    Inputs:
+    Lx: labeled feature data.
+    Ly: class labels for labeled data.
+    Ux: unlabeled feature data.
+    Uy: labels for unlabeled data (all labels should be -1).
+    filename: filename for saved plot.
+        The file must be saved with extension .joblib.
+        Added to filename if not included as input.
+    n: number of singular values to include in PCA analysis.
+    '''
+
+    plt.rcParams.update({'font.size': 20})
+    # only saving colors for binary classification with unlabeled instances
+    col_dict = {-1: 'tab:gray', 0: 'tab:orange', 1: 'tab:blue'}
+
+    pcadata = np.append(Lx, Ux, axis=0)
+    normalizer = StandardScaler()
+    x = normalizer.fit_transform(pcadata)
+    print(np.mean(pcadata), np.std(pcadata))
+    print(np.mean(x), np.std(x))
+
+    n = 2
+    pca = PCA(n_components=n)
+    principalComponents = pca.fit_transform(x)
+    print(pca.explained_variance_ratio_)
+    print(pca.singular_values_)
+    print(pca.components_)
+
+    alph = ["A", "B", "C", "D", "E", "F", "G", "H",
+            "I", "J", "K", "L", "M", "N", "O", "P",
+            "Q", "R", "S", "T", "U", "V", "W", "X",
+            "Y", "Z"]
+    jobs = alph[:n]
+
+    fig, axes = plt.subplots(n, n, figsize=(15, 15))
+
+    for row in range(axes.shape[0]):
+        for col in range(axes.shape[1]):
+            ax = axes[row, col]
+            if row == col:
+                ax.tick_params(
+                    axis='both', which='both',
+                    bottom='off', top='off',
+                    labelbottom='off',
+                    left='off', right='off',
+                    labelleft='off'
+                )
+                ax.text(0.5, 0.5, jobs[row], horizontalalignment='center')
+            else:
+                for idx, color in col_dict.items():
+                    indices = np.where(np.append(Ly, Uy, axis=0) == idx)[0]
+                    ax.scatter(principalComponents[indices, row],
+                               principalComponents[indices, col],
+                               c=color,
+                               label='class '+str(idx))
+    fig.tight_layout()
+    if filename[-4:] != '.png':
+        filename += '.png'
+    fig.savefig(filename)
+
+
+def plot_cf(testy, predy, title, filename):
+    '''
+    Uses sklearn metric to compute a confusion matrix for visualization
+    Inputs:
+    testy: array/vector with ground-truth labels for test/evaluation set
+    predy: array/vector with predicted sample labels from trained model
+    title: string title for plot
+    filename: string with extension for confusion matrix file
+    '''
+
+    cf_matrix = confusion_matrix(testy, predy)
+    ax = sns.heatmap(cf_matrix, annot=True, cmap='Blues')
+
+    ax.set_title(title)
+    ax.set_xlabel('\nPredicted Values')
+    ax.set_ylabel('Actual Values ')
+
+    # Ticket labels - List must be in alphabetical order
+    ax.xaxis.set_ticklabels(['0(SNM)', '1(other)'])
+    ax.yaxis.set_ticklabels(['0(SNM)', '1(other)'])
+    # Save the visualization of the Confusion Matrix.
+    plt.savefig(filename)
diff --git a/tests/test_BackgroundEstimator.py b/tests/test_BackgroundEstimator.py
index 2d10c89..efc1299 100644
--- a/tests/test_BackgroundEstimator.py
+++ b/tests/test_BackgroundEstimator.py
@@ -77,7 +77,6 @@ def test_write():
     bckg.write(ofilename=ofilename)
 
     results = np.loadtxt(fname=ofilename+'.csv', delimiter=',')
-    print(results)
 
     # the resulting observation should be:
     #   counts * integration / live-time
diff --git a/tests/test_models.py b/tests/test_models.py
new file mode 100644
index 0000000..e3fb086
--- /dev/null
+++ b/tests/test_models.py
@@ -0,0 +1,491 @@
+# diagnostics
+import numpy as np
+from datetime import datetime, timedelta
+# testing models
+from sklearn.model_selection import train_test_split
+from sklearn.preprocessing import StandardScaler
+import tests.test_data as test_data
+# hyperopt
+from hyperopt.pyll.base import scope
+from hyperopt import hp
+# testing utils
+import scripts.utils as utils
+# models
+from models.LogReg import LogReg
+from models.SSML.CoTraining import CoTraining
+from models.SSML.LabelProp import LabelProp
+from models.SSML.ShadowNN import ShadowNN
+from models.SSML.ShadowCNN import ShadowCNN
+# testing write
+import joblib
+import os
+
+# initialize sample data
+start_date = datetime(2019, 2, 2)
+delta = timedelta(seconds=1)
+timestamps = np.arange(start_date,
+                       start_date + (test_data.timesteps * delta),
+                       delta).astype('datetime64[s]').astype('float64')
+
+live = np.full((len(timestamps),), test_data.livetime)
+sample_val = 1.0
+spectra = np.full((len(timestamps), test_data.energy_bins),
+                  np.full((1, test_data.energy_bins), sample_val))
+# setting up for rejected null hypothesis
+rejected_H0_time = np.random.choice(spectra.shape[0],
+                                    test_data.timesteps//2,
+                                    replace=False)
+spectra[rejected_H0_time] = 100.0
+
+labels = np.full((spectra.shape[0],), 0)
+labels[rejected_H0_time] = 1
+
+
+def test_utils():
+    X, Ux, y, Uy = train_test_split(spectra,
+                                    labels,
+                                    test_size=0.5,
+                                    random_state=0)
+    Uy = np.full_like(Uy, -1)
+
+    # test cross validation for supervised data using LogReg
+    params = {'max_iter': 2022, 'tol': 0.5, 'C': 5.0}
+    model = LogReg(params=params)
+    max_acc_model = utils.cross_validation(model=model,
+                                           X=X,
+                                           y=y,
+                                           params=params)
+    assert max_acc_model['accuracy'] >= 0.5
+
+    # test cross validation for supervised data and StratifiedKFold with LogReg
+    params = {'max_iter': 2022, 'tol': 0.5, 'C': 5.0}
+    model = LogReg(params=params)
+    max_acc_model = utils.cross_validation(model=model,
+                                           X=X,
+                                           y=y,
+                                           params=params,
+                                           stratified=True)
+    assert max_acc_model['accuracy'] >= 0.5
+
+    # test cross validation for SSML with LabelProp
+    params = {'gamma': 10, 'n_neighbors': 15, 'max_iter': 2022, 'tol': 0.5}
+    model = LabelProp(params=params)
+    max_acc_model = utils.cross_validation(model=model,
+                                           X=np.append(X, Ux, axis=0),
+                                           y=np.append(y, Uy, axis=0),
+                                           params=params,
+                                           stratified=True)
+    assert max_acc_model['accuracy'] >= 0.5
+
+    # data split for data visualization
+    X_train, X_test, y_train, y_test = train_test_split(X,
+                                                        y,
+                                                        test_size=0.2,
+                                                        random_state=0)
+
+    filename = 'test_pca'
+    utils.pca(X_train, y_train, Ux, np.full_like(Uy, -1), filename)
+    os.remove(filename+'.png')
+
+    filename = 'test_multiD_pca'
+    utils.multiD_pca(X_train, y_train, Ux, np.full_like(Uy, -1), filename, n=5)
+    os.remove(filename+'.png')
+
+    # normalization
+    normalizer = StandardScaler()
+    normalizer.fit(X_train)
+
+    X_train = normalizer.transform(X_train)
+    X_test = normalizer.transform(X_test)
+
+    # default behavior
+    model = LogReg(params=None, random_state=0)
+    model.train(X_train, y_train)
+
+    # testing train and predict methods
+    pred, acc = model.predict(X_test, y_test)
+
+    filename = 'test_cf'
+    utils.plot_cf(y_test, pred, title=filename, filename=filename)
+    os.remove(filename+'.png')
+
+
+def test_LogReg():
+    # test saving model input parameters
+    params = {'max_iter': 2022, 'tol': 0.5, 'C': 5.0}
+    model = LogReg(params=params)
+
+    assert model.model.max_iter == params['max_iter']
+    assert model.model.tol == params['tol']
+    assert model.model.C == params['C']
+
+    X_train, X_test, y_train, y_test = train_test_split(spectra,
+                                                        labels,
+                                                        test_size=0.2,
+                                                        random_state=0)
+
+    # normalization
+    normalizer = StandardScaler()
+    normalizer.fit(X_train)
+
+    X_train = normalizer.transform(X_train)
+    X_test = normalizer.transform(X_test)
+
+    # default behavior
+    model = LogReg(params=None, random_state=0)
+    model.train(X_train, y_train)
+
+    # testing train and predict methods
+    pred, acc = model.predict(X_test, y_test)
+
+    assert acc > 0.7
+    np.testing.assert_equal(pred, y_test)
+
+    # testing hyperopt optimize methods
+    space = {'max_iter': scope.int(hp.quniform('max_iter',
+                                               10,
+                                               10000,
+                                               10)),
+             'tol': hp.loguniform('tol', 1e-5, 1e-1),
+             'C': hp.uniform('C', 0.001, 1000.0)
+             }
+    data_dict = {'trainx': X_train,
+                 'testx': X_test,
+                 'trainy': y_train,
+                 'testy': y_test
+                 }
+    model.optimize(space, data_dict, max_evals=2, verbose=True)
+
+    assert model.best['accuracy'] >= model.worst['accuracy']
+    assert model.best['status'] == 'ok'
+
+    # testing model write to file method
+    filename = 'test_LogReg'
+    ext = '.joblib'
+    model.save(filename)
+    model_file = joblib.load(filename+ext)
+    assert model_file.best['params'] == model.best['params']
+
+    os.remove(filename+ext)
+
+
+def test_CoTraining():
+    # test saving model input parameters
+    params = {'max_iter': 2022, 'tol': 0.5, 'C': 5.0}
+    model = CoTraining(params=params)
+
+    assert model.model1.max_iter == params['max_iter']
+    assert model.model1.tol == params['tol']
+    assert model.model1.C == params['C']
+
+    assert model.model2.max_iter == params['max_iter']
+    assert model.model2.tol == params['tol']
+    assert model.model2.C == params['C']
+
+    X, Ux, y, Uy = train_test_split(spectra,
+                                    labels,
+                                    test_size=0.5,
+                                    random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X,
+                                                        y,
+                                                        test_size=0.2,
+                                                        random_state=0)
+
+    # normalization
+    normalizer = StandardScaler()
+    normalizer.fit(X_train)
+
+    X_train = normalizer.transform(X_train)
+    X_test = normalizer.transform(X_test)
+    Ux = normalizer.transform(Ux)
+
+    # default behavior
+    model = CoTraining(params=None, random_state=0)
+    model.train(X_train, y_train, Ux)
+
+    # testing train and predict methods
+    pred, acc, *_ = model.predict(X_test, y_test)
+
+    assert acc > 0.7
+    np.testing.assert_equal(pred, y_test)
+
+    # testing hyperopt optimize methods
+    space = {'max_iter': scope.int(hp.quniform('max_iter',
+                                               10,
+                                               10000,
+                                               10)),
+             'tol': hp.loguniform('tol', 1e-5, 1e-3),
+             'C': hp.uniform('C', 1.0, 1000.0),
+             'n_samples': scope.int(hp.quniform('n_samples',
+                                                1,
+                                                20,
+                                                1)),
+             'seed': 0
+             }
+    data_dict = {'trainx': X_train,
+                 'testx': X_test,
+                 'trainy': y_train,
+                 'testy': y_test,
+                 'Ux': Ux
+                 }
+    model.optimize(space, data_dict, max_evals=2, verbose=True)
+
+    assert model.best['accuracy'] >= model.worst['accuracy']
+    assert model.best['status'] == 'ok'
+
+    # testing model plotting method
+    filename = 'test_plot'
+    model.plot_cotraining(model1_accs=model.best['model1_acc_history'],
+                          model2_accs=model.best['model2_acc_history'],
+                          filename=filename)
+    os.remove(filename+'.png')
+
+    # testing model write to file method
+    filename = 'test_LogReg'
+    ext = '.joblib'
+    model.save(filename)
+    model_file = joblib.load(filename+ext)
+    assert model_file.best['params'] == model.best['params']
+
+    os.remove(filename+ext)
+
+
+def test_LabelProp():
+    # test saving model input parameters
+    params = {'gamma': 10, 'n_neighbors': 15, 'max_iter': 2022, 'tol': 0.5}
+    model = LabelProp(params=params)
+
+    assert model.model.gamma == params['gamma']
+    assert model.model.n_neighbors == params['n_neighbors']
+    assert model.model.max_iter == params['max_iter']
+    assert model.model.tol == params['tol']
+
+    # there should be no normalization on LabelProp data
+    # since it depends on the distances between samples
+    X, Ux, y, Uy = train_test_split(spectra,
+                                    labels,
+                                    test_size=0.5,
+                                    random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X,
+                                                        y,
+                                                        test_size=0.2,
+                                                        random_state=0)
+
+    # default behavior
+    model = LabelProp(params=None, random_state=0)
+    model.train(X_train, y_train, Ux)
+
+    # testing train and predict methods
+    pred, acc = model.predict(X_test, y_test)
+
+    # the default n_neighbors(=7) from sklearn is too large
+    # for the size of this dataset
+    # therefore the accuracy is expected to be poor
+    # a better value for this dataset would be n_neighbors=2
+    # (tested when specifying params in LabelProp.__init__)
+    assert acc >= 0.5
+    # uninteresting test if LabelProp predicts all one class
+    # TODO: make the default params test meaningful
+    assert np.count_nonzero(pred == y_test) > 0
+
+    # testing hyperopt optimize methods
+    space = {'max_iter': scope.int(hp.quniform('max_iter',
+                                               10,
+                                               10000,
+                                               10)),
+             'tol': hp.loguniform('tol', 1e-6, 1e-4),
+             'gamma': hp.uniform('gamma', 1, 50),
+             'n_neighbors': scope.int(hp.quniform('n_neighbors',
+                                                  1,
+                                                  X_train.shape[0],
+                                                  1))
+             }
+    data_dict = {'trainx': X_train,
+                 'testx': X_test,
+                 'trainy': y_train,
+                 'testy': y_test,
+                 'Ux': Ux
+                 }
+    model.optimize(space, data_dict, max_evals=2, verbose=True)
+
+    assert model.best['accuracy'] >= model.worst['accuracy']
+    assert model.best['status'] == 'ok'
+
+    # testing model write to file method
+    filename = 'test_LogReg'
+    ext = '.joblib'
+    model.save(filename)
+    model_file = joblib.load(filename+ext)
+    assert model_file.best['params'] == model.best['params']
+
+    os.remove(filename+ext)
+
+
+def test_ShadowNN():
+    # check default parameter settings
+    model = ShadowNN()
+    assert model.params == {'binning': 1}
+    assert model.eaat is not None
+    assert model.eaat_opt is not None
+    assert model.xEnt is not None
+    assert model.input_length == 1000
+
+    X, Ux, y, Uy = train_test_split(spectra,
+                                    labels,
+                                    test_size=0.5,
+                                    random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X,
+                                                        y,
+                                                        test_size=0.2,
+                                                        random_state=0)
+
+    # normalization
+    normalizer = StandardScaler()
+    normalizer.fit(X_train)
+
+    X_train = normalizer.transform(X_train)
+    X_test = normalizer.transform(X_test)
+    Ux = normalizer.transform(Ux)
+
+    params = {'hidden_layer': 10,
+              'alpha': 0.1,
+              'xi': 1e-3,
+              'eps': 1.0,
+              'lr': 0.1,
+              'momentum': 0.9,
+              'binning': 20}
+    # default behavior
+    model = ShadowNN(params=params, random_state=0)
+    acc_history = model.train(X_train, y_train, Ux, X_test, y_test)
+
+    # testing train and predict methods
+    pred, acc = model.predict(X_test, y_test)
+
+    # test for agreement between training and testing
+    # (since the same data is used for diagnostics in this test)
+    assert acc_history[-1] == acc
+
+    # Shadow/PyTorch reports accuracies as percentages
+    # rather than decimals
+    # uninteresting test if Shadow predicts all one class
+    # TODO: make the default params test meaningful
+    assert np.count_nonzero(pred == y_test) > 0
+
+    # testing hyperopt optimize methods
+    space = {'hidden_layer': 10,
+             'alpha': 0.1,
+             'xi': 1e-3,
+             'eps': 1.0,
+             'lr': 0.1,
+             'momentum': 0.9,
+             'binning': scope.int(hp.quniform('binning',
+                                              10,
+                                              20,
+                                              1))
+             }
+    data_dict = {'trainx': X_train,
+                 'testx': X_test,
+                 'trainy': y_train,
+                 'testy': y_test,
+                 'Ux': Ux
+                 }
+    model.optimize(space, data_dict, max_evals=2, verbose=True)
+
+    assert model.best['accuracy'] >= model.worst['accuracy']
+    assert model.best['status'] == 'ok'
+
+    # testing model write to file method
+    filename = 'test_LogReg'
+    ext = '.joblib'
+    model.save(filename)
+    model_file = joblib.load(filename+ext)
+    assert model_file.best['params'] == model.best['params']
+
+    os.remove(filename+ext)
+
+
+def test_ShadowCNN():
+    # check default parameter settings
+    model = ShadowCNN()
+    assert model.params == {'binning': 1, 'batch_size': 1}
+    assert model.model is not None
+    assert model.eaat is not None
+    assert model.optimizer is not None
+
+    X, Ux, y, Uy = train_test_split(spectra,
+                                    labels,
+                                    test_size=0.5,
+                                    random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X,
+                                                        y,
+                                                        test_size=0.2,
+                                                        random_state=0)
+
+    # normalization
+    normalizer = StandardScaler()
+    normalizer.fit(X_train)
+
+    X_train = normalizer.transform(X_train)
+    X_test = normalizer.transform(X_test)
+    Ux = normalizer.transform(Ux)
+
+    params = {'layer1': 2,
+              'kernel': 3,
+              'alpha': 0.1,
+              'xi': 1e-3,
+              'eps': 1.0,
+              'lr': 0.1,
+              'momentum': 0.9,
+              'binning': 20,
+              'batch_size': 4,
+              'drop_rate': 0.1}
+
+    # default behavior
+    model = ShadowCNN(params=params, random_state=0)
+    losscurve, evalcurve = model.train(X_train, y_train, Ux, X_test, y_test)
+
+    # testing train and predict methods
+    pred, acc = model.predict(X_test, y_test)
+
+    # test for agreement between training and testing
+    # (since the same data is used for diagnostics in this test)
+    assert evalcurve[-1] == acc
+
+    # Shadow/PyTorch reports accuracies as percentages
+    # rather than decimals
+    # uninteresting test if Shadow predicts all one class
+    # TODO: make the default params test meaningful
+    assert np.count_nonzero(pred == y_test) > 0
+
+    # testing hyperopt optimize methods
+    space = params
+    space['binning'] = scope.int(hp.quniform('binning',
+                                 10,
+                                 20,
+                                 1))
+    data_dict = {'trainx': X_train,
+                 'testx': X_test,
+                 'trainy': y_train,
+                 'testy': y_test,
+                 'Ux': Ux
+                 }
+    model.optimize(space, data_dict, max_evals=2, verbose=True)
+
+    assert model.best['accuracy'] >= model.worst['accuracy']
+    assert model.best['status'] == 'ok'
+
+    # testing model plotting method
+    filename = 'test_plot'
+    model.plot_training(losscurve=model.best['losscurve'],
+                        evalcurve=model.best['evalcurve'],
+                        filename=filename)
+    os.remove(filename+'.png')
+
+    # testing model write to file method
+    filename = 'test_LogReg'
+    ext = '.joblib'
+    model.save(filename)
+    model_file = joblib.load(filename+ext)
+    assert model_file.best['params'] == model.best['params']
+
+    os.remove(filename+ext)