Source code for magni.cs.phase_transition._simulation

    Copyright (c) 2014-2017, Magni developers.
    All rights reserved.
    See LICENSE.rst for further information.

Module providing the actual simulation functionality.

Routine listings
run(algorithm, path, label)
     Simulate a reconstruction algorithm.

See Also
magni.cs.phase_transition._config : Configuration options.

The results of the simulation are backed up throughout the simulation. In case
the simulation is interrupted during execution, the simulation will resume from
the last backup point when run again.


from __future__ import division
import os
import random
import sys
import time

import numpy as np

from magni.cs.phase_transition import _backup
from magni.cs.phase_transition import config as _conf
from magni.cs.phase_transition import _data
from magni.utils.multiprocessing import File as _File
from magni.utils.multiprocessing import process as _process
from magni.utils import split_path as _split_path

if sys.version_info[0] == 2:
    iter_range = xrange(2**30)
    iter_range = range(2**32)

[docs]def run(algorithm, path, label, pre_simulation_hook=None): """ Simulate a reconstruction algorithm. The simulation results are stored in a HDF5 database rather than returned by the function. Parameters ---------- algorithm : function A function handle to the reconstruction algorithm. path : str The path of the HDF5 database where the results should be stored. label : str The label assigned to the simulation results pre_simumlation_hook : callable A handle to a callable which should be run *just* before the call to the reconstruction algorithm (the default is None, which implies that no pre hook is run). """ tmp_dir = _split_path(path)[0] + '.tmp' + os.sep tmp_file = tmp_dir + label.replace('/', '#') + '.hdf5' if not os.path.isdir(tmp_dir): os.mkdir(tmp_dir) if not os.path.isfile(tmp_file): _backup.create(tmp_file) done = _backup.get(tmp_file) shape = [len(_conf['delta']), len(_conf['rho'])] random.seed(_conf['seed']) seeds = np.array( random.sample(iter_range, shape[0] * shape[1])).reshape(shape) tasks = [(algorithm, (i, j), seeds[i, j], tmp_file, pre_simulation_hook) for i in range(shape[0]) for j in range(shape[1]) if not done[i, j]] _process(_simulate, args_list=tasks, maxtasks=_conf['maxpoints']) with _File(tmp_file, 'r') as f: stat_time = f.root.time[:] stat_dist = f.root.dist[:] stat_mse = f.root.mse[:] stat_norm = f.root.norm[:] with _File(path, 'a') as f: f.create_array('/' + label, 'time', stat_time, createparents=True) f.create_array('/' + label, 'dist', stat_dist, createparents=True) f.create_array('/' + label, 'mse', stat_mse, createparents=True) f.create_array('/' + label, 'norm', stat_norm, createparents=True) os.remove(tmp_file) if len(os.listdir(tmp_dir)) == 0: os.removedirs(tmp_dir)
[docs]def _simulate(algorithm, ij_tuple, seed, path, pre_simulation_hook=None): """ Run a number of monte carlo simulations in a single delta-rho point. The result of a simulation is the simulation error distance, i.e., the ratio between the energy of the coefficient residual and the energy of the coefficient vector. The time of the simulation is the execution time of the reconstruction attempt. Parameters ---------- algorithm : function A function handle to the reconstruction algorithm. ij_tuple : tuple A tuple (i, j) containing the parameters i, j as listed below. i : int The delta-index of the point in the delta-rho grid. j : int The rho-index of the point in the delta-rho grid. seed : int The seed to use in the random number generator when generating the problem suite instances. path : str The path of the HDF5 backup database. pre_simulation_hook : callable A handle to a callable which should be run *just* before the call to the reconstruction algorithm (the default is None, which implies that no pre hook is run). See Also -------- magni.cs.phase_transition._config: Configuration options. magni.cs.phase_transition._data.generate_matrix : Matrix generation. magni.cs.phase_transition._data.generate_vector : Coefficient vector generation. Notes ----- The `pre_simulation_hook` may be used to setup the simulation to match the specfic simulation parameters, e.g. if an oracle estimator is used in the reconstruction algorithm. The `pre_simulation_hook` takes one argument which is the locals() dict. The following reconstruction statistics are computed: * time: Measured algorithm run time in seconds. * dist: Normalised mean squared error (NMSE) - ( ||alpha_hat - alpha|| / ||alpha||)^2 * mse: Mean squared error (MSE) - 1/n * ||alpha_hat - alpha||^2 * norm: True vector norm - ||alpha|| """ i, j = ij_tuple n = _conf['problem_size'] m = int(np.round(n * _conf['delta'][i])) k = int(np.round(m * _conf['rho'][j])) noise = _conf['noise'] stat_time = np.zeros(_conf['monte_carlo'], dtype=np.float64) stat_dist = stat_time.copy() stat_mse = stat_time.copy() stat_norm = stat_time.copy() if k > 0: np.random.seed(seed) for l in range(_conf['monte_carlo']): A = _data.generate_matrix(m, n) alpha = _data.generate_vector(n, k) if noise is not None: e = _data.generate_noise(m, n, k) y = + e else: y = if pre_simulation_hook is not None: pre_simulation_hook(locals()) start = time.time() alpha_hat = algorithm(y, A, **_conf['algorithm_kwargs']) end = time.time() error_norm = np.linalg.norm(alpha_hat - alpha) alpha_norm = np.linalg.norm(alpha) stat_time[l] = end - start stat_dist[l] = (error_norm / alpha_norm)**2 stat_mse[l] = 1.0/n * error_norm**2 stat_norm[l] = alpha_norm _backup.set(path, (i, j), stat_time, stat_dist, stat_mse, stat_norm)