From 52d300269c395f3afac63a22de1de971118ace63 Mon Sep 17 00:00:00 2001 From: cameronosmith Date: Thu, 16 Apr 2020 20:31:17 -0700 Subject: [PATCH 1/5] No significant code, just example class structure and basic project layout. --- src/AuxFunc.cc | 379 -------------------- src/AuxFunc.h | 52 --- src/CCM.cc | 808 ------------------------------------------- src/Common.cc | 192 ---------- src/Common.h | 372 +------------------- src/DateTime.h | 32 -- src/DateTimeUtil.cc | 172 --------- src/EDM.cc | 21 ++ src/EDM.h | 35 ++ src/EDM_Functions.cc | 39 +++ src/Embed.cc | 183 ---------- src/Embed.h | 36 -- src/Eval.cc | 687 ------------------------------------ src/Interface.cc | 7 - src/Multiview.cc | 769 ---------------------------------------- src/Neighbors.cc | 277 --------------- src/Neighbors.h | 31 -- src/Parameter.cc | 625 --------------------------------- src/Parameter.h | 122 ------- src/SMap.cc | 513 --------------------------- src/Simplex.cc | 243 +------------ src/Simplex.h | 21 ++ src/Version.h | 22 -- src/makefile | 60 +--- src/makefile.windows | 74 ---- 25 files changed, 134 insertions(+), 5638 deletions(-) delete mode 100644 src/AuxFunc.cc delete mode 100644 src/AuxFunc.h delete mode 100644 src/CCM.cc delete mode 100644 src/Common.cc delete mode 100644 src/DateTime.h delete mode 100644 src/DateTimeUtil.cc create mode 100644 src/EDM.cc create mode 100644 src/EDM.h create mode 100644 src/EDM_Functions.cc delete mode 100644 src/Embed.cc delete mode 100644 src/Embed.h delete mode 100644 src/Eval.cc delete mode 100644 src/Interface.cc delete mode 100644 src/Multiview.cc delete mode 100644 src/Neighbors.cc delete mode 100644 src/Neighbors.h delete mode 100644 src/Parameter.cc delete mode 100644 src/Parameter.h delete mode 100644 src/SMap.cc create mode 100644 src/Simplex.h delete mode 100644 src/Version.h delete mode 100644 src/makefile.windows diff --git a/src/AuxFunc.cc b/src/AuxFunc.cc deleted file mode 100644 index 52f2f62..0000000 --- a/src/AuxFunc.cc +++ /dev/null @@ -1,379 +0,0 @@ - -#include "AuxFunc.h" -#include "DateTime.h" - -namespace EDM_AuxFunc { - std::mutex mtx; -} - -//--------------------------------------------------------------- -// Common code for Simplex and Smap: -// 1) Extract or Embed() data into dataBlock -// 2) Get target (library) vector -// 3) DeletePartialDataRows() -// 4) Adjust param.library and param.prediction indices -// 5) FindNeighbors() -// -// NOTE: time column is not returned in the embedding dataBlock. -// -// NOTE: If data is embedded by Embed(), the returned dataBlock -// has tau * (E-1) fewer rows than data. Since data is -// included in the returned DataEmbedNN struct, the first -// (or last) tau * (E-1) data rows are deleted to match -// dataBlock. The target vector is also reduced. -// -// NOTE: If rows are deleted, then the library and prediction -// vectors in Parameters are updated to reflect this. -//--------------------------------------------------------------- -DataEmbedNN EmbedNN( DataFrame *data, - Parameters ¶m, - bool checkDataRows ) -{ - DataFrame &dataIn = std::ref( *data ); - - if ( checkDataRows ) { - CheckDataRows( param, dataIn, "EmbedNN: Input data" ); - } - - //---------------------------------------------------------- - // Extract or embed dataIn into dataBlock - //---------------------------------------------------------- - DataFrame dataBlock; // Multivariate or embedded DataFrame - - if ( param.embedded ) { - // dataIn is a multivariable block, no embedding needed - // Select the specified columns into dataBlock - if ( param.columnNames.size() ) { - dataBlock = dataIn.DataFrameFromColumnNames( param.columnNames ); - } - else if ( param.columnIndex.size() ) { - dataBlock = dataIn.DataFrameFromColumnIndex( param.columnIndex ); - } - else { - throw std::runtime_error( "EmbedNN(): colNames and " - " colIndex are empty.\n" ); - } - } - else { - // embedded = false: Create the embedding dataBlock via Embed() - // dataBlock will have tau * (E-1) fewer rows than dataIn - dataBlock = Embed( dataIn, param.E, param.tau, - param.columns_str, param.verbose ); - } - - //---------------------------------------------------------- - // Get target (library) vector - //---------------------------------------------------------- - std::valarray targetIn; - if ( param.targetIndex ) { - targetIn = dataIn.Column( param.targetIndex ); - } - else if ( param.targetName.size() ) { - targetIn = dataIn.VectorColumnName( param.targetName ); - } - else { - // Default to first column - targetIn = dataIn.Column( 0 ); - } - - //------------------------------------------------------------ - // embedded = false: Embed() was called on dataIn - // Remove target, dataIn rows as needed - // Adjust param.library and param.prediction indices - //------------------------------------------------------------ - if ( not param.embedded ) { - - if ( param.E < 1 ) { - std::stringstream errMsg; - errMsg << "EmbedNN(): E = " << param.E << " is invalid.\n" ; - throw std::runtime_error( errMsg.str() ); - } - - size_t shift = abs( param.tau ) * ( param.E - 1 ); - - // Copy targetIn excluding partial data into targetEmbed - std::valarray targetEmbed( dataIn.NRows() - shift ); - - // Bogus cast to ( std::valarray ) for MSVC - // as it doesn't export its own slice_array applied to [] - if ( param.tau < 0 ) { - targetEmbed = ( std::valarray ) - targetIn[ std::slice( shift, targetIn.size() - shift, 1 ) ]; - } - else { - targetEmbed = ( std::valarray ) - targetIn[ std::slice( 0, targetIn.size() - shift, 1 ) ]; - } - - // Resize targetIn to ignore partial data rows - targetIn.resize( targetEmbed.size() ); - - // Copy target without partial data into resized targetIn - std::slice targetEmbed_i = std::slice( 0, targetEmbed.size(), 1 ); - targetIn[ targetEmbed_i ] = ( std::valarray ) - targetEmbed[ targetEmbed_i ]; - - // Delete dataIn top or bottom rows of partial data - if ( not dataIn.PartialDataRowsDeleted() ) { - // Not thread safe - std::lock_guard lck( EDM_AuxFunc::mtx ); - - dataIn.DeletePartialDataRows( shift, param.tau ); - } - - // Adjust param.library and param.prediction vectors of indices - if ( shift > 0 ) { - param.DeleteLibPred( shift ); - } - - // Check boundaries again since rows were removed - if ( checkDataRows ) { - CheckDataRows( param, dataIn, "EmbedNN: Embedded data" ); - } - } - - //---------------------------------------------------------- - // Nearest neighbors - //---------------------------------------------------------- - Neighbors neighbors = FindNeighbors( dataBlock, param ); - - // Create struct to return the objects - DataEmbedNN dataEmbedNN = DataEmbedNN( &dataIn, dataBlock, - targetIn, neighbors ); - return dataEmbedNN; -} - -//---------------------------------------------------------- -// Common code for Simplex and Smap output generation -//---------------------------------------------------------- -DataFrame FormatOutput( Parameters param, - std::valarray predictions, - std::valarray const_predictions, - std::valarray variance, - std::valarray target_vec, - std::vector time, - std::string timeName ) -{ - - //---------------------------------------------------- - // Time vector with additional Tp points - //---------------------------------------------------- - size_t N_time = time.size(); - size_t N_row = predictions.size(); - - // Populate vector of time strings for output - std::vector timeOut( N_row + param.Tp ); - - if ( N_time ) { - FillTimes( param, time, std::ref( timeOut ) ); - } - - //---------------------------------------------------- - // Observations: add Tp nan at end - //---------------------------------------------------- - std::valarray observations( N_row + param.Tp ); - std::slice pred_i = std::slice( param.prediction[0], N_row, 1 ); - - observations[ std::slice( 0, N_row, 1 ) ] = - ( std::valarray ) target_vec[ pred_i ]; - - for ( size_t i = N_row; i < N_row + param.Tp; i++ ) { - observations[ i ] = NAN; - } - - //---------------------------------------------------- - // Predictions & variance: insert Tp nan at start - //---------------------------------------------------- - std::valarray predictionsOut( N_row + param.Tp ); - for ( size_t i = 0; i < param.Tp; i++ ) { - predictionsOut[ i ] = NAN; - } - predictionsOut[ std::slice(param.Tp, N_row, 1) ] = predictions; - - std::valarray constPredictionsOut( N_row + param.Tp ); - if ( param.const_predict ) { - for ( size_t i = 0; i < param.Tp; i++ ) { - constPredictionsOut[ i ] = NAN; - } - constPredictionsOut[ std::slice(param.Tp, N_row, 1) ] = - const_predictions; - } - - std::valarray varianceOut( N_row + param.Tp ); - for ( size_t i = 0; i < param.Tp; i++ ) { - varianceOut[ i ] = NAN; - } - varianceOut[ std::slice(param.Tp, N_row, 1) ] = variance; - - //---------------------------------------------------- - // Create output DataFrame - //---------------------------------------------------- - size_t dataFrameColumms = param.const_predict ? 4 : 3; - - DataFrame dataFrame( N_row + param.Tp, dataFrameColumms ); - - if ( param.const_predict ) { - dataFrame.ColumnNames() = { "Observations", "Predictions", - "Pred_Variance", "Const_Predictions" }; - } - else { - dataFrame.ColumnNames() = {"Observations","Predictions","Pred_Variance"}; - } - - if ( N_time ) { - dataFrame.TimeName() = timeName; - dataFrame.Time() = timeOut; - } - - dataFrame.WriteColumn( 0, observations ); - dataFrame.WriteColumn( 1, predictionsOut ); - dataFrame.WriteColumn( 2, varianceOut ); - - if ( param.const_predict ) { - dataFrame.WriteColumn( 3, constPredictionsOut ); - } - -#ifdef DEBUG_ALL - std::cout << "FormatOutput() time " << timeOut.size() - << " pred " << predictionsOut.size() - << " obs " << observations.size() << std::endl; - std::cout << "FormatOutput() dataFrame -------------------" << std::endl; - std::cout << dataFrame; -#endif - - return dataFrame; -} - -//---------------------------------------------------------- -// Copy strings of time values into timeOut. -// If prediction times exceed times from the data, -// create new entries for the additional times. -//---------------------------------------------------------- -void FillTimes( Parameters param, - std::vector time, - std::vector &timeOut ) -{ - size_t N_time = time.size(); - size_t N_row = param.prediction.size(); - size_t max_pred_i = param.prediction[ N_row - 1 ]; - - if ( max_pred_i >= N_time ) { - // If tau > 0 end rows were deleted. max_pred_i might exceed time bounds - max_pred_i = N_time - 1; - } - - if ( timeOut.size() != N_row + param.Tp ) { - std::stringstream errMsg; - errMsg << "FillTimes(): timeOut vector length " << timeOut.size() - << " is not equal to the number of predictions + Tp " - << N_row + param.Tp << std::endl; - throw std::runtime_error( errMsg.str() ); - } - - // Fill in times guaranteed to be in param.prediction indices - for ( auto i = 0; i < N_row; i++ ) { - size_t pred_i = param.prediction[ i ]; - if ( pred_i < N_time ) { - timeOut[ i ] = time[ pred_i ]; - } - } - - // Now fill in times beyond param.prediction indices - if ( max_pred_i + param.Tp < N_time ) { - // All prediction times are available in time, get the rest - for ( auto i = 0; i < param.Tp; i++ ) { - timeOut[ N_row + i ] = time[ max_pred_i + i + 1 ]; - } - } - else { - // Tp introduces time values beyond the range of time - bool time_format_warning_printed = false; - - // Try to parse the last time vector string as a date or datetime - // if dtinfo.unrecognized_fmt = true; it is not a date or datetime - datetime_info dtinfo = parse_datetime( time[ max_pred_i ] ); - - for ( auto i = 0; i < param.Tp; i++ ) { - std::stringstream tss; - - if ( dtinfo.unrecognized_fmt ) { - // Numeric so add Tp - tss << std::stod( time[ max_pred_i ] ) + i + 1; - } - else { - int time_delta = i + 1; - // Get last two datetimes to compute time diff to add time delta - std::string time_new( time[ max_pred_i ] ); - std::string time_old( time[ max_pred_i - 1 ] ); - std::string new_time = - increment_datetime_str( time_old, time_new, time_delta ); - - // Add +ti if not a recognized format (datetime util returns "") - if ( new_time.size() ) { - tss << new_time; - } - else { - tss << time[ max_pred_i ] << " +" << i + 1; - - if ( not time_format_warning_printed ) { - std::cout << "FillTimes(): " - << "time column is unrecognized time format." - << "\n\tManually adding + tp to the last" - << " time column available." << std::endl; - time_format_warning_printed = true; - } - } - } - - timeOut[ N_row + i ] = tss.str(); - } - } -} - -//---------------------------------------------------------- -// Validate dataFrameIn rows against lib and pred indices -//---------------------------------------------------------- -void CheckDataRows( Parameters param, - DataFrame &dataFrameIn, - std::string call ) -{ - // param.prediction & library have been zero-offset in Validate() - // to convert from user specified data row to array indicies - size_t prediction_max_i = param.prediction[ param.prediction.size() - 1 ]; - size_t library_max_i = param.library [ param.library.size() - 1 ]; - - size_t shift; - if ( param.embedded ) { - shift = 0; - } - else { - if ( param.E < 1 ) { - std::stringstream errMsg; - errMsg << "CheckDataRows(): E = " << param.E << " is invalid.\n" ; - throw std::runtime_error( errMsg.str() ); - } - - shift = abs( param.tau ) * ( param.E - 1 ); - } - - if ( dataFrameIn.NRows() <= prediction_max_i ) { - std::stringstream errMsg; - errMsg << "CheckDataRows(): " << call - << ": The prediction index " - << prediction_max_i + 1 - << " exceeds the number of data rows " - << dataFrameIn.NRows(); - throw std::runtime_error( errMsg.str() ); - } - - if ( dataFrameIn.NRows() <= library_max_i + shift ) { - std::stringstream errMsg; - errMsg << "CheckDataRows(): " << call - << ": The library index " << library_max_i + 1 - << " + tau(E-1) " << shift << " = " - << library_max_i + 1 + shift - << " exceeds the number of data rows " - << dataFrameIn.NRows(); - throw std::runtime_error( errMsg.str() ); - } -} diff --git a/src/AuxFunc.h b/src/AuxFunc.h deleted file mode 100644 index 63c7ad6..0000000 --- a/src/AuxFunc.h +++ /dev/null @@ -1,52 +0,0 @@ -#ifndef AUXFUNC -#define AUXFUNC - -#include -#include // std::ref - -#include "Common.h" -#include "Neighbors.h" -#include "Embed.h" - -//---------------------------------------------------------------- -// Data Input, embedding and NN structure to accommodate -// common initial processing in Simplex and Smap -//---------------------------------------------------------------- -struct DataEmbedNN { - DataFrame *dataIn; - DataFrame dataFrame; - std::valarray targetVec; - Neighbors neighbors; - - // Constructors - DataEmbedNN() {} - - DataEmbedNN( DataFrame *dataIn, - DataFrame dataFrame, - std::valarray targetVec, - Neighbors neighbors ) : - dataIn( dataIn ), dataFrame( dataFrame ), targetVec( targetVec ), - neighbors( neighbors ) {} -}; - -// Prototypes -DataEmbedNN EmbedNN( DataFrame *dataIn, - Parameters ¶m, - bool checkDataRows = true ); - -DataFrame FormatOutput( Parameters param, - std::valarray predictions, - std::valarray const_predictions, - std::valarray variance, - std::valarray target_vec, - std::vector time, - std::string timeName ); - -void FillTimes( Parameters param, - std::vector time, - std::vector &timeOut ); - -void CheckDataRows( Parameters param, - DataFrame &dataFrameIn, - std::string call ); -#endif diff --git a/src/CCM.cc b/src/CCM.cc deleted file mode 100644 index 4d8edbd..0000000 --- a/src/CCM.cc +++ /dev/null @@ -1,808 +0,0 @@ - -#include -#include -#include -#include -#include - -#ifdef CCM_THREADED // Defined in makefile -// Two explicit CrossMap() threads are invoked. -// One for forward mapping, one for inverse mapping. -#include -#endif - -#include "Common.h" -#include "Embed.h" -#include "AuxFunc.h" - -namespace EDM_CCM { - std::mutex mtx; - std::mutex q_mtx; - std::queue< std::exception_ptr > exceptionQ; - // Define the initial maximum distance for neigbors to avoid sort() - // DBL_MAX is a Macro equivalent to: std::numeric_limits::max() - // The issue with std::sort is that it ignores ties... - double DistanceMax = std::numeric_limits::max(); - double DistanceLimit = std::numeric_limits::max() - 1; -} - -//---------------------------------------------------------------- -// forward declarations -//---------------------------------------------------------------- -void CrossMap( Parameters param, - DataFrame< double > dataFrameIn, - bool includeData, - const CrossMapValues &crossMapValues ); - -DataFrame< double > CCMDistances( const DataFrame< double > &dataBlock, - Parameters param ); - -Neighbors CCMNeighbors( const DataFrame< double > &Distances, - std::vector< size_t > lib_i, - Parameters param ); - -DataFrame SimplexProjection( Parameters param, - DataEmbedNN embedNN, - bool checkDataRows ); - -//---------------------------------------------------------------- -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: -// which is a wrapper for CrossMap() -//---------------------------------------------------------------- -CCMValues CCM( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - int E, - int Tp, - int knn, - int tau, - std::string columns, - std::string target, - std::string libSizes_str, - int sample, - bool random, - bool replacement, - unsigned seed, - bool includeData, - bool verbose ) -{ - //---------------------------------------------------------- - // Load data to dataFrameIn - //---------------------------------------------------------- - DataFrame< double > dataFrameIn( pathIn, dataFile ); - - CCMValues ccmValues = CCM( dataFrameIn, - pathOut, - predictFile, - E, - Tp, - knn, - tau, - columns, - target, - libSizes_str, - sample, - random, - replacement, - seed, - includeData, - verbose ); - return ccmValues; -} - -//---------------------------------------------------------------- -// API Overload 2: DataFrame passed in -// Implemented a wrapper for CrossMap() -//---------------------------------------------------------------- -CCMValues CCM( DataFrame< double > dataFrameIn, - std::string pathOut, - std::string predictFile, - int E, - int Tp, - int knn, - int tau, - std::string columns, - std::string target, - std::string libSizes_str, - int sample, - bool random, - bool replacement, - unsigned seed, - bool includeData, - bool verbose ) -{ - if ( not columns.size() ) { - throw std::runtime_error("CCM() must specify the column to embed."); - } - if ( not target.size() ) { - throw std::runtime_error("CCM() must specify the target."); - } - if ( not libSizes_str.size() ) { - throw std::runtime_error("CCM() must specify library sizes."); - } - - Parameters param = Parameters( Method::CCM, - "", // pathIn - "", // dataFile - pathOut, // - predictFile, // - "", // lib_str - "", // pred_str - E, // - Tp, // - knn, // - tau, // - 0, // theta - 0, // exclusionRadius - columns, // - target, // - false, // embedded - false, // const_predict - verbose, // - "", // SmapFile - "", // blockFile - "", // derivatives_str - 0, // svdSig - 0, // tikhonov - 0, // elasticNet - 0, // multi - libSizes_str, // - sample, // - random, // - replacement, // - seed ); // - - if ( param.columnNames.size() > 1 ) { - std::cout << "WARNING: CCM() Only the first column will be mapped.\n"; - } - - // Create Parameters object that switches column[0] and target - // for the inverse mapping - Parameters inverseParam( param ); // copy constructor - std::string newTarget( param.columns_str ); - inverseParam.columns_str = param.target_str; - inverseParam.target_str = newTarget; - - // Validate converts column_str & target_str to columnNames, targetName - inverseParam.Validate(); - -#ifdef DEBUG_ALL - std::cout << "CCM() params:\n"; - std::cout << param; - std::cout << "CCM() inverseParams:\n"; - std::cout << inverseParam; -#endif - - //------------------------------------------------------------ - // Setup DataFrames for output CrossMapValues structs - //------------------------------------------------------------ - DataFrame LibStats1( param.librarySizes.size(), 4, - "LibSize rho RMSE MAE" ); - DataFrame LibStats2( param.librarySizes.size(), 4, - "LibSize rho RMSE MAE" ); - size_t maxSamples; - if ( param.randomLib ) { - // Random samples from library - maxSamples = param.librarySizes.size() * param.subSamples; - } - else { - // Contiguous samples up to the size of the library - maxSamples = param.librarySizes.size(); - } - - DataFrame PredictionStats1( maxSamples, 8, - "N E nn tau LibSize rho RMSE MAE" ); - DataFrame PredictionStats2( maxSamples, 8, - "N E nn tau LibSize rho RMSE MAE" ); - - // Instantiate CrossMapValues output structs and insert DataFrames - CrossMapValues col_to_target = CrossMapValues(); - CrossMapValues target_to_col = CrossMapValues(); - - col_to_target.LibStats = LibStats1; - target_to_col.LibStats = LibStats2; - - if ( includeData ) { - col_to_target.PredictStats = PredictionStats1; - target_to_col.PredictStats = PredictionStats2; - } - -#ifdef CCM_THREADED - std::thread CrossMapColTarget( CrossMap, param, dataFrameIn, includeData, - std::ref( col_to_target ) ); - - std::thread CrossMapTargetCol( CrossMap, inverseParam, dataFrameIn, - includeData, std::ref( target_to_col ) ); - CrossMapColTarget.join(); - CrossMapTargetCol.join(); - - // If thread threw exception, get from queue and rethrow - if ( not EDM_CCM::exceptionQ.empty() ) { - std::lock_guard lck( EDM_CCM::q_mtx ); - - // Take the first exception in the queue - std::exception_ptr exceptionPtr = EDM_CCM::exceptionQ.front(); - - // Unroll all other exception from the thread/loops - while( not EDM_CCM::exceptionQ.empty() ) { - // JP When do these exception_ptr get deleted? Is it a leak? - EDM_CCM::exceptionQ.pop(); - } - std::rethrow_exception( exceptionPtr ); - } -#else - CrossMap( param, dataFrameIn, includeData, std::ref( col_to_target)); - CrossMap( inverseParam, dataFrameIn, includeData, std::ref( target_to_col)); -#endif - - //----------------------------------------------------------------- - // Output - //----------------------------------------------------------------- - // Create unified column names of output DataFrame - std::stringstream libRhoNames; - libRhoNames << "LibSize " - << param.columnNames[0] << ":" << param.targetName << " " - << param.targetName << ":" << param.columnNames[0]; - - // Unified LibStats output DataFrame - DataFrame PredictLibRho( param.librarySizes.size(), 3, - libRhoNames.str() ); - - PredictLibRho.WriteColumn( 0, col_to_target.LibStats.Column( 0 ) ); - PredictLibRho.WriteColumn( 1, col_to_target.LibStats.Column( 1 ) ); - PredictLibRho.WriteColumn( 2, target_to_col.LibStats.Column( 1 ) ); - - if ( param.predictOutputFile.size() ) { - // Write to disk - PredictLibRho.WriteData( param.pathOut, param.predictOutputFile ); - } - - // Output struct - CCMValues ccmValues; - ccmValues.AllLibStats = PredictLibRho; - - if ( includeData ) { - // Now handle data for each prediction instance - ccmValues.CrossMap1.PredictStats = col_to_target.PredictStats; - ccmValues.CrossMap2.PredictStats = target_to_col.PredictStats; - ccmValues.CrossMap1.Predictions = col_to_target.Predictions; - ccmValues.CrossMap2.Predictions = target_to_col.Predictions; - } - - return ccmValues; -} - -//---------------------------------------------------------------- -// CrossMap() -// Worker function for CCM. -// Return DataFrame of rho, RMSE, MAE values for param.librarySizes -// -// NOTE: This is a bit of a kludge at the moment... would be nice -// to pass in a reference to dataFrameIn, however... Embed() -// returns a dataBlock with (E-1)*tau fewer rows since -// partial data vectors are removed. To maintain proper -// alignment with the target in the dataFrameIn, -// dataFrameIn.DeletePartialDataRows() is called. If we use -// a reference to dataFrameIn, then the second thread will -// receive a reference to dataFrameIn that has rows deleted -// and Embed() will then not be correct. So at the moment -// we'll pass a copy of the dataFrameIn to each thread. -//---------------------------------------------------------------- -void CrossMap( Parameters paramCCM, - DataFrame< double > dataFrameIn, - bool includeData, - const CrossMapValues &crossMapValuesIn ) { - - // Get a local reference for CrossMapValues - CrossMapValues &crossMapValues = - const_cast< CrossMapValues & >( crossMapValuesIn ); - - if ( paramCCM.verbose ) { - std::lock_guard lck( EDM_CCM::mtx ); - std::stringstream msg; - msg << "CrossMap(): Simplex cross mapping from " - << paramCCM.columnNames[0] - << " to " << paramCCM.targetName << " E=" << paramCCM.E - << " knn=" << paramCCM.knn << " Library range: [" - << paramCCM.libSizes_str << "] "; - for ( size_t i = 0; i < paramCCM.librarySizes.size(); i++ ) { - msg << paramCCM.librarySizes[ i ] << " "; - } msg << std::endl << std::endl; - std::cout << msg.str(); - } - - try { - //------------------------------------------------------------ - // Generate embedding on data to be cross mapped (-c column) - // dataBlock will have tau * (E-1) fewer rows than dataFrameIn - // JP: Should this be allocated on the heap? - //------------------------------------------------------------ - DataFrame dataBlock = Embed( dataFrameIn, - paramCCM.E, - paramCCM.tau, - paramCCM.columnNames[0], - paramCCM.verbose ); - - size_t N_row = dataBlock.NRows(); - - // NOTE: No need to adjust param.library and param.prediction indices - // [param.DeleteLibPred( shift );] since pred will be created - // below based on N_row of dataBlock. - - //-------------------------------------------------------------- - // Remove dataFrameIn rows to match embedded dataBlock with - // partial data rows ignored: CrossMap() -> Embed() -> MakeBlock() - // This removal of partial data rows is also done in EmbedNN() - //-------------------------------------------------------------- - if ( paramCCM.E < 1 ) { - std::stringstream errMsg; - errMsg << "CrossMap(): E = " << paramCCM.E << " is invalid.\n" ; - throw std::runtime_error( errMsg.str() ); - } - - size_t shift = abs( paramCCM.tau ) * ( paramCCM.E - 1 ); - { - std::lock_guard lck( EDM_CCM::mtx ); - if ( not dataFrameIn.PartialDataRowsDeleted() ) { - // Not thread safe. - dataFrameIn.DeletePartialDataRows( shift, paramCCM.tau ); - } - } - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << ">>>> CrossMap() dataFrameIn-----------------------\n"; - std::cout << dataFrameIn; - std::cout << "<<<< dataFrameIn----------------------------------\n"; - std::cout << ">>>> dataBlock------------------------------------\n"; - std::cout << dataBlock; - std::cout << "<<<< dataBlock------------------------------------\n"; - } -#endif - - //----------------------------------------------------------------- - // Create Parameters for SimplexProjection - // Add library and prediction indices for the entire library - //----------------------------------------------------------------- - std::stringstream ss; - ss << "1 " << N_row; - paramCCM.lib_str = ss.str(); - paramCCM.pred_str = ss.str(); - // Validate converts lib_str, pred_str to library & prediction vectors - paramCCM.Validate(); - - //----------------------------------------------------------------- - // Set number of samples - //----------------------------------------------------------------- - size_t maxSamples; - if ( paramCCM.randomLib ) { - // Random samples from library - maxSamples = paramCCM.subSamples; - } - else { - // Contiguous samples up to the size of the library - maxSamples = 1; - } - - //----------------------------------------------------------------- - // Create random number generator: DefaultRandEngine - //----------------------------------------------------------------- - if ( paramCCM.randomLib ) { - if ( paramCCM.seed == 0 ) { - // Select a random seed - typedef std::chrono::high_resolution_clock CCMclock; - CCMclock::time_point beginning = CCMclock::now(); - CCMclock::duration duration = CCMclock::now() - beginning; - paramCCM.seed = duration.count(); - } - } - std::default_random_engine DefaultRandomEngine( paramCCM.seed ); - - //----------------------------------------------------------------- - // Distance for all possible pred : lib E-dimensional vector pairs - // Distances is a square Matrix of all row to to row distances - //----------------------------------------------------------------- - DataFrame< double > Distances = CCMDistances( std::ref( dataBlock ), - paramCCM ); - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "CrossMap() " << paramCCM.columnNames[0] << " to " - << paramCCM.targetName << " Distances: " << Distances.NRows() - << " x " << Distances.NColumns() << std::endl; - } -#endif - - //---------------------------------------------------------- - // Predictions - //---------------------------------------------------------- - size_t predictionCount = 0; - // Loop for library sizes - for ( size_t lib_size_i = 0; - lib_size_i < paramCCM.librarySizes.size(); lib_size_i++ ) { - - size_t lib_size = paramCCM.librarySizes[ lib_size_i ]; - - // Create random RNG sampler for this lib_size - std::uniform_int_distribution distribution( 0, N_row - 1 ); - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "lib_size: " << lib_size - << " ------------------------------------------\n"; - } -#endif - - std::valarray< double > rho ( maxSamples ); - std::valarray< double > RMSE( maxSamples ); - std::valarray< double > MAE ( maxSamples ); - - // Loop for subsamples - for ( size_t n = 0; n < maxSamples; n++ ) { - - // Vector of row indices to include in this lib_size evaluation - std::vector< size_t > lib_i( lib_size ); - - if ( paramCCM.randomLib ) { - // Uniform random sample of rows - - if ( paramCCM.replacement ) { - // With replacement - for ( size_t i = 0; i < lib_size; i++ ) { - lib_i[ i ] = distribution( DefaultRandomEngine ); - } - } - else { - // Without replacement lib_size elements from [0, N_row-1] - // Robert W. Floyd's algorithm - // NOTE: c++17 has the sample() function in - if ( lib_size >= N_row ) { - std::stringstream errMsg; - errMsg << "CrossMap(): lib_size=" << lib_size - << " must be less than N_row=" << N_row - << " for random sample without replacement."; - throw std::runtime_error( errMsg.str() ); - } - - // unordered set to store samples - std::unordered_set samples; - - // Sample and insert values into samples - size_t r = 0; - while( samples.size() < lib_size ) { - size_t v = distribution( DefaultRandomEngine ); - if ( not samples.insert( v ).second ) { - samples.insert( r ); - r++; - } - } - - // Copy samples into result - std::vector result(samples.begin(), samples.end()); - - // Shuffle result - std::shuffle( result.begin(), result.end(), - DefaultRandomEngine ); - - lib_i = result; // Copy result to lib_i - } - } - else { - // Not random samples, contiguous samples increasing size - if ( lib_size >= N_row ) { - // library size exceeded, back down - lib_i.resize( N_row ); - std::iota( lib_i.begin(), lib_i.end(), 0 ); - lib_size = N_row; - - if ( paramCCM.verbose ) { - std::stringstream msg; - msg << "CCM(): Sequential library samples," - << " max lib_size is " << N_row - << ", lib_size has been limited.\n"; - std::cout << msg.str(); - } - } - else { - // Contiguous blocks up to N_rows = maxSamples - if ( n + lib_size < N_row ) { - std::iota( lib_i.begin(), lib_i.end(), n ); - } - else { - // n + lib_size > N_row, wrap around to data origin - std::vector< size_t > lib_start( N_row - n ); - std::iota( lib_start.begin(), lib_start.end(), n ); - - size_t max_i = std::min( lib_size-(N_row - n), N_row ); - std::vector< size_t > lib_wrap( max_i ); - std::iota( lib_wrap.begin(), lib_wrap.end(), 0 ); - - // Build new lib_i - lib_i = std::vector< size_t > ( lib_start ); - lib_i.insert( lib_i.end(), - lib_wrap.begin(), - lib_wrap.end() ); - lib_size = lib_i.size(); - } - } - } - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "lib_i: (" << lib_i.size() << ") "; - for ( size_t i = 0; i < lib_i.size(); i++ ) { - std::cout << lib_i[i] << " "; - } std::cout << std::endl; - } -#endif - - //---------------------------------------------------------- - // Nearest neighbors : Local CCMNeighbors() function - //---------------------------------------------------------- - Neighbors neighbors = CCMNeighbors( Distances, lib_i, paramCCM ); - - //---------------------------------------------------------- - // Subset dataFrameIn to lib_i - //---------------------------------------------------------- - DataFrame< double > dataFrameLib_i( lib_i.size(), - dataFrameIn.NColumns(), - dataFrameIn.ColumnNames() ); - - for ( size_t i = 0; i < lib_i.size(); i++ ) { - dataFrameLib_i.WriteRow( i, dataFrameIn.Row( lib_i[ i ] ) ) ; - } - - std::valarray targetVec = - dataFrameLib_i.VectorColumnName( paramCCM.targetName ); - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "dataFrameLib_i -------------------------------\n"; - std::cout << dataFrameLib_i; - } -#endif - - //---------------------------------------------------------- - // Pack embedding, target, neighbors for SimplexProjection - //---------------------------------------------------------- - DataEmbedNN embedNN = DataEmbedNN( &dataFrameLib_i, dataBlock, - targetVec, neighbors ); - - //---------------------------------------------------------- - // Simplex Projection: lib_str & pred_str set from N_row - //---------------------------------------------------------- - DataFrame S = SimplexProjection( paramCCM, embedNN, false ); - - VectorError ve = ComputeError( - S.VectorColumnName( "Observations" ), - S.VectorColumnName( "Predictions" ) ); - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "CCM Simplex ---------------------------------\n"; - S.MaxRowPrint() = S.NRows(); - std::cout << S; - std::cout << "rho " << ve.rho << " RMSE " << ve.RMSE - << " MAE " << ve.MAE << std::endl; - } -#endif - // Record values for these samples - rho [ n ] = ve.rho; - RMSE[ n ] = ve.RMSE; - MAE [ n ] = ve.MAE; - - if ( includeData ) { - // Save stats for this prediction - std::valarray< double > predOutVec( 8 ); - predOutVec[ 0 ] = predictionCount + 1; // N - predOutVec[ 1 ] = paramCCM.E; // E - predOutVec[ 2 ] = paramCCM.knn; // nn - predOutVec[ 3 ] = paramCCM.tau; // tau - predOutVec[ 4 ] = lib_size; // LibSize - predOutVec[ 5 ] = ve.rho; // rho - predOutVec[ 6 ] = ve.RMSE; // RMSE - predOutVec[ 7 ] = ve.MAE; // MAE - - crossMapValues.PredictStats.WriteRow(predictionCount,predOutVec); - - // Save predictions - crossMapValues.Predictions.push_front( S ); - } - - predictionCount++; - } // for ( n = 0; n < maxSamples; n++ ) - - std::valarray< double > statVec( 4 ); - statVec[ 0 ] = lib_size; - statVec[ 1 ] = rho.sum() / maxSamples; - statVec[ 2 ] = RMSE.sum() / maxSamples; - statVec[ 3 ] = MAE.sum() / maxSamples; - - crossMapValues.LibStats.WriteRow( lib_size_i, statVec ); - } // for ( lib_size : param.librarySizes ) - - } // try - catch(...) { - // push exception pointer onto queue for main thread to catch - std::lock_guard lck( EDM_CCM::q_mtx ); - EDM_CCM::exceptionQ.push( std::current_exception() ); - } -} - -//--------------------------------------------------------------------- -// Note that for CCM the library and prediction rows are the same. -// Note that dataBlock does NOT have the time in column 0. -// -// Return Distances: a square matrix with distances. -// Matrix elements D[i,j] hold the distance between the E-dimensional -// phase space point (vector) between rows (observations) i and j. -//--------------------------------------------------------------------- -DataFrame< double > CCMDistances( const DataFrame< double > &dataBlock, - Parameters param ) { - - size_t N_row = dataBlock.NRows(); - - size_t E = param.E; - - DataFrame< double > D = DataFrame< double >( N_row, N_row ); - - // Initialise D to DistanceMax to avoid sort() : Add init constructor? - std::valarray< double > row_init( EDM_CCM::DistanceMax, N_row ); - for ( size_t row = 0; row < N_row; row++ ) { - D.WriteRow( row, row_init ); - } - - for ( size_t row = 0; row < N_row; row++ ) { - // Get E-dimensional vector from this library row - std::valarray< double > v1_ = dataBlock.Row( row ); - // The first column (i=0) is NOT time, use it - std::valarray< double > v1 = v1_[ std::slice( 0, E, 1 ) ]; - - // Only compute upper triangular D, the diagonal and - // lower left are redundant: (col < N_row); row >= col - for ( size_t col = 0; col < N_row; col++ ) { - // Avoid redundant computations - if ( row >= col ) { - continue; // Computed in upper triangle, copied below - } - - // Find distance between vector (v) and other library vector - std::valarray< double > v2_ = dataBlock.Row( col ); - // The first column (i=0) is NOT time, use it - std::valarray< double > v2 = v2_[ std::slice( 0, E, 1 ) ]; - - D( row, col ) = Distance( v1, v2, DistanceMetric::Euclidean ); - - // Insert degenerate values since D[i,j] = D[j,i] - D( col, row ) = D( row, col ); - } - } - return D; -} - -//--------------------------------------------------------------------- -// Return Neighbors { neighbors, distances }. neighbors is a matrix of -// row indices in the library matrix. Each neighbors row represents one -// prediction vector. Columns are the indices of knn nearest neighbors -// for the prediction vector (phase-space point) in the library matrix. -// distances is a matrix with the same shape as neighbors holding the -// corresponding distance values in each row. -// -// Note that the indices in neighbors are not the original indices in -// the libraryMatrix rows (observations), but are with respect to the -// distances subset defined by the list of rows lib_i, and so have values -// from 0 to len(lib_i)-1. -// -//--------------------------------------------------------------------- -Neighbors CCMNeighbors( const DataFrame< double > &DistancesIn, - std::vector< size_t > lib_i, - Parameters param ) { - - size_t N_row = lib_i.size(); - size_t knn = param.knn; - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "CCMNeighbors Distances\n"; - for ( size_t r = 0; r < 5; r++ ) { - for ( int c = 0; c < 5; c++ ) { - std::cout << DistancesIn(r,c) << " "; - } std::cout << std::endl; - } - std::cout << "lib_i N_row: " << N_row - << " DistancesIn NRow: " << DistancesIn.NRows() << std::endl; - } -#endif - - // Matrix to hold libraryMatrix row indices - // One row for each prediction vector, knn columns for each index - DataFrame< size_t > neighbors( N_row, knn ); - - // Matrix to hold libraryMatrix knn distance values - // One row for each prediction vector, k_NN columns for each index - DataFrame< double > distances( N_row, knn ); - - // For each prediction vector (row in predictionMatrix) find the list - // of library indices that are the closest knn points - size_t row = 0; - std::valarray< double > knn_distances( knn ); - std::valarray< size_t > knn_neighbors( knn ); - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "CCMNeighbors lib_i: "; - for ( size_t i = 0; i < lib_i.size(); i++ ) { - std::cout << lib_i[i] << " "; - } std::cout << std::endl << std::flush; - } -#endif - - size_t shift0 = abs( param.tau ) * param.E; - - for ( auto row_i : lib_i ) { - - // Take Distances( row, col ) a row at a time - // col represent the other row distance - std::valarray< double > dist_row = DistancesIn.Row( row_i ); - - // These new column indices are with respect to the lib_i vector - // not the original Distances with all other columns - - // Reset the neighbor and distance vectors for this pred row - for ( size_t i = 0; i < knn; i++ ) { - knn_neighbors[ i ] = 0; - // This avoids the need to sort the distances of this row - knn_distances[ i ] = EDM_CCM::DistanceMax; - } - - for ( size_t col_i = 0; col_i < N_row; col_i++ ) { - - if ( col_i > N_row - shift0 ) { - continue; - } - - double d_i = dist_row[ lib_i[col_i] ]; - // If d_i is less than values in knn_distances, add to list - auto max_it = std::max_element( begin( knn_distances ), - end ( knn_distances ) ); - if ( d_i < *max_it ) { - size_t max_i = std::distance( begin(knn_distances), max_it ); - knn_neighbors[ max_i ] = col_i; // Save the index - knn_distances[ max_i ] = d_i; // Save the value - } - } - - neighbors.WriteRow( row, knn_neighbors ); - distances.WriteRow( row, knn_distances ); - - row = row + 1; - } - - Neighbors ccmNeighbors = Neighbors(); - ccmNeighbors.neighbors = neighbors; - ccmNeighbors.distances = distances; - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_CCM::mtx ); - std::cout << "CCMNeighbors knn_neighbors\n"; - for ( size_t r = 0; r < 5; r++ ) { - for ( int c = 0; c < ccmNeighbors.neighbors.NColumns(); c++ ) { - std::cout << ccmNeighbors.neighbors(r,c) << " "; - } std::cout << std::endl; - } - } -#endif - - return ccmNeighbors; -} diff --git a/src/Common.cc b/src/Common.cc deleted file mode 100644 index ebdd311..0000000 --- a/src/Common.cc +++ /dev/null @@ -1,192 +0,0 @@ - -#include -#include - -#include "Common.h" - -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -std::string ToLower( std::string str ) { - - std::string lowerStr( str ); - std::transform( lowerStr.begin(), lowerStr.end(), - lowerStr.begin(), ::tolower ); - - return lowerStr; -} - -//------------------------------------------------------------------- -// Called in two different contexts: -// DataFrame.h ReadData() Are data columns numeric or labels? -// Parameter.cc Validate() Are columns/target integer index or label? -// -// param integer : true: only digits, false: digits plus '-', '.' -// return -// true : str has only numeric characters -// false : str has non-numeric characters -//------------------------------------------------------------------- -bool OnlyDigits( std::string str, bool integer ) { - - if ( not str.size() ) { - throw std::runtime_error( "OnlyDigits(): String is empty.\n" ); - } - - // Remove whitespace - std::string str_( str ); - str_.erase( std::remove_if( str_.begin(), str_.end(), ::isspace ), - str_.end() ); - - std::string digits; - if ( integer ) { - digits = "0123456789"; - } - else { - digits = "-.0123456789"; - } - - // Is str_ purely numeric characters? - bool onlyDigits = strspn(str_.c_str(), digits.c_str()) == str_.size(); - - return onlyDigits; -} - -//---------------------------------------------------------------- -// SplitString -// -// Purpose: like Python string.split() -// -// Arguments: inString : string to be split -// delimeters : string of delimeters -// -// Note: A typical delimeter string: delimeters = " \t,\n;" -// -// Return: vector of tokens -//---------------------------------------------------------------- -std::vector SplitString( std::string inString, - std::string delimeters ) { - - size_t pos = 0; - size_t eos = 0; - size_t wordStart = 0; - size_t wordEnd = 0; - - bool foundStart = false; - bool foundEnd = false; - - std::vector splitString; - - std::string word; - - eos = inString.length(); - - while ( pos <= eos ) { - if ( not foundStart ) { - if ( delimeters.find( inString[pos] ) == delimeters.npos ) { - // this char (inString[pos]) is not a delimeter - wordStart = pos; - foundStart = true; - pos++; - continue; - } - } - if ( foundStart and not foundEnd ) { - if ( delimeters.find( inString[pos] ) != delimeters.npos - or pos == eos ) { - // this char (inString[pos]) is a delimeter or - // at the end of the string - wordEnd = pos; - foundEnd = true; - } - } - if ( foundStart and foundEnd ) { - foundStart = false; - foundEnd = false; - - word = inString.substr( wordStart, wordEnd - wordStart ); - - // remove whitespace - word.erase( std::remove_if( word.begin(), word.end(), ::isspace ), - word.end() ); - - splitString.push_back( word ); - } - if ( pos == eos ) { - break; - } - pos++; - } - - return splitString; -} - -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -VectorError ComputeError( std::valarray< double > obsIn, - std::valarray< double > predIn ) { - - // Check for nan in vectors - size_t nanObs = 0; - size_t nanPred = 0; - for ( auto o : obsIn ) { if ( std::isnan( o ) ) { nanObs++; } } - for ( auto p : predIn ) { if ( std::isnan( p ) ) { nanPred++; } } - - if ( nanObs != nanPred ) { - std::stringstream errMsg; - errMsg << "ComputeError(): obs has " << nanObs << " nan, pred has " - << nanPred << " nan. ComputeError result invalid.\n"; - std::cout << errMsg.str() << std::flush; - } - - // JP: Assume that nan are at the beginning of predictions and - // at the end of observations... probably not a robust assumption - // but this is the case if the data were prepared from Embed() - // and there were initially no nans - std::valarray pred( predIn.size() - 2 * nanPred ); - std::valarray obs ( predIn.size() - 2 * nanPred ); - if ( nanPred > 0 ) { - // ignore nanPred initial pred, and nanObs end obs - pred = predIn[ std::slice( nanPred, pred.size(), 1 ) ]; - obs = obsIn [ std::slice( nanPred, pred.size(), 1 ) ]; - } - else { - pred = std::valarray( predIn ); - obs = std::valarray( obsIn ); - } - - size_t N = pred.size(); - - std::valarray< double > two( 2, N ); // Vector of 2's for squaring - - double sumPred = pred.sum(); - double sumObs = obs.sum(); - double meanPred = sumPred / N; - double meanObs = sumObs / N; - double sumSqrPred = pow( pred, two ).sum(); - double sumSqrObs = pow( obs, two ).sum(); - double sumErr = abs( obs - pred ).sum(); - double sumSqrErr = pow( obs - pred, two ).sum(); - double sumProd = ( obs * pred ).sum(); - - double rho; // Pearson correlation coefficient - - if ( sumSqrPred * N == sumSqrPred ) { - rho = 0; - } - else { - rho = ( sumProd - N * meanObs * meanPred ) / - ( std::sqrt( ( sumSqrObs - N * pow( meanObs, 2 ) ) ) * - std::sqrt( ( sumSqrPred - N * pow( meanPred, 2 ) ) ) ); - } - - VectorError vectorError = VectorError(); - - vectorError.RMSE = sqrt( sumSqrErr / N ); - - vectorError.MAE = sumErr / N; - - vectorError.rho = rho; - - return vectorError; -} diff --git a/src/Common.h b/src/Common.h index 3d8edd6..4f9f3d6 100644 --- a/src/Common.h +++ b/src/Common.h @@ -1,3 +1,5 @@ +// File for declarations common to all EDM code + #ifndef COMMON_H #define COMMON_H @@ -11,375 +13,15 @@ #include #ifdef _MSC_VER -#include // macro constants for MSVC C++ operators not in ISO646 -#endif - -#include "DataFrame.h" // has #include Common.h - -// forward declaration for SMap solver -std::valarray < double > SVD( DataFrame < double > A, - std::valarray< double > B ); - -// Enumerations -enum class Method { None, Embed, Simplex, SMap, CCM }; -enum class DistanceMetric { Euclidean, Manhattan }; - -//--------------------------------------------------------- -// Data structs -//--------------------------------------------------------- -struct VectorError { - double rho; - double RMSE; - double MAE; -}; - -struct SMapValues { - DataFrame< double > predictions; - DataFrame< double > coefficients; -}; - -// Return object for CrossMap() worker function -struct CrossMapValues { - DataFrame< double > LibStats; // mean libsize, rho, RMSE, MAE - DataFrame< double > PredictStats; // each predict libsize, rho, RMSE, MAE - std::forward_list< DataFrame< double > > Predictions; -}; - -// Return object for CCM() with two CrossMapValues -struct CCMValues { - DataFrame< double > AllLibStats; // unified mean libsize, rho, RMSE, MAE - CrossMapValues CrossMap1; - CrossMapValues CrossMap2; -}; - -struct MultiviewValues { - DataFrame< double > Combo_rho; // col_i..., rho, MAE, RMSE - DataFrame< double > Predictions; - std::vector< std::string > Combo_rho_table; // includes column names - -#ifdef MULTIVIEW_VALUES_OVERLOAD - // Don't define constructors for the setuptools module build on Windows - // The MSVC compiler with pybind11 does not handle overloads easily... - // https://pybind11.readthedocs.io/en/stable/classes.html - - // Constructors - MultiviewValues(); - - MultiviewValues( DataFrame< double > combo_rho, - DataFrame< double > predictions, - std::vector< std::string > combo_rho_table ): - Combo_rho( combo_rho ), Predictions( predictions ), - Combo_rho_table( combo_rho_table ) {} +#include // Macro constants for MSVC C++ operators not in ISO646 #endif -}; - -//------------------------------------------------------------- -// Prototypes -// Primary API functions generally have two call-signatures. -// The first takes a (path, file name) pair specifying the data -// file image on disk to be loaded and converted to a data frame. -// The second replaces these two arguments with a DataFrame object. -// -// NOTE: These are the first declarations seen by the compiler -// for the API and provide default argument values -//------------------------------------------------------------- -std::string ToLower ( std::string str ); -bool OnlyDigits( std::string str, bool integerOnly ); - -std::vector SplitString( std::string inString, - std::string delimeters ); - -VectorError ComputeError( std::valarray< double > obs, - std::valarray< double > pred ); - -std::string increment_datetime_str( std::string datetime1, - std::string datetime2, - int tp ); - -// API functions Embed() and MakeBlock() are in Embed.h Embed.cc - -DataFrame Simplex( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - int exclusionRadius = 0, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool const_predict = false, - bool verbose = true ); - -DataFrame Simplex( DataFrame< double > &dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - int exclusionRadius = 0, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool const_predict = false, - bool verbose = true ); - -// SMap is a special case since it can be called with a function pointer -// to the SVD solver. This is done so that interfaces such as pybind11 -// can provide their own object for the solver. -// 1) Data path/file with default SVD (LAPACK) assigned in Smap.cc 2) -SMapValues SMap( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - double theta = 0, - int exclusionRadius = 0, - std::string columns = "", - std::string target = "", - std::string smapFile = "", - std::string derivatives = "", - bool embedded = false, - bool const_predict = false, - bool verbose = true ); - -// 2) DataFrame with default SVD (LAPACK) assigned in Smap.cc 2) -SMapValues SMap( DataFrame< double > &dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - double theta = 0, - int exclusionRadius = 0, - std::string columns = "", - std::string target = "", - std::string smapFile = "", - std::string derivatives = "", - bool embedded = false, - bool const_predict = false, - bool verbose = true ); - -// 3) Data path/file with external solver object, init to default SVD -SMapValues SMap( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - double theta = 0, - int exclusionRadius = 0, - std::string columns = "", - std::string target = "", - std::string smapFile = "", - std::string derivatives = "", - std::valarray (*solver)(DataFrame < double >, - std::valarray < double >) = &SVD, - bool embedded = false, - bool const_predict = false, - bool verbose = true ); - -// 4) DataFrame with external solver object, init to default SVD -SMapValues SMap( DataFrame< double > &dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - double theta = 0, - int exclusionRadius = 0, - std::string columns = "", - std::string target = "", - std::string smapFile = "", - std::string derivatives = "", - std::valarray (*solver)(DataFrame < double >, - std::valarray < double >) = &SVD, - bool embedded = false, - bool const_predict = false, - bool verbose = true ); - -CCMValues CCM( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - int E = 0, - int Tp = 0, - int knn = 0, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - std::string libSizes_str = "", - int sample = 0, - bool random = true, - bool replacement = false, - unsigned seed = 0, // seed=0: use RNG - bool includeData = false, - bool verbose = true ); - -CCMValues CCM( DataFrame< double > dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - int E = 0, - int Tp = 0, - int knn = 0, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - std::string libSizes_str = "", - int sample = 0, - bool random = true, - bool replacement = false, - unsigned seed = 0, // seed=0: use RNG - bool includeData = false, - bool verbose = true ); - -MultiviewValues Multiview( std::string pathIn = "./", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int D = 0, - int E = 1, - int Tp = 1, - int knn = 0, - int tau = -1, - std::string columns = "", - std::string target = "", - int multiview = 0, - int exclusionRadius = 0, - bool trainLib = true, - bool verbose = false, - unsigned nThreads = 4 ); - -MultiviewValues Multiview( DataFrame< double >, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int D = 0, - int E = 1, - int Tp = 1, - int knn = 0, - int tau = -1, - std::string columns = "", - std::string target = "", - int multiview = 0, - int exclusionRadius = 0, - bool trainLib = true, - bool verbose = false, - unsigned nThreads = 4 ); - -DataFrame EmbedDimension( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int maxE = 10, - int Tp = 1, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool verbose = true, - unsigned nThreads = 4 ); -DataFrame EmbedDimension( DataFrame< double > &dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int maxE = 10, - int Tp = 1, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool verbose = true, - unsigned nThreads = 4 ); +#include "DataFrame.h" // Has #include Common.h -DataFrame PredictInterval( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int maxTp = 10, - int E = 0, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool verbose = true, - unsigned nThreads = 4 ); +// Define EDM function and class signatures in its own namespace -DataFrame PredictInterval( DataFrame< double > &dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - int maxTp = 10, - int E = 0, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool verbose = true, - unsigned nThreads = 4 ); +namespace cppEDM { -DataFrame PredictNonlinear( std::string pathIn = "./data/", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - std::string theta = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool verbose = true, - unsigned nThreads = 4 ); +} -DataFrame PredictNonlinear( DataFrame< double > &dataFrameIn, - std::string pathOut = "./", - std::string predictFile = "", - std::string lib = "", - std::string pred = "", - std::string theta = "", - int E = 0, - int Tp = 1, - int knn = 0, - int tau = -1, - std::string colNames = "", - std::string targetName = "", - bool embedded = false, - bool verbose = true, - unsigned nThreads = 4 ); #endif diff --git a/src/DateTime.h b/src/DateTime.h deleted file mode 100644 index 9084864..0000000 --- a/src/DateTime.h +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef DATETIMEUTIL_H -#define DATETIMEUTIL_H - -#include // for testing -#include // time formatting -#include -#include -#include -#include -#include - -const int iso_start_year = 1900; -const int iso_start_month = 1; - -struct datetime_info { - struct tm time = {}; - std::string datetime_fmt; - bool unrecognized_fmt = false; -}; - -// Prototypes -void parse_datetime_str ( struct tm & time_obj, - std::string datetime_str, - bool date_fmt ); - -datetime_info parse_datetime ( std::string datetime ); - -std::string increment_datetime_str ( std::string datetime1, - std::string datetime2, - int tp ); - -#endif diff --git a/src/DateTimeUtil.cc b/src/DateTimeUtil.cc deleted file mode 100644 index 6dfc69b..0000000 --- a/src/DateTimeUtil.cc +++ /dev/null @@ -1,172 +0,0 @@ -#include "DateTime.h" - -// Provide some utility for parsing datetime std::strings -// to add some tp increment to a datetime std::string past the given -// range the time column - -//--------------------------------------------------------------------- -// TIME FORMATS supported: -// YYYY-MM-DD -// HH:MM:SS -// YYYY-MM-DDTHH:MM:SS (2019-06-30T10:26:10) -// hh:mm:ss.sss -//--------------------------------------------------------------------- -// regex's used in parsing and their time formats. in pair for easier checking -std::regex regEx_yyyymmdd ("\\d{4}-\\d{2}-\\d{2}"); -std::string fmt_yyyymmdd ("%Y-%m-%d"); -std::regex regEx_hhmmss ("\\d{2}:\\d{2}:\\d{2}"); -std::string fmt_hhmmss ("%H:%M:%S"); -std::regex regEx_yymmddthhmmss ("\\d{4}-\\d{2}-\\d{2}T\\d{2}:\\d{2}:\\d{2}"); -std::string fmt_yymmddthhmmss ("%Y-%m-%dT%H:%M:%S"); -std::regex regEx_yymmddhhmmss ("\\d{4}-\\d{2}-\\d{2} \\d{2}:\\d{2}:\\d{2}"); -std::string fmt_yymmddhhmmss ("%Y-%m-%d %H:%M:%S"); -std::regex regEx_hhmmsssss ("\\d{2}:\\d{2}:\\d{2}\\.\\d{3}"); -std::string fmt_hhmmsssss ("%H:%M:%S"); - -//---------------------------------------------------------------------- -// Parse a date or time string into a tm obj -// @param tm : the tm object to populate -// @param datetime_str : the date or time string to populate -// @param date_fmt : true if this is a date object -// @return : none, just populates the tm obj -//---------------------------------------------------------------------- -void parse_datetime_str ( struct tm & time_obj, - std::string datetime_str, - bool date_fmt ) { - // parsing delim is different for date or time - char parse_delim = date_fmt ? '-' : ':'; - - // parse the string for it's tokens - std::stringstream parseable_str ( datetime_str ); - std::string token; - std::vector tokens; - - while( getline( parseable_str, token, parse_delim ) ) { - tokens.push_back( token ); - } - - // populate the date time obj - if ( date_fmt ) { - time_obj.tm_mday = stod(tokens[2]); - time_obj.tm_mon = stod(tokens[1]) - iso_start_month; - time_obj.tm_year = stod(tokens[0]) - iso_start_year; - } - else { - time_obj.tm_sec = stod(tokens[2]); - time_obj.tm_min = stod(tokens[1]); - time_obj.tm_hour = stod(tokens[0]); - } - - int err = mktime( &time_obj ); - - if ( err < 0 ) { - std::stringstream errMsg; - errMsg << "parse_datetime_str() mktime failed on " - << datetime_str; - throw( errMsg.str() ); - } -} - -//---------------------------------------------------------------------- -// Parse the datetime std::string into a struct tm -// @param datetime : the datetime to parse -// @return datetime : the datetime to parse -//---------------------------------------------------------------------- -datetime_info parse_datetime ( std::string datetime ) { - - datetime_info output; - - // check which time format we have. populate output with each case - if ( std::regex_match( datetime, regEx_yyyymmdd ) ) { - output.datetime_fmt = fmt_yyyymmdd; - parse_datetime_str( output.time, datetime, true ); - } - else if ( std::regex_match( datetime, regEx_hhmmss )) { - output.datetime_fmt = fmt_hhmmss; - parse_datetime_str( output.time, datetime, false ); - } - else if ( std::regex_match( datetime, regEx_yymmddhhmmss )) { - output.datetime_fmt = fmt_yymmddhhmmss; - // split by " ", then split first by - second by : - int delim_pos = datetime.find(' '); - std::string date = datetime.substr(0, delim_pos); - std::string time = datetime.substr(delim_pos+1, datetime.size()); - parse_datetime_str( output.time, date, true ); - parse_datetime_str( output.time, time, false ); - } - else if ( std::regex_match( datetime, regEx_yymmddthhmmss )) { - output.datetime_fmt = fmt_yymmddthhmmss; - // split by T, then split first by - second by : - int delim_pos = datetime.find('T'); - std::string date = datetime.substr(0, delim_pos); - std::string time = datetime.substr(delim_pos+1, datetime.size()); - parse_datetime_str( output.time, date, true ); - parse_datetime_str( output.time, time, false ); - } - else if ( std::regex_match( datetime, regEx_hhmmsssss )) { - output.datetime_fmt = fmt_hhmmsssss; - // trim the milliseconds off and parse - datetime = datetime.substr(0,datetime.size()-4); - parse_datetime_str( output.time, datetime, false ); - } - else { - output.unrecognized_fmt = true; - } - return output; -} - -//---------------------------------------------------------------------- -// Generate a new datetime + delta past the range of given -//---------------------------------------------------------------------- -// -// @params datetime1/2 : the two last time std::strings -// to compute the delta unit -// we increment from datetime2 -// @param tp : the amount to increment the time diff by -// @return : the new incremented timestd::string -//---------------------------------------------------------------------- -std::string increment_datetime_str ( std::string datetime1, - std::string datetime2, int tp ) { - //parse datetimes - datetime_info dtinfo1 = parse_datetime( datetime1 ); - datetime_info dtinfo2 = parse_datetime( datetime2 ); - - if ( dtinfo1.unrecognized_fmt or dtinfo2.unrecognized_fmt ) { - // return empty string - return std::string(); - } - - // get the delta unit between two datetimes in the time col - size_t seconds_diff = difftime( mktime( &dtinfo2.time ), - mktime( &dtinfo1.time ) ); - - if ( seconds_diff == 0 ) { - seconds_diff = 1; //if millisec, want some update - } - - // increment the time and format - dtinfo2.time.tm_sec += tp * seconds_diff; - - int err = mktime( &dtinfo2.time ); - - if ( err < 0 ) { - std::stringstream errMsg; - errMsg << "increment_datetime_str() mktime failed on " - << datetime2; - throw( errMsg.str() ); - } - - // format incremented time - char tmp_buffer [ BUFSIZ ]; - - size_t n_char = strftime( tmp_buffer, BUFSIZ, - dtinfo2.datetime_fmt.c_str(), &dtinfo2.time ); - if ( n_char == 0 ) { - std::stringstream errMsg; - errMsg << "increment_datetime_str(): Failed on " - << datetime1 << ", " << datetime2 << " tp = " << tp; - throw( errMsg.str() ); - } - - return std::string( tmp_buffer ); -} diff --git a/src/EDM.cc b/src/EDM.cc new file mode 100644 index 0000000..8a7e536 --- /dev/null +++ b/src/EDM.cc @@ -0,0 +1,21 @@ +// Contains class definitions for EDM object. See EDM.h header for +// description of EDM class and how to use the methods. + +#include "EDM.h" + +//---------------------------------------------------------------- +// EDM() : Constructor +// +// data : Input dataframe containing the time series to model. +// embedded : Flag on whether the input data is already embedded. +// E : The embedding dimension to use when performing embedding. +// tau : The steps between each index in the embedding. Negative +// tau is probably what you intend to use; positive tau yields +// an embedding (E_t,E_t+tau...) (future forward embedding). +//---------------------------------------------------------------- +EDM::EDM ( DataFrame & data, bool embedded, int E, int tau ) : + data( data ) { + + // check parameters + +} diff --git a/src/EDM.h b/src/EDM.h new file mode 100644 index 0000000..fca7c9b --- /dev/null +++ b/src/EDM.h @@ -0,0 +1,35 @@ +// Contains class declaration for EDM object. The EDM object is a general +// data-processing class which holds the data object and its embedding, +// and will have its specific projection methods defined by Simplex and SMap. + +#ifndef EDM_H +#define EDM_H + +#include "Common.h" + +//---------------------------------------------------------------- +// EDM Class +// EDM maintains the central data object used for prediction and +// defines data iterating patterns common to Simplex and SMap. +// It should be treated as an abstract class as it performs no +// prediction, but for sake of compiler simplicity we are choosing +// to leave it as a virtual class for now +//---------------------------------------------------------------- +class EDM { + + // The input time series (potentially not embedded) + DataFrame & data; + + // The target column this class should model for + std::string targetName; + + public: + + //---------------------------------------------------------------- + // EDM() : Constructor + //---------------------------------------------------------------- + EDM ( DataFrame & data, bool embedded, int E, int tau ); + +}; + +#endif diff --git a/src/EDM_Functions.cc b/src/EDM_Functions.cc new file mode 100644 index 0000000..a6f8068 --- /dev/null +++ b/src/EDM_Functions.cc @@ -0,0 +1,39 @@ +// These functions are function-wrappers for the EDM algorithm classes +// such that the user does not have to deal with the OOD. + +#include "Common.h" + +//---------------------------------------------------------------- +// Simplex prediction algorithm with data input as DataFrame +// See EDM class for parameter definitions. +//---------------------------------------------------------------- +DataFrame Simplex( DataFrame & data, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ) { + + // Create Simplex object and get its projection + +} + +//---------------------------------------------------------------- +// Simplex prediction algorithm with data input as filepath +// See EDM class for parameter definitions. +//---------------------------------------------------------------- +DataFrame Simplex( std::string pathIn, std::string dataFile, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ) { + + // Just load in dataframe and delegate to Simplex with DataFrame input + + DataFrame data (pathIn, dataFile); + + return Simplex ( data, pathOut, predictFile, lib, pred, E, Tp, knn, tau, + exclusionRadius, columns, target, embedded, const_predict, verbose ); +} + diff --git a/src/Embed.cc b/src/Embed.cc deleted file mode 100644 index c1fba55..0000000 --- a/src/Embed.cc +++ /dev/null @@ -1,183 +0,0 @@ - -#include "Embed.h" - -// NOTE: The returned data block does NOT have the time column - -//---------------------------------------------------------------- -// API Overload 1: Explicit data file path/name -// Embed DataFrame columns in E dimensions. -// Data is read from path/dataFile -// Implemented as a wrapper to API Overload 2: -// which is a wrapper for MakeBlock() -// -// NOTE: Truncates data by tau * (E-1) rows to remove -// nan values (partial data rows) -//---------------------------------------------------------------- -DataFrame< double > Embed( std::string path, - std::string dataFile, - int E, // embedding dimension - int tau, // time step delay - std::string columns, // column names or indices - bool verbose ) { - - DataFrame< double > dataFrame( path, dataFile ); - DataFrame< double > embedded = Embed(dataFrame, E, tau, columns, verbose); - return embedded; -} - -//---------------------------------------------------------------- -// API Overload 2: DataFrame provided -// Implemented as a wrapper for MakeBlock() -// Note: dataFrame must have the columnNameToIndex map -//---------------------------------------------------------------- -DataFrame< double > Embed( DataFrame< double > dataFrameIn, - int E, - int tau, - std::string columns, - bool verbose ) { - - // Parameter.Validate will convert columns into a vector of names - // or a vector of column indices - Parameters param = Parameters( Method::Embed, "", "", "", "", - "1 1", "1 1", E, 0, 0, tau, 0, 0, - columns, "", false, false, verbose ); - - if ( not param.columnIndex.size() and - dataFrameIn.ColumnNameToIndex().empty() ) { - throw std::runtime_error("Embed(DataFrame): columnNameIndex empty.\n"); - } - - // If columns provided, validate they are in dataFrameIn - for ( auto colName : param.columnNames ) { - auto ci = find( dataFrameIn.ColumnNames().begin(), - dataFrameIn.ColumnNames().end(), colName ); - - if ( ci == dataFrameIn.ColumnNames().end() ) { - std::stringstream errMsg; - errMsg << "Embed(DataFrame): Failed to find column " - << colName << " in dataFrame with columns: [ "; - for ( auto col : dataFrameIn.ColumnNames() ) { - errMsg << col << " "; - } errMsg << " ]\n"; - throw std::runtime_error( errMsg.str() ); - } - } - - // Get column names for MakeBlock - std::vector< std::string > colNames; - if ( param.columnNames.size() ) { - // column names are strings use as-is - colNames = param.columnNames; - } - else if ( param.columnIndex.size() ) { - // columns are indices : Create column names for MakeBlock - for ( size_t i = 0; i < param.columnIndex.size(); i++ ) { - std::stringstream ss; - ss << "V" << param.columnIndex[i]; - colNames.push_back( ss.str() ); - } - } - else { - throw std::runtime_error( "Embed(DataFrame): columnNames and " - " columnIndex are empty.\n" ); - } - - // Extract the specified columns (sub)DataFrame from dataFrameIn - DataFrame< double > dataFrame; - - if ( param.columnNames.size() ) { - dataFrame = dataFrameIn.DataFrameFromColumnNames( param.columnNames ); - } - else if ( param.columnIndex.size() ) { - // already have column indices - // Note there will be no column names transferred - dataFrame = dataFrameIn.DataFrameFromColumnIndex( param.columnIndex ); - } - - DataFrame< double > embedding = MakeBlock( dataFrame, E, tau, - colNames, verbose ); - - return embedding; -} - -//------------------------------------------------------------------------ -// MakeBlock from dataFrame -// Ignores the first (or last) tau * (E-1) dataFrame rows of partial data. -// Does not validate parameters or columns, use Embed() -//------------------------------------------------------------------------ -DataFrame< double > MakeBlock( DataFrame< double > dataFrame, - int E, - int tau, - std::vector columnNames, - bool verbose ) { - - if ( columnNames.size() != dataFrame.NColumns() ) { - std::stringstream errMsg; - errMsg << "MakeBlock: The number of columns in the dataFrame (" - << dataFrame.NColumns() << ") is not equal to the number " - << "of columns specified (" << columnNames.size() << ").\n";; - throw std::runtime_error( errMsg.str() ); - } - - if ( E < 1 ) { - std::stringstream errMsg; - errMsg << "MakeBlock(): E = " << E << " is invalid.\n" ; - throw std::runtime_error( errMsg.str() ); - } - - size_t NRows = dataFrame.NRows(); // number of input rows - size_t NColOut = dataFrame.NColumns() * E; // number of output columns - size_t NPartial = abs( tau ) * (E-1); // rows to shift & delete - - // Create embedded data frame column names X(t-0) X(t-1)... - std::vector< std::string > newColumnNames( NColOut ); - size_t newCol_i = 0; - for ( size_t col = 0; col < columnNames.size(); col ++ ) { - for ( size_t e = 0; e < E; e++ ) { - std::stringstream ss; - if ( tau < 0 ) { - ss << columnNames[ col ] << "(t-" << e << ")"; - } - else { - ss << columnNames[ col ] << "(t+" << e << ")"; - } - newColumnNames[ newCol_i ] = ss.str(); - newCol_i++; - } - } - - // Ouput data frame with tau * E-1 fewer rows - DataFrame< double > embedding( NRows - NPartial, NColOut, newColumnNames ); - - // To keep track of where to insert column in new data frame - size_t colCount = 0; - - // Slice to ignore rows with partial data - std::slice slice_i; - if ( tau < 0 ) { - slice_i = std::slice( NPartial, NRows - NPartial, 1 ); - } - else { - slice_i = std::slice( 0, NRows - NPartial, 1 ); - } - - // Shift column data and write to embedding data frame - for ( size_t col = 0; col < dataFrame.NColumns(); col++ ) { - // for each embedding dimension - for ( size_t e = 0; e < E; e++ ) { - - std::valarray< double > column = dataFrame.Column( col ); - - // Returns a copy of the valarray object with its elements - // shifted left n spaces (or right if n is negative). - std::valarray< double > tmp = column.shift( e * tau ); - - // Write shifted columns to the output embedding DataFrame - embedding.WriteColumn( colCount, tmp[ slice_i ] ); - - colCount++; - } - } - - return embedding; -} diff --git a/src/Embed.h b/src/Embed.h deleted file mode 100644 index 3bccba0..0000000 --- a/src/Embed.h +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef EMBED_H -#define EMBED_H - -#include "Common.h" -#include "Parameter.h" - -//---------------------------------------------------------------- -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: -// which is a wrapper for MakeBlock() -//---------------------------------------------------------------- -DataFrame< double > Embed ( std::string path = "", - std::string dataFile = "", - int E = 0, - int tau = 0, - std::string columns = "", - bool verbose = false ); - -//---------------------------------------------------------------- -// API Overload 2: DataFrame provided -// Implemented as a wrapper for MakeBlock() -//---------------------------------------------------------------- -DataFrame< double > Embed ( DataFrame< double > dataFrame, - int E = 0, - int tau = 0, - std::string columns = "", - bool verbose = false ); - -//---------------------------------------------------------------- -//---------------------------------------------------------------- -DataFrame< double > MakeBlock ( DataFrame< double > dataFrame, - int E, - int tau, - std::vector columnNames, - bool verbose ); -#endif diff --git a/src/Eval.cc b/src/Eval.cc deleted file mode 100644 index d974463..0000000 --- a/src/Eval.cc +++ /dev/null @@ -1,687 +0,0 @@ - -#include -#include -#include -#include - -#include "Common.h" - -namespace EDM_Eval { - // Thread Work Queue : Vector of int - typedef std::vector< int > WorkQueue; - - // Thread exception_ptr queue - std::queue< std::exception_ptr > embedDimExceptQ; - std::queue< std::exception_ptr > predictIntExceptQ; - std::queue< std::exception_ptr > predictNLExceptQ; - - // atomic counters for all threads - std::atomic tp_count_i (0); // initialize to 0 - std::atomic embed_count_i(0); // initialize to 0 - std::atomic smap_count_i (0); // initialize to 0 - - std::mutex mtx; - std::mutex q_mtx; -} - -//---------------------------------------------------------------- -// Forward declaration: -// Worker thread for EmbedDimension() -//---------------------------------------------------------------- -void EmbedThread( EDM_Eval::WorkQueue &workQ, - DataFrame< double > &data, - DataFrame< double > &E_rho, - std::string lib, - std::string pred, - int Tp, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose ); - -//---------------------------------------------------------------- -// Forward declaration: -// Worker thread for PredictInterval() -//---------------------------------------------------------------- -void PredictIntervalThread( EDM_Eval::WorkQueue &workQ, - DataFrame< double > &data, - DataFrame< double > &Tp_rho, - std::string lib, - std::string pred, - int E, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose ); - -//---------------------------------------------------------------- -// Forward declaration: -// Worker thread for PredictNonLinear() -//---------------------------------------------------------------- -void SMapThread( EDM_Eval::WorkQueue &workQ, - DataFrame< double > &data, - DataFrame< double > &Theta_rho, - std::vector ThetaValues, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose ); - -//---------------------------------------------------------------- -// EmbedDimension() : Evaluate Simplex rho vs. dimension E -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: -//---------------------------------------------------------------- -DataFrame EmbedDimension( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int maxE, - int Tp, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose, - unsigned nThreads ) { - - // Create DataFrame (constructor loads data) - DataFrame< double > *dataFrameIn = - new DataFrame< double > ( pathIn, dataFile ); - - DataFrame E_rho = EmbedDimension( std::ref( *dataFrameIn ), - pathOut, - predictFile, - lib, - pred, - maxE, - Tp, - tau, - colNames, - targetName, - embedded, - verbose, - nThreads ); - - delete dataFrameIn; - - return E_rho; -} - -//---------------------------------------------------------------- -// EmbedDimension() : Evaluate Simplex rho vs. dimension E -// API Overload 2: DataFrame provided -//---------------------------------------------------------------- -DataFrame EmbedDimension( DataFrame< double > &data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int maxE, - int Tp, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose, - unsigned nThreads ) { - - // Container for results - DataFrame E_rho( maxE, 2, "E rho" ); - - // Build work queue - EDM_Eval::WorkQueue workQ( maxE ); - - // Insert dimension values into work queue - for ( auto i = 0; i < maxE; i++ ) { - workQ[ i ] = i + 1; - } - - unsigned maxThreads = std::thread::hardware_concurrency(); - if ( maxThreads < nThreads ) { nThreads = maxThreads; } - if ( nThreads > maxE ) { nThreads = maxE; } - - // thread container - std::vector< std::thread > threads; - for ( unsigned i = 0; i < nThreads; i++ ) { - - threads.push_back( std::thread( EmbedThread, - std::ref( workQ ), - std::ref( data ), - std::ref( E_rho ), - lib, - pred, - Tp, - tau, - colNames, - targetName, - embedded, - verbose ) ); - } - - // join threads - for ( auto &thrd : threads ) { - thrd.join(); - } - - // If thread threw exception, get from queue and rethrow - if ( not EDM_Eval::embedDimExceptQ.empty() ) { - std::lock_guard lck( EDM_Eval::q_mtx ); - - // Take the first exception in the queue - std::exception_ptr exceptionPtr = EDM_Eval::embedDimExceptQ.front(); - - // Unroll all other exception from the thread/loops - while( not EDM_Eval::embedDimExceptQ.empty() ) { - // JP When do these exception_ptr get deleted? Is it a leak? - EDM_Eval::embedDimExceptQ.pop(); - } - std::rethrow_exception( exceptionPtr ); - } - - if ( predictFile.size() ) { - E_rho.WriteData( pathOut, predictFile ); - } - - return E_rho; -} - -//---------------------------------------------------------------- -// Worker thread for EmbedDimension() -//---------------------------------------------------------------- -void EmbedThread( EDM_Eval::WorkQueue &workQ, - DataFrame< double > &data, - DataFrame< double > &E_rho, - std::string lib, - std::string pred, - int Tp, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose ) -{ - - std::size_t i = std::atomic_fetch_add( &EDM_Eval::embed_count_i, - std::size_t(1) ); - - while( i < workQ.size() ) { - - // WorkQueue stores E - int E = workQ[ i ]; - - // Simplex() -> EmbedNN() -> Embed() -> DeletePartialDataRows() - // In a multthreaded application we need to pass a unique copy of - // the data so that DeletePartialDataRows() is not recursively - // applied to the same data frame. - DataFrame< double > localData( data ); - - try { - DataFrame S = Simplex( std::ref( localData ), - "", // pathOut, - "", // predictFile, - lib, - pred, - E, - Tp, - 0, // knn - tau, - 0, // exclusionRadius - colNames, - targetName, - embedded, - false, // const_predict - verbose ); - - VectorError ve = ComputeError( S.VectorColumnName("Observations"), - S.VectorColumnName("Predictions")); - - E_rho.WriteRow( i, std::valarray({ (double) E, ve.rho })); - - if ( verbose ) { - std::lock_guard lck( EDM_Eval::mtx ); - std::cout << "EmbedThread() workQ[" << workQ[i] << "] E " << E - << " rho " << ve.rho << " RMSE " << ve.RMSE - << " MAE " << ve.MAE << std::endl << std::endl; - } - } - catch(...) { - // push exception pointer onto queue for main thread to catch - std::lock_guard lck( EDM_Eval::q_mtx ); - EDM_Eval::embedDimExceptQ.push( std::current_exception() ); - } - - i = std::atomic_fetch_add(&EDM_Eval::embed_count_i, std::size_t(1)); - } - - // Reset counter - std::atomic_store( &EDM_Eval::embed_count_i, std::size_t(0) ); -} - -//----------------------------------------------------------------- -// PredictInterval() : Evaluate Simplex rho vs. predict interval Tp -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: -//----------------------------------------------------------------- -DataFrame PredictInterval( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int maxTp, - int E, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose, - unsigned nThreads ) { - - // Create DataFrame (constructor loads data) - DataFrame< double > *dataFrameIn = - new DataFrame< double > ( pathIn, dataFile ); - - DataFrame Tp_rho = PredictInterval( std::ref( *dataFrameIn ), - pathOut, - predictFile, - lib, - pred, - maxTp, - E, - tau, - colNames, - targetName, - embedded, - verbose ); - delete dataFrameIn; - - return Tp_rho; -} - -//----------------------------------------------------------------- -// PredictInterval() : Evaluate Simplex rho vs. predict interval Tp -// API Overload 2: DataFrame provided -//----------------------------------------------------------------- -DataFrame PredictInterval( DataFrame< double > &data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int maxTp, - int E, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose, - unsigned nThreads ) { - - // Container for results - DataFrame Tp_rho( maxTp, 2, "Tp rho" ); - - // Build work queue - EDM_Eval::WorkQueue workQ( maxTp ); - - // Insert Tp values into work queue - for ( auto i = 0; i < maxTp; i++ ) { - workQ[ i ] = i + 1; - } - - unsigned maxThreads = std::thread::hardware_concurrency(); - if ( maxThreads < nThreads ) { nThreads = maxThreads; } - if ( nThreads > maxTp ) { nThreads = maxTp; } - - // thread container - std::vector< std::thread > threads; - for ( unsigned i = 0; i < nThreads; ++i ) { - threads.push_back( std::thread( PredictIntervalThread, - std::ref( workQ ), - std::ref( data ), - std::ref( Tp_rho ), - lib, - pred, - E, - tau, - colNames, - targetName, - embedded, - verbose ) ); - } - - // join threads - for ( auto &thrd : threads ) { - thrd.join(); - } - - // If thread threw exception, get from queue and rethrow - if ( not EDM_Eval::predictIntExceptQ.empty() ) { - std::lock_guard lck( EDM_Eval::q_mtx ); - - // Take the first exception in the queue - std::exception_ptr exceptionPtr = EDM_Eval::predictIntExceptQ.front(); - - // Unroll all other exception from the thread/loops - while( not EDM_Eval::predictIntExceptQ.empty() ) { - // JP When do these exception_ptr get deleted? Is it a leak? - EDM_Eval::predictIntExceptQ.pop(); - } - std::rethrow_exception( exceptionPtr ); - } - - if ( predictFile.size() ) { - Tp_rho.WriteData( pathOut, predictFile ); - } - - return Tp_rho; -} - -//---------------------------------------------------------------- -// Worker thread for PredictInterval() -//---------------------------------------------------------------- -void PredictIntervalThread( EDM_Eval::WorkQueue &workQ, - DataFrame< double > &data, - DataFrame< double > &Tp_rho, - std::string lib, - std::string pred, - int E, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose ) -{ - std::size_t i = std::atomic_fetch_add( &EDM_Eval::tp_count_i, - std::size_t(1) ); - - while( i < workQ.size() ) { - - // WorkQueue stores Tp - int Tp = workQ[ i ]; - - // In a multthreaded application we need to pass a unique copy of - // the data so that DeletePartialDataRows() is not recursively - // applied to the same data frame. - DataFrame< double > localData( data ); - - try { - DataFrame S = Simplex( std::ref( localData ), - "", // pathOut, - "", // predictFile, - lib, - pred, - E, - Tp, - 0, // knn - tau, - 0, // exclusionRadius - colNames, - targetName, - embedded, - false, // const_pred - verbose ); - - VectorError ve = ComputeError( S.VectorColumnName("Observations"), - S.VectorColumnName("Predictions")); - - Tp_rho.WriteRow( i, std::valarray({ (double) Tp, ve.rho })); - - if ( verbose ) { - std::lock_guard lck( EDM_Eval::mtx ); - std::cout << "PredictIntervalThread() workQ[" << workQ[i] - << "] Tp " << Tp - << " rho " << ve.rho << " RMSE " << ve.RMSE - << " MAE " << ve.MAE << std::endl << std::endl; - } - } - catch(...) { - // push exception pointer onto queue for main thread to catch - std::lock_guard lck( EDM_Eval::q_mtx ); - EDM_Eval::predictIntExceptQ.push( std::current_exception() ); - } - - i = std::atomic_fetch_add( &EDM_Eval::tp_count_i, std::size_t(1) ); - } - - // Reset counter - std::atomic_store( &EDM_Eval::tp_count_i, std::size_t(0) ); -} - -//---------------------------------------------------------------- -// PredictNonlinear() : Smap rho vs. localisation parameter theta -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: -//---------------------------------------------------------------- -DataFrame PredictNonlinear( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - std::string theta, - int E, - int Tp, - int knn, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose, - unsigned nThreads ) { - - // Create DataFrame (constructor loads data) - DataFrame< double > *dataFrameIn = - new DataFrame< double > ( pathIn, dataFile ); - - DataFrame< double > Theta_rho = PredictNonlinear( std::ref( *dataFrameIn ), - pathOut, - predictFile, - lib, - pred, - theta, - E, - Tp, - knn, - tau, - colNames, - targetName, - embedded, - verbose ); - delete dataFrameIn; - - return Theta_rho; -} - -//---------------------------------------------------------------- -// PredictNonlinear() : Smap rho vs. localisation parameter theta -// API Overload 2: DataFrame provided -//---------------------------------------------------------------- -DataFrame PredictNonlinear( DataFrame< double > &data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - std::string theta, - int E, - int Tp, - int knn, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose, - unsigned nThreads ) { - - std::vector ThetaValues( { 0.01, 0.1, 0.3, 0.5, 0.75, 1, - 1.5, 2, 3, 4, 5, 6, 7, 8, 9 } ); - - if ( theta.size() ) { - // Use theta values passed in as parameter string - ThetaValues.clear(); - - std::vector theta_vec = SplitString( theta, " \t,\n" ); - - try { - for ( auto ci = theta_vec.begin(); ci != theta_vec.end(); ++ci ) { - ThetaValues.push_back( std::stod( *ci ) ); - } - } - catch ( const std::invalid_argument &ia ) { - std::stringstream errMsg; - errMsg << "PredictNonlinear(): Unable to convert theta [" - << ia.what() << "] to numeric."; - throw std::runtime_error( errMsg.str() ); - } - } - - // Container for results - DataFrame Theta_rho( ThetaValues.size(), 2, "Theta rho" ); - - // Build work queue - EDM_Eval::WorkQueue workQ( ThetaValues.size() ); - - // Insert ThetaValues indexes into work queue - for ( auto i = 0; i < ThetaValues.size(); i++ ) { - workQ[ i ] = i; - } - - unsigned maxThreads = std::thread::hardware_concurrency(); - if ( maxThreads < nThreads ) { nThreads = maxThreads; } - if ( nThreads > ThetaValues.size() ) { nThreads = ThetaValues.size(); } - - // thread container - std::vector< std::thread > threads; - for ( unsigned i = 0; i < nThreads; ++i ) { - threads.push_back( std::thread( SMapThread, - std::ref( workQ ), - std::ref( data ), - std::ref( Theta_rho ), - ThetaValues, - lib, - pred, - E, - Tp, - knn, - tau, - colNames, - targetName, - embedded, - verbose ) ); - } - - // join threads - for ( auto &thrd : threads ) { - thrd.join(); - } - - // If thread threw exception, get from queue and rethrow - if ( not EDM_Eval::predictNLExceptQ.empty() ) { - std::lock_guard lck( EDM_Eval::q_mtx ); - - // Take the first exception in the queue - std::exception_ptr exceptionPtr = EDM_Eval::predictNLExceptQ.front(); - - // Unroll all other exception from the thread/loops - while( not EDM_Eval::predictNLExceptQ.empty() ) { - // JP When do these exception_ptr get deleted? Is it a leak? - EDM_Eval::predictNLExceptQ.pop(); - } - std::rethrow_exception( exceptionPtr ); - } - - if ( predictFile.size() ) { - Theta_rho.WriteData( pathOut, predictFile ); - } - - return Theta_rho; -} - -//---------------------------------------------------------------- -// Worker thread for PredictNonlinear() -//---------------------------------------------------------------- -void SMapThread( EDM_Eval::WorkQueue &workQ, - DataFrame< double > &data, - DataFrame< double > &Theta_rho, - std::vector ThetaValues, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - std::string colNames, - std::string targetName, - bool embedded, - bool verbose ) -{ - - std::size_t i = - std::atomic_fetch_add( &EDM_Eval::smap_count_i, std::size_t(1) ); - - while( i < workQ.size() ) { - - double theta = ThetaValues[ workQ[ i ] ]; - - // In a multthreaded application we need to pass a unique copy of - // the data so that DeletePartialDataRows() is not recursively - // applied to the same data frame. - DataFrame< double > localData( data ); - - try { - SMapValues S = SMap( std::ref( localData ), - "", - "", // predictFile - lib, - pred, - E, - Tp, - knn, - tau, - theta, - 0, // exclusionRadius - colNames, - targetName, - "", // smapFile - "", // derivatives - embedded, - false, // const_predict - verbose ); - - DataFrame< double > predictions = S.predictions; - DataFrame< double > coefficients = S.coefficients; - - VectorError ve = ComputeError( - predictions.VectorColumnName( "Observations" ), - predictions.VectorColumnName( "Predictions" ) ); - - Theta_rho.WriteRow( i, std::valarray({ theta, ve.rho })); - - if ( verbose ) { - std::lock_guard lck( EDM_Eval::mtx ); - std::cout << "Theta " << theta - << " rho " << ve.rho << " RMSE " << ve.RMSE - << " MAE " << ve.MAE << std::endl << std::endl; - } - } - catch(...) { - // push exception pointer onto queue for main thread to catch - std::lock_guard lck( EDM_Eval::q_mtx ); - EDM_Eval::predictNLExceptQ.push( std::current_exception() ); - } - - i = std::atomic_fetch_add( &EDM_Eval::smap_count_i, std::size_t(1) ); - } - - // Reset counter - std::atomic_store( &EDM_Eval::smap_count_i, std::size_t(0) ); -} diff --git a/src/Interface.cc b/src/Interface.cc deleted file mode 100644 index 17c2671..0000000 --- a/src/Interface.cc +++ /dev/null @@ -1,7 +0,0 @@ -#include "Common.h" - -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -void Interface() { -} diff --git a/src/Multiview.cc b/src/Multiview.cc deleted file mode 100644 index 9050b35..0000000 --- a/src/Multiview.cc +++ /dev/null @@ -1,769 +0,0 @@ - -//-------------------------------------------------------------------- -// Data input requires columns to specify timeseries columns -// that will be embedded by Embed(), and target for predictions. -// -// D represents the number of variables to combine for each -// assessment, if not specified, it is the number of columns. -// E is the embedding dimension of each variable. -// If E = 1, no time delay embedding is done, but the variables -// in the embedding are named X(t-0), Y(t-0)... -// -// Parameters.Validate() with Method::Simplex sets knn equal to E+1 -// if knn not specified, so we need to explicitly set knn to D + 1. -// -// multiview is the number of top-ranked D-dimensional predictions -// to "average" for the final prediction. Corresponds to parameter -// k in Ye & Sugihara with default k = sqrt(m) where m is the -// number of combinations C(n,D) available from the n = D * E -// columns taken D at-a-time. -// -// Ye H., and G. Sugihara, 2016. Information leverage in -// interconnected ecosystems: Overcoming the curse of dimensionality. -// Science 353:922–925. -//-------------------------------------------------------------------- -// -// NOTE: Multiview evaluates the top projections using in-sample -// library predictions. It can be shown that highly accurate -// in-sample predictions can be made from arbitrary non- -// constant, non-oscillatory vectors. Therefore, some attention -// may be warranted to filter prospective embedding vectors. -//-------------------------------------------------------------------- - -#include -#include -#include -#include - -#include "Common.h" -#include "AuxFunc.h" - -namespace EDM_Multiview { - // Thread Work Queue : Vector of combos indices - typedef std::vector< int > WorkQueue; - - // Thread exception_ptr queue - std::queue< std::exception_ptr > exceptionQ; - - // atomic counter for all threads - std::atomic eval_i(0); // initialize to 0 - - std::mutex mtx; - std::mutex q_mtx; -} - -//---------------------------------------------------------------- -// forward declarations -//---------------------------------------------------------------- -std::vector< std::vector< size_t > > Combination( int n, int k ); - -DataFrame SimplexProjection( Parameters param, - DataEmbedNN embedNN, - bool checkDataRows = true ); - -void EvalComboThread( Parameters param, - EDM_Multiview::WorkQueue workQ, - std::vector< std::vector< size_t > > combos, - DataFrame< double > &embedding, - std::valarray< double > &targetVec, - DataFrame< double > &combos_rho, - std::vector< DataFrame< double > > &prediction ); - -std::vector< std::string > ComboRhoTable( DataFrame combosRho, - std::vector< std::string > colNames ); - -//---------------------------------------------------------------- -// Multiview() : Evaluate Simplex rho vs. dimension E -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: -//---------------------------------------------------------------- -MultiviewValues Multiview( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int D, - int E, - int Tp, - int knn, - int tau, - std::string columns, - std::string target, - int multiview, - int exclusionRadius, - bool trainLib, - bool verbose, - unsigned nThreads ) { - - // Create DataFrame (constructor loads data) - DataFrame< double > dataFrameIn( pathIn, dataFile ); - - MultiviewValues result = Multiview( dataFrameIn, - pathOut, - predictFile, - lib, - pred, - D, - E, - Tp, - knn, - tau, - columns, - target, - multiview, - exclusionRadius, - trainLib, - verbose, - nThreads ); - return result; -} - -//---------------------------------------------------------------- -// Multiview() -// API Overload 2: DataFrame provided -//---------------------------------------------------------------- -MultiviewValues Multiview( DataFrame< double > data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int D, - int E, - int Tp, - int knn, - int tau, - std::string columns, - std::string target, - int multiview, - int exclusionRadius, - bool trainLib, - bool verbose, - unsigned nThreads ) { - - // Require at least E = 1 - if ( E < 1 ) { - std::stringstream errMsg; - errMsg << " Multiview(): E = " << E << " is invalid.\n" ; - throw std::runtime_error( errMsg.str() ); - } - - // Create local Parameters struct. Note embedded = true - Parameters param = Parameters( Method::Simplex, "", "", - pathOut, predictFile, - lib, pred, E, Tp, knn, tau, 0, - exclusionRadius, columns, target, - true, // embedded true - false, verbose, - "", "", "", 0, 0, 0, multiview ); - - if ( not param.columnNames.size() ) { - throw std::runtime_error( "Multiview() requires column names." ); - } - if ( not param.targetName.size() ) { - throw std::runtime_error( "Multiview() requires target name." ); - } - // Ensure that params are validated so columnNames are populated - if ( not param.validated ) { - throw std::runtime_error( "Multiview() params not validated." ); - } - - // Validate that columns & target are in data - for ( auto colName : param.columnNames ) { - auto ci = find( data.ColumnNames().begin(), - data.ColumnNames().end(), colName ); - - if ( ci == data.ColumnNames().end() ) { - std::stringstream errMsg; - errMsg << "Multiview(): Failed to find column " - << colName << " in dataFrame with columns: [ "; - for ( auto col : data.ColumnNames() ) { - errMsg << col << " "; - } errMsg << " ]\n"; - throw std::runtime_error( errMsg.str() ); - } - } - auto ti = find( data.ColumnNames().begin(), - data.ColumnNames().end(), param.targetName ); - if ( ti == data.ColumnNames().end() ) { - std::stringstream errMsg; - errMsg << "Multiview(): Failed to find target " - << param.targetName << " in dataFrame with columns: [ "; - for ( auto col : data.ColumnNames() ) { - errMsg << col << " "; - } errMsg << " ]\n"; - throw std::runtime_error( errMsg.str() ); - } - - // Validate data rows against lib and pred indices - CheckDataRows( param, std::ref( data ), "Multiview()" ); - - //------------------------------------------------------------ - // Generate embedding on param.columns_str - // embedding will have tau * (E-1) fewer rows than data, - // the time column is not returned in embedding. - // Column names are changed from X to X(t-0), X(t-1),... - //------------------------------------------------------------ - DataFrame< double > embedding = Embed( data, - param.E, - param.tau, - param.columns_str, - param.verbose ); - - size_t shift = abs( param.tau ) * ( param.E - 1 ); - - // Delete data top rows of partial data - if ( not data.PartialDataRowsDeleted() ) { - // Not thread safe - std::lock_guard lck( EDM_Multiview::mtx ); - - data.DeletePartialDataRows( shift, param.tau ); - } - - // Adjust param.library and param.prediction vectors of indices - if ( shift > 0 ) { - param.DeleteLibPred( shift ); - } - - // Get target - std::valarray< double > targetVec = data.VectorColumnName(param.targetName); - - // Save param.predictOutputFile and reset so Simplex() does not write - std::string outputFile = param.predictOutputFile; - param.predictOutputFile = ""; - - // Establish the state-space dimension D - // default to the number of input columns (not embedded columns) - if ( D == 0 ) { - D = param.columnNames.size(); - } - if ( D > embedding.NColumns() ) { - std::stringstream msg; - msg << "WARNING: Multiview(): D = " << D - << " exceeds the number of columns in the embedding: " - << embedding.NColumns() << ". D set to " - << embedding.NColumns() << std::endl; - std::cout << msg.str(); - - D = param.columnNames.size(); - } - - // Combinations of possible embedding variables, D at-a-time - // Note that these combinations are not zero-offset, i.e. - // Combination( 3, 2 ) = [(1, 2), (1, 3), (2, 3)] - // These correspond to column indices +1 - std::vector< std::vector< size_t > > combos = - Combination( embedding.NColumns(), D ); - -#ifdef DEBUG_ALL - std::cout << "Multiview(): " << combos.size() << " combos:\n"; - for ( auto i = 0; i < combos.size(); i++ ) { - std::vector< size_t > combo_i = combos[i]; - std::cout << "["; - for ( auto j = 0; j < combo_i.size(); j++ ) { - std::cout << combo_i[j] << ","; - } - std::cout << "] "; - } std::cout << std::endl; -#endif - - // Establish number of ensembles if not specified - if ( not param.MultiviewEnsemble ) { - // Ye & Sugihara suggest sqrt( m ) as the number of embeddings to avg - param.MultiviewEnsemble = std::max(2, (int) std::sqrt(combos.size())); - - std::stringstream msg; - msg << "Multiview() Set view sample size to " - << param.MultiviewEnsemble << std::endl; - std::cout << msg.str(); - } - - // validate number of combinations - if ( param.MultiviewEnsemble > combos.size() ) { - std::stringstream msg; - msg << "WARNING: Multiview(): multiview ensembles " - << param.MultiviewEnsemble - << " exceeds the number of available combinations: " - << combos.size() << ". Set to " << combos.size() << std::endl; - std::cout << msg.str(); - - param.MultiviewEnsemble = combos.size(); - } - - //--------------------------------------------------------------- - // Evaluate variable combinations. - // Note that this is done within the library itself (in-sample). - //--------------------------------------------------------------- - // Save a copy of the specified prediction observation rows. - std::vector prediction = param.prediction; - - if ( trainLib ) { - // Override param.prediction for in-sample forecast skill evaluation - param.prediction = param.library; - } - - // This is not a good implementation... - // Replace param.E with the number of dimensions, recall embbeded = true - param.E = D; - - // Create column names for the results DataFrame - // One row for each combo: D columns (a combo), rho, MAE, RMSE - // JP: NOTE DataFrame is based on valarray, so we can't store non - // numeric values (column names) in the results DataFrame. - std::stringstream header; - for ( auto i = 1; i <= D; i++ ) { - header << "Col_" << i << " "; - } - header << "rho MAE RMSE"; - - // Results Data Frame: D columns (a combo), rho, mae, rmse - DataFrame combos_rho( combos.size(), D + 3, header.str() ); - - // Results vector of DataFrame's with prediction results - std::vector< DataFrame< double > > combos_prediction( combos.size() ); - - // Build work queue - EDM_Multiview::WorkQueue workQ( combos.size() ); - - // Insert combos index into work queue - for ( auto i = 0; i < combos.size(); i++ ) { - workQ[ i ] = i; - } - - unsigned maxThreads = std::thread::hardware_concurrency(); - if ( maxThreads < nThreads ) { nThreads = maxThreads; } - - // thread container - std::vector< std::thread > threads; - for ( unsigned i = 0; i < nThreads; ++i ) { - threads.push_back( std::thread( EvalComboThread, - param, - workQ, - combos, - std::ref( embedding ), - std::ref( targetVec ), - std::ref( combos_rho ), - std::ref( combos_prediction ) ) ); - } - - // join threads - for ( auto &thrd : threads ) { - thrd.join(); - } - - // If thread threw exception, get from queue and rethrow - if ( not EDM_Multiview::exceptionQ.empty() ) { - std::lock_guard lck( EDM_Multiview::q_mtx ); - - // Take the first exception in the queue - std::exception_ptr exceptionPtr = EDM_Multiview::exceptionQ.front(); - - // Unroll all other exception from the thread/loops - while( not EDM_Multiview::exceptionQ.empty() ) { - // JP When do these exception_ptr get deleted? Is it a leak? - EDM_Multiview::exceptionQ.pop(); - } - std::rethrow_exception( exceptionPtr ); - } - - //----------------------------------------------------------------- - // Rank forecasts. If trainLib true these are in-sample (library) - //----------------------------------------------------------------- - // Make pairs of row indices and rho - std::valarray< double > rho = combos_rho.VectorColumnName( "rho" ); - // vector of indices - std::valarray< size_t > indices( rho.size() ); - std::iota( begin( indices ), end( indices ), 0 ); - - // Ensure that rho is the first of the pair so sort will work - std::vector< std::pair< double, int > > combo_sort( rho.size() ); - for ( size_t i = 0; i < rho.size(); i++ ) { - combo_sort[ i ] = std::make_pair( rho[i], indices[i] ); - } - - // sort pairs and reverse for largest rho first - std::sort ( combo_sort.begin(), combo_sort.end() ); - std::reverse( combo_sort.begin(), combo_sort.end() ); - -#ifdef DEBUG_ALL - std::cout << "Multiview(): combos:\n" << combos_rho << std::endl; - std::cout << "Ranked combos:\n"; - for ( auto i = 0; i < combo_sort.size(); i++ ) { - std::cout << "("; - std::pair< double, int > combo_pair = combo_sort[ i ]; - std::cout << combo_pair.first << "," - << combo_pair.second << ") "; - } std::cout << std::endl; -#endif - - // --------------------------------------------------------------- - // Perform predictions with the top multiview embeddings - // --------------------------------------------------------------- - if ( trainLib ) { - // Reset the user specified prediction vector - param.prediction = prediction; - } - - // Get top param.MultiviewEnsemble combos - size_t nEnsemble = std::min( (int) combo_sort.size(), - param.MultiviewEnsemble ); - std::vector< std::pair< double, int > > - combo_best( combo_sort.begin(), combo_sort.begin() + nEnsemble ); - -#ifdef DEBUG_ALL - std::cout << "Multiview(): Best combos:\n"; - for ( auto i = 0; i < combo_best.size(); i++ ) { - std::pair< double, int > combo_pair = combo_best[ i ]; - std::vector< size_t > this_combo = combos[ combo_pair.second ]; - std::cout << "(" << combo_pair.first << " ["; - for ( auto j = 0; j < this_combo.size(); j++ ) { - std::cout << this_combo[j] << ","; - } std::cout << "]) "; - } std::cout << std::endl; -#endif - - // Create combos_best (vector of column numbers) from combo_best - std::vector< std::vector< size_t > > combos_best( param.MultiviewEnsemble ); - for ( auto i = 0; i < combo_best.size(); i++ ) { - std::pair< double, int > combo_pair = combo_best[ i ]; - std::vector< size_t > this_combo = combos[ combo_pair.second ]; - combos_best[ i ] = this_combo; - } - - // Results Data Frame: D columns (a combo), and rho mae rmse - DataFrame combos_rho_pred( param.MultiviewEnsemble, - D + 3, header.str() ); - - // Results vector of DataFrame's with prediction results - // Used to compute the multiview ensemble average prediction - std::vector< DataFrame< double > > - combos_rho_prediction( param.MultiviewEnsemble ); - - //-------------------------------------------------------------------- - // If trainLib false, no need to compute these projections - //-------------------------------------------------------------------- - if ( trainLib ) { - // Build work queue - EDM_Multiview::WorkQueue workQ_pred( param.MultiviewEnsemble ); - - // Insert combos index into work queue - for ( auto i = 0; i < param.MultiviewEnsemble; i++ ) { - workQ_pred[ i ] = i; - } - - // thread container - std::vector< std::thread > threads_pred; - for ( unsigned i = 0; i < nThreads; ++i ) { - threads_pred.push_back( - std::thread( EvalComboThread, - param, - workQ_pred, - combos_best, - std::ref( embedding ), - std::ref( targetVec ), - std::ref( combos_rho_pred ), - std::ref( combos_rho_prediction) ) ); - } - - // join threads - for ( auto &thrd : threads_pred ) { - thrd.join(); - } - - // If thread threw exception, get from queue and rethrow - if ( not EDM_Multiview::exceptionQ.empty() ) { - std::lock_guard lck( EDM_Multiview::q_mtx ); - - // Take the first exception in the queue - std::exception_ptr exceptionPtr = EDM_Multiview::exceptionQ.front(); - - // Unroll all other exception from the thread/loops - while( not EDM_Multiview::exceptionQ.empty() ) { - // JP When do these exception_ptr get deleted? Is it a leak? - EDM_Multiview::exceptionQ.pop(); - } - std::rethrow_exception( exceptionPtr ); - } -#ifdef DEBUG_ALL - for ( auto cpi = combos_rho_prediction.begin(); - cpi != combos_rho_prediction.end(); ++cpi ) { - std::cout << *cpi; - } - std::cout << combos_rho_pred; -#endif - } // if ( trainLib ) - //------------------------------------------------------------------- - // else: trainLib = false - // Projections were initally made with lib != pred - //------------------------------------------------------------------- - else { - // Insert top prediction results into combos_rho_pred - // Insert top predictions into combos_rho_prediction - int row; - for ( size_t row_i = 0; row_i < param.MultiviewEnsemble; row_i++ ) { - row = combo_sort[ row_i ].second; // row index of best rho - combos_rho_pred.WriteRow( row_i, combos_rho.Row( row ) ); - combos_rho_prediction[ row_i ] = combos_prediction[ row ]; - } - } - - //---------------------------------------------------------- - // Compute Multiview averaged prediction - // combos_rho_prediction is a vector of DataFrames with - // columns [ Observations, Predictions ] - //---------------------------------------------------------- - // Get copy of Observations - std::valarray< double > - Obs = combos_rho_prediction[0].VectorColumnName( "Observations" ); - - // Create ensemble average prediction vector - std::valarray< double > Predictions( 0., Obs.size() ); - - // Compute ensemble prediction (this seems clunky...) - // combos_rho_prediction is a vector of DataFrames with prediction - // from each combo - for ( auto i = 0; i < param.MultiviewEnsemble; i++ ) { - std::valarray< double > prediction_i = - combos_rho_prediction[ i ].VectorColumnName( "Predictions" ); - - // Accumulate prediction values - for ( auto j = 0; j < Predictions.size(); j++ ) { - Predictions[ j ] += prediction_i[ j ]; - } - } - // Mean of prediction values - for ( auto i = 0; i < Predictions.size(); i++ ) { - Predictions[ i ] /= param.MultiviewEnsemble; - } - - // Error of ensemble prediction - VectorError ve = ComputeError( Obs, Predictions ); - - // Output Prediction DataFrame - DataFrame< double > Prediction( Predictions.size(), 2, - "Observations Predictions" ); - // Output time vector - std::vector< std::string > predTime( param.prediction.size() + param.Tp ); - - FillTimes( param, data.Time(), std::ref( predTime ) ); - - Prediction.Time() = predTime; - Prediction.TimeName() = data.TimeName(); - Prediction.WriteColumn( 0, Obs ); - Prediction.WriteColumn( 1, Predictions ); - - if ( outputFile.size() ) { - Prediction.WriteData( param.pathOut, outputFile ); - } - - // Create combos_rho table with column names - std::vector< std::string > comboTable = - ComboRhoTable( combos_rho_pred, embedding.ColumnNames() ); - - if ( param.verbose ) { - std::cout << "Multiview(): rho " << ve.rho - << " MAE " << ve.MAE << " RMSE " << ve.RMSE << std::endl; - std::cout << std::endl << "Multiview Combinations:" << std::endl; - for ( auto tableRow : comboTable ) { - std::cout << tableRow << std::endl; - } std::cout << std::endl; - } - - struct MultiviewValues MV( combos_rho_pred, Prediction, comboTable ); - - return MV; -} - -//---------------------------------------------------------------- -// Worker thread -// Output: Write rho to combos_rho DataFrame, -// Simplex results to combos_prediction -//---------------------------------------------------------------- -void EvalComboThread( Parameters param, - EDM_Multiview::WorkQueue workQ, - std::vector< std::vector< size_t > > combos, - DataFrame< double > &embedding, - std::valarray< double > &targetVec, - DataFrame< double > &combos_rho, - std::vector< DataFrame< double > > &combos_prediction ) -{ - // atomic_fetch_add(): Adds val to the contained value and returns - // the value it had immediately before the operation. - std::size_t eval_i = - std::atomic_fetch_add( &EDM_Multiview::eval_i, std::size_t(1) ); - - while( eval_i < workQ.size() ) { - - // WorkQueue stores combo index in combos - size_t combo_i = workQ[ eval_i ]; - - // Get the combo for this thread - std::vector< size_t > combo = combos[ combo_i ]; - - try { - - // Local copy with combo column indices (zero-offset) - std::vector< size_t > combo_cols( combo ); - - // Zero offset combo column indices for dataFrame - for ( auto ci = combo_cols.begin(); ci != combo_cols.end(); ++ci ) { - *ci = *ci - 1; - } - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_Multiview::mtx ); - std::cout << "EvalComboThread() Thread [" - << std::this_thread::get_id() << "] "; - std::cout << "combo: ["; - for ( auto i = 0; i < combo.size(); i++ ) { - std::cout << combo[i] << ","; - } std::cout << "] rho = "; - } -#endif - - // Select combo columns from the data - DataFrame comboData = - embedding.DataFrameFromColumnIndex( combo_cols ); - - // Compute neighbors on comboData - Neighbors neighbors = FindNeighbors( comboData, param ); - - // Pack embedding, target, neighbors for SimplexProjection - DataEmbedNN embedNN = DataEmbedNN( &embedding, comboData, - targetVec, neighbors ); - - // combo prediction - // This is an embedded = true, E = D columns prediction - DataFrame S = SimplexProjection( param, embedNN ); - - // Write combo prediction DataFrame - combos_prediction[ eval_i ] = S; - - // Evaluate combo prediction - VectorError ve = ComputeError( S.VectorColumnName( "Observations" ), - S.VectorColumnName( "Predictions" ) ); - -#ifdef DEBUG_ALL - { - std::lock_guard lck( EDM_Multiview::mtx ); - std::cout << ve.rho << std::endl; - std::cout << "-------------- embedding -------------------\n"; - std::cout << embedding; - std::cout << "-------------- comboData -------------------\n"; - std::cout << comboData; - } -#endif - - // Write combo and rho to the Data Frame - // D columns (a combo), rho, MAE, RMSE - std::valarray< double > combo_row( combo.size() + 3 ); - for ( auto i = 0; i < combo.size(); i++ ) { - combo_row[ i ] = combo[ i ]; - } - combo_row[ combo.size() ] = ve.rho; - combo_row[ combo.size() + 1 ] = ve.MAE; - combo_row[ combo.size() + 2 ] = ve.RMSE; - - combos_rho.WriteRow( eval_i, combo_row ); - - } // try - catch(...) { - // push exception pointer onto queue for main thread to catch - std::lock_guard lck( EDM_Multiview::q_mtx ); - EDM_Multiview::exceptionQ.push( std::current_exception() ); - } - - eval_i = std::atomic_fetch_add(&EDM_Multiview::eval_i, std::size_t(1)); - } - - // Reset counter - std::atomic_store( &EDM_Multiview::eval_i, std::size_t(0) ); -} - -//---------------------------------------------------------------- -// Return combinations C(n,k) as a vector of vectors -//---------------------------------------------------------------- -std::vector< std::vector< size_t > > Combination( int n, int k ) { - - std::vector< bool > v(n); - for ( size_t i = 0; i < n; ++i ) { - v[i] = ( i >= (n - k) ); - } - - std::vector< std::vector< size_t > > combos; - - do { - std::vector< size_t > this_combo( k ); - size_t j = 0; - for ( size_t i = 0; i < n; ++i ) { - if ( v[i] ) { - this_combo[ j ] = i + 1; - j++; - } - } - // insert this tuple in the combos vector - combos.push_back( this_combo ); - - } while ( std::next_permutation( v.begin(), v.end() ) ); - - return combos; -} - -//---------------------------------------------------------------- -// Return combos_rho_pred DataFrame as a vector of strings -// with column names. -//---------------------------------------------------------------- -std::vector< std::string > ComboRhoTable( - DataFrame combos_rho_pred, - std::vector< std::string > columnNames ) -{ - - // combos_rho_pred has E + 3 columns: Col_1, ... Col_E, rho, MAE, RMSE - size_t nCol = combos_rho_pred.NColumns() - 3; // JP Hardcoded silliness! - - if ( nCol > columnNames.size() ) { - std::stringstream errMsg; - errMsg << "ComboRhoTable(): Combos_rho has " << nCol - << " columns, but the data embedding has " - << columnNames.size() << " elements."; - throw std::runtime_error( errMsg.str() ); - } - - std::vector< std::string > table; - - // Header - std::stringstream header; - for ( size_t col = 0; col < nCol; col++ ) { // column indices - header << "col_" << col + 1 << ", "; - } - for ( size_t col = 0; col < nCol; col++ ) { // column names - header << "name_" << col + 1 << ", "; - } - header << "rho, MAE, RMSE"; - table.push_back( header.str() ); - - // Process each row of combos_rho_pred - for ( size_t row = 0; row < combos_rho_pred.NRows(); row++ ) { - std::stringstream rowsstring; - rowsstring.precision( 4 ); - - std::valarray< double > rowValues = combos_rho_pred.Row( row ); - - for ( size_t col = 0; col < nCol; col++ ) { - rowsstring << std::setw(4) << rowValues[ col ] << ", "; - } - for ( size_t col = 0; col < nCol; col++ ) { - size_t col_i = (size_t) rowValues[ col ]; - rowsstring << columnNames[ col_i - 1 ] << ", "; - } - - rowsstring << std::setw(6) << rowValues[ nCol ] << ", "; // rho - rowsstring << std::setw(6) << rowValues[ nCol + 1 ] << ", "; // MAE - rowsstring << std::setw(6) << rowValues[ nCol + 2 ]; // RMSE - - table.push_back( rowsstring.str() ); - } - - return table; -} diff --git a/src/Neighbors.cc b/src/Neighbors.cc deleted file mode 100644 index 87595cb..0000000 --- a/src/Neighbors.cc +++ /dev/null @@ -1,277 +0,0 @@ - -#include "Neighbors.h" - -//---------------------------------------------------------------- -Neighbors:: Neighbors() {} -Neighbors::~Neighbors() {} - -namespace EDM_Neighbors { - // Define the initial maximum distance for neigbors to avoid sort() - // DBL_MAX is a Macro equivalent to: std::numeric_limits::max() - double DistanceMax = std::numeric_limits::max(); - double DistanceLimit = std::numeric_limits::max() - 1; -} - -//---------------------------------------------------------------- -// It is assumed that the data frame has only columns of data for -// which knn will be computed. The (time) column is not present. -//---------------------------------------------------------------- -Neighbors FindNeighbors( - DataFrame dataFrame, - Parameters parameters ) -{ - -#ifdef DEBUG_ALL - PrintDataFrameIn( dataFrame, parameters ); -#endif - - if ( not parameters.validated ) { - std::string errMsg("FindNeighbors(): Parameters not validated." ); - throw( std::runtime_error( errMsg ) ); - } - - if ( parameters.embedded and parameters.E > dataFrame.NColumns() ) { - std::stringstream errMsg; - errMsg << "WARNING: FindNeighbors() Multivariate data " - << "(embedded = true): The number of dataFrame columns (" - << dataFrame.NColumns() << ") is less than the embedding " - << "dimension E (" << parameters.E << ")\n"; - std::cout << errMsg.str(); - } - - size_t N_library_rows = parameters.library.size(); - size_t N_prediction_rows = parameters.prediction.size(); - size_t N_columns = dataFrame.NColumns(); - - auto max_lib_it = std::max_element( parameters.library.begin(), - parameters.library.end() ); - size_t max_lib_index = *max_lib_it; - - // Maximum column index. - // We assume that the dataFrame has been selected to the proper columns - size_t maxCol_i = N_columns - 1; - - if ( parameters.verbose ) { - // Identify degenerate library : prediction points by - // set_intersection() of lib & pred indices, needs a result vector - std::vector< double > result( N_library_rows + N_prediction_rows, 0 ); - - std::vector< double >::iterator ii = set_intersection ( - parameters.prediction.begin(), parameters.prediction.end(), - parameters.library.begin(), parameters.library.end(), - result.begin() ); - - if ( ii != result.begin() ) { - // Overlapping indices exist - std::stringstream msg; - msg << "WARNING: FindNeighbors(): Degenerate library and " - << " prediction data found. Overlap indices: "; - for ( auto ri = result.begin(); ri != ii; ++ri ) { - msg << *ri << " "; - } msg << std::endl; - std::cout << msg.str(); - } - } - - // Neighbors: struct on local stack to be returned by copy - Neighbors neighbors = Neighbors(); - neighbors.neighbors = DataFrame(N_prediction_rows, parameters.knn); - neighbors.distances = DataFrame(N_prediction_rows, parameters.knn); - - // Vectors to hold indices and values from each comparison - std::valarray k_NN_neighbors( parameters.knn ); - std::valarray k_NN_distances( parameters.knn ); - - //------------------------------------------------------------------- - // For each prediction vector (row in prediction DataFrame) find the - // list of library indices that are within k_NN points - //------------------------------------------------------------------- - for ( size_t row_i = 0; row_i < parameters.prediction.size(); row_i++ ) { - // Get the prediction vector for this pred_row index - size_t pred_row = parameters.prediction[ row_i ]; - std::valarray pred_vec = dataFrame.Row( pred_row ); - - // Reset the neighbor and distance vectors for this pred row - for ( size_t i = 0; i < parameters.knn; i++ ) { - k_NN_neighbors[ i ] = 0; - // JP: Used to avoid sort() - k_NN_distances[ i ] = EDM_Neighbors::DistanceMax; - } - - //-------------------------------------------------------------- - // Library Rows - //-------------------------------------------------------------- - for ( size_t row_j = 0; row_j < parameters.library.size(); row_j++ ) { - // Get the library vector for this lib_row index - size_t lib_row = parameters.library[ row_j ]; - std::valarray lib_vec = dataFrame.Row( lib_row ); - - // If the library point is degenerate with the prediction, - // ignore it. - if ( lib_row == pred_row ) { -#ifdef DEBUG_ALL - if ( parameters.verbose ) { - std::stringstream msg; - msg << "FindNeighbors(): Ignoring degenerate lib_row " - << lib_row << " and pred_row " << pred_row << std::endl; - std::cout << msg.str(); - } -#endif - continue; - } - - // Apply temporal exclusion radius: units are data rows, not time - if ( parameters.exclusionRadius ) { - int xrad = (int) lib_row - pred_row; - if ( std::abs( xrad ) <= parameters.exclusionRadius ) { - continue; - } - } - - // If this lib_row + args.Tp >= max_lib_index, then this neighbor - // would be outside the library, keep looking if noNeighborLimit - if ( not parameters.noNeighborLimit ) { - if ( lib_row + parameters.Tp > max_lib_index ) { - continue; - } - } - - // Find distance between the prediction vector - // and each of the library vectors - // The 1st column (j=0) of Time has been excluded above - double d_i = Distance( lib_vec, pred_vec, - DistanceMetric::Euclidean ); - - // If d_i is less than values in k_NN_distances, add to list - auto max_it = std::max_element( begin( k_NN_distances ), - end( k_NN_distances ) ); - if ( d_i < *max_it ) { - size_t max_i = std::distance( begin(k_NN_distances), max_it ); - k_NN_neighbors[ max_i ] = lib_row; // Save the index - k_NN_distances[ max_i ] = d_i; // Save the value - } - } // for ( row_j = 0; row_j < library.size(); row_j++ ) - - if ( *std::max_element( begin( k_NN_distances ), - end ( k_NN_distances ) ) > - EDM_Neighbors::DistanceLimit ) { - - std::stringstream errMsg; - errMsg << "FindNeighbors(): Failed to find " - << parameters.knn << " knn neighbors. The library " - << "may be too small." << std::endl; - throw std::runtime_error( errMsg.str() ); - } - - // Check for ties. JP: Need to address this, not just warning - // First sort a copy of k_NN_neighbors so unique() will work - std::valarray k_NN_neighborCopy( k_NN_neighbors ); - std::sort( begin( k_NN_neighborCopy ), end( k_NN_neighborCopy ) ); - - // ui is iterator to first non unique element - auto ui = std::unique( begin( k_NN_neighborCopy ), - end ( k_NN_neighborCopy ) ); - - if ( std::distance( begin( k_NN_neighborCopy ), ui ) != - k_NN_neighborCopy.size() ) { - std::cout << "WARNING: FindNeighbors(): Degenerate neighbors./n"; - } - - // Write the neighbor indices and distance values - neighbors.neighbors.WriteRow( row_i, k_NN_neighbors ); - neighbors.distances.WriteRow( row_i, k_NN_distances ); - - } // for ( row_i = 0; row_i < predictionRows->size(); row_i++ ) - -#ifdef DEBUG_ALL - const Neighbors &neigh = neighbors; - PrintNeighborsOut( neigh ); -#endif - - return neighbors; -} - -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -double Distance( const std::valarray &v1, - const std::valarray &v2, - DistanceMetric metric ) -{ - double distance = 0; - - // For efficiency sake, we forego the usual validation of v1 & v2. - - if ( metric == DistanceMetric::Euclidean ) { - double sum = 0; - double delta = 0; - for ( size_t i = 0; i < v1.size(); i++ ) { - delta = v2[i] - v1[i]; - sum += delta * delta; // avoid call to pow() - } - distance = sqrt( sum ); - - // Note: this implicit implementation is slower - // std::valarray delta = v2 - v1; - // distance = sqrt( (delta * delta).sum() ); - } - else if ( metric == DistanceMetric::Manhattan ) { - double sum = 0; - for ( size_t i = 0; i < v1.size(); i++ ) { - sum += abs( v2[i] - v1[i] ); - } - distance = sum; - } - else { - std::stringstream errMsg; - errMsg << "Distance() Invalid DistanceMetric: " - << static_cast( metric ); - throw std::runtime_error( errMsg.str() ); - } - - return distance; -} - -#ifdef DEBUG_ALL -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -void PrintDataFrameIn( const DataFrame &dataFrame, - const Parameters ¶meters ) -{ - std::cout << "FindNeighbors(): library:" << std::endl; - for ( size_t row = 0; row < parameters.library.size(); row++ ) { - size_t row_i = parameters.library[row]; - std::cout << "row " << row_i << " : "; - for ( size_t col = 0; col < dataFrame.NColumns(); col++ ) { - std::cout << dataFrame(row_i,col) << " "; - } std::cout << std::endl; - } - std::cout << "FindNeighbors(): prediction:" << std::endl; - for ( size_t row = 0; row < parameters.prediction.size(); row++ ) { - size_t row_i = parameters.prediction[row]; - std::cout << "row " << row_i << " : "; - for ( size_t col = 0; col < dataFrame.NColumns(); col++ ) { - std::cout << dataFrame(row_i,col) << " "; - } std::cout << std::endl; - } -} - -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -void PrintNeighborsOut( const Neighbors &neighbors ) -{ - std::cout << "FindNeighbors(): neighbors:distances" << std::endl; - //for ( size_t i = 0; i < neighbors.neighbors.NRows(); i++ ) { - for ( size_t i = 0; i < 5; i++ ) { - std::cout << "Row " << i << " | "; - for ( size_t j = 0; j < neighbors.neighbors.NColumns(); j++ ) { - std::cout << neighbors.neighbors( i, j ) << " "; - } std::cout << " : "; - for ( size_t j = 0; j < neighbors.neighbors.NColumns(); j++ ) { - std::cout << neighbors.distances( i, j ) << " "; - } std::cout << std::endl; - } -} -#endif diff --git a/src/Neighbors.h b/src/Neighbors.h deleted file mode 100644 index 8655c7b..0000000 --- a/src/Neighbors.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef NEIGHBORS_H -#define NEIGHBORS_H - -#include -#include - -#include "Common.h" -#include "Parameter.h" - -// Return structure of FindNeighbors() -struct Neighbors { - DataFrame neighbors; - DataFrame distances; - Neighbors(); - ~Neighbors(); -}; - -// Prototypes -Neighbors FindNeighbors( DataFrame dataFrame, - Parameters parameters ); - -void PrintDataFrameIn( const DataFrame &dataFrame, - const Parameters ¶meters ); - -void PrintNeighborsOut( const Neighbors &neighbors ); - -double Distance( const std::valarray &v1, - const std::valarray &v2, - DistanceMetric metric ); - -#endif diff --git a/src/Parameter.cc b/src/Parameter.cc deleted file mode 100644 index 4b1b30c..0000000 --- a/src/Parameter.cc +++ /dev/null @@ -1,625 +0,0 @@ - -#include "Parameter.h" - -//---------------------------------------------------------------- -// Constructor -// Default values set in Parameter.h -//---------------------------------------------------------------- -Parameters::Parameters( - Method method, - std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib_str, - std::string pred_str, - int E, - int Tp, - int knn, - int tau, - double theta, - int exclusionRadius, - - std::string columns_str, - std::string target_str, - - bool embedded, - bool const_predict, - bool verbose, - - std::string SmapFile, - std::string blockFile, - std::string derivatives_str, - - double svdSig, - double tikhonov, - double elasticNet, - - int multi, - std::string libSizes_str, - int sample, - bool random, - bool replacement, - unsigned rseed, - bool noNeigh - ) : - // Variable initialization from Parameters arguments - method ( method ), - pathIn ( pathIn ), - dataFile ( dataFile ), - pathOut ( pathOut ), - predictOutputFile( predictFile ), - lib_str ( lib_str ), - pred_str ( pred_str ), - E ( E ), - Tp ( Tp ), - knn ( knn ), - tau ( tau ), - theta ( theta ), - exclusionRadius ( exclusionRadius ), - - columns_str ( columns_str ), - target_str ( target_str ), - targetIndex ( 0 ), - - embedded ( embedded ), - const_predict ( const_predict ), - verbose ( verbose ), - - SmapOutputFile ( SmapFile ), - blockOutputFile ( blockFile ), - - derivatives_str ( derivatives_str ), - SVDSignificance ( svdSig ), - TikhonovAlpha ( tikhonov ), - ElasticNetAlpha ( elasticNet ), - - MultiviewEnsemble( multi ), - libSizes_str ( libSizes_str ), - subSamples ( sample ), - randomLib ( random ), - replacement ( replacement ), - seed ( rseed ), - noNeighborLimit ( noNeigh ), - - // Set validated flag and instantiate Version - validated ( false ), - version ( 1, 3, 5, "2020-04-16" ) -{ - // Constructor code - if ( method != Method::None ) { - - Validate(); - - if ( verbose ) { - version.ShowVersion(); - } - } -} - -//---------------------------------------------------------------- -// Destructor -//---------------------------------------------------------------- -Parameters::~Parameters() {} - -//---------------------------------------------------------------- -// -//---------------------------------------------------------------- -void Parameters::Load() {} - -//---------------------------------------------------------------- -// Index offsets, generate library and prediction indices, -// and parameter validation -//---------------------------------------------------------------- -void Parameters::Validate() { - - validated = true; - - if ( not embedded and tau == 0 ) { - std::string errMsg( "Parameters::Validate(): " - "tau must be non-zero.\n" ); - throw std::runtime_error( errMsg ); - } - - if ( Tp < 0 ) { - std::string errMsg( "Parameters::Validate(): " - "Tp must be positive.\n" ); - throw std::runtime_error( errMsg ); - } - - //-------------------------------------------------------------- - // Generate library indices: Apply zero-offset - //-------------------------------------------------------------- - if ( lib_str.size() ) { - // Parse lib_str into vector of strings - std::vector lib_vec = SplitString( lib_str, " \t," ); - if ( lib_vec.size() % 2 != 0 ) { - std::string errMsg( "Parameters::Validate(): " - "library must be even number of integers.\n" ); - throw std::runtime_error( errMsg ); - } - - // Generate vector of start, stop index pairs - std::vector< std::pair< size_t, size_t > > libPairs; - for ( size_t i = 0; i < lib_vec.size(); i = i + 2 ) { - libPairs.emplace_back( std::make_pair( std::stoi( lib_vec[i] ), - std::stoi( lib_vec[i+1] ) ) ); - } - - size_t nLib = 0; // Count of lib items - - // Get number of lib indices, validate end > start - for ( auto thisPair : libPairs ) { - size_t lib_start = thisPair.first; - size_t lib_end = thisPair.second; - - nLib += lib_end - lib_start + 1; - - // Validate end > stop indices - if ( method == Method::Simplex or method == Method::SMap ) { - // Don't check if method == None, Embed or CCM since default - // of "1 1" is used. - if ( lib_start >= lib_end ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): library start " - << lib_start << " exceeds end " << lib_end << ".\n"; - throw std::runtime_error( errMsg.str() ); - } - } - } - - // Create library vector of indices - library = std::vector< size_t >( nLib ); - size_t i = 0; - for ( auto thisPair : libPairs ) { - for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { - library[ i ] = li - 1; // apply zero-offset - i++; - } - } - } - - //-------------------------------------------------------------- - // Generate prediction indices: Apply zero-offset - //-------------------------------------------------------------- - if ( pred_str.size() ) { - // Parse pred_str into vector of strings - std::vector pred_vec = SplitString( pred_str, " \t," ); - if ( pred_vec.size() % 2 != 0 ) { - std::string errMsg( "Parameters::Validate(): " - "prediction must be even number of integers.\n"); - throw std::runtime_error( errMsg ); - } - - // Generate vector of start, stop index pairs - std::vector< std::pair< size_t, size_t > > predPairs; - for ( size_t i = 0; i < pred_vec.size(); i = i + 2 ) { - predPairs.emplace_back( std::make_pair( std::stoi( pred_vec[i] ), - std::stoi( pred_vec[i+1]))); - } - - size_t nPred = 0; // Count of pred items - - // Get number of pred indices, validate end > start - for ( auto thisPair : predPairs ) { - size_t pred_start = thisPair.first; - size_t pred_end = thisPair.second; - - nPred += pred_end - pred_start + 1; - - // Validate end > stop indices - if ( method == Method::Simplex or method == Method::SMap ) { - // Don't check if method == None, Embed or CCM since default - // of "1 1" is used. - if ( pred_start >= pred_end ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): prediction start " - << pred_start << " exceeds end " << pred_end << ".\n"; - throw std::runtime_error( errMsg.str() ); - } - } - } - - // Create prediction vector of indices - prediction = std::vector< size_t >( nPred ); - size_t i = 0; - for ( auto thisPair : predPairs ) { - for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { - prediction[ i ] = li - 1; // apply zero-offset - i++; - } - } - } - - if ( method == Method::Simplex or method == Method::SMap ) { - if ( not library.size() ) { - std::string errMsg( "Parameters::Validate(): " - "library indices not found.\n" ); - throw std::runtime_error( errMsg ); - } - if ( not prediction.size() ) { - std::string errMsg( "Parameters::Validate(): " - "prediction indices not found.\n" ); - throw std::runtime_error( errMsg ); - } - } - else { - // Defaults if Method is None, Embed or CCM - if ( not library.size() ) { - library = std::vector( 1, 0 ); - } - if ( not prediction.size() ) { - prediction = std::vector( 1, 0 ); - } - } - -#ifdef DEBUG_ALL - PrintIndices( library, prediction ); -#endif - - //-------------------------------------------------------------- - // Convert multi argument parameters from string to vectors - //-------------------------------------------------------------- - - // Columns - // If columns are purely integer, then populate vector columnIndex - // Otherwise fill in vector columnNames - if ( columns_str.size() ) { - - std::vector columns_vec = SplitString( columns_str, - " \t,\n" ); - - bool onlyDigits = false; - - for ( auto ci = columns_vec.begin(); ci != columns_vec.end(); ++ci ) { - onlyDigits = OnlyDigits( *ci, true ); - - if ( not onlyDigits ) { break; } - } - - if ( onlyDigits ) { - for ( auto ci = columns_vec.begin(); - ci != columns_vec.end(); ++ci ) { - columnIndex.push_back( std::stoi( *ci ) ); - } - } - else { - columnNames = columns_vec; - } - } - - if ( not columnIndex.size() and not columnNames.size() ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): Simplex/CCM: " - << " No valid columns found." << std::endl; - throw std::runtime_error( errMsg.str() ); - } - - // target - if ( target_str.size() ) { - bool onlyDigits = OnlyDigits( target_str, true ); - if ( onlyDigits ) { - targetIndex = std::stoi( target_str ); - } - else { - targetName = target_str; - } - } - - // Derivatives - if ( derivatives_str.size() > 0 ) { - std::vector der_vec = SplitString(derivatives_str," \t,"); - if ( der_vec.size() < 2 ) { - std::string errMsg( "Parameters::Validate(): " - "derivatives must be integer pairs.\n"); - throw std::runtime_error( errMsg ); - } - derivatives = std::vector( der_vec.size() ); - for ( size_t i = 0; i < der_vec.size(); i++ ) { - derivatives.push_back( std::stoi( der_vec[i] ) ); - } - } - - //-------------------------------------------------------------------- - // CCM sample not 0 if random is true - if ( method == Method::CCM ) { - if ( randomLib ) { - if ( subSamples < 1 ) { - std::string errMsg( "Parameters::Validate(): " - "CCM samples must be > 0.\n" ); - throw std::runtime_error( errMsg ); - } - } - } - - // CCM librarySizes - if ( libSizes_str.size() > 0 ) { - std::vector libsize_vec = SplitString(libSizes_str," \t,"); - if ( libsize_vec.size() != 3 ) { - std::string errMsg( "Parameters::Validate(): " - "CCM librarySizes must be three integers.\n" ); - throw std::runtime_error( errMsg ); - } - - size_t start = std::stoi( libsize_vec[0] ); - size_t stop = std::stoi( libsize_vec[1] ); - size_t increment = std::stoi( libsize_vec[2] ); - - if ( increment < 1 ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): " - << "CCM librarySizes increment " << increment - << " is invalid.\n"; - throw std::runtime_error( errMsg.str() ); - } - - if ( start > stop ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): " - << "CCM librarySizes start " << start - << " stop " << stop << " are invalid.\n"; - throw std::runtime_error( errMsg.str() ); - } - - size_t N_lib = std::floor( (stop-start)/increment + 1/increment ) + 1; - - if ( start < E ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): " - << "CCM librarySizes start < E = " << E << "\n"; - throw std::runtime_error( errMsg.str() ); - } - else if ( start < 3 ) { - std::string errMsg( "Parameters::Validate(): " - "CCM librarySizes start < 3.\n" ); - throw std::runtime_error( errMsg ); - } - - // Create the librarySizes vector - librarySizes = std::vector( N_lib, 0 ); - - // Fill in the sizes - size_t libSize = start; - for ( size_t i = 0; i < librarySizes.size(); i++ ) { - librarySizes[i] = libSize; - libSize = libSize + increment; - } - } - - //-------------------------------------------------------------------- - // Simplex and knn not specified: not embedded : knn set to E+1 - // embedded : knn set to size( columns ) - // S-Map require knn > E + 1, default is all neighbors. - if ( method == Method::Simplex or method == Method::CCM ) { - if ( knn < 1 ) { - if ( not embedded ) { - knn = E + 1; - if ( verbose ) { - std::stringstream msg; - msg << "Parameters::Validate(): Set knn = " << knn - << " (E+1) for Simplex. " << std::endl; - std::cout << msg.str(); - } - } - else { // embedded = true - if ( columnIndex.size() ) { - knn = columnIndex.size() + 1; - } - else if ( columnNames.size() ) { - knn = columnNames.size() + 1; - } - - if ( verbose ) { - std::stringstream msg; - msg << "Parameters::Validate(): Set knn = " << knn - << " for Simplex (embedded = true). " << std::endl; - std::cout << msg.str(); - } - } - } - if ( knn < E + 1 ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate(): Simplex knn of " << knn - << " is less than E+1 = " << E + 1 << std::endl; - throw std::runtime_error( errMsg.str() ); - } - } - else if ( method == Method::SMap ) { - if ( knn > 0 ) { - if ( knn < E + 1 ) { - std::stringstream errMsg; - errMsg << "Parameters::Validate() S-Map knn must be at least " - " E+1 = " << E + 1 << ".\n"; - throw std::runtime_error( errMsg.str() ); - } - } - else { - // default knn = 0, set knn value - knn = library.size() - Tp * (E + 1); - if ( verbose ) { - std::stringstream msg; - msg << "Parameters::Validate(): Set knn = " << knn - << " for SMap. " << std::endl; - std::cout << msg.str(); - } - } - if ( not embedded and columnNames.size() > 1 ) { - std::string msg( "Parameters::Validate() WARNING: " - "Multivariable S-Map should use " - "-e (embedded) data input to ensure " - "data/dimension correspondance.\n" ); - std::cout << msg; - } - - // S-Map coefficient columns for derivatives start at 1 since the 0th - // column is the S-Map linear prediction bias term - if ( derivatives.size() > 1 ) { - std::vector::iterator it = std::find( derivatives.begin(), - derivatives.end(), 0); - if ( it != derivatives.end() ) { - std::string errMsg( "Parameters::Validate() S-Map coefficient " - " columns for derivatives can not use column 0.\n"); - throw std::runtime_error( errMsg ); - } - if ( derivatives.size() % 2 ) { - std::string errMsg( "Parameters::Validate() S-Map coefficient " - " columns for derivatives must be in pairs.\n"); - throw std::runtime_error( errMsg ); - } - } - - // Tikhonov and ElasticNet are mutually exclusive - if ( TikhonovAlpha and ElasticNetAlpha ) { - std::string errMsg( "Parameters::Validate() Multiple S-Map solve " - "methods specified. Use one or none of: " - "tikhonov, elasticNet.\n"); - throw std::runtime_error( errMsg ); - } - - // Very small alphas don't make sense in elastic net - if ( ElasticNetAlpha < 0.01 ) { - std::cout << "Parameters::Validate() ElasticNetAlpha too small." - " Setting to 0.01."; - ElasticNetAlpha = 0.01; - } - if ( ElasticNetAlpha > 1 ) { - std::cout << "Parameters::Validate() ElasticNetAlpha too large." - " Setting to 1."; - ElasticNetAlpha = 1; - } - } - else if ( method == Method::Embed ) { - // no-op - } - else { - throw std::runtime_error( "Parameters::Validate() " - "Prediction method error.\n" ); - } -} - -//------------------------------------------------------------ -// Adjust lib/pred concordant with Embed() removal of tau(E-1) -// rows, and DeletePartialDataRow() -//------------------------------------------------------------ -void Parameters::DeleteLibPred( size_t shift ) { - - size_t library_len = library.size(); - size_t prediction_len = prediction.size(); - - // If [0, 1, ... shift] (negative tau) or - // [N-shift, ... N-1, N] (positive tau) are in library or prediction - // those rows were deleted, delete these index elements. - // First, create vectors of indices to delete. - std::vector< size_t > deleted_pred_elements( shift, 0 ); - std::vector< size_t > deleted_lib_elements ( shift, 0 ); - - if ( tau < 0 ) { - std::iota(deleted_pred_elements.begin(), deleted_pred_elements.end(),0); - std::iota(deleted_lib_elements.begin(), deleted_lib_elements.end(), 0); - } - else { - std::iota( deleted_pred_elements.begin(), - deleted_pred_elements.end(), prediction_len - shift ); - std::iota( deleted_lib_elements.begin(), - deleted_lib_elements.end(), library_len - shift ); - } - - // Erase elements of row indices that were deleted - for ( auto element = deleted_lib_elements.begin(); - element != deleted_lib_elements.end(); element++ ) { - - std::vector< size_t >::iterator it; - it = std::find( library.begin(), library.end(), *element ); - - if ( it != library.end() ) { - library.erase( it ); - } - } - - for ( auto element = deleted_pred_elements.begin(); - element != deleted_pred_elements.end(); element++ ) { - - std::vector< size_t >::iterator it; - it = std::find( prediction.begin(), - prediction.end(), *element ); - - if ( it != prediction.end() ) { - prediction.erase( it ); - } - } - - // Now offset all values by shift so that vectors indices - // in library and prediction refer to the same data rows - // before the deletion/shift. - if ( tau < 0 ) { - for ( auto li = library.begin(); li != library.end(); li++ ) { - *li = *li - shift; - } - for ( auto pi = prediction.begin(); pi != prediction.end(); pi++ ) { - *pi = *pi - shift; - } - } - // tau > 0 : Forward shifting: no adjustment needed from origin -} - -//------------------------------------------------------------------ -// Overload << to output to ostream -//------------------------------------------------------------------ -std::ostream& operator<< ( std::ostream &os, Parameters &p ) { - - // print info about the dataframe - os << "Parameters: -------------------------------------------\n"; - - std::string method("Unknown"); - if ( p.method == Method::Simplex ) { method = "Simplex"; } - else if ( p.method == Method::SMap ) { method = "SMap"; } - else if ( p.method == Method::CCM ) { method = "CCM"; } - else if ( p.method == Method::None ) { method = "None"; } - else if ( p.method == Method::Embed ) { method = "Embed"; } - - os << "Method: " << method - << " E=" << p.E << " Tp=" << p.Tp - << " knn=" << p.knn << " tau=" << p.tau << " theta=" << p.theta - << std::endl; - - if ( p.columnNames.size() ) { - os << "Column Names : [ "; - for ( auto ci = p.columnNames.begin(); - ci != p.columnNames.end(); ++ci ) { - os << *ci << " "; - } os << "]" << std::endl; - } - - if ( p.targetName.size() ) { - os << "Target: " << p.targetName << std::endl; - } - - os << "Library: [" << p.library[0] << " : " - << p.library[ p.library.size() - 1 ] << "] " - << "Prediction: [" << p.prediction[0] << " : " - << p.prediction[ p.prediction.size() - 1 ] - << "] " << std::endl; - - - os << "-------------------------------------------------------\n"; - - return os; -} - -#ifdef DEBUG_ALL -//------------------------------------------------------------------ -// -//------------------------------------------------------------------ -void Parameters::PrintIndices( std::vector library, - std::vector prediction ) -{ - std::cout << "Parameters(): library: "; - for ( auto li = library.begin(); li != library.end(); ++li ) { - std::cout << *li << " "; - } std::cout << std::endl; - std::cout << "Parameters(): prediction: "; - for ( auto pi = prediction.begin(); pi != prediction.end(); ++pi ) { - std::cout << *pi << " "; - } std::cout << std::endl; -} -#endif diff --git a/src/Parameter.h b/src/Parameter.h deleted file mode 100644 index 60125c3..0000000 --- a/src/Parameter.h +++ /dev/null @@ -1,122 +0,0 @@ -#ifndef PARAMETER_H -#define PARAMETER_H - -#include -#include - -#include "Common.h" -#include "Version.h" - -//------------------------------------------------------------ -// -//------------------------------------------------------------ -class Parameters { - -public: // Not protected with accessors. - Method method; // Simplex or SMap enum class - - std::string pathIn; // path for input dataFile - std::string dataFile; // input dataFile (assumed .csv) - std::string pathOut; // path for output files - std::string predictOutputFile;// - std::string lib_str; // multi argument parameters for library - std::string pred_str; // multi argument parameters for prediction - - std::vector library; // library row indices - std::vector prediction; // prediction row indices - - int E; // dimension - int Tp; // prediction interval - int knn; // k nearest neighbors - int tau; // block embedding delay - double theta; // S Map localization - int exclusionRadius; // temporal rows to ignore in predict - - std::string columns_str; - std::string target_str; - std::vector columnNames; // column names - std::vector columnIndex; // column indices - - std::string targetName; // target column name - size_t targetIndex; // target column index - - bool embedded; // true if data is already embedded/block - bool const_predict; // true to compute non "predictor" stats - bool verbose; - - std::string SmapOutputFile; // - std::string blockOutputFile; // Embed() output file - - std::string derivatives_str; - std::vector derivatives;// list of column indices for derivatives - - double SVDSignificance; // SVD singular value cutoff - double TikhonovAlpha; // Initial alpha parameter - double ElasticNetAlpha; // Initial alpha parameter - - int MultiviewEnsemble; // Number of ensembles in multiview - - std::string libSizes_str; - std::vector librarySizes;// CCM library sizes to evaluate - int subSamples; // CCM number of samples to draw - bool randomLib; // CCM randomly select subsets if true - bool replacement; // CCM random select with replacement if true - unsigned seed; // CCM random selection RNG seed - - bool noNeighborLimit; // Strictly forbid neighbors outside library - bool validated; - - Version version; // Version object, instantiated in constructor - - friend std::ostream& operator<<(std::ostream &os, Parameters ¶ms); - - // Constructor declaration and default arguments - Parameters( - Method method = Method::None, - std::string pathIn = "./", - std::string dataFile = "", - std::string pathOut = "./", - std::string predictFile = "", - std::string lib_str = "", - std::string pred_str = "", - int E = 0, - int Tp = 0, - int knn = 0, - int tau = -1, - double theta = 0, - int exclusionRadius = 0, - - std::string columns_str = "", - std::string target_str = "", - - bool embedded = false, - bool const_predict = false, - bool verbose = false, - - std::string SmapFile = "", - std::string blockFile = "", - std::string derivatives_str = "", - - double svdSig = 1E-5, - double tikhonov = 0, - double elasticNet = 0.1, - - int multi = 0, - std::string libSizes_str = "", - int sample = 0, - bool random = true, - bool replacement = false, - unsigned seed = 0, // 0: Generate random seed in CCM - bool noNeighbor = false - ); - - ~Parameters(); - - void Validate(); // Parameter validation and index offsets - void Load(); // Populate the parameters from arguments - void DeleteLibPred( size_t shift ); // Adjust for embedding - void PrintIndices( std::vector library, - std::vector prediction ); -}; - -#endif diff --git a/src/SMap.cc b/src/SMap.cc deleted file mode 100644 index f1c42a3..0000000 --- a/src/SMap.cc +++ /dev/null @@ -1,513 +0,0 @@ - -#include "Common.h" -#include "Parameter.h" -#include "Embed.h" -#include "Neighbors.h" -#include "AuxFunc.h" - -// forward declarations -std::valarray< double > Lapack_SVD( int m, - int n, - double *a, - double *b, - double rcond ); - -//--------------------------------------------------------------------- -// Overload 1: Explicit data file path/name with internal SVD (LAPACK) -// Implemented as a wrapper to API overload 2 -> 4 -//--------------------------------------------------------------------- -SMapValues SMap( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - double theta, - int exclusionRadius, - std::string columns, - std::string target, - std::string smapFile, - std::string derivatives, - bool embedded, - bool const_predict, - bool verbose ) -{ - // DataFrame constructor loads data - DataFrame< double > dataFrameIn( pathIn, dataFile ); - - // Call overload 2 with the dataFrameIn - SMapValues SMapOutput = SMap( dataFrameIn, pathOut, predictFile, - lib, pred, E, Tp, knn, tau, theta, - exclusionRadius, - columns, target, smapFile, derivatives, - embedded, const_predict, verbose ); - return SMapOutput; -} - -//---------------------------------------------------------------- -// Overload 2: DataFrame provided with internal SVD (LAPACK) -// Implemented as a wrapper to API overload 4 -//---------------------------------------------------------------- -SMapValues SMap( DataFrame< double > &data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - double theta, - int exclusionRadius, - std::string columns, - std::string target, - std::string smapFile, - std::string derivatives, - bool embedded, - bool const_predict, - bool verbose ) { - - // Call overload 4 with default SVD function (below) - SMapValues SMapOutput = SMap( data, pathOut, predictFile, - lib, pred, E, Tp, knn, tau, theta, - exclusionRadius, - columns, target, smapFile, derivatives, - &SVD, - embedded, const_predict, verbose); - - return SMapOutput; -} - -//---------------------------------------------------------------- -// Overload 3: Explicit data file path/name and solver -// Implemented as a wrapper to API overload 4 -//---------------------------------------------------------------- -SMapValues SMap( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - double theta, - int exclusionRadius, - std::string columns, - std::string target, - std::string smapFile, - std::string derivatives, - std::valarray (*solver) (DataFrame < double >, - std::valarray < double > ), - bool embedded, - bool const_predict, - bool verbose ) -{ - // DataFrame constructor loads data - DataFrame< double > dataFrameIn( pathIn, dataFile ); - - // Call overload 4 with dataFrameIn and solver object - SMapValues SMapOutput = SMap( dataFrameIn, pathOut, predictFile, - lib, pred, E, Tp, knn, tau, theta, - exclusionRadius, - columns, target, smapFile, derivatives, - solver, embedded, const_predict, verbose ); - return SMapOutput; -} - -//---------------------------------------------------------------- -// Overload 4: Solver & DataFrame provided -//---------------------------------------------------------------- -SMapValues SMap( DataFrame< double > &data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - double theta, - int exclusionRadius, - std::string columns, - std::string target, - std::string smapFile, - std::string derivatives, - std::valarray (*solver) (DataFrame < double >, - std::valarray < double > ), - bool embedded, - bool const_predict, - bool verbose ) -{ - - Parameters param = Parameters( Method::SMap, "", "", - pathOut, predictFile, - lib, pred, E, Tp, knn, tau, theta, - exclusionRadius, columns, target, - embedded, const_predict, verbose, - smapFile, "", derivatives ); - - //---------------------------------------------------------- - // Load data, Embed, compute Neighbors - //---------------------------------------------------------- - DataEmbedNN dataEmbedNN = EmbedNN( &data, std::ref( param ) ); - - // Unpack the dataEmbedNN for convenience - DataFrame *dataInRef = dataEmbedNN.dataIn; - DataFrame dataBlock = dataEmbedNN.dataFrame; - std::valarray target_vec = dataEmbedNN.targetVec; - Neighbors neighbors = dataEmbedNN.neighbors; - - DataFrame &dataIn = std::ref( *dataInRef ); - - //---------------------------------------------------------- - // SMap projection - //---------------------------------------------------------- - size_t library_N_row = param.library.size(); - size_t predict_N_row = param.prediction.size(); - size_t N_row = neighbors.neighbors.NRows(); - - auto max_lib_it = std::max_element( param.library.begin(), - param.library.end() ); - size_t max_lib_index = *max_lib_it; - - if ( predict_N_row != N_row ) { - std::stringstream errMsg; - errMsg << "SMap(): Number of prediction rows (" << predict_N_row - << ") does not match the number of neighbor rows (" - << N_row << ").\n"; - throw std::runtime_error( errMsg.str() ); - } - if ( neighbors.distances.NColumns() != param.knn ) { - std::stringstream errMsg; - errMsg << "SMap(): Number of neighbor columns (" - << neighbors.distances.NColumns() - << ") does not match knn (" << param.knn << ").\n"; - throw std::runtime_error( errMsg.str() ); - } - - std::valarray< double > predictions = std::valarray< double >( N_row ); - std::valarray< double > variance = std::valarray< double >( N_row ); - - // Init coefficients to NAN ? - DataFrame< double > coefficients = DataFrame< double >( N_row, - param.E + 1 ); - DataFrame< double > derivative; - DataFrame< double > tangents; - - //------------------------------------------------------------ - // Process each prediction row - //------------------------------------------------------------ - for ( size_t row = 0; row < N_row; row++ ) { - - double D_avg = neighbors.distances.Row( row ).sum() / param.knn; - - // Compute weight vector - std::valarray< double > w = std::valarray< double >( param.knn ); - if ( param.theta > 0 ) { - w = std::exp( (-param.theta/D_avg) * neighbors.distances.Row(row) ); - } - else { - w = std::valarray< double >( 1, param.knn ); - } - - DataFrame< double > A = DataFrame< double >(param.knn, param.E + 1); - std::valarray< double > B = std::valarray< double >( param.knn ); - - // Populate matrix A (exp weighted future prediction), and - // vector B (target BC's) for this row (observation). - size_t lib_row; - size_t lib_row_base; - - for ( size_t k = 0; k < param.knn; k++ ) { - lib_row_base = neighbors.neighbors( row, k ); - lib_row = lib_row_base + param.Tp; - - if ( lib_row > max_lib_index ) { - // The knn index + Tp is outside the library domain - // Can only happen if noNeighborLimit = true is used. - if ( param.verbose ) { - std::stringstream msg; - msg << "SMap() in row " << row << " libRow " << lib_row - << " exceeds library domain.\n"; - std::cout << msg.str(); - } - // Use the neighbor at the 'base' of the trajectory - B[ k ] = target_vec[ lib_row_base ]; - } - else { - B[ k ] = target_vec[ lib_row ]; - } - - //--------------------------------------------------------------- - // Linear system coefficient matrix - //--------------------------------------------------------------- - // NOTE: The matrix A has a (weighted) constant (1) first column - // to enable a linear intercept/bias term. - // NOTE: The dataBlock does not have a time vector, and only - // has columns from the embedding. So the coefficient - // matrix A has E+1 columns, while the dataBlock has E. - //--------------------------------------------------------------- - A( k, 0 ) = w[ k ]; // Intercept bias terms in column 0 (weighted) - - for ( size_t j = 1; j < param.E + 1; j++ ) { - A( k, j ) = w[k] * dataBlock( lib_row_base, j-1 ); - } - } - - B = w * B; // Weighted target vector - - // Estimate linear mapping of predictions A onto target B - std::valarray < double > C = solver( A, B ); - - // Prediction is local linear projection - double prediction = C[ 0 ]; // C[ 0 ] is the bias term - - for ( size_t e = 1; e < param.E + 1; e++ ) { - prediction = prediction + - C[ e ] * dataBlock( param.prediction[ row ], e-1 ); - } - - predictions[ row ] = prediction; - coefficients.WriteRow( row, C ); - - // "Variance" estimate assuming weights are probabilities - std::valarray< double > deltaSqr = std::pow(target_vec - predictions, 2); - variance[ row ] = ( w * deltaSqr ).sum() / w.sum(); - - } // for ( row = 0; row < predict_N_row; row++ ) - - // non "predictions" X(t+1) = X(t) if const_predict specified - std::valarray< double > const_predictions( 0., N_row ); - if ( param.const_predict ) { - std::slice pred_slice = - std::slice( param.prediction[ 0 ], param.prediction.size(), 1 ); - - const_predictions = target_vec[ pred_slice ]; - } - - //----------------------------------------------------- - // Derivatives - //----------------------------------------------------- - - - //---------------------------------------------------- - // Ouput - //---------------------------------------------------- - // Observations & predictions: Adjust rows/nan for Tp - DataFrame dataOut = FormatOutput( param, - predictions, - const_predictions, - variance, - target_vec, - dataIn.Time(), - dataIn.TimeName() ); - - // Coefficient output DataFrame: N_row + Tp rows - DataFrame< double > coefOut = DataFrame< double >( dataOut.NRows(), - param.E + 1 ); - - // Set time from: dataIn -> FormatOutput() -> FillTimes() -> dataOut - if ( dataOut.Time().size() ) { - coefOut.Time() = dataOut.Time(); - coefOut.TimeName() = dataOut.TimeName(); - } - // else { throw ? } JP - - // Populate coefOut column names: C0, C1, C2, ... - std::vector coefNames; - for ( size_t col = 0; col < coefficients.NColumns(); col++ ) { - std::stringstream coefName; - coefName << "C" << col; - coefNames.push_back( coefName.str() ); - } - coefOut.ColumnNames() = coefNames; - - // coefficients have N_row's; coefOut has N_row + Tp - // Create coefficient column vector with Tp nan rows at the - // beginning of coefOut as in FormatOutput() - std::valarray coefColumnVec( NAN, dataOut.NRows() ); - - // Copy coefficients vectors into coefOut - std::slice coef_i = std::slice( param.Tp, N_row, 1 ); - for ( size_t col = 0; col < coefOut.NColumns(); col++ ) { - coefColumnVec[ coef_i ] = coefficients.Column( col ); - coefOut.WriteColumn( col, coefColumnVec ); - } - - if ( param.predictOutputFile.size() ) { - // Write predictions to disk - dataOut.WriteData( param.pathOut, param.predictOutputFile ); - } - if ( param.SmapOutputFile.size() ) { - // Write Smap coefficients to disk - coefOut.WriteData( param.pathOut, param.SmapOutputFile ); - } - - SMapValues values = SMapValues(); - values.predictions = dataOut; - values.coefficients = coefOut; - - return values; -} - -//---------------------------------------------------------------- -// Singular Value Decomposition -//---------------------------------------------------------------- -std::valarray < double > SVD( DataFrame < double > A, - std::valarray< double > B ) { - - // NOTE: A elements are Row Major format - // Convert A to column major for LAPACK dgelss() - // a is the memory start location pointer to colMajorElements - std::valarray < double > colMajorElements = A.ColumnMajorData(); - double *a = &( colMajorElements[0] ); - - double *b = &( B[0] ); - - std::valarray < double > C = - Lapack_SVD( A.NRows(), // number of rows - A.NColumns(), // number of columns - a, // A - b, // b - 1.E-9 ); // rcond - -#ifdef DEBUG_ALL - std::cout << "SVD------------------------\n"; - std::cout << "A ----------\n"; - std::cout << A << std::endl; -#endif - - return C; -} - -//------------------------------------------------------------------------- -// subroutine dgelss() -//----------------------------------------------------------------------- -// DGELSS computes the minimum norm solution to a real linear least -// squares problem: -// -// Minimize 2-norm(| b - A*x |). -// -// using the singular value decomposition (SVD) of A. A is an M-by-N -// matrix which may be rank-deficient. -// -// Several right hand side vectors b and solution vectors x can be -// handled in a single call; they are stored as the columns of the -// M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X. -// -// The effective rank of A is determined by treating as zero those -// singular values which are less than RCOND times the largest singular -// value. -// -// INFO is INTEGER -// = 0: successful exit -// < 0: if INFO = -i, the i-th argument had an illegal value. -// > 0: the algorithm for computing the SVD failed to converge; -// if INFO = i, i off-diagonal elements of an intermediate -// bidiagonal form did not converge to zero. -//----------------------------------------------------------------------- -//------------------------------------------------------------------------- -// DOUBLE PRECISION = REAL*8 = c++ double -//------------------------------------------------------------------------- -extern "C" { - - void dgelss_( int *M, - int *N, - int *NRHS, - double *A, - int *LDA, - double *B, - int *LDB, - double *S, - double *RCOND, - int *RANK, - double *WORK, - int *LWORK, - int *INFO ); -} - -//----------------------------------------------------------------------- -// -//----------------------------------------------------------------------- -std::valarray< double > Lapack_SVD( int m, // number of rows in matrix - int n, // number of columns in matrix - double *a, // pointer to top-left corner - double *b, - double rcond ) -{ - int N_SingularValues = m < n ? m : n; - - // s to hold singular values - // MSVC BS: Have to use static const int size, or new - double *s = new double[ N_SingularValues ]; - - int lda = m; // LDA >= max(1,M) - int ldb = m; // LDB >= max(1,max(M,N)) - int nrhs = 1; - - // Workspace and info variables: - // MSVC BS: Have to use static const int size, or new - int *iwork = new int[ 8 * N_SingularValues ]; - - double workSize = 0; // To query optimal work size - int lwork = -1; // To query optimal work size - int info = 0; // return code - int rank = 0; - -#ifdef DEBUG_ALL - std::cout << "m.row=" << m << " n.col=" << n << " lda=" << lda - << " s.n=" << N_SingularValues << std::endl; - - for ( size_t i = 0; i < m*n; i++ ) { - std::cout << a[i] << " "; - } std::cout << std::endl; -#endif - - // Call dgelss with lwork = -1 to query optimal workspace size: - dgelss_( &m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, - &rank, &workSize, &lwork, &info ); - - if ( info ) { - throw std::runtime_error( "Lapack_SVD(): dgelss failed on query.\n" ); - } - -#ifdef DEBUG_ALL - std::cout << "Optimal work size is " << workSize << std::endl; -#endif - - // Optimal workspace size is returned in workSize. - // MSVC BS: Have to use static const int size, or new - double *work = new double[ (size_t) workSize ]; - - lwork = (int) workSize; - - // Call dgelss for SVD solution using lwork workSize: - dgelss_( &m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, - &rank, work, &lwork, &info ); - - if ( info ) { - throw std::runtime_error( "Lapack_SVD(): dgelss failed.\n" ); - } - -#ifdef DEBUG_ALL - std::cout << "Solution: [ "; - for ( auto i = 0; i < N_SingularValues; i++ ) { - std::cout << b[i] << " "; - } std::cout << "]" << std::endl; -#endif - - // Copy solution vector in b to C - std::valarray< double > C( b, N_SingularValues ); - - delete[] s; - delete[] work; - delete[] iwork; - - return C; -} diff --git a/src/Simplex.cc b/src/Simplex.cc index b68b0ad..abb3f91 100644 --- a/src/Simplex.cc +++ b/src/Simplex.cc @@ -1,244 +1,11 @@ +// Definitions for the Simplex class -#include "Common.h" -#include "Parameter.h" -#include "Neighbors.h" -#include "Embed.h" -#include "AuxFunc.h" +#include "Simplex.h" -// Forward declaration -DataFrame SimplexProjection( Parameters param, - DataEmbedNN embedNN, - bool checkDataRows = true ); - -//---------------------------------------------------------------- -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: //---------------------------------------------------------------- -DataFrame Simplex( std::string pathIn, - std::string dataFile, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - int exclusionRadius, - std::string columns, - std::string target, - bool embedded, - bool const_predict, - bool verbose ) { - - // DataFrame constructor loads data - DataFrame< double > *dataFrameIn = - new DataFrame< double > ( pathIn, dataFile ); - - // Pass data frame to Simplex - DataFrame< double > S = Simplex( std::ref( *dataFrameIn ), - pathOut, - predictFile, - lib, - pred, - E, - Tp, - knn, - tau, - exclusionRadius, - columns, - target, - embedded, - const_predict, - verbose ); - delete dataFrameIn; - - return S; -} - +// Simplex() : Constructor +// See EDM class for descriptions of parameters not described //---------------------------------------------------------------- -// API Overload 2: DataFrame provided -//---------------------------------------------------------------- -DataFrame Simplex( DataFrame< double > &data, - std::string pathOut, - std::string predictFile, - std::string lib, - std::string pred, - int E, - int Tp, - int knn, - int tau, - int exclusionRadius, - std::string columns, - std::string target, - bool embedded, - bool const_predict, - bool verbose ) { - - Parameters param = Parameters( Method::Simplex, "", "", - pathOut, predictFile, - lib, pred, E, Tp, knn, tau, 0, - exclusionRadius, - columns, target, embedded, - const_predict, verbose ); - - //---------------------------------------------------------- - // Embed, compute Neighbors - //---------------------------------------------------------- - DataEmbedNN embedNN = EmbedNN( &data, std::ref( param ) ); - - DataFrame S = SimplexProjection( param, embedNN ); - - return S; -} - -//---------------------------------------------------------------- -// Simplex Projection -//---------------------------------------------------------------- -DataFrame SimplexProjection( Parameters param, - DataEmbedNN embedNN, - bool checkDataRows ) { - - // Unpack the data, (embedding dataBlock not used), target & neighbors - DataFrame *dataIn = embedNN.dataIn; // used for output - std::valarray target_vec = embedNN.targetVec; - Neighbors neighbors = embedNN.neighbors; - - size_t library_N_row = param.library.size(); - size_t N_row = neighbors.neighbors.NRows(); - - auto max_lib_it = std::max_element( param.library.begin(), - param.library.end() ); - size_t max_lib_index = *max_lib_it; - -#ifdef DEBUG_ALL - std::cout << "SimplexProjection -------------------------\n"; - std::cout << "Neighbors: (" << neighbors.neighbors.NRows() << "x" - << neighbors.neighbors.NColumns() << ")\n"; - std::cout << neighbors.neighbors; - std::cout << "Target: (" << target_vec.size() << ")\n"; - for ( size_t row = 0; row < 10; row++ ) { - std::cout << target_vec[ row ] << " "; - } std::cout << std::endl; - std::cout << "-------------------------------------------\n\n"; -#endif - - if ( N_row != neighbors.distances.NRows() ) { - std::stringstream errMsg; - errMsg << "Simplex(): Number of neighbor rows " << N_row - << " doesn't match the number of distances rows " - << neighbors.distances.NRows() << std::endl; - throw std::runtime_error( errMsg.str() ); - } - - double minWeight = 1.E-6; - std::valarray predictions( 0., N_row ); - std::valarray variance ( 0., N_row ); - - // Process each prediction row in neighbors - for ( size_t row = 0; row < N_row; row++ ) { - - std::valarray distanceRow = neighbors.distances.Row( row ); - - // Establish exponential weight reference, the 'distance scale' - double minDistance = distanceRow.min(); - - // Compute weight (vector) for each k_NN - std::valarray weightedDistances( minWeight, param.knn ); - - if ( minDistance == 0 ) { - // Handle cases of distanceRow = 0 : can't divide by minDistance - for ( size_t i = 0; i < param.knn; i++ ) { - if ( distanceRow[i] > 0 ) { - weightedDistances[i] = exp( -distanceRow[i] / minDistance ); - } - else { - // Setting weight = 1 implies that the corresponding - // library target vector is the same as the observation - // so it will be given full-weight in the prediction. - weightedDistances[i] = 1; - } - } - } - else { - // exp() is a valarray<> overload (vectorized?) - weightedDistances = exp( -distanceRow / minDistance ); - } - - // weight vector - std::valarray weights( param.knn ); - for ( size_t i = 0; i < param.knn; i++ ) { - weights[i] = std::max( weightedDistances[i], minWeight ); - } - - // target library vector, one element for each knn - std::valarray libTarget( param.knn ); - - for ( size_t k = 0; k < param.knn; k++ ) { - size_t libRow = (size_t) neighbors.neighbors( row, k ) + param.Tp; - - if ( libRow > max_lib_index ) { - // The k_NN index + Tp is outside the library domain - // Can only happen if noNeighborLimit = true is used. - if ( param.verbose ) { - std::stringstream msg; - msg << "Simplex() in row " << row << " libRow " << libRow - << " exceeds library domain.\n"; - std::cout << msg.str(); - } - - // Use the neighbor at the 'base' of the trajectory - libTarget[ k ] = target_vec[ libRow - param.Tp ]; - } - else { - libTarget[ k ] = target_vec[ libRow ]; - } - } - - // Prediction is average of weighted library projections - predictions[ row ] = ( weights * libTarget ).sum() / weights.sum(); - - // "Variance" estimate assuming weights are probabilities - std::valarray< double > deltaSqr = std::pow(libTarget - predictions, 2); - variance[ row ] = ( weights * deltaSqr ).sum() / weights.sum(); - - } // for ( row = 0; row < N_row; row++ ) - - // non "predictions" X(t+1) = X(t) if const_predict specified - std::valarray< double > const_predictions( 0., N_row ); - if ( param.const_predict ) { - std::slice pred_slice = - std::slice( param.prediction[ 0 ], param.prediction.size(), 1 ); - - const_predictions = target_vec[ pred_slice ]; - } - - //---------------------------------------------------- - // Ouput - //---------------------------------------------------- - DataFrame dataFrame = FormatOutput( param, - predictions, - const_predictions, - variance, - target_vec, - dataIn->Time(), - dataIn->TimeName() ); - - if ( param.predictOutputFile.size() ) { - // Write to disk - dataFrame.WriteData( param.pathOut, param.predictOutputFile ); - } - -#ifdef DEBUG_ALL - std::cout << dataFrame; - VectorError ve = ComputeError( - dataFrame.VectorColumnName( "Observations" ), - dataFrame.VectorColumnName( "Predictions" ) ); - std::cout << "-------------------------------------------\n"; - std::cout << "rho " << ve.rho << " RMSE " << ve.RMSE - << " MAE " << ve.MAE << std::endl; - std::cout << "-------------------------------------------\n"; -#endif +Simplex::Simplex ( DataFrame & data, bool embedded, int E, int tau ) { - return dataFrame; } diff --git a/src/Simplex.h b/src/Simplex.h new file mode 100644 index 0000000..df3be5f --- /dev/null +++ b/src/Simplex.h @@ -0,0 +1,21 @@ +// Contains the Simplex class declarations +// The Simplex class provides Simplex-algorithm-specific definitions for the +// EDM class data processing methods + +#include "Common.h" + +//---------------------------------------------------------------- +// Simplex class +// The Simplex class inherits from the EDM class and defines the +// Simplex-specific projection methods for its parent EDM class +//---------------------------------------------------------------- +class Simplex { + + public: + + //---------------------------------------------------------------- + // Simplex() : Constructor + //---------------------------------------------------------------- + Simplex ( DataFrame & data, bool embedded, int E, int tau ); + +}; diff --git a/src/Version.h b/src/Version.h deleted file mode 100644 index b367bef..0000000 --- a/src/Version.h +++ /dev/null @@ -1,22 +0,0 @@ -#ifndef VERSION_H -#define VERSION_H - -//------------------------------------------------------------ -// Instantiated in Parameters() constructor -//------------------------------------------------------------ -class Version { - int Major; - int Minor; - int Micro; - std::string Date; - -public: - Version( int Major, int Minor, int Micro, std::string Date ) : - Major( Major ), Minor( Minor ), Micro( Micro ), Date ( Date ) {}; - - void ShowVersion() { - std::cout << "cppEDM Version " << Major << "." - << Minor << "." << Micro << " " << Date << std::endl; - } -}; -#endif diff --git a/src/makefile b/src/makefile index b22079f..23053ac 100644 --- a/src/makefile +++ b/src/makefile @@ -1,12 +1,11 @@ CC = g++ -OBJ = Common.o AuxFunc.o DateTimeUtil.o Parameter.o Embed.o Interface.o\ - Neighbors.o Simplex.o Eval.o CCM.o Multiview.o SMap.o +OBJ = EDM.o EDM_Functions.o Simplex.o LIB = libEDM.a -CFLAGS = -std=c++11 -DCCM_THREADED -DMULTIVIEW_VALUES_OVERLOAD -O3 -Wreorder # -g -DDEBUG -DDEBUG_ALL -LFLAGS = -L./ -lstdc++ -lEDM -lpthread # -llapacke -llapack -lblas +CFLAGS = -std=c++11 -O3 -Wreorder -M # -g -DDEBUG -DDEBUG_ALL +LFLAGS = -L./ -lstdc++ -lEDM -lpthread all: $(LIB) ar -rcs $(LIB) $(OBJ) @@ -20,41 +19,8 @@ distclean: $(LIB): $(OBJ) -AuxFunc.o: AuxFunc.cc - $(CC) -c AuxFunc.cc $(CFLAGS) - -DateTimeUtil.o: DateTimeUtil.cc - $(CC) -c DateTimeUtil.cc $(CFLAGS) - -Common.o: Common.cc - $(CC) -c Common.cc $(CFLAGS) - -Parameter.o: Parameter.cc - $(CC) -c Parameter.cc $(CFLAGS) - -Embed.o: Embed.cc - $(CC) -c Embed.cc $(CFLAGS) - -Interface.o: Interface.cc - $(CC) -c Interface.cc $(CFLAGS) - -Neighbors.o: Neighbors.cc - $(CC) -c Neighbors.cc $(CFLAGS) - -Simplex.o: Simplex.cc - $(CC) -c Simplex.cc $(CFLAGS) - -Eval.o: Eval.cc - $(CC) -c Eval.cc $(CFLAGS) - -CCM.o: CCM.cc - $(CC) -c CCM.cc $(CFLAGS) - -Multiview.o: Multiview.cc - $(CC) -c Multiview.cc $(CFLAGS) - -SMap.o: SMap.cc - $(CC) -c SMap.cc $(CFLAGS) +%.o : %.cc %.h + $(CC) -c $< -o $@ SRCS = `echo ${OBJ} | sed -e 's/.o /.cc /g'` depend: @@ -62,18 +28,4 @@ depend: makedepend -Y $(SRCS) # DO NOT DELETE -Common.o: Common.h DataFrame.h -AuxFunc.o: AuxFunc.h Common.h DataFrame.h Neighbors.h Parameter.h Version.h -AuxFunc.o: Embed.h DateTime.h -DateTimeUtil.o: DateTime.h -Parameter.o: Parameter.h Common.h DataFrame.h Version.h -Embed.o: Embed.h Common.h DataFrame.h Parameter.h Version.h -Interface.o: Common.h DataFrame.h -Neighbors.o: Neighbors.h Common.h DataFrame.h Parameter.h Version.h -Simplex.o: Common.h DataFrame.h Parameter.h Version.h Neighbors.h Embed.h -Simplex.o: AuxFunc.h -Eval.o: Common.h DataFrame.h -CCM.o: Common.h DataFrame.h Embed.h Parameter.h Version.h AuxFunc.h -CCM.o: Neighbors.h -Multiview.o: Common.h DataFrame.h AuxFunc.h Neighbors.h Parameter.h Version.h -Multiview.o: Embed.h +EDM.o: EDM.h Common.h DataFrame.h diff --git a/src/makefile.windows b/src/makefile.windows deleted file mode 100644 index 2d9de31..0000000 --- a/src/makefile.windows +++ /dev/null @@ -1,74 +0,0 @@ - -CC = cl -OBJ = Common.obj AuxFunc.obj DateTimeUtil.obj Parameter.obj Embed.obj\ - Interface.obj Neighbors.obj Simplex.obj Eval.obj CCM.obj\ - Multiview.obj SMap.obj - -LIB = EDM.lib - -CFLAGS = -DCCM_THREADED -DMULTIVIEW_VALUES_OVERLOAD /EHsc /MD # /MT -DDEBUG -DDEBUG_ALL - -all: $(LIB) - lib /NODEFAULTLIB:LIBCMT /NODEFAULTLIB:library /OUT:$(LIB) $(OBJ) - cp $(LIB) ..\lib - -clean: - del -f $(OBJ) $(LIB) - -distclean: - del -f $(OBJ) $(LIB) ../lib/$(LIB) *~ *.bak *.csv - -$(LIB): $(OBJ) - -AuxFunc.obj: AuxFunc.cc - $(CC) /c AuxFunc.cc $(CFLAGS) - -DateTimeUtil.obj: DateTimeUtil.cc - $(CC) /c DateTimeUtil.cc $(CFLAGS) - -Common.obj: Common.cc - $(CC) /c Common.cc $(CFLAGS) - -Parameter.obj: Parameter.cc - $(CC) /c Parameter.cc $(CFLAGS) - -Embed.obj: Embed.cc - $(CC) /c Embed.cc $(CFLAGS) - -Interface.obj: Interface.cc - $(CC) /c Interface.cc $(CFLAGS) - -Neighbors.obj: Neighbors.cc - $(CC) /c Neighbors.cc $(CFLAGS) - -Simplex.obj: Simplex.cc - $(CC) /c Simplex.cc $(CFLAGS) - -Eval.obj: Eval.cc - $(CC) /c Eval.cc $(CFLAGS) - -CCM.obj: CCM.cc - $(CC) /c CCM.cc $(CFLAGS) - -Multiview.obj: Multiview.cc - $(CC) /c Multiview.cc $(CFLAGS) - -SMap.obj: SMap.cc - $(CC) /c SMap.cc $(CFLAGS) - -# Depedencies from makedepend on Linux -Common.obj: Common.h DataFrame.h -AuxFunc.obj: AuxFunc.h Common.h DataFrame.h Neighbors.h Parameter.h Version.h -AuxFunc.obj: Embed.h DateTime.h -DateTimeUtil.obj: DateTime.h -Parameter.obj: Parameter.h Common.h DataFrame.h Version.h -Embed.obj: Embed.h Common.h DataFrame.h Parameter.h Version.h -Interface.obj: Common.h DataFrame.h -Neighbors.obj: Neighbors.h Common.h DataFrame.h Parameter.h Version.h -Simplex.obj: Common.h DataFrame.h Parameter.h Version.h Neighbors.h Embed.h -Simplex.obj: AuxFunc.h -Eval.obj: Common.h DataFrame.h -CCM.obj: Common.h DataFrame.h Embed.h Parameter.h Version.h AuxFunc.h -CCM.obj: Neighbors.h -Multiview.obj: Common.h DataFrame.h AuxFunc.h Neighbors.h Parameter.h Version.h -Multiview.obj: Embed.h From 23dc6132a94bedc5863c876135a600009e1a959a Mon Sep 17 00:00:00 2001 From: cameronosmith Date: Sat, 25 Apr 2020 14:43:27 -0700 Subject: [PATCH 2/5] EDM neighbors search added, runs but not numerically verified. --- src/Common.cc | 192 +++++++++++++ src/Common.h | 24 +- src/EDM.cc | 297 +++++++++++++++++++- src/EDM.h | 40 ++- src/EDM_Functions.cc | 64 +++-- src/Embed.cc | 182 +++++++++++++ src/Embed.h | 37 +++ src/Parameter.cc | 625 +++++++++++++++++++++++++++++++++++++++++++ src/Parameter.h | 122 +++++++++ src/Simplex.cc | 11 +- src/Simplex.h | 17 +- src/Version.h | 22 ++ src/makefile | 9 +- 13 files changed, 1593 insertions(+), 49 deletions(-) create mode 100644 src/Common.cc create mode 100644 src/Embed.cc create mode 100644 src/Embed.h create mode 100644 src/Parameter.cc create mode 100644 src/Parameter.h create mode 100644 src/Version.h diff --git a/src/Common.cc b/src/Common.cc new file mode 100644 index 0000000..ebdd311 --- /dev/null +++ b/src/Common.cc @@ -0,0 +1,192 @@ + +#include +#include + +#include "Common.h" + +//---------------------------------------------------------------- +// +//---------------------------------------------------------------- +std::string ToLower( std::string str ) { + + std::string lowerStr( str ); + std::transform( lowerStr.begin(), lowerStr.end(), + lowerStr.begin(), ::tolower ); + + return lowerStr; +} + +//------------------------------------------------------------------- +// Called in two different contexts: +// DataFrame.h ReadData() Are data columns numeric or labels? +// Parameter.cc Validate() Are columns/target integer index or label? +// +// param integer : true: only digits, false: digits plus '-', '.' +// return +// true : str has only numeric characters +// false : str has non-numeric characters +//------------------------------------------------------------------- +bool OnlyDigits( std::string str, bool integer ) { + + if ( not str.size() ) { + throw std::runtime_error( "OnlyDigits(): String is empty.\n" ); + } + + // Remove whitespace + std::string str_( str ); + str_.erase( std::remove_if( str_.begin(), str_.end(), ::isspace ), + str_.end() ); + + std::string digits; + if ( integer ) { + digits = "0123456789"; + } + else { + digits = "-.0123456789"; + } + + // Is str_ purely numeric characters? + bool onlyDigits = strspn(str_.c_str(), digits.c_str()) == str_.size(); + + return onlyDigits; +} + +//---------------------------------------------------------------- +// SplitString +// +// Purpose: like Python string.split() +// +// Arguments: inString : string to be split +// delimeters : string of delimeters +// +// Note: A typical delimeter string: delimeters = " \t,\n;" +// +// Return: vector of tokens +//---------------------------------------------------------------- +std::vector SplitString( std::string inString, + std::string delimeters ) { + + size_t pos = 0; + size_t eos = 0; + size_t wordStart = 0; + size_t wordEnd = 0; + + bool foundStart = false; + bool foundEnd = false; + + std::vector splitString; + + std::string word; + + eos = inString.length(); + + while ( pos <= eos ) { + if ( not foundStart ) { + if ( delimeters.find( inString[pos] ) == delimeters.npos ) { + // this char (inString[pos]) is not a delimeter + wordStart = pos; + foundStart = true; + pos++; + continue; + } + } + if ( foundStart and not foundEnd ) { + if ( delimeters.find( inString[pos] ) != delimeters.npos + or pos == eos ) { + // this char (inString[pos]) is a delimeter or + // at the end of the string + wordEnd = pos; + foundEnd = true; + } + } + if ( foundStart and foundEnd ) { + foundStart = false; + foundEnd = false; + + word = inString.substr( wordStart, wordEnd - wordStart ); + + // remove whitespace + word.erase( std::remove_if( word.begin(), word.end(), ::isspace ), + word.end() ); + + splitString.push_back( word ); + } + if ( pos == eos ) { + break; + } + pos++; + } + + return splitString; +} + +//---------------------------------------------------------------- +// +//---------------------------------------------------------------- +VectorError ComputeError( std::valarray< double > obsIn, + std::valarray< double > predIn ) { + + // Check for nan in vectors + size_t nanObs = 0; + size_t nanPred = 0; + for ( auto o : obsIn ) { if ( std::isnan( o ) ) { nanObs++; } } + for ( auto p : predIn ) { if ( std::isnan( p ) ) { nanPred++; } } + + if ( nanObs != nanPred ) { + std::stringstream errMsg; + errMsg << "ComputeError(): obs has " << nanObs << " nan, pred has " + << nanPred << " nan. ComputeError result invalid.\n"; + std::cout << errMsg.str() << std::flush; + } + + // JP: Assume that nan are at the beginning of predictions and + // at the end of observations... probably not a robust assumption + // but this is the case if the data were prepared from Embed() + // and there were initially no nans + std::valarray pred( predIn.size() - 2 * nanPred ); + std::valarray obs ( predIn.size() - 2 * nanPred ); + if ( nanPred > 0 ) { + // ignore nanPred initial pred, and nanObs end obs + pred = predIn[ std::slice( nanPred, pred.size(), 1 ) ]; + obs = obsIn [ std::slice( nanPred, pred.size(), 1 ) ]; + } + else { + pred = std::valarray( predIn ); + obs = std::valarray( obsIn ); + } + + size_t N = pred.size(); + + std::valarray< double > two( 2, N ); // Vector of 2's for squaring + + double sumPred = pred.sum(); + double sumObs = obs.sum(); + double meanPred = sumPred / N; + double meanObs = sumObs / N; + double sumSqrPred = pow( pred, two ).sum(); + double sumSqrObs = pow( obs, two ).sum(); + double sumErr = abs( obs - pred ).sum(); + double sumSqrErr = pow( obs - pred, two ).sum(); + double sumProd = ( obs * pred ).sum(); + + double rho; // Pearson correlation coefficient + + if ( sumSqrPred * N == sumSqrPred ) { + rho = 0; + } + else { + rho = ( sumProd - N * meanObs * meanPred ) / + ( std::sqrt( ( sumSqrObs - N * pow( meanObs, 2 ) ) ) * + std::sqrt( ( sumSqrPred - N * pow( meanPred, 2 ) ) ) ); + } + + VectorError vectorError = VectorError(); + + vectorError.RMSE = sqrt( sumSqrErr / N ); + + vectorError.MAE = sumErr / N; + + vectorError.rho = rho; + + return vectorError; +} diff --git a/src/Common.h b/src/Common.h index 4f9f3d6..eb4af8b 100644 --- a/src/Common.h +++ b/src/Common.h @@ -8,7 +8,7 @@ #include #include #include -#include +#include #include #include @@ -18,10 +18,30 @@ #include "DataFrame.h" // Has #include Common.h -// Define EDM function and class signatures in its own namespace +// Data structs +struct VectorError { + double rho; + double RMSE; + double MAE; +}; + +// Enumerations +enum class Method { None, Embed, Simplex, SMap, CCM }; +enum class DistanceMetric { Euclidean, Manhattan }; + +// Namespace with all functions user should have namespace cppEDM { + // Forward declarations of EDM Functions + + DataFrame Simplex( DataFrame & data, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ) ; + } #endif diff --git a/src/EDM.cc b/src/EDM.cc index 8a7e536..e11b413 100644 --- a/src/EDM.cc +++ b/src/EDM.cc @@ -3,19 +3,300 @@ #include "EDM.h" + +namespace EDM_Neighbors { + // Define the initial maximum distance for neigbors to avoid sort() + // DBL_MAX is a Macro equivalent to: std::numeric_limits::max() + double DistanceMax = std::numeric_limits::max(); + double DistanceLimit = std::numeric_limits::max() - 1; +} + +//---------------------------------------------------------------- +// Distance computation between two vectors for several metrics +//---------------------------------------------------------------- +double Distance( const std::valarray &v1, + const std::valarray &v2, + DistanceMetric metric ) +{ + double distance = 0; + + // For efficiency sake, we forego the usual validation of v1 & v2. + + if ( metric == DistanceMetric::Euclidean ) { + double sum = 0; + double delta = 0; + for ( size_t i = 0; i < v1.size(); i++ ) { + delta = v2[i] - v1[i]; + sum += delta * delta; // avoid call to pow() + } + distance = sqrt( sum ); + + // Note: this implicit implementation is slower + // std::valarray delta = v2 - v1; + // distance = sqrt( (delta * delta).sum() ); + } + else if ( metric == DistanceMetric::Manhattan ) { + double sum = 0; + for ( size_t i = 0; i < v1.size(); i++ ) { + sum += abs( v2[i] - v1[i] ); + } + distance = sum; + } + else { + std::stringstream errMsg; + errMsg << "Distance() Invalid DistanceMetric: " + << static_cast( metric ); + throw std::runtime_error( errMsg.str() ); + } + + return distance; +} + +//---------------------------------------------------------------- +// EDM() : Constructor +// Creates the data embedding and checks parameters +// +// data : Input dataframe containing the time series to model +// method : Enum on which method is using this EDM object (for parameters) +// E : The embedding dimension to use when performing embedding +// tau : The steps between each index in the embedding. Negative +// tau is probably what you intend to use; positive tau yields +// an embedding (E_t,E_t+tau...) (future forward embedding) +// columns : Columns for embedding +// targetName : Dimension to project onto for prediction +// embedded : Whether data supplied is already embedded or not +// verbose : Verbose information flag +// +//---------------------------------------------------------------- +EDM::EDM ( DataFrame & data, int E, int tau, + std::string columns, std::string targetName, + bool embedded, bool verbose ): data(data), targetName( targetName ){ + + // Check parameters + + // Create embedding + + embedding = embedded ? data : Embed( data, E, tau, columns, verbose ); + +} + +//---------------------------------------------------------------- +// EDM() : ComputeNeighbors +// Computes neighbors for every prediction index +// lib, pred : Library and prediction ranges +// verbose : Verbose information flag +// +// return : List of DF where first element is neighbors, second is distances +// +//---------------------------------------------------------------- +EDM::Neighbors EDM::ComputeNeighbors ( std::string lib, std::string pred, + int Tp, int knn, int exclusionRadius, bool verbose ){ + + //////////////////////////////////////////////////// + // Parse prediction and lib into vectors of indices + //////////////////////////////////////////////////// + + // Holds the vector of indices for both lib and pred + std::vector< std::vector< size_t > > ranges; + + // Process both range strings + for ( std::string rangeStr : {lib, pred} ) { + + // Validate that number of ranges is even + std::vector rangeVec = SplitString( rangeStr, " \t," ); + if ( rangeVec.size() % 2 != 0 ) { + std::string errMsg( "Parameters::Validate(): " + "disjoint range must be even number of integers.\n" ); + throw std::runtime_error( errMsg ); + } + + // Generate vector of start, stop index pairs + std::vector< std::pair< size_t, size_t > > rangePairs; + for ( size_t i = 0; i < rangeVec.size(); i = i + 2 ) { + rangePairs.emplace_back( std::make_pair(std::stoi( rangeVec[i]), + std::stoi( rangeVec[i+1]))); + } + + // Create library vector of indices + + std::vector rangeIndicesVec; + + for ( auto thisPair : rangePairs ) { + for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { + + rangeIndicesVec.push_back( li - 1 ); // apply zero-offset + + } + } + ranges.push_back( rangeIndicesVec ); + } + + std::vector libraryIndices = ranges[0]; + std::vector predIndices = ranges[1]; + + auto max_lib_it = std::max_element(libraryIndices.begin(), + libraryIndices.end() ); + size_t max_lib_index = *max_lib_it; + + // Check/Set knn. Note that SMap should set knn to -1 for full library + // If knn=-1, set to full library, if knn=0, set to E+1, if E>knn>-1, error + + if ( knn == -1 ) { + knn = libraryIndices.size() - Tp * (data.NColumns() + 1); + } + else if ( knn == 0 ) { + + knn = data.NColumns() + 1; + + if ( verbose ) { + std::stringstream msg; + msg << "Parameters::Validate(): Set knn = " << knn + << " " << std::endl; + std::cout << msg.str(); + } + } + else if ( knn < data.NColumns() + 1 ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): knn of " << knn + << " is less than E+1 = " << data.NColumns() + 1 << std::endl; + throw std::runtime_error( errMsg.str() ); + } + + // Neighbors: struct on local stack to be returned by copy + Neighbors neighbors; + neighbors.neighbors = DataFrame (predIndices.size(), knn); + neighbors.distances = DataFrame (predIndices.size(), knn); + + // Vectors to hold indices and values from each comparison + std::valarray k_NN_neighbors( knn ); + std::valarray k_NN_distances( knn ); + + //------------------------------------------------------------------- + // For each prediction vector (row in prediction DataFrame) find the + // list of library indices that are within k_NN points + //------------------------------------------------------------------- + for ( size_t pred_row_idx = 0; pred_row_idx < predIndices.size(); pred_row_idx++ ) { + + // Get the current query/pred row + size_t pred_row = predIndices[ pred_row_idx ]; + std::valarray pred_vec = data.Row( pred_row ); + + // Reset the neighbor and distance vectors for this pred row + for ( size_t i = 0; i < knn; i++ ) { + k_NN_neighbors[ i ] = 0; + // JP: Used to avoid sort() + k_NN_distances[ i ] = EDM_Neighbors::DistanceMax; + } + + //-------------------------------------------------------------- + // Compute distance on every library row to query/pred row + //-------------------------------------------------------------- + for ( size_t row_j = 0; row_j < libraryIndices.size(); row_j++ ) { + // Get the library vector for this lib_row index + size_t lib_row = libraryIndices[ row_j ]; + std::valarray lib_vec = data.Row( lib_row ); + + // If the library point is degenerate with the prediction, + // ignore it. + if ( lib_row == pred_row ) { +#ifdef DEBUG_ALL + if ( verbose ) { + std::stringstream msg; + msg << "FindNeighbors(): Ignoring degenerate lib_row " + << lib_row << " and pred_row " << pred_row << std::endl; + std::cout << msg.str(); + } +#endif + continue; + } + + // Apply temporal exclusion radius: units are data rows, not time + if ( exclusionRadius ) { + int xrad = (int) lib_row - pred_row; + if ( std::abs( xrad ) <= exclusionRadius ) { + continue; + } + } + // If this lib_row + args.Tp >= max_lib_index, then this neighbor + // would be outside the library, keep looking if noNeighborLimit + if ( not noNeighborLimit ) { + if ( lib_row + Tp > max_lib_index ) { + continue; + } + } + + // Find distance between the prediction vector + // and each of the library vectors + // The 1st column (j=0) of Time has been excluded above + double d_i = Distance( lib_vec, pred_vec, + DistanceMetric::Euclidean ); + + // If d_i is less than values in k_NN_distances, add to list + auto max_it = std::max_element( begin( k_NN_distances ), + end( k_NN_distances ) ); + if ( d_i < *max_it ) { + size_t max_i = std::distance( begin(k_NN_distances), max_it ); + k_NN_neighbors[ max_i ] = lib_row; // Save the index + k_NN_distances[ max_i ] = d_i; // Save the value + } + } // for ( row_j = 0; row_j < library.size(); row_j++ ) + + // Assert that at least knn were found + if ( *std::max_element( begin( k_NN_distances ), + end ( k_NN_distances ) ) > + EDM_Neighbors::DistanceLimit ) { + + std::stringstream errMsg; + errMsg << "FindNeighbors(): Failed to find " + << knn << " knn neighbors. The library " + << "may be too small." << std::endl; + throw std::runtime_error( errMsg.str() ); + } + + // Check for ties. JP: Need to address this, not just warning + // First sort a copy of k_NN_neighbors so unique() will work + std::valarray k_NN_neighborCopy( k_NN_neighbors ); + std::sort( begin( k_NN_neighborCopy ), end( k_NN_neighborCopy ) ); + + // ui is iterator to first non unique element + auto ui = std::unique( begin( k_NN_neighborCopy ), + end ( k_NN_neighborCopy ) ); + + if ( std::distance( begin( k_NN_neighborCopy ), ui ) != + k_NN_neighborCopy.size() ) { + std::cout << "WARNING: FindNeighbors(): Degenerate neighbors./n"; + } + + // Write the neighbor indices and distance values + neighbors.neighbors.WriteRow( pred_row_idx, k_NN_neighbors ); + neighbors.distances.WriteRow( pred_row_idx, k_NN_distances ); + + } // for ( pred_row_idx = 0; pred_row_idx < predictionRows->size(); pred_row_idx++ ) + + return neighbors; + +} + //---------------------------------------------------------------- -// EDM() : Constructor +// EDM() : Project +// Finds neighbors, performs weighting on neighbors, and +// projects onto pred range. // // data : Input dataframe containing the time series to model. -// embedded : Flag on whether the input data is already embedded. -// E : The embedding dimension to use when performing embedding. -// tau : The steps between each index in the embedding. Negative // tau is probably what you intend to use; positive tau yields // an embedding (E_t,E_t+tau...) (future forward embedding). //---------------------------------------------------------------- -EDM::EDM ( DataFrame & data, bool embedded, int E, int tau ) : - data( data ) { +std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, + int Tp, int knn, int exclusionRadius, bool verbose ) +{ + + // Find neighbors + ComputeNeighbors(lib,pred,Tp,knn,exclusionRadius,verbose); - // check parameters - + // Weight neighbors + + // Project onto target for prediction + + return std::list>(); } + diff --git a/src/EDM.h b/src/EDM.h index fca7c9b..9681560 100644 --- a/src/EDM.h +++ b/src/EDM.h @@ -6,6 +6,7 @@ #define EDM_H #include "Common.h" +#include "Embed.h" //---------------------------------------------------------------- // EDM Class @@ -17,18 +18,45 @@ //---------------------------------------------------------------- class EDM { - // The input time series (potentially not embedded) + // The input time series and its embedding DataFrame & data; + DataFrame embedding; - // The target column this class should model for + // The dimension to be project onto std::string targetName; + // Flag on whether to have no neighbor limit in neighbor search + bool noNeighborLimit = false; + + // Return structure of FindNeighbors() + struct Neighbors { + DataFrame neighbors; + DataFrame distances; + Neighbors(); + ~Neighbors(); + }; + public: - //---------------------------------------------------------------- - // EDM() : Constructor - //---------------------------------------------------------------- - EDM ( DataFrame & data, bool embedded, int E, int tau ); + //---------------------------------------------------------------- + // EDM() : Constructor + //---------------------------------------------------------------- + EDM ( DataFrame & data, int E, int tau, + std::string columns, std::string targetName, + bool embedded, bool verbose ); + + //---------------------------------------------------------------- + // EDM() : ComputeNeighbors + //---------------------------------------------------------------- + Neighbors ComputeNeighbors ( + std::string lib, std::string pred, int Tp, int knn, + int exclusionRadius, bool verbose ); + + //---------------------------------------------------------------- + // EDM() : Project + //---------------------------------------------------------------- + std::list< DataFrame > Project ( std::string lib, std::string pred, + int Tp, int knn, int exclusionRadius, bool verbose ); }; diff --git a/src/EDM_Functions.cc b/src/EDM_Functions.cc index a6f8068..21546a0 100644 --- a/src/EDM_Functions.cc +++ b/src/EDM_Functions.cc @@ -2,38 +2,48 @@ // such that the user does not have to deal with the OOD. #include "Common.h" +#include "Simplex.h" -//---------------------------------------------------------------- -// Simplex prediction algorithm with data input as DataFrame -// See EDM class for parameter definitions. -//---------------------------------------------------------------- -DataFrame Simplex( DataFrame & data, - std::string pathOut, std::string predictFile, - std::string lib, std::string pred, - int E, int Tp, int knn, int tau, int exclusionRadius, - std::string columns, std::string target, - bool embedded, bool const_predict, bool verbose ) { +namespace cppEDM { - // Create Simplex object and get its projection + //---------------------------------------------------------------- + // Simplex prediction algorithm with data input as DataFrame + // See EDM class for parameter definitions. + //---------------------------------------------------------------- + DataFrame Simplex( DataFrame & data, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ) { -} + // Create Simplex object and get its projection -//---------------------------------------------------------------- -// Simplex prediction algorithm with data input as filepath -// See EDM class for parameter definitions. -//---------------------------------------------------------------- -DataFrame Simplex( std::string pathIn, std::string dataFile, - std::string pathOut, std::string predictFile, - std::string lib, std::string pred, - int E, int Tp, int knn, int tau, int exclusionRadius, - std::string columns, std::string target, - bool embedded, bool const_predict, bool verbose ) { + SimplexMachine simplexManager( data, pathOut, predictFile, + lib, pred, E, Tp, knn, tau, exclusionRadius, columns, target, + embedded, const_predict, verbose ); - // Just load in dataframe and delegate to Simplex with DataFrame input + return DataFrame(); - DataFrame data (pathIn, dataFile); + } - return Simplex ( data, pathOut, predictFile, lib, pred, E, Tp, knn, tau, - exclusionRadius, columns, target, embedded, const_predict, verbose ); -} + //---------------------------------------------------------------- + // Simplex prediction algorithm with data input as filepath + // See EDM class for parameter definitions. + //---------------------------------------------------------------- + DataFrame Simplex( std::string pathIn, std::string dataFile, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ) { + + // Just load in dataframe and delegate to Simplex with DataFrame input + DataFrame data (pathIn, dataFile); + + return cppEDM::Simplex ( data, pathOut, predictFile, lib, pred, E, Tp, knn, tau, + exclusionRadius, columns, target, embedded, const_predict, verbose ); + } + +} diff --git a/src/Embed.cc b/src/Embed.cc new file mode 100644 index 0000000..ba1b854 --- /dev/null +++ b/src/Embed.cc @@ -0,0 +1,182 @@ + +// NOTE: The returned data block does NOT have the time column +#include "Embed.h" + +//---------------------------------------------------------------- +// API Overload 1: Explicit data file path/name +// Embed DataFrame columns in E dimensions. +// Data is read from path/dataFile +// Implemented as a wrapper to API Overload 2: +// which is a wrapper for MakeBlock() +// +// NOTE: Truncates data by tau * (E-1) rows to remove +// nan values (partial data rows) +//---------------------------------------------------------------- +DataFrame< double > Embed( std::string path, + std::string dataFile, + int E, // embedding dimension + int tau, // time step delay + std::string columns, // column names or indices + bool verbose ) { + + DataFrame< double > dataFrame( path, dataFile ); + DataFrame< double > embedded = Embed(dataFrame, E, tau, columns, verbose); + return embedded; +} + +//---------------------------------------------------------------- +// API Overload 2: DataFrame provided +// Implemented as a wrapper for MakeBlock() +// Note: dataFrame must have the columnNameToIndex map +//---------------------------------------------------------------- +DataFrame< double > Embed( DataFrame< double > dataFrameIn, + int E, + int tau, + std::string columns, + bool verbose ) { + + // Parameter.Validate will convert columns into a vector of names + // or a vector of column indices + Parameters param = Parameters( Method::Embed, "", "", "", "", + "1 1", "1 1", E, 0, 0, tau, 0, 0, + columns, "", false, false, verbose ); + + if ( not param.columnIndex.size() and + dataFrameIn.ColumnNameToIndex().empty() ) { + throw std::runtime_error("Embed(DataFrame): columnNameIndex empty.\n"); + } + + // If columns provided, validate they are in dataFrameIn + for ( auto colName : param.columnNames ) { + auto ci = find( dataFrameIn.ColumnNames().begin(), + dataFrameIn.ColumnNames().end(), colName ); + + if ( ci == dataFrameIn.ColumnNames().end() ) { + std::stringstream errMsg; + errMsg << "Embed(DataFrame): Failed to find column " + << colName << " in dataFrame with columns: [ "; + for ( auto col : dataFrameIn.ColumnNames() ) { + errMsg << col << " "; + } errMsg << " ]\n"; + throw std::runtime_error( errMsg.str() ); + } + } + + // Get column names for MakeBlock + std::vector< std::string > colNames; + if ( param.columnNames.size() ) { + // column names are strings use as-is + colNames = param.columnNames; + } + else if ( param.columnIndex.size() ) { + // columns are indices : Create column names for MakeBlock + for ( size_t i = 0; i < param.columnIndex.size(); i++ ) { + std::stringstream ss; + ss << "V" << param.columnIndex[i]; + colNames.push_back( ss.str() ); + } + } + else { + throw std::runtime_error( "Embed(DataFrame): columnNames and " + " columnIndex are empty.\n" ); + } + + // Extract the specified columns (sub)DataFrame from dataFrameIn + DataFrame< double > dataFrame; + + if ( param.columnNames.size() ) { + dataFrame = dataFrameIn.DataFrameFromColumnNames( param.columnNames ); + } + else if ( param.columnIndex.size() ) { + // already have column indices + // Note there will be no column names transferred + dataFrame = dataFrameIn.DataFrameFromColumnIndex( param.columnIndex ); + } + + DataFrame< double > embedding = MakeBlock( dataFrame, E, tau, + colNames, verbose ); + + return embedding; +} + +//------------------------------------------------------------------------ +// MakeBlock from dataFrame +// Ignores the first (or last) tau * (E-1) dataFrame rows of partial data. +// Does not validate parameters or columns, use Embed() +//------------------------------------------------------------------------ +DataFrame< double > MakeBlock( DataFrame< double > dataFrame, + int E, + int tau, + std::vector columnNames, + bool verbose ) { + + if ( columnNames.size() != dataFrame.NColumns() ) { + std::stringstream errMsg; + errMsg << "MakeBlock: The number of columns in the dataFrame (" + << dataFrame.NColumns() << ") is not equal to the number " + << "of columns specified (" << columnNames.size() << ").\n";; + throw std::runtime_error( errMsg.str() ); + } + + if ( E < 1 ) { + std::stringstream errMsg; + errMsg << "MakeBlock(): E = " << E << " is invalid.\n" ; + throw std::runtime_error( errMsg.str() ); + } + + size_t NRows = dataFrame.NRows(); // number of input rows + size_t NColOut = dataFrame.NColumns() * E; // number of output columns + size_t NPartial = abs( tau ) * (E-1); // rows to shift & delete + + // Create embedded data frame column names X(t-0) X(t-1)... + std::vector< std::string > newColumnNames( NColOut ); + size_t newCol_i = 0; + for ( size_t col = 0; col < columnNames.size(); col ++ ) { + for ( size_t e = 0; e < E; e++ ) { + std::stringstream ss; + if ( tau < 0 ) { + ss << columnNames[ col ] << "(t-" << e << ")"; + } + else { + ss << columnNames[ col ] << "(t+" << e << ")"; + } + newColumnNames[ newCol_i ] = ss.str(); + newCol_i++; + } + } + + // Ouput data frame with tau * E-1 fewer rows + DataFrame< double > embedding( NRows - NPartial, NColOut, newColumnNames ); + + // To keep track of where to insert column in new data frame + size_t colCount = 0; + + // Slice to ignore rows with partial data + std::slice slice_i; + if ( tau < 0 ) { + slice_i = std::slice( NPartial, NRows - NPartial, 1 ); + } + else { + slice_i = std::slice( 0, NRows - NPartial, 1 ); + } + + // Shift column data and write to embedding data frame + for ( size_t col = 0; col < dataFrame.NColumns(); col++ ) { + // for each embedding dimension + for ( size_t e = 0; e < E; e++ ) { + + std::valarray< double > column = dataFrame.Column( col ); + + // Returns a copy of the valarray object with its elements + // shifted left n spaces (or right if n is negative). + std::valarray< double > tmp = column.shift( e * tau ); + + // Write shifted columns to the output embedding DataFrame + embedding.WriteColumn( colCount, tmp[ slice_i ] ); + + colCount++; + } + } + + return embedding; +} diff --git a/src/Embed.h b/src/Embed.h new file mode 100644 index 0000000..dd08100 --- /dev/null +++ b/src/Embed.h @@ -0,0 +1,37 @@ +#ifndef EMBED_H +#define EMBED_H + +#include "Common.h" +#include "Parameter.h" + +//---------------------------------------------------------------- +// API Overload 1: Explicit data file path/name +// Implemented as a wrapper to API Overload 2: +// which is a wrapper for MakeBlock() +//---------------------------------------------------------------- +DataFrame< double > Embed ( std::string path = "", + std::string dataFile = "", + int E = 0, + int tau = 0, + std::string columns = "", + bool verbose = false ); + +//---------------------------------------------------------------- +// API Overload 2: DataFrame provided +// Implemented as a wrapper for MakeBlock() +//---------------------------------------------------------------- +DataFrame< double > Embed ( DataFrame< double > dataFrame, + int E = 0, + int tau = 0, + std::string columns = "", + bool verbose = false ); + +//---------------------------------------------------------------- +//---------------------------------------------------------------- +DataFrame< double > MakeBlock ( DataFrame< double > dataFrame, + int E, + int tau, + std::vector columnNames, + bool verbose ); +#endif + diff --git a/src/Parameter.cc b/src/Parameter.cc new file mode 100644 index 0000000..4b1b30c --- /dev/null +++ b/src/Parameter.cc @@ -0,0 +1,625 @@ + +#include "Parameter.h" + +//---------------------------------------------------------------- +// Constructor +// Default values set in Parameter.h +//---------------------------------------------------------------- +Parameters::Parameters( + Method method, + std::string pathIn, + std::string dataFile, + std::string pathOut, + std::string predictFile, + std::string lib_str, + std::string pred_str, + int E, + int Tp, + int knn, + int tau, + double theta, + int exclusionRadius, + + std::string columns_str, + std::string target_str, + + bool embedded, + bool const_predict, + bool verbose, + + std::string SmapFile, + std::string blockFile, + std::string derivatives_str, + + double svdSig, + double tikhonov, + double elasticNet, + + int multi, + std::string libSizes_str, + int sample, + bool random, + bool replacement, + unsigned rseed, + bool noNeigh + ) : + // Variable initialization from Parameters arguments + method ( method ), + pathIn ( pathIn ), + dataFile ( dataFile ), + pathOut ( pathOut ), + predictOutputFile( predictFile ), + lib_str ( lib_str ), + pred_str ( pred_str ), + E ( E ), + Tp ( Tp ), + knn ( knn ), + tau ( tau ), + theta ( theta ), + exclusionRadius ( exclusionRadius ), + + columns_str ( columns_str ), + target_str ( target_str ), + targetIndex ( 0 ), + + embedded ( embedded ), + const_predict ( const_predict ), + verbose ( verbose ), + + SmapOutputFile ( SmapFile ), + blockOutputFile ( blockFile ), + + derivatives_str ( derivatives_str ), + SVDSignificance ( svdSig ), + TikhonovAlpha ( tikhonov ), + ElasticNetAlpha ( elasticNet ), + + MultiviewEnsemble( multi ), + libSizes_str ( libSizes_str ), + subSamples ( sample ), + randomLib ( random ), + replacement ( replacement ), + seed ( rseed ), + noNeighborLimit ( noNeigh ), + + // Set validated flag and instantiate Version + validated ( false ), + version ( 1, 3, 5, "2020-04-16" ) +{ + // Constructor code + if ( method != Method::None ) { + + Validate(); + + if ( verbose ) { + version.ShowVersion(); + } + } +} + +//---------------------------------------------------------------- +// Destructor +//---------------------------------------------------------------- +Parameters::~Parameters() {} + +//---------------------------------------------------------------- +// +//---------------------------------------------------------------- +void Parameters::Load() {} + +//---------------------------------------------------------------- +// Index offsets, generate library and prediction indices, +// and parameter validation +//---------------------------------------------------------------- +void Parameters::Validate() { + + validated = true; + + if ( not embedded and tau == 0 ) { + std::string errMsg( "Parameters::Validate(): " + "tau must be non-zero.\n" ); + throw std::runtime_error( errMsg ); + } + + if ( Tp < 0 ) { + std::string errMsg( "Parameters::Validate(): " + "Tp must be positive.\n" ); + throw std::runtime_error( errMsg ); + } + + //-------------------------------------------------------------- + // Generate library indices: Apply zero-offset + //-------------------------------------------------------------- + if ( lib_str.size() ) { + // Parse lib_str into vector of strings + std::vector lib_vec = SplitString( lib_str, " \t," ); + if ( lib_vec.size() % 2 != 0 ) { + std::string errMsg( "Parameters::Validate(): " + "library must be even number of integers.\n" ); + throw std::runtime_error( errMsg ); + } + + // Generate vector of start, stop index pairs + std::vector< std::pair< size_t, size_t > > libPairs; + for ( size_t i = 0; i < lib_vec.size(); i = i + 2 ) { + libPairs.emplace_back( std::make_pair( std::stoi( lib_vec[i] ), + std::stoi( lib_vec[i+1] ) ) ); + } + + size_t nLib = 0; // Count of lib items + + // Get number of lib indices, validate end > start + for ( auto thisPair : libPairs ) { + size_t lib_start = thisPair.first; + size_t lib_end = thisPair.second; + + nLib += lib_end - lib_start + 1; + + // Validate end > stop indices + if ( method == Method::Simplex or method == Method::SMap ) { + // Don't check if method == None, Embed or CCM since default + // of "1 1" is used. + if ( lib_start >= lib_end ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): library start " + << lib_start << " exceeds end " << lib_end << ".\n"; + throw std::runtime_error( errMsg.str() ); + } + } + } + + // Create library vector of indices + library = std::vector< size_t >( nLib ); + size_t i = 0; + for ( auto thisPair : libPairs ) { + for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { + library[ i ] = li - 1; // apply zero-offset + i++; + } + } + } + + //-------------------------------------------------------------- + // Generate prediction indices: Apply zero-offset + //-------------------------------------------------------------- + if ( pred_str.size() ) { + // Parse pred_str into vector of strings + std::vector pred_vec = SplitString( pred_str, " \t," ); + if ( pred_vec.size() % 2 != 0 ) { + std::string errMsg( "Parameters::Validate(): " + "prediction must be even number of integers.\n"); + throw std::runtime_error( errMsg ); + } + + // Generate vector of start, stop index pairs + std::vector< std::pair< size_t, size_t > > predPairs; + for ( size_t i = 0; i < pred_vec.size(); i = i + 2 ) { + predPairs.emplace_back( std::make_pair( std::stoi( pred_vec[i] ), + std::stoi( pred_vec[i+1]))); + } + + size_t nPred = 0; // Count of pred items + + // Get number of pred indices, validate end > start + for ( auto thisPair : predPairs ) { + size_t pred_start = thisPair.first; + size_t pred_end = thisPair.second; + + nPred += pred_end - pred_start + 1; + + // Validate end > stop indices + if ( method == Method::Simplex or method == Method::SMap ) { + // Don't check if method == None, Embed or CCM since default + // of "1 1" is used. + if ( pred_start >= pred_end ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): prediction start " + << pred_start << " exceeds end " << pred_end << ".\n"; + throw std::runtime_error( errMsg.str() ); + } + } + } + + // Create prediction vector of indices + prediction = std::vector< size_t >( nPred ); + size_t i = 0; + for ( auto thisPair : predPairs ) { + for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { + prediction[ i ] = li - 1; // apply zero-offset + i++; + } + } + } + + if ( method == Method::Simplex or method == Method::SMap ) { + if ( not library.size() ) { + std::string errMsg( "Parameters::Validate(): " + "library indices not found.\n" ); + throw std::runtime_error( errMsg ); + } + if ( not prediction.size() ) { + std::string errMsg( "Parameters::Validate(): " + "prediction indices not found.\n" ); + throw std::runtime_error( errMsg ); + } + } + else { + // Defaults if Method is None, Embed or CCM + if ( not library.size() ) { + library = std::vector( 1, 0 ); + } + if ( not prediction.size() ) { + prediction = std::vector( 1, 0 ); + } + } + +#ifdef DEBUG_ALL + PrintIndices( library, prediction ); +#endif + + //-------------------------------------------------------------- + // Convert multi argument parameters from string to vectors + //-------------------------------------------------------------- + + // Columns + // If columns are purely integer, then populate vector columnIndex + // Otherwise fill in vector columnNames + if ( columns_str.size() ) { + + std::vector columns_vec = SplitString( columns_str, + " \t,\n" ); + + bool onlyDigits = false; + + for ( auto ci = columns_vec.begin(); ci != columns_vec.end(); ++ci ) { + onlyDigits = OnlyDigits( *ci, true ); + + if ( not onlyDigits ) { break; } + } + + if ( onlyDigits ) { + for ( auto ci = columns_vec.begin(); + ci != columns_vec.end(); ++ci ) { + columnIndex.push_back( std::stoi( *ci ) ); + } + } + else { + columnNames = columns_vec; + } + } + + if ( not columnIndex.size() and not columnNames.size() ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): Simplex/CCM: " + << " No valid columns found." << std::endl; + throw std::runtime_error( errMsg.str() ); + } + + // target + if ( target_str.size() ) { + bool onlyDigits = OnlyDigits( target_str, true ); + if ( onlyDigits ) { + targetIndex = std::stoi( target_str ); + } + else { + targetName = target_str; + } + } + + // Derivatives + if ( derivatives_str.size() > 0 ) { + std::vector der_vec = SplitString(derivatives_str," \t,"); + if ( der_vec.size() < 2 ) { + std::string errMsg( "Parameters::Validate(): " + "derivatives must be integer pairs.\n"); + throw std::runtime_error( errMsg ); + } + derivatives = std::vector( der_vec.size() ); + for ( size_t i = 0; i < der_vec.size(); i++ ) { + derivatives.push_back( std::stoi( der_vec[i] ) ); + } + } + + //-------------------------------------------------------------------- + // CCM sample not 0 if random is true + if ( method == Method::CCM ) { + if ( randomLib ) { + if ( subSamples < 1 ) { + std::string errMsg( "Parameters::Validate(): " + "CCM samples must be > 0.\n" ); + throw std::runtime_error( errMsg ); + } + } + } + + // CCM librarySizes + if ( libSizes_str.size() > 0 ) { + std::vector libsize_vec = SplitString(libSizes_str," \t,"); + if ( libsize_vec.size() != 3 ) { + std::string errMsg( "Parameters::Validate(): " + "CCM librarySizes must be three integers.\n" ); + throw std::runtime_error( errMsg ); + } + + size_t start = std::stoi( libsize_vec[0] ); + size_t stop = std::stoi( libsize_vec[1] ); + size_t increment = std::stoi( libsize_vec[2] ); + + if ( increment < 1 ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): " + << "CCM librarySizes increment " << increment + << " is invalid.\n"; + throw std::runtime_error( errMsg.str() ); + } + + if ( start > stop ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): " + << "CCM librarySizes start " << start + << " stop " << stop << " are invalid.\n"; + throw std::runtime_error( errMsg.str() ); + } + + size_t N_lib = std::floor( (stop-start)/increment + 1/increment ) + 1; + + if ( start < E ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): " + << "CCM librarySizes start < E = " << E << "\n"; + throw std::runtime_error( errMsg.str() ); + } + else if ( start < 3 ) { + std::string errMsg( "Parameters::Validate(): " + "CCM librarySizes start < 3.\n" ); + throw std::runtime_error( errMsg ); + } + + // Create the librarySizes vector + librarySizes = std::vector( N_lib, 0 ); + + // Fill in the sizes + size_t libSize = start; + for ( size_t i = 0; i < librarySizes.size(); i++ ) { + librarySizes[i] = libSize; + libSize = libSize + increment; + } + } + + //-------------------------------------------------------------------- + // Simplex and knn not specified: not embedded : knn set to E+1 + // embedded : knn set to size( columns ) + // S-Map require knn > E + 1, default is all neighbors. + if ( method == Method::Simplex or method == Method::CCM ) { + if ( knn < 1 ) { + if ( not embedded ) { + knn = E + 1; + if ( verbose ) { + std::stringstream msg; + msg << "Parameters::Validate(): Set knn = " << knn + << " (E+1) for Simplex. " << std::endl; + std::cout << msg.str(); + } + } + else { // embedded = true + if ( columnIndex.size() ) { + knn = columnIndex.size() + 1; + } + else if ( columnNames.size() ) { + knn = columnNames.size() + 1; + } + + if ( verbose ) { + std::stringstream msg; + msg << "Parameters::Validate(): Set knn = " << knn + << " for Simplex (embedded = true). " << std::endl; + std::cout << msg.str(); + } + } + } + if ( knn < E + 1 ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): Simplex knn of " << knn + << " is less than E+1 = " << E + 1 << std::endl; + throw std::runtime_error( errMsg.str() ); + } + } + else if ( method == Method::SMap ) { + if ( knn > 0 ) { + if ( knn < E + 1 ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate() S-Map knn must be at least " + " E+1 = " << E + 1 << ".\n"; + throw std::runtime_error( errMsg.str() ); + } + } + else { + // default knn = 0, set knn value + knn = library.size() - Tp * (E + 1); + if ( verbose ) { + std::stringstream msg; + msg << "Parameters::Validate(): Set knn = " << knn + << " for SMap. " << std::endl; + std::cout << msg.str(); + } + } + if ( not embedded and columnNames.size() > 1 ) { + std::string msg( "Parameters::Validate() WARNING: " + "Multivariable S-Map should use " + "-e (embedded) data input to ensure " + "data/dimension correspondance.\n" ); + std::cout << msg; + } + + // S-Map coefficient columns for derivatives start at 1 since the 0th + // column is the S-Map linear prediction bias term + if ( derivatives.size() > 1 ) { + std::vector::iterator it = std::find( derivatives.begin(), + derivatives.end(), 0); + if ( it != derivatives.end() ) { + std::string errMsg( "Parameters::Validate() S-Map coefficient " + " columns for derivatives can not use column 0.\n"); + throw std::runtime_error( errMsg ); + } + if ( derivatives.size() % 2 ) { + std::string errMsg( "Parameters::Validate() S-Map coefficient " + " columns for derivatives must be in pairs.\n"); + throw std::runtime_error( errMsg ); + } + } + + // Tikhonov and ElasticNet are mutually exclusive + if ( TikhonovAlpha and ElasticNetAlpha ) { + std::string errMsg( "Parameters::Validate() Multiple S-Map solve " + "methods specified. Use one or none of: " + "tikhonov, elasticNet.\n"); + throw std::runtime_error( errMsg ); + } + + // Very small alphas don't make sense in elastic net + if ( ElasticNetAlpha < 0.01 ) { + std::cout << "Parameters::Validate() ElasticNetAlpha too small." + " Setting to 0.01."; + ElasticNetAlpha = 0.01; + } + if ( ElasticNetAlpha > 1 ) { + std::cout << "Parameters::Validate() ElasticNetAlpha too large." + " Setting to 1."; + ElasticNetAlpha = 1; + } + } + else if ( method == Method::Embed ) { + // no-op + } + else { + throw std::runtime_error( "Parameters::Validate() " + "Prediction method error.\n" ); + } +} + +//------------------------------------------------------------ +// Adjust lib/pred concordant with Embed() removal of tau(E-1) +// rows, and DeletePartialDataRow() +//------------------------------------------------------------ +void Parameters::DeleteLibPred( size_t shift ) { + + size_t library_len = library.size(); + size_t prediction_len = prediction.size(); + + // If [0, 1, ... shift] (negative tau) or + // [N-shift, ... N-1, N] (positive tau) are in library or prediction + // those rows were deleted, delete these index elements. + // First, create vectors of indices to delete. + std::vector< size_t > deleted_pred_elements( shift, 0 ); + std::vector< size_t > deleted_lib_elements ( shift, 0 ); + + if ( tau < 0 ) { + std::iota(deleted_pred_elements.begin(), deleted_pred_elements.end(),0); + std::iota(deleted_lib_elements.begin(), deleted_lib_elements.end(), 0); + } + else { + std::iota( deleted_pred_elements.begin(), + deleted_pred_elements.end(), prediction_len - shift ); + std::iota( deleted_lib_elements.begin(), + deleted_lib_elements.end(), library_len - shift ); + } + + // Erase elements of row indices that were deleted + for ( auto element = deleted_lib_elements.begin(); + element != deleted_lib_elements.end(); element++ ) { + + std::vector< size_t >::iterator it; + it = std::find( library.begin(), library.end(), *element ); + + if ( it != library.end() ) { + library.erase( it ); + } + } + + for ( auto element = deleted_pred_elements.begin(); + element != deleted_pred_elements.end(); element++ ) { + + std::vector< size_t >::iterator it; + it = std::find( prediction.begin(), + prediction.end(), *element ); + + if ( it != prediction.end() ) { + prediction.erase( it ); + } + } + + // Now offset all values by shift so that vectors indices + // in library and prediction refer to the same data rows + // before the deletion/shift. + if ( tau < 0 ) { + for ( auto li = library.begin(); li != library.end(); li++ ) { + *li = *li - shift; + } + for ( auto pi = prediction.begin(); pi != prediction.end(); pi++ ) { + *pi = *pi - shift; + } + } + // tau > 0 : Forward shifting: no adjustment needed from origin +} + +//------------------------------------------------------------------ +// Overload << to output to ostream +//------------------------------------------------------------------ +std::ostream& operator<< ( std::ostream &os, Parameters &p ) { + + // print info about the dataframe + os << "Parameters: -------------------------------------------\n"; + + std::string method("Unknown"); + if ( p.method == Method::Simplex ) { method = "Simplex"; } + else if ( p.method == Method::SMap ) { method = "SMap"; } + else if ( p.method == Method::CCM ) { method = "CCM"; } + else if ( p.method == Method::None ) { method = "None"; } + else if ( p.method == Method::Embed ) { method = "Embed"; } + + os << "Method: " << method + << " E=" << p.E << " Tp=" << p.Tp + << " knn=" << p.knn << " tau=" << p.tau << " theta=" << p.theta + << std::endl; + + if ( p.columnNames.size() ) { + os << "Column Names : [ "; + for ( auto ci = p.columnNames.begin(); + ci != p.columnNames.end(); ++ci ) { + os << *ci << " "; + } os << "]" << std::endl; + } + + if ( p.targetName.size() ) { + os << "Target: " << p.targetName << std::endl; + } + + os << "Library: [" << p.library[0] << " : " + << p.library[ p.library.size() - 1 ] << "] " + << "Prediction: [" << p.prediction[0] << " : " + << p.prediction[ p.prediction.size() - 1 ] + << "] " << std::endl; + + + os << "-------------------------------------------------------\n"; + + return os; +} + +#ifdef DEBUG_ALL +//------------------------------------------------------------------ +// +//------------------------------------------------------------------ +void Parameters::PrintIndices( std::vector library, + std::vector prediction ) +{ + std::cout << "Parameters(): library: "; + for ( auto li = library.begin(); li != library.end(); ++li ) { + std::cout << *li << " "; + } std::cout << std::endl; + std::cout << "Parameters(): prediction: "; + for ( auto pi = prediction.begin(); pi != prediction.end(); ++pi ) { + std::cout << *pi << " "; + } std::cout << std::endl; +} +#endif diff --git a/src/Parameter.h b/src/Parameter.h new file mode 100644 index 0000000..60125c3 --- /dev/null +++ b/src/Parameter.h @@ -0,0 +1,122 @@ +#ifndef PARAMETER_H +#define PARAMETER_H + +#include +#include + +#include "Common.h" +#include "Version.h" + +//------------------------------------------------------------ +// +//------------------------------------------------------------ +class Parameters { + +public: // Not protected with accessors. + Method method; // Simplex or SMap enum class + + std::string pathIn; // path for input dataFile + std::string dataFile; // input dataFile (assumed .csv) + std::string pathOut; // path for output files + std::string predictOutputFile;// + std::string lib_str; // multi argument parameters for library + std::string pred_str; // multi argument parameters for prediction + + std::vector library; // library row indices + std::vector prediction; // prediction row indices + + int E; // dimension + int Tp; // prediction interval + int knn; // k nearest neighbors + int tau; // block embedding delay + double theta; // S Map localization + int exclusionRadius; // temporal rows to ignore in predict + + std::string columns_str; + std::string target_str; + std::vector columnNames; // column names + std::vector columnIndex; // column indices + + std::string targetName; // target column name + size_t targetIndex; // target column index + + bool embedded; // true if data is already embedded/block + bool const_predict; // true to compute non "predictor" stats + bool verbose; + + std::string SmapOutputFile; // + std::string blockOutputFile; // Embed() output file + + std::string derivatives_str; + std::vector derivatives;// list of column indices for derivatives + + double SVDSignificance; // SVD singular value cutoff + double TikhonovAlpha; // Initial alpha parameter + double ElasticNetAlpha; // Initial alpha parameter + + int MultiviewEnsemble; // Number of ensembles in multiview + + std::string libSizes_str; + std::vector librarySizes;// CCM library sizes to evaluate + int subSamples; // CCM number of samples to draw + bool randomLib; // CCM randomly select subsets if true + bool replacement; // CCM random select with replacement if true + unsigned seed; // CCM random selection RNG seed + + bool noNeighborLimit; // Strictly forbid neighbors outside library + bool validated; + + Version version; // Version object, instantiated in constructor + + friend std::ostream& operator<<(std::ostream &os, Parameters ¶ms); + + // Constructor declaration and default arguments + Parameters( + Method method = Method::None, + std::string pathIn = "./", + std::string dataFile = "", + std::string pathOut = "./", + std::string predictFile = "", + std::string lib_str = "", + std::string pred_str = "", + int E = 0, + int Tp = 0, + int knn = 0, + int tau = -1, + double theta = 0, + int exclusionRadius = 0, + + std::string columns_str = "", + std::string target_str = "", + + bool embedded = false, + bool const_predict = false, + bool verbose = false, + + std::string SmapFile = "", + std::string blockFile = "", + std::string derivatives_str = "", + + double svdSig = 1E-5, + double tikhonov = 0, + double elasticNet = 0.1, + + int multi = 0, + std::string libSizes_str = "", + int sample = 0, + bool random = true, + bool replacement = false, + unsigned seed = 0, // 0: Generate random seed in CCM + bool noNeighbor = false + ); + + ~Parameters(); + + void Validate(); // Parameter validation and index offsets + void Load(); // Populate the parameters from arguments + void DeleteLibPred( size_t shift ); // Adjust for embedding + void PrintIndices( std::vector library, + std::vector prediction ); +}; + +#endif diff --git a/src/Simplex.cc b/src/Simplex.cc index abb3f91..d45b3c3 100644 --- a/src/Simplex.cc +++ b/src/Simplex.cc @@ -6,6 +6,15 @@ // Simplex() : Constructor // See EDM class for descriptions of parameters not described //---------------------------------------------------------------- -Simplex::Simplex ( DataFrame & data, bool embedded, int E, int tau ) { +SimplexMachine::SimplexMachine ( DataFrame & data, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ): + EDM{data, E, tau, columns, target, embedded, verbose } { + + Project(lib,pred,Tp,knn,exclusionRadius,verbose); + } diff --git a/src/Simplex.h b/src/Simplex.h index df3be5f..1ade394 100644 --- a/src/Simplex.h +++ b/src/Simplex.h @@ -2,20 +2,31 @@ // The Simplex class provides Simplex-algorithm-specific definitions for the // EDM class data processing methods +#ifndef SIMPLEX_H +#define SIMPLEX_H + #include "Common.h" +#include "EDM.h" //---------------------------------------------------------------- // Simplex class // The Simplex class inherits from the EDM class and defines the // Simplex-specific projection methods for its parent EDM class //---------------------------------------------------------------- -class Simplex { +class SimplexMachine : public EDM { public: //---------------------------------------------------------------- - // Simplex() : Constructor + // Simplex() : Constructor with DataFrame input //---------------------------------------------------------------- - Simplex ( DataFrame & data, bool embedded, int E, int tau ); + SimplexMachine ( DataFrame &data, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ); }; + +#endif diff --git a/src/Version.h b/src/Version.h new file mode 100644 index 0000000..b367bef --- /dev/null +++ b/src/Version.h @@ -0,0 +1,22 @@ +#ifndef VERSION_H +#define VERSION_H + +//------------------------------------------------------------ +// Instantiated in Parameters() constructor +//------------------------------------------------------------ +class Version { + int Major; + int Minor; + int Micro; + std::string Date; + +public: + Version( int Major, int Minor, int Micro, std::string Date ) : + Major( Major ), Minor( Minor ), Micro( Micro ), Date ( Date ) {}; + + void ShowVersion() { + std::cout << "cppEDM Version " << Major << "." + << Minor << "." << Micro << " " << Date << std::endl; + } +}; +#endif diff --git a/src/makefile b/src/makefile index 23053ac..8f023ee 100644 --- a/src/makefile +++ b/src/makefile @@ -1,6 +1,6 @@ CC = g++ -OBJ = EDM.o EDM_Functions.o Simplex.o +OBJ = Common.o Embed.o EDM.o Simplex.o EDM_Functions.o Parameter.o LIB = libEDM.a @@ -28,4 +28,9 @@ depend: makedepend -Y $(SRCS) # DO NOT DELETE -EDM.o: EDM.h Common.h DataFrame.h +Common.o: Common.h DataFrame.h +Parameter.o: Common.o +Embed.o: Parameter.o +EDM.o: Embed.o +Simplex.o: EDM.o +EDM_Functions.o: Simplex.o From 23b9314b842b05a484e1a8bb0e1462dc5e2ed4c0 Mon Sep 17 00:00:00 2001 From: cameronosmith Date: Thu, 30 Apr 2020 21:37:46 -0700 Subject: [PATCH 3/5] Up to distance function added. Parameter object mostly removed from project, at least from EDM.cc as of now. Integrated all of Parameter into EDM class --- src/Common.h | 9 ++ src/EDM.cc | 281 ++++++++++++++++++++++++++------------------- src/EDM.h | 21 ++-- src/EDM_Helpers.cc | 238 ++++++++++++++++++++++++++++++++++++++ src/Simplex.cc | 2 +- src/makefile | 7 +- 6 files changed, 431 insertions(+), 127 deletions(-) create mode 100644 src/EDM_Helpers.cc diff --git a/src/Common.h b/src/Common.h index eb4af8b..334466f 100644 --- a/src/Common.h +++ b/src/Common.h @@ -4,6 +4,7 @@ #define COMMON_H #include +#include #include #include #include @@ -42,6 +43,14 @@ namespace cppEDM { std::string columns, std::string target, bool embedded, bool const_predict, bool verbose ) ; + DataFrame Simplex( std::string pathIn, std::string dataFile, + std::string pathOut, std::string predictFile, + std::string lib, std::string pred, + int E, int Tp, int knn, int tau, int exclusionRadius, + std::string columns, std::string target, + bool embedded, bool const_predict, bool verbose ) ; + + } #endif diff --git a/src/EDM.cc b/src/EDM.cc index e11b413..e10d788 100644 --- a/src/EDM.cc +++ b/src/EDM.cc @@ -2,7 +2,7 @@ // description of EDM class and how to use the methods. #include "EDM.h" - +#include "EDM_Helpers.cc" namespace EDM_Neighbors { // Define the initial maximum distance for neigbors to avoid sort() @@ -11,46 +11,9 @@ namespace EDM_Neighbors { double DistanceLimit = std::numeric_limits::max() - 1; } -//---------------------------------------------------------------- -// Distance computation between two vectors for several metrics -//---------------------------------------------------------------- -double Distance( const std::valarray &v1, - const std::valarray &v2, - DistanceMetric metric ) -{ - double distance = 0; - - // For efficiency sake, we forego the usual validation of v1 & v2. - - if ( metric == DistanceMetric::Euclidean ) { - double sum = 0; - double delta = 0; - for ( size_t i = 0; i < v1.size(); i++ ) { - delta = v2[i] - v1[i]; - sum += delta * delta; // avoid call to pow() - } - distance = sqrt( sum ); - - // Note: this implicit implementation is slower - // std::valarray delta = v2 - v1; - // distance = sqrt( (delta * delta).sum() ); - } - else if ( metric == DistanceMetric::Manhattan ) { - double sum = 0; - for ( size_t i = 0; i < v1.size(); i++ ) { - sum += abs( v2[i] - v1[i] ); - } - distance = sum; - } - else { - std::stringstream errMsg; - errMsg << "Distance() Invalid DistanceMetric: " - << static_cast( metric ); - throw std::runtime_error( errMsg.str() ); - } - - return distance; -} +#ifdef DEBUG_ALL +void PrintNeighbors( const EDM::Neighbors &neighbors ); +#endif //---------------------------------------------------------------- // EDM() : Constructor @@ -63,90 +26,166 @@ double Distance( const std::valarray &v1, // tau is probably what you intend to use; positive tau yields // an embedding (E_t,E_t+tau...) (future forward embedding) // columns : Columns for embedding -// targetName : Dimension to project onto for prediction // embedded : Whether data supplied is already embedded or not // verbose : Verbose information flag // //---------------------------------------------------------------- -EDM::EDM ( DataFrame & data, int E, int tau, - std::string columns, std::string targetName, - bool embedded, bool verbose ): data(data), targetName( targetName ){ +EDM::EDM ( DataFrame & data, int E, int tau, + std::string columns, std::string targetName, + bool embedded, bool verbose ): + data(data), targetName( targetName ), E(E), tau(tau), embedded(embedded) { + + ///////////////////////////////////////////// + // Validate parameters and create embedding + ///////////////////////////////////////////// + + if ( not embedded and tau == 0 ) { + std::string errMsg( "Parameters::Validate(): " + "tau must be non-zero.\n" ); + throw std::runtime_error( errMsg ); + } - // Check parameters + std::vector columnNames = ParseColumnNames( columns ); + std::vector columnIndex = ParseColumnIndices( columns ); - // Create embedding + if ( not columnIndex.size() and not columnNames.size() ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): No valid columns found." << std::endl; + throw std::runtime_error( errMsg.str() ); + } - embedding = embedded ? data : Embed( data, E, tau, columns, verbose ); + if ( embedded ) { + // dataIn is a multivariable block, no embedding needed + // Select the specified columns into dataBlock + embedding = columnNames.size() ? + data.DataFrameFromColumnNames( columnNames ) : + data.DataFrameFromColumnIndex( columnIndex ); + } + else { + // embedded = false: Create the embedding dataBlock via Embed() + // dataBlock will have tau * (E-1) fewer rows than dataIn + embedding = Embed( data, E, tau, columns, verbose ); + } } //---------------------------------------------------------------- -// EDM() : ComputeNeighbors -// Computes neighbors for every prediction index -// lib, pred : Library and prediction ranges -// verbose : Verbose information flag +// EDM() : Project +// Finds neighbors, performs weighting on neighbors, and +// projects onto pred range. // -// return : List of DF where first element is neighbors, second is distances -// +// targetName : Dimension to project onto for prediction +// +// data : Input dataframe containing the time series to model. +// tau is probably what you intend to use; positive tau yields +// an embedding (E_t,E_t+tau...) (future forward embedding). //---------------------------------------------------------------- -EDM::Neighbors EDM::ComputeNeighbors ( std::string lib, std::string pred, - int Tp, int knn, int exclusionRadius, bool verbose ){ +std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, + std::string target, int Tp, int knn, + int exclusionRadius, bool verbose ){ + + // Validate Tp + if ( Tp < 0 ) { + std::string errMsg( "Parameters::Validate(): " + "Tp must be positive.\n" ); + throw std::runtime_error( errMsg ); + } + + // Parse lib and pred range strings + std::vector libIndices = ParseRangeStr( lib ); + std::vector predIndices = ParseRangeStr( pred ); - //////////////////////////////////////////////////// - // Parse prediction and lib into vectors of indices - //////////////////////////////////////////////////// - - // Holds the vector of indices for both lib and pred - std::vector< std::vector< size_t > > ranges; - - // Process both range strings - for ( std::string rangeStr : {lib, pred} ) { - - // Validate that number of ranges is even - std::vector rangeVec = SplitString( rangeStr, " \t," ); - if ( rangeVec.size() % 2 != 0 ) { - std::string errMsg( "Parameters::Validate(): " - "disjoint range must be even number of integers.\n" ); - throw std::runtime_error( errMsg ); - } + CheckDataRows( data.NRows(), libIndices.back(), predIndices.back(), + E, tau, embedded ); + + //---------------------------------------------------------- + // Get target (library) vector + //---------------------------------------------------------- + + std::valarray targetIn; + + std::vector< size_t > colIndices = ParseColumnIndices(target); + std::vector< std::string > colNames = ParseColumnNames(target); + + // Note CS : check with JP to make sure target range selection appropriate + // If a column name or idx specified, extract that vector + if ( colIndices.size() or colNames.size() ) { + + // Get name of target column in embedding + + std::string targName; - // Generate vector of start, stop index pairs - std::vector< std::pair< size_t, size_t > > rangePairs; - for ( size_t i = 0; i < rangeVec.size(); i = i + 2 ) { - rangePairs.emplace_back( std::make_pair(std::stoi( rangeVec[i]), - std::stoi( rangeVec[i+1]))); + // If embedded, then column names are not formatted + if (embedded) { + targName = colNames.size()? colNames[0] : + embedding.ColumnNames()[ colIndices[0] ]; + } + // If we performed embedding, reproduce formatting of target col name + // as prefix "V" if idx specified, and append (t+-0) + else { + std::string originalName = colNames.size() ? colNames[0] : + "V" + std::to_string( colIndices[0] ); + std::string tauSign = tau > 0 ? "+" : ""; + targName = originalName+ "(t" + tauSign + std::to_string(tau) + ")"; } - // Create library vector of indices + targetIn = embedding.VectorColumnName( targName ); + } + // Default to first column + else { + targetIn = data.Column( 0 ); + } - std::vector rangeIndicesVec; + // Shift lib and pred by number of partial rows if we embedded - for ( auto thisPair : rangePairs ) { - for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { + if ( not embedded ) { - rangeIndicesVec.push_back( li - 1 ); // apply zero-offset + size_t shift = abs( tau ) * ( E - 1 ); - } + if ( shift > 0 ) { + libIndices = AdjustRange( shift, tau, libIndices); + predIndices = AdjustRange( shift, tau, predIndices); } - ranges.push_back( rangeIndicesVec ); + + // Check boundaries again since rows were potentially removed + CheckDataRows( data.NRows(), libIndices.back(), predIndices.back(), + E, tau, embedded ); } - std::vector libraryIndices = ranges[0]; - std::vector predIndices = ranges[1]; + // Find neighbors + ComputeNeighbors(libIndices, predIndices, Tp, knn, exclusionRadius, verbose); - auto max_lib_it = std::max_element(libraryIndices.begin(), - libraryIndices.end() ); - size_t max_lib_index = *max_lib_it; + // Weight neighbors + + // Project onto target for prediction + + return std::list>(); +} + +//---------------------------------------------------------------- +// EDM() : ComputeNeighbors +// Computes neighbors for every prediction index +// lib, pred : Library and prediction ranges +// verbose : Verbose information flag +// +// return : List of DF where first element is neighbors, second is distances +// +//---------------------------------------------------------------- +EDM::Neighbors EDM::ComputeNeighbors ( + std::vector libraryIndices, std::vector predIndices, + int Tp, int knn, int exclusionRadius, bool verbose ){ + //------------------------------------------------------------------------- // Check/Set knn. Note that SMap should set knn to -1 for full library // If knn=-1, set to full library, if knn=0, set to E+1, if E>knn>-1, error + //------------------------------------------------------------------------- if ( knn == -1 ) { - knn = libraryIndices.size() - Tp * (data.NColumns() + 1); + knn = libraryIndices.size() - Tp * (E + 1); } else if ( knn == 0 ) { - knn = data.NColumns() + 1; + knn = E + 1; if ( verbose ) { std::stringstream msg; @@ -155,13 +194,17 @@ EDM::Neighbors EDM::ComputeNeighbors ( std::string lib, std::string pred, std::cout << msg.str(); } } - else if ( knn < data.NColumns() + 1 ) { + else if ( knn < E + 1 ) { std::stringstream errMsg; errMsg << "Parameters::Validate(): knn of " << knn - << " is less than E+1 = " << data.NColumns() + 1 << std::endl; + << " is less than E+1 = " << E + 1 << std::endl; throw std::runtime_error( errMsg.str() ); } + //------------------------------------------------------------------------- + // Variables for holding neighbor results and definitions for convenience + //------------------------------------------------------------------------- + // Neighbors: struct on local stack to be returned by copy Neighbors neighbors; neighbors.neighbors = DataFrame (predIndices.size(), knn); @@ -171,15 +214,21 @@ EDM::Neighbors EDM::ComputeNeighbors ( std::string lib, std::string pred, std::valarray k_NN_neighbors( knn ); std::valarray k_NN_distances( knn ); + // Max index for out of bounds query in exclusion function + auto max_lib_it = std::max_element(libraryIndices.begin(), + libraryIndices.end() ); + size_t max_lib_index = *max_lib_it; + //------------------------------------------------------------------- // For each prediction vector (row in prediction DataFrame) find the // list of library indices that are within k_NN points //------------------------------------------------------------------- for ( size_t pred_row_idx = 0; pred_row_idx < predIndices.size(); pred_row_idx++ ) { + // Get the current query/pred row size_t pred_row = predIndices[ pred_row_idx ]; - std::valarray pred_vec = data.Row( pred_row ); + std::valarray pred_vec = embedding.Row( pred_row ); // Reset the neighbor and distance vectors for this pred row for ( size_t i = 0; i < knn; i++ ) { @@ -194,7 +243,7 @@ EDM::Neighbors EDM::ComputeNeighbors ( std::string lib, std::string pred, for ( size_t row_j = 0; row_j < libraryIndices.size(); row_j++ ) { // Get the library vector for this lib_row index size_t lib_row = libraryIndices[ row_j ]; - std::valarray lib_vec = data.Row( lib_row ); + std::valarray lib_vec = embedding.Row( lib_row ); // If the library point is degenerate with the prediction, // ignore it. @@ -273,30 +322,32 @@ EDM::Neighbors EDM::ComputeNeighbors ( std::string lib, std::string pred, } // for ( pred_row_idx = 0; pred_row_idx < predictionRows->size(); pred_row_idx++ ) + #ifdef DEBUG_ALL + const Neighbors &neigh = neighbors; + PrintNeighbors( neigh ); + #endif + return neighbors; } + //---------------------------------------------------------------- -// EDM() : Project -// Finds neighbors, performs weighting on neighbors, and -// projects onto pred range. -// -// data : Input dataframe containing the time series to model. -// tau is probably what you intend to use; positive tau yields -// an embedding (E_t,E_t+tau...) (future forward embedding). +// Debug method to print out all neighbors //---------------------------------------------------------------- -std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, - int Tp, int knn, int exclusionRadius, bool verbose ) +#ifdef DEBUG_ALL +void PrintNeighbors( const EDM::Neighbors &neighbors ) { - - // Find neighbors - ComputeNeighbors(lib,pred,Tp,knn,exclusionRadius,verbose); - - // Weight neighbors - - // Project onto target for prediction - - return std::list>(); + std::cout << "FindNeighbors(): neighbors:distances" << std::endl; + //for ( size_t i = 0; i < neighbors.neighbors.NRows(); i++ ) { + for ( size_t i = 0; i < 5; i++ ) { + std::cout << "Row " << i << " | "; + for ( size_t j = 0; j < neighbors.neighbors.NColumns(); j++ ) { + std::cout << neighbors.neighbors( i, j ) << " "; + } std::cout << " : "; + for ( size_t j = 0; j < neighbors.neighbors.NColumns(); j++ ) { + std::cout << neighbors.distances( i, j ) << " "; + } std::cout << std::endl; + } } - +#endif diff --git a/src/EDM.h b/src/EDM.h index 9681560..5ecb174 100644 --- a/src/EDM.h +++ b/src/EDM.h @@ -22,22 +22,25 @@ class EDM { DataFrame & data; DataFrame embedding; + // The current embedding dimension (meaning may change when implementing + // EmbedDimension) and tau + int E, tau; + bool embedded; + // The dimension to be project onto std::string targetName; // Flag on whether to have no neighbor limit in neighbor search bool noNeighborLimit = false; + public: + // Return structure of FindNeighbors() struct Neighbors { DataFrame neighbors; DataFrame distances; - Neighbors(); - ~Neighbors(); }; - public: - //---------------------------------------------------------------- // EDM() : Constructor //---------------------------------------------------------------- @@ -49,14 +52,16 @@ class EDM { // EDM() : ComputeNeighbors //---------------------------------------------------------------- Neighbors ComputeNeighbors ( - std::string lib, std::string pred, int Tp, int knn, - int exclusionRadius, bool verbose ); + std::vector libraryIndices, std::vector predIndices, + int Tp, int knn, int exclusionRadius, bool verbose ); + //---------------------------------------------------------------- // EDM() : Project //---------------------------------------------------------------- - std::list< DataFrame > Project ( std::string lib, std::string pred, - int Tp, int knn, int exclusionRadius, bool verbose ); + std::list< DataFrame > Project (std::string lib,std::string pred, + std::string target, int Tp, int knn, + int exclusionRadius, bool verbose ); }; diff --git a/src/EDM_Helpers.cc b/src/EDM_Helpers.cc new file mode 100644 index 0000000..e40dd3d --- /dev/null +++ b/src/EDM_Helpers.cc @@ -0,0 +1,238 @@ +// This file is for the helper functions (parsing ranges, validating ranges,...) +// that are only used in EDM which prevents us adding it to Common, but is +// extraneous to the actual EDM behavior +// Doesn't warrant a header file since only used in EDM.cc + +#include "Common.h" + +//---------------------------------------------------------------- +// Distance computation between two vectors for several metrics +//---------------------------------------------------------------- +double Distance( const std::valarray &v1, + const std::valarray &v2, + DistanceMetric metric ) +{ + double distance = 0; + + // For efficiency sake, we forego the usual validation of v1 & v2. + + if ( metric == DistanceMetric::Euclidean ) { + double sum = 0; + double delta = 0; + for ( size_t i = 0; i < v1.size(); i++ ) { + delta = v2[i] - v1[i]; + sum += delta * delta; // avoid call to pow() + } + distance = sqrt( sum ); + + // Note: this implicit implementation is slower + // std::valarray delta = v2 - v1; + // distance = sqrt( (delta * delta).sum() ); + } + else if ( metric == DistanceMetric::Manhattan ) { + double sum = 0; + for ( size_t i = 0; i < v1.size(); i++ ) { + sum += abs( v2[i] - v1[i] ); + } + distance = sum; + } + else { + std::stringstream errMsg; + errMsg << "Distance() Invalid DistanceMetric: " + << static_cast( metric ); + throw std::runtime_error( errMsg.str() ); + } + + return distance; +} + +//---------------------------------------------------------- +// Method to validate dataFrameIn rows against lib and pred indices +// prediction_max_i and library_max_i are the last time index for each range +// numRows : number of rows in the input dataframe +// maxPredIdx : the max index in the pred indices +// maxLibIdx : the max index in the lib indices +// return : none, throws runtime error if invalid ranges +//---------------------------------------------------------- +void CheckDataRows( size_t numRows, size_t maxPredIdx, size_t maxLibIdx, + int E, int tau, bool embedded ){ + + // param.prediction & library have been zero-offset in Validate() + // to convert from user specified data row to array indicies + + size_t shift = embedded ? 0 : abs( tau ) * ( E - 1 ); + + if ( numRows <= maxPredIdx ) { + std::stringstream errMsg; + errMsg << "CheckDataRows(): The prediction index " << maxPredIdx + 1 + << " exceeds the number of data rows " << numRows; + throw std::runtime_error( errMsg.str() ); + } + + if ( numRows <= maxLibIdx + shift ) { + std::stringstream errMsg; + errMsg << "CheckDataRows(): The library index " << maxLibIdx + 1 + << " + tau(E-1) " << shift << " = " << maxLibIdx + 1 + shift + << " exceeds the number of data rows " << numRows; + throw std::runtime_error( errMsg.str() ); + } + +} + +//---------------------------------------------------------------- +// Function to parse a potentially disjoint range string into its indices +// rangeStr : range string to parse +// return : vector of corresponding indices for given range string +//---------------------------------------------------------------- +std::vector ParseRangeStr ( std::string rangeStr ) { + + // Validate that number of ranges is even + std::vector rangeVec = SplitString( rangeStr, " \t," ); + if ( rangeVec.size() % 2 != 0 ) { + std::string errMsg( "Parameters::Validate(): " + "disjoint range must be even number of integers.\n" ); + throw std::runtime_error( errMsg ); + } + + // Generate vector of start, stop index pairs + std::vector< std::pair< size_t, size_t > > rangePairs; + for ( size_t i = 0; i < rangeVec.size(); i = i + 2 ) { + rangePairs.emplace_back( std::make_pair(std::stoi( rangeVec[i]), + std::stoi( rangeVec[i+1]))); + } + + // Create library vector of indices + + std::vector rangeIndicesVec; + + for ( auto thisPair : rangePairs ) { + for ( size_t li = thisPair.first; li <= thisPair.second; li++ ) { + + rangeIndicesVec.push_back( li - 1 ); // apply zero-offset + + } + } + + // Assert size of range + if ( not rangeIndicesVec.size() ) { + std::string errMsg( "Parameters::Validate(): " + "prediction or library range not found.\n" ); + throw std::runtime_error( errMsg ); + } + + return rangeIndicesVec; + +} + +//---------------------------------------------------------------- +// Function to parse multi-arg str into int vector (column indices) +// If column string contains only names, not indices, empty vector returned +//---------------------------------------------------------------- +std::vector ParseColumnIndices ( std::string indicesStr ) { + + // If columns are purely integer, then populate vector columnIndex + // Else return empty vector as no indices specified + if ( indicesStr.size() ) { + + std::vector columnSplit = SplitString( indicesStr, + " \t,\n" ); + bool onlyDigits = false; + + for ( auto ci = columnSplit.begin(); ci != columnSplit.end(); ++ci ) { + onlyDigits = OnlyDigits( *ci, true ); + + if ( not onlyDigits ) { break; } + } + + if ( onlyDigits ) { + std::vector columnIndices; + for ( auto columnIndex : columnSplit ) { + columnIndices.push_back( std::stoi( columnIndex ) ); + } + return columnIndices; + } + } + + // Either column str empty or contained alpha characters + return std::vector(); +} +//---------------------------------------------------------------- +// Function to parse multi-arg str into str vector (column names) +// If column string contains only indices, not names, empty vector returned +//---------------------------------------------------------------- +std::vector ParseColumnNames ( std::string columnsStr ) { + + // Fill columns into vector columnNames + if ( columnsStr.size() ) { + + std::vector columnsVec = SplitString( columnsStr, + " \t,\n" ); + + bool onlyDigits = false; + + for ( auto ci = columnsVec.begin(); ci != columnsVec.end(); ++ci ) { + onlyDigits = OnlyDigits( *ci, true ); + + if ( not onlyDigits ) { break; } + } + + if ( not onlyDigits ) { + return columnsVec; + } + } + + // Either column str empty or contained only indices + return std::vector(); +} +//------------------------------------------------------------ +// Adjust lib/pred concordant with Embed() removal of tau(E-1) +// rows, and DeletePartialDataRow() +// shift : the shift amount (number of deleted rows) +// tau : tau +// rangeVec : the range vector (lib or pred) +// return : the adjusted range vectors +//------------------------------------------------------------ +std::vector AdjustRange(size_t shift, int tau, + std::vector rangeVec){ + + // If [0, 1, ... shift] (negative tau) or + // [N-shift, ... N-1, N] (positive tau) are in library or prediction + // those rows were deleted, delete these index elements. + // First, create vectors of indices to delete. + std::vector< size_t > deleted_range_elements( shift, 0 ); + + if ( tau < 0 ) { + std::iota(deleted_range_elements.begin(),deleted_range_elements.end(),0); + } + else { + std::iota( deleted_range_elements.begin(), deleted_range_elements.end(), + rangeVec.size() - shift); + } + + // Erase elements of row indices that were deleted + for ( auto element = deleted_range_elements.begin(); + element != deleted_range_elements.end(); element++ ) { + + std::vector< size_t >::iterator it; + it = std::find( rangeVec.begin(), rangeVec.end(), *element ); + + if ( it != rangeVec.end() ) { + rangeVec.erase( it ); + } + } + + // Now offset all values by shift so that vectors indices + // in library and prediction refer to the same data rows + // before the deletion/shift. + if ( tau < 0 ) { + for ( auto ri = rangeVec.begin(); ri != rangeVec.end(); ri++ ) { + *ri = *ri - shift; + } + } + // tau > 0 : Forward shifting: no adjustment needed from origin + + return rangeVec; +} + + + diff --git a/src/Simplex.cc b/src/Simplex.cc index d45b3c3..16fad86 100644 --- a/src/Simplex.cc +++ b/src/Simplex.cc @@ -14,7 +14,7 @@ SimplexMachine::SimplexMachine ( DataFrame & data, bool embedded, bool const_predict, bool verbose ): EDM{data, E, tau, columns, target, embedded, verbose } { - Project(lib,pred,Tp,knn,exclusionRadius,verbose); + Project(lib,pred,target,Tp,knn,exclusionRadius,verbose); } diff --git a/src/makefile b/src/makefile index 8f023ee..a18f40b 100644 --- a/src/makefile +++ b/src/makefile @@ -1,10 +1,10 @@ CC = g++ -OBJ = Common.o Embed.o EDM.o Simplex.o EDM_Functions.o Parameter.o +OBJ = Common.o Embed.o EDM.o Simplex.o EDM_Functions.o Parameter.o EDM_Helpers.o LIB = libEDM.a -CFLAGS = -std=c++11 -O3 -Wreorder -M # -g -DDEBUG -DDEBUG_ALL +CFLAGS = -std=c++11 -O3 -Wreorder -M -g -DDEBUG -DDEBUG_ALL LFLAGS = -L./ -lstdc++ -lEDM -lpthread all: $(LIB) @@ -31,6 +31,7 @@ depend: Common.o: Common.h DataFrame.h Parameter.o: Common.o Embed.o: Parameter.o -EDM.o: Embed.o +EDM_Helpers.o: Common.o +EDM.o: Embed.o EDM_Helpers.o Simplex.o: EDM.o EDM_Functions.o: Simplex.o From 2a97ddef522f4ee225db9c3668fddeb222d3ffb4 Mon Sep 17 00:00:00 2001 From: cameronosmith Date: Thu, 30 Apr 2020 22:20:50 -0700 Subject: [PATCH 4/5] Typo in target column name selection --- src/EDM.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EDM.cc b/src/EDM.cc index e10d788..5afa58c 100644 --- a/src/EDM.cc +++ b/src/EDM.cc @@ -125,8 +125,8 @@ std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, else { std::string originalName = colNames.size() ? colNames[0] : "V" + std::to_string( colIndices[0] ); - std::string tauSign = tau > 0 ? "+" : ""; - targName = originalName+ "(t" + tauSign + std::to_string(tau) + ")"; + char tauSign = tau > 0 ? '+' : '-'; + targName = originalName+ "(t" + tauSign + "0)"; } targetIn = embedding.VectorColumnName( targName ); From 6e6a21037089dd085861a011d31f8e34f577ae6b Mon Sep 17 00:00:00 2001 From: SoftwareLiteracy Date: Tue, 19 May 2020 12:51:47 -0400 Subject: [PATCH 5/5] v2.0.0 Tweak --- src/EDM.cc | 63 +++++++++++++++++++++++++++--------------------- src/EDM.h | 47 ++++++++++++++++++++++-------------- src/Embed.cc | 6 +++-- src/Parameter.cc | 14 +++-------- src/Parameter.h | 10 +++----- src/Simplex.cc | 33 ++++++++++++++++--------- src/makefile | 35 ++++++++++++++++----------- 7 files changed, 119 insertions(+), 89 deletions(-) diff --git a/src/EDM.cc b/src/EDM.cc index 5afa58c..93ac486 100644 --- a/src/EDM.cc +++ b/src/EDM.cc @@ -30,14 +30,21 @@ void PrintNeighbors( const EDM::Neighbors &neighbors ); // verbose : Verbose information flag // //---------------------------------------------------------------- -EDM::EDM ( DataFrame & data, int E, int tau, - std::string columns, std::string targetName, - bool embedded, bool verbose ): - data(data), targetName( targetName ), E(E), tau(tau), embedded(embedded) { - - ///////////////////////////////////////////// +EDM::EDM ( DataFrame & data, + int E, + int tau, + bool embedded, + std::string columns, + std::string target, + bool noNeighborLimit, + bool verbose ): + data(data), E(E), tau(tau), embedded(embedded), + columns(columns), target(target), noNeighborLimit(noNeighborLimit), + verbose(verbose) { + + //------------------------------------------------------------- // Validate parameters and create embedding - ///////////////////////////////////////////// + //------------------------------------------------------------- if ( not embedded and tau == 0 ) { std::string errMsg( "Parameters::Validate(): " @@ -66,7 +73,6 @@ EDM::EDM ( DataFrame & data, int E, int tau, // dataBlock will have tau * (E-1) fewer rows than dataIn embedding = Embed( data, E, tau, columns, verbose ); } - } //---------------------------------------------------------------- @@ -74,22 +80,19 @@ EDM::EDM ( DataFrame & data, int E, int tau, // Finds neighbors, performs weighting on neighbors, and // projects onto pred range. // -// targetName : Dimension to project onto for prediction +// target : Dimension to project onto for prediction // // data : Input dataframe containing the time series to model. -// tau is probably what you intend to use; positive tau yields +// -tau is probably what you intend to use; positive tau yields // an embedding (E_t,E_t+tau...) (future forward embedding). //---------------------------------------------------------------- -std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, - std::string target, int Tp, int knn, - int exclusionRadius, bool verbose ){ - - // Validate Tp - if ( Tp < 0 ) { - std::string errMsg( "Parameters::Validate(): " - "Tp must be positive.\n" ); - throw std::runtime_error( errMsg ); - } +std::list< DataFrame > EDM::Project ( std::string lib, + std::string pred, + std::string target, + int Tp, + int knn, + int exclusionRadius, + bool verbose ) { // Parse lib and pred range strings std::vector libIndices = ParseRangeStr( lib ); @@ -153,7 +156,7 @@ std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, } // Find neighbors - ComputeNeighbors(libIndices, predIndices, Tp, knn, exclusionRadius, verbose); + FindNeighbors( libIndices, predIndices, Tp, knn, exclusionRadius, verbose ); // Weight neighbors @@ -163,7 +166,7 @@ std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, } //---------------------------------------------------------------- -// EDM() : ComputeNeighbors +// EDM() : FindNeighbors // Computes neighbors for every prediction index // lib, pred : Library and prediction ranges // verbose : Verbose information flag @@ -171,9 +174,12 @@ std::list< DataFrame > EDM::Project ( std::string lib, std::string pred, // return : List of DF where first element is neighbors, second is distances // //---------------------------------------------------------------- -EDM::Neighbors EDM::ComputeNeighbors ( - std::vector libraryIndices, std::vector predIndices, - int Tp, int knn, int exclusionRadius, bool verbose ){ +EDM::Neighbors EDM::FindNeighbors ( std::vector libraryIndices, + std::vector predIndices, + int Tp, + int knn, + int exclusionRadius, + bool verbose ) { //------------------------------------------------------------------------- // Check/Set knn. Note that SMap should set knn to -1 for full library @@ -223,7 +229,8 @@ EDM::Neighbors EDM::ComputeNeighbors ( // For each prediction vector (row in prediction DataFrame) find the // list of library indices that are within k_NN points //------------------------------------------------------------------- - for ( size_t pred_row_idx = 0; pred_row_idx < predIndices.size(); pred_row_idx++ ) { + for ( size_t pred_row_idx = 0; pred_row_idx < predIndices.size(); + pred_row_idx++ ) { // Get the current query/pred row @@ -231,7 +238,7 @@ EDM::Neighbors EDM::ComputeNeighbors ( std::valarray pred_vec = embedding.Row( pred_row ); // Reset the neighbor and distance vectors for this pred row - for ( size_t i = 0; i < knn; i++ ) { + for ( int i = 0; i < knn; i++ ) { k_NN_neighbors[ i ] = 0; // JP: Used to avoid sort() k_NN_distances[ i ] = EDM_Neighbors::DistanceMax; @@ -312,7 +319,7 @@ EDM::Neighbors EDM::ComputeNeighbors ( end ( k_NN_neighborCopy ) ); if ( std::distance( begin( k_NN_neighborCopy ), ui ) != - k_NN_neighborCopy.size() ) { + (long int) k_NN_neighborCopy.size() ) { std::cout << "WARNING: FindNeighbors(): Degenerate neighbors./n"; } diff --git a/src/EDM.h b/src/EDM.h index 5ecb174..7e3c14a 100644 --- a/src/EDM.h +++ b/src/EDM.h @@ -24,14 +24,15 @@ class EDM { // The current embedding dimension (meaning may change when implementing // EmbedDimension) and tau - int E, tau; + int E; + int tau; bool embedded; + + std::string columns; // Timeseries to be embedded, or the embedding + std::string target; // The dimension to be projected onto - // The dimension to be project onto - std::string targetName; - - // Flag on whether to have no neighbor limit in neighbor search - bool noNeighborLimit = false; + bool noNeighborLimit; + bool verbose; public: @@ -44,25 +45,35 @@ class EDM { //---------------------------------------------------------------- // EDM() : Constructor //---------------------------------------------------------------- - EDM ( DataFrame & data, int E, int tau, - std::string columns, std::string targetName, - bool embedded, bool verbose ); + EDM ( DataFrame & data, + int E, + int tau, + bool embedded, + std::string columns, + std::string target, + bool noNeighborLimit, + bool verbose ); //---------------------------------------------------------------- - // EDM() : ComputeNeighbors + // EDM() : FindNeighbors //---------------------------------------------------------------- - Neighbors ComputeNeighbors ( - std::vector libraryIndices, std::vector predIndices, - int Tp, int knn, int exclusionRadius, bool verbose ); + Neighbors FindNeighbors ( std::vector libraryIndices, + std::vector predIndices, + int Tp, + int knn, + int exclusionRadius, + bool verbose ); //---------------------------------------------------------------- // EDM() : Project //---------------------------------------------------------------- - std::list< DataFrame > Project (std::string lib,std::string pred, - std::string target, int Tp, int knn, - int exclusionRadius, bool verbose ); - + std::list< DataFrame > Project ( std::string lib, + std::string pred, + std::string target, + int Tp, + int knn, + int exclusionRadius, + bool verbose ); }; - #endif diff --git a/src/Embed.cc b/src/Embed.cc index ba1b854..1a7ce76 100644 --- a/src/Embed.cc +++ b/src/Embed.cc @@ -110,6 +110,8 @@ DataFrame< double > MakeBlock( DataFrame< double > dataFrame, std::vector columnNames, bool verbose ) { + if ( verbose ){} + if ( columnNames.size() != dataFrame.NColumns() ) { std::stringstream errMsg; errMsg << "MakeBlock: The number of columns in the dataFrame (" @@ -132,7 +134,7 @@ DataFrame< double > MakeBlock( DataFrame< double > dataFrame, std::vector< std::string > newColumnNames( NColOut ); size_t newCol_i = 0; for ( size_t col = 0; col < columnNames.size(); col ++ ) { - for ( size_t e = 0; e < E; e++ ) { + for ( int e = 0; e < E; e++ ) { std::stringstream ss; if ( tau < 0 ) { ss << columnNames[ col ] << "(t-" << e << ")"; @@ -163,7 +165,7 @@ DataFrame< double > MakeBlock( DataFrame< double > dataFrame, // Shift column data and write to embedding data frame for ( size_t col = 0; col < dataFrame.NColumns(); col++ ) { // for each embedding dimension - for ( size_t e = 0; e < E; e++ ) { + for ( int e = 0; e < E; e++ ) { std::valarray< double > column = dataFrame.Column( col ); diff --git a/src/Parameter.cc b/src/Parameter.cc index 4b1b30c..67d18ca 100644 --- a/src/Parameter.cc +++ b/src/Parameter.cc @@ -31,11 +31,8 @@ Parameters::Parameters( std::string blockFile, std::string derivatives_str, - double svdSig, - double tikhonov, - double elasticNet, - int multi, + bool trainLib, std::string libSizes_str, int sample, bool random, @@ -68,13 +65,10 @@ Parameters::Parameters( SmapOutputFile ( SmapFile ), blockOutputFile ( blockFile ), - derivatives_str ( derivatives_str ), - SVDSignificance ( svdSig ), - TikhonovAlpha ( tikhonov ), - ElasticNetAlpha ( elasticNet ), MultiviewEnsemble( multi ), + MultiviewTrainLib( trainLib ), libSizes_str ( libSizes_str ), subSamples ( sample ), randomLib ( random ), @@ -84,7 +78,7 @@ Parameters::Parameters( // Set validated flag and instantiate Version validated ( false ), - version ( 1, 3, 5, "2020-04-16" ) + version ( 2, 0, 0, "2020-05-30" ) { // Constructor code if ( method != Method::None ) { @@ -363,7 +357,7 @@ void Parameters::Validate() { size_t N_lib = std::floor( (stop-start)/increment + 1/increment ) + 1; - if ( start < E ) { + if ( start < (size_t) E ) { std::stringstream errMsg; errMsg << "Parameters::Validate(): " << "CCM librarySizes start < E = " << E << "\n"; diff --git a/src/Parameter.h b/src/Parameter.h index 60125c3..22c5a87 100644 --- a/src/Parameter.h +++ b/src/Parameter.h @@ -55,7 +55,8 @@ class Parameters { double ElasticNetAlpha; // Initial alpha parameter int MultiviewEnsemble; // Number of ensembles in multiview - + bool MultiviewTrainLib; // Use in sample prediction for k select + std::string libSizes_str; std::vector librarySizes;// CCM library sizes to evaluate int subSamples; // CCM number of samples to draw @@ -94,14 +95,11 @@ class Parameters { bool verbose = false, std::string SmapFile = "", - std::string blockFile = "", + std::string blockFile = "", std::string derivatives_str = "", - double svdSig = 1E-5, - double tikhonov = 0, - double elasticNet = 0.1, - int multi = 0, + bool trainLib = true, std::string libSizes_str = "", int sample = 0, bool random = true, diff --git a/src/Simplex.cc b/src/Simplex.cc index 16fad86..bb3d503 100644 --- a/src/Simplex.cc +++ b/src/Simplex.cc @@ -6,15 +6,26 @@ // Simplex() : Constructor // See EDM class for descriptions of parameters not described //---------------------------------------------------------------- -SimplexMachine::SimplexMachine ( DataFrame & data, - std::string pathOut, std::string predictFile, - std::string lib, std::string pred, - int E, int Tp, int knn, int tau, int exclusionRadius, - std::string columns, std::string target, - bool embedded, bool const_predict, bool verbose ): - EDM{data, E, tau, columns, target, embedded, verbose } { - - Project(lib,pred,target,Tp,knn,exclusionRadius,verbose); - - +SimplexMachine::SimplexMachine ( + DataFrame & data, + std::string pathOut, + std::string predictFile, + std::string lib, + std::string pred, + int E, + int Tp, + int knn, + int tau, + int exclusionRadius, + std::string columns, + std::string target, + bool embedded, + bool const_predict, + bool verbose ): + EDM{ data, E, tau, embedded, columns, target, true, verbose } { + + Project( lib, pred, target, Tp, knn, exclusionRadius, verbose ); + if ( pathOut.size() ) {} + if ( predictFile.size() ) {} + if ( const_predict ) {} } diff --git a/src/makefile b/src/makefile index a18f40b..00e8c82 100644 --- a/src/makefile +++ b/src/makefile @@ -1,14 +1,21 @@ +.PHONY: all clean distclean depend CC = g++ -OBJ = Common.o Embed.o EDM.o Simplex.o EDM_Functions.o Parameter.o EDM_Helpers.o +HEADERS = Common.h DataFrame.h EDM.h Embed.h Parameter.h Simplex.h Version.h + +SRCS = Common.cc EDM.cc EDM_Functions.cc EDM_Helpers.cc Embed.cc\ + Parameter.cc Simplex.cc + +OBJ = $(SRCS:%.cc=%.o) LIB = libEDM.a -CFLAGS = -std=c++11 -O3 -Wreorder -M -g -DDEBUG -DDEBUG_ALL -LFLAGS = -L./ -lstdc++ -lEDM -lpthread +CFLAGS += -std=c++11 -DMULTIVIEW_VALUES_OVERLOAD -O3\ + -Wpedantic -Wall -Wextra -Wreorder # -g -DDEBUG -DDEBUG_ALL + +LFLAGS = -L./ -lstdc++ -lEDM -lpthread # -llapacke -llapack -lblas all: $(LIB) - ar -rcs $(LIB) $(OBJ) cp $(LIB) ../lib/ clean: @@ -18,20 +25,20 @@ distclean: rm -f $(OBJ) $(LIB) ../lib/$(LIB) *~ *.bak *.csv $(LIB): $(OBJ) + ar -rcs $(LIB) $(OBJ) -%.o : %.cc %.h - $(CC) -c $< -o $@ +%.o : %.cc + $(CC) $(CFLAGS) -c $< -SRCS = `echo ${OBJ} | sed -e 's/.o /.cc /g'` depend: @echo ${SRCS} - makedepend -Y $(SRCS) + makedepend -Y $(SRCS) -w160 # DO NOT DELETE Common.o: Common.h DataFrame.h -Parameter.o: Common.o -Embed.o: Parameter.o -EDM_Helpers.o: Common.o -EDM.o: Embed.o EDM_Helpers.o -Simplex.o: EDM.o -EDM_Functions.o: Simplex.o +EDM.o: EDM.h Common.h DataFrame.h Embed.h Parameter.h Version.h EDM_Helpers.cc +EDM_Functions.o: Common.h DataFrame.h Simplex.h EDM.h Embed.h Parameter.h Version.h +EDM_Helpers.o: Common.h DataFrame.h +Embed.o: Embed.h Common.h DataFrame.h Parameter.h Version.h +Parameter.o: Parameter.h Common.h DataFrame.h Version.h +Simplex.o: Simplex.h Common.h DataFrame.h EDM.h Embed.h Parameter.h Version.h