1
2import algopy
3import numpy
4
5
6class Model:
7
8    # Function evaluations
9    def eval_f(self, x):
10        return (x[0]-10)**2 + 5*(x[1]-12)**2 + x[2]**4 + 3*(x[3]-11)**2 + 10*x[4]**6
11        + 7*x[5]**2 + x[6]**4 - 4*x[5]*x[6] - 10*x[5] - 8*x[6]
12
13    def eval_g(self, x):
14        out = algopy.zeros(3, dtype=x)
15        out[0] = 2*x[0]**2 + 3*x[1]**4 + x[2] + 4*x[3]**2 + 5*x[4]
16        out[1] = 7*x[0] + 3*x[1] + 10*x[2]**2 + x[3] - x[4]
17        out[2] = 23*x[0] + x[1]**2 + 6*x[5]**2 - 8*x[6]  -4*x[0]**2 - x[1]**2 + 3*x[0]*x[1] -2*x[2]**2 - 5*x[5]+11*x[6]
18        return out
19
20    def eval_Lagrangian(self, lam,x):
21        return self.eval_f(x) + algopy.dot(lam, self.eval_g(x))
22
23
24    # Forward Mode Derivative Evaluations
25    def eval_grad_f_forward(self, x):
26        x = algopy.UTPM.init_jacobian(x)
27        return algopy.UTPM.extract_jacobian(self.eval_f(x))
28
29    def eval_jac_g_forward(self, x):
30        x = algopy.UTPM.init_jacobian(x)
31        return algopy.UTPM.extract_jacobian(self.eval_g(x))
32
33    def eval_jac_vec_g_forward(self, x, v):
34        x = algopy.UTPM.init_jac_vec(x, v)
35        return algopy.UTPM.extract_jac_vec(self.eval_g(x))
36
37    def eval_grad_Lagrangian_forward(self, lam, x):
38        return self.eval_grad_f_forward(x) + algopy.dot(lam, self.eval_jac_g_forward(x))
39
40    def eval_hess_Lagrangian_forward(self, lam, x):
41        x = algopy.UTPM.init_hessian(x)
42        return algopy.UTPM.extract_hessian(x.size, self.eval_Lagrangian(lam, x))
43
44    def eval_vec_hess_g_forward(self, w, x):
45        x = algopy.UTPM.init_hessian(x)
46        tmp = algopy.dot(w, self.eval_g(x))
47        return algopy.UTPM.extract_hessian(x.size, tmp)
48
49    # Reverse Mode Derivative Evaluations
50    def trace_eval_f(self, x):
51        cg = algopy.CGraph()
52        x = algopy.Function(x)
53        y = self.eval_f(x)
54        cg.trace_off()
55        cg.independentFunctionList = [x]
56        cg.dependentFunctionList = [y]
57        self.cg = cg
58
59    def trace_eval_g(self, x):
60        cg2 = algopy.CGraph()
61        x = algopy.Function(x)
62        y = self.eval_g(x)
63        cg2.trace_off()
64        cg2.independentFunctionList = [x]
65        cg2.dependentFunctionList = [y]
66        self.cg2 = cg2
67
68    def eval_grad_f_reverse(self, x):
69        return self.cg.gradient(x)
70
71    def eval_jac_g_reverse(self, x):
72        return self.cg2.jacobian(x)
73
74    def eval_hess_f_reverse(self, x):
75        return self.cg.hessian(x)
76
77    def eval_hess_vec_f_reverse(self, x, v):
78        return self.cg.hess_vec(x,v)
79
80    def eval_vec_hess_g_reverse(self, w, x):
81        return self.cg2.vec_hess(w, x)
82
83    def eval_grad_Lagrangian_reverse(self, lam, x):
84        return self.cg.gradient(x) + self.cg2.vec_jac(lam, x)
85
86    def eval_hess_Lagrangian_reverse(self, lam, x):
87        return self.cg.hessian(x) + self.cg2.vec_hess(lam, x)
88
89
90lam = numpy.array([1,1,1],dtype=float)
91x = numpy.array([1,2,3,4,0,1,1],dtype=float)
92v = numpy.array([1,1,1,1,1,1,1],dtype=float)
93lagra = numpy.array([1,2,0],dtype=float)
94V = numpy.eye(7)
95
96m = Model()
97
98print('normal function evaluation')
99m.eval_f(x)
100m.eval_g(x)
101
102print('Forward Mode')
103
104grad_f_forward          = m.eval_grad_f_forward(x)
105jac_g_forward           = m.eval_jac_g_forward(x)
106jac_vec_g_forward       = m.eval_jac_vec_g_forward(x,[1,0,0,0,0,0,0])
107grad_Lagrangian_forward = m.eval_grad_Lagrangian_forward(lam, x)
108hess_Lagrangian_forward = m.eval_hess_Lagrangian_forward(lam, x)
109vec_hess_g_forward      = m.eval_vec_hess_g_forward(lagra, x)
110
111print(grad_f_forward)
112print(jac_g_forward)
113print(jac_vec_g_forward)
114print(grad_Lagrangian_forward)
115print(hess_Lagrangian_forward)
116print(vec_hess_g_forward)
117
118print('Reverse Mode')
119m.trace_eval_f(x)
120m.trace_eval_g(x)
121grad_f_reverse           = m.eval_grad_f_reverse(x)
122jac_g_reverse            = m.eval_jac_g_reverse(x)
123hess_f_reverse           = m.eval_hess_f_reverse(x)
124hess_vec_f_reverse       = m.eval_hess_vec_f_reverse(x,v)
125vec_hess_g_reverse       = m.eval_vec_hess_g_reverse(lagra, x)
126grad_Lagrangian_reverse  = m.eval_grad_Lagrangian_reverse(lam, x)
127hess_Lagrangian_reverse  = m.eval_hess_Lagrangian_reverse(lam, x)
128
129print(grad_f_reverse)
130print(jac_g_reverse)
131print(hess_f_reverse)
132print(hess_vec_f_reverse)
133print(vec_hess_g_reverse)
134print(grad_Lagrangian_reverse)
135
136from numpy.testing import assert_almost_equal
137assert_almost_equal(grad_f_forward, grad_f_reverse)
138assert_almost_equal(jac_g_forward, jac_g_reverse)
139assert_almost_equal(grad_Lagrangian_forward, grad_Lagrangian_reverse)
140assert_almost_equal(hess_Lagrangian_forward, hess_Lagrangian_reverse)
141assert_almost_equal(vec_hess_g_forward, vec_hess_g_reverse)