NSHM27 branch integration for upstream PR #269

src/main/java/org/opensha/nshmp/shaded/model/NshmErf.java

@@ -11,22 +11,28 @@
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Comparator;
+import java.util.EnumMap;
 import java.util.List;
+import java.util.Map;
 import java.util.Map.Entry;
 import java.util.Set;
 import java.util.stream.Collectors;
 import java.util.stream.IntStream;
 
+import org.apache.commons.numbers.core.Precision;
 import org.opensha.commons.data.TimeSpan;
 import org.opensha.sha.earthquake.AbstractERF;
 import org.opensha.sha.earthquake.ProbEqkSource;
 import org.opensha.sha.earthquake.param.IncludeBackgroundOption;
+import org.opensha.sha.imr.attenRelImpl.nshmp.NSHMP_GMM_Wrapper;
 import org.opensha.sha.util.TectonicRegionType;
 
+import com.google.common.base.Preconditions;
 import com.google.common.collect.Multimap;
 import com.google.common.collect.MultimapBuilder;
 
 import org.opensha.nshmp.shaded.data.NshmpIndexing;
+import org.opensha.nshmp.shaded.gmm.NshmpGmm;
 import org.opensha.nshmp.shaded.model.NshmpSystemRuptureSet.SystemRupture;
 import org.opensha.nshmp.shaded.tree.NshmpBranch;
 
@@ -43,6 +49,7 @@ public class NshmErf extends AbstractERF {
 	private final boolean volcanic;
 	private final boolean grid;
 	private final boolean faults;
+	private Set<TectonicRegionType> trts;
 
 	public NshmErf(Path path, Set<TectonicRegionType> trts,
 			IncludeBackgroundOption gridOption) {
@@ -52,6 +59,7 @@ public NshmErf(Path path, Set<TectonicRegionType> trts,
 	public NshmErf(NshmpHazardModel model, Set<TectonicRegionType> trts,
 			IncludeBackgroundOption gridOption) {
 		this.model = model;
+		this.trts = trts;
 		allSources = new ArrayList<>();
 		sourceMap = MultimapBuilder
 				.enumKeys(TectonicRegionType.class)
@@ -69,6 +77,8 @@ public NshmErf(NshmpHazardModel model, Set<TectonicRegionType> trts,
 
 		init();
 	}
+
+	private static final boolean PROCESS_TRT_OVERRIDES = true;
 
 	private void init() {
 
@@ -79,37 +89,134 @@ private void init() {
 
 		// nshmp-haz initializers
 		Multimap<NshmpTectonicSetting, NshmpSourceTree> trees = model.trees();
+		if (PROCESS_TRT_OVERRIDES) {
+			List<NshmpTRTBranch> trtMappedRuptureSets = parseTRTOverrides(trees);
+			for (NshmpTRTBranch trtBranch : trtMappedRuptureSets) {
+				if (!trts.isEmpty() && !trts.contains(trtBranch.trt))
+					continue;
+				List<NshmSource> sources = initBranch(trtBranch.branch);
+				sources.forEach(s -> s.setTectonicRegionType(trtBranch.trt));
+				allSources.addAll(sources);
+				sourceMap.putAll(trtBranch.trt, sources);
+			}
+			allSources.sort(new Comparator<NshmSource>() {
+				@Override
+				public int compare(NshmSource o1, NshmSource o2) {
+					return Integer.compare(o1.getNSHM_ID(), o2.getNSHM_ID());
+				}
+			});
+		} else {
+			for (Entry<NshmpTectonicSetting, NshmpSourceTree> entry : trees.entries()) {
+				NshmpTectonicSetting setting = entry.getKey();
+				NshmpSourceTree tree = entry.getValue();
+				NshmpSourceType type = tree.type();
+
+				TectonicRegionType trt = NshmUtil.tectonicSettingToType(setting, type);
+				if (!trts.isEmpty() && !trts.contains(trt))
+					continue;
+				List<NshmSource> sources = initTree(tree);
+				sources.forEach(s -> s.setTectonicRegionType(trt));
+				allSources.addAll(sources);
+				sourceMap.putAll(trt, sources);
+			}
+		}
+	}
+
+	private record NshmpTRTBranch(TectonicRegionType trt, NshmpBranch<NshmpRuptureSet> branch) {}
+
+	/**
+	 * Map each branch to a tectonic region type, also processing any GMM-overrides that Peter might sneak in there
+	 * to change the TRT within a different NshmpTectonicSetting
+	 * @param trees
+	 * @return
+	 */
+	private static List<NshmpTRTBranch> parseTRTOverrides(
+			Multimap<NshmpTectonicSetting, NshmpSourceTree> trees) {
+		List<NshmpTRTBranch> ret = new ArrayList<>();
+
 		for (Entry<NshmpTectonicSetting, NshmpSourceTree> entry : trees.entries()) {
-
 			NshmpTectonicSetting setting = entry.getKey();
 			NshmpSourceTree tree = entry.getValue();
 			NshmpSourceType type = tree.type();
-
-			if (setting == NshmpTectonicSetting.SUBDUCTION) {
-				if ((type == NshmpSourceType.INTERFACE || type == NshmpSourceType.INTERFACE_CLUSTER
-						|| type == NshmpSourceType.INTERFACE_GRID || type == NshmpSourceType.INTERFACE_SYSTEM) && !subInterface) {
+			TectonicRegionType origTRT = NshmUtil.tectonicSettingToType(setting, type);
+
+			System.out.println("Processing tree "+tree.id()+". "+tree.name()+" ("+origTRT.name()+")");
+
+			for (NshmpBranch<NshmpRuptureSet> branch : tree) {
+				if (origTRT == VOLCANIC) {
+					// copy over as is, no current GMM type for volcanic
+					ret.add(new NshmpTRTBranch(origTRT, branch));
 					continue;
 				}
-				if ((type == NshmpSourceType.SLAB || type == NshmpSourceType.INTRASLAB_GRID) && !subSlab) {
-					continue;
+				NshmpRuptureSet rs = branch.value();
+				System.out.println("\tProcessing RS "+rs.id()+". "+rs.name()+" ("+origTRT.name()+")");
+				NshmpGmmTree gmmTree = rs.gmmTree();
+				double weightOrig = 0d;
+				double weightOverrideSum = 0d;
+				Map<TectonicRegionType, Double> overrideWeights = new EnumMap<>(TectonicRegionType.class);
+				for (NshmpBranch<NshmpGmm> gmm : gmmTree.tree()) {
+					TectonicRegionType gmmTRT = NSHMP_GMM_Wrapper.trtForType(gmm.value().type());
+					if (gmmTRT != origTRT) {
+						// we have a GMM TRT override
+						weightOverrideSum += gmm.weight();
+						if (overrideWeights.containsKey(gmmTRT))
+							overrideWeights.put(gmmTRT, overrideWeights.get(gmmTRT)+gmm.weight());
+						else
+							overrideWeights.put(gmmTRT, gmm.weight());
+					} else {
+						weightOrig += gmm.weight();
+					}
+				}
+				if (weightOverrideSum > 0d) {
+					// we have overrides to process
+					Preconditions.checkState(Precision.equals(weightOverrideSum + weightOrig, 1d, 1e-4));
+					System.err.println("Detected a GMM TRT override for source "+rs.id()+". "+rs.name()+" with original TRT="+origTRT.name());
+					if (weightOrig > 0d) {
+						System.err.println("\t"+origTRT.name()+":\t"+weightOrig);
+						ret.add(new NshmpTRTBranch(origTRT, new WeightScaledBranch(branch, weightOrig, origTRT)));
+					}
+					for (TectonicRegionType trt : overrideWeights.keySet()) {
+						double weight = overrideWeights.get(trt);
+						System.err.println("\t"+trt.name()+":\t"+weight);
+						ret.add(new NshmpTRTBranch(trt, new WeightScaledBranch(branch, weight, trt)));
+					}
+				} else {
+					// copy over as is
+					ret.add(new NshmpTRTBranch(origTRT, branch));
 				}
 			}
-			if (setting == NshmpTectonicSetting.STABLE_CRUST && !stableCrust) {
-				continue;
-			}
-			if (setting == NshmpTectonicSetting.ACTIVE_CRUST && !activeCrust) {
-				continue;
-			}
-			if (setting == NshmpTectonicSetting.VOLCANIC && !volcanic) {
-				continue;
-			}
+		}
+
+		return ret;
+	}
+
+	private static class WeightScaledBranch implements NshmpBranch<NshmpRuptureSet> {
+
+		private NshmpBranch<NshmpRuptureSet> upstream;
+		private double scale;
+		private TectonicRegionType trt;
+
+		public WeightScaledBranch(NshmpBranch<NshmpRuptureSet> upstream, double scale, TectonicRegionType trt) {
+			this.upstream = upstream;
+			this.scale = scale;
+			this.trt = trt;
+		}
+
+		@Override
+		public String id() {
+			return upstream.id()+"-"+trt.name();
+		}
+
+		@Override
+		public NshmpRuptureSet value() {
+			return upstream.value();
+		}
 
-			TectonicRegionType trt = NshmUtil.tectonicSettingToType(setting, type);
-			List<NshmSource> sources = initTree(tree);
-			sources.forEach(s -> s.setTectonicRegionType(trt));
-			allSources.addAll(sources);
-			sourceMap.putAll(trt, sources);
+		@Override
+		public double weight() {
+			return upstream.weight()*scale;
 		}
+
 	}
 
 	public List<NshmSource> allSources() {
@@ -152,6 +259,14 @@ public int compare(NshmSource o1, NshmSource o2) {
 		return sources;
 	}
 
+	private List<NshmSource> initBranch(NshmpBranch<NshmpRuptureSet> branch) {
+		List<NshmSource> sources = new ArrayList<>();
+		double duration = getTimeSpan().getDuration();
+		sources.addAll(sourcesFromBranch(branch, duration));
+
+		return sources;
+	}
+
 	private List<NshmSource> sourcesFromBranch(
 			NshmpBranch<NshmpRuptureSet> branch,
 			double duration) {

src/main/java/org/opensha/nshmp/shaded/model/NshmSource.java

@@ -49,11 +49,7 @@ public RuptureSurface getSourceSurface() {
 		throw new UnsupportedOperationException();
 	}
 
-	public Object getDelegate() {
-		return delegate;
-	}
-
-	static class Fault extends NshmSource<NshmpIterableRuptureSet> {
+	public static class Fault extends NshmSource<NshmpIterableRuptureSet> {
 
 		private final List<NshmRupture> ruptures;
 
@@ -99,7 +95,7 @@ public void setDuration(double duration) {
 		}
 	}
 
-	static class Point extends NshmSource<NshmpGridSource> {
+	public static class Point extends NshmSource<NshmpGridSource> {
 
 		final double weight;
 		double duration;
@@ -142,7 +138,7 @@ public void setDuration(double duration) {
 		}
 	}
 
-	static class System extends NshmSource<SystemRupture> {
+	public static class System extends NshmSource<SystemRupture> {
 
 		final List<NshmSurface> surfaces;
 		final NshmRupture rupture;

src/main/java/org/opensha/nshmp/shaded/model/NshmSurface.java

@@ -11,6 +11,8 @@
 import org.opensha.sha.faultSurface.RuptureSurface;
 import org.opensha.sha.faultSurface.cache.CacheEnabledSurface;
 import org.opensha.sha.faultSurface.cache.SurfaceDistances;
+import org.opensha.sha.faultSurface.utils.ptSrcCorr.PointSourceDistanceCorrection;
+import org.opensha.sha.util.TectonicRegionType;
 
 import org.opensha.nshmp.shaded.fault.surface.NshmpDefaultGriddedSurface;
 import org.opensha.nshmp.shaded.fault.surface.NshmpGriddedSurface;
@@ -36,6 +38,64 @@ public class NshmSurface implements CacheEnabledSurface {
 	public NshmSurface(org.opensha.nshmp.shaded.fault.surface.NshmpRuptureSurface delegate) {
 		this.delegate = delegate;
 	}
+
+	/**
+	 * This creates a point surface that will work with existing OpenSHA point-source optimizations
+	 * @param delegate
+	 * @return
+	 */
+	public static org.opensha.sha.faultSurface.PointSurface buildPointSurface(
+			org.opensha.nshmp.shaded.fault.surface.NshmpRuptureSurface delegate) {
+		// this is the point surface
+		double len = 0d;
+		try {
+			len = delegate.length();
+		} catch (Exception e) {}
+		org.opensha.sha.faultSurface.PointSurface surf = new org.opensha.sha.faultSurface.PointSurface(
+				NshmUtil.toOpenShaLocation(delegate.centroid()), delegate.dip(), delegate.depth(),
+				delegate.depth() + delegate.width()*Math.sin(Math.toRadians(delegate.dip())), len);
+		return new org.opensha.sha.faultSurface.PointSurface.DistanceCorrecting(
+				surf, new DelegatePointSourceCorrection(delegate), null, Double.NaN);
+	}
+
+	/**
+	 * Delegate point source correction that passes through to NSHMP-haz. This is required for point source
+	 * optimizations to work with wrapped point sources. As part of that, equals/hashCode have to be constant
+	 * for all delegate corrections.
+	 */
+	private static class DelegatePointSourceCorrection implements PointSourceDistanceCorrection.Single {
+
+		private org.opensha.nshmp.shaded.fault.surface.NshmpRuptureSurface delegate;
+
+		private DelegatePointSourceCorrection(org.opensha.nshmp.shaded.fault.surface.NshmpRuptureSurface delegate) {
+			this.delegate = delegate;
+		}
+
+		@Override
+		public SurfaceDistances getCorrectedDistance(Location location, org.opensha.sha.faultSurface.PointSurface surf,
+				TectonicRegionType trt, double mag, double horzDist) {
+			NshmpDistance distance = delegate.distanceTo(NshmUtil.fromOpenShaLocation(location));
+			return new SurfaceDistances.Precomputed(location, distance.rRup, distance.rJB, distance.rX);
+		}
+
+		private static final int hashCode = DelegatePointSourceCorrection.class.hashCode();
+
+		@Override
+		public int hashCode() {
+			return hashCode;
+		}
+
+		@Override
+		public boolean equals(Object obj) {
+//			return obj instanceof DelegatePointSourceCorrection;
+			if (!(obj instanceof DelegatePointSourceCorrection))
+				return false;
+			DelegatePointSourceCorrection other = (DelegatePointSourceCorrection)obj;
+//			System.out.println("Delegate classes:\t"+delegate.getClass()+"\t"+other.delegate.getClass());
+			return other.delegate.getClass().equals(delegate.getClass());
+		}
+
+	}
 
 	// return nshmp-haz rupture centroid as OpenSHA location for
 	// use in computing min distance to a fault system subsection