From 677f22cfd4aaddd49421754feb25263397358141 Mon Sep 17 00:00:00 2001 From: James Brown <64858662+james-d-brown@users.noreply.github.com> Date: Mon, 23 Dec 2024 19:38:19 +0000 Subject: [PATCH] Add an event detection method based on Regina and Ogden (2021), #130 --- build.gradle | 1 + .../pipeline/pooling/EventsGenerator.java | 15 +- src/wres/pipeline/pooling/PoolFactory.java | 33 +- wres-config/nonsrc/schema.yml | 27 + .../config/yaml/DeclarationValidator.java | 54 ++ .../yaml/components/EventDetection.java | 22 +- .../components/EventDetectionParameters.java | 44 + .../config/yaml/components/FeatureGroups.java | 1 + .../deserializers/DurationDeserializer.java | 16 +- .../DurationIntervalDeserializer.java | 8 +- .../EventDetectionDeserializer.java | 46 +- .../deserializers/TimeScaleDeserializer.java | 2 +- .../config/yaml/DeclarationFactoryTest.java | 12 +- .../config/yaml/DeclarationValidatorTest.java | 6 + .../wres/datamodel/time/TimeSeriesSlicer.java | 14 +- .../EventDetectionException.java | 5 +- .../wres/eventdetection/EventDetector.java | 6 +- .../eventdetection/EventDetectorFactory.java | 9 +- .../ReginaOgdenEventDetector.java | 773 +++++++++++++++++- .../ReginaOgdenEventDetectorTest.java | 225 +++++ 20 files changed, 1266 insertions(+), 53 deletions(-) create mode 100644 wres-config/src/wres/config/yaml/components/EventDetectionParameters.java create mode 100644 wres-eventdetection/test/wres/eventdetection/ReginaOgdenEventDetectorTest.java diff --git a/build.gradle b/build.gradle index 241277b7da..b5e24d4068 100755 --- a/build.gradle +++ b/build.gradle @@ -660,6 +660,7 @@ project(':wres-eventdetection') { implementation project(':wres-datamodel') api project(':wres-statistics') implementation 'org.slf4j:slf4j-api:2.0.13' + implementation group: 'org.apache.commons', name: 'commons-math3', version: '3.6.1' runtimeOnly('ch.qos.logback:logback-classic:1.5.11') { // Not used at runtime, bloat diff --git a/src/wres/pipeline/pooling/EventsGenerator.java b/src/wres/pipeline/pooling/EventsGenerator.java index d42dcf4c62..6f5a01d6d8 100644 --- a/src/wres/pipeline/pooling/EventsGenerator.java +++ b/src/wres/pipeline/pooling/EventsGenerator.java @@ -17,7 +17,6 @@ import wres.config.yaml.components.EvaluationDeclaration; import wres.config.yaml.components.EventDetection; import wres.config.yaml.components.EventDetectionDataset; -import wres.config.yaml.components.EventDetectionMethod; import wres.datamodel.scale.TimeScaleOuter; import wres.datamodel.space.Feature; import wres.datamodel.space.FeatureGroup; @@ -27,7 +26,6 @@ import wres.datamodel.time.TimeWindowOuter; import wres.datamodel.time.TimeSeries; import wres.eventdetection.EventDetector; -import wres.eventdetection.EventDetectorFactory; import wres.io.project.Project; import wres.io.retrieving.RetrieverFactory; import wres.statistics.generated.TimeScale; @@ -40,14 +38,17 @@ * @param baselineUpscaler the upscaler for single-valued time-series with a baseline orientation * @param covariateUpscaler the upscaler for single-valued time-series with a covariate orientation * @param measurementUnit the measurement unit + * @param eventDetector the event detector * @author James Brown */ record EventsGenerator( TimeSeriesUpscaler leftUpscaler, TimeSeriesUpscaler rightUpscaler, TimeSeriesUpscaler baselineUpscaler, TimeSeriesUpscaler covariateUpscaler, - String measurementUnit ) + String measurementUnit, + EventDetector eventDetector ) { + /** * Construct and validate. * @@ -56,6 +57,7 @@ record EventsGenerator( TimeSeriesUpscaler leftUpscaler, * @param baselineUpscaler the upscaler for single-valued time-series with a baseline orientation * @param covariateUpscaler the upscaler for single-valued time-series with a covariate orientation * @param measurementUnit the measurement unit + * @param eventDetector the event detector */ EventsGenerator { @@ -64,6 +66,7 @@ record EventsGenerator( TimeSeriesUpscaler leftUpscaler, Objects.requireNonNull( baselineUpscaler ); Objects.requireNonNull( covariateUpscaler ); Objects.requireNonNull( measurementUnit ); + Objects.requireNonNull( eventDetector ); } /** Logger. */ @@ -356,11 +359,9 @@ private Set doEventDetection( TimeSeries timeSeries, + "metadata: {}.", timeSeries.getEvents() .size(), timeSeries.getMetadata() ); - // Get an event detector - EventDetector eventDetector = EventDetectorFactory.getEventDetector( EventDetectionMethod.DEFAULT ); - // Unbounded time window, placeholder - return eventDetector.detect( timeSeries, details.detection() ); + return this.eventDetector() + .detect( timeSeries ); } /** diff --git a/src/wres/pipeline/pooling/PoolFactory.java b/src/wres/pipeline/pooling/PoolFactory.java index 732d48b88f..bbdbfc6c29 100644 --- a/src/wres/pipeline/pooling/PoolFactory.java +++ b/src/wres/pipeline/pooling/PoolFactory.java @@ -78,6 +78,8 @@ import wres.datamodel.time.TimeSeriesUpscaler; import wres.datamodel.time.TimeWindowOuter; import wres.datamodel.baselines.PersistenceGenerator; +import wres.eventdetection.EventDetector; +import wres.eventdetection.EventDetectorFactory; import wres.io.project.Project; import wres.io.retrieving.CachingSupplier; import wres.io.retrieving.RetrieverFactory; @@ -2938,11 +2940,30 @@ private PoolFactory( Project project ) this.baselineEnsembleUpscaler = TimeSeriesOfEnsembleUpscaler.of( baselineLenient, this.getUnitMapper() .getUnitAliases() ); - this.eventsGenerator = new EventsGenerator( this.getLeftSingleValuedUpscaler(), - this.getRightSingleValuedUpscaler(), - this.getBaselineSingleValuedUpscaler(), - this.getCovariateSingleValuedUpscaler(), - this.getUnitMapper() - .getDesiredMeasurementUnitName() ); + + if ( Objects.nonNull( this.project.getDeclaration() + .eventDetection() ) ) + { + EventDetector eventDetector = EventDetectorFactory.getEventDetector( this.getProject() + .getDeclaration() + .eventDetection() + .method(), + this.getProject() + .getDeclaration() + .eventDetection() + .parameters() ); + + this.eventsGenerator = new EventsGenerator( this.getLeftSingleValuedUpscaler(), + this.getRightSingleValuedUpscaler(), + this.getBaselineSingleValuedUpscaler(), + this.getCovariateSingleValuedUpscaler(), + this.getUnitMapper() + .getDesiredMeasurementUnitName(), + eventDetector ); + } + else + { + this.eventsGenerator = null; + } } } \ No newline at end of file diff --git a/wres-config/nonsrc/schema.yml b/wres-config/nonsrc/schema.yml index 4dac11a53b..63abef4763 100644 --- a/wres-config/nonsrc/schema.yml +++ b/wres-config/nonsrc/schema.yml @@ -591,9 +591,36 @@ definitions: minItems: 1 method: - "$ref": "#/definitions/EventDetectionMethodEnum" + parameters: + anyOf: + - "$ref": "#/definitions/EventDetectionParametersReginaOgden" required: - dataset + EventDetectionParametersReginaOgden: + title: Time-series event detection parameters for the Regina-Ogden method + description: "Event detection parameters for the Regina-Ogden event + detection method." + type: object + additionalProperties: false + properties: + window_size: + type: integer + minimum: 1 + minimum_event_duration: + type: integer + minimum: 0 + half_life: + type: integer + minimum: 1 + start_radius: + type: integer + minimum: 0 + duration_unit: + "$ref": "#/definitions/DurationEnum" + required: + - duration_unit + CrossPair: title: Cross-pairing of time-series for consistency description: "Applies cross-pairing to the time-series within an evaluation diff --git a/wres-config/src/wres/config/yaml/DeclarationValidator.java b/wres-config/src/wres/config/yaml/DeclarationValidator.java index cf21a8b1f8..8bde9bd364 100644 --- a/wres-config/src/wres/config/yaml/DeclarationValidator.java +++ b/wres-config/src/wres/config/yaml/DeclarationValidator.java @@ -2399,6 +2399,60 @@ private static List eventDetectionIsValid( EvaluationDecl events.add( warn ); } + // Check parameters + List parameters = DeclarationValidator.eventDetectionParametersAreValid( declaration ); + events.addAll( parameters ); + + return Collections.unmodifiableList( events ); + } + + /** + * Checks that any event detection parameters are valid. + * @param declaration the evaluation declaration + * @return the validation events encountered + */ + private static List eventDetectionParametersAreValid( EvaluationDeclaration declaration ) + { + List events = new ArrayList<>(); + + // Parameters undefined for which estimates/defaults are speculative: warn + if ( Objects.isNull( declaration.eventDetection() + .parameters() ) + || Objects.isNull( declaration.eventDetection() + .parameters() + .windowSize() ) ) + { + EvaluationStatusEvent warn + = EvaluationStatusEvent.newBuilder() + .setStatusLevel( StatusLevel.WARN ) + .setEventMessage( "Event detection was declared, but the window size " + + "parameter was undefined. An attempt will be made to " + + "choose a reasonable default by inspecting the " + + "time-series data, but it is strongly recommended that " + + "you instead declare the 'window_size' explicitly as " + + "the default value may not be appropriate." ) + .build(); + events.add( warn ); + } + if ( Objects.isNull( declaration.eventDetection() + .parameters() ) + || Objects.isNull( declaration.eventDetection() + .parameters() + .halfLife() ) ) + { + EvaluationStatusEvent warn + = EvaluationStatusEvent.newBuilder() + .setStatusLevel( StatusLevel.WARN ) + .setEventMessage( "Event detection was declared, but the half-life " + + "parameter was undefined. An attempt will be made to " + + "choose a reasonable default by inspecting the " + + "time-series data, but it is strongly recommended that " + + "you instead declare the 'half_life' explicitly as the " + + "default value may not be appropriate." ) + .build(); + events.add( warn ); + } + return Collections.unmodifiableList( events ); } diff --git a/wres-config/src/wres/config/yaml/components/EventDetection.java b/wres-config/src/wres/config/yaml/components/EventDetection.java index 4cf84b2fe6..08aafb851a 100644 --- a/wres-config/src/wres/config/yaml/components/EventDetection.java +++ b/wres-config/src/wres/config/yaml/components/EventDetection.java @@ -11,16 +11,23 @@ /** * Used for the detection of discrete events within time-series. - * @param datasets the datasets to use for event detection + * + * @param datasets the datasets to use for event detection, required + * @param method the event detection method, optional + * @param parameters the event detection parameters, optional unless required by a particular method */ @RecordBuilder @JsonDeserialize( using = EventDetectionDeserializer.class ) public record EventDetection( @JsonProperty( "dataset" ) Set datasets, - @JsonProperty( "method ") EventDetectionMethod method ) + @JsonProperty( "method " ) EventDetectionMethod method, + EventDetectionParameters parameters ) { /** * Creates an instance. - * @param datasets the datasets to use for event detection + * + * @param datasets the datasets to use for event detection, required + * @param method the event detection method, optional + * @param parameters the event detection parameters, optional unless required by a particular method */ public EventDetection { @@ -31,9 +38,16 @@ public record EventDetection( @JsonProperty( "dataset" ) Set geometryGroups, Map> +public class DurationIntervalDeserializer extends JsonDeserializer> { /** The name of the first duration within the interval. */ private final String firstName; @@ -25,15 +25,15 @@ public class DurationIntervalDeserializer extends JsonDeserializer deserialize( JsonParser jp, DeserializationContext context ) + public Pair deserialize( JsonParser jp, DeserializationContext context ) throws IOException { Objects.requireNonNull( jp ); ObjectReader mapper = ( ObjectReader ) jp.getCodec(); JsonNode node = mapper.readTree( jp ); - Duration first = DurationDeserializer.getDuration( mapper, node, this.firstName ); - Duration second = DurationDeserializer.getDuration( mapper, node, this.secondName ); + Duration first = DurationDeserializer.getDuration( mapper, node, this.firstName, "unit" ); + Duration second = DurationDeserializer.getDuration( mapper, node, this.secondName, "unit" ); return Pair.of( first, second ); } diff --git a/wres-config/src/wres/config/yaml/deserializers/EventDetectionDeserializer.java b/wres-config/src/wres/config/yaml/deserializers/EventDetectionDeserializer.java index bb7ddcbceb..54bf9d4526 100644 --- a/wres-config/src/wres/config/yaml/deserializers/EventDetectionDeserializer.java +++ b/wres-config/src/wres/config/yaml/deserializers/EventDetectionDeserializer.java @@ -1,6 +1,7 @@ package wres.config.yaml.deserializers; import java.io.IOException; +import java.time.Duration; import java.util.Collections; import java.util.Objects; import java.util.Set; @@ -16,6 +17,7 @@ import wres.config.yaml.components.EventDetectionBuilder; import wres.config.yaml.components.EventDetectionDataset; import wres.config.yaml.components.EventDetectionMethod; +import wres.config.yaml.components.EventDetectionParametersBuilder; /** * Custom deserializer for a {@link EventDetection}. @@ -24,6 +26,10 @@ */ public class EventDetectionDeserializer extends JsonDeserializer { + + /** Duration unit name. */ + private static final String DURATION_UNIT = "duration_unit"; + @Override public EventDetection deserialize( JsonParser jp, DeserializationContext context ) throws IOException @@ -70,10 +76,48 @@ else if ( datasetNode.isArray() ) method = reader.readValue( methodNode, EventDetectionMethod.class ); } + EventDetectionParametersBuilder parameters = EventDetectionParametersBuilder.builder(); + + if ( node.has( "parameters" ) ) + { + JsonNode parametersNode = node.get( "parameters" ); + if ( parametersNode.has( "window_size" ) ) + { + Duration windowSize = DurationDeserializer.getDuration( reader, + parametersNode, + "window_size", + DURATION_UNIT ); + parameters.windowSize( windowSize ); + } + if ( parametersNode.has( "start_radius" ) ) + { + Duration windowSize = DurationDeserializer.getDuration( reader, + parametersNode, + "start_radius", + DURATION_UNIT ); + parameters.windowSize( windowSize ); + } + if ( parametersNode.has( "half_life" ) ) + { + Duration windowSize = DurationDeserializer.getDuration( reader, + parametersNode, + "half_life", + DURATION_UNIT ); + parameters.windowSize( windowSize ); + } + if ( parametersNode.has( "minimum_event_duration" ) ) + { + Duration windowSize = DurationDeserializer.getDuration( reader, + parametersNode, + "minimum_event_duration", + DURATION_UNIT ); + parameters.windowSize( windowSize ); + } + } return EventDetectionBuilder.builder() .datasets( datasets ) .method( method ) + .parameters( parameters.build() ) .build(); } - } \ No newline at end of file diff --git a/wres-config/src/wres/config/yaml/deserializers/TimeScaleDeserializer.java b/wres-config/src/wres/config/yaml/deserializers/TimeScaleDeserializer.java index 7e59e2dfa5..f7480cbe4f 100644 --- a/wres-config/src/wres/config/yaml/deserializers/TimeScaleDeserializer.java +++ b/wres-config/src/wres/config/yaml/deserializers/TimeScaleDeserializer.java @@ -46,7 +46,7 @@ public wres.config.yaml.components.TimeScale deserialize( JsonParser jp, Deseria if ( node.has( "period" ) && node.has( "unit" ) ) { - java.time.Duration duration = DurationDeserializer.getDuration( mapper, node, "period" ); + java.time.Duration duration = DurationDeserializer.getDuration( mapper, node, "period", "unit" ); Duration protoDuration = Duration.newBuilder() .setSeconds( duration.getSeconds() ) .setNanos( duration.getNano() ) diff --git a/wres-config/test/wres/config/yaml/DeclarationFactoryTest.java b/wres-config/test/wres/config/yaml/DeclarationFactoryTest.java index 950b57c44a..b9fb4e83d9 100644 --- a/wres-config/test/wres/config/yaml/DeclarationFactoryTest.java +++ b/wres-config/test/wres/config/yaml/DeclarationFactoryTest.java @@ -49,6 +49,8 @@ import wres.config.yaml.components.EventDetectionBuilder; import wres.config.yaml.components.EventDetectionDataset; import wres.config.yaml.components.EventDetectionMethod; +import wres.config.yaml.components.EventDetectionParameters; +import wres.config.yaml.components.EventDetectionParametersBuilder; import wres.config.yaml.components.FeatureGroups; import wres.config.yaml.components.FeatureGroupsBuilder; import wres.config.yaml.components.FeatureService; @@ -2423,7 +2425,7 @@ void testDeserializeWithEventDetection() throws IOException } @Test - void testDeserializeWithEventDetectionUsingExplicitDatasetAndMethod() throws IOException + void testDeserializeWithEventDetectionUsingExplicitDatasetAndMethodWithParameters() throws IOException { String yaml = """ observed: some_file.csv @@ -2431,13 +2433,21 @@ void testDeserializeWithEventDetectionUsingExplicitDatasetAndMethod() throws IOE event_detection: dataset: observed method: regina-ogden + parameters: + window_size: 1 + duration_unit: hours """; EvaluationDeclaration actual = DeclarationFactory.from( yaml ); + EventDetectionParameters parameters = + EventDetectionParametersBuilder.builder() + .windowSize( java.time.Duration.ofHours( 1 ) ) + .build(); EventDetection eventDetection = EventDetectionBuilder.builder() .datasets( Set.of( EventDetectionDataset.OBSERVED ) ) .method( EventDetectionMethod.REGINA_OGDEN ) + .parameters( parameters ) .build(); EvaluationDeclaration expected = EvaluationDeclarationBuilder.builder() diff --git a/wres-config/test/wres/config/yaml/DeclarationValidatorTest.java b/wres-config/test/wres/config/yaml/DeclarationValidatorTest.java index bb0b33cafe..e92e0bd1c4 100644 --- a/wres-config/test/wres/config/yaml/DeclarationValidatorTest.java +++ b/wres-config/test/wres/config/yaml/DeclarationValidatorTest.java @@ -2689,6 +2689,12 @@ void testExplicitPoolsAndEventDetectionProducesErrorAndWarnings() () -> assertTrue( DeclarationValidatorTest.contains( events, "Event detection was declared alongside reference date " + "pools", + StatusLevel.WARN ) ), + () -> assertTrue( DeclarationValidatorTest.contains( events, + "it is strongly recommended that you instead declare the 'window_size'", + StatusLevel.WARN ) ), + () -> assertTrue( DeclarationValidatorTest.contains( events, + "it is strongly recommended that you instead declare the 'half_life'", StatusLevel.WARN ) ) ); } diff --git a/wres-datamodel/src/wres/datamodel/time/TimeSeriesSlicer.java b/wres-datamodel/src/wres/datamodel/time/TimeSeriesSlicer.java index be759a9b63..438ea98b79 100644 --- a/wres-datamodel/src/wres/datamodel/time/TimeSeriesSlicer.java +++ b/wres-datamodel/src/wres/datamodel/time/TimeSeriesSlicer.java @@ -1738,8 +1738,9 @@ public static TimeSeries applyOffsetToValidTimes( TimeSeries toTransfo } /** - * Snips the input series to the prescribed time window. Only snips lead durations with respect to reference times - * with the type {@link ReferenceTimeType#T0}. + * Snips the input series to the prescribed time window using a right-closed interval for each time dimension, + * i.e., the upper bound is included, the lower bound is excluded. Only snips lead durations with respect to + * reference times with the type {@link ReferenceTimeType#T0}. * * @param the time-series event value type * @param toSnip the time-series to snip @@ -1764,18 +1765,19 @@ public static TimeSeries snip( TimeSeries toSnip, TimeWindowOuter snip snipTo.getLatestValidTime() ); TimeWindowOuter partialSnip = TimeWindowOuter.of( inner ); - LOGGER.debug( "Snipping paired time-series {} to the time window of {}.", + LOGGER.debug( "Snipping time-series {} to the time window of {}.", toSnip.hashCode(), partialSnip ); returnMe = TimeSeriesSlicer.filter( returnMe, partialSnip ); // For all other reference time types, filter the datetimes only - if ( toSnip.getReferenceTimes().containsKey( ReferenceTimeType.T0 ) + if ( toSnip.getReferenceTimes() + .containsKey( ReferenceTimeType.T0 ) && !snipTo.bothLeadDurationsAreUnbounded() ) { - LOGGER.debug( "Additionally snipping paired time-series {} to lead durations ({},{}] for the reference " - + "time type of {}.", + LOGGER.debug( "Additionally snipping time-series {} to lead durations ({},{}] for the reference time " + + "type of {}.", toSnip.hashCode(), snipTo.getEarliestLeadDuration(), snipTo.getLatestLeadDuration(), diff --git a/wres-eventdetection/src/wres/eventdetection/EventDetectionException.java b/wres-eventdetection/src/wres/eventdetection/EventDetectionException.java index c88bfe4dfa..4c40d06e22 100644 --- a/wres-eventdetection/src/wres/eventdetection/EventDetectionException.java +++ b/wres-eventdetection/src/wres/eventdetection/EventDetectionException.java @@ -11,11 +11,10 @@ public class EventDetectionException extends RuntimeException * Constructs a {@link EventDetectionException} with the specified message. * * @param message the message. - * @param cause the cause of the exception */ - public EventDetectionException( final String message, final Throwable cause ) + public EventDetectionException( String message ) { - super( message, cause ); + super( message ); } } diff --git a/wres-eventdetection/src/wres/eventdetection/EventDetector.java b/wres-eventdetection/src/wres/eventdetection/EventDetector.java index 73825e8b62..7aaf99afe9 100644 --- a/wres-eventdetection/src/wres/eventdetection/EventDetector.java +++ b/wres-eventdetection/src/wres/eventdetection/EventDetector.java @@ -2,7 +2,6 @@ import java.util.Set; -import wres.config.yaml.components.EventDetection; import wres.datamodel.time.TimeSeries; import wres.datamodel.time.TimeWindowOuter; @@ -19,11 +18,10 @@ public interface EventDetector * Performs event detection. * * @param timeSeries the time-series data - * @param parameters the event detection parameters * @return the detected events - * @throws NullPointerException if either input is null + * @throws NullPointerException if the input is null * @throws EventDetectionException if the event detection fails for any other reason */ - Set detect( TimeSeries timeSeries, EventDetection parameters ); + Set detect( TimeSeries timeSeries ); } diff --git a/wres-eventdetection/src/wres/eventdetection/EventDetectorFactory.java b/wres-eventdetection/src/wres/eventdetection/EventDetectorFactory.java index 5ba4b42fcb..022573c8f5 100644 --- a/wres-eventdetection/src/wres/eventdetection/EventDetectorFactory.java +++ b/wres-eventdetection/src/wres/eventdetection/EventDetectorFactory.java @@ -1,7 +1,9 @@ package wres.eventdetection; import java.util.Objects; + import wres.config.yaml.components.EventDetectionMethod; +import wres.config.yaml.components.EventDetectionParameters; /** * A factory class that generates {@link EventDetector} for detecting events in time-series data. @@ -16,18 +18,21 @@ public class EventDetectorFactory * Returns an {@link EventDetector} for the named {@link EventDetectionMethod}. * * @param method the method, required + * @param parameters the event detection parameters * @return the detector * @throws NullPointerException if the method is null * @throws IllegalArgumentException if the method is unsupported */ - public static EventDetector getEventDetector( EventDetectionMethod method ) + public static EventDetector getEventDetector( EventDetectionMethod method, + EventDetectionParameters parameters ) { Objects.requireNonNull( method ); + Objects.requireNonNull( parameters ); if ( method == EventDetectionMethod.REGINA_OGDEN || method == EventDetectionMethod.DEFAULT ) { - return ReginaOgdenEventDetector.of(); + return ReginaOgdenEventDetector.of( parameters ); } throw new IllegalArgumentException( "Unsupported event detection method: " + method + "." ); diff --git a/wres-eventdetection/src/wres/eventdetection/ReginaOgdenEventDetector.java b/wres-eventdetection/src/wres/eventdetection/ReginaOgdenEventDetector.java index 60016b2fd6..14575850e8 100644 --- a/wres-eventdetection/src/wres/eventdetection/ReginaOgdenEventDetector.java +++ b/wres-eventdetection/src/wres/eventdetection/ReginaOgdenEventDetector.java @@ -1,20 +1,40 @@ package wres.eventdetection; +import java.time.Duration; +import java.time.Instant; +import java.util.Arrays; +import java.util.Collections; +import java.util.Comparator; +import java.util.Iterator; import java.util.Objects; import java.util.Set; +import java.util.SortedSet; +import java.util.TreeSet; +import java.util.stream.Collectors; +import com.google.protobuf.Timestamp; +import org.apache.commons.lang3.ArrayUtils; +import org.apache.commons.math3.stat.descriptive.rank.Median; +import org.apache.commons.math3.stat.ranking.NaNStrategy; import org.slf4j.Logger; import org.slf4j.LoggerFactory; -import wres.config.yaml.components.EventDetection; +import wres.config.yaml.components.EventDetectionParameters; +import wres.config.yaml.components.EventDetectionParametersBuilder; +import wres.datamodel.scale.TimeScaleOuter; +import wres.datamodel.time.Event; import wres.datamodel.time.TimeSeries; +import wres.datamodel.time.TimeSeriesSlicer; import wres.datamodel.time.TimeWindowOuter; +import wres.statistics.MessageFactory; +import wres.statistics.generated.TimeWindow; /** * An implementation of the method described in * Regina and Ogden (2021). Ported from the original Python code, available here: * HydroTools. * + * @author Jason Regina * @author James Brown */ public class ReginaOgdenEventDetector implements EventDetector @@ -22,33 +42,770 @@ public class ReginaOgdenEventDetector implements EventDetector /** Logger. */ private static final Logger LOGGER = LoggerFactory.getLogger( ReginaOgdenEventDetector.class ); + /** Median calculator. */ + private static final Median MEDIAN = new Median().withNaNStrategy( NaNStrategy.REMOVED ); + + /** The detection parameters. */ + private final EventDetectionParameters parameters; + @Override - public Set detect( TimeSeries timeSeries, EventDetection parameters ) + public Set detect( TimeSeries timeSeries ) { Objects.requireNonNull( timeSeries ); - Objects.requireNonNull( parameters ); LOGGER.debug( "Performing event detection using the Regina-Ogden method, which is described here: " + "https://doi.org/10.1002/hyp.14405" ); - // Unbounded time window, placeholder - return Set.of( TimeWindowOuter.of( wres.statistics.MessageFactory.getTimeWindow() ) ); + if ( timeSeries.getEvents() + .size() < 2 ) + { + LOGGER.debug( "Cannot detect events in a time-series that contains fewer than two event values: {}.", + timeSeries.getEvents() + .size() ); + + return Set.of(); + } + + return this.detectEvents( timeSeries, + this.getParameters() ); } /** * Creates an instance. + * + * @param parameters the event detection parameters * @return a detector instance + * @throws NullPointerException if the parameters are null + * @throws IllegalArgumentException if the parameters are invalid + */ + + static ReginaOgdenEventDetector of( EventDetectionParameters parameters ) + { + return new ReginaOgdenEventDetector( parameters ); + } + + /** + * @return the event detection parameters + */ + private EventDetectionParameters getParameters() + { + return this.parameters; + } + + /** + * Performs event detection with the supplied parameter. Equivalent to the list_events method in the + * original codebase. Models the trend in a time series by taking the maximum of two rolling minimum filters + * applied in a forward and backward fashion. Remove the trend and residual components. The method aims to produce + * a detrended time series with a median of 0.0. It assumes any residual components less than twice the detrended + * median are random noise. + * + * @param timeSeries the time-series + * @param parameters the event detection parameters + * @return the events + */ + + private Set detectEvents( TimeSeries timeSeries, + EventDetectionParameters parameters ) + { + parameters = this.setDefaultParameterValues( parameters, timeSeries ); + + // Smooth the time-series + TimeSeries smoothed = this.exponentialMovingAverage( timeSeries, parameters.halfLife() ); + + // Detrend the time-series + TimeSeries detrended = this.detrend( smoothed, parameters.halfLife(), parameters.windowSize() ); + + // Mark the events based on the smoothed, detrended flow + TimeSeries mapped = TimeSeriesSlicer.transform( detrended, + t -> t > 0.0, + m -> m ); + + // Filter events by duration + TimeSeries filtered = this.filterEventsByDuration( mapped, parameters.minimumEventDuration() ); + + // Get the provisional events + Set provisional = this.getEventBoundaries( filtered ); + + // Refine the events + return this.refineEvents( detrended, provisional, parameters.startRadius() ); + } + + /** + * Sets the default parameter values, inspecting the timeseries where useful. + * @param parameters parameters + * @param timeSeries the time-series */ - static ReginaOgdenEventDetector of() + private EventDetectionParameters setDefaultParameterValues( EventDetectionParameters parameters, + TimeSeries timeSeries ) { - return new ReginaOgdenEventDetector(); + EventDetectionParametersBuilder builder = EventDetectionParametersBuilder.builder( parameters ); + if ( Objects.isNull( builder.minimumEventDuration() ) + && Objects.nonNull( builder.halfLife() ) ) + { + LOGGER.warn( "When setting the event detection parameters, discovered that the half-life was defined and " + + "the minimum event duration was undefined. Setting the minimum event duration to the " + + "half-life of {}.", builder.halfLife() ); + builder.minimumEventDuration( builder.halfLife() ); + } + + if ( Objects.isNull( builder.startRadius() ) ) + { + builder.startRadius( Duration.ZERO ); + } + + if ( Objects.isNull( builder.minimumEventDuration() ) ) + { + builder.minimumEventDuration( Duration.ZERO ); + } + + this.setSeriesSpecificParameterDefaults( builder, timeSeries, parameters ); + + return builder.build(); + } + + /** + * Sets the time-series specific event detection parameter estimates. + * + * @param builder the builder + * @param timeSeries the time-series + * @param parameters the declared parameters + */ + private void setSeriesSpecificParameterDefaults( EventDetectionParametersBuilder builder, + TimeSeries timeSeries, + EventDetectionParameters parameters ) + { + // Calculate the series-specific parameter defaults + if ( Objects.isNull( builder.halfLife() ) + || Objects.isNull( builder.windowSize() ) ) + { + Duration averageTimestep = this.getAverageTimestep( timeSeries ); + + if ( Objects.isNull( builder.halfLife() ) ) + { + if ( Objects.nonNull( builder.windowSize() ) ) + { + builder.halfLife( builder.windowSize() + .dividedBy( 10 ) ); + LOGGER.warn( "When performing event detection with the Regina-Ogden method, the half life was " + + "undefined. However, the window size was defined. The default half life is {}, " + + "which is one tenth of the window size. This default may not be appropriate and it " + + "is strongly recommended that you set the half life explicitly using the " + + "'half_life' parameter.", + builder.halfLife() ); + } + else + { + builder.halfLife( averageTimestep.multipliedBy( 10 ) ); + LOGGER.warn( "When performing event detection with the Regina-Ogden method, the half life was " + + "undefined. The default half life is {}, which is ten times the average time-step " + + "associated with the time-series used for event detection. This default may not be " + + "appropriate and it is strongly recommended that you set the half life explicitly using " + + "the 'half_life' parameter.", + builder.halfLife() ); + } + } + + if ( Objects.isNull( builder.windowSize() ) ) + { + // Check original halflife, not current estimate + if ( Objects.nonNull( parameters.halfLife() ) ) + { + builder.windowSize( parameters.halfLife() + .multipliedBy( 10 ) ); + LOGGER.warn( "When performing event detection with the Regina-Ogden method, the window wize for " + + "smoothing and detecting trends was undefined. However, the half life was defined. " + + "The default window size is {}, which is ten times the half-life. This default may " + + "not be appropriate and it is strongly recommended that you set the window size " + + "explicitly using the 'window_size' parameter.", + builder.windowSize() ); + } + else + { + builder.windowSize( averageTimestep.multipliedBy( 100 ) ); + LOGGER.warn( "When performing event detection with the Regina-Ogden method, the window wize for " + + "smoothing and detecting trends was undefined. The default window size is {}, which is " + + "one hundred times the average time-step associated with the time-series used for event " + + "detection. This default may not be appropriate and it is strongly recommended that you " + + "set the window size explicitly using the 'window_size' parameter.", + builder.windowSize() ); + } + } + } + } + + /** + * Derives the event boundaries from the supplied dichotomous events. + * + * @param dichotomousEvents the event markers + * @return the event boundaries + * @throws NullPointerException if the input is null + */ + private Set getEventBoundaries( TimeSeries dichotomousEvents ) + { + Objects.requireNonNull( dichotomousEvents ); + + // Short circuit + if ( dichotomousEvents.getEvents() + .isEmpty() ) + { + return Set.of(); + } + + Set timeWindows = new TreeSet<>(); + + TimeWindow.Builder currentEvent = null; + + Event lastEvent = null; + for ( Event nextEvent : dichotomousEvents.getEvents() ) + { + if ( Boolean.TRUE.equals( nextEvent.getValue() ) ) + { + if ( Objects.isNull( currentEvent ) ) + { + Timestamp nextTime = MessageFactory.getTimestamp( nextEvent.getTime() ); + currentEvent = MessageFactory.getTimeWindow() + .toBuilder() + .setEarliestValidTime( nextTime ); + } + } + else + { + if ( Objects.nonNull( currentEvent ) ) + { + Timestamp nextTime = MessageFactory.getTimestamp( lastEvent.getTime() ); + currentEvent.setLatestValidTime( nextTime ); + TimeWindow timeWindow = currentEvent.build(); + TimeWindowOuter wrapped = TimeWindowOuter.of( timeWindow ); + timeWindows.add( wrapped ); + currentEvent = null; + } + } + + lastEvent = nextEvent; + } + + // Mop up an ongoing event + if ( Objects.nonNull( currentEvent ) ) + { + Timestamp nextTime = MessageFactory.getTimestamp( lastEvent.getTime() ); + currentEvent.setLatestValidTime( nextTime ); + TimeWindow timeWindow = currentEvent.build(); + TimeWindowOuter wrapped = TimeWindowOuter.of( timeWindow ); + timeWindows.add( wrapped ); + } + + return Collections.unmodifiableSet( timeWindows ); + } + + /** + * Removes the trend component from the prescribed time-series. + * + * @param timeSeries the time-series. + * @param halfLife the half-life + * @param smoothingWindow the smoothing window duration + * @return the detrended series + * @throws EventDetectionException if the smoothing window spans fewer than two time-steps, on average + */ + private TimeSeries detrend( TimeSeries timeSeries, Duration halfLife, Duration smoothingWindow ) + { + Objects.requireNonNull( timeSeries ); + Objects.requireNonNull( halfLife ); + + double[] eventValues = timeSeries.getEvents() + .stream() + .mapToDouble( Event::getValue ) + .toArray(); + + // Calculate the window size in integer steps + long windowSize = this.getWindowSizeFromDuration( timeSeries, smoothingWindow ); + + if ( windowSize < 2 ) + { + throw new EventDetectionException( "The window size for event detection must be at least two time-steps, " + + "as this is used to determine the trend in the time-series data. " + + "However, the smoothing duration of " + + parameters.windowSize() + + " spans only " + + windowSize + + " time-steps, on average, which is insufficient. Please increase the " + + "window size for smoothing and try again." ); + } + + double[] trend = this.rollingMinimum( eventValues, windowSize ); + + double[] detrended = new double[eventValues.length]; + + // Subtract the trend and remove the residuals + for ( int i = 0; i < detrended.length; i++ ) + { + detrended[i] = eventValues[i] - trend[i]; + } + + double residual = MEDIAN.evaluate( detrended ) * 2.0; + + SortedSet> detrendedEvents = new TreeSet<>(); + Iterator> iterator = timeSeries.getEvents() + .iterator(); + int current = 0; + while ( iterator.hasNext() ) + { + Event nextValue = iterator.next(); + + double value = Math.max( 0.0, detrended[current] - residual ); + + Event detrendedEvent = Event.of( nextValue.getTime(), value ); + detrendedEvents.add( detrendedEvent ); + current++; + } + + return TimeSeries.of( timeSeries.getMetadata(), detrendedEvents ); + } + + /** + * Smoothes a time-series with an exponential moving average. + * + * @param timeSeries the time-series + * @param halfLife the half life + * @return the smoothed series + */ + private TimeSeries exponentialMovingAverage( TimeSeries timeSeries, + Duration halfLife ) + { + Objects.requireNonNull( timeSeries ); + Objects.requireNonNull( halfLife ); + + if ( timeSeries.getEvents() + .size() < 2 ) + { + LOGGER.debug( "Cannot exponentially smooth a time-series with fewer than two values: {}.", timeSeries ); + return timeSeries; + } + + Iterator> iterator = timeSeries.getEvents() + .iterator(); + + Event firstEvent = iterator.next(); + double lastValue = firstEvent.getValue(); + + double weight = 1.0; + + SortedSet> smoothed = new TreeSet<>(); + smoothed.add( firstEvent ); // Initialize with first value + Instant lastTime = firstEvent.getTime(); + long halfLifeMillis = halfLife.toMillis(); + while ( iterator.hasNext() ) + { + Event currentValue = iterator.next(); + + // Calculate alpha based on the current timestep + Duration timestep = Duration.between( lastTime, currentValue.getTime() ); + long timestepMillis = timestep.toMillis(); + + double alpha = 1.0 - Math.exp( -Math.log( 2 ) / halfLifeMillis * timestepMillis ); + + // Handle NaN + if ( Double.isNaN( currentValue.getValue() ) ) + { + // First event, then insert NaN + if ( currentValue == firstEvent ) + { + Event smoothedEvent = Event.of( currentValue.getTime(), Double.NaN ); + smoothed.add( smoothedEvent ); + } + // Carry forward previous value + else + { + Event smoothedEvent = Event.of( currentValue.getTime(), lastValue ); + smoothed.add( smoothedEvent ); + } + + // Decay the weight + weight = weight * ( 1.0 - alpha ); + } + else + { + double smoothedValue = alpha * currentValue.getValue() + + ( 1.0 - alpha ) * lastValue * weight; + lastValue = smoothedValue; + Event smoothedEvent = Event.of( currentValue.getTime(), smoothedValue ); + smoothed.add( smoothedEvent ); + weight = weight * ( 1.0 - alpha ) + alpha; + } + + lastTime = currentValue.getTime(); + } + + return TimeSeries.of( timeSeries.getMetadata(), smoothed ); + } + + /** + * Calculates the average duration between time-steps and returns the whole number of times this average fits inside + * the smoothing window. + * + * @param timeSeries the time-series to inspect + * @param smoothingWindow the smoothing window duration + */ + + private long getWindowSizeFromDuration( TimeSeries timeSeries, Duration smoothingWindow ) + { + Duration averageTimestep = this.getAverageTimestep( timeSeries ); + + return smoothingWindow.dividedBy( averageTimestep ); + } + + /** + * Calculates the average time-step in the series. + * + * @param timeSeries the time-series to inspect + */ + + private Duration getAverageTimestep( TimeSeries timeSeries ) + { + // Sum the durations between times + Instant lastTime = null; + Duration sum = Duration.ZERO; + for ( Event next : timeSeries.getEvents() ) + { + Instant nextTime = next.getTime(); + if ( Objects.nonNull( lastTime ) ) + { + Duration between = Duration.between( lastTime, nextTime ); + sum = sum.plus( between ); + } + lastTime = nextTime; + } + + return sum.dividedBy( timeSeries.getEvents() + .size() ); + } + + /** + * Calculates a rolling minimum. + * + * @param values the time-series values in time order + * @param windowSize the window size + * @return the rolling minimum + */ + private double[] rollingMinimum( double[] values, long windowSize ) + { + // Forward pass + double[] forwardFiltered = this.rollingMinimumInner( values, windowSize ); + + // Reverse pass (flip the time series and values) + double[] reversedValues = new double[forwardFiltered.length]; + System.arraycopy( values, 0, reversedValues, 0, values.length ); + ArrayUtils.reverse( reversedValues ); + + // Backward pass on reversed data + double[] backwardFiltered = this.rollingMinimumInner( reversedValues, windowSize ); + + // Reverse the backward result to align with the original time series + ArrayUtils.reverse( backwardFiltered ); + + // Combine forward and backward passes by taking the maximum + double[] finalResult = new double[values.length]; + for ( int i = 0; i < values.length; i++ ) + { + finalResult[i] = Math.max( forwardFiltered[i], backwardFiltered[i] ); + } + + return finalResult; + } + + /** + * Calculates a rolling minimum value from the input. + * @param values the values + * @param windowSize the window size + * @return the rolling minimum + */ + private double[] rollingMinimumInner( double[] values, long windowSize ) + { + double[] result = new double[values.length]; + Arrays.fill( result, Double.NaN ); + + for ( int i = 0; i < values.length; i++ ) + { + long start = Math.max( 0, ( i + 1 ) - windowSize ); + double min = values[i]; + + for ( int j = i; j >= start; j-- ) + { + min = Math.min( min, values[j] ); + } + result[i] = min; + } + return result; + } + + /** + * Finds the time associated with the minimum value in the series using a search window. + * @param origin the origin + * @param radius the search radius + * @param timeSeries the time-series + * @return the time of the local minimum + */ + private Instant findLocalMinimum( Instant origin, + Duration radius, + TimeSeries timeSeries ) + { + Objects.requireNonNull( origin ); + Objects.requireNonNull( radius ); + Objects.requireNonNull( timeSeries ); + + // Convert radius to duration and subtract one instant from the left bound because the snip method is + // right-closed + Instant left = origin.minus( radius ) + .minus( Duration.ofNanos( 1 ) ); + Instant right = origin.plus( radius ); + + TimeWindow window = MessageFactory.getTimeWindow( left, right ); + TimeWindowOuter wrapped = TimeWindowOuter.of( window ); + + // Snip the series to the bounds + TimeSeries snipped = TimeSeriesSlicer.snip( timeSeries, wrapped ); + + if ( snipped.getEvents() + .isEmpty() ) + { + throw new IllegalArgumentException( "Could not find any time-series values within the search window: " + + wrapped + + "." ); + } + + // Sort the events by value + SortedSet> sorted = new TreeSet<>( Comparator.comparingDouble( Event::getValue ) ); + sorted.addAll( snipped.getEvents() ); + + return sorted.iterator() + .next() + .getTime(); + } + + /** + * Retains only those events that span a minimum duration. + * + * @param dichotomousEvents the time-series of event markers + * @param minimumEventDuration the minimum event duration + * @return the adjusted series with only those events included that span the minimum duration + * @throws NullPointerException if any input is null + */ + private TimeSeries filterEventsByDuration( TimeSeries dichotomousEvents, + Duration minimumEventDuration ) + { + Objects.requireNonNull( dichotomousEvents ); + Objects.requireNonNull( minimumEventDuration ); + + if ( minimumEventDuration == Duration.ZERO ) + { + LOGGER.debug( "Events were not filtered by duration because the minimum duration was {}.", Duration.ZERO ); + + return dichotomousEvents; + } + + // Create a series with only events + TimeSeries eventSeries = TimeSeriesSlicer.filterByEvent( dichotomousEvents, Event::getValue ); + + // Short-circuit + if ( eventSeries.getEvents() + .isEmpty() ) + { + LOGGER.debug( "Events were not filtered by duration because no events were detected." ); + + return dichotomousEvents; + } + + // Determine whether the events meet the minimum span + SortedSet> adjustedEvents = new TreeSet<>(); + Iterator> eventIterator = dichotomousEvents.getEvents() + .iterator(); + + // If the timescale period is non-instantaneous, then add to the event duration + Duration scaleAdjustment = this.getTimeScaleAdjustment( dichotomousEvents.getTimeScale() ); + + while ( eventIterator.hasNext() ) + { + Event start = eventIterator.next(); + + // Event is beginning, find the end and the span + if ( Boolean.TRUE.equals( start.getValue() ) ) + { + LOGGER.debug( "Event started at {}", start ); + + SortedSet> gathered = new TreeSet<>(); + gathered.add( start ); + Event last = start; + while ( eventIterator.hasNext() ) // NOSONAR only one conceptual break when event ends + { + Event next = eventIterator.next(); + + // Event is formally ending + if ( Boolean.FALSE.equals( next.getValue() ) ) + { + SortedSet> validEvents = this.getValidEvents( start.getTime(), + last.getTime(), + scaleAdjustment, + minimumEventDuration, + gathered ); + // Add the new state for all values within the event period + adjustedEvents.addAll( validEvents ); + + // Add the current non-event + adjustedEvents.add( next ); + + break; + } + // Ending because no more values, so next is last + else if ( !eventIterator.hasNext() ) + { + SortedSet> validEvents = this.getValidEvents( start.getTime(), + next.getTime(), + scaleAdjustment, + minimumEventDuration, + gathered ); + adjustedEvents.addAll( validEvents ); + + // Add the current event + adjustedEvents.add( next ); + + break; + } + else + { + gathered.add( next ); + } + + last = next; + } + } + else + { + adjustedEvents.add( start ); + } + } + + return TimeSeries.of( dichotomousEvents.getMetadata(), adjustedEvents ); + } + + /** + * Marks events that are within the minimum duration as events, otherwise non-events. + * @param startTime the event start time + * @param endTime the event end time + * @param scaleAdjustment the timescale adjustment + * @param minimumEventDuration the minimum event duration + * @param gathered the gathered data to adjust + * @return the valid events + */ + private SortedSet> getValidEvents( Instant startTime, + Instant endTime, + Duration scaleAdjustment, + Duration minimumEventDuration, + SortedSet> gathered ) + { + // Are the gathered values within an event of appropriate span? + Duration duration = Duration.between( startTime, endTime ) + .plus( scaleAdjustment ); + + boolean valid = duration.compareTo( minimumEventDuration ) >= 0; + + return gathered.stream() + .map( e -> Event.of( e.getTime(), valid ) ) + .collect( Collectors.toCollection( TreeSet::new ) ); + } + + /** + * @param timeScale the timescale + * @return {@link TimeScaleOuter#getPeriod()} if the input is non-instantaneous, else {@link Duration#ZERO} + */ + + private Duration getTimeScaleAdjustment( TimeScaleOuter timeScale ) + { + // If the timescale period is non-instantaneous, then add to the event duration + Duration scalePeriod = Duration.ZERO; + + if ( Objects.nonNull( timeScale ) + && !timeScale.isInstantaneous() ) + { + scalePeriod = timeScale.getPeriod(); + } + + return scalePeriod; + } + + /** + * Refines the event points, returning a revised series of events. + * + * @param timeSeries the time-series + * @param events the provisional events to refine + * @param startRadius the start radius for searching + * @return the refined event markers + * @throws NullPointerException if any input is null + */ + private Set refineEvents( TimeSeries timeSeries, + Set events, + Duration startRadius ) + { + Objects.requireNonNull( timeSeries ); + Objects.requireNonNull( events ); + Objects.requireNonNull( startRadius ); + + // Short-circuit + if ( events.isEmpty() ) + { + return events; + } + + Set adjustedWindows = new TreeSet<>(); + for ( TimeWindowOuter nextEvent : events ) + { + Instant eventStart = nextEvent.getEarliestValidTime(); + Instant refinedStart = this.findLocalMinimum( eventStart, + startRadius, + timeSeries ); + + TimeWindow adjusted = nextEvent.getTimeWindow() + .toBuilder() + .setEarliestValidTime( MessageFactory.getTimestamp( refinedStart ) ) + .build(); + + adjustedWindows.add( TimeWindowOuter.of( adjusted ) ); + } + + return Collections.unmodifiableSet( adjustedWindows ); } /** * Hidden constructor. + * @param parameters the event detection parameters */ - private ReginaOgdenEventDetector() + private ReginaOgdenEventDetector( EventDetectionParameters parameters ) { + Objects.requireNonNull( parameters ); + + if ( !Objects.equals( parameters.halfLife(), parameters.minimumEventDuration() ) ) + { + LOGGER.warn( "When performing event detection with the Regina-Ogden method, discovered a half-life " + + "of {} and a minimum event duration of {}. In general, the minimum event duration should " + + "match the half-life to reduce the number of false positives caused by measurement " + + "variability.", parameters.halfLife(), parameters.minimumEventDuration() ); + } + + if ( Objects.nonNull( parameters.minimumEventDuration() ) + && parameters.minimumEventDuration() + .isNegative() ) + { + throw new IllegalArgumentException( "The minimum event duration for event detection cannot be negative: " + + parameters.minimumEventDuration() ); + } + + if ( Objects.nonNull( parameters.startRadius() ) + && parameters.startRadius() + .isNegative() ) + { + throw new IllegalArgumentException( "The start radius for event detection cannot be a negative duration: " + + parameters.startRadius() ); + } + + this.parameters = parameters; } } diff --git a/wres-eventdetection/test/wres/eventdetection/ReginaOgdenEventDetectorTest.java b/wres-eventdetection/test/wres/eventdetection/ReginaOgdenEventDetectorTest.java new file mode 100644 index 0000000000..6c89693fc8 --- /dev/null +++ b/wres-eventdetection/test/wres/eventdetection/ReginaOgdenEventDetectorTest.java @@ -0,0 +1,225 @@ +package wres.eventdetection; + +import java.time.Duration; +import java.time.Instant; +import java.util.Map; +import java.util.Set; + +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +import wres.config.yaml.components.EventDetectionParameters; +import wres.config.yaml.components.EventDetectionParametersBuilder; +import wres.datamodel.scale.TimeScaleOuter; +import wres.datamodel.space.Feature; +import wres.datamodel.time.Event; +import wres.datamodel.time.TimeSeries; +import wres.datamodel.time.TimeSeriesMetadata; +import wres.datamodel.time.TimeWindowOuter; +import wres.statistics.MessageFactory; +import wres.statistics.generated.TimeWindow; + +/** + * Tests the {@link ReginaOgdenEventDetector}. + * + * @author James Brown + */ + +class ReginaOgdenEventDetectorTest +{ + @Test + void testDetectWithZeroTrendAndTwoEvents() + { + EventDetectionParameters parameters = EventDetectionParametersBuilder.builder() + .windowSize( Duration.ofHours( 6 ) ) + .minimumEventDuration( Duration.ZERO ) + .halfLife( Duration.ofHours( 2 ) ) + .build(); + + ReginaOgdenEventDetector detector = ReginaOgdenEventDetector.of( parameters ); + + Event one = Event.of( Instant.parse( "2079-12-03T00:00:00Z" ), 5.0 ); + Event two = Event.of( Instant.parse( "2079-12-03T01:00:00Z" ), 5.0 ); + Event three = Event.of( Instant.parse( "2079-12-03T02:00:00Z" ), 24.0 ); + Event four = Event.of( Instant.parse( "2079-12-03T03:00:00Z" ), 25.0 ); + Event five = Event.of( Instant.parse( "2079-12-03T04:00:00Z" ), 5.0 ); + Event six = Event.of( Instant.parse( "2079-12-03T05:00:00Z" ), 5.0 ); + Event seven = Event.of( Instant.parse( "2079-12-03T06:00:00Z" ), 5.0 ); + Event eight = Event.of( Instant.parse( "2079-12-03T07:00:00Z" ), 84.0 ); + Event nine = Event.of( Instant.parse( "2079-12-03T08:00:00Z" ), 85.0 ); + Event ten = Event.of( Instant.parse( "2079-12-03T09:00:00Z" ), 87.0 ); + Event eleven = Event.of( Instant.parse( "2079-12-03T10:00:00Z" ), 5.0 ); + Event twelve = Event.of( Instant.parse( "2079-12-03T11:00:00Z" ), 5.0 ); + + TimeSeriesMetadata metadata = TimeSeriesMetadata.of( Map.of(), + TimeScaleOuter.of(), + "foo", + Feature.of( MessageFactory.getGeometry( "bar" ) ), + "baz" ); + TimeSeries timeSeries = + new TimeSeries.Builder() + .addEvent( one ) + .addEvent( two ) + .addEvent( three ) + .addEvent( four ) + .addEvent( five ) + .addEvent( six ) + .addEvent( seven ) + .addEvent( eight ) + .addEvent( nine ) + .addEvent( ten ) + .addEvent( eleven ) + .addEvent( twelve ) + .setMetadata( metadata ) + .build(); + + Set actual = detector.detect( timeSeries ); + + Instant startOne = Instant.parse( "2079-12-03T03:00:00Z" ); + Instant endOne = Instant.parse( "2079-12-03T03:00:00Z" ); + + TimeWindow expectedInnerOne = MessageFactory.getTimeWindow() + .toBuilder() + .setEarliestValidTime( MessageFactory.getTimestamp( startOne ) ) + .setLatestValidTime( MessageFactory.getTimestamp( endOne ) ) + .build(); + + Instant startTwo = Instant.parse( "2079-12-03T08:00:00Z" ); + Instant endTwo = Instant.parse( "2079-12-03T10:00:00Z" ); + + TimeWindow expectedInnerTwo = MessageFactory.getTimeWindow() + .toBuilder() + .setEarliestValidTime( MessageFactory.getTimestamp( startTwo ) ) + .setLatestValidTime( MessageFactory.getTimestamp( endTwo ) ) + .build(); + + Set expected = Set.of( TimeWindowOuter.of( expectedInnerOne ), + TimeWindowOuter.of( expectedInnerTwo ) ); + + + assertEquals( expected, actual ); + } + + @Test + void testDetectWithTrendAndOneEvent() + { + TimeSeries testSeries = this.createTestSeries( false ); + + EventDetectionParameters parameters = EventDetectionParametersBuilder.builder() + .windowSize( Duration.ofDays( 7 ) ) + .minimumEventDuration( Duration.ZERO ) + .halfLife( Duration.ofHours( 6 ) ) + .build(); + + ReginaOgdenEventDetector detector = ReginaOgdenEventDetector.of( parameters ); + + Set actual = detector.detect( testSeries ); + + Instant start = Instant.parse( "2018-01-22T02:00:00Z" ); + Instant end = Instant.parse( "2018-01-29T00:00:00Z" ); + + TimeWindow expectedInner = MessageFactory.getTimeWindow() + .toBuilder() + .setEarliestValidTime( MessageFactory.getTimestamp( start ) ) + .setLatestValidTime( MessageFactory.getTimestamp( end ) ) + .build(); + + Set expected = Set.of( TimeWindowOuter.of( expectedInner ) ); + + // Assert that one event is detected + assertEquals( expected, actual ); + } + + @Test + void testDetectWithTrendAndOneEventAndNoise() + { + TimeSeries testSeries = this.createTestSeries( true ); + + EventDetectionParameters parameters = EventDetectionParametersBuilder.builder() + .halfLife( Duration.ofHours( 6 ) ) + .windowSize( Duration.ofDays( 7 ) ) + .minimumEventDuration( Duration.ofHours( 6 ) ) + .startRadius( Duration.ofHours( 7 ) ) + .build(); + + ReginaOgdenEventDetector detector = ReginaOgdenEventDetector.of( parameters ); + + Set actual = detector.detect( testSeries ); + + Instant start = Instant.parse( "2018-01-21T19:00:00Z" ); + Instant end = Instant.parse( "2018-01-29T00:00:00Z" ); + + TimeWindow expectedInner = MessageFactory.getTimeWindow() + .toBuilder() + .setEarliestValidTime( MessageFactory.getTimestamp( start ) ) + .setLatestValidTime( MessageFactory.getTimestamp( end ) ) + .build(); + + Set expected = Set.of( TimeWindowOuter.of( expectedInner ) ); + + // Assert that one event is detected + assertEquals( expected, actual ); + } + + /** + * @param addNoise is true to add noise, false otherwise + * @return the test series + */ + private TimeSeries createTestSeries( boolean addNoise ) + { + // Create a trend array + double[] trend = new double[1000]; + for ( int i = 0; i < 1000; i++ ) + { + trend[i] = 100.0 * Math.exp( -0.002 * ( i + 1 ) * 1.0 ); + } + + // Create the event flows, which span 500 times + double[] eventFlows = new double[500]; + for ( int i = 0; i < 500; i++ ) + { + eventFlows[i] = 1000.0 * Math.exp( -0.004 * ( i + 1 ) * 1.0 ); + } + + // Concatenate the event to zero flows, which span the first 500 times + double[] flow = new double[1000]; + System.arraycopy( eventFlows, 0, flow, 500, 500 ); // Second part: event flows + + // Add a noise component if needed + double[] noise = new double[1000]; + if ( addNoise ) + { + for ( int i = 0; i < 1000; i++ ) + { + noise[i] = i * 334 * Math.PI / 1000; // Simulating np.linspace(0.0, 334*np.pi, 1000) + noise[i] = 5.0 * Math.sin( noise[i] ); // noise + } + } + + // Sum the trend, flow and noise arrays to get total flow + double[] totalFlow = new double[1000]; + for ( int i = 0; i < 1000; i++ ) + { + totalFlow[i] = trend[i] + flow[i] + noise[i]; + } + + // Create a time-series + TimeSeriesMetadata metadata = TimeSeriesMetadata.of( Map.of(), + TimeScaleOuter.of(), + "foo", + Feature.of( MessageFactory.getGeometry( "bar" ) ), + "baz" ); + TimeSeries.Builder timeSeries = new TimeSeries.Builder().setMetadata( metadata ); + + Instant startDate = Instant.parse( "2018-01-01T00:00:00Z" ); + for ( int i = 0; i < totalFlow.length; i++ ) + { + Instant nextDate = startDate.plus( Duration.ofHours( i ) ); + Event nextEvent = Event.of( nextDate, totalFlow[i] ); + timeSeries.addEvent( nextEvent ); + } + + return timeSeries.build(); + } +} \ No newline at end of file