Open yugao0912 opened 2 months ago
ahaselsteiner
commented
2 months ago Hi @yugao0912,
This looks indeed strange. Can you post the script that you used?
Yes, you can take the same model and provide it as a parameter to IFORMContour.
This example script in which we compute different contours for the same model might help you: https://github.com/virocon-organization/virocon/blob/main/examples/compare_contour_methods.py
-Andreas
yugao0912
commented
2 months ago Sure! I did not modify any other codes except for the data file. It seems that I cannot upload the python file. I have attached the data file. I am unsure whether you could see them.
Additionally, may I confirm I just need to change the codes to IFORMCoutour here:
[cid:c9de9864-0317-4f21-9526-18e2e8c3bf6e]
The codes please see below:
""" A comprehensive example that shows the whole workflow of 1) Loading data 2) Defining the model structure for a joint distribution 3) Estimating the parameter values of the model ("Fitting") 4) Computing a 3D environmental contour
"""
import numpy as np import matplotlib.pyplot as plt import pandas as pd
from virocon import ( GlobalHierarchicalModel, LogNormalDistribution, ExponentiatedWeibullDistribution, DependenceFunction, WidthOfIntervalSlicer, HighestDensityContour, plot_marginal_quantiles, plot_dependence_functions, )
data = pd.read_csv("datasets/coastDat2_oneyear.csv", sep=";", skipinitialspace=True) data.index = pd.to_datetime(data.pop(data.columns[0]), format="%Y-%m-%d-%H")
def _power3(x, a, b, c): return a + b * x**c
def _exp3(x, a, b, c): return a + b np.exp(c x)
def _alpha3(x, a, b, c, d_of_x): return (a + b * xc) / 2.0445 (1 / d_of_x(x))
def _logistics4(x, a=1, b=1, c=-1, d=1): return a + b / (1 + np.exp(c * (x - d)))
bounds = [(0, None), (0, None), (None, None)] logistics_bounds = [(0, None), (0, None), (None, 0), (0, None)]
power3 = DependenceFunction(_power3, bounds, latex="$a + b x^c$") exp3 = DependenceFunction(_exp3, bounds, latex="$a + b \exp(c x)$") logistics4 = DependenceFunction( _logistics4, logistics_bounds, weights=lambda x, y: y, latex="$a + b / (1 + \exp[c (x -d)])$", ) alpha3 = DependenceFunction( _alpha3, bounds, d_of_x=logistics4, weights=lambda x, y: y, latex="$(a + b * x^c) / 2.0445^{1 / F()}$", )
dist_description_0 = { "distribution": ExponentiatedWeibullDistribution(), "intervals": WidthOfIntervalSlicer(2, min_n_points=50), }
dist_description_1 = { "distribution": ExponentiatedWeibullDistribution(f_delta=5), "intervals": WidthOfIntervalSlicer(0.5), "conditional_on": 0, "parameters": { "alpha": alpha3, "beta": logistics4, }, }
dist_description_2 = { "distribution": LogNormalDistribution(), "conditional_on": 1, "parameters": {"mu": power3, "sigma": exp3}, }
model = GlobalHierarchicalModel( [dist_description_0, dist_description_1, dist_description_2] )
semantics = { "names": ["Wind speed", "Significant wave height", "Zero-up-crossing period"], "symbols": ["V", "H_s", "T_z"], "units": ["m/s", "m", "s"], }
model.fit(data)
print(model)
fig1, axs = plt.subplots(1, 3, figsize=[20, 7], dpi=300) plot_marginal_quantiles(model, data, semantics, axes=axs) fig2, axs = plt.subplots(1, 4, figsize=[20, 7], dpi=300) plot_dependence_functions(model, semantics, axes=axs)
state_duration = 3 # hours return_period = 50 # years alpha = state_duration / (return_period 365.25 24) HDC = HighestDensityContour(model, alpha, limits=[(0, 50), (0, 25), (0, 25)])
rng = np.random.default_rng(42) n_contour_points = len(HDC.coordinates) random_points = rng.choice( n_contour_points, int(0.05 * n_contour_points), replace=False )
Xs = HDC.coordinates[random_points, 0] Ys = HDC.coordinates[random_points, 1] Zs = HDC.coordinates[random_points, 2]
fig = plt.figure(dpi=300, figsize=[15, 7]) ax = fig.add_subplot(111, projection="3d") ax.scatter(Xs, Ys, Zs, c="#004488") ax.set_xlabel("Wind speed (m/s)") ax.set_ylabel("Significant wave height (m)") ax.set_zlabel("Zero-up-crossing period (s)")
Best wishes, Yu
发件人: Andreas Haselsteiner @.> 发送时间: 2024年7月17日 20:36 收件人: virocon-organization/virocon @.> 抄送: Yu Gao @.>; Mention @.> 主题: Re: [virocon-organization/virocon] Problems with the results from the 3D_contour.py (Issue #219)
Hi @yugao0912https://github.com/yugao0912,
This looks indeed strange. Can you post the script that you used?
Yes, you can take the same model and provide it as a parameter to IFORMContour. This example script in which we compute different contours for the same model might help you: https://github.com/virocon-organization/virocon/blob/main/examples/compare_contour_methods.py
― Reply to this email directly, view it on GitHubhttps://github.com/virocon-organization/virocon/issues/219#issuecomment-2234092890, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A56V74GCRF35XQSCHH3N2RDZM3BUDAVCNFSM6AAAAABLBESGTSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMZUGA4TEOBZGA. You are receiving this because you were mentioned.Message ID: @.***>
ahaselsteiner
commented
2 months ago Thanks. The data file is unfortunately not uploaded. Could you give the upload another try? I'd like to reproduce the problem on my machine. The figures of the model fit look good to me such that I would expect that contour calculation should work.
It appears that the contour calculation algorithm does not find 1-alpha probability within the given limits and therefore takes up the complete space.
For debugging, you could increase the alpha value (e.g. to 0.1 or 0.5) to see if contour calculation works then:
HDC = HighestDensityContour(model, 0.1, limits=[(0, 50), (0, 25), (0, 25)])
or you could calculate an IFORM contour instead as discussed.
yugao0912
commented
2 months ago [https://res.public.onecdn.static.microsoft/assets/mail/file-icon/png/generic_16x16.png]3D_contour.pyhttps://1drv.ms/u/c/93aa9071f6020a74/EXnm0ZL_ii1Mqkt03TEuAKwB4TxAMebmMU3O-SzOMIqYhw
Hi Andreas,
Thank you for your suggestion. I have uploaded the file to the OneDrive. I am not sure whether it could work or not. If it can not still work, would you mind provide your personal email so that I could send you the file. Otherwise, I can only attach the codes.
Best wishes, Yu
发件人: Andreas Haselsteiner @.> 发送时间: 2024年7月18日 22:28 收件人: virocon-organization/virocon @.> 抄送: Yu Gao @.>; Mention @.> 主题: Re: [virocon-organization/virocon] Problems with the results from the 3D_contour.py (Issue #219)
Thanks. The data file is unfortunately not uploaded. Could you give the upload another try? I'd like to reproduce the problem on my machine. The figures of the model fit look good to me such that I would expect that contour calculation should work.
It appears that the contour calculation algorithm does not find 1-alpha probability within the given limits and therefore takes up the complete space.
For debugging, you could increase the alpha value (e.g. to 0.1 or 0.5) to see if contour calculation works then: HDC = HighestDensityContour(model, 0.1, limits=[(0, 50), (0, 25), (0, 25)])
or you could calculate an IFORM contour instead as discussed.
― Reply to this email directly, view it on GitHubhttps://github.com/virocon-organization/virocon/issues/219#issuecomment-2237647446, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A56V74BYOISBUR6L7J3G3B3ZNAXOHAVCNFSM6AAAAABLBESGTSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMZXGY2DONBUGY. You are receiving this because you were mentioned.Message ID: @.***>
yugao0912
commented
1 month ago Thanks. The data file is unfortunately not uploaded. Could you give the upload another try? I'd like to reproduce the problem on my machine. The figures of the model fit look good to me such that I would expect that contour calculation should work.
It appears that the contour calculation algorithm does not find 1-alpha probability within the given limits and therefore takes up the complete space.
For debugging, you could increase the alpha value (e.g. to 0.1 or 0.5) to see if contour calculation works then:
HDC = HighestDensityContour(model, 0.1, limits=[(0, 50), (0, 25), (0, 25)])or you could calculate an IFORM contour instead as discussed.
Hello Andreas,
I am aware that this script may be still not available. I would copy the codes from the scripts here.
And I also have another question about how to plot different contours in one figure. For example, I would like to plot the 1-year contour and 50-year contour based on the same dataset in one figure. But it seems that the original method does not work for this.
These codes are directly pasted from the script. I did not change anything but the path of the file.
Looking forward to your reply.
Best wishes, Yu
ahaselsteiner
commented
1 month ago Hi Yu,
I had a look: The upper limit for Tz was too low, that's why the algorithm used the complete probability space within the given limits.
I increased the limit to 50 seconds. Now a contour with 1 - alpha probability can be found:
This is the line of code I changed:
HDC = HighestDensityContour(model, alpha, limits=[(0, 20), (0, 10), (0, 50)]) # dimensions: wind speed, hs, tz
And here is the complete code for easy reference:
![https://res.public.onecdn.static.microsoft/assets/mail/file-icon/png/generic_16x16.png]3D_contour.py](https://res.public.onecdn.static.microsoft/assets/mail/file-icon/png/generic_16x16.png]3D_contour.py){kind=link}
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from virocon import (
GlobalHierarchicalModel,
LogNormalDistribution,
ExponentiatedWeibullDistribution,
DependenceFunction,
WidthOfIntervalSlicer,
HighestDensityContour,
plot_marginal_quantiles,
plot_dependence_functions,
)
# Load sea state measurements.
#data = pd.read_csv("datasets/coastDat2_oneyear.csv", sep=";", skipinitialspace=True)
data = pd.read_csv("3D.csv", sep=";", skipinitialspace=True)
data.index = pd.to_datetime(data.pop(data.columns[0]), format="%Y-%m-%d-%H")
def _power3(x, a, b, c):
return a + b * x**c
def _exp3(x, a, b, c):
return a + b * np.exp(c * x)
def _alpha3(x, a, b, c, d_of_x):
return (a + b * x**c) / 2.0445 ** (1 / d_of_x(x))
def _logistics4(x, a=1, b=1, c=-1, d=1):
return a + b / (1 + np.exp(c * (x - d)))
bounds = [(0, None), (0, None), (None, None)]
logistics_bounds = [(0, None), (0, None), (None, 0), (0, None)]
power3 = DependenceFunction(_power3, bounds, latex="$a + b * x^c$")
exp3 = DependenceFunction(_exp3, bounds, latex="$a + b * \exp(c * x)$")
logistics4 = DependenceFunction(
_logistics4,
logistics_bounds,
weights=lambda x, y: y,
latex="$a + b / (1 + \exp[c * (x -d)])$",
)
alpha3 = DependenceFunction(
_alpha3,
bounds,
d_of_x=logistics4,
weights=lambda x, y: y,
latex="$(a + b * x^c) / 2.0445^{1 / F()}$",
)
# Wind speed.
dist_description_0 = {
"distribution": ExponentiatedWeibullDistribution(),
"intervals": WidthOfIntervalSlicer(2, min_n_points=50),
}
# Wave height.
dist_description_1 = {
"distribution": ExponentiatedWeibullDistribution(f_delta=5),
"intervals": WidthOfIntervalSlicer(0.5),
"conditional_on": 0,
"parameters": {
"alpha": alpha3,
"beta": logistics4,
},
}
# Zero-up-crossing period.
dist_description_2 = {
"distribution": LogNormalDistribution(),
"conditional_on": 1,
"parameters": {"mu": power3, "sigma": exp3},
}
model = GlobalHierarchicalModel(
[dist_description_0, dist_description_1, dist_description_2]
)
# Define a dictionary that describes the model.
semantics = {
"names": ["Wind speed", "Significant wave height", "Zero-up-crossing period"],
"symbols": ["V", "H_s", "T_z"],
"units": ["m/s", "m", "s"],
}
# Fit the model to the data (estimate the model's parameter values).
model.fit(data)
# Print the estimated parameter values.
print(model)
# Create plots to inspect the model's goodness-of-fit.
fig1, axs = plt.subplots(1, 3)
plot_marginal_quantiles(model, data, semantics, axes=axs)
fig2, axs = plt.subplots(1, 4)
plot_dependence_functions(model, semantics, axes=axs)
# Calculate 3D Contour.
state_duration = 3 # hours
return_period = 50 # years
alpha = state_duration / (return_period * 365.25 * 24)
HDC = HighestDensityContour(model, alpha, limits=[(0, 20), (0, 10), (0, 50)])
# randomly select only 5% of the contour's points to increase performance
rng = np.random.default_rng(42)
n_contour_points = len(HDC.coordinates)
random_points = rng.choice(
n_contour_points, int(0.05 * n_contour_points), replace=False
)
Xs = HDC.coordinates[random_points, 0]
Ys = HDC.coordinates[random_points, 1]
Zs = HDC.coordinates[random_points, 2]
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(Xs, Ys, Zs, c="#004488")
ax.set_xlabel("Wind speed (m/s)")
ax.set_ylabel("Significant wave height (m)")
ax.set_zlabel("Zero-up-crossing period (s)")
plt.show()
For showing two contours in one plot: I'd use slices through the "3D contour" (the surface) instead of a 3D plot.
-Andreas
Hello,
I was trying to run the _3Dcontour.py file. It can be run successfully. But the results seem to be wired. I am unsure whether it is because the maximum of my wind speed is 5 m/s.
Besides, I would like to ask what should I do if I want to change the method to compute the contour. Do I just need to use the IFORMContour to replace the HighestDensityContour?
Looking forward to your reply.
Best wishes, Yu