magnetic_field_factory.py 6.26 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# IMAGINE is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18 19 20 21 22

import numpy as np

from keepers import Loggable

Theo Steininger's avatar
Theo Steininger committed
23
from imagine.carrier_mapper import unity_mapper
24

25 26
from nifty import FieldArray, RGSpace

27 28 29 30 31 32
from magnetic_field import MagneticField


class MagneticFieldFactory(Loggable, object):

    def __init__(self, box_dimensions, resolution):
Theo Steininger's avatar
Theo Steininger committed
33
        self.logger.debug("Setting up MagneticFieldFactory.")
34 35
        self.box_dimensions = box_dimensions
        self.resolution = resolution
36 37 38
        self._parameter_defaults = self._initial_parameter_defaults
        self._variable_to_parameter_mappings = \
            self._initial_variable_to_parameter_mappings
39

40 41 42 43 44 45 46 47 48 49 50 51
        distances = np.array(self.box_dimensions) / np.array(self.resolution)
        self._grid_space = RGSpace(shape=self.resolution,
                                   distances=distances)
        self._vector = FieldArray(shape=(3,))
        self._ensemble_cache = {}

    def _get_ensemble(self, ensemble_size):
        if ensemble_size not in self._ensemble_cache:
            self._ensemble_cache[ensemble_size] = \
                                FieldArray(shape=(ensemble_size,))
        return self._ensemble_cache[ensemble_size]

52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
    @property
    def box_dimensions(self):
        return self._box_dimensions

    @box_dimensions.setter
    def box_dimensions(self, box_dimensions):
        dim = tuple(np.array(box_dimensions, dtype=np.float))
        if len(dim) != 3:
            raise ValueError("Input of box_dimensions must have length three.")
        self._box_dimensions = dim

    @property
    def resolution(self):
        return self._resolution

    @resolution.setter
    def resolution(self, resolution):
        resolution = tuple(np.array(resolution, dtype=np.int))
        if len(resolution) != 3:
            raise ValueError("Input for resolution must have length three.")
        self._resolution = resolution

74 75 76
    @property
    def magnetic_field_class(self):
        return MagneticField
77

78 79 80
    @property
    def _initial_parameter_defaults(self):
        return {}
81

82 83 84
    @property
    def _initial_variable_to_parameter_mappings(self):
        return {}
85 86 87

    @staticmethod
    def _interval(mean, sigma, n):
88
        return [mean-n*sigma, mean+n*sigma]
89 90 91

    @staticmethod
    def _positive_interval(mean, sigma, n):
92
        return [max(0, mean-n*sigma), mean+n*sigma]
93 94 95 96 97 98 99 100 101 102 103

    @property
    def parameter_defaults(self):
        return self._parameter_defaults

    @parameter_defaults.setter
    def parameter_defaults(self, new_defaults):
        self._parameter_defaults.update((str(k), np.float(v))
                                        for k, v in new_defaults.items()
                                        if k in self._parameter_defaults)

104 105 106 107 108 109 110 111 112
    @property
    def variable_defaults(self):
        variable_defaults = {}
        for parameter in self.parameter_defaults:
            low, high = self.variable_to_parameter_mappings[parameter]
            default = self.parameter_defaults[parameter]
            variable_defaults[parameter] = (default - low)/(high - low)
        return variable_defaults

113 114 115 116 117 118 119 120 121
    @property
    def variable_to_parameter_mappings(self):
        return self._variable_to_parameter_mappings

    @variable_to_parameter_mappings.setter
    def variable_to_parameter_mappings(self, new_mapping):
        """
        The parameter-mapping must be a dictionary with
        key: parameter-name
122
        value: [min, max]
123 124 125 126
        """
        for k, v in new_mapping.items():
            if k in self._variable_to_parameter_mappings:
                key = str(k)
Theo Steininger's avatar
Theo Steininger committed
127 128
                value = [np.float(v[0]), np.float(v[1])]
                self._variable_to_parameter_mappings.update({key: value})
129 130 131 132 133 134 135 136
                self.logger.debug("Updated variable_to_parameter_mapping %s "
                                  "to %s" % (key, str(value)))

    def _map_variables_to_parameters(self, variables):
        parameter_dict = {}
        for variable_name in variables:
            if variable_name in self.variable_to_parameter_mappings:
                mapping = self.variable_to_parameter_mappings[variable_name]
Theo Steininger's avatar
Theo Steininger committed
137 138
                mapped_variable = unity_mapper(variables[variable_name],
                                               a=mapping[0],
139
                                               b=mapping[1])
Theo Steininger's avatar
Theo Steininger committed
140 141 142 143
#                mapped_variable = carrier_mapper(variables[variable_name],
#                                                 a=mapping[0],
#                                                 m=mapping[1],
#                                                 b=mapping[2])
144 145 146 147 148
            else:
                mapped_variable = np.float(variables[variable_name])
            parameter_dict[variable_name] = mapped_variable
        return parameter_dict

149
    def generate(self, variables={}, ensemble_size=1, random_seed=None):
150 151 152 153
        mapped_variables = self._map_variables_to_parameters(variables)
        work_parameters = self.parameter_defaults.copy()
        work_parameters.update(mapped_variables)

154 155 156
        domain = (self._get_ensemble(ensemble_size),
                  self._grid_space,
                  self._vector)
157

158
        result_magnetic_field = self.magnetic_field_class(
159 160
                                              domain=domain,
                                              parameters=work_parameters,
161 162
                                              distribution_strategy='equal',
                                              random_seed=random_seed)
Theo Steininger's avatar
Theo Steininger committed
163 164
        self.logger.debug("Generated magnetic field with work-parameters %s" %
                          work_parameters)
165
        return result_magnetic_field