Skip to content

Commit

Permalink
update generating initial parameter realizations
Browse files Browse the repository at this point in the history
use the same seed for all realizations. only shuffle before writing the file.
  • Loading branch information
martinvonk committed Sep 29, 2024
1 parent 34e8bc9 commit e56a22c
Showing 1 changed file with 19 additions and 21 deletions.
40 changes: 19 additions & 21 deletions pastas_plugins/pest/solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -406,7 +406,9 @@ def __init__(
self.computername = get_computername()
copy_file(self.exe_agent, self.temp_ws) # copy agent executable

def solve(self, silent: bool = False, **kwargs) -> Tuple[bool, np.ndarray, np.ndarray]:
def solve(
self, silent: bool = False, **kwargs
) -> Tuple[bool, np.ndarray, np.ndarray]:
"""
Solve the optimization problem using the pest_hp solver.
Expand Down Expand Up @@ -571,6 +573,7 @@ def run_ensembles(
The correlation coefficient of the observation noise, by default 0.0.
ies_parameter_ensemble_method : Optional[Literal["norm", "truncnorm", "uniform"]], optional
The method to distribution of the prior for the parameter ensemble, by default None.
If None the parameter distribution is drwan by pestpp-ies itself.
pestpp_options : Optional[Dict], optional
Additional PEST++ options, by default None.
Returns
Expand Down Expand Up @@ -637,8 +640,7 @@ def parameter_distribution(
pmax: float,
par_sigma_range: float,
method: Literal["norm", "truncnorm", "uniform"],
seed: int = pyemu.en.SEED,
) -> np.array:
) -> np.ndarray[float]:
"""Generate a distribution of parameter values based on the specified method.
Parameters
Expand All @@ -657,8 +659,6 @@ def parameter_distribution(
Method to use for generating the distribution. 'norm' generates a
normal distribution, 'truncnorm' generates a truncated normal
distribution, and 'uniform' generates a uniform distribution.
seed : int, optional
Random seed for reproducibility, by default pyemu.en.SEED.
Returns
-------
Expand All @@ -667,9 +667,7 @@ def parameter_distribution(
"""
if method == "norm":
scale = min(initial - pmin, pmax - initial) / (par_sigma_range / 2)
rvs = norm(loc=initial, scale=scale).rvs(
size=ies_num_reals, random_state=seed
)
rvs = np.sort(norm(loc=initial, scale=scale).rvs(size=ies_num_reals))
rvs[rvs < pmin] = pmin
rvs[rvs > pmax] = pmax
elif method == "truncnorm":
Expand All @@ -688,17 +686,15 @@ def parameter_distribution(
right_ies_num_reals = int(
np.ceil((pmax - initial) / (pmax - pmin) * ies_num_reals)
)
rvs_left = tnorm_left.rvs(size=left_ies_num_reals, random_state=seed)
rvs_right = tnorm_right.rvs(size=right_ies_num_reals, random_state=seed)
rvs = np.append(rvs_left, rvs_right)[:ies_num_reals]
rvs_left = tnorm_left.rvs(size=left_ies_num_reals)
rvs_right = tnorm_right.rvs(size=right_ies_num_reals)
rvs = np.sort(np.append(rvs_left, rvs_right)[:ies_num_reals])
rvs[rvs < pmin] = pmin
rvs[rvs > pmax] = pmax
elif method == "uniform":
# rvs = uniform(loc=pmin, scale=pmax).rvs(ies_num_reals, random_state=pyemu.en.SEED)
rvs = np.linspace(
pmin, pmax, ies_num_reals
start=pmin, stop=pmax, num=ies_num_reals
) # linspace ensures pmin and pmax are in the ensembles
np.random.default_rng(seed=seed).shuffle(rvs)
else:
raise ValueError(f"{method=} should be 'norm', 'truncnorm' or 'uniform'.")
return rvs
Expand All @@ -710,7 +706,7 @@ def generate_observation_noise(
standard_deviation: float,
correlation_coefficient: float = 0.0,
seed: int = pyemu.en.SEED,
) -> np.array:
) -> np.ndarray[float]:
"""Generate a matrix of normally distributed and optionally correlated noise
Parameters
Expand Down Expand Up @@ -745,6 +741,7 @@ def write_ensemble_parameter_distribution(
method: Literal["norm", "truncnorm", "uniform"] = "norm",
par_sigma_range: float = 4.0,
ies_add_base: bool = True,
seed: int = pyemu.en.SEED,
) -> None:
"""
Generate and write an ensemble of parameter distributions to a CSV file.
Expand All @@ -762,6 +759,8 @@ def write_ensemble_parameter_distribution(
ies_add_base : bool, optional
If True, add the base parameter values to the ensemble. Default is
True.
seed : int, optional
Random seed for reproducibility, by default pyemu.en.SEED.
Returns
-------
Expand All @@ -771,7 +770,6 @@ def write_ensemble_parameter_distribution(
par_df = pd.DataFrame(
index=pd.Index(range(self.ies_num_reals)), columns=pst.parameter_data.index
)
seed = pyemu.en.SEED
for pname, pdata in pst.parameter_data.iterrows():
rvs = PestIesSolver.parameter_distribution(
ies_num_reals=self.ies_num_reals,
Expand All @@ -780,10 +778,12 @@ def write_ensemble_parameter_distribution(
pmax=pdata.at["parubnd"],
par_sigma_range=par_sigma_range,
method=method,
seed=seed,
)
seed += 1
par_df[pname] = rvs
# shuffle each column with the initial parameters independently
par_df.loc[:, :] = np.random.default_rng(seed=seed).permuted(
par_df.values, axis=0
)
if ies_add_base:
par_df.loc[self.ies_num_reals - 1] = pst.parameter_data.loc[
:, "parval1"
Expand Down Expand Up @@ -1112,9 +1112,7 @@ def __init__(
)

def start(
self,
pestpp_options: Optional[Dict] = None,
silent: bool = False
self, pestpp_options: Optional[Dict] = None, silent: bool = False
) -> None:
"""
Start the PESTPP-SEN analysis.
Expand Down

0 comments on commit e56a22c

Please sign in to comment.