Source code for pymls.solver

#! /usr/bin/env python
# -*- coding:utf8 -*-
#
# solver.py
#
# This file is part of pymls, a software distributed under the MIT license.
# For any question, please contact one of the authors cited below.
#
# Copyright (c) 2017
# 	Olivier Dazel <olivier.dazel@univ-lemans.fr>
# 	Mathieu Gaborit <gaborit@kth.se>
# 	Peter Göransson <pege@kth.se>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#

import numpy as np
from numpy.lib.scimath import sqrt

from pymls.analysis import Analysis
from pymls.interface.utils import generic_interface, rigid_interface
from pymls.layers import generic_layer, StochasticLayer
import pymls.backing as backing
from mediapack import Air


[docs]class IncompleteDefinitionError(Exception): """ Exception raised when attempting to solve an incomplete system Either missing layers definition or no backing will produce such an error.""" def __init__(self, msg='The definition is incomplete and no analysis can be performed'): super().__init__(msg)
[docs]class Solver(object): """ Stores a system to solve and parameters for the analysis. Performs analysis and gives back raw unmodified/cleaned results. All post-processing should be done *out* of this class. Parameters ---------- layers : list of Layer/StochasticLayer instances The right most layer appears last in the list. backing : function reference from `pymls.backing` Describe the type of backing condition media : list of Media subclasses, optional Stores all media used in the system for later reference analyses : list of Analysis instances, optional If only one instance is provided with a list, the constructor will wrap it into a list. Attributes ---------- media : list of Media subclasses, optional layers : list of Layer/StochasticLayer instances backing : function reference from `pymls.backing` analyses : list of Analysis instances, optional resultset : list of dict Contains the results for all analysis and metadata Methods ------- solve(frequencies, angles, n_draws, prng_state) : list of dict Starts the solving process w/w stochastic parameters. check_is_complete() : bool Check that all required data has been provided and gathers media. """ def __init__(self, media=None, analyses=None, layers=None, backing=None): self.media = media if media is not None else [] self.layers = layers if layers is not None else [] self.backing = backing if type(analyses) == Analysis: self.analyses = [analyses] elif type(analyses) == list: self.analyses = analyses else: self.analyses = [] self.resultset = []
[docs] def check_is_complete(self): """ Check that all required data has been provided and gathers media. Returns ------- bool True if the described is complete and ready to be solved. Raises ------ IncompleteDefinitionError If the system is incomplete (missing layer or backing) """ # not empty layer list if not self.layers: raise IncompleteDefinitionError("Empty layer list") # if media are missing, grab them from the layers list missing_media = {l.medium for l in self.layers} - set(self.media) self.media += list(missing_media) if self.backing not in [backing.transmission, backing.rigid]: raise IncompleteDefinitionError('No backing provided') return True
[docs] def solve(self, frequencies=None, angles=0, n_draws=1000, prng_state=None): """ Starts the solving process w/w stochastic parameters. The function looks for `StochasticLayer` instances in the `layers` list and flag them. It creates an `Analysis` if all parameters are provided upon call and runs the corresponding solver functions for all registered analysis, gathering results in `resultset`. Parameters ---------- frequencies : list, optional list of frequency where to compute the analysis. If it isn't provided and no `Analysis` has been registered before hand, the function will return nothing. angles : optional Defaults to 0. Can be a list or anything `Analysis` can parse to an iterable. n_draws : int Number of draws for the stochastic analyses. prng_state : tuple Saved state for Numpy's pseudo random number generator (see `numpy.random.get_state`) .. _numpy.random.get_state: https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.get_state.html Returns ------- resultset : dict or list of dict Set of all computed results and relevant metadata provided as a dict for easy serialisation. """ self.check_is_complete() self.resultset = [] self.n_draws = n_draws self.stochastic_layers = list(filter( lambda _: type(_[1]) == StochasticLayer, enumerate(self.layers) )) if self.stochastic_layers and prng_state is None: self.prng_state = np.random.get_state() if frequencies is not None: self.analyses = [Analysis( 'auto', frequencies, angles, len(self.stochastic_layers) > 0 )] self.n_analyses = 0 for a in self.analyses: if a.enable_stochastic: partial_resultset = self.__run_stochastic_analysis(a) self.resultset += partial_resultset else: result = self.__run__analysis(a) self.resultset.append(result) if len(self.resultset) == 1: return self.resultset[0] else: return self.resultset
def __run_stochastic_analysis(self, a): """ Runs a stochastic solver for `Analysis` `a` with stochastic layers Note that it will run the analysis for each stochastic layer separately, re initialising the PRNG state between two runs. Parameters ---------- a : Analysis The one to be run. Returns ------- partial_resultset : list of dict Return a list of all result dicts computed Notes ----- The result dict has the following structure: .. code-block:: python result = { 'name': str_var, 'enable_stochastic': bool_var, 'stochastics': { 'layer': int_id, # index of the layer simulated 'param': str_param, # variated parameter 'values': [], # values of the parameters }, 'f': freqs_vect, 'angle': angles_vect, 'R': [], # reflection coefficient per freq 'T': [], # transmission coefficient per freq } `R` and `T` are matrices (rows corresponding to the freqs/angles axis, columns to the draws) See Also -------- layers.layer.StochasticLayer For stochastic layer definition """ partial_resultset = [] analysis_size = len(a.freqs)*len(a.angles) for l_id, l in self.stochastic_layers: self.n_analyses += 1 self.__reinit_stochastic_solver() result = { 'name': a.name, 'enable_stochastic': a.enable_stochastic, 'stochastics': { 'layer': l_id, 'param': l.stochastic_param, 'values': [], }, 'f': a.freqs, 'angle': a.angles, 'R': [[] for _ in range(analysis_size)], 'T': [[] for _ in range(analysis_size)], } for i_draw in range(self.n_draws): draw = l.new_draw() result['stochastics']['values'].append(draw) for analysis_point, (f, angle) in enumerate(a): (R, T) = self.__solve_one_frequency(f, angle) result['R'][analysis_point].append(R) if T is not None: result['T'][analysis_point].append(T) partial_resultset.append(result) return partial_resultset def __run__analysis(self, a): """ Runs a solver for `Analysis` `a` Parameters ---------- a : Analysis The one to be run. Returns ------- result : dict Return the result of the computation and metadata in a dict format Notes ----- The result dict has the following structure: .. code-block:: python result = { 'name': str_var, 'enable_stochastic': bool_var, 'f': freqs_vect, 'angle': angles_vect, 'R': [], # reflection coefficient per freq 'T': [], # transmission coefficient per freq } """ self.n_analyses += 1 result = { 'name': a.name, 'enable_stochastic': a.enable_stochastic, 'f': a.freqs, 'angle': a.angles, 'R': [], 'T': [], } for f, angle in a: (R, T) = self.__solve_one_frequency(f, angle) result['R'].append(R) if T is not None: result['T'].append(T) return result def __reinit_stochastic_solver(self): """ Utility function to re-initialise the solver between two stochastic analyses. """ np.random.set_state(self.prng_state) for _,l in self.stochastic_layers: l.reinit() def __solve_one_frequency(self, frequency, theta_inc): """ Solve for one frequency `frequency` and one angle `theta_inc` Parameters ---------- frequency : Frequency for the resolution theta_inc : Angle of incidence for the resolution Returns ------- reflx_coefficient : complex Reflection coefficient trans_coefficient : complex or None Transmission coefficient """ omega = frequency*2*np.pi for L in self.layers: L.update_frequency(omega) # compute k_x k_x = omega/Air.c*np.sin(theta_inc*np.pi/180) k_air = omega*sqrt(Air.rho/Air.K) k_z = sqrt(k_air**2-k_x**2) # load the backing vector to initiate recursion Omega_plus = self.backing(omega, k_x) # go backward (from last to first layer) and compute successive # Omega_plus/minus back_prop = np.eye(1) for invertedi_L, L in enumerate(self.layers[::-1]): i_L = len(self.layers)-invertedi_L-1 if invertedi_L == 0: # right-most layer if self.backing == backing.transmission: # check if the last layer is identical to the transmission medium if L.medium.MODEL == 'fluid' and L.medium.c == Air.c and L.medium.rho == Air.rho: Omega_plus *= np.exp(-1j*k_z*L.thickness) continue interface_func = generic_interface(L.medium, Air) else: interface_func = rigid_interface(L.medium) else: interface_func = generic_interface( L.medium, self.layers[i_L+1].medium ) if interface_func is not None: (Omega_minus, tau) = interface_func(Omega_plus) else: Omega_minus = Omega_plus tau = np.eye(int(len(Omega_minus)/2)) layer_func = generic_layer(L.medium) (Omega_plus, xi) = layer_func(Omega_minus, omega, k_x, L.medium, L.thickness) if self.backing == backing.transmission: back_prop = back_prop.dot(tau).dot(xi) # last interface interface_func = generic_interface(Air, self.layers[0].medium) if interface_func is not None: (Omega_minus, tau) = interface_func(Omega_plus) else: Omega_minus = Omega_plus tau = np.eye(int(len(Omega_minus)/2)) if self.backing == backing.transmission: back_prop = back_prop.dot(tau) # Solve for the first layer u_z = 1j*k_z/(Air.rho*omega**2) Omega_0_fluid = np.array([ [-u_z], [-1] ], dtype=np.complex) S_fluid = np.array([ [-u_z], [1] ], dtype=np.complex) temp = np.array([ [Omega_minus[0,0], Omega_0_fluid[0,0]], [Omega_minus[1,0], Omega_0_fluid[1,0]] ]) X = np.linalg.inv(temp).dot(S_fluid) reflx_coefficient = X[1,0] X_0_minus = X[0,0] if self.backing == backing.transmission: trans_coefficient = back_prop*X_0_minus trans_coefficient = trans_coefficient[0,0] else: trans_coefficient = None return (reflx_coefficient, trans_coefficient)
[docs] def compute_fields(self, layer_id, frequency, theta_inc): """ Returns the backpropagation matrix from the first interface to layer num. `layer_id`. Parameters ---------- frequency : float frequency at which backprop is computed theta_inc : angle of incidence layer_id : int id of the layer up to which backpropagate Returns ------- layer_func : callable Function to get the propagation in the layer Raises ------ ValueError : if the id is invalid """ if layer_id >= len(self.layers): raise ValueError("Supplied layer's id is out of range.") omega = frequency*2*np.pi for L in self.layers: L.update_frequency(omega) # compute k_x k_x = omega/Air.c*np.sin(theta_inc*np.pi/180) # load the backing vector to initiate recursion Omega_plus = self.backing(omega, k_x) # goes backward (from last to first layer) and compute successive # Omega_plus/minus back_prop = None for invertedi_L, L in enumerate(self.layers[::-1]): i_L = len(self.layers)-invertedi_L-1 if invertedi_L == 0: # right-most layer if self.backing == backing.transmission: interface_func = generic_interface(L.medium, Air) else: interface_func = rigid_interface(L.medium) else: interface_func = generic_interface(L.medium, self.layers[i_L+1].medium) if interface_func is not None: (Omega_minus, tau) = interface_func(Omega_plus) else: Omega_minus = Omega_plus tau = np.eye(int(len(Omega_minus)/2)) layer_func = generic_layer(L.medium) (Omega_plus, xi) = layer_func(Omega_minus, omega, k_x, L.medium, L.thickness) if i_L < layer_id: # Store the backprop to the desired layer back_prop = tau.dot(xi) if back_prop is None else back_prop.dot(tau).dot(xi) if i_L == layer_id: omega_plus_target = Omega_plus # last interface interface_func = generic_interface(Air, self.layers[0].medium) if interface_func is not None: (Omega_minus, tau) = interface_func(Omega_plus) else: Omega_minus = Omega_plus tau = np.eye(int(len(Omega_minus)/2)) back_prop = tau if layer_id == 0 else back_prop.dot(tau) # Solve for the first layer k_air = omega*sqrt(Air.rho/Air.K) k_z = sqrt(k_air**2-k_x**2) u_z = 1j*k_z/(Air.rho*omega**2) Omega_0_fluid = np.array([ [-u_z], [-1] ], dtype=np.complex) S_fluid = np.array([ [-u_z], [1] ], dtype=np.complex) temp = np.array([ [Omega_minus[0,0], Omega_0_fluid[0,0]], [Omega_minus[1,0], Omega_0_fluid[1,0]] ]) X = np.linalg.inv(temp).dot(S_fluid) S = omega_plus_target.dot(back_prop)*X[0,0] return S.reshape((max(S.shape,)))