GMCA_optimisation_problem/callback_routine.py at main · vj2Sheffield/GMCA_optimisation_problem · 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
119
120
121
122
123
124
125
126
127
128
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 22 13:58:48 2025

@author: Victoria Johnson
"""

import numpy as np
from scipy.optimize import fmin_bfgs
from scipy import optimize
import matplotlib.pyplot as plt
import pandas as pd
from GMCA_admin_functions import get_clusters
from GMCA_functions_beta_gamma import x_to_u


def MDA_new(y, u, cluster_dict):
    A = cluster_dict["A"]
    B = cluster_dict["B"]
    intercept = cluster_dict["intercept"]
    JPC = cluster_dict["baseline_JPC"]

    node_outputs = [intercept[0] + B[0][0]*u[0] + B[0][1]*JPC + B[0][2]*u[1] - y[0],
    intercept[1] + B[1][0]*u[0] + B[1][1]*JPC + A[1][0]*y[0] - y[1],
    intercept[2] + B[2][0]*u[0] + A[2][1]*y[1] + A[2][4]*y[4] + A[2][5]*y[5] + A[2][7]*y[7] - y[2],
    intercept[3] + B[3][0]*u[0] + B[3][1]*JPC + A[3][0]*y[0] + A[3][1]*y[1] + A[3][6]*y[6] - y[3],
    intercept[4] + B[4][0]*u[0] + A[4][1]*y[1] + A[4][3]*y[3] + A[4][5]*y[5] + A[4][6]*y[6] + A[4][7]*y[7] - y[4],
    intercept[5] + B[5][1]*JPC + A[5][0]*y[0] + A[5][1]*y[1] + A[5][3]*y[3] +A[5][7]*y[7] - y[5],
    intercept[6] + B[6][1]*JPC + A[6][1]*y[1] - y[6],
    intercept[7] + B[7][0]*u[0] + B[7][1]*JPC + A[7][1]*y[1] + A[7][3]*y[3] - y[7]]

    return node_outputs

Nfeval = 1
fx_vals = []
x_vals = []
def callback(x, fx):
    global Nfeval, fx_vals, x_vals
    fx_vals.append(fx)
    x_vals.append(x)

    Nfeval += 1


def run_MDA_new(u, cluster_dict):
# =============================================================================
#     Inputs:
#       u: a 1 x 20 array - the inputs to the MDA
#       cluster_dict: A dictionary containing:
#           A_matrix: All the A matrices for the model (list of np.ndarray)
#           B_matrix: All the B matrices for the model (list of np.ndarray)
#           int_matrix: All the intercepts for the model (list of np.ndarray)
#           baseline_JPC: All the baseline JPC values for all LAs (np.ndarray)
#           LA_names: The LA names (list of str)
#
#     Outputs:
#       sol_array: An array of size 10 x 8 containing the node values for each LA
# =============================================================================
    u = np.array(u).reshape((10, 2))
    sol_array = np.zeros([10, 8])
    for m in range(1,2):
        u_in = u[m, :]
        cluster_dict_m = {
            "A": cluster_dict["A_matrix"][m],
            "B": cluster_dict["B_matrix"][m],
            "intercept": cluster_dict["intercept_matrix"][m],
            "baseline_JPC": cluster_dict["baseline_JPC"][m]}
        sol = optimize.root(MDA_new, [0, 0, 0, 0, 0, 0, 0, 0], args = (u_in, cluster_dict_m), method='broyden1', callback = callback)
        sol_array[m, :] = sol.x
    print(sol)

    return sol_array

# Load in data
pops = np.array(pd.read_csv('Model_Data/populations.csv', header = None))[:, 1]
beta = np.array(pd.read_csv('Model_Data/beta.csv', header = None))
gamma = np.array(pd.read_csv('Model_Data/gamma.csv', header = None))
LA_params = pd.read_csv('Model_Data/LA_params.csv')
model_params = pd.read_csv('Model_Data/model_params.csv')
current_funding = np.array(LA_params["current_funding"])
baseline_JPC = np.array(LA_params["baseline_JPC"])
mu = np.array(model_params["mu"])
sigma = np.array(model_params["sigma"])

# Private settings
#G1 = 30.382003
G1 = sum(current_funding)
N = 22
lb = np.zeros(N)
lb[1:20:2]=-current_funding
lb[21] = -sum(current_funding)
ub = G1*np.ones(N)

cluster_A, cluster_B, cluster_int = get_clusters()

cluster_dict = {
  "A_matrix": cluster_A,
  "B_matrix": cluster_B,
  "intercept_matrix": cluster_int,
  "baseline_JPC": baseline_JPC,
  "LA_names": ['BOL', 'BUR', 'MAN', 'OLD', 'ROC', 'SAL', 'STO', 'TAM', 'TRA', 'WIG']
}


X = np.array(pd.read_csv('Results/rep15_X.csv'))
u = x_to_u(X[0,:], beta, gamma, mu, sigma)
run_MDA_new(u, cluster_dict)
a=np.concatenate( fx_vals, axis=0 )
a=a.reshape((len(fx_vals),8))
b=np.concatenate(x_vals, axis=0 )
b=b.reshape((len(x_vals),8))

fig, ([ax1, ax2]) = plt.subplots(2, 1, figsize =(8, 10))
ax1.plot(np.arange(0, len(fx_vals)), a)
ax1.set_title("Residuals")
ax1.set_xlabel("Solver iteration")
ax1.set_ylabel("Residuals")


ax2.plot(np.arange(0, len(x_vals)), b)
ax2.set_title("Node value")
ax2.set_xlabel("Solver iteration")
ax2.set_ylabel("Node value (z-score)")
ax1.legend(["EPR", "PVT", "FPR", "DHI", "EIQ", "MCS", "PMR", "HLE"])

plt.savefig('Plots/solver_progress_BOL_sol_0.eps', dpi = 200)
plt.show