Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP Allow -E+n for modelnorm output #7358

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
Open

WIP Allow -E+n for modelnorm output #7358

wants to merge 9 commits into from

Conversation

PaulWessel
Copy link
Member

See forum post for the origin of this PR.

Needs much work still - such as matching what the formum post seems to do?

See forum post for the origin of this PR.
@PaulWessel PaulWessel added the new feature PR that implements a new feature or capability in GMT label Apr 4, 2023
@PaulWessel PaulWessel added this to the 6.5.0 milestone Apr 4, 2023
@PaulWessel PaulWessel self-assigned this Apr 4, 2023
@PaulWessel
Copy link
Member Author

Can any of @GenericMappingTools/core see whaT I am missing repeatedly?

@wasjabloch
Copy link

wasjabloch commented Apr 12, 2023

I have successfully extracted the model variance using this branch. In order to plot an L-curve to determine the optimal number of singular values to use, I computed the data variance from the output grids as (python syntax):

datnorm = sqrt(sum(((pre - obs) / uncert) ** 2))

with predicted data values pre, observed data values obs and 1-sigma uncertainties (as supplied in the 4th column via -W) uncert.

It would be convenient to also output the data norm as an additional column into misfitfile.

During the course of my work, it turned out that the L-curve is not as easy to interpret as I wished to, so I ended up computing Akaikes Information Criterion (AIC)

aic = 2*k - 2*ln(L)

with number of singular values used k and the model likelihood L (of k). I assumed Gaussian uncertainties and computed L as:

L = product( exp(-(obs - pred)**2 /(2 * uncert)**2) / sqrt(2 * pi * uncert**2) )

The extra steps needed (apart from calculating the model norm modnorm) can all be done manually outside GMT. However, writing out all the grids and extracting pred can be quite time consuming, so it would be really nice if these metrics could be provided in dry-run mode, that is in concert with -C+n. I am unsure, however, if -C+n currently does allow to compute pred.

One glitch in the current implementation: modnorm should be sqrt(modnorm).

@PaulWessel
Copy link
Member Author

Might you be able to submit a branch based off this one that implements what you describe?

@wasjabloch
Copy link

I could try getting started, but an important prerequisite would be to get the file output working. Until now, modnorm does not end up in misfitfile, but only 0s are written to the last column.

@PaulWessel
Copy link
Member Author

OK, I have a fixed a bunch of issues:

  1. modnorm is now sort (mod norm) and is written properly in the correct column.
  2. The column headers are now written to output file

@PaulWessel
Copy link
Member Author

I also merged in master so the branch is up to date with master plus the specific changes above.

@seisman seisman modified the milestones: 6.5.0, 6.6.0 Jan 7, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
new feature PR that implements a new feature or capability in GMT
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants