In [1]:

    

import numpy as np
import openturns as ot
from depimpact import ConservativeEstimate
from dependence import quantile_func
import dask
from depimpact.utils import get_grid_sample
from depimpact.tests.test_functions import func_cum_sum_weight
from depimpact.tests.test_functions import func_spec
from depimpact.dependence_plot import set_style_paper

import matplotlib.pyplot as plt
%matplotlib inline

%load_ext autoreload
%autoreload 2


    

In [2]:

    

K = 200
n = 10000
dim = 4
alpha = 0.1
tau_max = 0.8
q_func = quantile_func(alpha)
template_function = func_spec
families = np.tril(np.ones((dim, dim), dtype=int), k=-1)
corr_dim = int(dim*(dim-1)/2)


    

In [3]:

    

def objective(hyperparams, verbose=False):
    margin_params = hyperparams[:dim*2]
    model_params = hyperparams[dim*2:]
    
    margins = []
    for a, b in zip(margin_params[:-1:2], margin_params[1::2]):
        marginal = ot.Uniform(a, b)
        if verbose:
            print(marginal)
        margins.append(marginal)
        
    func = lambda x: template_function(x, a=model_params)
    quant_estimate = ConservativeEstimate(model_func=func, margins=margins, families=families)
    
    indep_quant = quant_estimate.independence(n, q_func=q_func, keep_input_sample=False)
    grid_result = quant_estimate.gridsearch(K, n, q_func=q_func, grid_type='lhs', keep_input_samples=False)
    
    indep_quant = indep_quant.quantity
    min_result = grid_result.min_result   
    min_quantity = min_result.quantity
    min_kendall = min_result.kendall_tau    
    deviation = indep_quant - min_quantity
    
    constraint = 0.
    for kendall in min_kendall:
        constraint += max(0., abs(kendall) - tau_max)
    
    print('Mean of min kendall:', np.mean(np.abs(min_kendall)))
    print('Deviation:', deviation)
    return deviation**2 * (100*constraint - 1.)


    

In [ ]:

    

from skopt import gp_minimize

space_margin_params = []
for i in range(dim):
    space_margin_params.append((-1., -0.01))
    space_margin_params.append((0.01, 1.))
    
space_model_params = [(-1., 1.)]*6
space = space_margin_params + space_model_params


    

In [ ]:

    

n_calls = 100
n_random_starts = 80

res = gp_minimize(objective, space, n_calls=n_calls, verbose=True, n_random_starts=n_random_starts)


    

    




Iteration No: 1 started. Evaluating function at random point.
('Mean of min kendall:', 0.55275000000000007)
('Deviation:', 0.3054226148523935)
Iteration No: 1 ended. Evaluation done at random point.
Time taken: 24.2852
Function value obtained: -0.0933
Current minimum: -0.0933
Iteration No: 2 started. Evaluating function at random point.
('Mean of min kendall:', 0.73919999999999997)
('Deviation:', 0.821418883110292)
Iteration No: 2 ended. Evaluation done at random point.
Time taken: 23.1057
Function value obtained: 20.9571
Current minimum: -0.0933
Iteration No: 3 started. Evaluating function at random point.
('Mean of min kendall:', 0.65669999999999995)
('Deviation:', 0.13927403788219994)
Iteration No: 3 ended. Evaluation done at random point.
Time taken: 22.7569
Function value obtained: 0.2184
Current minimum: -0.0933
Iteration No: 4 started. Evaluating function at random point.
('Mean of min kendall:', 0.61875000000000002)
('Deviation:', 0.24025739380184863)
Iteration No: 4 ended. Evaluation done at random point.
Time taken: 22.9538
Function value obtained: 0.6676
Current minimum: -0.0933
Iteration No: 5 started. Evaluating function at random point.
('Mean of min kendall:', 0.71445000000000025)
('Deviation:', 0.5271181750039338)
Iteration No: 5 ended. Evaluation done at random point.
Time taken: 22.9858
Function value obtained: 11.8463
Current minimum: -0.0933
Iteration No: 6 started. Evaluating function at random point.
('Mean of min kendall:', 0.79199999999999993)
('Deviation:', 0.3570503714098834)
Iteration No: 6 ended. Evaluation done at random point.
Time taken: 22.6487
Function value obtained: 7.0372
Current minimum: -0.0933
Iteration No: 7 started. Evaluating function at random point.
('Mean of min kendall:', 0.6583500000000001)
('Deviation:', 0.2341021630385529)
Iteration No: 7 ended. Evaluation done at random point.
Time taken: 22.9217
Function value obtained: 0.9426
Current minimum: -0.0933
Iteration No: 8 started. Evaluating function at random point.
('Mean of min kendall:', 0.44055000000000005)
('Deviation:', 1.0666525160340223)
Iteration No: 8 ended. Evaluation done at random point.
Time taken: 22.8660
Function value obtained: 7.5262
Current minimum: -0.0933
Iteration No: 9 started. Evaluating function at random point.
('Mean of min kendall:', 0.60555000000000003)
('Deviation:', 0.5453617149927612)
Iteration No: 9 ended. Evaluation done at random point.
Time taken: 22.4849
Function value obtained: 2.4656
Current minimum: -0.0933
Iteration No: 10 started. Evaluating function at random point.
('Mean of min kendall:', 0.81179999999999997)
('Deviation:', 0.4366438012477776)
Iteration No: 10 ended. Evaluation done at random point.
Time taken: 23.0413
Function value obtained: 5.8637
Current minimum: -0.0933
Iteration No: 11 started. Evaluating function at random point.
('Mean of min kendall:', 0.69299999999999995)
('Deviation:', 0.7453713153551842)
Iteration No: 11 ended. Evaluation done at random point.
Time taken: 22.8814
Function value obtained: 9.9365
Current minimum: -0.0933
Iteration No: 12 started. Evaluating function at random point.
('Mean of min kendall:', 0.42404999999999998)
('Deviation:', 0.09659168474079116)
Iteration No: 12 ended. Evaluation done at random point.
Time taken: 23.2036
Function value obtained: 0.0894
Current minimum: -0.0933
Iteration No: 13 started. Evaluating function at random point.
('Mean of min kendall:', 0.67980000000000007)
('Deviation:', 0.5556500264446773)
Iteration No: 13 ended. Evaluation done at random point.
Time taken: 23.0504
Function value obtained: 8.2729
Current minimum: -0.0933
Iteration No: 14 started. Evaluating function at random point.
('Mean of min kendall:', 0.46200000000000002)
('Deviation:', 0.05083090370487353)
Iteration No: 14 ended. Evaluation done at random point.
Time taken: 23.1904
Function value obtained: -0.0026
Current minimum: -0.0933
Iteration No: 15 started. Evaluating function at random point.
('Mean of min kendall:', 0.56759999999999999)
('Deviation:', 0.18246244880515183)
Iteration No: 15 ended. Evaluation done at random point.
Time taken: 24.0936
Function value obtained: 0.8033
Current minimum: -0.0933
Iteration No: 16 started. Evaluating function at random point.
('Mean of min kendall:', 0.67154999999999998)
('Deviation:', 0.39352650275611945)
Iteration No: 16 ended. Evaluation done at random point.
Time taken: 23.5520
Function value obtained: 3.0291
Current minimum: -0.0933
Iteration No: 17 started. Evaluating function at random point.
('Mean of min kendall:', 0.62534999999999996)
('Deviation:', 0.8290347607648338)
Iteration No: 17 ended. Evaluation done at random point.
Time taken: 23.4424
Function value obtained: 11.1411
Current minimum: -0.0933
Iteration No: 18 started. Evaluating function at random point.
('Mean of min kendall:', 0.64185000000000003)
('Deviation:', 0.514973020359379)
Iteration No: 18 ended. Evaluation done at random point.
Time taken: 23.4900
Function value obtained: 2.1985
Current minimum: -0.0933
Iteration No: 19 started. Evaluating function at random point.
('Mean of min kendall:', 0.6946500000000001)
('Deviation:', 0.8562404195088538)
Iteration No: 19 ended. Evaluation done at random point.
Time taken: 23.3342
Function value obtained: 18.4167
Current minimum: -0.0933
Iteration No: 20 started. Evaluating function at random point.
('Mean of min kendall:', 0.65670000000000006)
('Deviation:', 0.534698049938902)
Iteration No: 20 ended. Evaluation done at random point.
Time taken: 22.6837
Function value obtained: 2.0871
Current minimum: -0.0933
Iteration No: 21 started. Evaluating function at random point.
('Mean of min kendall:', 0.49830000000000024)
('Deviation:', 0.7836149974595696)
Iteration No: 21 ended. Evaluation done at random point.
Time taken: 23.5209
Function value obtained: 10.7490
Current minimum: -0.0933
Iteration No: 22 started. Evaluating function at random point.
('Mean of min kendall:', 0.52800000000000014)
('Deviation:', 0.7470466941859177)
Iteration No: 22 ended. Evaluation done at random point.
Time taken: 23.2047
Function value obtained: 7.9415
Current minimum: -0.0933
Iteration No: 23 started. Evaluating function at random point.

In [ ]:

    

K = 500
margin_params = res.x[:dim*2]
model_params = res.x[dim*2:]
res_margins = []
for a, b in zip(margin_params[:-1:2], margin_params[1::2]):
    res_margins.append(ot.Uniform(a, b))
res_func = lambda x: template_function(x, a=model_params)
res_quant_estimate = ConservativeEstimate(model_func=res_func, margins=res_margins, families=families)


    

In [ ]:

    

grid_result_lhs = res_quant_estimate.gridsearch(K, n, q_func=q_func, grid_type='lhs', keep_input_samples=False)
print('Min quantile: {0}\nMin Kendall: {1}'.format(grid_result_lhs.min_result.quantity, grid_result_lhs.min_result.kendall_tau))


    

In [ ]:

    

grid_result_vertices = res_quant_estimate.gridsearch(K, n, q_func=q_func, grid_type='vertices', keep_input_samples=False)
print('Min quantile: {0}\nMin Kendall: {1}'.format(grid_result_vertices.min_result.quantity, grid_result_vertices.min_result.kendall_tau))


    

In [ ]:

    

indep_result = res_quant_estimate.independence(n, q_func=q_func, keep_input_sample=False)
print('Min quantile: {0} at independence'.format(indep_result.quantity))


    

In [ ]:

    

kendalls_lhs = grid_result_lhs.kendalls
kendalls_vertices = grid_result_vertices.kendalls
dev_kendall_lhs = abs(kendalls_lhs).mean(axis=1)
dev_kendall_vertices = abs(kendalls_vertices).mean(axis=1)
quantities_lhs = grid_result_lhs.quantities
quantities_vertices = grid_result_vertices.quantities
quantity_indep = indep_result.quantity


    

In [ ]:

    

min_kendall_lhs = grid_result_lhs.min_result.kendall_tau
min_kendall_vertices = grid_result_vertices.min_result.kendall_tau
min_dev_kendall_lhs = np.abs(min_kendall_lhs).mean()
min_dev_kendall_vertices = np.abs(min_kendall_vertices).mean()
min_quantity_lhs = grid_result_lhs.min_result.quantity
min_quantity_vertices = grid_result_vertices.min_result.quantity


    

In [ ]:

    

set_style_paper()
fig, ax = plt.subplots(figsize=(7, 4))    
ax.plot(dev_kendall_lhs, quantities_lhs, 'g.', label='EGS K=%d' % (K))
ax.plot(min_dev_kendall_lhs, min_quantity_lhs, 'go', label='Min EGS')
ax.plot(dev_kendall_vertices, quantities_vertices, 'r.', label='BEGS K=%d' % (min(K, 3**dim-1)))
ax.plot(min_dev_kendall_vertices, min_quantity_vertices, 'ro', label='Min BEGS')
ax.plot(0., quantity_indep, 'bo', label='Independence')
ax.legend(loc=0)
ax.set_xlabel('Kendall coefficient deviation')
ax.set_ylabel('Quantile at $\\alpha = %.2f$' % (alpha))
fig.tight_layout()
fig.savefig('./output/optim/non_monotonic_multidim_test_quantile_dim_%d_K_%d.pdf' % (dim, K))
fig.savefig('./output/optim/non_monotonic_multidim_test_quantile_dim_%d_K_%d.png' % (dim, K))


    

In [ ]:

    

fig, axes = plt.subplots(dim-1, dim-1, figsize=(3*dim, 2.5*dim), sharex=True, sharey=True)

k = 0
for i in range(dim-1):
    for j in range(i+1):
        ax = axes[i, j] if dim > 2 else axes
        ax.plot(kendalls_lhs[:, k], quantities_lhs, 'g.', label='EGS K=%d' % (K))
        ax.plot(min_kendall_lhs[k], min_quantity_lhs, 'go', label='Min EGS')
        ax.plot(kendalls_vertices[:, k], quantities_vertices, 'r.', label='BEGS K=%d' % (min(K, 3**dim-1)))
        ax.plot(min_kendall_vertices[k], min_quantity_vertices, 'ro', label='Min BEGS')
        ax.plot(0., quantity_indep, 'bo', label='Independence')
        k += 1
        if i == dim-2:
            ax.set_xlabel('Kendall coefficient')
        if j == 0:
            ax.set_ylabel('Quantile at $\\alpha = %.2f$' % (alpha))
        if i == j+1:
            ax.legend(loc=0)
            
fig.tight_layout()
fig.savefig('./output/optim/matrix_plot_non_monotonic_multidim_test_quantile_dim_%d_K_%d.pdf' % (dim, K))
fig.savefig('./output/optim/matrix_plot_non_monotonic_multidim_test_quantile_dim_%d_K_%d.png' % (dim, K))