diff --git a/src/CCM.cc b/src/CCM.cc index e074858..7d28c14 100644 --- a/src/CCM.cc +++ b/src/CCM.cc @@ -5,7 +5,7 @@ #include #include -#ifdef CCM_THREADED // Defined in makefile +#ifdef CCM_THREADED // Two explicit CrossMap() threads are invoked. One for forward mapping, // one for inverse mapping. The call signature of CrossMap() is // dependent on which path is used. This should probably be unified @@ -134,6 +134,7 @@ DataFrame CCM( DataFrame< double > dataFrameIn, tau, // 0, // theta 0, // exclusionRadius + DataFrame(),// exclusionMatrix columns, // target, // false, // embedded diff --git a/src/Common.h b/src/Common.h index 7dbf0ff..9801781 100644 --- a/src/Common.h +++ b/src/Common.h @@ -99,6 +99,7 @@ DataFrame Simplex( std::string pathIn = "./data/", int knn = 0, int tau = 1, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), std::string colNames = "", std::string targetName = "", bool embedded = false, @@ -115,6 +116,7 @@ DataFrame Simplex( DataFrame< double > &dataFrameIn, int knn = 0, int tau = 1, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), std::string colNames = "", std::string targetName = "", bool embedded = false, @@ -133,6 +135,7 @@ SMapValues SMap( std::string pathIn = "./data/", int tau = 1, double theta = 0, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), std::string columns = "", std::string target = "", std::string smapFile = "", @@ -152,6 +155,7 @@ SMapValues SMap( DataFrame< double > &dataFrameIn, int tau = 1, double theta = 0, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), std::string columns = "", std::string target = "", std::string smapFile = "", @@ -207,6 +211,7 @@ MultiviewValues Multiview( std::string pathIn = "./", std::string target = "", int multiview = 0, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), bool verbose = false, unsigned nThreads = 4 ); @@ -223,6 +228,7 @@ MultiviewValues Multiview( DataFrame< double >, std::string target = "", int multiview = 0, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), bool verbose = false, unsigned nThreads = 4 ); diff --git a/src/DataFrame.h b/src/DataFrame.h index 6eb8ab0..62f22e5 100644 --- a/src/DataFrame.h +++ b/src/DataFrame.h @@ -62,7 +62,7 @@ class DataFrame { //----------------------------------------------------------------- // Constructors //----------------------------------------------------------------- - DataFrame() {} + DataFrame(): DataFrame(0,0) {} //----------------------------------------------------------------- // Load data from CSV file path/fileName, populate DataFrame diff --git a/src/Embed.cc b/src/Embed.cc index 48e9fab..3412c60 100644 --- a/src/Embed.cc +++ b/src/Embed.cc @@ -39,8 +39,8 @@ DataFrame< double > Embed( DataFrame< double > dataFrameIn, // Parameter.Validate will convert columns into a vector of names // or a vector of column indices Parameters param = Parameters( Method::Embed, "", "", "", "", - "1 1", "1 1", E, 0, 0, tau, 0, 0, - columns, "", false, false, verbose ); + "1 1", "1 1", E, 0, 0, tau, 0, 0, DataFrame(), + columns, "", false, false, verbose ); if ( not param.columnIndex.size() and dataFrameIn.ColumnNameToIndex().empty() ) { diff --git a/src/Eval.cc b/src/Eval.cc index d974463..bdee18a 100644 --- a/src/Eval.cc +++ b/src/Eval.cc @@ -237,6 +237,7 @@ void EmbedThread( EDM_Eval::WorkQueue &workQ, 0, // knn tau, 0, // exclusionRadius + DataFrame(),// exclusionMatrix colNames, targetName, embedded, @@ -425,6 +426,7 @@ void PredictIntervalThread( EDM_Eval::WorkQueue &workQ, 0, // knn tau, 0, // exclusionRadius + DataFrame(),// exclusionMatrix colNames, targetName, embedded, @@ -649,6 +651,7 @@ void SMapThread( EDM_Eval::WorkQueue &workQ, tau, theta, 0, // exclusionRadius + DataFrame(),// exclusionMatrix colNames, targetName, "", // smapFile diff --git a/src/Multiview.cc b/src/Multiview.cc index a163f18..34c4d3a 100644 --- a/src/Multiview.cc +++ b/src/Multiview.cc @@ -82,6 +82,7 @@ MultiviewValues Multiview( std::string pathIn, std::string target, int multiview, int exclusionRadius, + const DataFrame & exclusionMatrix, bool verbose, unsigned nThreads ) { @@ -101,6 +102,7 @@ MultiviewValues Multiview( std::string pathIn, target, multiview, exclusionRadius, + exclusionMatrix, verbose, nThreads ); return result; @@ -123,6 +125,7 @@ MultiviewValues Multiview( DataFrame< double > data, std::string target, int multiview, int exclusionRadius, + const DataFrame & exclusionMatrix, bool verbose, unsigned nThreads ) { @@ -130,7 +133,7 @@ MultiviewValues Multiview( DataFrame< double > data, Parameters param = Parameters( Method::Simplex, "", "", pathOut, predictFile, lib, pred, E, Tp, knn, tau, 0, - exclusionRadius, columns, target, + exclusionRadius, exclusionMatrix, columns, target, true, false, verbose, "", "", "", 0, 0, 0, multiview ); diff --git a/src/Neighbors.cc b/src/Neighbors.cc index 5727b97..906bfe1 100644 --- a/src/Neighbors.cc +++ b/src/Neighbors.cc @@ -92,26 +92,45 @@ Neighbors FindNeighbors( //-------------------------------------------------------------- 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; } + //skip excluded neighbors in exclusion matrix if provided + if ( parameters.exclusionMatrix.NRows() and + parameters.method != Method::CCM) { + + bool exclude_row = false; + + //if neighb idx OR any idx's in its lag are excluded, skip + + int lead_time_idx = (parameters.E-1)*parameters.tau + lib_row; + int pred_time_row = pred_row + parameters.tau*(parameters.E-1); + + for ( int curr_row = lib_row; curr_row <= lead_time_idx; + curr_row += parameters.tau ) { + if( parameters.exclusionMatrix( pred_time_row , curr_row ) ){ + exclude_row = true; + } + } + + if ( exclude_row ) { + 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; } diff --git a/src/Parameter.cc b/src/Parameter.cc index 0735688..6ec89e6 100644 --- a/src/Parameter.cc +++ b/src/Parameter.cc @@ -20,6 +20,8 @@ Parameters::Parameters( float theta, int exclusionRadius, + const DataFrame &exclusionMatrix, + std::string columns_str, std::string target_str, @@ -58,6 +60,8 @@ Parameters::Parameters( tau ( tau ), theta ( theta ), exclusionRadius ( exclusionRadius ), + + exclusionMatrix ( exclusionMatrix ), columns_str ( columns_str ), target_str ( target_str ), @@ -142,6 +146,7 @@ void Parameters::Validate() { library = std::vector( lib_end - lib_start + 1 ); std::iota ( library.begin(), library.end(), lib_start - 1 ); + } //-------------------------------------------------------------- @@ -169,6 +174,19 @@ void Parameters::Validate() { prediction = std::vector( pred_end - pred_start + 1 ); std::iota ( prediction.begin(), prediction.end(), pred_start - 1 ); + + //also check exclusion matrix size while we have pred string parsed + + if ( exclusionMatrix.NRows() and method != Method::CCM and ( + exclusionMatrix.NRows() < pred_end or + exclusionMatrix.NColumns() < pred_end ) ){ + + std::string errMsg( "Parameters::Validate(): " + "The range of rows in the Exclusion Matrix " + "is smaller than range predicting on.\n" ); + throw std::runtime_error( errMsg ); + } + } if ( method == Method::Simplex or method == Method::SMap ) { @@ -193,6 +211,7 @@ void Parameters::Validate() { } } + #ifdef DEBUG_ALL PrintIndices( library, prediction ); #endif @@ -253,6 +272,7 @@ void Parameters::Validate() { } } + // CCM librarySizes if ( libSizes_str.size() > 0 ) { std::vector libsize_vec = SplitString(libSizes_str," \t,"); diff --git a/src/Parameter.h b/src/Parameter.h index ab2aca5..3c3e209 100644 --- a/src/Parameter.h +++ b/src/Parameter.h @@ -21,6 +21,8 @@ class Parameters { int knn; // k nearest neighbors int tau; // block embedding delay int exclusionRadius; // temporal rows to ignore in predict + + const DataFrame &exclusionMatrix; float theta; // S Map localization float SVDSignificance; // SVD singular value cutoff @@ -85,6 +87,8 @@ class Parameters { float theta = 0, int exclusionRadius = 0, + const DataFrame &exclusionMatrix = DataFrame(), + std::string columns_str = "", std::string target_str = "", diff --git a/src/SMap.cc b/src/SMap.cc index 51ef93c..fa423af 100644 --- a/src/SMap.cc +++ b/src/SMap.cc @@ -31,6 +31,7 @@ SMapValues SMap( std::string pathIn, int tau, double theta, int exclusionRadius, + const DataFrame & exclusionMatrix, std::string columns, std::string target, std::string smapFile, @@ -44,7 +45,7 @@ SMapValues SMap( std::string pathIn, SMapValues SMapOutput = SMap( dataFrameIn, pathOut, predictFile, lib, pred, E, Tp, knn, tau, theta, - exclusionRadius, + exclusionRadius, exclusionMatrix, columns, target, smapFile, derivatives, embedded, const_predict, verbose ); return SMapOutput; @@ -64,6 +65,7 @@ SMapValues SMap( DataFrame< double > &data, int tau, double theta, int exclusionRadius, + const DataFrame & exclusionMatrix, std::string columns, std::string target, std::string smapFile, @@ -76,7 +78,7 @@ SMapValues SMap( DataFrame< double > &data, Parameters param = Parameters( Method::SMap, "", "", pathOut, predictFile, lib, pred, E, Tp, knn, tau, theta, - exclusionRadius, columns, target, + exclusionRadius, exclusionMatrix, columns, target, embedded, const_predict, verbose, smapFile, "", derivatives ); diff --git a/src/Simplex.cc b/src/Simplex.cc index f6a936a..eb705fb 100644 --- a/src/Simplex.cc +++ b/src/Simplex.cc @@ -25,6 +25,7 @@ DataFrame Simplex( std::string pathIn, int knn, int tau, int exclusionRadius, + const DataFrame & exclusionMatrix, std::string columns, std::string target, bool embedded, @@ -46,6 +47,7 @@ DataFrame Simplex( std::string pathIn, knn, tau, exclusionRadius, + exclusionMatrix, columns, target, embedded, @@ -69,6 +71,7 @@ DataFrame Simplex( DataFrame< double > &data, int knn, int tau, int exclusionRadius, + const DataFrame &exclusionMatrix, std::string columns, std::string target, bool embedded, @@ -78,7 +81,7 @@ DataFrame Simplex( DataFrame< double > &data, Parameters param = Parameters( Method::Simplex, "", "", pathOut, predictFile, lib, pred, E, Tp, knn, tau, 0, - exclusionRadius, + exclusionRadius, exclusionMatrix, columns, target, embedded, const_predict, verbose ); @@ -134,12 +137,13 @@ DataFrame SimplexProjection( Parameters param, 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 @@ -220,7 +224,6 @@ DataFrame SimplexProjection( Parameters param, } #ifdef DEBUG_ALL - std::cout << dataFrame; VectorError ve = ComputeError( dataFrame.VectorColumnName( "Observations" ), dataFrame.VectorColumnName( "Predictions" ) ); diff --git a/tests/CustomExclusionTest.cc b/tests/CustomExclusionTest.cc new file mode 100644 index 0000000..22da02c --- /dev/null +++ b/tests/CustomExclusionTest.cc @@ -0,0 +1,29 @@ +// custom exclusion data test + +#include "TestCommon.h" +#include +#include "Embed.cc" + +using namespace std; + +int main () { + + DataFrame < double > twospike_data( "data/", "twospike_data.csv" ); + + DataFrame< double > exclusion_mat( twospike_data.NRows(), + twospike_data.NRows() ); + + // want to exclude indices 20-32 when predicting indices 72-84 (non-1 vals) + + for ( int j = 72; j < 83; j++ ) + for ( int i = 20; i < 32; i++ ) + exclusion_mat(j,i) = 1; + + // try predicting on data now without non-1 points + + cout<< Simplex( twospike_data, "","","1 53","75 90",3,1,4,3, + 0, exclusion_mat, + "wave","wave", false,false,false); + + +} diff --git a/tests/MultiviewTest.cc b/tests/MultiviewTest.cc index 5a013a9..9d41138 100644 --- a/tests/MultiviewTest.cc +++ b/tests/MultiviewTest.cc @@ -29,6 +29,7 @@ int main( int argc, char *argv[] ) { "x_t", // target, 0, // multiview 0, // exclusionRadius + DataFrame(), false, // verbose, 1 ); // nThreads diff --git a/tests/SMapTest.cc b/tests/SMapTest.cc index 30481ff..e1b8d22 100644 --- a/tests/SMapTest.cc +++ b/tests/SMapTest.cc @@ -19,7 +19,7 @@ int main () { SMapValues smapVals = SMap ( circleData, "./data/", "Smap_circle_cppEDM.csv", - " 1 100 ", "101 198", 2, 1, 0, 1, 4, 0, + " 1 100 ", "101 198", 2, 1, 0, 1, 4, 0,DataFrame(), "x y", "x", "", "", true, false, false ); DataFrame < double > cppOutput = smapVals.predictions; @@ -37,7 +37,7 @@ int main () { smapVals = SMap ( "../data/", "block_3sp.csv", "./data/", "Smap_embd_block_3sp_cppEDM.csv", - " 1 99 ", "100 198", 3, 1, 0, 1, 2, 0, + " 1 99 ", "100 198", 3, 1, 0, 1, 2, 0,DataFrame(), "x_t y_t z_t", "x_t", "", "", true, false, false ); cppOutput = smapVals.predictions; diff --git a/tests/SimplexTest.cc b/tests/SimplexTest.cc index bb5b7c3..91cd39b 100644 --- a/tests/SimplexTest.cc +++ b/tests/SimplexTest.cc @@ -19,7 +19,7 @@ int main () { // Generate cpp output cppOutput = Simplex ( "../data/", "block_3sp.csv", "./data/", "Smplx_embd_block_3sp_cppEDM.csv", - "1 99","100 198", 3, 1, 0, 1, 0, + "1 99","100 198", 3, 1, 0, 1, 0, DataFrame(), "x_t y_t z_t", "x_t", true, false, false ); // Comparison MakeTest ( "block_3sp.csv embedded data test", pyOutput, cppOutput ); @@ -35,7 +35,7 @@ int main () { // Generate cpp output cppOutput = Simplex ( "../data/", "block_3sp.csv", "./data/", "Smplx_E3_block_3sp_cppEDM.csv", - "1 100", "101 195", 3, 1, 0, 1, 0, + "1 100", "101 195", 3, 1, 0, 1, 0,DataFrame(), "x_t", "x_t", false, false, false ); // Comparison MakeTest ( "block_3sp.csv dynamic embedding test", pyOutput, cppOutput ); @@ -53,7 +53,7 @@ int main () { // Generate cpp output cppOutput = Simplex ( "../data/", "S12CD-S333-SumFlow_1980-2005.csv", "./data/", "Smplx_S12CD_E3_cppEDM.csv", - "1 800", "801 1375", 3, 1, 0, 1, 0, + "1 800", "801 1375", 3, 1, 0, 1, 0,DataFrame(), "S12.C.D.S333", "S12.C.D.S333", false, false, false ); cppOutput.MaxRowPrint() = 5; // Set number of rows to print diff --git a/tests/makefile b/tests/makefile index 7a3c247..37e4e02 100644 --- a/tests/makefile +++ b/tests/makefile @@ -1,19 +1,24 @@ CC = g++ -EXE = SimplexTest TestCommonTest SMapTest CCMTest MultiviewTest DateTimeTest +EXE = SimplexTest TestCommonTest SMapTest CCMTest MultiviewTest DateTimeTest \ + CustomExclusionTest OBJ = $(EXE:=.o) TestCommon.o CFLAGS = -std=c++11 -D PRINT_DIFFERENCE_IN_RESULTS LFLAGS = -lstdc++ -L../lib/ -I../src/ -lEDM -lpthread -llapack -all: $(EXE) +all: $(EXE) # Need to figure out how to build TestCommon.o automatically SimplexTest: SimplexTest.cc $(CC) TestCommon.cc -c $(CFLAGS) $(LFLAGS) $(CC) $@.cc -o $@ $(CFLAGS) $(LFLAGS) TestCommon.o +CustomExclusionTest: CustomExclusionTest.cc + $(CC) CustomExclusionTest.cc -c $(CFLAGS) $(LFLAGS) + $(CC) $@.cc -o $@ $(CFLAGS) $(LFLAGS) TestCommon.o + DateTimeTest: DateTimeTest.cc ../src/DateTimeUtil.cc $(CC) DateTimeTest.cc -c $(CFLAGS) $(LFLAGS) $(CC) $@.cc -o $@ $(CFLAGS) $(LFLAGS) TestCommon.o