vespuccidataset.cpp

Go to the documentation of this file.

 /*******************************************************************************
    Copyright (C) 2014-2016 Wright State University - All Rights Reserved
    Daniel P. Foose - Maintainer/Lead Developer

    This file is part of Vespucci.

    Vespucci is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Vespucci is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Vespucci.  If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "Data/Analysis/abstractdataanalyzer.h"
#include "Data/Dataset/vespuccidataset.h"
#include "Data/Import/textimport.h"
#include "Math/Clustering/agglomerativeclustering.h"
#include <H5Cpp.h>
#include "Math/Baseline/rollingball.h"
#include <Math/Smoothing/denoise.h>
#include <Math/Clustering/kmeanswrapper.h>
using namespace arma;
using namespace std;


VespucciDataset::~VespucciDataset()
{
    for (auto key : MapKeys())
        workspace_->GetMap(name_, key)->HideMapWindow();

}

bool VespucciDataset::Save(QString filename)
{
    using namespace H5;
    QStringList core_names = CoreMatrixKeys();
    QStringList results_names = AnalysisResultsKeys();
    try{
        H5File file(filename.toStdString(), H5F_ACC_TRUNC);
        DataSpace str_dataspace(H5S_SCALAR);
        StrType str_type(PredType::C_S1, H5T_VARIABLE);

        Attribute str_attr;
        if (file.attrExists("Name")) str_attr = file.openAttribute("Name");
        else str_attr = file.createAttribute("Name", str_type, str_dataspace);
        str_attr.write(str_type, name_.toStdString());

        if (file.attrExists("x-Axis Description")) str_attr = file.openAttribute("x-Axis Description");
        else str_attr = file.createAttribute("x-Axis Description", str_type, str_dataspace);
        str_attr.write(str_type, x_axis_description_.toStdString());

        if (file.attrExists("y-Axis Description")) str_attr = file.openAttribute("y-Axis Description");
        else str_attr = file.createAttribute("y-Axis Description", str_type, str_dataspace);
        str_attr.write(str_type, y_axis_description_.toStdString());

        for (auto core_name: core_names){
            hsize_t dims[2];
            dims[0] = GetCoreMatrix(core_name).n_rows;
            dims[1] = GetCoreMatrix(core_name).n_cols;
            DataSpace dataspace(2, dims);
            DataSet ds(file.createDataSet(core_name.toStdString(),
                                          PredType::NATIVE_DOUBLE, dataspace));
            ds.write(GetCoreMatrix(core_name).memptr(),
                     PredType::NATIVE_DOUBLE, dataspace);
            ds.close();
        }

        for (auto results_name: results_names){
            QSharedPointer<AnalysisResults> results = GetAnalysisResult(results_name);
            QStringList matrix_names = results->KeyList();
            Group results_group(file.createGroup(results_name.toStdString()));
            if (results_group.attrExists("type")) str_attr = results_group.openAttribute("type");
            else str_attr = results_group.createAttribute("type", str_type, str_dataspace);
            std::string results_type = results->type().toStdString();
            str_attr.write(str_type, results_type);
            for (auto matrix_name: matrix_names){
                hsize_t dims[2];
                dims[0] = results->GetMatrix(matrix_name).n_rows;
                dims[1] = results->GetMatrix(matrix_name).n_cols;
                DataSpace dataspace(2, dims);
                DataSet ds(results_group.createDataSet(matrix_name.toStdString(),
                                                       PredType::NATIVE_DOUBLE,
                                                       dataspace));
                ds.write(results->GetMatrix(matrix_name).memptr(),
                         PredType::NATIVE_DOUBLE);
                ds.close();
            }
            results_group.close();
        }
        file.close();
    }
    catch(FileIException error){
        QString msg = QString::fromStdString(error.getDetailMsg());
        main_window_->DisplayWarning("HDF5 FileIException", msg);
        return false;
    }
    catch(DataSetIException error){
        QString msg = QString::fromStdString(error.getDetailMsg());
        main_window_->DisplayWarning("HDF5 DataSetIException", msg);
        return false;
    }
    catch(AttributeIException error){
        QString msg = QString::fromStdString(error.getDetailMsg());
        main_window_->DisplayWarning("HDF5 AttributeIException", msg);
        return false;
    }

    last_save_filename_ = filename;
    saved_ = true;
    state_changed_ = false;
    return true;
}

bool VespucciDataset::SaveSpectrum(QString filename, uword column, file_type type)
{
    std::string filename_stdstring = filename.toStdString();
    vec spectrum;
    bool success;
    try{
        spectrum = spectra_.col(column);
        success = spectrum.save(filename_stdstring, type);
    }
    catch(exception e){
        main_window_->DisplayExceptionWarning("VespucciDataset::SaveSpectrum", e);
        success = false;
    }
    return success;

}

void VespucciDataset::SetOldCopies()
{
    spectra_old_ = spectra_;
    x_old_ = x_;
    y_old_ = y_;
    abscissa_old_ = abscissa_;
}

bool VespucciDataset::Contains(const QString &key) const
{
    return CoreMatrixKeys().contains(key) || AnalysisResultsKeys().contains(key)
            || AuxiliaryMatrixKeys().contains(key) || key == "Auxiliary Matrices";
}

bool VespucciDataset::IsValid() const
{
    bool has_spatial = x_.n_rows;
    //x and y have same size
    bool valid_spatial = has_spatial && (x_.n_rows == y_.n_rows);
    //x and y have same number of rows as spectra has columns
    bool spatial_match = valid_spatial && (x_.n_rows == spectra_.n_cols);
    //spectra and abscissa have same number of rows
    bool spectral_match = spectra_.n_elem && abscissa_.n_rows == spectra_.n_rows;
    return spectral_match && (spatial_match || !has_spatial);
}

bool VespucciDataset::state_changed() const
{
    return state_changed_;
}

QString VespucciDataset::filename() const
{
    return filename_;
}

QString VespucciDataset::last_save_filename() const
{
    return last_save_filename_;
}

bool VespucciDataset::saved() const
{
    return saved_;
}


VespucciDataset::VespucciDataset(const QString &h5_filename,
                                 MainWindow *main_window,
                                 QSharedPointer<VespucciWorkspace> ws)
    : auxiliary_matrices_(new AnalysisResults("Auxiliary Matrices", "Auxiliary Matrices"))
{
    saved_ = true;
    filename_ = h5_filename;
    last_save_filename_ = h5_filename;
    workspace_ = ws;
    QString fallback_name = "Dataset";
    directory_ = workspace_->directory_ptr();
    int counter = 1;
    while (workspace_->dataset_names().contains(fallback_name))
        fallback_name = "Dataset " + QString::number(counter++);
    QString abs_label = workspace_->settings()->value("absLabel").toString();
    QString abs_units = workspace_->settings()->value("absUnits").toString();
    QString ord_label = workspace_->settings()->value("ordLabel").toString();
    QString ord_units = workspace_->settings()->value("ordUnits").toString();
    QString fallback_x_description = abs_label + "(" + abs_units + ")";
    QString fallback_y_description = ord_label + "(" + ord_units + ")";
    main_window_ = main_window;
    if (!QFile::exists(h5_filename)){
        constructor_canceled_ = true;
        return;
    }
    using namespace H5;

    try{
        H5File file(h5_filename.toStdString(), H5F_ACC_RDONLY);

        if (file.attrExists("Name")){
            Attribute attr = file.openAttribute("Name");
            H5std_string buf("");
            DataType type = attr.getDataType();
            attr.read(type, buf);
            name_ = QString::fromStdString(buf);
        }
        else{
            name_ = fallback_name;
        }
        if (file.attrExists("x-Axis Description")){
            Attribute attr = file.openAttribute("x-Axis Description");
            H5std_string buf("");
            DataType type = attr.getDataType();
            attr.read(type, buf);
            x_axis_description_ = QString::fromStdString(buf);
        }
        else{
            x_axis_description_ = fallback_x_description;
        }
        if (file.attrExists("y-Axis Description")){
            Attribute attr = file.openAttribute("y-Axis Description");
            H5std_string buf("");
            DataType type = attr.getDataType();
            attr.read(type, buf);
            y_axis_description_ = QString::fromStdString(buf);
        }
        else{
            y_axis_description_ = fallback_y_description;
        }

        for (hsize_t i = 0; i < file.getNumObjs(); ++i){
            H5G_obj_t obj_type = file.getObjTypeByIdx(i);
            string obj_name = file.getObjnameByIdx(i);

            if (obj_type == H5G_DATASET){
                DataSet ds = file.openDataSet(obj_name);
                DataSpace dspace = ds.getSpace();
                hsize_t dims[2];
                dspace.getSimpleExtentDims(dims);
                mat matrix(dims[0], dims[1]);
                ds.read(matrix.memptr(), PredType::NATIVE_DOUBLE);
                AddMatrix(QString::fromStdString(obj_name), matrix);
                ds.close();
            }
            if (obj_type == H5G_GROUP){
                Group group = file.openGroup(obj_name);
                H5std_string result_type("");
                if (group.attrExists("type")){
                    Attribute attr = group.openAttribute("type");
                    DataType type = attr.getDataType();
                    attr.read(type, result_type);
                }
                QSharedPointer<AnalysisResults> result(new AnalysisResults(QString::fromStdString(obj_name), QString::fromStdString(result_type)));
                for (hsize_t j = 0; j < group.getNumObjs(); ++ j){
                    H5G_obj_t subobj_type = group.getObjTypeByIdx(j);
                    string subobj_name = group.getObjnameByIdx(j);
                    if (subobj_type == H5G_DATASET){
                        DataSet ds = group.openDataSet(subobj_name);
                        DataSpace dspace = ds.getSpace();
                        hsize_t dims[2];
                        dspace.getSimpleExtentDims(dims);
                        mat matrix(dims[0], dims[1]);
                        ds.read(matrix.memptr(), PredType::NATIVE_DOUBLE);
                        result->AddMatrix(QString::fromStdString(subobj_name),
                                          matrix);
                        ds.close();
                    }//if (subobj_type == H5G_DATASET
                }//for objects in group
                AddAnalysisResult(result);
                group.close();
            }//if (obj_type == H5G_GROUP)
        }//for objects in file

        file.close();
        if (AnalysisResultsKeys().contains("Auxiliary Matrices")){
            auxiliary_matrices_ = GetAnalysisResult("Auxiliary Matrices");
        }
        else{
            auxiliary_matrices_ = QSharedPointer<AnalysisResults>(new AnalysisResults("Auxiliary Matrices", "Auxiliary Matrices"));
            AddAnalysisResult(auxiliary_matrices_);
        }
        constructor_canceled_ = false;
    }//try
    catch(H5::Exception error){
        cerr << error.getDetailMsg() << "\n";
        constructor_canceled_ = IsValid();
    }
    catch(std::exception e){
        cerr << e.what() << "\n";
        constructor_canceled_ = IsValid();
    }

    if (!constructor_canceled_){
        if (AnalysisResultsKeys().contains("Auxiliary Matrices")){
            auxiliary_matrices_ = GetAnalysisResult("Auxiliary Matrices");
        }
        else{
            AddAnalysisResult(auxiliary_matrices_);
        }
    }

}

VespucciDataset::VespucciDataset(QString text_filename,
                                 MainWindow *main_window,
                                 QString *directory,
                                 QString name,
                                 QString x_axis_description,
                                 QString y_axis_description,
                                 bool swap_spatial,
                                 std::string format)
    : auxiliary_matrices_(new AnalysisResults("Auxiliary Matrices", "Auxiliary Matrices"))
{
    saved_ = false;
    filename_ = text_filename;
    workspace_ = main_window->workspace_ptr();
    QDateTime datetime = QDateTime::currentDateTimeUtc();
    state_changed_ = true;
    AddAnalysisResult(auxiliary_matrices_);

    non_spatial_ = false;
    meta_ = false;
    //Set up variables unrelated to hyperspectral data:
    //map_list_model_ = new MapListModel(main_window, this);
    map_loading_count_ = 0;
    z_scores_calculated_ = false;
    directory_ = directory;
    flipped_ = swap_spatial;

    QProgressDialog progress("Loading Dataset...", "Cancel", 0, 100, NULL);
    vec indices_temp;


    if (format == "WideTabDel"){
        try{
            constructor_canceled_ = TextImport::ImportWideText(text_filename.toStdString(),
                                       spectra_,
                                       abscissa_,
                                       x_, y_, swap_spatial);
            indices_ = linspace<vec>(0, x_.n_elem - 1, x_.n_elem);
        }
        catch(exception e){
            string str = "Text import constructor:" +  string(e.what());
            throw std::runtime_error(str);
        }
    }
    else if (format == "WideCSV"){
        try{
            constructor_canceled_ = TextImport::ImportWideText(text_filename.toStdString(),
                                       spectra_,
                                       abscissa_,
                                       x_, y_, swap_spatial);
            indices_ = linspace<vec>(0, x_.n_elem - 1, x_.n_elem);
        }
        catch(exception e){
            string str = "Text import constructor: " + string(e.what());
            throw std::runtime_error(str);
        }
    }
    else if (format == "LongTabDel" || format == "LongCSV"){
        try{

            constructor_canceled_ = TextImport::ImportLongText(text_filename.toStdString(),
                                                               spectra_,
                                                               abscissa_,
                                                               x_, y_,
                                                               swap_spatial);
            indices_.set_size(x_.n_elem);
            for (uword i = 0; i < indices_.n_elem; ++i)
                indices_(i) = i;
        }
        catch(exception e){
            string str = "Text import constructor: " + string(e.what());
            throw std::runtime_error(str);
        }
    }
    else{
        throw std::runtime_error("Invalid file format for text import constructor");
    }

    constructor_canceled_ = false;
    name_ = name;
    x_axis_description_ = x_axis_description;
    y_axis_description_ = y_axis_description;
    main_window_ = main_window;
}

VespucciDataset::VespucciDataset(map<pair<int, int>, string> text_filenames,
                                 MainWindow *main_window,
                                 QString *directory,
                                 QString name,
                                 QString x_axis_description,
                                 QString y_axis_description,
                                 int rows, int cols)
     : auxiliary_matrices_(new AnalysisResults("Auxiliary Matrices", "Auxiliary Matrices"))
{
    saved_ = false;
    text_filenames_ = text_filenames;
    workspace_ = main_window->workspace_ptr();
    AddAnalysisResult(auxiliary_matrices_);
    QDateTime datetime = QDateTime::currentDateTimeUtc();
    non_spatial_ = false;
    meta_ = false;
    //Set up variables unrelated to hyperspectral data:
    //map_list_model_ = new MapListModel(main_window, this);
    map_loading_count_ = 0;
    z_scores_calculated_ = false;
    directory_ = directory;
    flipped_ = false;
    try{
        constructor_canceled_ = TextImport::ImportMultiplePoints(text_filenames,
                                                                 rows, cols,
                                                                 spectra_,
                                                                 abscissa_,
                                                                 x_,
                                                                 y_);

    }
    catch(exception e){
        string reason = "Multi import constructor: " + string(e.what());
        throw std::runtime_error(reason.c_str());
    }
    name_ = name;
    x_axis_description_ = x_axis_description;
    y_axis_description_ = y_axis_description;
    main_window_ = main_window;
}

VespucciDataset::VespucciDataset(QString name,
                                 MainWindow *main_window,
                                 QString *directory,
                                 QSharedPointer<VespucciDataset> original,
                                 uvec indices)
    : auxiliary_matrices_(new AnalysisResults("Auxiliary Matrices", "Auxiliary Matrices"))

{
    saved_ = false;
    workspace_ = main_window->workspace_ptr();
    AddAnalysisResult(auxiliary_matrices_);
    state_changed_ = true;
    non_spatial_ = true;
    meta_ = original->meta();
    map_loading_count_ = 0;
    //principal_components_calculated_ = false;
    //mlpack_pca_calculated_ = false;
    //partial_least_squares_calculated_ = false;
    //vertex_components_calculated_ = false;
    z_scores_calculated_ = false;
    directory_ = directory;
    vec parent_indices;

    try{
        spectra_ = original->spectra(indices);
        abscissa_ = original->wavelength();
        x_ = original->x(indices);
        y_ = original->y(indices);
        parent_indices = original->indices();
        indices_ = parent_indices(indices);
    }
    catch(exception e){
        string str = "Extraction constructor: " + string(e.what());
        throw std::runtime_error(str);
    }

    name_ = name;
    main_window_ = main_window;
    directory_ = directory;
}

VespucciDataset::VespucciDataset(QString name,
                                 MainWindow *main_window,
                                 QString *directory)
    : auxiliary_matrices_(new AnalysisResults("Auxiliary Matrices", "Auxiliary Matrices"))
{
    workspace_ = main_window->workspace_ptr();
    AddAnalysisResult(auxiliary_matrices_);
    state_changed_ = true;
    non_spatial_ = true;
    meta_ = true;
    map_loading_count_ = 0;
    z_scores_calculated_ = false;
    directory_ = directory;
    name_ = name;
    main_window_ = main_window;
    directory_ = directory;
    saved_ = false;
}

// PRE-PROCESSING FUNCTIONS //
void VespucciDataset::Undo()
{

    try{
        mat spectra_buffer = spectra_;
        vec x_buffer = x_;
        vec y_buffer = y_;
        vec abscissa_buffer = abscissa_;

        spectra_.swap(spectra_old_);
        x_.swap(x_old_);
        y_.swap(y_old_);
        abscissa_.swap(abscissa_old_);

        //spectra_ = spectra_old_;
        //x_ = x_old_;
        //y_ = y_old_;
        //abscissa_ = abscissa_old_;

        //spectra_old_ = spectra_buffer;
        //x_old_ = x_buffer;
        //y_old_ = y_buffer;
        //abscissa_old_ = abscissa_buffer;


    }
    catch(exception e){
        string str = "Undo: " + string(e.what());
        throw std::runtime_error(str);
    }

    state_changed_ = true;
    last_operation_ = "Undo";
}


void VespucciDataset::CropSpectra(double x_min, double x_max,
                                  double y_min, double y_max,
                                  double wl_min, double wl_max)
{

    SetOldCopies();
    if (!std::isnan(x_min) && !std::isnan(x_max) && !std::isnan(y_min) && !std::isnan(y_max)){
        uvec valid_indices = find ((x_ >= x_min) && (x_ <= x_max));
        try{
            spectra_ = spectra_.cols(valid_indices);
            y_ = y_.rows(valid_indices);
            x_ = x_.rows(valid_indices);
            valid_indices = find((y_ >= y_min) && (y_ <= y_max));
            y_ = y_.rows(valid_indices);
            x_ = x_.rows(valid_indices);
            spectra_ = spectra_.cols(valid_indices);
        }catch(exception e){
            string str = "CropSpectra: " + string(e.what());
            throw std::runtime_error(str);
        }
    }
    if (!std::isnan(wl_min) && !std::isnan(wl_max)){
        uvec valid_indices = find(abscissa_ >= wl_min && abscissa_ <= wl_max);
        try{
            spectra_ = spectra_.rows(valid_indices);
            abscissa_ = abscissa_.rows(valid_indices);

        }catch(exception e){
            string str = "CropSpectra: " + string(e.what());
            throw std::runtime_error(str);
        }
    }

    last_operation_ = "crop";
    operations_ << "CropSpectra("
                   + QString::number(x_min) + ", "
                   + QString::number(x_max) + ", "
                   + QString::number(y_min) + ", "
                   + QString::number(y_max) + ", "
                   + QString::number(wl_min) + ", "
                   + QString::number(wl_max) + ")";
    state_changed_ = true;
    workspace_->UpdateModel();
}


void VespucciDataset::MinMaxNormalize()
{
    try{
        SetOldCopies();
        vec spectrum_buffer;
        vec extrema_buffer;
        for (uword i = 0; i < spectra_.n_cols; ++i){
            spectrum_buffer = spectra_.col(i);
            extrema_buffer = spectrum_buffer.min()*ones<vec>(spectra_.n_rows);
            spectrum_buffer -= extrema_buffer;
            extrema_buffer = spectrum_buffer.max()*ones<vec>(spectra_.n_rows);
            spectrum_buffer /= extrema_buffer;
            spectra_.col(i) = spectrum_buffer;
        }
    }
    catch(exception e){
        string str = "MinMaxNormalize: " + string(e.what());
        throw std::runtime_error(str);
    }
    last_operation_ = "min/max normalize";
    operations_ << "MinMaxNormalize()";
    state_changed_ = true;
}

void VespucciDataset::VectorNormalize(uword norm)
{
    if (!norm || norm > 2) return;
    state_changed_ = true;

    SetOldCopies();
    try{
        spectra_ = normalise(spectra_, norm);
    }
    catch(exception e){
        string str = "VectorNormalize()" + string(e.what());
        throw std::runtime_error(str);
    }
    operations_ << "VectorNormalize(" + QString::number(norm) + ")";
    last_operation_ = "vector normalize";
}

void VespucciDataset::MeanCenter()
{
    state_changed_ = true;
    SetOldCopies();
    try{
        vec mean_intensity = mean(spectra_, 1);
        spectra_.each_col() -= mean_intensity;
    }
    catch(exception e){
        throw std::runtime_error("MeanCenter");
    }
    last_operation_ = "mean centering";
    operations_ << "MeanCenter()";
}

void VespucciDataset::PeakIntensityNormalize(double left_bound, double right_bound)
{
    SetOldCopies();
    vec positions;
    vec peak_maxes = Vespucci::Math::Quantification::FindPeakMaxMat(spectra_, abscissa_, left_bound, right_bound, positions);
    for (uword j = 0; j < spectra_.n_cols; ++j){
        spectra_.col(j) /= peak_maxes(j);
    }
    last_operation_ = "Peak intensity normalize";
    operations_ << "PeakIntensityNormalize("
                   + QString::number(left_bound) + ", "
                   + QString::number(right_bound) + ")";
}

void VespucciDataset::Booleanize(double min, double max, bool keep_inside, bool oneify)
{
    SetOldCopies();
    last_operation_ = "Booleanize";

    //set all values on the outside (or inside) of the range to zero
    try{
        if (keep_inside)
            spectra_.elem(find(spectra_ <= min && spectra_ >= max) ).zeros();
        else
            spectra_.elem(find(spectra_ >= min && spectra_ <= max) ).zeros();


        //set all non-zero values equal to one
        if (oneify)
            spectra_.elem( find(spectra_) ).ones();

    }
    catch (exception e){
        cerr << e.what();
        runtime_error f("Booleanize");
        main_window_->DisplayExceptionWarning(f);
    }
    operations_ << "Booleanize(" + QString::number(min) + ", "
                   + QString::number(max) + ", "
                   + (keep_inside ? "true" : "false") + ", "
                   + (oneify ? "true" : "false") + ")";
}

void VespucciDataset::Clamp(double min, double max)
{
    SetOldCopies();
    last_operation_ = "Clamp";
    state_changed_ = true;


    try{
        spectra_ = clamp(spectra_, min, max);
    }
    catch(exception e){
        cerr << e.what();
        runtime_error f("Clamp");
        main_window_->DisplayExceptionWarning(f);
    }
    operations_ << "Clamp("
                   + QString::number(min) + ", "
                   + QString::number(max) + ")";
}

void VespucciDataset::ShedZeroSpectra()
{
    SetOldCopies();
    vec current_col;
    uvec indices;
    uvec indices_buffer;
    for (uword i = 0; i < spectra_.n_cols; ++i){
        current_col = spectra_.col(i);
        indices_buffer = find(current_col, 1);
        if (indices_buffer.n_elem > 0)
            indices << i;
    }
    try{
        spectra_ = spectra_.cols(indices);
        x_ = x_.rows(indices);
        y_ = y_.rows(indices);
    }
    catch(exception e){
        main_window_->DisplayExceptionWarning("VespucciDataset::ShedZeroSpectra", e);
    }

    workspace_->UpdateModel();
    operations_ << "ShedZeroSpectra()";
}

void VespucciDataset::ShedZeroWavelengths()
{
    state_changed_ = true;

    SetOldCopies();
    vec current_row;
    uvec indices, indices_buffer;
    for (uword i = 0; i < spectra_.n_rows; ++i){
        current_row = spectra_.row(i);
        indices_buffer = find(current_row, 1);
        if (indices_buffer.n_elem > 0)
            indices << i;
    }
    try{
        spectra_ = spectra_.rows(indices);
        abscissa_ = abscissa_.rows(indices);
    }
    catch(exception e){
        main_window_->DisplayExceptionWarning("VespucciDataset::ShedZeroWavelengths", e);
    }
    workspace_->UpdateModel();
    operations_ << "ShedZeroWavelengths()";
}

mat VespucciDataset::ZScoreNormCopy()
{

    mat normalized;
    //mat normalized_copy(num_rows, num_cols);
    try{
        normalized = Vespucci::Math::Normalization::StandardScore(spectra_);
    }
    catch(exception e){
        string str = "ZScoreNormCopy: " + string(e.what());
        throw std::runtime_error(str);
    }
    return trans(normalized);

}

void VespucciDataset::ZScoreNormalize()
{
    state_changed_ = true;

    SetOldCopies();

    try{
        spectra_ = Vespucci::Math::Normalization::StandardScoreMat(spectra_);
    }
    catch(exception e){
        string str = "ZScoreNormalize: " + string(e.what());
        throw std::runtime_error(str);
    }

    z_scores_calculated_ = true;
    last_operation_ = "Z-score normalize";
    operations_ << "ZScoreNormalize()";

}


void VespucciDataset::SNVNormalize(double offset, bool center)
{
    state_changed_ = true;

    SetOldCopies();
    try{
        spectra_ = Vespucci::Math::Normalization::SNVNorm(spectra_, offset, center);
    }
    catch(exception e){
        string str = "SNVNormalize: " + string(e.what());
        throw std::runtime_error(str);
    }

    last_operation_ = "SNVNormalize";
    operations_ << "SNVNormalize()";
}

void VespucciDataset::AbsoluteValue()
{
    state_changed_ = true;

    SetOldCopies();

    try{
        spectra_ = arma::abs(spectra_);
    }
    catch(exception e){
        string str =
        "AbsoluteValue: " + string(e.what());
        throw std::runtime_error(str);
    }

    last_operation_ = "Absolute value";
    operations_ << "AbsoluteValue()";
}

void VespucciDataset::SubtractBackground(const QStringList &data_keys)
{
    mat background = workspace_->GetMatrix(data_keys);
    state_changed_ = true;
    SetOldCopies();
    if (background.n_rows != spectra_.n_rows){
        QMessageBox::warning(0,
                             "Improper Dimensions!",
                             "The background spectrum has a different number of"
                             " points than the map data."
                             " No subtraction can be performed");
        return;
    }
    else{
        try{
            spectra_.each_col() -= background.col(0);
        }
        catch(exception e){
            string str = "SubtractBackground: " + string(e.what());
            throw std::runtime_error(str);
        }
    }
    last_operation_ = "background correction";
}

void VespucciDataset::MFBaseline(int window_size, int iterations)
{
    state_changed_ = true;
    SetOldCopies();
    try{
        baselines_ = spectra_;
        for (uword i = 0; i < uword(iterations); ++i){
            baselines_ = Vespucci::Math::Smoothing::MedianFilterMat(baselines_, window_size);
        }
        spectra_ -= baselines_;
    }
    catch(exception e){
        string str = "MFBaseline: " + string(e.what());
        throw std::runtime_error(str);
    }

    last_operation_ = "baseline correction (median filter)";
    operations_ << "MFBaseline("
                   + QString::number(window_size) + ", "
                   + QString::number(iterations) + ")";
}

void VespucciDataset::RollingBallBaseline(size_t wm, size_t ws)
{
    state_changed_ = true;
    SetOldCopies();
    try{
        mat baselines;
        spectra_ = Vespucci::Math::Baseline::RollingBallBaselineMat(spectra_,
                                                                    baselines,
                                                                    wm, ws);
        AddAuxiliaryMatrix("Rolling Ball Baselines", baselines);
    }catch(exception e){
        string str = "RollingBallBaseline " + string(e.what());
        throw std::runtime_error(str);
    }
    last_operation_ = "baseline correction (rolling ball)";
    operations_ << "RollingBallBaseline("
                   + QString::number(wm) + ", "
                   + QString::number(ws) + ")";
}

void VespucciDataset::CWTBaseline(int lambda, int penalty_order, double SNR_threshold, double peak_shape_threshold)
{

}

void VespucciDataset::IModPolyBaseline(const uword poly_order, const uword max_it, double threshold)
{
    state_changed_ = true;
    SetOldCopies();
    mat baselines(spectra_.n_rows, spectra_.n_cols);
    vec baseline, corrected;
    QProgressDialog *progress = new QProgressDialog(main_window_);
    progress->setMinimum(0);
    progress->setMaximum(spectra_.n_cols-1);
    progress->setWindowTitle("Correcting Baseline");
    progress->setWindowModality(Qt::WindowModal);
    QString colcount = " of " + QString::number(spectra_.n_cols);
    progress->setLabelText("Spectrum 0" + colcount);
    progress->show();
    double err;
    try{
        for (uword i = 0; i < spectra_.n_cols; ++i){
            progress->setValue(i);
            progress->setLabelText("Spectrum " + QString::number(i+1) + colcount);
            Vespucci::Math::LinLeastSq::IModPoly(spectra_.col(i),
                                                 abscissa_, baseline,
                                                 corrected, err,
                                                 poly_order, max_it, threshold);
            baselines.col(i) = baseline;
            spectra_.col(i) = corrected;
            if (progress->wasCanceled()){
                progress->close();
                break;
                Undo();
            }
        }
    }catch(exception e){
          main_window_->DisplayExceptionWarning("VespucciDataset::IModPolyBaseline", e);
    }
    progress->close();
    if (!progress->wasCanceled()){
        //AddAnalysisResult("IModPoly Baselines", baselines);
    }
    operations_ << "IModPolyBaseline("
                   + QString::number(poly_order) + ", "
                   + QString::number(max_it) + ", "
                   + QString::number(threshold) + ")";
    last_operation_ = "IModPoly baseline";
}

void VespucciDataset::RemoveClippedSpectra(double threshold)
{
    SetOldCopies();
    rowvec spectra_max = max(spectra_, 0);

    uvec valid_indices = find(spectra_max < threshold);
    if (valid_indices.n_elem == 0)
        return;

    try{
        spectra_ = spectra_.cols(valid_indices);
        x_ = x_.rows(valid_indices);
        y_ = y_.rows(valid_indices);
    }
    catch(exception e){
        string str = "RemoveClipped Spectra: " + string(e.what());
        throw std::runtime_error(str);
    }

    if (spectra_.n_rows != spectra_old_.n_rows)
        non_spatial_ = true;

    workspace_->UpdateModel();
    last_operation_ = "threshold";
    operations_ << "RemovedClippedSpectra(" + QString::number(threshold) + ")";
}

void VespucciDataset::RemoveFlatSpectra(double threshold)
{
    SetOldCopies();
    rowvec spectra_max = max(spectra_, 0);

    uvec valid_indices = find(spectra_max > threshold);
    if (valid_indices.n_elem == 0)
        return;

    try{
        spectra_ = spectra_.cols(valid_indices);
        x_ = x_.rows(valid_indices);
        y_ = y_.rows(valid_indices);
    }
    catch(exception e){
        string str = "RemoveClipped Spectra: " + string(e.what());
        throw std::runtime_error(str);
    }

    if (spectra_.n_rows != spectra_old_.n_rows)
        non_spatial_ = true;

    workspace_->UpdateModel();
    last_operation_ = "threshold";
    operations_ << "RemoveFlatSpectra(" + QString::number(threshold) + ")";
}

void VespucciDataset::ZeroClippedSpectra(double threshold)
{
    SetOldCopies();
    rowvec spectra_max = max(spectra_, 0);
    uvec invalid_indices = find(spectra_max > threshold);
    if (invalid_indices.n_elem == 0)
        return;

    try{
        for (uword i = 0; i < invalid_indices.n_rows; ++i)
            spectra_.col(invalid_indices(i)) = zeros(spectra_.n_rows);
    }
    catch(exception e){
        string str = "RemoveClipped Spectra: " + string(e.what());
        throw std::runtime_error(str);
    }

    if (spectra_.n_rows != spectra_old_.n_rows)
        non_spatial_ = true;

    workspace_->UpdateModel();
    last_operation_ = "threshold";
    operations_ << "ZeroClippedSpectra(" + QString::number(threshold) + ")";
}

void VespucciDataset::ZeroFlatSpectra(double threshold)
{
    SetOldCopies();
    rowvec spectra_max = max(spectra_, 0);
    uvec invalid_indices = find(spectra_max < threshold);
    if (invalid_indices.n_elem == 0)
        return;

    try{
        for (uword i = 0; i < invalid_indices.n_rows; ++i)
            spectra_.col(invalid_indices(i)) = zeros(spectra_.n_rows);
    }
    catch(exception e){
        string str = "RemoveClipped Spectra: " + string(e.what());
        throw std::runtime_error(str);
    }

    if (spectra_.n_rows != spectra_old_.n_rows)
        non_spatial_ = true;

    workspace_->UpdateModel();
    last_operation_ = "threshold";
    operations_ << "ZeroFlatSpectra(" + QString::number(threshold) + ")";
}

//Filtering functions

void VespucciDataset::MedianFilter(unsigned int window_size)
{
    state_changed_ = true;
    mat processed;
    SetOldCopies();
    try{
        processed = Vespucci::Math::Smoothing::MedianFilterMat(spectra_, window_size);
        spectra_ = processed;
    }
    catch(exception e){
        string str = "MedianFilter: " + string(e.what());
        throw std::runtime_error(str);
    }

    last_operation_ = "median filter";
    operations_ << "MedianFilter(" + QString::number(window_size) + ")";
}


void VespucciDataset::LinearMovingAverage(unsigned int window_size)
{
    state_changed_ = true;
    SetOldCopies();
    try{
        vec filter = Vespucci::Math::Smoothing::CreateMovingAverageFilter(window_size);
        for (uword j = 0; j < spectra_.n_cols; ++j){
            spectra_.col(j) = Vespucci::Math::Smoothing::ApplyFilter(spectra_.col(j), filter);
        }
    }
    catch(exception e){
        string str = "LinearMovingAverage: " + string(e.what());
        throw std::runtime_error(str);
    }

    last_operation_ = "moving average filter";
    operations_ << "LinearMovingAverage(" +QString::number(window_size) + ")";
}

void VespucciDataset::SingularValue(unsigned int singular_values)
{
    state_changed_ = true;
    SetOldCopies();
    mat U;
    vec s;
    mat V;
    try{
        spectra_ = Vespucci::Math::Smoothing::SVDDenoise(spectra_, singular_values, U, s, V);
    }
    catch(exception e){
        string str = "SingularValue: " + string(e.what());
        throw std::runtime_error(str);
    }

    vec SVD_vec({double(singular_values)});
    QSharedPointer<AnalysisResults> results(new AnalysisResults("QUIC-SVD", "QUIC-SVD"));
    results->AddMatrix("Left Singular Vectors", U);
    results->AddMatrix("Right Singular Vectors", V);
    results->AddMatrix("Singular Values", s);
    results->AddMatrix("Rank", SVD_vec);
    AddAnalysisResult(results);

    last_operation_ = "truncated SVD de-noise";
    operations_ << "SingularValue(" + QString::number(singular_values) + ")";
}


int VespucciDataset::QUIC_SVD(double epsilon)
{
    SetOldCopies();
    state_changed_ = true;

    uword rank;
    mat U, V;
    vec s;
    try{
        spectra_ = Vespucci::Math::Smoothing::QUICSVDDenoise(spectra_, epsilon, U, s, V, rank);
    }catch(exception e){
        string str = "QUIC_SVD: " + string(e.what());
        throw std::runtime_error(str);
    }

    vec SVD_vec({double(rank)});

    QSharedPointer<AnalysisResults> results(new AnalysisResults("QUIC-SVD", "QUIC-SVD"));
    results->AddMatrix("Left Singualr Vectors", U);
    results->AddMatrix("Right Singualr Vectors", V);
    results->AddMatrix("Singualr Values", s);
    results->AddMatrix("Rank", SVD_vec);

    last_operation_ = "QUIC_SVD de-noise";
    operations_ << "QUIC_SVD(" + QString::number(epsilon) +  " )";
    return rank;
}

void VespucciDataset::SavitzkyGolay(unsigned int derivative_order,
                                    unsigned int polynomial_order,
                                 unsigned int window_size)
{
    state_changed_ = true;
    SetOldCopies();

    try{
        spectra_ = Vespucci::Math::Smoothing::sgolayfilt(spectra_,
                                        polynomial_order,
                                        window_size,
                                        derivative_order,
                                        1);
    }
    catch(exception e){
        string str = "SavitzkyGolay: " + string(e.what());
        throw std::runtime_error(str);
    }

    last_operation_ = "Savitzky-Golay filtering";
    operations_ << "SavitzkyGolay("
                   + QString::number(derivative_order) + ", "
                   + QString::number(polynomial_order) + ", "
                   + QString::number(window_size) + ")";
}

void VespucciDataset::Scale(double scaling_factor)
{
    state_changed_ = true;
    SetOldCopies();

    try{
        spectra_ = spectra_ * scaling_factor;
    }catch(exception e){
     string str = "Scale(): " + string(e.what());
     throw std::runtime_error(str);
    }

    last_operation_ = "Scaling";
    operations_ << "Scale(" + QString::number(scaling_factor) + ")";
}

void VespucciDataset::ShedSpectrum(const uword index)
{
    state_changed_ = true;

    SetOldCopies();
    try{
        spectra_.shed_col(index);
        x_.shed_row(index);
        y_.shed_row(index);
        indices_.shed_row(index);
    }
    catch(exception e){
        cout << e.what();
        std::runtime_error exc("VespucciDataset::ShedSpectrum");
        main_window_->DisplayExceptionWarning(exc);
    }
    cout << "spectra_ columns (post) = " << spectra_.n_cols;
    operations_ << "ShedSpectrum(" + QString::number(index) + ")";
    workspace_->UpdateModel();
}

void VespucciDataset::ZeroSpectrum(const uword index)
{
    state_changed_ = true;
    SetOldCopies();
    try{
        spectra_.col(index).fill(0);
    }
    catch(exception e){
        cout << e.what();
        std::runtime_error exc("VespucciDataset::ShedSpectrum");
        main_window_->DisplayExceptionWarning(exc);
    }
    operations_ << "ZeroSpectrum(" + QString::number(index) + ")";
    workspace_->UpdateModel();
}

int VespucciDataset::HySime()
{
    cout << "VespucciDataset::HySime\n";

    wall_clock timer;
    mat noise, noise_correlation, subspace;
    cout << "call EstimateAdditiveNoise\n";

    timer.tic();
    Vespucci::Math::DimensionReduction::EstimateAdditiveNoise(noise, noise_correlation, spectra_);
    cout << "took " << timer.toc() << " seconds.\n";

    cout << "Call HySime\n";

    timer.tic();
    int k = Vespucci::Math::DimensionReduction::HySime(spectra_, noise, noise_correlation, subspace);
    cout << "Took " << timer.toc() << " seconds.\n";

    return k;
}

void VespucciDataset::TransformAbscissa(QString input_units, double input_factor, QString output_units, double output_factor, QString description)
{
    state_changed_ = true;
    if (input_units == "Wavelength"){
        if (output_units == "Energy"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::WavelengthToEnergy(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavenumber"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::WavelengthToWavenumber(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Frequency"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::WavelengthToFrequency(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavelength"){
            SetOldCopies();
            abscissa_ = (input_factor * output_factor) * abscissa_;
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else{return;}
    }
    else if (input_units == "Energy"){
        if (output_units == "Energy"){
            SetOldCopies();
            abscissa_ = (input_factor * output_factor) * abscissa_;
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavenumber"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::EnergyToWavenumber(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Frequency"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::EnergyToFrequency(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavelength"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::EnergyToWavelength(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else{return;}
    }
    else if (input_units == "Wavenumber"){
        if (output_units == "Energy"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::WavenumberToEnergy(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavenumber"){
            SetOldCopies();
            abscissa_ = (input_factor * output_factor) * abscissa_;
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Frequency"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::WavenumberToFrequency(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavelength"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::WavenumberToWavelength(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else{return;}
    }
    else if (input_units == "Frequency"){
        if (output_units == "Energy"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::FrequencyToEnergy(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavenumber"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::FrequencyToWavenumber(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Frequency"){
            SetOldCopies();
            abscissa_ = (input_factor * output_factor) * abscissa_;
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else if (output_units == "Wavelength"){
            SetOldCopies();
            abscissa_ = Vespucci::Math::FrequencyToWavelength(abscissa_, output_factor, input_factor);
            x_axis_description_ = description;
            last_operation_ = "Abscissa Transform";
        }
        else{return;}//do nothing for invalid input
    }
    else{return;}
}

void VespucciDataset::InterpolateToNewAbscissa(const vec &new_abscissa, unsigned int poly_order, unsigned int window_size)
{
    state_changed_ = true;
    mat new_spectra;
    try{
        new_spectra = Vespucci::Math::Smoothing::InterpolateToNewAbscissa(spectra_, abscissa_, new_abscissa, window_size, poly_order);
    }catch(exception e){
        string str = "InterpolateToNewAbscissa(vec, uint, uint): " + string(e.what());
        throw std::runtime_error(str);
    }
    SetOldCopies();
    abscissa_ = new_abscissa;
    spectra_ = new_spectra;
    last_operation_ = "Abscissa Interpolation";
    workspace_->UpdateModel();
}

void VespucciDataset::InterpolateToNewAbscissa(const vec &new_abscissa)
{
    state_changed_ = true;
    mat new_spectra;
    try{
        new_spectra = Vespucci::Math::Smoothing::InterpolateToNewAbscissa(spectra_, abscissa_, new_abscissa);
    }catch(exception e){
        string str = "InterpolateToNewAbscissa(vec): " + string(e.what());
        throw std::runtime_error(str);
    }
    SetOldCopies();
    abscissa_ = new_abscissa;
    spectra_ = new_spectra;
    last_operation_ = "Abscissa Interpolation";
    workspace_->UpdateModel();
}

void VespucciDataset::FourierTransform(int n)
{
    state_changed_ = true;
    cx_mat f_spectra(spectra_.n_rows, spectra_.n_cols);
    vec f_abscissa(abscissa_.n_rows);
    try{
        Vespucci::Math::Transform::fft_mat(spectra_, abscissa_,
                                           f_spectra, f_abscissa,
                                           n);
    }catch(exception e){
        string str = "FourierTransform: " + string(e.what());
        throw std::runtime_error(str);
    }
    SetOldCopies();
    spectra_ = real(f_spectra);
    spectra_imag_ = imag(f_spectra);
    abscissa_ = f_abscissa;
    x_axis_description_ = "Frequency (Hz)";
}

void VespucciDataset::InverseFourierTransform(int n)
{
    state_changed_ = true;
    cx_mat t_spectra(spectra_.n_rows, spectra_.n_cols);
    vec t_abscissa(abscissa_.n_rows);
    try{
        Vespucci::Math::Transform::ifft_mat(cx_spectra(),
                                            abscissa_,
                                            t_spectra,
                                            t_abscissa,
                                            n);
    }catch(exception e){
        string str = "InverseFourierTransform: " + string(e.what());
        throw std::runtime_error(str);
    }
    SetOldCopies();
    spectra_ = real(t_spectra);
    spectra_imag_ = imag(t_spectra);
    abscissa_ = t_abscissa;
    x_axis_description_ = "Time (s)";
}

void VespucciDataset::ApplyFTWeight(QString type,
                                    double param)
{
    state_changed_ = true;
    SetOldCopies();
    string shape;
    if (type == "Exponential"){shape = "exp";}
    else{shape = "gaus";}
    try{
        spectra_ = Vespucci::Math::Transform::ApplyWeights(spectra_, abscissa_, shape, param);
    }catch(exception e){
        main_window_->DisplayExceptionWarning("VespucciDataset::ApplyFTWeight", e);
    }
}

void VespucciDataset::ApplyFTWeight(double start_offset,
                                    double end_offset,
                                    double power)
{
    state_changed_ = true;
    SetOldCopies();
    try{
        spectra_ = Vespucci::Math::Transform::ApplySBWeights(spectra_,
                                                             abscissa_,
                                                             start_offset,
                                                             end_offset,
                                                             power);
    }catch(exception e){
        main_window_->DisplayExceptionWarning("VespucciDataset::ApplyFTWeight", e);
    }

}


void VespucciDataset::Univariate(const QString &name,
                                 double &left_bound, double &right_bound,
                                 uword bound_window)
{
    state_changed_ = true;
    QSharedPointer<UnivariateData> univariate_data(new UnivariateData(name));
    try{
        univariate_data->Apply(left_bound, right_bound, bound_window, spectra_, abscissa_);
    }catch(exception e){
        main_window_->DisplayExceptionWarning(e);
        return;
    }
    AddAnalysisResult(univariate_data);
    workspace_->UpdateModel();

    operations_ << "Univariate("
                   + name + ", "
                   + QString::number(left_bound) + ", "
                   + QString::number(right_bound) + ", "
                   + QString::number(bound_window) + ")";

}

void VespucciDataset::FitPeak(const QString &name, const QString &peak_shape, double &left_bound, double &right_bound)
{
    state_changed_ = true;
    QSharedPointer<UnivariateData> univariate_data(new UnivariateData(name));
    try{
        univariate_data->Apply(peak_shape, left_bound, right_bound, spectra_, abscissa_);
    }catch (exception e){
        main_window_->DisplayExceptionWarning(e);
    }
    AddAnalysisResult(univariate_data);
    workspace_->UpdateModel();
}

void VespucciDataset::BandRatio(const QString &name, double &first_left_bound, double &first_right_bound, double &second_left_bound, double &second_right_bound, uword bound_window)
{
    state_changed_ = true;
    QSharedPointer<UnivariateData> univariate_data(new UnivariateData(name));
    try{
        univariate_data->Apply(first_left_bound, first_right_bound,
                               second_left_bound, second_right_bound,
                               bound_window,
                               spectra_, abscissa_);
    }catch(exception e){
        main_window_->DisplayExceptionWarning(e);
        return;
    }
    AddAnalysisResult(univariate_data);
    workspace_->UpdateModel();
    operations_ << "BandRatio("
                   + name + ", "
                   + QString::number(first_left_bound) + ", "
                   + QString::number(first_right_bound) + ", "
                   + QString::number(second_left_bound) + ", "
                   + QString::number(second_right_bound) + ", "
                   + QString::number(bound_window) + ")";


}


void VespucciDataset::PrincipalComponents(const QString &name, bool scale_data)
{
    QSharedPointer<MlpackPCAData> mlpack_pca_data(new MlpackPCAData(name));
    state_changed_ = true;

    try{
        mlpack_pca_data->Apply(spectra_, scale_data);
    }
    catch(exception e){
        throw std::runtime_error("VespucciDataset::PrincipalComponents");
    }
    AddAnalysisResult(mlpack_pca_data);
    workspace_->UpdateModel();

}

void VespucciDataset::PrincipalComponents(const QString &name)
{
    state_changed_ = true;

    QSharedPointer<PrincipalComponentsData> pca_data(new PrincipalComponentsData(name));
    try{
        pca_data->Apply(spectra_);
    }catch(exception e){
        string str = "PrincipalComponents: " + string(e.what());
        throw std::runtime_error(str);
    }
    AddAnalysisResult(pca_data);
    workspace_->UpdateModel();
    operations_ << "PrincipalComponents(" + name + ")";


}

void VespucciDataset::FindPeaks(const QString &name, double sel, double threshold, uword poly_order, uword window_size)
{/*

    mat peak_magnitudes;
    mat peak_positions;


    try{
        peak_positions =
                Vespucci::Math::PeakFinding::FindPeaksMat(spectra_,
                                       sel, threshold,
                                       poly_order, window_size,
                                       peak_magnitudes);
    }catch(exception e){
        cerr << e.what();
        throw std::runtime_error("FindPeaks");
    }
    QSharedPointer<AnalysisResults> peak_mag(new AnalysisResults(peak_magnitudes));
    QSharedPointer<AnalysisResults> peak_pos(new AnalysisResults(peak_positions));
    analysis_results_.insert("Peak Positions", peak_pos);
    analysis_results_.insert("Peak Magnitudes", peak_mag);

*/}

void VespucciDataset::AgglomerativeClustering(const QString &name, const QString &linkage, const QString &metric)
{
    QSharedPointer<AnalysisResults> ahca_results(new AnalysisResults(name, "AHCA"));
    mat assignments;
    Vespucci::Math::Clustering::AHCA ahca;
    try{
        ahca.SetLinkage(linkage.toStdString());
        ahca.SetMetric(metric.toStdString());
        ahca.Link(spectra_);
        assignments = ahca.Cluster(spectra_.n_cols);
    }catch (exception e){
        string str = "AgglomerativeClustering: " + string(e.what());
        throw runtime_error(str);
    }

    ahca_results->AddMatrix("Assignments", assignments);
    ahca_results->AddMatrix("Spectrum Distances", ahca.dist());
    ahca_results->AddMatrix("Cluster Distances", ahca.merge_data());
    AddAnalysisResult(ahca_results);
    state_changed_ = true;
    operations_ << "AgglomerativeClustering(" + name + ", " + linkage + ", " + metric + ")";
    workspace_->UpdateModel();
}

void VespucciDataset::CalculateRepresentativeSpectrum(const QString &name, QString statistic, QString metric)
{
    uword index;
    vec rep;
    try{
        rep = Vespucci::Math::RepresentativeSpectrum(spectra_, index, metric.toStdString(), statistic.toStdString());
    }catch(exception e){
        string str = "CalculateRepresentativeSpectrum: " + string(e.what());
        throw runtime_error(str);
    }
    QString matrix_name = name;
    int i = 1;
    while (auxiliary_matrices_->HasMatrix(matrix_name))
        matrix_name = name + " (" + QString::number(i++) + ")";
    auxiliary_matrices_->AddMatrix(matrix_name, rep);
    workspace_->UpdateModel();
}


void VespucciDataset::VertexComponents(const QString &name, uword endmembers)
{
    state_changed_ = true;

    if(endmembers == 0){
        QMessageBox alert;
        alert.setText("The HySime algorithm will take a long time.");
        alert.setInformativeText("OK to continue");

        alert.setStandardButtons(QMessageBox::Ok | QMessageBox::Cancel);
        alert.setWindowTitle("Principal Components Analysis");
        alert.setIcon(QMessageBox::Question);

        int ret = alert.exec();


        if (ret == QMessageBox::Cancel){
            return;
        }
        endmembers = HySime();
    }

    QSharedPointer<VCAData> vertex_components_data(new VCAData(name));

    try{
        vertex_components_data->Apply(spectra_, endmembers);
    }catch(exception e){
        cerr << "VespucciDataset::VertexComponents()\n";

        throw std::runtime_error("VertexComponents: " + string(e.what()));
    }
    AddAnalysisResult(vertex_components_data);
    workspace_->UpdateModel();
    operations_ << "VertexComponents("
                   + name + ", "
                   + QString::number(endmembers) + ")";
}




void VespucciDataset::PartialLeastSquares(const QString &name, uword components)
{
    state_changed_ = true;
    if(components == 0){
        components = HySime();
    }

    QSharedPointer<PLSData> pls_data(new PLSData(name));


    try{
        pls_data->Classify(spectra_, abscissa_, components);
    }catch(exception e){
        string str = "PartialLeastSquares: " + string(e.what());
        throw std::runtime_error(str);
    }

    AddAnalysisResult(pls_data);
    workspace_->UpdateModel();
    operations_ << "PartialLeastSquares("
                   + name + ", "
                   + QString::number(components) + ")";

}

void VespucciDataset::PLSCalibration(const QString &name, const QStringList &control_keys)
{
    state_changed_ = true;
    QSharedPointer<PLSData> pls_data(new PLSData(name));
    mat controls = workspace_->GetMatrix(control_keys);

    try{
        pls_data->Calibrate(spectra_, controls);
    }catch(exception e){
        string str = "PartialLeastSquares: " + string(e.what());
        throw std::runtime_error(str);
    }

    AddAnalysisResult(pls_data);
    workspace_->UpdateModel();
}

void VespucciDataset::TrainPLSDA(const QString &name, const QStringList &label_keys)
{
    state_changed_ = true;
    QSharedPointer<PLSData> pls_data(new PLSData(name));
    mat labels = workspace_->GetMatrix(label_keys);

    try{
        pls_data->Discriminate(spectra_, labels);
    }catch(exception e){
        string str = "PartialLeastSquares: " + string(e.what());
        throw std::runtime_error(str);
    }

    AddAnalysisResult(pls_data);
    workspace_->UpdateModel();
}

void VespucciDataset::CorrelationAnalysis(const QString &control_key, QString name)
{

    if (!auxiliary_matrices_->HasMatrix(control_key)){
        main_window_->DisplayWarning("Object does not exit", "The object at the specified key does not exist!");
        return;
    }
    vec control;
    try{
        control = auxiliary_matrices_->GetMatrix(control_key).col(0);
    }catch(exception e){
        main_window_->DisplayExceptionWarning("VespucciDataset::CorrelationAnalysis", e);
        return;
    }

    QSharedPointer<UnivariateData> univariate_data(new UnivariateData(name));
    try{
        univariate_data->ApplyCorrelation(spectra_, control);
    }catch(exception e){
        main_window_->DisplayExceptionWarning("UnivariateData::apply", e);
        return;
    }

    state_changed_ = true;
    AddAnalysisResult(univariate_data);
    workspace_->UpdateModel();
}


void VespucciDataset::KMeans(const QString &name, const QString &metric_text, const QString &partition_policy, bool allow_empty, size_t clusters)
{
    mat centroids;
    vec assignments;
    try{
        Vespucci::Math::KMeansWrapper k(partition_policy.toStdString(), metric_text.toStdString(), allow_empty);
        assignments = k.Cluster(spectra_, clusters, centroids);
    }catch(exception e){
        main_window_->DisplayExceptionWarning("KMeansWrapper: ", e);
        return;
    }

    QSharedPointer<AnalysisResults> results(new AnalysisResults(name, "k-Means Analysis"));
    results->AddMatrix("Assignments", assignments);
    results->AddMatrix("Centroids", centroids);
    AddAnalysisResult(results);
    workspace_->UpdateModel();
    operations_ << "KMeans("
                   + name + ", "
                   + QString::number(clusters) + ", "
                   + metric_text + ", "
                   + partition_policy + ", "
                   + (allow_empty ? "true" : "false" ) + ")";
}

void VespucciDataset::ClassicalLeastSquares(const QString &name, const QStringList &reference_keys)
{
    state_changed_ = true;
    mat reference = workspace_->GetMatrix(reference_keys);
    mat coefs;
    try{
        coefs = Vespucci::Math::LinLeastSq::OrdinaryLeastSquares(reference, spectra_);
    }catch(exception e){
        string str = "Vespucci::Math::LinLeastSq::OrdinaryLeastSquares" + string(e.what());
        throw runtime_error(str);
    }
    coefs = coefs.t(); //to make coefficients mappable
    QSharedPointer<AnalysisResults> results(new AnalysisResults(name, "CLS Analysis"));
    results->AddMatrix("Coefficients", coefs);
    AddAnalysisResult(results);
    workspace_->UpdateModel();
    QString operation = "ClassicalLeastSquares(" + name;
    for (auto param: reference_keys)
        operation = operation + ", " + param;
    operation = operation + ")";
    operations_ << operation;
}


// HELPER FUNCTIONS //

uvec VespucciDataset::FindRange(double start, double end) const
{
    uvec indices(2);
    indices(0) = FindIndex(start);
    indices(1) = FindIndex(end);
    return indices;
}

uword VespucciDataset::FindOrigin() const
{
    double delta = std::max(std::abs((x_(1) - x(0))), std::abs((y_(1) - y_(0))));
    uvec zero_x = find(((0-delta) <= x_) && (0 + delta) >= x_);
    vec sub_y = y_.elem(zero_x);
    uvec zero_y = find(((0-delta) <= sub_y) && ((0+delta) >= sub_y));

    // If for some reason this doesn't work, also find the point halfway down.
    uword mid = ((x_.n_rows % 2 == 0) ? x_.n_rows : (2*(x_.n_rows + 1)) / 2);
    return (zero_y.n_rows > 0 ? zero_y(0) : mid);
}

vec VespucciDataset::PointSpectrum(uword index) const
{
    return spectra_.col((index < spectra_.n_cols ? index : spectra_.n_cols - 1));
}

QVector<double> VespucciDataset::WavelengthQVector() const
{
    std::vector<double> abscissa_stdvector =
            conv_to< std::vector<double> >::from(abscissa_);

    QVector<double> abscissa_qvector =
            QVector<double>::fromStdVector(abscissa_stdvector);

    return abscissa_qvector;
}

uword VespucciDataset::FindIndex(double abscissa_value) const
{
    double delta = std::fabs(abscissa_(1) - abscissa_(0));
    uvec indices = find((abscissa_value - delta) < abscissa_ <= (abscissa_value + delta));
    return indices(0);
}

QCPRange VespucciDataset::ValueRange() const
{
    double lower = y_.min();
    double upper = y_.max();
    QCPRange range(upper, lower);
    return range;
}
QCPRange VespucciDataset::KeyRange() const
{
    double lower = x_.min();
    double upper = x_.max();
    QCPRange range(upper, lower);
    return range;
}

int VespucciDataset::KeySize() const
{
    uword i;
    uword x_count=1;
    double x_buf;

    //loop through x until a value different then the first is met.
    if (!flipped_){
        x_count = 1; //this counts the first entry in x_
        x_buf = x_(0);
        for(i=0; i<x_.n_elem; ++i){
            if(x_(i)!=x_buf){
                ++x_count;
                x_buf=x_(i);
            }
        }
    } else{
        x_count = 0;
        for (i=0; i<x_.n_elem; ++i){
            if(y_(i)!=y_(0)){
                break;
            } else{
                ++x_count;
            }
        }
    }

    return x_count;
}

int VespucciDataset::ValueSize() const
{

    uword i = 0;
    uword y_count;


    //long-text files hold x constant and vary y
    //until x is different, count y
    //reverse if flipped
    if (!flipped_){
        y_count = 0;
        for (i=0; i<x_.n_elem; ++i){
            if(x_(i)!=x_(0)){
                break;
            }
            else{
                ++y_count;
            }
        }
    } else{
        y_count = 1;
        double y_buf = y_(0);
        for(i=0; i<y_.n_elem; ++i){
            if(y_(i)!=y_buf){
                ++y_count;
                y_buf=y_(i);
            }
        }
    }


    return y_count;
}


// MEMBER ACCESS FUNCTIONS //
vec VespucciDataset::wavelength() const
{
    return abscissa_;
}

vec VespucciDataset::abscissa() const
{
    return abscissa_;
}

vec VespucciDataset::wavelength(uvec indices) const
{
    return abscissa_.rows(indices);
}

vec VespucciDataset::x() const
{
    return x_;
}

vec VespucciDataset::indices() const
{
    return indices_;
}

mat *VespucciDataset::indices_ptr()
{
    return (mat *) &indices_;
}

QStringList VespucciDataset::operations()
{
    return operations_;
}

void VespucciDataset::SetIndices(vec indices)
{
    indices_ = indices;
}

vec VespucciDataset::x(uvec indices) const
{
    return x_(indices);
}

double VespucciDataset::x(uword index) const
{
    if (index >= x_.n_rows)
        return x_(x_.n_rows - 1);
    else
        return x_(index);
}



vec VespucciDataset::y() const
{
    return y_;
}

vec VespucciDataset::y(uvec indices) const
{
    return y_(indices);
}

double VespucciDataset::y(uword index) const
{
    if (index >= y_.n_rows)
        return y_(y_.n_rows - 1);
    else
        return y_(index);
}

mat VespucciDataset::spectra() const
{
    return spectra_;
}

cx_mat VespucciDataset::cx_spectra() const
{
    if (spectra_imag_.n_rows == spectra_.n_rows && spectra_imag_.n_cols == spectra_.n_cols)
        return cx_mat(spectra_, spectra_imag_);
    else
        return cx_mat(spectra_, zeros(spectra_.n_rows, spectra_.n_cols));
}

cx_mat VespucciDataset::cx_spectra(uvec indices) const
{
    if (spectra_imag_.n_rows == spectra_.n_rows && spectra_imag_.n_cols == spectra_.n_cols)
        return cx_mat(spectra_.cols(indices), spectra_imag_.cols(indices));
    else
        return cx_mat(spectra_.cols(indices), zeros(spectra_.n_rows, indices.n_rows));
}

mat VespucciDataset::spectra(uvec indices) const
{
    return spectra_.cols(indices);
}

const QString VespucciDataset::name() const
{
    return name_;
}

void VespucciDataset::SetName(QString new_name)
{
    name_ = new_name;
}

void VespucciDataset::SetData(const mat &spectra, const vec &wavelength, const vec &x, const vec &y)
{
    spectra_ = spectra;
    abscissa_ = wavelength;
    x_ = x;
    y_ = y;
}

//MAP HANDLING FUNCTIONS
//QStringList VespucciDataset::map_names()
//{
//    return map_names_;
//}

int VespucciDataset::map_loading_count() const
{
    return map_loading_count_;
}


void VespucciDataset::AddMap(QSharedPointer<MapData> map)
{
    QString basename = map->name();
    int i = 0;
    while(MapKeys().contains(map->name()))
        map->SetName(basename + " (" + QString::number(++i) + ")", map->type());
    maps_.append(map);
    workspace_->UpdateModel();
}

void VespucciDataset::RemoveMap(const QString &name)
{
    if (MapKeys().contains(name)){
        int i = maps_.size();
        while (i--){
            if (maps_[i]->name() == name)
                maps_.removeAt(i); //might be safe to break the while here
        }
    }
}

QCPRange VespucciDataset::WavelengthRange() const
{
    double min = abscissa_.min();
    double max = abscissa_.max();
    QCPRange range(min, max);
    return range;
}

QCPRange VespucciDataset::PointSpectrumRange(int i) const
{
    vec row = spectra_.col(i);
    double min = row.min();
    double max = row.max();

    QCPRange range(min, max);
    return range;
}

QCPColorGradient VespucciDataset::GetGradient(int gradient_number) const
{
    switch (gradient_number)
    {
    case 0: return QCPColorGradient::cbBuGn;
    case 1: return QCPColorGradient::cbBuPu;
    case 2: return QCPColorGradient::cbGnBu;
    case 3: return QCPColorGradient::cbOrRd;
    case 4: return QCPColorGradient::cbPuBu;
    case 5: return QCPColorGradient::cbPuBuGn;
    case 6: return QCPColorGradient::cbPuRd;
    case 7: return QCPColorGradient::cbRdPu;
    case 8: return QCPColorGradient::cbYlGn;
    case 9: return QCPColorGradient::cbYlGnBu;
    case 10: return QCPColorGradient::cbYlOrBr;
    case 11: return QCPColorGradient::cbYlOrRd;
    case 12: return QCPColorGradient::cbBlues;
    case 13: return QCPColorGradient::cbGreens;
    case 14: return QCPColorGradient::cbOranges;
    case 15: return QCPColorGradient::cbPurples;
    case 16: return QCPColorGradient::cbReds;
    case 17: return QCPColorGradient::cbGreys;
    case 18: return QCPColorGradient::gpGrayscale;
    case 19: return QCPColorGradient::gpNight;
    case 20: return QCPColorGradient::gpCandy;
    case 21: return QCPColorGradient::gpIon;
    case 22: return QCPColorGradient::gpThermal;
    case 23: return QCPColorGradient::gpPolar;
    case 24: return QCPColorGradient::gpSpectrum;
    case 25: return QCPColorGradient::gpJet;
    case 26: return QCPColorGradient::gpHues;
    case 27: return QCPColorGradient::gpHot;
    case 28: return QCPColorGradient::gpCold;
    case 29: return QCPColorGradient::cbBrBG;
    case 30: return QCPColorGradient::cbPiYG;
    case 31: return QCPColorGradient::cbPRGn;
    case 32: return QCPColorGradient::cbPuOr;
    case 33: return QCPColorGradient::cbRdBu;
    case 34: return QCPColorGradient::cbRdGy;
    case 35: return QCPColorGradient::cbRdYlBu;
    case 36: return QCPColorGradient::cbRdYlGn;
    case 37: return QCPColorGradient::cbSpectral;
    case 38: return QCPColorGradient::vSpectral;
    default: return QCPColorGradient::gpCold;
    }
}

QCPColorGradient VespucciDataset::GetClusterGradient(int clusters) const
{
    switch (clusters)
    {
    case 2: return QCPColorGradient::cbCluster2;
    case 3: return QCPColorGradient::cbCluster3;
    case 4: return QCPColorGradient::cbCluster4;
    case 5: return QCPColorGradient::cbCluster5;
    case 6: return QCPColorGradient::cbCluster6;
    case 7: return QCPColorGradient::cbCluster7;
    case 8: return QCPColorGradient::cbCluster8;
    case 9: return QCPColorGradient::cbCluster9;
    default: return QCPColorGradient::gpJet; //not as good, but should still be different
    }
}

bool VespucciDataset::ConstructorCancelled() const
{
    return constructor_canceled_;
}

mat VespucciDataset::AverageSpectrum(bool stats) const
{
    mat spec_mean = mean(spectra_, 1);
    vec spec_stddev;
    //insert stddevs on next line if requested
    if (stats){
        spec_stddev = stddev(spectra_, 0, 1);
        spec_mean.insert_cols(1, spec_stddev);
    }
    return spec_mean;
}



const QString VespucciDataset::x_axis_description() const
{
    return x_axis_description_;
}

void VespucciDataset::SetXDescription(QString description)
{
    x_axis_description_ = description;
}

void VespucciDataset::SetYDescription(QString description)
{
    y_axis_description_ = description;
}

const QString VespucciDataset::y_axis_description() const
{
    return y_axis_description_;
}


mat* VespucciDataset::spectra_ptr()
{
    return &spectra_;
}

const mat &VespucciDataset::spectra_ref()
{
    return spectra_;
}

const vec &VespucciDataset::abscissa_ref()
{
    return abscissa_;
}

const vec &VespucciDataset::x_ref()
{
    return x_;
}

const vec &VespucciDataset::y_ref()
{
    return y_;
}

size_t VespucciDataset::columns() const
{
    return spectra_.n_cols;
}

double VespucciDataset::AbscissaMin() const
{
    return abscissa_.min();
}

double VespucciDataset::AbscissaMax() const
{
    return abscissa_.max();
}


bool VespucciDataset::Undoable() const
{
    return (spectra_old_.n_elem > 0 ? true : false);
}

uword VespucciDataset::UniqueX() const
{
    vec unique_x = unique(x_);
    return unique_x.n_rows;
}

uword VespucciDataset::UniqueY() const
{
    vec unique_y = unique(y_);
    return unique_y.n_rows;
}


int VespucciDataset::UnivariateCount() const
{
    return analysis_results_.size();
}

mat* VespucciDataset::abscissa_ptr()
{
    return &abscissa_;
}

mat* VespucciDataset::wavelength_ptr()
{
    return &abscissa_;
}

mat* VespucciDataset::x_ptr()
{
    return &x_;
}

mat* VespucciDataset::y_ptr()
{
    return &y_;
}

bool VespucciDataset::non_spatial() const
{
    return non_spatial_;
}

bool VespucciDataset::meta() const
{
    return meta_;
}

void VespucciDataset::SetParentDatasetIndices(mat parent_dataset_indices)
{
    parent_dataset_indices_ = parent_dataset_indices;
}

mat *VespucciDataset::parent_dataset_indices()
{
    return &parent_dataset_indices_;
}

const QString VespucciDataset::last_operation() const
{
    return last_operation_;
}


void VespucciDataset::AddAnalysisResult(QSharedPointer<AnalysisResults> analysis_result)
{
    QString name = analysis_result->name();
    QString new_name = name;
    int i = 0;
    while (AnalysisResultsKeys().contains(new_name))
        new_name = name + " " + QString::number(++i);
    analysis_result->SetName(new_name);
    analysis_result->AddParent(name_, 0, spectra_.n_cols);
    analysis_results_.append(analysis_result);
    state_changed_ = true;
    workspace_->UpdateModel();
}

void VespucciDataset::AddAnalysisResult(QSharedPointer<AnalysisResults> analysis_result, uword start_row, uword end_row)
{
    QString name = analysis_result->name();
    QString new_name = name;
    int i = 0;
    while (AnalysisResultsKeys().contains(new_name))
        new_name = name + " " + QString::number(++i);
    analysis_result->SetName(new_name);
    analysis_result->AddParent(name_, start_row, end_row);
    analysis_results_.append(analysis_result);
    state_changed_ = true;
    workspace_->UpdateModel();
}


QStringList VespucciDataset::AnalysisResultsKeys() const
{
    QStringList keys;
    for (auto result: analysis_results_)
        keys << result->name();
    return keys;
}

QMap<QString, QStringList> VespucciDataset::AnalysisResultsTreeStructure() const
{
    QMap<QString, QStringList> tree_structure;
    for (auto result: analysis_results_)
        tree_structure[result->name()] = result->KeyList();
    return tree_structure;
}

void VespucciDataset::ImportAuxiliaryMatrix(const QString &name, const QString &filename)
{
    mat matrix;
    bool ok = matrix.load(filename.toStdString());
    if (ok)
        auxiliary_matrices_->AddMatrix(name, matrix);
    else
        main_window_->DisplayWarning("Matrix Not Loaded", "The file " + name + " could not be loaded");
}

void VespucciDataset::AddAuxiliaryMatrix(const QString &name, mat &matrix)
{
    QString basename = name;
    QString new_name = name;
    int i = 0;
    while (auxiliary_matrices_->HasMatrix(new_name))
        new_name = basename + " (" + QString::number(++i) + ")";
    auxiliary_matrices_->AddMatrix(new_name, matrix);
    workspace_->UpdateModel();
    state_changed_ = true;
}

void VespucciDataset::AddMatrix(const QString &name, mat &matrix)
{
    if (name == "Spectra")
        spectra_ = matrix;
    else if (name == "x")
        x_ = matrix.col(0);
    else if (name == "y")
        y_ = matrix.col(0);
    else if (name == "Spectral Abscissa")
        abscissa_ = matrix.col(0);
    else
        AddAuxiliaryMatrix(name, matrix);
    state_changed_ = true;
    workspace_->UpdateModel();
}

QStringList VespucciDataset::AuxiliaryMatrixKeys() const
{
    return auxiliary_matrices_->KeyList();
}

QStringList VespucciDataset::CoreMatrixKeys() const
{
    return {"Spectra", "Spectral Abscissa", "x", "y"};
}

const mat & VespucciDataset::GetAnalysisResultMatrix(const QString &results_key, const QString &matrix_key) const
{
    for (auto result: analysis_results_){
        if (result->name() == results_key)
            if (result->KeyList().contains(matrix_key))
                return result->GetMatrix(matrix_key);
    }

    return empty_matrix_;
}

QSharedPointer<AnalysisResults> VespucciDataset::GetAnalysisResult(const QString &key)
{

    for (auto result: analysis_results_){
        if (result->name() == key)
            return result;
    }

    QSharedPointer<AnalysisResults> dummy(new AnalysisResults("Dummy", "Dummy"));
    return dummy;
}

const mat & VespucciDataset::GetAuxiliaryMatrix(const QString &key) const
{
    if (auxiliary_matrices_->HasMatrix(key)) return auxiliary_matrices_->GetMatrix(key);
    else
        return empty_matrix_;
}

const mat &VespucciDataset::GetCoreMatrix(const QString &key) const
{
    if (key == "Spectra"){return spectra_;}
    else if (key == "Spectral Abscissa"){return abscissa_;}
    else if (key == "x"){return x_;}
    else if (key == "y"){return y_;}
    else{return empty_matrix_;}
}

bool VespucciDataset::IsCoreMatrix(const QString &key) const
{
    return (key == "Spectra" || key == "Spectral Abscissa" || key == "x" || key == "y");
}

QSharedPointer<MapData> VespucciDataset::GetMapData(const QString &key)
{
    for (auto map: maps_){
        if (map->name() == key) return map;
    }
    return QSharedPointer<MapData>(0); //return null mapdata
}

void VespucciDataset::CreateMap(const QString &map_name,
                                const QString &results_key,
                                const QString &matrix_key,
                                uword column,
                                QCPColorGradient gradient)
{
    if (AnalysisResultsKeys().contains(results_key)
            && GetAnalysisResult(results_key)->HasMatrix(matrix_key)){
        QStringList data_keys = {name_, results_key, matrix_key};
        QString map_type = GetAnalysisResult(results_key)->type() + " Color Map";
        QSharedPointer<MapData> new_map(new MapData(map_name,
                                                    map_type,
                                                    data_keys,
                                                    column,
                                                    workspace_));

        new_map->setGradient(gradient);
        AddMap(new_map);
        workspace_->UpdateModel();
    }
}

void VespucciDataset::CreateMap(const QString &map_name,
                                const QString &matrix_key,
                                uword column,
                                QCPColorGradient gradient)
{
    if (auxiliary_matrices_->HasMatrix(matrix_key)){
        QStringList data_keys = {name_, matrix_key};
        QSharedPointer<MapData> new_map(new MapData(map_name,
                                                    QString(),
                                                    data_keys,
                                                    column,
                                                    workspace_));
        new_map->setGradient(gradient);
        AddMap(new_map);
        workspace_->UpdateModel();
    }
}

bool VespucciDataset::ShowMapViewer(const QString &map_key, bool show)
{
    if (!MapKeys().contains(map_key)) return false;
    bool ok = MapKeys().contains(map_key);
        if (ok){
        for (auto map: maps_){
            if (map->name() == map_key){
                map->ShowMapWindow(show);
            }
        }
    }
    return ok;
}

QStringList VespucciDataset::MapKeys() const
{
    QStringList map_keys;
    for (auto map: maps_)
        map_keys << map->name();
    return map_keys;
}