broil/machine_replacement_experiment_stress_test_numrewards.py at main · dsbrown1331/broil · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
import bayesian_irl
import mdp_worlds
import utils
import mdp
import numpy as np
import scipy
import random
import generate_efficient_frontier

def generate_reward_sample(n_states):
    #costs for no-op are -N(0,10-4) except last state that is -N(100,800)
    r_noop = []
    locs = 1/2
    scales = [i*10 for i in range(n_states)]
    for i in range(n_states):
        r_noop.append(-np.random.gamma(locs, scales[i], 1)[0])
    r_noop = np.array(r_noop)
    #print(r_noop)

    #costs for repair are -N(130,1) for all but last state where it is -N(130,20)
    r_repair = -100 + -1 * np.random.randn(n_states)
    #print(r_repair)
    return np.concatenate((r_noop, r_repair))

# def generate_reward_sample():
#     #costs for no-op are -N(0,10-4) except last state that is -N(100,800)
#     r_noop = []
#     shapes = [1,1,1,3]
#     scales = [20,30,40,50]
#     for i in range(num_states):
#         r_noop.append(-np.random.gamma(shapes[i], scales[i], 1)[0])
#     r_noop = np.array(r_noop)
#     #print(r_noop)

#     #costs for repair are -N(130,1) for all but last state where it is -N(130,20)
#     r_repair = -100 + -1 * np.random.randn(4)
#     #print(r_repair)
#     return np.concatenate((r_noop, r_repair))

# def generate_reward_mean():
#     #costs for no-op are -N(0,10-4) except last state that is -N(100,800)
#     r_noop = []
#     loc = 1/4
#     scales = [20, 40,80,300]
#     for i in range(num_states):
#         r_noop.append(-loc * scales[i])
#     r_noop = np.array(r_noop)
#     #print(r_noop)

#     #costs for repair are -N(130,1) for all but last state where it is -N(130,20)
#     r_repair = -100 * np.ones(4)
#     #print(r_repair)
#     return np.concatenate((r_noop, r_repair))

def generate_posterior_samples(num_samples, n_states):
    #print("samples")
    all_samples = []
    for i in range(num_samples):
        r_sample = generate_reward_sample(n_states)
        all_samples.append(r_sample)
        # r_string = ""
        # for r in r_sample:
        #     r_string += "{:.1f}\t".format(r)
        # print(r_string)

    #print("mean of posterior from samples")
    #print(np.mean(all_samples, axis=0))
    #print(generate_reward_mean(3))

    posterior = np.array(all_samples)

    return posterior.transpose()  #each column is a reward sample


if __name__=="__main__":
    seed = 1234
    np.random.seed(seed)
    scipy.random.seed(seed)
    random.seed(seed)

    num_states = 100
    num_reps = 20
    num_samples_list = [10,50,100,200, 500,1000,2000,5000]
    gamma = 0.95
    alpha = 0.95
    lamda = 0.5


    run_times = np.zeros((num_reps, len(num_samples_list)))
    for i,num_samples in enumerate(num_samples_list):
        print("num rewards", num_samples)
        for rep in range(num_reps):
            print(rep)
            posterior = generate_posterior_samples(num_samples, num_states)


            r_sa = np.mean(posterior, axis=1)
            #print("rsa", r_sa)
            init_distribution = np.ones(num_states)/num_states  #uniform distribution
            mdp_env = mdp.MachineReplacementMDP(num_states, r_sa, gamma, init_distribution)
            #run CVaR optimization, just the robust version since we don't have demos
            u_expert = np.zeros(mdp_env.num_actions * mdp_env.num_states)

            # print("solving for CVaR optimal policy")
            posterior_probs = np.ones(num_samples) / num_samples  #uniform dist since samples from MCMC
            import time
            t = time.time()
            cvar_opt_usa, cvar, exp_ret = mdp.solve_max_cvar_policy(mdp_env, u_expert, posterior, posterior_probs, alpha, False, lamda)
            run_times[rep,i]  = time.time() - t


print(run_times)
print(np.mean(run_times, axis=0))
print(np.std(run_times, axis=0))
import os
if not os.path.exists('./results/stress_test/'):
    os.makedirs('./results/stress_test/')
np.savetxt("./results/stress_test/machine_replace_rewards.csv", run_times, delimiter=",")