Skip to content

Commit

Permalink
improvements to residual plots; XYZ now with log axis
Browse files Browse the repository at this point in the history
  • Loading branch information
Kevin Milner committed Nov 7, 2024
1 parent ee413ad commit cf9a330
Show file tree
Hide file tree
Showing 4 changed files with 82 additions and 45 deletions.
34 changes: 29 additions & 5 deletions src/main/java/scratch/kevin/quantum/ToyRuptureSetBuilder.java
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,9 @@
import java.io.File;
import java.io.FileOutputStream;
import java.io.IOException;
import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Date;
import java.util.List;

import org.opensha.commons.data.CSVWriter;
Expand Down Expand Up @@ -35,10 +37,27 @@ public class ToyRuptureSetBuilder {
public static void main(String[] args) throws IOException {
boolean rebuildRupSet = false;

File outputDir = new File("/home/kevin/markdown/inversions/2024_10_31-quantum_test_rup_set_problem");
String dirName = new SimpleDateFormat("yyyy_MM_dd").format(new Date());
// String dirName = "2024_11_06";

dirName += "-quantum_test_rup_set_problem";

// String[] nameSearches = { "Likely" };
// dirName += "-70_sects";

// String[] nameSearches = { "Cedar", "Mountain", "Mahogany" };
// dirName += "-124_sects";

String[] nameSearches = { "San", "Felipe", "proxy" };
dirName += "-213_sects";

System.out.println(dirName);

File outputDir = new File("/home/kevin/markdown/inversions/"+dirName);
Preconditions.checkState(outputDir.exists() || outputDir.mkdir());

File rsFile = new File(outputDir, "orig_rup_set.zip");
File exhaustiveFile = new File(outputDir, "exhaustive_rup_set.zip");
FaultSystemRupSet origRupSet;
if (rsFile.exists() && !rebuildRupSet) {
origRupSet = FaultSystemRupSet.load(rsFile);
Expand All @@ -51,7 +70,7 @@ public static void main(String[] args) throws IOException {
System.out.println("Building connectivity clusters");
ConnectivityClusters clusters = ConnectivityClusters.build(origRupSet);

int targetParentID = FaultSectionUtils.findParentSectionID(origRupSet.getFaultSectionDataList(), "Likely");
int targetParentID = FaultSectionUtils.findParentSectionID(origRupSet.getFaultSectionDataList(), nameSearches);
FaultSection targetSect = null;
for (FaultSection sect : origRupSet.getFaultSectionDataList()) {
if (sect.getParentSectionId() == targetParentID) {
Expand Down Expand Up @@ -116,7 +135,14 @@ public static void main(String[] args) throws IOException {
// if true, allow splays (using default settings)
rsConfig.setSplays(true); // enabled

FaultSystemRupSet rupSet = rsConfig.build(FaultSysTools.defaultNumThreads());
FaultSystemRupSet rupSet;
if (exhaustiveFile.exists() && !rebuildRupSet) {
rupSet = FaultSystemRupSet.load(exhaustiveFile);
} else {
rupSet = rsConfig.build(FaultSysTools.defaultNumThreads());

rupSet.write(exhaustiveFile);
}
System.out.println("Original had "+origRupSet.getNumRuptures());
System.out.println("New has "+rupSet.getNumRuptures());

Expand Down Expand Up @@ -171,8 +197,6 @@ public static void main(String[] args) throws IOException {
// now write the actual ruptures
writeRuptures(rupSet.getSectionIndicesForAllRups(), new File(outputDir, "orig_ruptures.csv"), subSects.size());
writeRuptures(rupSet.getSectionIndicesForAllRups(), new File(outputDir, "exhaustive_ruptures.csv"), subSects.size());

rupSet.write(new File(outputDir, "exhaustive_rup_set.zip"));
}

private static void writeMatrix(double[][] mat, File file) throws IOException {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -634,14 +634,18 @@ public void generateGMPE_Page(File outputDir, List<String> headerLines, AttenRel
lines.addAll(table.build());
lines.add("");

if (!siteBundles.isEmpty() && sites.size() > 1) {
System.out.println("All sites aggregated");
if (!siteBundles.isEmpty() && this.sites.size() > 1) {
System.out.println("Bundled sites aggregated");
lines.add("### Bundled z-scores");
lines.add(topLink); lines.add("");

for (int s=0; s<siteBundles.size(); s++) {
String name = siteBundleNames.get(s);
String prefix = "site_bundle_"+name.replaceAll("\\W+", "_")+"_std_norm";
String sanitizedName = name.replace("<=", "_le_").replace("<", "_less_")
.replace(">=", "_ge_").replace(">", "_greater_").replaceAll("\\W+", "_");
while (sanitizedName.contains("__"))
sanitizedName = sanitizedName.replace("__", "_");
String prefix = "site_bundle_"+sanitizedName+"_std_norm";

System.out.println("Processing site bundle: "+name);

Expand Down Expand Up @@ -1054,9 +1058,12 @@ public void generateGMPE_Page(File outputDir, List<String> headerLines, AttenRel
List<ResidualType> residualTypes = new ArrayList<>();
residualTypes.add(ResidualType.MAG);
ScalarIMR tempGMPE = checkOutGMPE(gmpeRef);
if (tempGMPE.getPropagationEffectParams().containsParameter(DistanceJBParameter.NAME))
ScalarIMR subGMPE = tempGMPE instanceof MultiIMR_Averaged_AttenRel ? ((MultiIMR_Averaged_AttenRel)tempGMPE).getIMRs().get(0) : null;
if (tempGMPE.getPropagationEffectParams().containsParameter(DistanceJBParameter.NAME)
|| (subGMPE != null && subGMPE.getPropagationEffectParams().containsParameter(DistanceJBParameter.NAME)))
residualTypes.add(ResidualType.DIST_JB);
if (tempGMPE.getPropagationEffectParams().containsParameter(DistanceRupParameter.NAME))
if (tempGMPE.getPropagationEffectParams().containsParameter(DistanceRupParameter.NAME)
|| (subGMPE != null && subGMPE.getPropagationEffectParams().containsParameter(DistanceRupParameter.NAME)))
residualTypes.add(ResidualType.DIST_RUP);
checkInGMPE(gmpeRef, tempGMPE);
sites = new ArrayList<>(this.sites); // we previously added a null element for all sites above
Expand Down Expand Up @@ -1389,7 +1396,7 @@ public Double getValue(RuptureComparison<?> comp, Site site) {
return comp.getDistanceJB(site);
}
},
VS30("Vs30", "m/s", Vs30_Param.NAME, true, 200d, Double.NaN) {
VS30("Vs30", "m/s", Vs30_Param.NAME, false, 200d, Double.NaN) {
@Override
public Double getValue(RuptureComparison<?> comp, Site site) {
return site.getParameter(Double.class, Vs30_Param.NAME).getValue();
Expand Down
64 changes: 35 additions & 29 deletions src/main/java/scratch/kevin/simCompare/ResidualScatterPlot.java
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,7 @@ else if (minX > myMinX+halfDelta)
else
xRange = new Range(minX, maxX);

SimpleRegression linearRegression = getRegression(xy, logX, xRange);
SimpleRegression linearRegression = getRegression(xy, logX, xRange, true, true);
regressionIntercept = linearRegression.getIntercept();
regressionLinearSlope = linearRegression.getSlope();
System.out.println("Linear regression yInt: "+linearRegression.getIntercept());
Expand All @@ -168,7 +168,7 @@ else if (minX > myMinX+halfDelta)
Range regressRange = new Range(binCenter - 0.5*deltaX, binCenter + 0.5*deltaX);
if (logX)
regressRange = new Range(Math.pow(10, regressRange.getLowerBound()), Math.pow(10, regressRange.getUpperBound()));
regressionBinned[i] = getRegressionFit(xy, logX, regressRange);
regressionBinned[i] = getRegressionFit(xy, logX, regressRange, i == 0, i == numRegressionBins-1);
if (regressionBinned[i] != null)
dataSigmaBinned[i] = calcSigma(regressionBinned[i], xy);
}
Expand Down Expand Up @@ -239,11 +239,11 @@ public void plotScatter(File outputDir, String prefix) throws IOException {
funcs.add(zeroLine);
chars.add(new PlotCurveCharacterstics(PlotLineType.DASHED, 3f, Color.DARK_GRAY));

addResidualFuncs(funcs, chars, true);
addResidualFuncs(funcs, chars);

PlotSpec spec = new PlotSpec(funcs, chars, title, xAxisLabel, residualLabel);
spec.setLegendVisible(true);
List<XYTextAnnotation> anns = buildAnnotations(false);
List<XYTextAnnotation> anns = buildAnnotations();
if (isFiltered) {
// add filtered note
Font font = new Font(Font.SANS_SERIF, Font.BOLD, 20);
Expand Down Expand Up @@ -355,8 +355,9 @@ else if (!Doubles.isFinite(minZ))
ArrayList<PlotCurveCharacterstics> chars = new ArrayList<>();

Range xRange = this.xRange;
if (logX)
xRange = new Range(Math.log10(xRange.getLowerBound()), Math.log10(xRange.getUpperBound()));
// now plotting with log axis; XYZ is still in log spacing, but range in linear
// if (logX)
// xRange = new Range(Math.log10(xRange.getLowerBound()), Math.log10(xRange.getUpperBound()));

DefaultXY_DataSet zeroLine = new DefaultXY_DataSet();
zeroLine.set(xRange.getLowerBound(), 0d);
Expand All @@ -365,7 +366,7 @@ else if (!Doubles.isFinite(minZ))
funcs.add(0, zeroLine);
chars.add(0, new PlotCurveCharacterstics(PlotLineType.DASHED, 3f, Color.DARK_GRAY));

addResidualFuncs(funcs, chars, false);
addResidualFuncs(funcs, chars);

// add custom y axis guides
for (int i=(int)Math.round(yRange.getLowerBound()); i<=(int)Math.round(yRange.getUpperBound()); i++) {
Expand All @@ -381,7 +382,7 @@ else if (!Doubles.isFinite(minZ))

xyzSpec.setXYElems(funcs);
xyzSpec.setXYChars(chars);
xyzSpec.setPlotAnnotations(buildAnnotations(true));
xyzSpec.setPlotAnnotations(buildAnnotations());

PlotPreferences plotPrefs = PlotPreferences.getDefault();
plotPrefs.setTickLabelFontSize(18);
Expand All @@ -390,7 +391,7 @@ else if (!Doubles.isFinite(minZ))
plotPrefs.setBackgroundColor(Color.WHITE);
HeadlessGraphPanel xyzGP = new HeadlessGraphPanel(plotPrefs);
System.out.println("Plotting with xRange: "+xRange);
xyzGP.drawGraphPanel(xyzSpec, false, false, xRange, yRange);
xyzGP.drawGraphPanel(xyzSpec, logX, false, xRange, yRange);
// write plot
xyzGP.getYAxis().setStandardTickUnits(getYTickUnits());
xyzGP.getChartPanel().setSize(plotWidth, plotWidth-100);
Expand All @@ -399,7 +400,7 @@ else if (!Doubles.isFinite(minZ))
xyzGP.saveAsPDF(new File(outputDir, prefix+".pdf").getAbsolutePath());
}

private List<XYTextAnnotation> buildAnnotations(boolean doLog) {
private List<XYTextAnnotation> buildAnnotations() {
List<XYTextAnnotation> anns = new ArrayList<>();

double yPos = yRange.getUpperBound() * 0.95;
Expand All @@ -412,7 +413,8 @@ private List<XYTextAnnotation> buildAnnotations(boolean doLog) {
xMin = Math.log10(xMin);
}
double xPos = xMin + 0.95*(xMax - xMin);
if (logX && !doLog)
// if (logX && !doLog)
if (logX)
xPos = Math.pow(10, xPos);

Font font = new Font(Font.SANS_SERIF, Font.BOLD, 24);
Expand Down Expand Up @@ -448,7 +450,8 @@ private List<XYTextAnnotation> buildAnnotations(boolean doLog) {
if (plotLinearResidual) {
yPos = yRange.getUpperBound() * 0.95;
xPos = xMin + 0.05*(xMax - xMin);
if (logX && !doLog)
// if (logX && !doLog)
if (logX)
xPos = Math.pow(10, xPos);

ann = new XYTextAnnotation("Linear LSQ Fit", xPos, yPos);
Expand Down Expand Up @@ -478,7 +481,7 @@ private List<XYTextAnnotation> buildAnnotations(boolean doLog) {
private static final DecimalFormat annDF = new DecimalFormat("0.00");
private static final DecimalFormat regressionDF = new DecimalFormat("0.00E0");

private void addResidualFuncs(List<XY_DataSet> funcs, List<PlotCurveCharacterstics> chars, boolean scatter) {
private void addResidualFuncs(List<XY_DataSet> funcs, List<PlotCurveCharacterstics> chars) {
Color linearColor = Color.BLACK;
Color binnedColor = Color.BLACK;
Color gmpeColor = new Color(35, 72, 132);
Expand All @@ -489,9 +492,9 @@ private void addResidualFuncs(List<XY_DataSet> funcs, List<PlotCurveCharactersti
// linear fit
if (plotLinearResidual) {
regressionLinear.setName("Linear LS Fit");
if (!scatter && logX)
funcs.add(getInLogX(regressionLinear));
else
// if (!scatter && logX)
// funcs.add(getInLogX(regressionLinear));
// else
funcs.add(regressionLinear);
chars.add(new PlotCurveCharacterstics(PlotLineType.SOLID, residualThickness, linearColor));
}
Expand Down Expand Up @@ -530,11 +533,11 @@ private void addResidualFuncs(List<XY_DataSet> funcs, List<PlotCurveCharactersti
tickEnd = binMiddle + binFractEachSide*delta;
}

if (!scatter && logX) {
binMiddle = Math.log10(binMiddle);
tickStart = Math.log10(tickStart);
tickEnd = Math.log10(tickEnd);
}
// if (!scatter && logX) {
// binMiddle = Math.log10(binMiddle);
// tickStart = Math.log10(tickStart);
// tickEnd = Math.log10(tickEnd);
// }

binnedPoints.set(binMiddle, binY);
addSigmaFuncs(funcs, chars, binMiddle, binY, tickStart, tickEnd, dataSigmaBinned[i],
Expand Down Expand Up @@ -608,23 +611,26 @@ private static TickUnits getYTickUnits() {
return tus;
}

private static SimpleRegression getRegression(XY_DataSet xy, boolean logX, Range xRange) {
private static SimpleRegression getRegression(XY_DataSet xy, boolean logX, Range xRange, boolean includeBelow, boolean includeAbove) {
// regression
SimpleRegression regression = new SimpleRegression();
for (Point2D pt : xy) {
if (!xRange.contains(pt.getX()))
continue;
if (logX)
pt = new Point2D.Double(Math.log10(pt.getX()), pt.getY());
regression.addData(pt.getX(), pt.getY());
boolean contains = xRange.contains(pt.getX());
contains |= includeBelow && pt.getX() <= xRange.getLowerBound();
contains |= includeAbove && pt.getX() >= xRange.getUpperBound();
if (contains) {
if (logX)
pt = new Point2D.Double(Math.log10(pt.getX()), pt.getY());
regression.addData(pt.getX(), pt.getY());
}
}
if (regression.getN() == 0l)
return null;
return regression;
}

private static ArbitrarilyDiscretizedFunc getRegressionFit(XY_DataSet xy, boolean logX, Range xRange) {
SimpleRegression regression = getRegression(xy, logX, xRange);
private static ArbitrarilyDiscretizedFunc getRegressionFit(XY_DataSet xy, boolean logX, Range xRange, boolean includeBelow, boolean includeAbove) {
SimpleRegression regression = getRegression(xy, logX, xRange, includeBelow, includeAbove);
if (regression == null)
return null;
return getRegressionFit(regression, xy, logX, xRange);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -636,16 +636,16 @@ public static void main(String[] args) throws ZipException, IOException {
// RuptureIdentifier loadIden = FocalMechIden.builder().forRake(-105, -75).forDip(35, 55).build();
// String loadIdenPrefix = "mech_normal";

boolean replotScatters = true;
boolean replotZScores = true;
boolean replotCurves = true;
boolean replotResiduals = true;

// boolean replotScatters = true;
// boolean replotZScores = true;
// boolean replotCurves = true;
// boolean replotResiduals = true;

boolean replotScatters = false;
boolean replotZScores = false;
boolean replotCurves = false;
boolean replotResiduals = true;

if (timeScale != 1d)
doRotD = false;

Expand Down

0 comments on commit cf9a330

Please sign in to comment.