autoPickWaveforms(final Collection frequencyBandParameters) {
return velocityMeasurements.parallelStream().filter(vel -> vel.getWaveform() != null).filter(vel -> vel.getWaveform().getAssociatedPicks() != null).map(vel -> {
SharedFrequencyBandParameters params = frequencyBandParameters.get(new FrequencyBand(vel.getWaveform().getLowFrequency(), vel.getWaveform().getHighFrequency()));
- Optional pick = vel.getWaveform().getAssociatedPicks().stream().filter(p -> PICK_TYPES.AP.name().equalsIgnoreCase(p.getPickType())).findFirst();
+ Optional pick = vel.getWaveform()
+ .getAssociatedPicks()
+ .stream()
+ .filter(p -> p.getPickType() != null && PICK_TYPES.AP.name().equalsIgnoreCase(p.getPickType().trim()))
+ .findFirst();
+
+ Optional endPick = vel.getWaveform()
+ .getAssociatedPicks()
+ .stream()
+ .filter(p -> p.getPickType() != null && PICK_TYPES.F.name().equalsIgnoreCase(p.getPickType().trim()))
+ .findFirst();
- Optional endPick = vel.getWaveform().getAssociatedPicks().stream().filter(p -> PICK_TYPES.F.name().equalsIgnoreCase(p.getPickType())).findFirst();
if ((!endPick.isPresent() || pick.isPresent()) && params != null) {
vel.getWaveform().getAssociatedPicks().forEach(p -> p.setWaveform(null));
vel.getWaveform().getAssociatedPicks().clear();
@@ -158,31 +164,43 @@ private Collection autoPickWaveforms(final Collection maxlength) {
+ stopTime = maxlength;
+ }
+
+ WaveformPick autoPick = new WaveformPick().setPickType(PICK_TYPES.F.name())
+ .setPickName(PICK_TYPES.F.getPhase())
+ .setWaveform(vel.getWaveform())
+ .setPickTimeSecFromOrigin((float) stopTime);
WaveformPick startPick = new WaveformPick().setPickType(PICK_TYPES.AP.name())
- .setPickName(PICK_TYPES.AP.name())
+ .setPickName(PICK_TYPES.AP.getPhase())
.setWaveform(vel.getWaveform())
.setPickTimeSecFromOrigin((float) startTime.subtractD(originTime));
- if (autoPick.getPickTimeSecFromOrigin() > startPick.getPickTimeSecFromOrigin()) {
- vel.getWaveform().getAssociatedPicks().add(autoPick);
- vel.getWaveform().getAssociatedPicks().add(startPick);
- }
+ vel.getWaveform().getAssociatedPicks().add(autoPick);
+ vel.getWaveform().getAssociatedPicks().add(startPick);
+
}
+
return vel;
}).filter(Objects::nonNull).collect(Collectors.toList());
}
diff --git a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/SiteCalibrationServiceImpl.java b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/SiteCalibrationServiceImpl.java
index 99007aa8..959901b0 100644
--- a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/SiteCalibrationServiceImpl.java
+++ b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/SiteCalibrationServiceImpl.java
@@ -18,6 +18,9 @@
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
+import java.util.Set;
+import java.util.TreeMap;
+import java.util.function.Function;
import java.util.stream.Collectors;
import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
@@ -73,8 +76,10 @@ public Map> measureSite
//Input
Map>> freqBandEvidStaMeasurementsMap = mapToEventAndStation(dataByFreqBand);
+ Map, Map>> weightFunctionMapByEvent = new HashMap<>();
+
//Step 1
- Map> staFreqBandSiteCorrectionMapMdac = new HashMap<>();
+ Map> staFreqBandSiteCorrectionMapReferenceEvents = new HashMap<>();
//Step 2
Map> averageMapByEvent = new HashMap<>();
//Step 3
@@ -82,29 +87,33 @@ public Map> measureSite
//Result
Map> siteCorrections = new HashMap<>();
+ //0) Determine if we have a RefMW in the dataset with a GT stress drop
+ //0-1A) If yes then omit everything above the corner frequency for MDAC model only events
+ //0-1B) Add it to the GT spectra event map
+ //0-1C) Weight the Mw fit for that event 1.0 at all bands
+ //0-2A) If no then weight the corner frequency highly, the low middle, and the high low
+
//1) Generate spectra for reference events and get the site correction for each station that saw it.
for (Entry>> evidStaMap : freqBandEvidStaMeasurementsMap.entrySet()) {
FrequencyBand freqBand = evidStaMap.getKey();
Double lowFreq = freqBand.getLowFrequency();
Double highFreq = freqBand.getHighFrequency();
double centerFreq = (highFreq + lowFreq) / 2;
-
for (Entry> staMwMap : evidStaMap.getValue().entrySet()) {
Event evid = staMwMap.getKey();
if (refMws.containsKey(evid.getEventId())) {
-
List refMwsParams = refMws.get(evid.getEventId());
ReferenceMwParameters refMw = refMwsParams.stream().findFirst().orElse(null);
if (refMw != null) {
double mw = refMw.getRefMw();
- MdacParametersFI mdacFiEntry = mdacFI;
+ MdacParametersFI mdacFiEntry = new MdacParametersFI(mdacFI);
if (refMw.getStressDropInMpa() != null && refMw.getStressDropInMpa() != 0.0) {
// If we know a MPA stress drop for this reference
// event we want to use stress instead of apparent
- // stress
- // so we set Psi == 0.0 to use Sigma as stress drop
+ // stress so we set Psi == 0.0 to use Sigma as stress drop
mdacFiEntry.setSigma(refMw.getStressDropInMpa());
mdacFiEntry.setPsi(0.0);
+ weightFunctionMapByEvent.put(evid, this::noWeights);
}
double[] refSpectra = mdac.calculateMdacSourceSpectra(psRows, mdacFiEntry, centerFreq, mw, UNUSED_DISTANCE);
@@ -112,18 +121,18 @@ public Map> measureSite
double amp = staMwEntry.getValue().getPathCorrected();
- if (!staFreqBandSiteCorrectionMapMdac.containsKey(staMwEntry.getKey())) {
- staFreqBandSiteCorrectionMapMdac.put(staMwEntry.getKey(), new HashMap());
+ if (!staFreqBandSiteCorrectionMapReferenceEvents.containsKey(staMwEntry.getKey())) {
+ staFreqBandSiteCorrectionMapReferenceEvents.put(staMwEntry.getKey(), new HashMap());
}
- if (!staFreqBandSiteCorrectionMapMdac.get(staMwEntry.getKey()).containsKey(freqBand)) {
- staFreqBandSiteCorrectionMapMdac.get(staMwEntry.getKey()).put(freqBand, new SummaryStatistics());
+ if (!staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).containsKey(freqBand)) {
+ staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).put(freqBand, new SummaryStatistics());
}
// Output should be Dyne-cm
double refAmp = Math.log10(refSpectra[1]) + DYNE_LOG10_ADJUSTMENT;
double ampDiff = refAmp - amp;
- staFreqBandSiteCorrectionMapMdac.get(staMwEntry.getKey()).get(freqBand).addValue(ampDiff);
+ staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).get(freqBand).addValue(ampDiff);
}
}
}
@@ -133,12 +142,10 @@ public Map> measureSite
//2) For every station with a site correction measured apply it to every other event and get average site term for every frequency band
for (Entry>> evidStaMap : freqBandEvidStaMeasurementsMap.entrySet()) {
FrequencyBand freqBand = evidStaMap.getKey();
-
- for (Entry> staMwMap : evidStaMap.getValue().entrySet()) {
- Event evid = staMwMap.getKey();
- for (Entry staMwEntry : staMwMap.getValue().entrySet()) {
- if (staFreqBandSiteCorrectionMapMdac.containsKey(staMwEntry.getKey()) && staFreqBandSiteCorrectionMapMdac.get(staMwEntry.getKey()).containsKey(freqBand)) {
-
+ for (Entry> evidStaEntry : evidStaMap.getValue().entrySet()) {
+ Event evid = evidStaEntry.getKey();
+ for (Entry staMwEntry : evidStaEntry.getValue().entrySet()) {
+ if (staFreqBandSiteCorrectionMapReferenceEvents.containsKey(staMwEntry.getKey()) && staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).containsKey(freqBand)) {
double amp = staMwEntry.getValue().getPathCorrected();
if (!averageMapByEvent.containsKey(evid)) {
averageMapByEvent.put(evid, new HashMap());
@@ -147,22 +154,23 @@ public Map> measureSite
averageMapByEvent.get(evid).put(freqBand, new SummaryStatistics());
}
- double refAmp = amp + staFreqBandSiteCorrectionMapMdac.get(staMwEntry.getKey()).get(freqBand).getMean();
-
+ double refAmp = amp + staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).get(freqBand).getMean();
averageMapByEvent.get(evid).get(freqBand).addValue(refAmp);
}
}
}
}
- //3) Offset the original measurements by the average site term for each station/frequency band to get the final site terms
+ //3) For all measurements offset by the average site term for each station/frequency band to get the final site terms
for (Entry>> evidStaMap : freqBandEvidStaMeasurementsMap.entrySet()) {
FrequencyBand freqBand = evidStaMap.getKey();
-
- for (Entry> staMwMap : evidStaMap.getValue().entrySet()) {
- Event evid = staMwMap.getKey();
+ for (Entry> evidStaEntry : evidStaMap.getValue().entrySet()) {
+ Event evid = evidStaEntry.getKey();
if (averageMapByEvent.containsKey(evid) && averageMapByEvent.get(evid).containsKey(freqBand)) {
- for (Entry staMwEntry : staMwMap.getValue().entrySet()) {
+ for (Entry staMwEntry : evidStaEntry.getValue().entrySet()) {
+ if (!weightFunctionMapByEvent.containsKey(evid)) {
+ weightFunctionMapByEvent.put(evid, this::lowerFreqHigherWeights);
+ }
double amp = staMwEntry.getValue().getPathCorrected();
if (!staFreqBandSiteCorrectionMapAverage.containsKey(staMwEntry.getKey())) {
staFreqBandSiteCorrectionMapAverage.put(staMwEntry.getKey(), new HashMap());
@@ -170,11 +178,15 @@ public Map> measureSite
if (!staFreqBandSiteCorrectionMapAverage.get(staMwEntry.getKey()).containsKey(freqBand)) {
staFreqBandSiteCorrectionMapAverage.get(staMwEntry.getKey()).put(freqBand, new SummaryStatistics());
}
-
- double refAmp = averageMapByEvent.get(evid).get(freqBand).getMean();
- double ampDiff = refAmp - amp;
-
- staFreqBandSiteCorrectionMapAverage.get(staMwEntry.getKey()).get(freqBand).addValue(ampDiff);
+ //We want to keep the reference events as-is
+ if (staFreqBandSiteCorrectionMapReferenceEvents.containsKey(staMwEntry.getKey())
+ && staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).containsKey(freqBand)) {
+ staFreqBandSiteCorrectionMapAverage.get(staMwEntry.getKey()).put(freqBand, staFreqBandSiteCorrectionMapReferenceEvents.get(staMwEntry.getKey()).get(freqBand));
+ } else {
+ double refAmp = averageMapByEvent.get(evid).get(freqBand).getMean();
+ double ampDiff = refAmp - amp;
+ staFreqBandSiteCorrectionMapAverage.get(staMwEntry.getKey()).get(freqBand).addValue(ampDiff);
+ }
}
}
}
@@ -184,12 +196,10 @@ public Map> measureSite
averageMapByEvent.clear();
for (Entry>> evidStaMap : freqBandEvidStaMeasurementsMap.entrySet()) {
FrequencyBand freqBand = evidStaMap.getKey();
-
- for (Entry> staMwMap : evidStaMap.getValue().entrySet()) {
- Event evid = staMwMap.getKey();
- for (Entry staMwEntry : staMwMap.getValue().entrySet()) {
+ for (Entry> evidStaEntry : evidStaMap.getValue().entrySet()) {
+ Event evid = evidStaEntry.getKey();
+ for (Entry staMwEntry : evidStaEntry.getValue().entrySet()) {
if (staFreqBandSiteCorrectionMapAverage.containsKey(staMwEntry.getKey()) && staFreqBandSiteCorrectionMapAverage.get(staMwEntry.getKey()).containsKey(freqBand)) {
-
double amp = staMwEntry.getValue().getPathCorrected();
if (!averageMapByEvent.containsKey(evid)) {
averageMapByEvent.put(evid, new HashMap());
@@ -204,14 +214,17 @@ public Map> measureSite
}
}
+ //TODO: Weight bands by number of points and std-dev instead
+
// 5) Measure the MW values per event
- List measuredMws = spectraCalc.measureMws(averageMapByEvent, PICK_TYPES.LG);
+ List measuredMws = spectraCalc.measureMws(averageMapByEvent, weightFunctionMapByEvent, PICK_TYPES.LG);
//Return results as a set of Site corrections
for (Entry> freqBandCorrectionMap : staFreqBandSiteCorrectionMapAverage.entrySet()) {
Station station = freqBandCorrectionMap.getKey();
+ Set> bandEntries = freqBandCorrectionMap.getValue().entrySet();
- for (Entry freqBandEntry : freqBandCorrectionMap.getValue().entrySet()) {
+ for (Entry freqBandEntry : bandEntries) {
SiteFrequencyBandParameters siteParam = new SiteFrequencyBandParameters();
siteParam.setStation(station);
siteParam.setHighFrequency(freqBandEntry.getKey().getHighFrequency());
@@ -232,6 +245,24 @@ public Map> measureSite
return siteCorrections;
}
+ private Map noWeights(Map frequencies) {
+ Map weightMap = new TreeMap<>();
+ for (Double frequency : frequencies.keySet()) {
+ weightMap.put(frequency, 1d);
+ }
+ return weightMap;
+ }
+
+ private Map lowerFreqHigherWeights(Map frequencies) {
+ Map weightMap = new TreeMap<>(frequencies);
+ int count = 0;
+ for (Entry entry : weightMap.entrySet()) {
+ entry.setValue(count == 0 || count == 2 ? 0.5 : count == 1 ? 1.0 : count == 3 ? 0.25 : 0.1);
+ count++;
+ }
+ return weightMap;
+ }
+
public static Map>> mapToEventAndStation(Map> dataByFreqBand) {
Map>> data = new HashMap<>();
diff --git a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CalibrationCurveFitter.java b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CalibrationCurveFitter.java
index e9f37e3b..55c0e2c8 100755
--- a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CalibrationCurveFitter.java
+++ b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CalibrationCurveFitter.java
@@ -91,7 +91,7 @@ public class CalibrationCurveFitter {
protected static final double YVV_MAX = 6.01;
protected static final double V_DIST_MAX = 1600;
- protected static final double YBB_MIN = -3.0E-2;
+ protected static final double YBB_MIN = -12.0E-2;
protected static final double YBB_MAX = 0.0005;
protected static final double B_DIST_MAX = 1550;
diff --git a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CodaSNREndTimePicker.java b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CodaSNREndTimePicker.java
index 4c6724a4..d6925592 100644
--- a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CodaSNREndTimePicker.java
+++ b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/CodaSNREndTimePicker.java
@@ -25,9 +25,9 @@ public class CodaSNREndTimePicker implements EndTimePicker {
private static final double BAD_PICK = -100.0;
@Override
- public double getEndTime(float[] waveform, double sampleRate, double startTimeEpochSeconds, int startOffset, double minLengthSec, double maxLengthSec, double minimumSnr) {
+ public double getEndTime(float[] waveform, double sampleRate, double startTimeEpochSeconds, int startOffset, double minLengthSec, double maxLengthSec, double minimumSnr, double noise) {
- Double snrPick = getSnrEndPick(waveform, sampleRate, startOffset, minLengthSec, maxLengthSec, minimumSnr, 20);
+ Double snrPick = getSnrEndPick(waveform, sampleRate, startOffset, minLengthSec, maxLengthSec, minimumSnr, noise, 20);
Double overallPick;
if (!Double.isNaN(snrPick)) {
@@ -39,14 +39,10 @@ public double getEndTime(float[] waveform, double sampleRate, double startTimeEp
return startTimeEpochSeconds + overallPick;
}
- private Double getSnrEndPick(final float[] waveform, final double sampleRate, int startOffset, final double minLengthSec, final double maxLengthSec, final double minimumSnr,
+ private Double getSnrEndPick(final float[] waveform, final double sampleRate, int startOffset, final double minLengthSec, final double maxLengthSec, final double minimumSnr, final double noise,
final int windowSize) {
- DescriptiveStatistics stats = new DescriptiveStatistics();
- for (float value : waveform) {
- stats.addValue(value);
- }
- Double noise = stats.getPercentile(50);
- DescriptiveStatistics obs = new DescriptiveStatistics((int) (windowSize * sampleRate));
+ int obsWindow = (int) (windowSize * sampleRate);
+ DescriptiveStatistics obs = new DescriptiveStatistics(obsWindow);
DescriptiveStatistics spike = new DescriptiveStatistics((int) (5 * sampleRate));
double snrTimePick = BAD_PICK;
@@ -55,15 +51,26 @@ private Double getSnrEndPick(final float[] waveform, final double sampleRate, in
if (waveform.length > startOffset && waveform.length - startOffset > minSamples) {
int stopIdx = waveform.length > maxSamples ? maxSamples : waveform.length;
- for (int i = startOffset; i < stopIdx; i++) {
- obs.addValue(waveform[i]);
- spike.addValue(waveform[i]);
- if (obs.getN() >= windowSize) {
- //TODO: .5 is a WAG on a data set, should try to find a better formulation that takes into account things like sample rate, window, and overall SNR
- if (obs.getMean() - minimumSnr <= noise || spike.getMean() > (obs.getMean() + .5d)) {
- break;
- } else {
- snrTimePick = (i - spike.getN()) / sampleRate;
+ if (waveform[startOffset] - minimumSnr >= noise) {
+ for (int i = startOffset; i < stopIdx; i++) {
+ obs.addValue(waveform[i]);
+ spike.addValue(waveform[i]);
+ if (obs.getN() >= windowSize) {
+ if (obs.getMean() <= minimumSnr + noise) {
+ for (int j = i - obsWindow; j < i; j++) {
+ if (j > 0 && waveform[j] <= minimumSnr + noise) {
+ snrTimePick = j / sampleRate;
+ break;
+ }
+ }
+ break;
+
+ //TODO: 1.1 is a WAG on a data set, should try to find a better formulation that takes into account things like sample rate, window, and overall SNR
+ } else if (spike.getMean() > (obs.getMean() * 1.1d)) {
+ break;
+ } else {
+ snrTimePick = (i - spike.getN()) / sampleRate;
+ }
}
}
}
diff --git a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/MaxVelocityCalculator.java b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/MaxVelocityCalculator.java
index 5050ca27..16a9175e 100644
--- a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/MaxVelocityCalculator.java
+++ b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/MaxVelocityCalculator.java
@@ -51,7 +51,6 @@ public MaxVelocityCalculator(WaveformToTimeSeriesConverter converter) {
public Collection computeMaximumVelocity(List waveforms) {
return waveforms.stream().parallel().map(rawWaveform -> {
TimeSeries waveform = converter.convert(rawWaveform);
- double median = waveform.getMedian();
double distance = EModel.getDistanceWGS84(
rawWaveform.getEvent().getLatitude(),
rawWaveform.getEvent().getLongitude(),
@@ -77,7 +76,7 @@ public Collection computeMaximumVelocity(List
double min_length = 25.;
if (endtime.subtract(starttime).getEpochTime() < min_length) {
- log.info("Coda window length too short: {}", endtime.subtract(starttime).getEpochTime());
+ log.debug("Coda window length too short: {}", endtime.subtract(starttime).getEpochTime());
}
// peakS[0] time in seconds for
@@ -92,12 +91,14 @@ public Collection computeMaximumVelocity(List
double velocity = distance / peakS[0];
+ double noise = WaveformUtils.getNoiseFloor(rawWaveform.getSegment());
+
// the envelope noise is in log10 units.
- double snrPeak = peakS[1] - median;
+ double snrPeak = peakS[1] - noise;
return new PeakVelocityMeasurement().setWaveform(rawWaveform)
.setNoiseStartSecondsFromOrigin(0d)
.setNoiseEndSecondsFromOrigin(20d)
- .setNoiseLevel(median)
+ .setNoiseLevel(noise)
.setSnr(snrPeak)
.setVelocity(velocity)
.setDistance(distance)
diff --git a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/SpectraCalculator.java b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/SpectraCalculator.java
index 46ab4038..8014338c 100644
--- a/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/SpectraCalculator.java
+++ b/calibration-service/calibration-service-impl/src/main/java/gov/llnl/gnem/apps/coda/calibration/service/impl/processing/SpectraCalculator.java
@@ -22,8 +22,8 @@
import java.util.Map;
import java.util.Map.Entry;
import java.util.Objects;
-import java.util.TreeMap;
import java.util.concurrent.ThreadLocalRandom;
+import java.util.function.Function;
import java.util.stream.Collectors;
import org.apache.commons.math3.analysis.MultivariateFunction;
@@ -50,7 +50,6 @@
import gov.llnl.gnem.apps.coda.calibration.model.domain.SiteFrequencyBandParameters;
import gov.llnl.gnem.apps.coda.calibration.model.domain.Spectra;
import gov.llnl.gnem.apps.coda.calibration.model.domain.SpectraMeasurement;
-import gov.llnl.gnem.apps.coda.calibration.service.api.EndTimePicker;
import gov.llnl.gnem.apps.coda.calibration.service.api.MdacParametersFiService;
import gov.llnl.gnem.apps.coda.calibration.service.api.MdacParametersPsService;
import gov.llnl.gnem.apps.coda.common.model.domain.Event;
@@ -77,17 +76,15 @@ public class SpectraCalculator {
private WaveformToTimeSeriesConverter converter;
private SyntheticCodaModel syntheticCodaModel;
- private EndTimePicker endTimePicker;
private MdacCalculatorService mdacService;
private MdacParametersFiService mdacFiService;
private MdacParametersPsService mdacPsService;
@Autowired
- public SpectraCalculator(WaveformToTimeSeriesConverter converter, SyntheticCodaModel syntheticCodaModel, EndTimePicker endTimePicker, MdacCalculatorService mdacService,
- MdacParametersFiService mdacFiService, MdacParametersPsService mdacPsService) {
+ public SpectraCalculator(WaveformToTimeSeriesConverter converter, SyntheticCodaModel syntheticCodaModel, MdacCalculatorService mdacService, MdacParametersFiService mdacFiService,
+ MdacParametersPsService mdacPsService) {
this.converter = converter;
this.syntheticCodaModel = syntheticCodaModel;
- this.endTimePicker = endTimePicker;
this.mdacService = mdacService;
this.mdacFiService = mdacFiService;
this.mdacPsService = mdacPsService;
@@ -146,14 +143,15 @@ private SpectraMeasurement measureAmplitudeForSynthetic(SyntheticCoda synth, Map
siteCorrection = frequencyBandSiteParameterMap.get(frequencyBand).get(station).getSiteTerm();
}
- double eshCorrection = log10ESHcorrection(synth.getSourceWaveform().getLowFrequency(),
- synth.getSourceWaveform().getHighFrequency(),
- params.getS1(),
- params.getS2(),
- params.getXc(),
- params.getXt(),
- params.getQ(),
- distance);
+ double eshCorrection = log10ESHcorrection(
+ synth.getSourceWaveform().getLowFrequency(),
+ synth.getSourceWaveform().getHighFrequency(),
+ params.getS1(),
+ params.getS2(),
+ params.getXc(),
+ params.getXt(),
+ params.getQ(),
+ distance);
double minlength = params.getMinLength();
double maxlength = params.getMaxLength();
@@ -161,22 +159,18 @@ private SpectraMeasurement measureAmplitudeForSynthetic(SyntheticCoda synth, Map
TimeT endTime = null;
if (synth.getSourceWaveform().getAssociatedPicks() != null) {
- endTime = synth.getSourceWaveform().getAssociatedPicks().stream().filter(p -> PICK_TYPES.F.name().equalsIgnoreCase(p.getPickType())).findFirst().map(pick -> {
- double startpick = startTime.subtractD(originTime);
- TimeT result = originTime;
- result = pick.getPickTimeSecFromOrigin() > (startpick + maxlength) ? result.add(startpick + maxlength) : result.add(pick.getPickTimeSecFromOrigin());
- return result;
- }).orElse(null);
- }
-
- if (endTime == null && autoPickingEnabled) {
- endTime = new TimeT(endTimePicker.getEndTime(envSeis.getData(),
- envSeis.getSamprate(),
- startTime.getEpochTime(),
- envSeis.getIndexForTime(startTime.getEpochTime()),
- minlength,
- maxlength,
- params.getMinSnr()));
+ endTime = synth.getSourceWaveform()
+ .getAssociatedPicks()
+ .stream()
+ .filter(p -> p.getPickType() != null && PICK_TYPES.F.name().equalsIgnoreCase(p.getPickType().trim()))
+ .findFirst()
+ .map(pick -> {
+ double startpick = startTime.subtractD(originTime);
+ TimeT result = originTime;
+ result = pick.getPickTimeSecFromOrigin() > (startpick + maxlength) ? result.add(startpick + maxlength) : result.add(pick.getPickTimeSecFromOrigin());
+ return result;
+ })
+ .orElse(null);
}
if (endTime == null) {
@@ -227,7 +221,7 @@ private SpectraMeasurement measureAmplitudeForSynthetic(SyntheticCoda synth, Map
correctedAmp = pathCorrectedAmp + siteCorrection;
if (endTime.subtract(startTime).getEpochTime() < minlength) {
- log.info("Coda window length too short for envelope: {}", synth.getSourceWaveform());
+ log.debug("Coda window length too short for envelope: {}", synth.getSourceWaveform());
return null;
}
@@ -440,14 +434,15 @@ public Spectra computeFitSpectra(MeasuredMwParameters event, List
* UCRL-ID-146882
*
*/
- public List measureMws(final Map> evidMap, final PICK_TYPES selectedPhase) {
+ public List measureMws(final Map> evidMap, Map, Map>> eventWeights,
+ final PICK_TYPES selectedPhase) {
final MdacParametersFI mdacFi = mdacFiService.findFirst();
final MdacParametersPS mdacPs = mdacPsService.findMatchingPhase(selectedPhase.getPhase());
List measuredMws = new ArrayList<>();
for (Entry> entry : evidMap.entrySet()) {
Map measurements = entry.getValue();
- double[] MoMw = fitMw(measurements, selectedPhase, mdacFi, mdacPs);
+ double[] MoMw = fitMw(measurements, selectedPhase, mdacFi, mdacPs, eventWeights.get(entry.getKey()));
if (MoMw == null) {
log.warn("MoMw calculation returned null value");
continue;
@@ -469,10 +464,11 @@ public List measureMws(final Map measurements, final PICK_TYPES phase, final MdacParametersFI mdacFi, final MdacParametersPS mdacPs) {
+ public double[] fitMw(final Map measurements, final PICK_TYPES phase, final MdacParametersFI mdacFi, final MdacParametersPS mdacPs,
+ Function, Map> weightFunction) {
double[] result = new double[5];
- final Map weightMap = new TreeMap<>();
+ final Map frequencyBands = new HashMap<>();
final double minMW = 0.5;
final double maxMW = 8.0;
final double minMPA = 0.00001;
@@ -481,20 +477,16 @@ public double[] fitMw(final Map measurements,
for (Entry meas : measurements.entrySet()) {
double logAmplitude = meas.getValue().getMean();
- double lowFreq = meas.getKey().getLowFrequency();
- double highFreq = meas.getKey().getHighFrequency();
- double centerFreq = lowFreq + 0.25 * (highFreq - lowFreq);
- weightMap.put(centerFreq, 0.0);
if (logAmplitude > 0.0) {
+ double lowFreq = meas.getKey().getLowFrequency();
+ double highFreq = meas.getKey().getHighFrequency();
+ double centerFreq = (highFreq + lowFreq) / 2.0;
+ frequencyBands.put(centerFreq, logAmplitude);
dataCount++;
}
}
- int count = 0;
- for (Entry entry : weightMap.entrySet()) {
- entry.setValue(count == 0 || count == 2 ? 0.5 : count == 1 ? 1.0 : count == 3 ? 0.25 : 0.1);
- count++;
- }
+ final Map weightMap = weightFunction.apply(frequencyBands);
MultivariateFunction mdacFunction = new MultivariateFunction() {
@@ -509,7 +501,7 @@ public double value(double[] point) {
double logAmplitude = meas.getValue().getMean();
double lowFreq = meas.getKey().getLowFrequency();
double highFreq = meas.getKey().getHighFrequency();
- double centerFreq = lowFreq + 0.25 * (highFreq - lowFreq);
+ double centerFreq = (highFreq + lowFreq) / 2.0;
if (mdacPs != null) {
// Note this is in dyne-cm to match Kevin
@@ -528,14 +520,14 @@ public double value(double[] point) {
ConvergenceChecker convergenceChecker = new SimplePointChecker<>(0.001, 0.001, 100000);
CMAESOptimizer optimizer = new CMAESOptimizer(1000000, 0, true, 0, 10, new MersenneTwister(), true, convergenceChecker);
- PointValuePair optimizerResult = optimizer.optimize(new MaxEval(1000000),
- new ObjectiveFunction(mdacFunction),
- GoalType.MINIMIZE,
- new SimpleBounds(new double[] { minMW, minMPA }, new double[] { maxMW, maxMPA }),
- new InitialGuess(new double[] { ThreadLocalRandom.current().nextDouble(minMW, maxMW),
- ThreadLocalRandom.current().nextDouble(minMPA, maxMPA) }),
- new CMAESOptimizer.Sigma(new double[] { 0.05, 1.0 }),
- new CMAESOptimizer.PopulationSize(50));
+ PointValuePair optimizerResult = optimizer.optimize(
+ new MaxEval(1000000),
+ new ObjectiveFunction(mdacFunction),
+ GoalType.MINIMIZE,
+ new SimpleBounds(new double[] { minMW, minMPA }, new double[] { maxMW, maxMPA }),
+ new InitialGuess(new double[] { ThreadLocalRandom.current().nextDouble(minMW, maxMW), ThreadLocalRandom.current().nextDouble(minMPA, maxMPA) }),
+ new CMAESOptimizer.Sigma(new double[] { 0.05, 1.0 }),
+ new CMAESOptimizer.PopulationSize(50));
// converted back into dyne-cm to match Kevin's format
double testMw = optimizerResult.getPoint()[0];
diff --git a/calibration-standalone/pom.xml b/calibration-standalone/pom.xml
index 4c97c762..f6a28b9d 100644
--- a/calibration-standalone/pom.xml
+++ b/calibration-standalone/pom.xml
@@ -85,17 +85,7 @@