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 //#include <thread> // for sleep_for
16 
17 
18 //KRD modified this starting from shubert.cpp
main(int argc,char ** argv)19 int main(int argc, char** argv)
20 {
21 
22   std::ifstream fin(argv[1]);
23   if (!fin) {
24     std::cerr << "\nError: failure opening " << argv[1] << std::endl;
25     exit(-1);
26   }
27   size_t i, j, k, num_vars, num_fns, num_deriv_vars;
28   std::string vars_text, fns_text, dvv_text;
29 
30   // Get the parameter std::vector and ignore the labels
31   fin >> num_vars >> vars_text;
32   std::vector<double> x(num_vars);
33   for (i=0; i<num_vars; i++) {
34     fin >> x[i];
35     fin.ignore(256, '\n');
36   }
37 
38   // Get the ASV std::vector and ignore the labels
39   fin >> num_fns >> fns_text;
40   if(num_fns!=1) {
41     std::cerr << "\nError: Smoothed Herbie doesn't support constraints but you said #constraints=" << num_fns-1 << std::endl;
42     exit(-1);
43   }
44   std::vector<int> ASV(num_fns);
45   fin >> ASV[0];
46   fin.ignore(256, '\n');
47 
48   // Get the DVV std::vector and ignore the labels
49   fin >> num_deriv_vars >> dvv_text;
50   std::vector<int> DVV(num_deriv_vars);
51   for (i=0; i<num_deriv_vars; i++) {
52     fin >> DVV[i];
53     fin.ignore(256, '\n');
54   }
55 
56   //srand ( (unsigned int) (time(NULL)/x[0]) );
57   //std::this_thread::sleep_for
58   //  (std::chrono::seconds((int)(3.0*((double)rand()/RAND_MAX))));
59 
60   //std::this_thread::sleep_for(std::chrono::seconds(5));
61 
62   // Compute the results and output them directly to argv[2] (the NO_FILTER
63   // option is used).  Response tags are now optional; output them for ease
64   // of results readability.
65   std::ofstream fout(argv[2]);
66   if (!fout) {
67     std::cerr << "\nError: failure creating " << argv[2] << std::endl;
68     exit(-1);
69   }
70   fout.precision(15); // 16 total digits
71   fout.setf(std::ios::scientific);
72   fout.setf(std::ios::right);
73 
74   // **** f:
75   // f(\underline{x})=-\prod_{i=1}^num_vars w(x_i) //latex notation
76   // w(x_i) = \exp(-(x_i-1)^2)+\exp(-0.8*(x_i+1)^2) //w for smoothed herbie
77   std::vector<double> w(num_vars);
78   if (ASV[0] >= 1)
79     for (i=0; i<num_vars; i++) {
80       double dtemp1=x[i]-1.0;
81       double dtemp2=x[i]+1.0;
82       w[i]=std::exp(-dtemp1*dtemp1)+std::exp(-0.8*dtemp2*dtemp2);
83     }
84 
85   if(ASV[0] & 1) {
86     double val = 1.0;
87     for (i=0; i<num_vars; i++)
88       val*=w[i];
89     fout << "                     " << -val << " f\n";
90   }
91 
92   // **** df/dx:
93   // df/dx_i=-(\prod_{j=1}^{i-1} w(x_j)) * dw(x_i)/dx_i * (\prod_{j=i+1}^{num_vars} w(x_j))
94   // dw/dx_i= -2*(x_i-1)*\exp(-(x_i-1)^2)-0.8*2*(x_i+1)*\exp(-0.8*(x_i+1)^2) //for smoothed herbie
95   std::vector<double> d1w(num_vars);
96   if (ASV[0] >= 2) {
97     for (i=0; i<num_deriv_vars; i++) {
98       int ii=DVV[i]-1;
99       double dtemp1=x[ii]-1.0;
100       double dtemp2=x[ii]+1.0;
101       d1w[ii]=-2.0*dtemp1*std::exp(-dtemp1*dtemp1)-1.6*dtemp2*std::exp(-0.8*dtemp2*dtemp2);
102     }
103   }
104   double d1;
105   if (ASV[0] & 2) {
106     fout << "[ ";
107     for (i=0; i<num_deriv_vars; i++) {
108       int ii=DVV[i]-1;
109       d1=d1w[ii];
110       for (j=0; j<ii; ++j)
111 	d1*=w[j];
112       for (j=ii+1; j<num_vars; ++j)
113 	d1*=w[j];
114       fout << -d1 << ' ';
115     }
116     fout << "]\n";
117   }
118 
119   // **** d^2f/dx^2:
120   // if i<k   d^2f/(dx_i*dx_k)=-(\prod_{j=1)^{i-1} w(x_j))   * dw(x_i)/dx_i ...
121   //                           *(\prod_{j=i+1)^{k-1} w(x_j)) * dw(x_k)/dx_k ...
122   //                           *(\prod_{j=k+1)^{num_vars} w(x_j))
123   // if i==k   d^2f/dx_i^2 =(\prod_{j=1)^{i-1} w(x_j)) * d^2w(x_i)/dx_i^2 * (\prod_{j=i+1)^{num_vars} w(x_j))
124   // d^2w(x_i)/dx_i^2 = -2*\exp(-(x_i-1)^2)+4*(x_i-1)^2*\exp(-(x_i-1)^2) ...
125   //                    -1.6*\exp(-0.8*(x_i+1)^2)+2.56*(x_i+1)^2*\exp(-0.8*(x_i+1)^2) //for smoothed herbie
126   if (ASV[0] & 4) {
127     fout << "[[ ";
128     double d2;
129     for (i=0; i<num_deriv_vars; ++i) {
130       for (k=0; k<num_deriv_vars; ++k) {
131 	int ii=DVV[i]-1;
132 	int kk=DVV[k]-1;
133 	if (ii==kk) {
134 	  double dtemp1=x[ii]-1.0; dtemp1*=dtemp1;
135 	  double dtemp2=x[ii]+1.0; dtemp2*=dtemp2;
136 	  d2=(-2.0+4.0*dtemp1)*std::exp(-dtemp1)
137 	    +(-1.6+2.56*dtemp2)*std::exp(-0.8*dtemp2);
138 	  for (j=0; j<ii; ++j)
139 	    d2*=w[j];
140 	  for (j=ii+1; j<num_vars; ++j)
141 	    d2*=w[j];
142 	}
143 	else {
144 	  d2=d1w[ii]*d1w[kk];
145 	  if(kk<ii) {
146 	    j=ii;
147 	    ii=kk;
148 	    kk=j;
149 	  }
150 	  for (j=0; j<ii; ++j)
151 	    d2*=w[j];
152 	  for (j=ii+1; j<kk; ++j)
153 	    d2*=w[j];
154 	  for (j=kk+1; j<num_vars; ++j)
155 	    d2*=w[j];
156 	}
157 	//	for (j=0; j<num_vars; ++j)
158 	//	  if((j!=ii)&&(j!=kk))
159 	//	    d2*=w[j];
160 	fout << -d2 << ' ';
161       }
162     }
163     fout << "]]\n";
164   }
165 
166   fout.flush();
167   fout.close();
168   return 0;
169 }
170