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.h b/src/Common.h index 3d8edd6..334466f 100644 --- a/src/Common.h +++ b/src/Common.h @@ -1,385 +1,56 @@ +// File for declarations common to all EDM code + #ifndef COMMON_H #define COMMON_H #include +#include #include #include #include #include -#include +#include #include #include #ifdef _MSC_VER -#include // macro constants for MSVC C++ operators not in ISO646 +#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 }; +#include "DataFrame.h" // Has #include Common.h -//--------------------------------------------------------- // 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 ) {} -#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 ); +// Enumerations +enum class Method { None, Embed, Simplex, SMap, CCM }; +enum class DistanceMetric { Euclidean, Manhattan }; -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 ); +// Namespace with all functions user should have +namespace cppEDM { -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 ); + // Forward declarations of EDM Functions -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 ); + 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 ) ; -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 ); + 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 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 ); -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..93ac486 --- /dev/null +++ b/src/EDM.cc @@ -0,0 +1,360 @@ +// Contains class definitions for EDM object. See EDM.h header for +// 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() + // DBL_MAX is a Macro equivalent to: std::numeric_limits::max() + double DistanceMax = std::numeric_limits::max(); + double DistanceLimit = std::numeric_limits::max() - 1; +} + +#ifdef DEBUG_ALL +void PrintNeighbors( const EDM::Neighbors &neighbors ); +#endif + +//---------------------------------------------------------------- +// 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 +// embedded : Whether data supplied is already embedded or not +// verbose : Verbose information flag +// +//---------------------------------------------------------------- +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(): " + "tau must be non-zero.\n" ); + throw std::runtime_error( errMsg ); + } + + std::vector columnNames = ParseColumnNames( columns ); + std::vector columnIndex = ParseColumnIndices( columns ); + + 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() ); + } + + 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() : Project +// Finds neighbors, performs weighting on neighbors, and +// projects onto pred range. +// +// 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 +// 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 ) { + + // Parse lib and pred range strings + std::vector libIndices = ParseRangeStr( lib ); + std::vector predIndices = ParseRangeStr( pred ); + + 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; + + // 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] ); + char tauSign = tau > 0 ? '+' : '-'; + targName = originalName+ "(t" + tauSign + "0)"; + } + + targetIn = embedding.VectorColumnName( targName ); + } + // Default to first column + else { + targetIn = data.Column( 0 ); + } + + // Shift lib and pred by number of partial rows if we embedded + + if ( not embedded ) { + + size_t shift = abs( tau ) * ( E - 1 ); + + if ( shift > 0 ) { + libIndices = AdjustRange( shift, tau, libIndices); + predIndices = AdjustRange( shift, tau, predIndices); + } + + // Check boundaries again since rows were potentially removed + CheckDataRows( data.NRows(), libIndices.back(), predIndices.back(), + E, tau, embedded ); + } + + // Find neighbors + FindNeighbors( libIndices, predIndices, Tp, knn, exclusionRadius, verbose ); + + // Weight neighbors + + // Project onto target for prediction + + return std::list>(); +} + +//---------------------------------------------------------------- +// EDM() : FindNeighbors +// 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::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 + // 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 * (E + 1); + } + else if ( knn == 0 ) { + + knn = E + 1; + + if ( verbose ) { + std::stringstream msg; + msg << "Parameters::Validate(): Set knn = " << knn + << " " << std::endl; + std::cout << msg.str(); + } + } + else if ( knn < E + 1 ) { + std::stringstream errMsg; + errMsg << "Parameters::Validate(): knn of " << knn + << " 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); + 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 ); + + // 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 = embedding.Row( pred_row ); + + // Reset the neighbor and distance vectors for this pred row + for ( int 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 = embedding.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 ) != + (long int) 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++ ) + + #ifdef DEBUG_ALL + const Neighbors &neigh = neighbors; + PrintNeighbors( neigh ); + #endif + + return neighbors; + +} + + +//---------------------------------------------------------------- +// Debug method to print out all neighbors +//---------------------------------------------------------------- +#ifdef DEBUG_ALL +void PrintNeighbors( const EDM::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/EDM.h b/src/EDM.h new file mode 100644 index 0000000..7e3c14a --- /dev/null +++ b/src/EDM.h @@ -0,0 +1,79 @@ +// 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" +#include "Embed.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 and its embedding + DataFrame & data; + DataFrame embedding; + + // The current embedding dimension (meaning may change when implementing + // EmbedDimension) and 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 + + bool noNeighborLimit; + bool verbose; + + public: + + // Return structure of FindNeighbors() + struct Neighbors { + DataFrame neighbors; + DataFrame distances; + }; + + //---------------------------------------------------------------- + // EDM() : Constructor + //---------------------------------------------------------------- + EDM ( DataFrame & data, + int E, + int tau, + bool embedded, + std::string columns, + std::string target, + bool noNeighborLimit, + bool verbose ); + + //---------------------------------------------------------------- + // EDM() : FindNeighbors + //---------------------------------------------------------------- + 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 ); +}; +#endif diff --git a/src/EDM_Functions.cc b/src/EDM_Functions.cc new file mode 100644 index 0000000..21546a0 --- /dev/null +++ b/src/EDM_Functions.cc @@ -0,0 +1,49 @@ +// 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" +#include "Simplex.h" + +namespace cppEDM { + + //---------------------------------------------------------------- + // 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 + + SimplexMachine simplexManager( data, pathOut, predictFile, + lib, pred, E, Tp, knn, tau, exclusionRadius, columns, target, + embedded, const_predict, verbose ); + + return DataFrame(); + + } + + //---------------------------------------------------------------- + // 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/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/Embed.cc b/src/Embed.cc index c1fba55..1a7ce76 100644 --- a/src/Embed.cc +++ b/src/Embed.cc @@ -1,7 +1,6 @@ -#include "Embed.h" - // NOTE: The returned data block does NOT have the time column +#include "Embed.h" //---------------------------------------------------------------- // API Overload 1: Explicit data file path/name @@ -111,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 (" @@ -133,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 << ")"; @@ -164,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/Embed.h b/src/Embed.h index 3bccba0..dd08100 100644 --- a/src/Embed.h +++ b/src/Embed.h @@ -34,3 +34,4 @@ DataFrame< double > MakeBlock ( DataFrame< double > dataFrame, 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 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/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..bb3d503 100644 --- a/src/Simplex.cc +++ b/src/Simplex.cc @@ -1,244 +1,31 @@ +// Definitions for the Simplex class -#include "Common.h" -#include "Parameter.h" -#include "Neighbors.h" -#include "Embed.h" -#include "AuxFunc.h" - -// Forward declaration -DataFrame SimplexProjection( Parameters param, - DataEmbedNN embedNN, - bool checkDataRows = true ); +#include "Simplex.h" //---------------------------------------------------------------- -// API Overload 1: Explicit data file path/name -// Implemented as a wrapper to API Overload 2: +// Simplex() : Constructor +// See EDM class for descriptions of parameters not described //---------------------------------------------------------------- -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; +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 } { - return S; -} - -//---------------------------------------------------------------- -// 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 - - return dataFrame; + Project( lib, pred, target, Tp, knn, exclusionRadius, verbose ); + if ( pathOut.size() ) {} + if ( predictFile.size() ) {} + if ( const_predict ) {} } diff --git a/src/Simplex.h b/src/Simplex.h new file mode 100644 index 0000000..1ade394 --- /dev/null +++ b/src/Simplex.h @@ -0,0 +1,32 @@ +// Contains the Simplex class declarations +// 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 SimplexMachine : public EDM { + + public: + + //---------------------------------------------------------------- + // Simplex() : Constructor with DataFrame input + //---------------------------------------------------------------- + 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/makefile b/src/makefile index b22079f..00e8c82 100644 --- a/src/makefile +++ b/src/makefile @@ -1,15 +1,21 @@ +.PHONY: all clean distclean depend 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 +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 -DCCM_THREADED -DMULTIVIEW_VALUES_OVERLOAD -O3 -Wreorder # -g -DDEBUG -DDEBUG_ALL +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: @@ -19,61 +25,20 @@ distclean: rm -f $(OBJ) $(LIB) ../lib/$(LIB) *~ *.bak *.csv $(LIB): $(OBJ) + ar -rcs $(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 + $(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 -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 +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 -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 +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 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