popt.misc_tools.optim_tools

  1"""
  2Collection of tools that can be used in optimization schemes. Only put tools here that are so general that they can
  3be used by several optimization schemes. If some method is only applicable to the update scheme you are
  4implementing, leave it in that class.
  5"""
  6import numpy as np
  7from scipy.linalg import block_diag
  8import os
  9from datetime import datetime
 10from scipy.optimize import OptimizeResult
 11
 12
 13def aug_optim_state(state, list_state):
 14    """
 15    Augment the state variables to get one augmented array.
 16
 17    Input:
 18            - state:                Dictionary of state variable for optimization. OBS: 1D arrays!
 19            - list_state:           Fixed list of keys in state dict.
 20
 21    Output:
 22            - aug_state:            Augmented 1D array of state variables.
 23
 24    """
 25    # Start with ensemble of first state variable
 26    aug = state[list_state[0]]
 27
 28    # Loop over the next states (if exists)
 29    for i in range(1, len(list_state)):
 30        aug = np.hstack((aug, state[list_state[i]]))
 31
 32    # Return the augmented array
 33    return aug
 34
 35
 36def update_optim_state(aug_state, state, list_state):
 37    """
 38    Extract the separate state variables from an augmented state array. It is assumed that the augmented state
 39    array is made in aug_optim_state method, hence this is the reverse method.
 40
 41    Input:
 42            - aug_state:                Augmented state array. OBS: 1D array.
 43            - state:                    Dictionary of state variable for optimization.
 44            - list_state:               Fixed list of keys in state dict.
 45
 46    Output:
 47            - state:                    State dictionary updated with aug_state.
 48
 49    """
 50
 51    # Loop over all entries in list_state and extract an array with same number of rows as the key in state
 52    # determines from aug and replace the values in state[key].
 53    # Init. a variable to keep track of which row in 'aug' we start from in each loop
 54    aug_row = 0
 55    for _, key in enumerate(list_state):
 56        # Find no. rows in state[key] to determine how many rows from aug to extract
 57        no_rows = state[key].shape[0]
 58
 59        # Extract the rows from aug and update 'state[key]'
 60        state[key] = aug_state[aug_row:aug_row + no_rows]
 61
 62        # Update tracking variable for row in 'aug'
 63        aug_row += no_rows
 64
 65    # Return
 66    return state
 67
 68
 69def corr2BlockDiagonal(state, corr):
 70    """
 71    Makes the correlation matrix block diagonal. The blocks are the state varible types.
 72
 73    Parameters
 74    ----------
 75    state: dict
 76        Current control state, including state names
 77
 78    corr : array_like
 79        Correlation matrix, of shape (d, d)
 80
 81    Returns
 82    -------
 83    corr_blocks : list
 84        block matrices, one for each variable type
 85
 86    """
 87
 88    statenames = list(state.keys())
 89    corr_blocks = []
 90    for name in statenames:
 91        dim = state[name].size
 92        corr_blocks.append(corr[:dim, :dim])
 93        corr = corr[dim:, dim:]
 94    return corr_blocks
 95
 96
 97def time_correlation(a, state, n_timesteps, dt=1.0):
 98    """
 99    Constructs correlation matrix with time correlation
100    using an autoregressive model.
101
102    .. math::
103        Corr(t_1, t_2) = a^{|t_1 - t_2|}
104    
105    Assumes that each varaible in state is time-order such that
106    `x = [x1, x2,..., xi,..., xn]`, where `i` is the time index, 
107    and `xi` is d-dimensional.
108
109    Parameters
110    -------------------------------------------------------------
111    a : float
112        Correlation coef, in range (0, 1).
113
114    state : dict
115        Control state (represented in a dict).
116    
117    n_timesteps : int
118        Number of time-steps to correlate for each component.
119    
120    dt : float or int
121        Duration between each time-step. Default is 1.
122
123    Returns
124    -------------------------------------------------------------
125    out : numpy.ndarray
126        Correlation matrix with time correlation    
127    """
128    dim_states = [int(state[name].size/n_timesteps) for name in list(state.keys())]
129    blocks     = []
130
131    # Construct correlation matrix
132    # m: variable type index  
133    # i: first time index
134    # j: second time index
135    # k: first dim index
136    # l: second dim index
137    for m in dim_states:
138        corr_single_block = np.zeros((m*n_timesteps, m*n_timesteps))
139        for i in range(n_timesteps):
140            for j in range(n_timesteps):
141                for k in range(m):
142                    for l in range(m):
143                       corr_single_block[i*m + k, j*m + l] = (k==l)*a**abs(dt*(i-j))
144        blocks.append(corr_single_block)
145
146    return block_diag(*blocks)
147
148
149def cov2corr(cov):
150    """
151    Transfroms a covaraince matrix to a correlation matrix
152
153    Parameters
154    -------------
155    cov : array_like
156        The covaraince matrix, of shape (d,d).
157
158    Returns
159    -------------
160    out : numpy.ndarray
161        The correlation matrix, of shape (d,d)
162    """
163    std  = np.sqrt(np.diag(cov))
164    corr = np.divide(cov, np.outer(std, std))
165    return corr
166
167
168def corr2cov(corr, std):
169    """
170    Transfroms a correlation matrix to a covaraince matrix
171
172    Parameters
173    ----------
174    corr : array_like
175        The correlation matrix, of shape (d,d).
176
177    std : array_like
178        Array of the standard deviations, of shape (d, ).
179
180    Returns
181    -------
182    out : numpy.ndarray
183        The covaraince matrix, of shape (d,d)
184    """
185    cov = np.multiply(corr, np.outer(std, std))
186    return cov
187
188
189def get_sym_pos_semidef(a):
190    """
191    Force matrix to positive semidefinite
192
193    Parameters
194    ----------
195    a : array_like
196        The input matrix, of shape (d,d)
197
198    Returns
199    -------
200    a : numpy.ndarray
201        The positive semidefinite matrix, of shape (d,d)
202    """
203
204    rtol = 1e-05
205    if not isinstance(a, int):
206        S, U = np.linalg.eigh(a)
207        if not np.all(S > 0):
208            S = np.clip(S, rtol, None)
209            a = (U * S) @ U.T
210    else:
211        a = np.maximum(a, rtol)
212    return a
213
214
215def clip_state(x, bounds):
216    """
217    Clip a state vector according to the bounds
218
219    Parameters
220    ----------
221    x : array_like
222        The input state
223
224    bounds : array_like
225        (min, max) pairs for each element in x. None is used to specify no bound.
226
227    Returns
228    -------
229    x : numpy.ndarray
230        The state after truncation
231    """
232
233    any_not_none = any(any(item) for item in bounds)
234    if any_not_none:
235        lb = np.array(bounds)[:, 0]
236        lb = np.where(lb is None, -np.inf, lb)
237        ub = np.array(bounds)[:, 1]
238        ub = np.where(ub is None, -np.inf, ub)
239        x = np.clip(x, lb, ub)
240    return x
241
242
243def get_optimize_result(obj):
244    """
245    Collect optimize results based on requested
246
247    Parameters
248    ----------
249    obj : popt.loop.optimize.Optimize
250        An instance of an optimization class
251
252    Returns
253    -------
254    save_dict : scipy.optimize.OptimizeResult
255        The requested optimization results
256    """
257
258    # Initialize dictionary of variables to save
259    save_dict = OptimizeResult({'success': True, 'x': obj.mean_state, 'fun': np.mean(obj.obj_func_values),
260                                'nit':  obj.iteration, 'nfev': obj.nfev, 'njev': obj.njev})
261    if hasattr(obj, 'epf') and obj.epf:
262        save_dict['epf_iteration'] = obj.epf_iteration
263
264    if 'savedata' in obj.options:
265
266        # Make sure "SAVEDATA" gives a list
267        if isinstance( obj.options['savedata'], list):
268            savedata = obj.options['savedata']
269        else:
270            savedata = [ obj.options['savedata']]
271      
272        # Loop over variables to store in save list
273        for save_typ in savedata:
274            if 'mean_state' in save_typ:
275                continue  # mean_state is alwaysed saved as 'x'
276            if save_typ in locals():
277                save_dict[save_typ] = eval('{}'.format(save_typ))
278            elif hasattr( obj, save_typ):
279                save_dict[save_typ] = eval(' obj.{}'.format(save_typ))
280            else:
281                print(f'Cannot save {save_typ}!\n\n')
282
283    if 'save_folder' in obj.options:
284        save_dict['save_folder'] = obj.options['save_folder']
285
286    return save_dict
287
288
289def save_optimize_results(intermediate_result):
290    """
291    Save optimize results
292
293    Parameters
294    ----------
295    intermediate_result : scipy.optimize.OptimizeResult
296        An instance of an OptimizeResult class
297    """
298    # Cast to OptimizeResult if a ndarray is passed as argument
299    if type(intermediate_result) is np.ndarray:
300        intermediate_result = OptimizeResult({'x': intermediate_result})
301
302    # Make folder (if it does not exist)
303    if 'save_folder' in intermediate_result:
304        save_folder = intermediate_result['save_folder']
305        if not os.path.exists(save_folder):
306            os.makedirs(save_folder)
307    else:
308        save_folder = './'
309
310    if 'nit' in intermediate_result:
311        suffix = str(intermediate_result['nit'])
312    else:
313        now = datetime.now()  # current date and time
314        suffix = now.strftime("%m_%d_%Y_%H_%M_%S")
315
316    # Save the variables
317    if 'epf_iteration' in intermediate_result:
318        np.savez(save_folder + '/optimize_result_{0}_{1}'.format(suffix, str(intermediate_result['epf_iteration'])), **intermediate_result)
319    else:
320        np.savez(save_folder + '/optimize_result_{0}'.format(suffix), **intermediate_result)