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 <string>
13 #include <vector>
14 #include <map>
15 #include <algorithm>
16 #include <cctype>
17 
18 enum var_t { X1, X2 };
19 
20 
main(int argc,char ** argv)21 int main(int argc, char** argv)
22 {
23   // The Rosenbrock function may be solved as either a general minimization
24   // problem with Objective function = 100.*(x1-x0^2)^2 + (1-x0)^2
25   // or a least squares problem with Term1 = 10.*(x1-x0^2) and Term2 = (1-x0).
26   // See p. 95 in Practical Optimization by Gill, Murray, and Wright.
27 
28   // This application program reads and writes parameter and response data
29   // directly so no input/output filters are needed.
30   std::ifstream fin(argv[1]);
31   if (!fin) {
32     std::cerr << "\nError: failure opening " << argv[1] << std::endl;
33     exit(-1);
34   }
35   size_t i, j, num_vars, num_fns, num_deriv_vars;
36   std::string vars_text, fns_text, dvv_text;
37 
38   // define the std::string to enumeration map
39   std::map<std::string, var_t> var_t_map;
40   var_t_map["x1"] = X1;
41   var_t_map["x2"] = X2;
42 
43   // Get the parameter std::vector and ignore the labels
44   fin >> num_vars >> vars_text;
45   std::map<var_t, double> vars;
46   std::vector<var_t> labels(num_vars);
47   double var_i; std::string label_i; var_t v_i;
48   std::map<std::string, var_t>::iterator v_iter;
49   for (i=0; i<num_vars; i++) {
50     fin >> var_i >> label_i;
51     transform(label_i.begin(), label_i.end(), label_i.begin(),
52 	      (int(*)(int))tolower);
53     v_iter = var_t_map.find(label_i);
54     if (v_iter == var_t_map.end()) {
55       std::cerr << "Error: label \"" << label_i << "\" not supported in analysis "
56 	   << "driver." << std::endl;
57       exit(-1);
58     }
59     else
60       v_i = v_iter->second;
61     vars[v_i] = var_i;
62     labels[i] = v_i;
63   }
64 
65   // Get the ASV std::vector and ignore the labels
66   fin >> num_fns >> fns_text;
67   std::vector<short> ASV(num_fns);
68   for (i=0; i<num_fns; i++) {
69     fin >> ASV[i];
70     fin.ignore(256, '\n');
71   }
72 
73   // Get the DVV std::vector and ignore the labels
74   fin >> num_deriv_vars >> dvv_text;
75   std::vector<var_t> DVV(num_deriv_vars);
76   unsigned int dvv_i;
77   for (i=0; i<num_deriv_vars; i++) {
78     fin >> dvv_i;
79     fin.ignore(256, '\n');
80     DVV[i] = labels[dvv_i-1];
81   }
82 
83   if (num_vars != 2) {
84     std::cerr << "Wrong number of variables for the rosenbrock problem\n";
85     exit(-1);
86   }
87   if (num_fns < 1 || num_fns > 2) { // 1 fn -> opt, 2 fns -> least sq
88     std::cerr << "Wrong number of functions in rosenbrock problem\n";
89     exit(-1);
90   }
91 
92   // Compute and output responses
93   bool least_sq_flag = (num_fns > 1) ? true : false;
94   double x1 = vars[X1], x2 = vars[X2];
95   double f1 = x2-x1*x1;
96   double f2 = 1.-x1;
97 
98   std::ofstream fout(argv[2]);
99   if (!fout) {
100     std::cerr << "\nError: failure creating " << argv[2] << std::endl;
101     exit(-1);
102   }
103   fout.precision(15); // 16 total digits
104   fout.setf(std::ios::scientific);
105   fout.setf(std::ios::right);
106 
107   if (least_sq_flag) {
108     if (ASV[0] & 1) // **** Residual R1:
109       fout << "                     " << 10.*f1  << " f1\n";
110     if (ASV[1] & 1) // **** Residual R2:
111       fout << "                     " << f2  << " f2\n";
112 
113     if (ASV[0] & 2) { // **** dR1/dx:
114       fout << "[ ";
115       for (i=0; i<num_deriv_vars; i++)
116 	switch (DVV[i]) {
117 	case X1: fout << -20.*x1 << ' '; break;
118 	case X2: fout <<  10.    << ' '; break;
119 	}
120       fout << "]\n ";
121     }
122     if (ASV[1] & 2) { // **** dR2/dx:
123       fout << "[ ";
124       for (i=0; i<num_deriv_vars; i++)
125 	switch (DVV[i]) {
126 	case X1: fout << -1. << ' '; break;
127 	case X2: fout <<  0. << ' '; break;
128 	}
129       fout << "]\n";
130     }
131 
132     if (ASV[0] & 4) { // **** d^2R1/dx^2:
133       fout << "[[ ";
134       for (i=0; i<num_deriv_vars; i++)
135 	for (j=0; j<num_deriv_vars; j++)
136 	  if (DVV[i] == X1 && DVV[j] == X1)
137 	    fout <<  -20. << ' ';
138 	  else
139 	    fout <<    0. << ' ';
140       fout << "]]\n";
141     }
142     if (ASV[1] & 4) { // **** d^2R2/dx^2:
143       fout << "[[ ";
144       for (i=0; i<num_deriv_vars; i++)
145 	for (j=0; j<num_deriv_vars; j++)
146 	  fout << 0. << ' ';
147       fout << "]]\n";
148     }
149   }
150   else {
151     if (ASV[0] & 1) { // **** f:
152       double f = 100.*f1*f1+f2*f2;
153       bool fail = (f > 250.); // fail away from min value (less problematic)
154       //bool fail = (f < 5.); // fail near min value (more problematic)
155       if (fail) fout << "                     nan f\n";
156       else      fout << "                     " << f << " f\n";
157     }
158 
159     if (ASV[0] & 2) { // **** df/dx:
160       fout << "[ ";
161       for (i=0; i<num_deriv_vars; i++)
162 	switch (DVV[i]) {
163 	case X1: { // failures for high gradients
164 	  double df_dx1 = -400.*f1*x1 - 2.*f2;
165 	  if (df_dx1 > 800.) fout << "nan ";
166 	  else               fout << df_dx1 << ' ';
167 	  break;
168 	}
169 	case X2: {
170 	  double df_dx2 = 200.*f1;
171 	  if (df_dx2 > 200.) fout << "nan ";
172 	  else               fout << df_dx2 << ' ';
173 	  break;
174 	}
175 	}
176       fout << "]\n";
177     }
178 
179     if (ASV[0] & 4) { // **** d^2f/dx^2:
180       fout << "[[ ";
181       for (i=0; i<num_deriv_vars; i++)
182 	for (j=0; j<num_deriv_vars; j++)
183 	  if (DVV[i] == X1 && DVV[j] == X1)
184 	    fout << -400.*(x2 - 3.*x1*x1) + 2. << ' ';
185 	  else if ( (DVV[i] == X1 && DVV[j] == X2) ||
186 		    (DVV[i] == X2 && DVV[j] == X1) )
187 	    fout << -400.*x1 << ' ';
188 	  else if (DVV[i] == X2 && DVV[j] == X2)
189 	    fout <<  200. << ' ';
190       fout << "]]\n";
191     }
192   }
193 
194   fout.flush();
195   fout.close();
196   return 0;
197 }
198