Set points in equadratures polynomial?

Questions
1

Is there a way to set the points in an equadratures polynomial object?

I did a bunch of model runs (I’m using the COAWST model) using the points output from equadratures:

order_parameters = 2
wave_height = eq.Parameter(distribution=‘uniform’, lower=0.001, upper=0.5, order=order_parameters)
wave_steepness = eq.Parameter(distribution=‘uniform’, lower=0.025, upper=0.055, order=order_parameters)
veg_density_powerscale = eq.Parameter(distribution=‘uniform’, lower=-0.3, upper=0.3, order=order_parameters)
veg_ht_powerscale = eq.Parameter(distribution=‘uniform’, lower=-0.3, upper=0.3, order=order_parameters)
veg_stem_powerscale = eq.Parameter(distribution=‘uniform’, lower=-0.3, upper=0.3, order=order_parameters)
start_location = eq.Parameter(distribution=‘uniform’, lower=-0.3, upper=0.3, order=order_parameters)
sediment_ws_power = eq.Parameter(distribution=‘uniform’, lower=-1.30103, upper=0.69897, order=order_parameters)
sediment_M_power = eq.Parameter(distribution=‘uniform’, lower=-5, upper=-4, order=order_parameters)
sediment_taucrit = eq.Parameter(distribution=‘uniform’, lower=0.01, upper=0.2, order=order_parameters)
sediment_supply_powerscale = eq.Parameter(distribution=‘uniform’, lower=-1, upper=1, order=order_parameters)
parameters = [wave_height,wave_steepness,veg_density_powerscale,veg_ht_powerscale,veg_stem_powerscale,start_location,sediment_ws_power,sediment_M_power,sediment_taucrit,sediment_supply_powerscale]

mybasis = eq.Basis(‘total-order’)
mypoly = eq.Poly(parameters, mybasis, method=‘least-squares’,
sampling_args={‘mesh’:‘tensor-grid’, ‘subsampling-algorithm’:‘qr’, ‘sampling-ratio’:1.0})

pts = mypoly.get_points()

(Note: mypoly.basis.cardinality is 66)

And now I want to find the Sobol’ indices of the model outputs. I load in some outputs of the model runs, eg. dissip_veg_5tc_marshflooded_allveg, and then:

mypoly.set_model(dissip_veg_5tc_marshflooded_allveg)
mean, var = mypoly.get_mean_and_variance()
mypoly.get_summary(‘SobolOut_DissipVeg.txt’)

But, I’ve been re-running the code listed above before I start the Sobol index analysis (starting with order_parameters = 2), and I’m getting parameter values to use in a different order than the order I got and used for the model runs. This results in computing large Sobol indices for parameters that don’t matter!

Is there a way to either load parameter values in to a poly, or to make sure that i’m generating a poly in the same order as before?

Thank you!

2

Hi @rmallen86, do you by chance have the old points and weights? If so, this might be the easiest fix:

mypoly._quadrature_points = old_pts
mypoly._quadrature_weights = old_wts

In case you don’t have the weights, then we could use the ordering in the points. I’m still not entirely sure what would cause a change in the ordering. I shall have a further ponder.

3

Thanks @psesh ,

I don’t have the old weights. Could you spell out a little more what it means to “use the ordering in the points”?

Thank you!
Rachel

4

Hi @rmallen86, as there has essentially been some permutation of the points (new vs old), I think the easiest fix would be to figure out what this permutation is, and then go from there. A crude way to do this would be the following. I hope this works :slight_smile: .

import numpy as np

# old_points = quadrature_points_from_before
# new_points = quadrature_points_now

indices = []
for row_new in new_points:
    counter = 0
    for row_old in old_points:
        if np.array_equiv(row_new, row_old):
            indices.append(counter)
        counter += 1

mypoly._quadrature_points = old_points[indices,:]
mypoly._quadrature_weights = mypoly._quadrature_weights[indices]
5

Thanks @psesh. I tried this, but it looks like only a handful of the new points were in the set of old points, so just matching up the sets doesn’t work.

I was looking through the sets of points, and I noticed that the weight of that set depends on the number of parameter values that match with the “base case”. That is, when all 10 parameter values match the base case, the weight is 8.78161992e-02, when 7 match the base case the weight is 2.14395017e-02, when 6 match it is 1.33996886e-02, and when 5 match it is 8.37480537e-03. Do you think I can assign the weights based on the number of parameter values that match the base case?

Thank you!
Rachel

6

Hi @rmallen86, may I trouble you to drop your old points and function evaluations, and new points as CSV files please? Hmm weights are typically function of the parameter values and the number of points, so they are specific to the points chosen. Can you also confirm whether you are using the same numpy / scipy version as before? Cheers!