colosseum.hardness.analysis

33@dataclass()
34class SingleInstanceHyperOptConfig:
35    mdp_class: Type["BaseMDP"]
36    """The class of the MDP with respect to whom we are performing the parameters optimization."""
37    mdp_parameters: Dict[str, Any]
38    """The dictionary containing the parameters for the MDP class with respect to whom we are performing the parameters optimization."""
39    n_seed: int
40    """The number of seed the agent/MDP interaction is repeated."""
41    optimization_horizon: int
42    """The optimization horizon of the agent/MDP interaction."""
43    num_samples: int
44    """The number of samples of the agent parameters."""
45    max_interaction_s: float
46    """The maximum amount of seconds allocated to training the agent."""
47    log_every: int
48    """The length of the interval between logging of the performance indicators."""
49    episodic_near_optimal_agent_class: Type["BaseAgent"] = PSRLEpisodic
50    """The class of the near optimal agent for the episodic setting. By default, it is ``PSRLEpisodic``."""
51    continuous_near_optimal_agent_class: Type["BaseAgent"] = UCRL2Continuous
52    """The class of the near optimal agent for the continuous setting. By default, it is ``UCRL2Continuous``."""

 20@gin.configurable
 21class PSRLEpisodic(BaseAgent):
 22    """
 23    The posterior sampling for reinforcement learning algorithm.
 24
 25    Osband, Ian, Daniel Russo, and Benjamin Van Roy. "(More) efficient reinforcement learning via posterior sampling."
 26    Advances in Neural Information Processing Systems 26 (2013).
 27    """
 28
 29    def step_update(
 30        self,
 31        ts_t: dm_env.TimeStep,
 32        a_t: "ACTION_TYPE",
 33        ts_tp1: dm_env.TimeStep,
 34        time: int,
 35    ):
 36        super(PSRLEpisodic, self).step_update(ts_t, a_t, ts_tp1, time)
 37
 38    @staticmethod
 39    def is_emission_map_accepted(emission_map: "EmissionMap") -> bool:
 40        return emission_map.is_tabular
 41
 42    @staticmethod
 43    def produce_gin_file_from_parameters(parameters: Dict[str, Any], index: int = 0):
 44        return (
 45            "from colosseum.agent.mdp_models import bayesian_models\n"
 46            f"prms_{index}/PSRLEpisodic.reward_prior_model = %bayesian_models.RewardsConjugateModel.N_NIG\n"
 47            f"prms_{index}/PSRLEpisodic.transitions_prior_model = %bayesian_models.TransitionsConjugateModel.M_DIR\n"
 48            f"prms_{index}/PSRLEpisodic.rewards_prior_prms = [{parameters['prior_mean']}, 1, 1, 1]\n"
 49            f"prms_{index}/PSRLEpisodic.transitions_prior_prms = [{parameters['transition_prior']}]"
 50        )
 51
 52    @staticmethod
 53    def is_episodic() -> bool:
 54        return True
 55
 56    @staticmethod
 57    def get_hyperparameters_search_spaces() -> Dict[str, tune.sample.Domain]:
 58        return dict(
 59            prior_mean=tune.uniform(0.001, 2.0), transition_prior=tune.uniform(0.001, 2)
 60        )
 61
 62    @staticmethod
 63    def get_agent_instance_from_parameters(
 64        seed: int,
 65        optimization_horizon: int,
 66        mdp_specs: MDPSpec,
 67        parameters: Dict[str, Any],
 68    ) -> "BaseAgent":
 69        return PSRLEpisodic(
 70            mdp_specs=mdp_specs,
 71            seed=seed,
 72            optimization_horizon=optimization_horizon,
 73            reward_prior_model=RewardsConjugateModel.N_NIG,
 74            transitions_prior_model=TransitionsConjugateModel.M_DIR,
 75            rewards_prior_prms=[parameters["prior_mean"], 1, 1, 1],
 76            transitions_prior_prms=[parameters["transition_prior"]],
 77        )
 78
 79    @property
 80    def current_optimal_stochastic_policy(self) -> np.ndarray:
 81        T_map, R_map = self._mdp_model.get_map_estimate()
 82        Q, _ = episodic_value_iteration(self._time_horizon, T_map, R_map)
 83        return get_policy_from_q_values(Q, True)
 84
 85    def __init__(
 86        self,
 87        seed: int,
 88        mdp_specs: MDPSpec,
 89        optimization_horizon: int,
 90        # MDP model parameters
 91        reward_prior_model: RewardsConjugateModel = None,
 92        transitions_prior_model: TransitionsConjugateModel = None,
 93        rewards_prior_prms=None,
 94        transitions_prior_prms=None,
 95        # Actor parameters
 96        epsilon_greedy: Union[float, Callable] = None,
 97        boltzmann_temperature: Union[float, Callable] = None,
 98    ):
 99        """
100
101        Parameters
102        ----------
103        seed : int
104            The random seed.
105        mdp_specs : MDPSpec
106            The full specification of the MDP.
107        optimization_horizon : int
108            The total number of interactions that the agent is expected to have with the MDP.
109        reward_prior_model : RewardsConjugateModel, optional
110            The reward priors.
111        transitions_prior_model : TransitionsConjugateModel, optional
112            The transitions priors.
113        rewards_prior_prms : Any
114            The reward prior parameters.
115        transitions_prior_prms : Any
116            The transitions prior parameters.
117        epsilon_greedy : Union[float, Callable], optional
118            The probability of selecting an action at random. It can be provided as a float or as a function of the
119            total number of interactions. By default, the probability is set to zero.
120        boltzmann_temperature : Union[float, Callable], optional
121            The parameter that controls the Boltzmann exploration. It can be provided as a float or as a function of
122            the total number of interactions. By default, Boltzmann exploration is disabled.
123        """
124
125        mdp_model = BayesianMDPModel(
126            seed,
127            mdp_specs,
128            reward_prior_model=reward_prior_model,
129            transitions_prior_model=transitions_prior_model,
130            rewards_prior_prms=rewards_prior_prms,
131            transitions_prior_prms=transitions_prior_prms,
132        )
133        actor = QValuesActor(seed, mdp_specs, epsilon_greedy, boltzmann_temperature)
134
135        super(PSRLEpisodic, self).__init__(
136            seed,
137            mdp_specs,
138            mdp_model,
139            actor,
140            optimization_horizon,
141        )
142
143    def episode_end_update(self):
144        Q, _ = episodic_value_iteration(self._time_horizon, *self._mdp_model.sample())
145        self._actor.set_q_values(Q)
146
147    def before_start_interacting(self):
148        self.episode_end_update()

 34@gin.configurable
 35class UCRL2Continuous(BaseAgent):
 36    """
 37    The second version of upper confidence for reinforcement learning algorithm.
 38
 39    Auer, Peter, Thomas Jaksch, and Ronald Ortner. "Near-optimal regret bounds for reinforcement learning." Advances in
 40    neural information processing systems 21 (2008).
 41
 42    Fruit, Ronan, Matteo Pirotta, and Alessandro Lazaric. "Improved analysis of ucrl2 with empirical bernstein inequality."
 43    arXiv preprint arXiv:2007.05456 (2020).
 44    """
 45
 46    @staticmethod
 47    def is_emission_map_accepted(emission_map: "EmissionMap") -> bool:
 48        return emission_map.is_tabular
 49
 50    @staticmethod
 51    def produce_gin_file_from_parameters(parameters: Dict[str, Any], index: int = 0):
 52        string = f"prms_{index}/UCRL2Continuous.bound_type_p='bernstein'\n"
 53        for k, v in parameters.items():
 54            string += f"prms_{index}/UCRL2Continuous.{k} = {v}\n"
 55        return string[:-1]
 56
 57    @staticmethod
 58    def is_episodic() -> bool:
 59        return False
 60
 61    @staticmethod
 62    def get_hyperparameters_search_spaces() -> Dict[str, tune.sample.Domain]:
 63        return {"alpha_p": tune.uniform(0.1, 3), "alpha_r": tune.uniform(0.1, 3)}
 64
 65    @staticmethod
 66    def get_agent_instance_from_parameters(
 67        seed: int,
 68        optimization_horizon: int,
 69        mdp_specs: "MDPSpec",
 70        parameters: Dict[str, Any],
 71    ) -> "BaseAgent":
 72        return UCRL2Continuous(
 73            mdp_specs=mdp_specs,
 74            seed=seed,
 75            optimization_horizon=optimization_horizon,
 76            alpha_p=parameters["alpha_p"],
 77            alpha_r=parameters["alpha_r"],
 78            bound_type_p="bernstein",
 79        )
 80
 81    @property
 82    def current_optimal_stochastic_policy(self) -> np.ndarray:
 83        Q, _ = discounted_value_iteration(self.P, self.estimated_rewards)
 84        return get_policy_from_q_values(Q, True)
 85
 86    def __init__(
 87        self,
 88        seed: int,
 89        mdp_specs: "MDPSpec",
 90        optimization_horizon: int,
 91        # MDP model parameters
 92        alpha_r=1.0,
 93        alpha_p=1.0,
 94        bound_type_p="_chernoff",
 95        bound_type_rew="_chernoff",
 96        # Actor parameters
 97        epsilon_greedy: Union[float, Callable] = None,
 98        boltzmann_temperature: Union[float, Callable] = None,
 99    ):
100        r"""
101        Parameters
102        ----------
103        seed : int
104            The random seed.
105        mdp_specs : MDPSpec
106            The full specification of the MDP.
107        optimization_horizon : int
108            The total number of interactions that the agent is expected to have with the MDP.
109        alpha_r : float
110            The :math:`\alpha` parameter for the rewards. By default, it is set to one.
111        alpha_p : float
112            The :math:`\alpha` parameter for the transitions. By default, it is set to one.
113        bound_type_p : str
114            The upper confidence bound type for the transitions. It can either be '_chernoff' or 'bernstein'. By default,
115            it is set to '_chernoff'.
116        bound_type_rew : str
117            The upper confidence bound type for the rewards. It can either be '_chernoff' or 'bernstein'. By default,
118            it is set to '_chernoff'.
119        epsilon_greedy : Union[float, Callable], optional
120            The probability of selecting an action at random. It can be provided as a float or as a function of the
121            total number of interactions. By default, the probability is set to zero.
122        boltzmann_temperature : Union[float, Callable], optional
123            The parameter that controls the Boltzmann exploration. It can be provided as a float or as a function of
124            the total number of interactions. By default, Boltzmann exploration is disabled.
125        """
126
127        n_states = self._n_states = mdp_specs.observations.num_values
128        n_actions = self._n_actions = mdp_specs.actions.num_values
129        self.reward_range = mdp_specs.rewards_range
130
131        assert bound_type_p in ["_chernoff", "bernstein"]
132        assert bound_type_rew in ["_chernoff", "bernstein"]
133
134        self.alpha_p = alpha_p
135        self.alpha_r = alpha_r
136
137        # initialize matrices
138        self.policy = np.zeros((n_states,), dtype=np.int_)
139        self.policy_indices = np.zeros((n_states,), dtype=np.int_)
140
141        # initialization
142        self.iteration = 0
143        self.episode = 0
144        self.delta = 1.0  # confidence
145        self.bound_type_p = bound_type_p
146        self.bound_type_rew = bound_type_rew
147
148        self.P = np.ones((n_states, n_actions, n_states), np.float32) / n_states
149
150        self.estimated_rewards = (
151            np.ones((n_states, n_actions), np.float32) * mdp_specs.rewards_range[1]
152        )
153        self.variance_proxy_reward = np.zeros((n_states, n_actions), np.float32)
154        self.estimated_holding_times = np.ones((n_states, n_actions), np.float32)
155        self.N = np.zeros((n_states, n_actions, n_states), dtype=np.int32)
156
157        self.current_state = None
158        self.artificial_episode = 0
159        self.episode_reward_data = dict()
160        self.episode_transition_data = dict()
161
162        super(UCRL2Continuous, self).__init__(
163            seed,
164            mdp_specs,
165            None,
166            QValuesActor(seed, mdp_specs, epsilon_greedy, boltzmann_temperature),
167            optimization_horizon,
168        )
169
170    def is_episode_end(
171        self,
172        ts_t: dm_env.TimeStep,
173        a_t: "ACTION_TYPE",
174        ts_tp1: dm_env.TimeStep,
175        time: int,
176    ) -> bool:
177        nu_k = len(self.episode_transition_data[ts_t.observation, a_t])
178        return nu_k >= max(1, self.N[ts_t.observation, a_t].sum() - nu_k)
179
180    def episode_end_update(self):
181        self.episode += 1
182        self.delta = 1 / math.sqrt(self.iteration + 1)
183
184        new_sp = self.solve_optimistic_model()
185        if new_sp is not None:
186            self.span_value = new_sp / self.reward_range[1]
187
188        if len(self.episode_transition_data) > 0:
189            self.model_update()
190            self.episode_reward_data = dict()
191            self.episode_transition_data = dict()
192
193    def before_start_interacting(self):
194        self.episode_end_update()
195
196    def step_update(
197        self, ts_t: dm_env.TimeStep, a_t: "ACTION_TYPE", ts_tp1: dm_env.TimeStep, h: int
198    ):
199        self.N[ts_t.observation, a_t, ts_tp1.observation] += 1
200
201        if (ts_t.observation, a_t) in self.episode_reward_data:
202            self.episode_reward_data[ts_t.observation, a_t].append(ts_tp1.reward)
203            if not ts_tp1.last():
204                self.episode_transition_data[ts_t.observation, a_t].append(
205                    ts_tp1.observation
206                )
207        else:
208            self.episode_reward_data[ts_t.observation, a_t] = [ts_tp1.reward]
209            if not ts_tp1.last():
210                self.episode_transition_data[ts_t.observation, a_t] = [
211                    ts_tp1.observation
212                ]
213
214    def model_update(self):
215        """
216        updates the model given the transitions obtained during the artificial episode.
217        """
218        for (s_tm1, action), r_ts in self.episode_reward_data.items():
219            # updated observations
220            scale_f = self.N[s_tm1, action].sum()
221            for r in r_ts:
222                # update the number of total iterations
223                self.iteration += 1
224
225                # update reward and variance estimate
226                scale_f += 1
227                old_estimated_reward = self.estimated_rewards[s_tm1, action]
228                self.estimated_rewards[s_tm1, action] *= scale_f / (scale_f + 1.0)
229                self.estimated_rewards[s_tm1, action] += r / (scale_f + 1.0)
230                self.variance_proxy_reward[s_tm1, action] += (
231                    r - old_estimated_reward
232                ) * (r - self.estimated_rewards[s_tm1, action])
233
234                # update holding time
235                self.estimated_holding_times[s_tm1, action] *= scale_f / (scale_f + 1.0)
236                self.estimated_holding_times[s_tm1, action] += 1 / (scale_f + 1)
237
238        for (s_tm1, action) in set(self.episode_transition_data.keys()):
239            self.P[s_tm1, action] = self.N[s_tm1, action] / self.N[s_tm1, action].sum()
240
241    def beta_r(self, nb_observations) -> np.ndarray:
242        """
243        calculates the confidence bounds on the reward.
244        Returns
245        -------
246        np.array
247            The vector of confidence bounds on the reward function (|S| x |A|)
248        """
249        S = self._n_states
250        A = self._n_actions
251        if self.bound_type_rew != "bernstein":
252            ci = _chernoff(
253                it=self.iteration,
254                N=nb_observations,
255                range=self.reward_range[1],
256                delta=self.delta,
257                sqrt_C=3.5,
258                log_C=2 * S * A,
259            )
260            return self.alpha_r * ci
261        else:
262            N = np.maximum(1, nb_observations)
263            Nm1 = np.maximum(1, nb_observations - 1)
264            var_r = self.variance_proxy_reward / Nm1
265            log_value = 2.0 * S * A * (self.iteration + 1) / self.delta
266            beta = bernstein(
267                scale_a=14 * var_r / N,
268                log_scale_a=log_value,
269                scale_b=49.0 * self.r_max / (3.0 * Nm1),
270                log_scale_b=log_value,
271                alpha_1=math.sqrt(self.alpha_r),
272                alpha_2=self.alpha_r,
273            )
274            return beta
275
276    def beta_p(self, nb_observations) -> np.ndarray:
277        """
278        calculates the confidence bounds on the transition probabilities.
279        Returns
280        -------
281        np.array
282            The vector of confidence bounds on the reward function (|S| x |A|)
283        """
284        S = self._n_states
285        A = self._n_actions
286        if self.bound_type_p != "bernstein":
287            beta = _chernoff(
288                it=self.iteration,
289                N=nb_observations,
290                range=1.0,
291                delta=self.delta,
292                sqrt_C=14 * S,
293                log_C=2 * A,
294            )
295            return self.alpha_p * beta.reshape([S, A, 1])
296        else:
297            N = np.maximum(1, nb_observations)
298            Nm1 = np.maximum(1, nb_observations - 1)
299            var_p = self.P * (1.0 - self.P)
300            log_value = 2.0 * S * A * (self.iteration + 1) / self.delta
301            beta = bernstein(
302                scale_a=14 * var_p / N[:, :, np.newaxis],
303                log_scale_a=log_value,
304                scale_b=49.0 / (3.0 * Nm1[:, :, np.newaxis]),
305                log_scale_b=log_value,
306                alpha_1=math.sqrt(self.alpha_p),
307                alpha_2=self.alpha_p,
308            )
309            return beta
310
311    def solve_optimistic_model(self) -> Union[None, float]:
312        """
313        solves the optimistic value iteration.
314        Returns
315        -------
316        Union[None, float]
317            The span value of the estimates from the optimistic value iteration or None if no solution has been found.
318        """
319        nb_observations = self.N.sum(-1)
320
321        beta_r = self.beta_r(nb_observations)  # confidence bounds on rewards
322        beta_p = self.beta_p(
323            nb_observations
324        )  # confidence bounds on transition probabilities
325
326        T = self.P
327        estimated_rewards = self.estimated_rewards
328
329        assert np.isclose(T.sum(-1), 1.0).all()
330        try:
331            res = extended_value_iteration(
332                T, estimated_rewards, beta_r, beta_p, self.reward_range[1]
333            )
334        except SystemError:
335            # Debug logs if the optimistic value iteration fails
336            os.makedirs(f"tmp{os.sep}error_ext_vi", exist_ok=True)
337            for i in range(100):
338                if not os.path.isfile(f"tmp{os.sep}error_ext_vi{os.sep}T{i}.npy"):
339                    np.save(f"tmp{os.sep}error_ext_vi{os.sep}T{i}.npy", T)
340                    np.save(
341                        f"tmp{os.sep}error_ext_vi{os.sep}estimated_rewards{i}.npy",
342                        estimated_rewards,
343                    )
344                    np.save(f"tmp{os.sep}error_ext_vi{os.sep}beta_r.npy{i}", beta_r)
345                    np.save(f"tmp{os.sep}error_ext_vi{os.sep}beta_p.npy{i}", beta_p)
346                    break
347            res = None
348
349        if res is not None:
350            span_value_new, self.Q, self.V = res
351            span_value = span_value_new
352            self._actor.set_q_values(self.Q)
353
354            assert span_value >= 0, "The span value cannot be lower than zero"
355            assert np.abs(span_value - span_value_new) < 1e-8
356
357            return span_value
358        return None

55@dataclass()
56class HardnessAnalysisParams:
57    mdp_class: Type["BaseMDP"]
58    """The class of the MDP whose hardness we are studying."""
59    varying_params_name: str
60    """The name of the parameter being varied."""
61    varying_params_values: Iterable
62    """The values of the parameter being varied."""
63    fixed_params: Dict[str, Any]
64    """The dictionary containing the names and values for the parameters being kept fixed."""
65    n_seeds_mdp: int
66    """The number of seeds used when instantiating the MDP."""
67    hardness_measures: Collection[Union[str, Callable[[Type["BaseMDP"]], float]]] = (
68        "diameter",
69        "value_norm",
70    )
71    """An iterable containing either the code name of an available measure of hardness or a function that takes an MDP
72    object as input and returns a value."""
73    near_optimal_agent_hyperopt_config: SingleInstanceHyperOptConfig = None
74    """The parameters optimization configuration for the near optimal agent. By default, it is None, which means
75    that the regret of the near optimal agent with tuned parameters is not computed and used as proxy for a 
76    complete measure of hardness."""
77    varying_params_name_clean: str = None
78    """The name of the parameter being varied in a clean format."""
79    retrieve_from_cache: bool = True
80    """If ture, the ``config.get_hardness_measures_cache_folder()`` is searched for a cached value of the measure.
81    By default, it is True."""
82
83    @property
84    def clean_varying_prm_name(self) -> str:
85        """
86        Returns
87        -------
88        str
89            A nicely formatted name for the varying parameter.
90        """
91        if self.varying_params_name_clean is None:
92            return self.varying_params_name
93        return self.varying_params_name_clean

 96def run_scenario_analysis(
 97    hap: HardnessAnalysisParams,
 98    ax=None,
 99):
100    """
101    runs a hardness analysis scenario.
102
103    Parameters
104    ----------
105    hap : HardnessAnalysisParams
106        The hardness analysis scenario to run.
107    ax : plt.Axes
108        The ax object where the plot will be put. By default, a new axis is created.
109    """
110
111    show = ax is None
112    if ax is None:
113        fig, ax = plt.subplots()
114
115    dfs = get_varying_parameter_dfs(hap, normalize_measures=True)
116
117    for measure_name, df in dfs.items():
118        sns.lineplot(
119            x=hap.clean_varying_prm_name,
120            y=measure_name,
121            data=df,
122            ax=ax,
123            label=measure_name,
124        )
125    plt.ylabel("Hardness measure value")
126    plt.legend()
127
128    if show:
129        plt.show()

132def get_varying_parameter_dfs(
133    hap: HardnessAnalysisParams,
134    normalize_measures: bool = False,
135) -> Dict[str, pd.DataFrame]:
136    """
137    computes the hardness measures for the scenarios and stored them in a `pd.DataFrame`s.
138
139    Parameters
140    ----------
141    hap : HardnessAnalysisParams
142        The hardness analysis scenario to run.
143    normalize_measures : bool
144        If True, the values of the hardness measures are normalized.
145
146    Returns
147    -------
148    Dict[str, pd.DataFrame]
149        A dictionary that associated the name of a hardness measure to its corresponding pd.DataFrame.
150    """
151
152    if hap.near_optimal_agent_hyperopt_config is not None:
153        raise NotImplementedError(
154            "The computation of the regret of the near optimal agent as a proxy for a complete measure of hardness is "
155            "being refactored at the moment."
156        )
157
158    measure_results = dict()
159    if config.get_available_cores() > 1:
160        inputs = [
161            (
162                hap.mdp_class,
163                hap.fixed_params,
164                hap.varying_params_name,
165                varying_value,
166                seed,
167                measure,
168                hap.retrieve_from_cache,
169            )
170            for varying_value in hap.varying_params_values
171            for measure in hap.hardness_measures
172            for seed in range(hap.n_seeds_mdp)
173        ]
174
175        with Pool(processes=config.get_available_cores()) as p:
176            for measure_name, varying_value, seed, measure_value in p.starmap_async(
177                _compute_hardness_measure, inputs
178            ).get():
179                _add_result(
180                    measure_results, measure_name, varying_value, seed, measure_value
181                )
182    else:
183        for seed in range(hap.n_seeds_mdp):
184            for measure in hap.hardness_measures:
185                for varying_value in hap.varying_params_values:
186                    out = compute_hardness_measure_for_varying_prm(
187                        hap.mdp_class,
188                        hap.fixed_params,
189                        hap.varying_params_name,
190                        varying_value,
191                        seed,
192                        measure,
193                        force_single_core=True,
194                        retrieve_from_cache=hap.retrieve_from_cache,
195                        return_n_states=hap.varying_params_name == "size",
196                    )
197                    if hap.varying_params_name == "size":
198                        measure_name, measure_value, n_states = out
199                    else:
200                        measure_name, measure_value = out
201
202                    _add_result(
203                        measure_results,
204                        measure_name,
205                        varying_value
206                        if hap.varying_params_name != "size"
207                        else n_states,
208                        seed,
209                        measure_value,
210                    )
211
212    for measure_name in measure_results:
213        # Create a Pandas DataFrame
214        df = pd.DataFrame.from_dict(measure_results[measure_name])
215
216        # Normalize the values if required
217        if normalize_measures:
218            min_value = df.loc[:, measure_name].min()
219            max_value = df.loc[:, measure_name].max()
220            if max_value > min_value + 1e-4:
221                df.loc[:, measure_name] = (df.loc[:, measure_name] - min_value) / (
222                    max_value - min_value
223                )
224            else:
225                df.loc[:, measure_name] = 0.5
226
227        # Clean the varying parameter name
228        df = df.rename(columns={"Varying value": hap.clean_varying_prm_name}).set_index(
229            hap.clean_varying_prm_name
230        )
231
232        measure_results[measure_name] = df
233
234    return measure_results