1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 #include <cstdlib>
10 #include <iostream>
11 #include <fstream>
12 #include <vector>
13 #include <string>
14 #include <cmath>
15
16
main(int argc,char ** argv)17 int main(int argc, char** argv)
18 {
19
20 // The illumination example in Boyd as a general minimization problem
21 // Objective function = ...
22
23 // This application program reads and writes parameter and response data
24 // directly so that the NO_FILTER option of dakota may be used.
25
26 std::ifstream fin(argv[1]);
27 if (!fin) {
28 std::cerr << "\nError: failure opening " << argv[1] << std::endl;
29 exit(-1);
30 }
31 size_t i, j, num_vars, num_fns;
32 std::string vars_text, fns_text;
33
34 // Get the parameter std::vector and ignore the labels
35 fin >> num_vars >> vars_text;
36 std::vector<double> x(num_vars);
37 for (i=0; i<num_vars; i++) {
38 fin >> x[i];
39 fin.ignore(256, '\n');
40 }
41
42 // Get the ASV std::vector and ignore the labels
43 fin >> num_fns >> fns_text;
44 std::vector<int> ASV(num_fns);
45 for (i=0; i<num_fns; i++) {
46 fin >> ASV[i];
47 fin.ignore(256, '\n');
48 }
49
50 if (num_vars != 7) {
51 std::cerr << "Wrong number of variables for the illumination problem" << std::endl;
52 exit(-1);
53 }
54 if (num_fns != 1) {
55 std::cerr << "Wrong number of functions for the illumination problem" << std::endl;
56 exit(-1);
57 }
58
59 // compute function and gradient values
60 double A[11][7] ={
61 { 0.347392, 0.205329, 0.191987, 0.077192, 0.004561, 0.024003, 0.000000},
62 { 0.486058, 0.289069, 0.379202, 0.117711, 0.006667, 0.032256, 0.000000},
63 { 0.752511, 0.611283, 2.417907, 0.701700, 0.473047, 0.285597, 0.319187},
64 { 0.303582, 0.364620, 1.898185, 0.693173, 0.607718, 0.328582, 0.437394},
65 { 0.540946, 0.411549, 1.696545, 0.391735, 0.177832, 0.110119, 0.083817},
66 { 0.651840, 0.540687, 3.208793, 0.639020, 0.293811, 0.156842, 0.128499},
67 { 0.098008, 0.245771, 0.742564, 0.807976, 0.929739, 0.435144, 0.669797},
68 { 0.000000, 0.026963, 0.000000, 0.246606, 0.414657, 0.231777, 0.372202},
69 { 0.285597, 0.320457, 0.851227, 0.584677, 0.616436, 0.341447, 0.477329},
70 { 0.324622, 0.306394, 0.991904, 0.477744, 0.376266, 0.158288, 0.198745},
71 { 0.000000, 0.050361, 0.000000, 0.212042, 0.434397, 0.286455, 0.462731} };
72
73 // KRD found bugs in previous computation; this Hessian no longer used:
74 double harray[7][7] ={
75 { 1.929437, 1.572662, 6.294004, 1.852205, 1.222324, 0.692036, 0.768564},
76 { 1.572662, 1.354287, 5.511537, 1.787932, 1.320048, 0.724301, 0.870382},
77 { 6.294004, 5.511537, 25.064512, 7.358494, 5.133563, 2.791970, 3.257364},
78 { 1.852205, 1.787932, 7.358494, 2.883178, 2.497491, 1.321922, 1.747230},
79 { 1.222324, 1.320048, 5.133563, 2.497491, 2.457733, 1.295927, 1.816568},
80 { 0.692036, 0.724301, 2.791970, 1.321922, 1.295927, 0.694642, 0.968982},
81 { 0.768564, 0.870382, 3.257364, 1.747230, 1.816568, 0.968982, 1.385357} };
82
83
84 std::ofstream fout(argv[2]);
85 if (!fout) {
86 std::cerr << "\nError: failure creating " << argv[2] << std::endl;
87 exit(-1);
88 }
89 fout.precision(15); // 16 total digits
90 fout.setf(std::ios::scientific);
91 fout.setf(std::ios::right);
92
93 // Calculation of grad(f) = df/dx_I and hess(f) = d^2f/(dx_I dx_J):
94 //
95 // U = sum_{i=0:10}{ ( 1 - sum_{j=0:6}{ A[i][j]*x[j] } )^2 }
96 // dU/dx_I = sum_{i=0:10}{ 2.0*( 1 - sum_{j=0:6}{ A[i][j]*x[j] } )*-A[i][I] }
97 // d^2U/(dx_I dx_J) = sum_{i=0:10){ 2.0 * A[i][I] * A[i][J] }
98 //
99 // f = sqrt(U)
100 // df/dx_I = 0.5*(dU/dx_I)/sqrt(U)
101 // = 0.5*(dU/dx_I)/f
102 // d^2f/(dx_I dx_J)
103 // = -0.25 * U^(-1.5) * dU/dx_I * dU/dx_J + 0.5/sqrt(U) * d2U/(dx_I dx_J)
104 // = ( -df/dx_I * df/dx_J + 0.5 * d2U/(dx_I dx_J) ) / f
105 //
106 // so to (efficiently) compute grad(f) we precompute f
107 // and to (efficiently) compute hess(f) we precompute f and grad(f)
108
109 int Ider, Jder;
110 double grad[7];
111 for(Ider=0; Ider<num_vars; ++Ider)
112 grad[Ider] = 0.0;
113
114 // **** f: (f is required to calculate any derivative of f; perform always)
115 double U = 0.0;
116 for (i=0; i<11; i++) {
117 double dtmp = 0.0;
118 for (j=0; j<num_vars; j++)
119 dtmp += A[i][j] * x[j];
120 dtmp = 1.0 - dtmp;
121 U += dtmp*dtmp;
122
123 // precompute grad(U) unconditionally as it might be needed (cheap)
124 for(Ider=0; Ider<num_vars; ++Ider)
125 grad[Ider] -= dtmp*2.0*A[i][Ider]; //this is grad(U)
126 }
127 double fx = sqrt(U);
128 if (ASV[0] & 1)
129 fout << " " << fx << " f\n";
130
131
132 // **** df/dx:
133 if (ASV[0] & 6) { // gradient required for itself and to calcuate the Hessian
134 for(Ider=0; Ider<num_vars; ++Ider)
135 grad[Ider] *= (0.5/fx); // this is now updated to be grad(f)
136
137 if (ASV[0] & 2) { // if ASV demands gradient, print it
138 fout << "[ ";
139 for (int Ider=0; Ider<num_vars; ++Ider)
140 fout << grad[Ider] << " ";
141 fout << "]\n";
142 }
143 }
144
145 // **** the hessian ddf/dxIdxJ
146 if (ASV[0] & 4) {
147 double hess[7][7];
148 fout << "[[ ";
149 for(Ider=0; Ider<num_vars; ++Ider) {
150 for(Jder=Ider; Jder<num_vars; ++Jder) {
151 // define triangular part of hess(U)
152 hess[Ider][Jder] = 0.0;
153 for(i=0; i<11; ++i)
154 hess[Ider][Jder] += A[i][Ider]*A[i][Jder]; // this is 0.5*hess(U)
155
156 hess[Jder][Ider]= hess[Ider][Jder] // this is hess(f)
157 = (hess[Ider][Jder]- grad[Ider]*grad[Jder])/fx;
158
159 }
160 for(Jder=0; Jder<num_vars; ++Jder)
161 fout << hess[Ider][Jder] << " ";
162 fout << "\n";
163 }
164 fout <<" ]] \n";
165 }
166
167 fout.flush();
168 fout.close();
169 return 0;
170 }
171