text:
# Exercice 13) Deep Deterministic Policy Gradients and Proximal Policy Optimization

In this exercise we will investigate two state-of-the-art algorithms: deep deterministic policy gradient (DDPG) and proximal policy optimization (PPO).

We will examine their performance on [Goddard's rocket problem](https://github.com/osannolik/gym-goddard).
This environment comes prepackaged in this notebook's folder, so it can be just imported.

```
First formulated by R. H. Goddard around 1910, this is a classical problem within dynamic optimization and optimal control. The task is simply to find the optimal thrust profile for a vertically ascending rocket in order for it to reach the maximum possible altitude, given that its mass decreases as the fuel is spent and that it is subject to varying drag and gravity.
```

The gym's observation space is the rocket's vertical position, velocity and mass, which are normalized for you to be within $[0, 1]$.

The rocket engine is assumed to be throttled such that the thrust can be continuously controlled between 0 to some maximum limit, which translates to an action space $\mathcal U \in [0, 1]$.

During the episode the reward is set to zero.
If the episode is terminating, the maximum hight reached by the rocket during the experiment is used as reward.
If the rocket has not depleted its tank after $300$ steps a penalty of $-1$ is applied and the episode terminates.

(Photo by SpaceX on Unsplash)

In order for the full notebook to run through, you will also need the stable-baslines3 package.

```
pip install stable-baselines3[extra]
```

```python
from rocket_env import GoddardEnv
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import torch
import torch.nn as nn
from torch.optim import Adam
from stable_baselines3 import PPO, DDPG, A2C
from stable_baselines3.common.callbacks import BaseCallback
from stable_baselines3.common.monitor import Monitor
from stable_baselines3.common.noise import NormalActionNoise
from copy import deepcopy

OPTIMAL_CONTROL = 0.0122079818367078
RANDOM_AVG = 0.00995
RANDOM_BEST = 0.01148
```

## 1) DDPG

Deep Deterministic Policy Gradient (DDPG) was first introduced [2015](https://arxiv.org/abs/1509.02971).
It uses four neural networks to interact on the continuous state and action space.
The first two we already know from the last exercise introduced Actor-Critic method: a critic to estimate the q-value function and an actor which deterministicly - completely trained - should give us the best action for a specific state with respect to the q-value function.
Additionally, DDPG provides two target networks which are (in the beginning) copies of the actor and critic. These target networks are updated time-delayed in a low-pass filter manner.
This enhances stability during the learning process because the networks used to calculate the target are not directly updated on the fly.
For more information see for example [here](https://spinningup.openai.com/en/latest/algorithms/ddpg.html).

During the learning process the critic is updated based on minimizing (!) the following loss function:

\begin{equation}
\mathcal L(w, {\mathcal D}) = \underset{(x,u,r,x',d) \sim {\mathcal D}}{{\mathrm E}}\left[
    \Bigg(  \left(r + \gamma \hat{q}_{}(x',\pi(x', \theta^-), w^-) \right) - \hat{q}(x,u, w) \Bigg)^2
    \right].
\end{equation}

Here, the target is calculated using the target networks. $w$ and $w^-$ are the parameters of the critic and critic-target networks $\hat{q}$, respectively and $\theta^-$ define the parameters of the actor-target network.

The policy/actor network is updated based on the idea to maximize (!) the expected return

\begin{align}
    \max_{\theta} \underset{x \sim {\mathcal D}}{{\mathrm E}}\left[\hat{q}(x, \pi_{}(x, \theta), w) \right],
\end{align}

where $\theta$ are the parameters of the policy network $\pi$.

The updates are performed off-policy by sampling from an experience replay buffer ${\mathcal D}$.
Since we are dealing with an deterministic actor, exploration during training is achieved by adding noise to the sampled actions.

### Task: Implement the DDPG with PyTorch
Write a DDPG algorithm using the Algo. 13.1 in lecture slides 13!
A lot of the functions and classes are already given because we switch our framework here from Tensorflow to [PyTorch](https://pytorch.org/).
There are a lot of tutorials on the PyTorch website to get familiar with how you can perform the update steps which will be most of your work in the follwoing cells.

Execute the following cell to make use of the defined multi layer perceptron and the plot function!

```python
def mlp(sizes, activation, output_activation=nn.Identity):
    """
    Defines a multi layer perceptron using pytorch layers and activation funtions
    """
    layers = []
    for j in range(len(sizes)-1):
        act = activation if j < len(sizes)-2 else output_activation
        layers += [nn.Linear(sizes[j], sizes[j+1]), act()]
    return nn.Sequential(*layers)

def plot_reward_trends(logs):
    repeats = len(logs['envs'])
    ncols = 2
    nrows = int(np.ceil(repeats / ncols))
    fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(ncols*8, nrows*4), sharex=True, sharey=True)
    for ax, _env, _model in zip(axes.flatten(), logs['envs'], logs['models']):
        ax.plot(_env.get_episode_rewards(), label='rewards', color='green')

        ax.set_ylabel('Reward')
        ax.set_xlabel('Episode')
        ax.set_ylim(0, 0.015)
        ax.axhline(OPTIMAL_CONTROL, ls='--', color='red', label='optimal')
        ax.axhline(RANDOM_BEST, ls='--', color='orange', label='best rnd actions')
        ax.axhline(RANDOM_AVG, ls='--', color='yellow', label='average rnd actions')
        ax.legend(loc='lower left')
        ax2 = ax.twinx()
        ax2.plot(_env.get_episode_lengths(), label='ep lengths')
        ax2.legend()
        ax2.set_ylabel('Steps')
    fig.tight_layout()

def test_agent(logs):
    """Deterministic test of the agent in the given env"""
    num_test_episodes = 400
    for i, (env, agent) in enumerate(zip(logs['envs'], logs['models'])):
        episode_reward = 0
        episode_len = 0
        state = env.reset()
        for j in tqdm(range(num_test_episodes)):
            episode_len += 1
            action = agent.decide(state, deterministic=True)
            next_state, reward, done, _ = env.step(action)
            episode_reward += reward
            env.render()
            state = next_state
            if done:
                state = env.reset()
                break

        print(f'Reward during test ({i}): {episode_reward}')
        env.close()
    print('Optimal control: 0.0122')

```

Execute the following cells to make use of the predefined actor & critic and the replay buffer!

```python
class Actor(nn.Module):

    def __init__(self, obs_dim, act_dim, hidden_sizes, activation=nn.ReLU, act_limit=1):
        super().__init__()
        pi_sizes = [obs_dim] + list(hidden_sizes) + [act_dim]
        self.pi = mlp(pi_sizes, activation, nn.Sigmoid)
        self.act_limit = act_limit

    def forward(self, state):
        return self.act_limit * self.pi(state)

    def act(self, state):
        # HINT: toch.no_grad -> result is not needed for backpropagation
        with torch.no_grad():
            return self.act_limit * self.pi(state).numpy()

class Critic(nn.Module):

    def __init__(self, obs_dim, act_dim, hidden_sizes, activation=nn.ReLU):
        super().__init__()
        self.q = mlp([obs_dim + act_dim] + list(hidden_sizes) + [1], activation)

    def forward(self, state, action):
        # torch.cat concatenates action to state

        q = self.q(torch.cat([state, action], dim=-1))
        return torch.squeeze(q, -1)  # To ensure q has right shape.

class ReplayBuffer:

        def __init__(self, obs_dim, action_dim, buffer_size):

            self.state_buf = np.zeros((buffer_size, obs_dim), dtype=np.float32)
            self.next_state_buf = np.zeros((buffer_size, obs_dim), dtype=np.float32)
            self.action_buf = np.zeros((buffer_size, action_dim), dtype=np.float32)
            self.reward_buf = np.zeros(buffer_size, dtype=np.float32)
            self.done_buf = np.zeros(buffer_size, dtype=np.float32)
            self.ptr, self.size, self.max_size = 0, 0, buffer_size

        def push(self, state, action, reward, next_state, done):

            self.state_buf[self.ptr] = state
            self.next_state_buf[self.ptr] = next_state
            self.action_buf[self.ptr] = action
            self.reward_buf[self.ptr] = reward
            self.done_buf[self.ptr] = done
            self.ptr = (self.ptr + 1) % self.max_size
            self.size = min(self.size + 1, self.max_size)

        def fetch(self, batch_size=32):

            idxs = np.random.randint(0, self.size, size=batch_size)

            batch = dict(state=self.state_buf[idxs],
                         next_state=self.next_state_buf[idxs],
                         action=self.action_buf[idxs],
                         reward=self.reward_buf[idxs],
                         done=self.done_buf[idxs])

            return {k: torch.as_tensor(v, dtype=torch.float32) for k, v in batch.items()}
```

Here, fill in the following code template to write a DDPG agent!

```python
class DDPG_agent:
    """Reference:
    https://github.com/openai/spinningup/blob/master/spinup/algos/pytorch/ddpg/ddpg.py
    """

    def __init__(self, env, actor_hidden_size, actor_number_layers, critic_hidden_size, critic_number_layers,
                 buffer_size, actor_lr, critic_lr, gamma, batch_size, learning_starts, tau):#, action_noise):
        self.actor_lr = actor_lr
        self.critic_lr = critic_lr
        self.batch_size = batch_size
        self.learning_starts = learning_starts
        self.gamma = gamma
        self.tau = tau

        obs_dim = env.observation_space.shape[0]
        act_dim = env.action_space.shape[0]

        # define A&C and the replaybuffer
        self.actor = None
        self.critic = None
        self.replay_buffer = None

        ### BEGIN SOLUTION
        self.actor = Actor(obs_dim, act_dim, [actor_hidden_size] * actor_number_layers, nn.ReLU)
        self.critic = Critic(obs_dim, act_dim, [critic_hidden_size] * critic_number_layers, nn.ReLU)
        self.replay_buffer = ReplayBuffer(obs_dim, act_dim, buffer_size)
        ### END SOLUTION

        # defines target networks
        self.actor_target = deepcopy(self.actor)
        self.critic_target = deepcopy(self.critic)

        # Uses Adam optimizer (see to ex12 for more explanation)
        self.actor_optimizer = Adam(self.actor.parameters(), lr=actor_lr)
        self.critic_optimizer = Adam(self.critic.parameters(), lr=critic_lr)

        # define action noise space for clipping
        self.action_space_high = env.action_space.high
        self.action_space_low = env.action_space.low

    def q_loss(self, data):
        ### BEGIN SOLUTION

        state, action, rerwad, next_state, done = data['state'], data['action'], data['reward'], data['next_state'], data['done']

        # Gradient descent for Q-function
        self.critic_optimizer.zero_grad()
        q = self.critic(state, action)
        with torch.no_grad():
            q_pi_targ = self.critic_target(next_state, self.actor_target(next_state))
            backup = rerwad + self.gamma * (1 - done) * q_pi_targ

        # MSE loss against Bellman backup
        return ((q - backup) ** 2).mean()

        ### END SOLUTION

    def pi_loss(self, data):
        """
        Calulate the loss for a gradient ascent(!) step
        """
        ### BEGIN SOLUTION

        state = data['state']

        self.actor_optimizer.zero_grad()
        return -self.critic(state, self.actor(state)).mean()

        ### END SOLUTION

    def deliberate(self, number_updates_per_step=1):
        """
        Fetches number_updates_per_step-times from replay_buffer, computes targets and calculates losses to
        update the q-function using gradient descent and the policy function using gradient ascent
        In the end updates the target networks
        """

        for _ in range(number_updates_per_step):
            ### BEGIN SOLUTION

            data = self.replay_buffer.fetch(self.batch_size)

            # Gradient descent for Q-function
            loss_q = self.q_loss(data)
            loss_q.backward()
            self.critic_optimizer.step()

            # Freeze Q-network so you don't waste computational effort
            # computing gradients for it during the policy learning step.
            for p in self.critic.parameters():
                p.requires_grad = False

            # Gradient ascent step for pi.
            loss_pi = self.pi_loss(data)
            loss_pi.backward()
            self.actor_optimizer.step()

            # Unfreeze Q-network so you can optimize it at next DDPG step.
            for p in self.critic.parameters():
                p.requires_grad = True

            # Finally, update target networks by polyak averaging.
            with torch.no_grad():
                # update Critc-target
                for p, p_targ in zip(self.critic.parameters(), self.critic_target.parameters()):
                    # uses in-place operations "mul_", "add_" to update target
                    # params, as opposed to "mul" and "add", which would make new tensors.
                    p_targ.data.mul_(1 - self.tau)
                    p_targ.data.add_((self.tau) * p.data)

                # update Actor-target
                for p, p_targ in zip(self.actor.parameters(), self.actor_target.parameters()):
                    # uses in-place operations "mul_", "add_" to update target
                    # params, as opposed to "mul" and "add", which would make new tensors.
                    p_targ.data.mul_(1 - self.tau)
                    p_targ.data.add_((self.tau) * p.data)

            ### END SOLUTION

    def decide(self, state, deterministic=False, noise_scale=0.1):
        """
        Returns action as nd-array depending on the state, adds scaled noise
        and clips the action depending on the action space
        """

        ### BEGIN SOLUTION
        action = self.actor.act(torch.as_tensor(state, dtype=torch.float32))

        if deterministic:
            return np.clip(action, self.action_space_low, self.action_space_high)
        else:
            return np.clip(action + noise_scale * np.random.randn(1), self.action_space_low, self.action_space_high)
        ### END SOLUTION

```

Run the following cell to train your agent. Note that, in contrast to classic RL applications, in standard libraries you have often the possibility to update the agent after full episodes instead of stepwise.
This could be helpful in our application due to the sparse reward signal of the rocket environment instead of after each environment step.
To overcome this issue we introduce the variable `learning_starts` and run here first 500 steps in the environment to be sure that we have at least one full episode in the replay buffer to learn from.

To compare episode-wise vs step-wise updating, both is executed 4 times in the following cell.
The learning could take some minutes. If you want to reduce it you can do only 1 run for each update for example (reduce `repeats`) or run less steps at all (reduce `total_timesteps`).

```python
train_logs_ddpg = {'episode': {'envs': [], 'models': []}, 'step_wise': {'envs': [], 'models': []}}
repeats = 4

total_timesteps = 20000

for train_freq, logs_ddpg in train_logs_ddpg.items():

    for rep in range(repeats):

        env = Monitor(GoddardEnv())

        myDDPG_agent = DDPG_agent(env=env, actor_hidden_size=8, actor_number_layers=1,
                               critic_hidden_size=8, critic_number_layers=1, buffer_size=int(1e6),
                               actor_lr=1e-4, critic_lr=1e-4, gamma=0.999, batch_size=256, learning_starts=500, tau=0.005)

        number_updates_per_step = 0
        episode_reward = 0
        episode_len = 0
        rewards = []
        episode_len_vec = []

        state = env.reset()

        for j in tqdm(range(total_timesteps)):
            episode_len += 1
            number_updates_per_step += 1
            action = myDDPG_agent.decide(state)
            next_state, reward, done, _ = env.step(action)
            episode_reward += reward

            myDDPG_agent.replay_buffer.push(state, action, reward, next_state, done)
            #env.render()
            if train_freq == 'step_wise' and j > myDDPG_agent.learning_starts:
                # Train every step using myDDPG_agent.batch_size batch
                myDDPG_agent.deliberate()

            state = next_state

            if done:
                if train_freq == 'episode' and j > myDDPG_agent.learning_starts:
                # Train after each episode number_updates_per_step times using myDDPG_agent.batch_size batch
                    myDDPG_agent.deliberate(number_updates_per_step)
                state = env.reset()
                rewards.append(episode_reward)
                episode_len_vec.append(episode_len)
                episode_reward = 0
                episode_len = 0
                number_updates_per_step = 0

        logs_ddpg['envs'] += [env]
        logs_ddpg['models'] += [myDDPG_agent]
        env.close()
```

    /home/webbah/anaconda3/envs/OMG/lib/python3.8/site-packages/gym/logger.py:30: UserWarning: WARN: Box bound precision lowered by casting to float32
      warnings.warn(colorize('%s: %s'%('WARN', msg % args), 'yellow'))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    /home/webbah/anaconda3/envs/OMG/lib/python3.8/site-packages/torch/autograd/__init__.py:130: UserWarning: CUDA initialization: Found no NVIDIA driver on your system. Please check that you have an NVIDIA GPU and installed a driver from http://www.nvidia.com/Download/index.aspx (Triggered internally at  /pytorch/c10/cuda/CUDAFunctions.cpp:100.)
      Variable._execution_engine.run_backward(

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))

```python
plot_reward_trends(train_logs_ddpg['episode'])
```

```python
plot_reward_trends(train_logs_ddpg['step_wise'])
```

Use the following cell to test your agent on the env using deterministic actions

```python
test_agent(train_logs_ddpg['episode'])
```

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (0): 0.00810419263634432

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (1): 0.011636630039107532

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (2): 0.01150983366270042

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (3): 0.011475247557406032
    Optimal control: 0.0122

```python
test_agent(train_logs_ddpg['step_wise'])
```

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (0): 0.011718857018671036

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (1): 0.011878249167981858

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (2): 0.011886240975892859

    HBox(children=(FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (3): 0.011885732739862931
    Optimal control: 0.0122

### Demo: StableBaselines3 usage

Alternatively, you can use readily available Python packages such as [stable-baselines3](https://github.com/DLR-RM/stable-baselines3) for employment of state-of-the-art algorithms.
In what follows below, the DDPG algorithm as utilized by stable-baselines3 is showcased.

```python
total_timesteps = int(20e3)

class Pbar(BaseCallback):
    def __init__(self):
        super().__init__()
        self.pbar = tqdm(desc='Training', total=total_timesteps)

    def _on_rollout_end(self) -> None:
        """This event is triggered before updating the policy."""
        # self.pbar.update(self.)

    def _on_step(self) -> bool:
        self.pbar.update()
        return True

policy_kwargs = dict(activation_fn=torch.nn.ReLU, net_arch=[dict(pi=[8], vf=[8])])

repeats = 4
logs = {'envs': [], 'models': []}
for _ in range(repeats):

    env = Monitor(GoddardEnv())

    #n_actions = env.action_space.shape[-1]
    #noise_var = 4
    #noise_theta = 25  # stiffness of OU
    #action_noise = OrnsteinUhlenbeckActionNoise(mean=np.zeros(n_actions), theta=noise_theta * np.ones(n_actions),
    #                                        sigma=noise_var * np.ones(n_actions))
    n_actions = env.action_space.shape[-1]
    action_noise = NormalActionNoise(mean=np.zeros(n_actions), sigma=0.1 * np.ones(n_actions))

    policy_kwargs = dict(activation_fn=torch.nn.ReLU, net_arch=dict(pi=[8] * 1, qf=[8] * 1))
    model = DDPG("MlpPolicy", env, policy_kwargs=policy_kwargs, batch_size=256, device='cpu',
             learning_rate=1e-4, gamma=0.999, learning_starts=500)
    model.learn(total_timesteps=total_timesteps, callback=Pbar())
    logs['envs'] += [env]
    logs['models'] += [model]
```

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=20000.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=20000.0), HTML(value='')))

```python
plot_reward_trends(logs)
```

Use the following cell to execute the last trained agent

```python
n_episodes = 1
max_steps_per_episode = 500
tst_logs = {'rewards': []}
for ep in range(n_episodes):
    state = env.reset()
    done = np.zeros(1, dtype=np.bool)
    cum_rew = 0
    k = 0
    while not np.all(done) and k < max_steps_per_episode:
        action, _ = model.predict(state, deterministic=True)
        state, reward, done, info = env.step(action)
        env.render()
        cum_rew += reward
        k += 1

    tst_logs['rewards'].append(cum_rew)
rewards = np.array(tst_logs['rewards']).ravel()
env.close()
```

    /home/wilhelmk/tools/anaconda3/envs/tf-cpu/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
    Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations

## 2) PPO

The [original paper from 2017](https://arxiv.org/abs/1707.06347) for the PPO came up with an idea to combine A3C (having multiple workers) and TRPO (using a trust region to improve the actor).
The PPO algorithm achieves this by hard clipping gradients in order to ensure that new policies won't be too far away from old ones.

In contrast to DDPG, PPO is an on-policy algorithm. In order to still apply mini-batch training, there is a so-called roll-out-buffer that is filled up with the current policy, on whose base a gradient ascent update would be done.

There are two variants on the PPO algorithm, from which we will implement the clip variant.

The actor (policy) update is computed according to

\begin{align}
\theta_{k+1} &= \arg \max_{\theta} \underset{\langle x,u \rangle \sim \pi_{\theta_k}}{{\mathbb E}}\left[\mathcal L(X,U,\theta_k, \theta)\right], \\
\mathcal L(x,u,\theta_k,\theta) &= \min\left(
\frac{\pi_{\theta}(u|x)}{\pi_{\theta_k}(u|x)}  A^{\pi_{\theta_k}}(x,u), \;\;
\text{clip}\left(\frac{\pi_{\theta}(u|x)}{\pi_{\theta_k}(u|x)}, 1 - \epsilon, 1+\epsilon \right) A^{\pi_{\theta_k}}(x, u)
\right),
\end{align}
where the expectation operator denotes the empirical average across the roll-out-buffer, $\theta$ the possible next weights of the actor, $\epsilon$ denotes the threshold within the next update is allowed, and $\mathcal L$ the long-term return.

On the opposite side, the critic's weights $\omega_k$ are updated through the plain stochastic mini-batch gradient descent on the mean squared error $\mathcal C$ between immediately seen rewards $r_{k+1}$ and corresponding estimated value $v_{w_k}(x_k)$ across the roll-out-buffer:
\begin{align}
w_{k+1} &\leftarrow \omega_k - \alpha \nabla\mathcal C(X,R_{k+1},\omega_k), \\
\mathcal C(x,r_{k+1}, \omega_k) &= \left(v_{\omega_k} (x_k)- r_{k+1}\right)^2
\end{align}

### Task: Implement the PPO clip variant with PyTorch
Fill in the below code to run the PPO-clip variant on the rocket environment.

For simplicity, do not contemplate a vectorized environment for synchronous training (like in A2C or A3C).
Moreover, calculate the generalized advantage estimate with $\lambda = 1$.

```python
class StochActor(nn.Module):
    """This stochastic actor learns a Gaussian distribution. While the mean value
    is learnt by a full-fledged MLP, the standard deviation is denoted by a
    single trainable weight. With Pytorch's distribution package, probabilities
    given a certain distribution and an action can be calculated."""
    def __init__(self, state_dim, action_dim):
        super().__init__()
        # The standard deviation is just one trainable weight
        self.log_std = torch.nn.Parameter(
                        torch.as_tensor(
                            -.5 * np.ones(action_dim, dtype=np.float32)))
        # the mean value is estimated by a full MLP
        self.mu_net = mlp(list(state_dim) + [8] + list(action_dim), nn.Tanh, nn.Sigmoid)

    def _distribution(self, state):
        return torch.distributions.normal.Normal(
                self.mu_net(state), torch.exp(self.log_std))

    def forward(self, state, action=None):
        pi = self._distribution(state)
        # if action is None, logp_a will be, too
        if action is None:
            logp_a = None
        else:
            logp_a = pi.log_prob(action).sum(axis=-1)
        return pi, logp_a

class PPOAgent:
    """Reference:
    https://github.com/openai/spinningup/blob/master/spinup/algos/pytorch/ppo/ppo.py
    """

    class RolloutBuffer:
        def __init__(self, size, action_dim, state_dim):
            self.action_buf = np.zeros((size, action_dim[0]), dtype=np.float32)
            self.state_buf = np.zeros((size, state_dim[0]), dtype=np.float32)
            self.rew_buf = np.zeros(size+1, dtype=np.float32)
            self.val_buf = np.zeros(size+1, dtype=np.float32)
            self.logp_buf = np.zeros(size, dtype=np.float32)
            self.i = 0
            self.max_size = size

        def push(self, state, action, reward, value, logp):
            """Append a sample to the buffer"""
            ### BEGIN SOLUTION
            self.state_buf[self.i] = state
            self.action_buf[self.i] = action
            self.rew_buf[self.i] = reward
            self.val_buf[self.i] = value
            self.logp_buf[self.i] = logp
            self.i += 1
            ### END SOLUTION

        def fetch(self, last_val=None, gamma=0.999):
            """Get all data from the rollout buffer.
            Returns a dictionary for state, action, rewards-to-go, advantages, and logp values.
            Entries might be of different length across multiple calls to this function as episodes
            potentially have different lengths. Rewards-to-go and advantages need to be computed from
            rewards, values, and gamma, as outlined in https://arxiv.org/abs/1506.02438 (or as done
            in the spinning up implementation). Normalize the advantage batch by standard scaling."""
            ### BEGIN SOLUTION
            if last_val is not None:
                self.rew_buf[self.i] = last_val
                self.val_buf[self.i] = last_val
            rewards = self.rew_buf[:self.i+1]
            values = self.val_buf[:self.i+1]

            # Advantage calculaton
            w = np.arange(rewards.size - 1)
            w = w - w[:, np.newaxis]
            w = np.triu(gamma ** w.clip(min=0)).T
            advantages = rewards[:-1] + gamma * values[1:] - values[:-1]
            advantages = (advantages.reshape(-1, 1) * w).sum(axis=0)

            # reward
            rewards_to_go = ((gamma*rewards[:-1]).reshape(-1, 1) * w).sum(axis=0)

            # normalize advantages
            advantages = (advantages - advantages.mean()) / np.std(advantages)

            data =  {'state': self.state_buf[:self.i],
                     'action': self.action_buf[:self.i],
                     'rewards_to_go': rewards_to_go.copy(),  # make stride positive
                     'advantages': advantages,
                     'logp': self.logp_buf[:self.i]}
            self.i = 0

            return {k: torch.as_tensor(v, dtype=torch.float32) for k, v in data.items()}
            ### END SOLUTION

    def __init__(self, buffer_size, env):
        self.buf = self.RolloutBuffer(buffer_size, action_dim=env.action_space.shape,
                                 state_dim=env.observation_space.shape)
        # use the below attributes for the actor-critic agent
        self.actor = None
        self.critic = None
        self.n_epochs = 8
        self.clip_ratio = 0.1
        ### BEGIN SOLUTION

        self.actor = StochActor(action_dim=env.action_space.shape, state_dim=env.observation_space.shape)
        self.critic = mlp(list(env.observation_space.shape) + [8, 1], nn.Tanh, nn.Tanh)
        self.pi_opt = Adam(self.actor.parameters(), lr=1e-3)
        self.v_opt = Adam(self.critic.parameters(), lr=1e-3)

        ### END SOLUTION

    def v_loss(self, data):
        ### BEGIN SOLUTION
        state, rewards_to_go = data['state'], data['rewards_to_go']
        return ((self.critic(state).squeeze(-1) - rewards_to_go)**2).mean()
        ### END SOLUTION

    def pi_loss(self, data):
        ### BEGIN SOLUTION
        state, action, advantage, logp_old = [data[k] for k
                                              in ('state', 'action', 'advantages', 'logp')]

        pi, logp = self.actor(state, action)
        ratio = torch.exp(logp - logp_old)
        clip_adv = torch.clamp(ratio, 1-self.clip_ratio, 1+self.clip_ratio) * advantage
        loss_pi = -(torch.min(ratio * advantage, clip_adv)).mean()

        return loss_pi
        ### END SOLUTION

    def deliberate(self, last_value=None):
        """Fetch the rollout buffer for all samples.
        Train the actor and critic for n_epochs by computing the pi and v
        loss iteratively.
        The rewards and values buffer needs a zero appended for terminal states,
        or the critic's estimate appended in case of intermediate states
        that result from timeout. This is denoted by last_value."""
        ### BEGIN SOLUTION
        batch = self.buf.fetch(last_value)
        loss_pi_old = self.pi_loss(batch).item()
        loss_v_old = self.v_loss(batch).item()

        # update policy n_epoch times
        for _ in range(self.n_epochs):
            self.pi_opt.zero_grad()
            loss_pi = self.pi_loss(batch)
            loss_pi.backward()
            self.pi_opt.step()

            self.v_opt.zero_grad()
            loss_v = self.v_loss(batch)
            loss_v.backward()
            self.v_opt.step()

        ### END SOLUTION

    def decide(self, state, deterministic=False):
        """Receive a state as torch tensor and return a tuple of numpy ndarrays
        in the form (action, value, log-probability-of-action-under-pi)"""
        ### BEGIN SOLUTION
        state = torch.as_tensor(state, dtype=torch.float32)
        with torch.no_grad():
            pi = self.actor._distribution(state)
            if deterministic:
                return pi.mean.numpy()
            else:
                a = pi.sample()
            logp_a = pi.log_prob(a).sum(axis=-1)
            v = self.critic(state)
        return a.numpy(), v.numpy(), logp_a.numpy()
        ### END SOLUTION

```

```python
# main training loop

n_episodes = 250
max_steps_per_episode = 500
train_logs_ppo_custom = {'envs': [], 'models': []}

def train_ppo(repeats):
    for rep in range(repeats):
        env = Monitor(GoddardEnv())
        agent = PPOAgent(max_steps_per_episode, env)

        for ep in tqdm(range(n_episodes), desc='Training'):
            state = env.reset()
            done = False
            k = 0
            while not done and k < max_steps_per_episode:
                action, value, logp_a = agent.decide(state)
                next_state, reward, done, info = env.step(action)

                # track experience
                agent.buf.push(state, action, reward, value, logp_a)

                state = next_state
                k += 1

            # if episode finished with done signal, the long term value should be 0
            #  for this terminal state
            if done:
                last_value = None
            else:
                _, last_value, _ = self.decide(state)
            agent.deliberate(last_value)

        train_logs_ppo_custom['envs'] += [env]
        train_logs_ppo_custom['models'] += [agent]
        env.close()
train_ppo(4)
```

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

```python
plot_reward_trends(train_logs_ppo_custom)
```

```python
test_agent(train_logs_ppo_custom)
```

    HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (0): 0.01202398780966707

    HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (1): 0.011901507527816824

    HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (2): 0.01186353153602493

    HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=400.0), HTML(value='')))

    Reward during test (3): 0.011754158952783422
    Optimal control: 0.0122

```python
# more samples
train_ppo(46)
```

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

    HBox(children=(HTML(value='Training'), FloatProgress(value=0.0, max=250.0), HTML(value='')))

```python
# plot average and standard deviation of 50 runs
repeats = len(train_logs_ppo_custom['envs'])
rewards_mat = np.asarray([_env.get_episode_rewards()[:250] for _env in train_logs_ppo_custom['envs']])
rewards_average = np.median(rewards_mat, axis=0)
rewards_std = rewards_mat.std(axis=0)

plt.figure()
ax = plt.gca()
ax.plot(rewards_average, color='green', label='median reward')
ax.plot(rewards_mat.max(axis=0), color='cyan', label='best seen reward')
ax.set_ylabel('Reward')
ax.set_xlabel('Episode')
ax.set_ylim(0, 0.015)
ax.axhline(OPTIMAL_CONTROL, ls='--', color='red', label='optimal')
ax.axhline(RANDOM_BEST, ls='--', color='orange', label='best rnd actions')
ax.axhline(RANDOM_AVG, ls='--', color='yellow', label='average rnd actions')
ax.legend(loc='lower left')
ax.set_title('Summary of 50 experiments')
plt.tight_layout()
```

Very strong performance is sometimes achieved even within 10 training episodes. Here, the 'optimal' performance is based on an analytical solution using a continuous-time model (planning) which is applied to an approximate discrete-time rocket simulation. Hence, there is a systematic deviation between the model of the analytical solution and the actual simulation behavior (i.e., it is actual non-optimal w.r.t. the simplified simulation).

The low median reward curve shows that several repetitions of an algorithm are inevitable for a robust assessment.

### Demo: StableBaselines3 usage

Alternatively you can use readily available Python packages such as [stable-baselines3](https://github.com/DLR-RM/stable-baselines3) for employment of state-of-the-art algorithms.
In what follows below, the PPO algorithm as utilized by stable-baselines3 is showcased.

```python
class Pbar(BaseCallback):
    """This is a callback that helps monitoring the training progress"""
    def __init__(self, total, n_rollout_steps):
        super().__init__()
        self.n_steps = n_rollout_steps
        self.pbar = tqdm(desc='Training', total=total)

    def _on_rollout_end(self) -> None:
        """This event is triggered before updating the policy."""
        self.pbar.update(self.n_steps)

    def _on_step(self) -> bool:
        return True

policy_kwargs = dict(activation_fn=torch.nn.Sigmoid, net_arch=[dict(pi=[8], vf=[8])])

n_steps = 256  # an episode takes around ~250 steps
total_timesteps = n_steps*250  # steps_per_episode * episodes

repeats = 4  # repeat experiment x times to assess scatter through random init
train_logs_ppo_sb3 = {'envs': [], 'models': []}
for _ in range(repeats):

    env = Monitor(GoddardEnv())
    model = PPO('MlpPolicy', env, n_steps=n_steps, n_epochs=8,  verbose=0, device='cpu', learning_rate=1e-3, batch_size=n_steps, policy_kwargs=policy_kwargs,
               clip_range=0.1, clip_range_vf=0.1, gamma=0.999)
    model = model.learn(total_timesteps=total_timesteps, callback=Pbar(total_timesteps, n_steps))
    train_logs_ppo_sb3['envs'] += [env]
    train_logs_ppo_sb3['models'] += [model]
```

```python
plot_reward_trends(train_logs_ppo_sb3)
```

Obviously, the chosen architecture and hyper parameters do well often, but can also fail abruptly.
What might be the reason for these failures?


text:
# Financial Asset Recommender: Profitability Predication

## Introduction

This notebook is aimed to serve as an introduction to the creation of a recommender that operates on an open dataset of stock market prices in order to recommend profitable assets. It covers the data download and processing steps, the calculation of the features fed into the prediction model, and finally a simple model that predicts the profitability of assets to rank them.

Before we get started, we need to confgure where we are to store the dataset and models produced. If you are running this notebook locally, then any folder on your machine should be fine. If you are working within a container, you may need to change the directory to a mounted writable file system.

```python
# Local Mode
storageDIR = "dataset" # creates a dataset directory in the same folder as the notebook

# Container Mode
#storageDIR = "/tmp/dataset" # creates a dataset directory in the /tmp/ directory for the container
```

## Dataset

Different types of financial asset recommendation system use different types of data to prodice their recommendations. This approach is known as Profitability Prediction, where assets that are predicted to gain significant value over the following year are recommended. This type of approach uses past pricing data, i.e. the price for different assets over time to identify pricing trends and hence future profitable assets. Hence, as input we need a dataset that for a range of assets contains their price history over time.

For illustration, in this notebook we will use an open dataset, although you may wish to swap this for data from your market of choice. The dataset used here is the Huge Stock Market Dataset compiled by Boris Marjanovic, and is publically available here.

The dataset comprises historical price and volume data for all US-based stocks and ETFs trading on the NYSE, NASDAQ, and NYSE markets, and runs up to the last quarter of 2017. For each financial asset (stock or ETF), it contains a series of price entries describing the market price of that asset on different days. Each entry is comprised of:
 - Date: The date of the pricing data
 - Open: Opening price for that day
 - High: The maximum price for that day
 - Low: The minimum price for that day
 - Close: The closing price for that day
 - Volume: The amount of the asset that is traded
 - OpenInt: The total number of outstanding contracts held by market participants

## Downloading the Dataset

This dataset can be downloaded through the Kaggle API. Users should make an account on the Kaggle website and download an API token in order to access this dataset on their local machine. This is available from the Account section of the Kaggle user profile.

To use the Kaggle API, we first install the package through pip. We then export the username and the API key from the aforementioned token, and use the download command to fetch and unzip the dataset. For our current experiment we will use only the set of Stocks.

(The below commands can also be entered directly into your terminal, without the ! prefix.)

```python
import os
os.environ["KAGGLE_USERNAME"] = "yourKaggleUsername"
os.environ["KAGGLE_KEY"] = "yourKaggleAPIKey"

!pip install kaggle
from kaggle.api.kaggle_api_extended import KaggleApi
api = KaggleApi()
api.authenticate()
api.dataset_download_files('borismarjanovic/price-volume-data-for-all-us-stocks-etfs', path=storageDIR)
print("Download Complete")
```

    Requirement already satisfied: kaggle in c:\programdata\anaconda3\lib\site-packages (1.5.12)
    Requirement already satisfied: certifi in c:\programdata\anaconda3\lib\site-packages (from kaggle) (2019.11.28)
    Requirement already satisfied: urllib3 in c:\programdata\anaconda3\lib\site-packages (from kaggle) (1.25.8)
    Requirement already satisfied: requests in c:\programdata\anaconda3\lib\site-packages (from kaggle) (2.22.0)
    Requirement already satisfied: tqdm in c:\programdata\anaconda3\lib\site-packages (from kaggle) (4.42.1)
    Requirement already satisfied: python-slugify in c:\programdata\anaconda3\lib\site-packages (from kaggle) (5.0.2)
    Requirement already satisfied: python-dateutil in c:\programdata\anaconda3\lib\site-packages (from kaggle) (2.8.1)
    Requirement already satisfied: six>=1.10 in c:\programdata\anaconda3\lib\site-packages (from kaggle) (1.14.0)
    Requirement already satisfied: idna<2.9,>=2.5 in c:\programdata\anaconda3\lib\site-packages (from requests->kaggle) (2.8)
    Requirement already satisfied: chardet<3.1.0,>=3.0.2 in c:\programdata\anaconda3\lib\site-packages (from requests->kaggle) (3.0.4)
    Requirement already satisfied: text-unidecode>=1.3 in c:\programdata\anaconda3\lib\site-packages (from python-slugify->kaggle) (1.3)
    Download Complete

## Extracting and Loading the Dataset

After downloading the dataset we need to extract the files and then convert it into a Pandas DataFrame, which is somewhat like a large data table that makes the raw data easier to analyse.

```python
import zipfile
import pandas as pd
import numpy as np
import glob, os, random, math

# Unzip the Dataset
with zipfile.ZipFile(storageDIR+"/price-volume-data-for-all-us-stocks-etfs.zip", 'r') as zip_ref:
    zip_ref.extractall(storageDIR)

# Replace this with your dataset path
path = storageDIR+'/Stocks/'
all_files = glob.glob(os.path.join(path, "*.us.txt"))
dfs = []

# Iterating through files and only using non-empty files
for f in all_files:
    if os.path.getsize(f) > 0:
        df = pd.read_csv(f)
        df['Stock'] = f.split('/')[-1].split('.')[0]
        dfs.append(df)
full_kaggle_df = pd.concat(dfs)

print("Dataset Extraction and Loading as Dataframe Complete")
```

    Dataset Extraction and Loading as Dataframe Complete

## Filtering the Dataset

Pandas allows us to perform manipulations on the pricing data so that we can extract only what we need for training the model. For the purpose of our illustration, we will only use pricing data from 2016 and 2017, where we consider 2016 as the 'past' and hence we can use that period for learning what assets are likely to be profitable, and 2017 as the future (which we will use to test how well our model performs later). In effect, for our experiment here we can consider an the date of asset recommendation to be the 1st of January 2017.

Lets first filter the dataset to only hold data from the dates we care about:

```python
# We only require the date, stock name, and closing price for that day
full_kaggle_df = full_kaggle_df[['Date', 'Close', 'Stock']]
full_kaggle_df['Date'] = pd.to_datetime(full_kaggle_df['Date'])
full_kaggle_df['year'] = full_kaggle_df['Date'].dt.year
full_kaggle_df['month'] = full_kaggle_df['Date'].dt.month
full_kaggle_df['day'] = full_kaggle_df['Date'].dt.day

# Selecting only that data from either 2016 or 2017
full_kaggle_df = full_kaggle_df[(full_kaggle_df['year'] >= 2016)].drop(columns=['year', 'month', 'day'])
full_kaggle_df['Date'] = full_kaggle_df['Date'].dt.strftime('%Y-%m-%d')

print("Dataset Filtering Complete")
```

## Weekly Price Change Calculation

Now that we have our data, we need to convert it into a form that makes it easier to work with. In this case, we can convert pricing data over time into weekly price changes for each asset. This makes the data more closely match what we want as output, while also enabling data normalization over the daily price points (which can help avoid outliers and in cases where particular days are missing in the data). In particular, we find the closing price for each stock as per the Monday of each week. We calculate the weekly return of each stock after aggregating close prices weekly. This is defined as the percentage change in the weekly prices and is represented as follows:

\begin{equation}
WeeklyChange = \frac{Current \ closing \ price \ - \ Preceding \ closing \ price}{Preceding \ closing \ price}
\end{equation}

Note that this step can take some time...

```python
stocks = full_kaggle_df['Stock'].unique().tolist()
pricedfs = []
for s in stocks:
    df = full_kaggle_df[full_kaggle_df['Stock'] == s]
    df['Date'] = pd.to_datetime(df['Date']) - pd.to_timedelta(7, unit='d')

    #Grouping close prices weekly
    df = df.groupby(['Stock', pd.Grouper(key='Date', freq='W-MON')])['Close'].nth(-1).reset_index().sort_values('Date')

    #Calculating percentage change in asset prices as returns
    df['returns'] = df['Close'].pct_change()
    df['returns'] = df['returns'].replace([np.inf, -np.inf], np.nan)
    df['returns'] = df['returns'].fillna(df['returns'].mean())

    df = df.iloc[1:]
    pricedfs.append(df)

print("Dataset Price Changes Extracted")
```

    C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:5: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      """

Below we show a sample DataFrame with weekly returns calculated.

```python
pricedfs[0].head()
```

      Stock
      Date
      Close
      returns

      1
      Stocks\a
      2016-01-04
      37.278
      -0.067560

      2
      Stocks\a
      2016-01-11
      36.463
      -0.021863

      3
      Stocks\a
      2016-01-18
      36.952
      0.013411

      4
      Stocks\a
      2016-01-25
      37.031
      0.002138

      5
      Stocks\a
      2016-02-01
      34.193
      -0.076638

We then filter the dataframes with 2016 pricing data for training purposes and those with 9 months of 2017 pricing data for testing, and take a random sample of 1000 assets from the filtered subset to speed up the experiment.

```python
# Filtering those dataframes with 2016 pricing data for training
newdfs_training = []

for d in pricedfs:
    d['Date'] = pd.to_datetime(d['Date'])
    d['year'] = d['Date'].dt.year
    d['month'] = d['Date'].dt.month
    d['day'] = d['Date'].dt.day

    d = d[(d['year'] == 2016)].drop(columns=['year', 'month', 'day'])
    d['Date'] = d['Date'].dt.strftime('%Y-%m-%d')
    if d.shape[0] >= 52:
        newdfs_training.append(d)

#Randomly sampling assets
random.seed(20)
newdfs_training = random.sample(newdfs_training, 1000)

# ==============================================================
# Filtering those dataframes with 2017 pricing data for testing

newdfs_test = []

for d in pricedfs:
    d['Date'] = pd.to_datetime(d['Date'])
    d['year'] = d['Date'].dt.year
    d['month'] = d['Date'].dt.month
    d['day'] = d['Date'].dt.day

    d = d[(d['year'] == 2017) & ((d['month'] != 10) & (d['month'] != 11) &
                                 (d['month'] != 12))].drop(columns=['year', 'month', 'day'])
    d['Date'] = d['Date'].dt.strftime('%Y-%m-%d')
    if d.shape[0] >= 39:
        newdfs_test.append(d)

#Randomly sampling assets
random.seed(20)
newdfs_test = random.sample(newdfs_test, 1000)

print("Dataset Divided into Training and Test Sets")
```

## Feature Creation for the Model

Now that we have the pricing data in a more useful form, we can now convert that data into additional indicators that a machine learned model can use for identifying patterns/trends. In effect, we want to capture how the price for an asset changed in the recent past, for use as indicators for future performance (of course past performance is not always a good indicator, and more advanced approaches may mix in other sources of evidence here). We convert the pricing data into three main indicator (feature) types:

 Average closing price

 Return on Investment

 Volatility

These are calculated over the last 3, 6, and 9 months from the point in time at which recommendations are requested.

Returns are calculated as follows over each specified period:

\begin{equation}
Return On Investment = \frac{Final \ closing \ price \ - \ Initial \ closing \ price}{Initial \ closing \ price}
\end{equation}

Volatility represent the risk of a stock as expressed by its fluctuations, and is expressed as the standard deviation of the returns of the stock. As our returns are calculated weekly, we also ensure our volatilities are converted into 3-month, 6-month, and 9-month estimates.

\begin{equation}
Volatility = \sqrt{\frac{1}{N-1} \sum_{i=1}^N (Return_i - \overline{Return})^2} * \sqrt{Number \ of \ weeks}
\end{equation}

```python

def calculate_features(newdfs):
    dfs = []
    for d in newdfs:
        d['Date'] = pd.to_datetime(d['Date'])
        df3m = d.set_index('Date').last('3M')
        d['avg_price_3m'] = df3m['Close'].mean()
        d['volatility_3m'] = df3m['returns'].std() * math.sqrt(13.0357)
        d['returns_3m'] = (df3m.tail(1)['Close'].values[0] - df3m.head(1)['Close'].values[0])/df3m.head(1)['Close'].values[0]

        df6m = d.set_index('Date').last('6M')
        d['avg_price_6m'] = df6m['Close'].mean()
        d['volatility_6m'] = df6m['returns'].std() * math.sqrt(26.0715)
        d['returns_6m'] = (df6m.tail(1)['Close'].values[0] - df6m.head(1)['Close'].values[0])/df6m.head(1)['Close'].values[0]

        df9m = d.set_index('Date').last('9M')
        d['avg_price_9m'] = df9m['Close'].mean()
        d['volatility_9m'] = df9m['returns'].std() * math.sqrt(39.1072)
        d['returns_9m'] = (df9m.tail(1)['Close'].values[0] - df9m.head(1)['Close'].values[0])/df9m.head(1)['Close'].values[0]

        d = d.head(1)
        dfs.append(d)

    all_feature_dfs = pd.concat(dfs)
    all_feature_dfs = all_feature_dfs.drop(columns=['Date', 'Close', 'returns'])
    all_feature_dfs.dropna(inplace=True)
    return all_feature_dfs

training_data = calculate_features(newdfs_training)

# Enter your training dataset path here
training_data.to_csv(storageDIR+'/Kaggle-2016-Training-Final.csv')

test_data = calculate_features(newdfs_test)

# Enter your test dataset path here
test_data.to_csv(storageDIR+'/Kaggle-2017-Test-Final.csv')

print("Dataset Pricing Trend Features Extracted")

```

    Dataset Pricing Trend Features Extracted

## Profitability prediction model - Random forest regression

After having curated our training and testing sets, we proceed to use a simple machine learning model to predict the 9-month profitability of assets from the test dataset. The model we choose is a Random Forest model. This is imported from the Scikit-Learn library, which is a commonly used package for introductory machine learning experiments and houses a wide variety of models. A random forest regression model is an ensemble model, which means that it combines the predictions of multiple machine learning models in order to improve prediction accuracy and avoid overfitting. A random forest model aggregates decision trees and takes the average of their predictions or probabilities for classes; decision trees are supervised models that attempt to infer simple decision rules from the features and use these to make predictions.

Random forests introduce randomness from two other sources, namely that each tree is built with a random sample from the training set, and that the splitting of nodes during tree construction is either found from all features or from a random subset of the maximum permitted number of features. These sources of randomness are intended to minimize the impact of each singular decision tree's sensitivity.

Here, we use the training and test data we just prepared and apply a scaler to them so that their values are standardized for prediction. We then apply a random forest regressor model with 500 estimators (decision trees) to calculate predicted returns over a 9-month period.

```python
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler

# Using precomputed data
training_data = pd.read_csv(storageDIR+'/Kaggle-2016-Training-Final.csv')
test_data = pd.read_csv(storageDIR+'/Kaggle-2017-Test-Final.csv')

train_instrument_ids = training_data['Stock']
returns_9m = training_data['returns_9m']
training_data = training_data.drop(columns=['Stock', 'returns_9m'])

test_instrument_ids = test_data['Stock']
test_returns_9m = test_data['returns_9m']
test_vol_9m = test_data['volatility_9m']
test_data = test_data.drop(columns=['Stock', 'returns_9m'])

# Scaling data values
scaler = StandardScaler()
scaled_data = scaler.fit_transform(training_data)

scaler = StandardScaler()
scaled_test = scaler.fit_transform(test_data)

# Prediction
model = RandomForestRegressor(n_estimators=500)
model.fit(scaled_data, returns_9m)
pred = model.predict(scaled_test)

pred_df = pd.DataFrame(pred)
pred_df = pred_df.rename(columns={0: 'Predicted returns'})
pred_df['Stock'] = pd.Series(test_instrument_ids)
pred_df['Actual returns (9M)'] = test_returns_9m
pred_df['Actual Volatility (9M)'] = test_vol_9m
```

## Model Effectiveness Evaluation

After the prediction of returns, we proceed onto the evaluation of our model. We can do so by means of more classical metrics that assess the difference between actual and predicted values; we can also look at the aggregate profitability and volatility of the recommended assets over the test period as a measure of suitability.

```python
pred_df = pred_df.sort_values(by='Predicted returns', ascending=False)
pred_df.head(10)
```

      Predicted returns
      Stock
      Actual returns (9M)
      Actual Volatility (9M)

      298
      2.704236
      Stocks\zgnx
      2.050980
      1.940495

      224
      2.287230
      Stocks\abeo
      1.933884
      0.860671

      506
      2.263269
      Stocks\arwr
      1.354497
      0.631877

      429
      2.251566
      Stocks\sppi
      2.081023
      0.620670

      398
      2.226679
      Stocks\wprt
      1.488189
      0.726520

      402
      2.224329
      Stocks\cyrx
      1.966049
      0.636090

      150
      2.199163
      Stocks\ensv
      0.037037
      0.849338

      705
      2.107315
      Stocks\xxii
      2.140152
      0.738352

      772
      2.102616
      Stocks\qure
      0.761017
      0.604679

      513
      2.093745
      Stocks\sify
      0.942658
      0.496918

In the above table, we are showing the top 10 stocks that were predicted to be profitable. The last two columns report the actual return on investment after 9 months and asset volitility. Note that a return value of 1.5 means a 150% return on investment.

We notice that predicted returns for the top stocks are exceedingly high, which is not ordinary. However, we can also see that the actual returns for these stocks are similar for several of these instances, i.e. the model is not wrong in predicting these as profitable investments in the short term. However, we can also see that the volitility fo these stocks is very high, i.e. these are 'high-risk' assets that may subsequently crash in price.

We can also analyse the statistics for this predictions across the dataset.

```python
pred_df.describe()
```

      Predicted returns
      Actual returns (9M)
      Actual Volatility (9M)

      count
      1000.000000
      1000.000000
      1000.000000

      mean
      0.291667
      0.117710
      0.279499

      std
      0.545248
      0.440275
      0.270453

      min
      -0.858322
      -0.998589
      0.008209

      25%
      0.019777
      -0.027597
      0.108639

      50%
      0.199883
      0.073051
      0.210807

      75%
      0.463422
      0.214930
      0.364167

      max
      2.704236
      7.232558
      3.707981

Here we are able to see the average returns and volatility of the entire test set.

We can also examine the predicted returns and volatility over the top 1, 5, 10, and 20 stocks, ranked by their predicted profitability. This will serve to provide an idea of how financially feasible our recommendations are, and whether the recommended stocks are an improvement to the average returns of the test set.

```python

ret_vol = pd.DataFrame({'Number of assets': [1, 5, 10, 20], 'Mean returns (9M), in %': [pred_df.head(1)['Actual returns (9M)'].mean(),
                        pred_df.head(5)['Actual returns (9M)'].mean(),
                        pred_df.head(10)['Actual returns (9M)'].mean(),
                        pred_df.head(20)['Actual returns (9M)'].mean()],
                        'Mean volatility (9M), in %': [pred_df.head(1)['Actual Volatility (9M)'].mean(),
                                                pred_df.head(5)['Actual Volatility (9M)'].mean(),
                                                pred_df.head(10)['Actual Volatility (9M)'].mean(),
                                                pred_df.head(20)['Actual Volatility (9M)'].mean()]
                       })
ret_vol['Mean returns (9M), in %'] *= 100
ret_vol['Mean volatility (9M), in %'] *= 100

ret_vol
```

      Number of assets
      Mean returns (9M), in %
      Mean volatility (9M), in %

      0
      1
      205.098039
      194.049514

      1
      5
      178.171489
      95.604667

      2
      10
      147.554883
      81.056093

      3
      20
      123.460761
      77.108316

The returns and volatility for the top stocks, ranked by predicted returns, are far higher than their averages across the test set. This indicates that ranking assets by their predicted returns can produce some highly profitable but risk-laden investment recommendations, which might be suitable for aggressive investors. However, it remains to be seen how much of this is owed to fluctuations and outliers in the data, and perhaps even if there are better ways to capture the returns and volatility of the dataset.

Next, we look at the differences between the actual and predicted returns.

```python
import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(10,10))
sns.regplot(pred_df['Actual returns (9M)'], pred_df['Predicted returns'])
plt.show()
```

```python
plt.figure(figsize=(10,10))
plt.plot(pred_df['Actual returns (9M)'] - pred_df['Predicted returns'] ,marker='o', linestyle='', c='crimson')
plt.xlabel('Observations', fontsize=15)
plt.ylabel('Difference between true and predicted returns', fontsize=15)
plt.show()
```

Lastly, we can examine the mean absolute error and mean squared error of the predictions. As these can be quite dependent on the dataset and problem in question, we also assume a simple baseline, by taking the median of all stock returns from the test dataset. We then compare the results of applying these metrics to the baseline and our predictor model.

```python
from sklearn.metrics import mean_absolute_error, mean_squared_error

testlist = pd.Series([test_returns_9m.median()]*1000)

mae_mse = pd.DataFrame({'Metric': ['MAE', 'MSE'],
                        'Baseline (median return)': [mean_absolute_error(pred, testlist),
                                                     mean_squared_error(pred, testlist)],
                       'Random forest regression': [mean_absolute_error(pred, test_returns_9m),
                                                   mean_squared_error(pred, test_returns_9m)]
                       })

mae_mse['Percentage difference'] = (mae_mse['Random forest regression'] -
                                    mae_mse['Baseline (median return)'])/mae_mse['Baseline (median return)']

mae_mse
```

      Metric
      Baseline (median return)
      Random forest regression
      Percentage difference

      0
      MAE
      0.394475
      0.259137
      -0.343083

      1
      MSE
      0.344792
      0.200128
      -0.419568

We can see from this that the random forest model presents an improvement (reduction) in both MAE and MSE.

```python

```

ext: ipynb


text:
```python
# h3: 3D Homogeneous Coordinates

import sympy
from sympy import *    # Make SymPy available to this program
sympy.init_printing()

from galgebra.printer import Fmt, GaPrinter # Need all these parameters
from galgebra.gprinter import gFormat, gprint
gFormat()

from galgebra.ga import *  # Make GAlgebra available to this program.
from galgebra.mv import *
from galgebra.printer import Fmt, GaPrinter, Format
    # Fmt:       sets the way that a multivector's basis expansion is output.
    # GaPrinter: makes GA output a little more readable.
    # Format:    turns on latex printer.

from galgebra.gprinter import gFormat, gprint
gFormat()
```

$\displaystyle \DeclareMathOperator{\Tr}{Tr}$$
$$\DeclareMathOperator{\Adj}{Adj}$$
$$\newcommand{\bfrac}[2]{\displaystyle\frac{#1}{#2}}$$
$$\newcommand{\lp}{\left (}$$
$$\newcommand{\rp}{\right )}$$
$$\newcommand{\paren}[1]{\lp {#1} \rp}$$
$$\newcommand{\half}{\frac{1}{2}}$$
$$\newcommand{\llt}{\left <}$$
$$\newcommand{\rgt}{\right >}$$
$$\newcommand{\abs}[1]{\left |{#1}\right | }$$
$$\newcommand{\pdiff}[2]{\bfrac{\partial {#1}}{\partial {#2}}}$$
$$\newcommand{\npdiff}[3]{\bfrac{\partial^{#3} {#1}}{\partial {#2}^{#3}}}$$
$$\newcommand{\lbrc}{\left \{}$$
$$\newcommand{\rbrc}{\right \}}$$
$$\newcommand{\W}{\wedge}$$
$$\newcommand{\prm}[1]{{#1}}$$
$$\newcommand{\ddt}[1]{\bfrac{d{#1}}{dt}}$$
$$\newcommand{\R}{\dagger}$$
$$\newcommand{\deriv}[3]{\bfrac{d^{#3}#1}{d{#2}^{#3}}}$$
$$\newcommand{\grade}[2]{\left < {#1} \right >_{#2}}$$
$$\newcommand{\f}[2]{{#1}\lp {#2} \rp}$$
$$\newcommand{\eval}[2]{\left . {#1} \right |_{#2}}$$
$$\newcommand{\bs}[1]{\boldsymbol{#1}}$$
$$\newcommand{\grad}{\bs{\nabla}}$

```python
h3coords = (x,y,z,e) = symbols('x y z e', real=True)
# h3 = Ga('e_x e_y e_z e_e', g=[1,1,1,1], coords=h3coords)
h3 = Ga('\mathbf{e}', g=[1,1,1,1], coords=h3coords)
(ex,ey,ez,ee) = h3.mv()
```

```python

```


text:
    Foundations of Data Science 7COM1073

        University of Hertfordshire

Assessment 1 – Data Classification (75 marks)

Mohammed Saaduddin Ahmed, MSc (SW) Data Science and Analytics, 19001506

This assessment includes two parts and counts for 75% of the overall assessment for this module.
The first part is 20 marks and should be submitted by noon on 22/11/2019 and the second one is 55 marks and should be handed by 20/12/2019.

The programming language you should use to finish this assessment is Python (in version 3 and above). In particular, you can use functions from the following packages: Numpy, Pandas, Matplotlib, Seaborn and Sklearn.

All Python skills needed to do this assessment have been covered in the practical sessions –practical notes are available on Canvas.

The information of the dataset which you will work with can be viewed in the following link: https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/
 Note that you should work on the red wine dataset (winequality-red.csv) only, which can be downloaded from ‘Data Folder’ in the link given above.

Part one (20 marks)

***

### **Task 1:** Data pre-processing and data exploration (15 marks)

- a. Use Pandas to load data
- b. Merge all the data with “quality” labels between 6-10 into Class 1 and similarly form Class 2 for the data with “quality” labels between 1-5.
- c. Report the number of features and number of rows in each class
- d. Choose an attribute and generate a boxplot for the two pre-defined classes.
- e. Show one scatter plot, that is, one feature against another feature. It is your choice to show which two features you want to use.

***

##### *a. Use Pandas to load data*

```python
import pandas as pd
wine = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", delimiter= ";")
#wine = pd.read_csv("winequality-red.csv")
print("Shape of the red wine dataset is :",wine.shape)
wine.head()
#wine.dtypes
```

    Shape of the red wine dataset is : (1599, 12)

      fixed acidity
      volatile acidity
      citric acid
      residual sugar
      chlorides
      free sulfur dioxide
      total sulfur dioxide
      density
      pH
      sulphates
      alcohol
      quality

      0
      7.4
      0.70
      0.00
      1.9
      0.076
      11.0
      34.0
      0.9978
      3.51
      0.56
      9.4
      5

      1
      7.8
      0.88
      0.00
      2.6
      0.098
      25.0
      67.0
      0.9968
      3.20
      0.68
      9.8
      5

      2
      7.8
      0.76
      0.04
      2.3
      0.092
      15.0
      54.0
      0.9970
      3.26
      0.65
      9.8
      5

      3
      11.2
      0.28
      0.56
      1.9
      0.075
      17.0
      60.0
      0.9980
      3.16
      0.58
      9.8
      6

      4
      7.4
      0.70
      0.00
      1.9
      0.076
      11.0
      34.0
      0.9978
      3.51
      0.56
      9.4
      5

##### *b. Merge all the data with “quality” labels between 6-10 into Class 1 and similarly form Class 2 for the data with “quality” labels between 1-5.*

```python
Class1 = wine.loc[(wine['quality'] >=6) & (wine['quality'] <=10)]
Class2 = wine.loc[(wine['quality'] >=0) & (wine['quality'] <=5)]
print("Number of observations in Class 1 are",Class1['quality'].count(),"and Number of observations in Class 2 are",Class2['quality'].count())

wine["quality"] = wine["quality"].replace([6,7,8,9,10],'Class 1')
wine["quality"] = wine["quality"].replace([1,2,3,4,5],'Class 2')
```

    Number of observations in Class 1 are 855 and Number of observations in Class 2 are 744

##### *c. Report the number of features and number of rows in each class*

```python
#c Report the number of features and number of rows in each class
print("Class 1 consists of",Class1.shape[0],"rows and",Class1.shape[1],"features")
print("Class 2 consists of",Class2.shape[0],"rows and",Class2.shape[1],"features")
```

    Class 1 consists of 855 rows and 12 features
    Class 2 consists of 744 rows and 12 features

##### *d. Choose an attribute and generate a boxplot for the two pre-defined classes*

```python
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline

data_to_plot = [Class1['density'], Class2['density']]

fig = plt.figure(1, figsize=(9, 6))
ax = fig.add_subplot(111)
ax.boxplot(data_to_plot)
plt.title(' Boxplot of density in Class 1 and Class 2 ',fontsize=15)
plt.xlabel("Class 1 and Class 2")
plt.ylabel("density "r'$( g/cm^3 )$')

fig.savefig('Boxplot of density in Class 1 and Class 2.png', bbox_inches='tight')
```

##### *e. Show one scatter plot, that is, one feature against another feature. It is your choice to show which two features you want to use.*

```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

fig = plt.figure(1, figsize=(9, 6))
ax = fig.add_subplot(111)
area = np.pi*7

plt.scatter(wine['free sulfur dioxide'], wine['total sulfur dioxide'], s = area , alpha=0.3)
plt.title(' Scatter plot of Free sulfur diaoxide and Total sulfur dioxide',size = 15)
plt.xlabel("Free sulfur dioxide "r'$(mg/dm^3)$')
plt.ylabel("Total sulfur dioxide "r'$(mg/dm^3)$')
plt.show()

fig.savefig('Scatter plot of Free sulfur diaoxide and Total sulfur dioxide.png', bbox_inches='tight')
```

***

### **Task 2:**  Computing probabilities using Python code for the given red wine dataset (5 marks)

- f. Prior probability:
        i. What is the probability of a wine classified as Class 1 (P(Class 1))?
       ii. What is the probability of a wine classified as Class 2 (P(Class 2))?
- g. Conditional probability:
        i. What is the probability of a wine having a pH value greater than 3.6 given it is classified as Class 1 (P(pH>3.6|Class 1))?
- h. Posterior probability:
        i.  What is the probability of a wine classified as Class 1 when it has a pH value greater than 3.6?

f. Prior probability:

   + *i. What is the probability of a wine classified as Class 1 (P(Class 1))?*
   + *ii. What is the probability of a wine classified as Class 2 (P(Class 2))?*

```python
# f. Prior probablitiy
import numpy as np
import sympy as sp
from sympy import FiniteSet

SS = SampleSpace = wine['quality'].count()
C1 = Class1['quality'].count()
C2 = Class2['quality'].count()

PC1 = C1/SampleSpace
PC2 = C2/SampleSpace

PC1_round = round(PC1,2)
PC2_round = round(PC2,2)

# i. What is the probability of a wine classified as Class 1 (P(Class 1))?
print("The probability of a wine classified as Class 1 =",PC1_round)

# ii. What is the probability of a wine classified as Class 2 (P(Class 2))?
print("The probability of a wine classified as Class 2 =",PC2_round)
```

    The probability of a wine classified as Class 1 = 0.53
    The probability of a wine classified as Class 2 = 0.47

Note :   Results are rounded upto two decimal places.

Number of observations in: SampleSpace = 1599, Class1 = 855, Class2 = 744

\begin{equation*}
P(Class 1) = \frac{855}{1599} = 0.53 \\ P(Class 2) = \frac{744}{1599} = 0.47
\end{equation*}

g. Conditional probability:

+ *i. What is the probability of a wine having a pH value greater than 3.6 given it is classified as Class 1 (P(pH>3.6|Class 1))?*

```python
#g. Conditional Probablity
#i. What is the probability of a wine having a pH value greater than 3.6 given it is classified as Class 1 (P(pH>3.6|Class 1))?
import pandas as pd

DataFrame = pd.DataFrame(Class1)
subsetDataFrame = DataFrame[DataFrame['pH'] > 3.6]

# P_P = P(pH>3.6 and Class 1) =>  probability of a wine having a pH value greater than 3.6 and are classified as Class 1
# P_C = P(Class1) => probabilty of a wine being classified as Class 1
# P_probability => (P(pH>3.6|Class 1))

P_P = subsetDataFrame['pH'].count()
print("The count of values which have pH > 3.6 and are classified as Class 1 :", P_P)
P_P = P_P/SampleSpace
P_C = Class1['pH'].count()
P_C = P_C/SampleSpace
P_probability = (P_P / P_C)

P_probability_round = round(P_probability,2)

print("\nThe probability of a wine having a pH value greater than 3.6 given it is classified as Class 1 (P(pH>3.6|Class 1)) :",P_probability_round)

```

    The count of values which have pH > 3.6 and are classified as Class 1 : 25

    The probability of a wine having a pH value greater than 3.6 given it is classified as Class 1 (P(pH>3.6|Class 1)) : 0.03

Note :  Results are rounded upto two decimal places.

Number of observations in: SampleSpace = 1599, Class1 = 855, Class2 = 744`

\begin{equation*}
                Conditional\;Probability: P(A|B) = \frac{P(A \cap B)}{P(B)}
\end{equation*}

$$ P(pH>3.6|Class\;1) = \frac{\frac{25}{1599}}{\frac{855}{1599}} = \frac{25}{855} = 0.03 $$

\begin{equation*}
Result : P(pH>3.6|Class\;1) = 0.03
\end{equation*}

h. Posterior probability:

+ *i.  What is the probability of a wine classified as Class 1 when it has a pH value greater than 3.6 (P(Class 1|pH>3.6))?*

```python
#h. Posterior probability:
#i.  What is the probability of a wine classified as Class 1 when it has a pH value greater than 3.6 P(Class 1|pH>3.6)?
import pandas as pd

DataFrame = pd.DataFrame(wine)
subsetDataFrame = DataFrame[DataFrame['pH'] > 3.6]

# P_P = P(pH>3.6) =>  probability of a wine having a pH value greater than 3.6
# P_C = P(Class1) => probabilty of a wine being classified as Class 1
# P_probability => (P(Class 1|pH>3.6))

P_P = subsetDataFrame['pH'].count()
print("The count of values which have pH > 3.6 :", P_P)
P_P = P_P/SampleSpace
P_C = Class1['pH'].count()
P_C = P_C/SampleSpace
P_probability_posterior = ((P_probability)*P_C)/P_P
P_probability_round = round(P_probability_posterior,2)

print("\nThe probability of a wine classified as Class 1 when it has a pH value greater than 3.6 (P(Class 1|pH>3.6)) :",P_probability_round)

```

    The count of values which have pH > 3.6 : 48

    The probability of a wine classified as Class 1 when it has a pH value greater than 3.6 (P(Class 1|pH>3.6)) : 0.52

Note :  Results are rounded upto two decimal places.

Number of observations in: SampleSpace = 1599, Class1 = 855, Class2 = 744, pH>3.6 = 48, pH>3.6∩Class1 = 25

\begin{equation*}
                Posterior\;Probability: P(B|A) = \frac{P(A|B)P(B)}{P(A)}
\end{equation*}

\begin{equation*}
P(Class\;1|pH>3.6) = \frac{P(pH>3.6|Class\;1)P(Class\;1)}{P(pH>3.6)}
\end{equation*}

\begin{equation*}
P(Class\;1|pH>3.6) = \frac{\frac{25}{85}\times{\frac{855}{1599}}}{\frac{48}{1599}} = \frac{25}{48} = 0.52
\end{equation*}

\begin{equation*}
Result : (Class\;1|pH>3.6) = 0.52
\end{equation*}

****

max_stars_repo_path: Assingments/Portfolio Assessment/Assessment 1 - Data Classification.ipynb

max_stars_repo_name: simplesaad/FoundationsOfDataScience

max_issues_repo_name: simplesaad/FoundationsOfDataScience

max_forks_repo_name: simplesaad/FoundationsOfDataScience


text:
```python
import numpy as np
import sympy as sym
import numba
import pydae.build as db
```

```python

```

\begin{eqnarray}
\dot \delta &=& \Omega_b \left(\omega - \omega_s\right)\\
\dot \omega &=& 1/(2 H) \left(p_m - p_e - D (\omega - \omega_s) \right)\\
\end{eqnarray}

$$ \sf
\Omega_{b} \left(\sf \omega - \omega_{s}\right)
$$

## System definition

```python
params = {'X_d':1.81,'X1d':0.3,'T1d0':8.0,   # synnchronous machine d-axis parameters
          'X_q':1.76,'X1q':0.65,'T1q0':1.0,  # synnchronous machine q-axis parameters
          'R_a':0.003,'X_l': 0.05,
          'H':6.0,'D':1.0,'S_n':1.0,
          'Omega_b':2*np.pi*50,'omega_s':1.0,
          'K_a':100, 'T_r':0.1, 'v_pss':0.0,
          'Droop':0.05,'T_m':5.0,'K_sec':0.001,
          'K_delta':0.01}

u_ini_dict = {'P_l':0.8, 'Q_l':0.2,'P_cig':0.1,'Q_cig':0.0,'v_ref':1.0}  # for the initialization problem
u_run_dict = {'P_l':0.8, 'Q_l':0.2,'P_cig':0.1,'Q_cig':0.0,'v_ref':1.0}  # for the running problem (here initialization and running problem are the same)

x_list = ['delta','omega','e1q','e1d','v_c','p_m','xi_m']    # [inductor current, PI integrator]
y_ini_list = ['i_d','i_q','v_1','theta_1','P_t','Q_t','v_f','p_m_ref','v_0','theta_0'] # for the initialization problem
y_run_list = ['i_d','i_q','v_1','theta_1','P_t','Q_t','v_f','p_m_ref','v_0','theta_0'] # for the running problem (here initialization and running problem are the same)

sys_vars = {'params':params,
            'u_list':u_run_dict,
            'x_list':x_list,
            'y_list':y_run_list}

exec(db.sym_gen_str())  # exec to generate the required symbolic varables and constants
```

```python
v_d = v_1*sin(delta - theta_1)
v_q = v_1*cos(delta - theta_1)

p_e = i_d*(v_d + R_a*i_d) + i_q*(v_q + R_a*i_q)

ddelta = Omega_b*(omega - omega_s) - K_delta*delta
domega = 1/(2*H)*(p_m - p_e - D*(omega - omega_s))
de1q = 1/T1d0*(-e1q - (X_d - X1d)*i_d + v_f)
de1d = 1/T1q0*(-e1d + (X_q - X1q)*i_q)
dv_c =   (v_1 - v_c)/T_r
dp_m =   (p_m_ref - p_m)/T_m
dxi_m =   omega - 1

g_1  = v_q + R_a*i_q + X1d*i_d - e1q
g_2  = v_d + R_a*i_d - X1q*i_q - e1d
g_3  = P_t*S_n + P_cig - (v_1*v_0*sin(theta_1 - theta_0))/X_l
g_4  = Q_t*S_n + Q_cig + (v_1*v_0*cos(theta_1 - theta_0))/X_l - v_1**2/X_l
g_5  = i_d*v_d + i_q*v_q - P_t
g_6  = i_d*v_q - i_q*v_d - Q_t
g_7  = K_a*(v_ref - v_c + v_pss) - v_f
g_8  = -p_m_ref - K_sec*xi_m - 1/Droop*(omega - 1)
g_9  = -P_l - (v_0*v_1*sin(theta_0 - theta_1))/X_l
g_10 = -Q_l + (v_0*v_1*cos(theta_0 - theta_1))/X_l - v_0**2/X_l

sys = {'name':'iso_milano_ex8p1_4ord_avr_gov',
       'params_dict':params,
       'f_list':[ddelta,domega,de1q,de1d,dv_c,dp_m,dxi_m],
       'g_list':[g_1,g_2,g_3,g_4,g_5,g_6,g_7,g_8,g_9,g_10],
       'x_list':x_list,
       'y_ini_list':y_ini_list,
       'y_run_list':y_run_list,
       'u_run_dict':u_run_dict,
       'u_ini_dict':u_ini_dict,
       'h_dict':{'p_m':p_m,'P_l':P_l,'Q_l':Q_l,'P_cig':P_cig,'Q_cig':Q_cig}}

sys = db.system(sys)
db.sys2num(sys)

```

```python

```

```python

```

```python

```

max_stars_repo_name: pydae/pydae

max_issues_repo_name: pydae/pydae

max_forks_repo_name: pydae/pydae


text:
# Inaugural Project

Imports and set magics:

```python
import numpy as np
import math
from scipy import optimize

# autoreload modules when code is run
%load_ext autoreload
%autoreload 2

# local modules
import inauguralproject
```

## Question 1

We start out by setting up the consumer mazimization problem. This is done by setting up a function returning the utility function and the objective function. These functions are set up in the inauguralproject.py file, to give a better overview of this file.

We then call upon a solver that can solve the maximization problem given the parameters.

```python
# call solver
def sol_case1(w, epsilon, m, nu, tau0, tau1, kappa):
    sol_case1_opt = optimize.minimize_scalar(inauguralproject.value_of_choice,method='bounded',
        bounds=(0,1),args=(w, epsilon, m, nu, tau0, tau1, kappa))
    l_opt = sol_case1_opt.x
    return l_opt
    l=l_opt
    u = inauguralproject.u_func(l, w, epsilon, m, nu, tau0, tau1, kappa)

```

## Question 2

Plot $l^{*}$and $c^{*}$ as a function of w in the range 0.5 to 1.5.

We create vectors for wage, labour and consumption, where the labour and consumptions vectors are empty. The labour and concumption is calculated for every value of w in a for loop.
From this we are able to plot the optimal values of labour and consumption

```python
N = 100
# vector of wages
w_vec = np.linspace(0.5,1.5,N)
l_vec = np.empty(N)
c_vec = np.empty(N)

# calculates l and c for every value of w
for i,w in enumerate(w_vec):
    w = w_vec[i]
    find_l=sol_case1(w, epsilon=0.3, m=1, nu=10, tau0=0.4, tau1=0.1, kappa=0.4)
    l_vec[i] = find_l
    l=l_vec[i]
    c_vec[i] = inauguralproject.c_func(l, w)

# creates the figures
import matplotlib.pyplot as plt
plt.style.use("seaborn")

fig = plt.figure(figsize=(10,4))

fig_w_l = fig.add_subplot(1,2,1)
fig_w_l.plot(w_vec,l_vec)

fig_w_l.set_title('Optimal $l$ as a funtion of $w$')
fig_w_l.set_xlabel('$w$')
fig_w_l.set_ylabel('$l*$')
fig_w_l.grid(True)

fig_w_c = fig.add_subplot(1,2,2)
fig_w_c.plot(w_vec,c_vec)

fig_w_c.set_title('Optimal $c$ as a funtion of $w$')
fig_w_c.set_xlabel('$w$')
fig_w_c.set_ylabel('$c*$')
fig_w_c.grid(True)
```

## Question 3

Condsider a population with $N=10.000$ individuals indexed by i.
Assume the distribution of wages is uniform such that:

\begin{equation}
w_{i}=\bigcup(0.5,1.5)
\end{equation}

Denote the optimal choices of individual i by $l_{i}^{*}$and $c_{i}^{*}.$

Calculate the total tax revenue given by:

\begin{equation}
T=\sum_{i=1}^{N}\left[\tau_{0}w_{i}l_{i}^{*}+\tau_{1}max\{wl-\kappa,0\}]\right]
\end{equation}

We set seed to be able to draw the same uniform numbers everytime we run the code. This time we compute the vector of wages by uniformly drawn numbers in the range 0.5 to 1.5 for N2 entries.

Then we compute two functions, which are loacted in the py.file. The first function is a function calculating the taxes paid by each individual while the second function returns the total tax revenue. Thus, the total taxrevenue is given by:

```python
np.random.seed(42)
N2 = 100 #OBS! Was supposed to be 10K, but code took forever to run in Q5
w_vec2 = np.random.uniform(0.5,1.5,N2)
l_vec2 = np.empty(N2)
t_vec = np.empty(N2)

print(f'The total tax revenue is given by {inauguralproject.tax_rev(epsilon=0.3, m=1, nu=10, tau0=0.4, tau1=0.1, kappa=0.4):.2f} ')

```

    The total tax revenue is given by 12.25

## Question 4

What would the tax revenue be instead of $\varepsilon=0.1$
We specify a new level of epsilon and call the total tax revnue function again:

```python
print(f'Given the new value of epsilon the tax revenue becomes {inauguralproject.tax_rev(epsilon=0.1, m=1, nu=10, tau0=0.4, tau1=0.1, kappa=0.4):.2f} ')
```

    Given the new value of epsilon the tax revenue becomes 28.18

## Question 5

Consider a politician who wishes to maximize the tax revenue.

Which $\tau_{0},\,\tau_{1}$ and $\kappa$ would you suggest her to
implement? Report the tax revenue you expect to obtain.

To solve this we define a new function that returns the objective function, that we want to maximize. To be able to solve this optimization problem we set some bounds and initial guesses.

From this we can optimize the objective function and use the optimal values of $\tau_{0},\,\tau_{1}$ and $\kappa$ to calculate the taxrevenue that we expect the politician to obtain

```python
# Setting bounds and initial guess
# Tax rates cannot be higher than 100 pct., kappa cannot be higher than max(l)*max(w)
bounds_tax = ((0,1),(0,1),(0,1.5))
initval = [0.4,0.1,0.4]

# Solving by calling the objective function
sol_opt_tax = optimize.minimize(inauguralproject.obj_func, initval, method='SLSQP',bounds=bounds_tax)

# Unpack solution
tau0_opt = sol_opt_tax.x[0]
tau1_opt = sol_opt_tax.x[1]
kappa_opt = sol_opt_tax.x[2]

print(f'Optimal tau0 is {tau0_opt:.3f}, optimal tau1 is {tau1_opt:.3f} and optimal kappa is {kappa_opt:.3f}')
print(f'Total tax revenue in this case will be {inauguralproject.tax_rev(0.1,1,10,tau0_opt,tau1_opt,kappa_opt):.3f}')
```

    Optimal tau0 is 0.914, optimal tau1 is 0.000 and optimal kappa is 0.000
    Total tax revenue in this case will be 50.551

# Conclusion

We have constructed a function that maximizes the consumer problem. We have further computed a function that can return the total tax revenue and which futrher allows ud to adjust the level of $\epsilon$.

Finally we have helped the politician by suggesting the values of $\tau_{0},\,\tau_{1}$ and $\kappa$, which yields the expected taxrevenue of 50.551.


text:
# Quantization of Signals

*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*

## Characteristic of a Linear Uniform Quantizer

The characteristics of a quantizer depend on the mapping functions $f(\cdot)$, $g(\cdot)$ and the rounding operation $\lfloor \cdot \rfloor$ introduced in the [previous section](introduction.ipynb). A linear quantizer bases on linear mapping functions $f(\cdot)$ and $g(\cdot)$. A uniform quantizer splits the mapped input signal into quantization steps of equal size. Quantizers can be described by their nonlinear in-/output characteristic $x_Q[k] = \mathcal{Q} \{ x[k] \}$, where $\mathcal{Q} \{ \cdot \}$ denotes the quantization process. For linear uniform quantization it is common to differentiate between two characteristic curves, the so called mid-tread and mid-rise.

### Mid-Tread Chacteristic Curve

The in-/output relation of the mid-tread quantizer is given as

\begin{equation}
x_Q[k] = Q \cdot \underbrace{\left\lfloor \frac{x[k]}{Q} + \frac{1}{2} \right\rfloor}_{index}
\end{equation}

where $Q$ denotes the quantization step size and $\lfloor \cdot \rfloor$ the [floor function](https://en.wikipedia.org/wiki/Floor_and_ceiling_functions) which maps a real number to the largest integer not greater than its argument. Without restricting $x[k]$ in amplitude, the resulting quantization indexes are [countable infinite](https://en.wikipedia.org/wiki/Countable_set). For a finite number of quantization indexes, the input signal has to be restricted to a minimal/maximal amplitude $x_\text{min} < x[k] < x_\text{max}$ before quantization. The resulting quantization characteristic of a linear uniform mid-tread quantizer is shown in the following

The term mid-tread is due to the fact that small values $|x[k]| < \frac{Q}{2}$ are mapped to zero.

#### Example

The quantization of one period of a sine signal $x[k] = A \cdot \sin[\Omega_0\,k]$ by a mid-tread quantizer is simulated. $A$ denotes the amplitude of the signal, $x_\text{min} = -1$ and $x_\text{max} = 1$ are the smallest and largest output values of the quantizer, respectively.

```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

A = 1.2  # amplitude of signal
Q = 1/10  # quantization stepsize
N = 2000  # number of samples

def uniform_midtread_quantizer(x, Q):
    # limiter
    x = np.copy(x)
    idx = np.where(np.abs(x) >= 1)
    x[idx] = np.sign(x[idx])
    # linear uniform quantization
    xQ = Q * np.floor(x/Q + 1/2)

    return xQ

def plot_signals(x, xQ):
    e = xQ - x
    plt.figure(figsize=(10,6))
    plt.plot(x, label=r'discrete signal $x[k]$')
    plt.plot(xQ, label=r'quantized signal $x_Q[k]$')
    plt.plot(e, label=r'quantization error $e[k]$')
    plt.xlabel(r'$k$')
    plt.axis([0, N, -1.1*A, 1.1*A])
    plt.legend()
    plt.grid()

# generate signal
x = A * np.sin(2*np.pi/N * np.arange(N))
# quantize signal
xQ = uniform_midtread_quantizer(x, Q)
# plot signals
plot_signals(x, xQ)
```

**Exercise**

* Change the quantization stepsize `Q` and the amplitude `A` of the signal. Which effect does this have on the quantization error?

### Mid-Rise Chacteristic Curve

The in-/output relation of the mid-rise quantizer is given as

\begin{equation}
x_Q[k] = Q \cdot \Big( \underbrace{\left\lfloor\frac{ x[k] }{Q}\right\rfloor}_{index}  + \frac{1}{2} \Big)
\end{equation}

where $\lfloor \cdot \rfloor$ denotes the floor function. The quantization characteristic of a linear uniform mid-rise quantizer is illustrated in the following

The term mid-rise copes for the fact that $x[k] = 0$ is not mapped to zero. Small positive/negative values around zero are rather mapped to $\pm \frac{Q}{2}$.

#### Example

The example from above is now evaluated for the mid-rise characteristic

```python
A = 1.2  # amplitude of signal
Q = 1/10  # quantization stepsize
N = 2000  # number of samples

def uniform_midrise_quantizer(x, Q):
    # limiter
    x = np.copy(x)
    idx = np.where(np.abs(x) >= 1)
    x[idx] = np.sign(x[idx])
    # linear uniform quantization
    xQ = Q * (np.floor(x/Q) + .5)

    return xQ

# generate signal
x = A * np.sin(2*np.pi/N * np.arange(N))
# quantize signal
xQ = uniform_midrise_quantizer(x, Q)
# plot signals
plot_signals(x, xQ)
```

**Exercise**

* What are the differences between the mid-tread and the mid-rise characteristic curves for the given example?

**Copyright**

This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2017*.

max_stars_repo_name: swchao/digitalSignalProcessingLecture