Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
4 : /* */
5 : /* Copyright 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 :
9 : #include "NeuralNetFreeEnergyBase.h"
10 : #include "MooseEnum.h"
11 : #include <fstream>
12 :
13 : InputParameters
14 38 : NeuralNetFreeEnergyBase::validParams()
15 : {
16 38 : auto params = ADMaterial::validParams();
17 38 : params.addClassDescription(
18 : "Evaluates a fitted deep neural network to obtain a free energy and its derivatives.");
19 :
20 76 : MooseEnum fileFormatEnum("MAGPIE GENANN");
21 76 : params.template addParam<MooseEnum>(
22 : "file_format", fileFormatEnum, "Weights and biases file format");
23 :
24 76 : params.template addParam<FileName>(
25 : "file_name",
26 : "Data file containing the weights and biasses for a fully connected deep neural network");
27 76 : params.addCoupledVar("inputs", "Coupled Variables that are inputs for the neural network");
28 76 : params.template addParam<std::vector<MaterialPropertyName>>(
29 : "prop_names", {}, "list of material properties fed from the outputs of the neural network");
30 76 : params.template addParam<bool>(
31 76 : "debug", false, "Tabulate the NN to a file for debugging purposes");
32 38 : return params;
33 38 : }
34 :
35 30 : NeuralNetFreeEnergyBase::NeuralNetFreeEnergyBase(const InputParameters & parameters)
36 : : ADMaterial(parameters),
37 60 : _file_format(getParam<MooseEnum>("file_format").template getEnum<FileFormat>()),
38 30 : _file_name(getParam<FileName>("file_name")),
39 90 : _output_name(getParam<std::vector<MaterialPropertyName>>("prop_names")),
40 30 : _n_output(_output_name.size()),
41 30 : _output(_n_output),
42 30 : _n_input(coupledComponents("inputs")),
43 30 : _input(_n_input),
44 60 : _d_output(_n_input * _n_output)
45 : {
46 : // open the NN data file
47 30 : std::ifstream ifile;
48 30 : ifile.open(_file_name);
49 30 : if (!ifile)
50 0 : paramError("file_name", "Unable to open file");
51 :
52 : // call the reader for the requested format
53 30 : switch (_file_format)
54 : {
55 15 : case FileFormat::MAGPIE:
56 15 : loadMagpieNet(ifile);
57 : break;
58 :
59 15 : case FileFormat::GENANN:
60 15 : loadGenANN(ifile);
61 : break;
62 :
63 0 : default:
64 0 : paramError("file_format", "Unknown file format");
65 : }
66 :
67 : // close parameter file
68 30 : ifile.close();
69 :
70 : // validate network properties
71 30 : if (_weight.size() != _bias.size())
72 0 : paramError("file_name", "Inconsistent layer numbers in datafile");
73 30 : if (_z[_n_layer - 1].size() != _n_output)
74 0 : paramError("prop_names",
75 : "Number of supplied property names must match the number of output nodes ",
76 : _z[_n_layer - 1].size(),
77 : " of the neural net");
78 30 : if (_n_input != _activation[0].size())
79 0 : paramError("inputs",
80 : "Number of supplied variables must match the number of input nodes ",
81 : _activation[0].size(),
82 : " of the neural net");
83 :
84 : // initialize storage for derivatives
85 30 : _prod.resize(_n_layer);
86 30 : _diff.resize(_n_layer);
87 30 : _diff[0] = _weight[0]; // _prod[0] is not used
88 105 : for (std::size_t i = 1; i < _n_layer; ++i)
89 : {
90 75 : _prod[i] = _weight[i];
91 75 : _diff[i] = DenseMatrix<ADReal>(_z[i].size(), _n_input);
92 : }
93 :
94 : // get coupled variables for the inputs
95 90 : for (std::size_t j = 0; j < _n_input; ++j)
96 60 : _input[j] = &adCoupledValue("inputs", j);
97 :
98 : // create material properties
99 : std::size_t k = 0;
100 60 : for (std::size_t i = 0; i < _n_output; ++i)
101 : {
102 30 : _output[i] = &declareADProperty<Real>(_output_name[i]);
103 90 : for (std::size_t j = 0; j < _n_input; ++j)
104 60 : _d_output[k++] = &declareADProperty<Real>(
105 120 : derivativePropertyNameFirst(_output_name[i], this->getVar("inputs", j)->name()));
106 : }
107 30 : }
108 :
109 : void
110 30 : NeuralNetFreeEnergyBase::initialSetup()
111 : {
112 60 : if (getParam<bool>("debug"))
113 0 : debugDump();
114 30 : }
115 :
116 : void
117 0 : NeuralNetFreeEnergyBase::debugDump()
118 : {
119 0 : std::ofstream ofile;
120 0 : ofile.open(_name + "_debug.dat");
121 0 : for (Real T = 0.0; T <= 1.0; T += 0.01)
122 0 : for (Real c = 0.0; c <= 1.0; c += 0.01)
123 : {
124 0 : ofile << T << ' ' << c;
125 :
126 : // set inputs
127 0 : _activation[0](0) = T;
128 0 : _activation[0](1) = c;
129 :
130 : // evaluate network
131 0 : evaluate();
132 :
133 : // output values
134 0 : for (std::size_t i = 0; i < _n_output; ++i)
135 0 : ofile << ' ' << _z[_n_layer - 1](i);
136 :
137 0 : ofile << '\n';
138 : }
139 0 : }
140 :
141 : void
142 15 : NeuralNetFreeEnergyBase::loadGenANN(std::ifstream & ifile)
143 : {
144 : std::size_t x, y;
145 :
146 : // read layer layout
147 : unsigned int n_inputs, n_hidden_layers, n_hidden, n_outputs;
148 30 : if (!(ifile >> n_inputs >> n_hidden_layers >> n_hidden >> n_outputs))
149 0 : mooseError("Error reading genann file header from file ", _file_name);
150 15 : _n_layer = n_hidden_layers + 1;
151 15 : _bias.resize(_n_layer);
152 15 : _weight.resize(_n_layer);
153 15 : _z.resize(_n_layer);
154 15 : _activation.resize(_n_layer);
155 15 : _d_activation.resize(_n_layer);
156 :
157 : // read weights and biases
158 45 : for (std::size_t i = 0; i < _n_layer; ++i)
159 : {
160 : // negative biases come first
161 30 : x = i < _n_layer - 1 ? n_hidden : n_outputs;
162 :
163 30 : _z[i] = DenseVector<ADReal>(x);
164 30 : _bias[i] = DenseVector<Real>(x);
165 420 : for (std::size_t j = 0; j < x; ++j)
166 390 : if (!(ifile >> _bias[i](j)))
167 0 : mooseError("Error reading biases from file ", _file_name);
168 : else
169 390 : _bias[i](j) *= -1.0;
170 :
171 : // next come the weights
172 30 : y = i > 0 ? n_hidden : n_inputs;
173 : // initialize weight matrix
174 30 : _weight[i] = DenseMatrix<Real>(x, y);
175 :
176 : // initialize computation buffers (including input and output)
177 30 : _activation[i] = DenseVector<ADReal>(y);
178 30 : _d_activation[i] = DenseVector<ADReal>(y);
179 30 : _z[i] = DenseVector<ADReal>(x);
180 :
181 435 : for (std::size_t k = 0; k < y; ++k)
182 1530 : for (std::size_t j = 0; j < x; ++j)
183 1125 : if (!(ifile >> _weight[i](j, k)))
184 0 : mooseError("Error reading weights from file ", _file_name);
185 : }
186 15 : }
187 :
188 : void
189 15 : NeuralNetFreeEnergyBase::loadMagpieNet(std::ifstream & ifile)
190 : {
191 : std::size_t x, y;
192 :
193 : // read weights (first the number of layers)
194 15 : ifile >> _n_layer;
195 15 : _weight.resize(_n_layer);
196 15 : _z.resize(_n_layer);
197 15 : _activation.resize(_n_layer);
198 15 : _d_activation.resize(_n_layer);
199 90 : for (std::size_t i = 0; i < _n_layer; ++i)
200 : {
201 75 : if (!(ifile >> y >> x))
202 0 : mooseError("Error reading file ", _file_name);
203 :
204 : // initialize weight matrix
205 75 : _weight[i] = DenseMatrix<Real>(x, y);
206 :
207 : // initialize computation buffers (including input and output)
208 75 : _activation[i] = DenseVector<ADReal>(y);
209 75 : _d_activation[i] = DenseVector<ADReal>(y);
210 75 : _z[i] = DenseVector<ADReal>(x);
211 :
212 1305 : for (std::size_t k = 0; k < y; ++k)
213 20130 : for (std::size_t j = 0; j < x; ++j)
214 18900 : if (!(ifile >> _weight[i](j, k)))
215 0 : mooseError("Error reading weights from file ", _file_name);
216 : }
217 :
218 : // read biases
219 15 : _bias.resize(_n_layer);
220 90 : for (std::size_t i = 0; i < _n_layer; ++i)
221 : {
222 75 : if (!(ifile >> x))
223 0 : mooseError("Error reading file ", _file_name);
224 :
225 75 : _z[i] = DenseVector<ADReal>(x);
226 75 : _bias[i] = DenseVector<Real>(x);
227 :
228 1290 : for (std::size_t j = 0; j < x; ++j)
229 1215 : if (!(ifile >> _bias[i](j)))
230 0 : mooseError("Error reading biases from file ", _file_name);
231 : }
232 15 : }
233 :
234 : void
235 4800 : NeuralNetFreeEnergyBase::computeQpProperties()
236 : {
237 : // set input nodes
238 14400 : for (std::size_t j = 0; j < _n_input; ++j)
239 9600 : _activation[0](j) = (*_input[j])[_qp];
240 :
241 : // evaluate network
242 4800 : evaluate();
243 :
244 : // copy back output values
245 : std::size_t k = 0;
246 9600 : for (std::size_t i = 0; i < _n_output; ++i)
247 : {
248 4800 : (*_output[i])[_qp] = _z[_n_layer - 1](i);
249 14400 : for (std::size_t j = 0; j < _n_input; ++j)
250 9600 : (*_d_output[k++])[_qp] = _diff[_n_layer - 1](i, j);
251 : }
252 4800 : }
253 :
254 : void
255 12000 : NeuralNetFreeEnergyBase::multiply(DenseMatrix<ADReal> & M1,
256 : const DenseMatrix<ADReal> & M2,
257 : const DenseMatrix<ADReal> & M3)
258 : {
259 : // Assertions to make sure we have been
260 : // passed matrices of the correct dimension.
261 : libmesh_assert_equal_to(M1.m(), M2.m());
262 : libmesh_assert_equal_to(M1.n(), M3.n());
263 : libmesh_assert_equal_to(M2.n(), M3.m());
264 :
265 : const unsigned int m_s = M2.m();
266 : const unsigned int p_s = M2.n();
267 : const unsigned int n_s = M1.n();
268 :
269 36000 : for (unsigned int j = 0; j < n_s; j++)
270 321600 : for (unsigned int i = 0; i < m_s; i++)
271 : {
272 297600 : M1.el(i, j) = 0.0;
273 6273600 : for (unsigned int k = 0; k < p_s; k++)
274 5976000 : M1.el(i, j) += M2.el(i, k) * M3.el(k, j);
275 : }
276 12000 : }
277 :
278 : void
279 4800 : NeuralNetFreeEnergyBase::evaluate()
280 : {
281 4800 : _layer = 0;
282 : while (true)
283 : {
284 : // apply weights and biases
285 16800 : _weight[_layer].vector_mult(_z[_layer], _activation[_layer]);
286 16800 : _z[_layer] += _bias[_layer];
287 :
288 : // derivatives
289 16800 : if (_layer > 0)
290 : {
291 : // prepare product of weights and activation function derivative (previous layer)
292 160800 : for (std::size_t j = 0; j < _weight[_layer].m(); ++j)
293 3136800 : for (std::size_t k = 0; k < _weight[_layer].n(); ++k)
294 5976000 : _prod[_layer](j, k) = _weight[_layer](j, k) * _d_activation[_layer](k);
295 :
296 : // multiply progressive Jacobian
297 12000 : multiply(_diff[_layer], _prod[_layer], _diff[_layer - 1]);
298 : }
299 :
300 : // bail to avoid applying activation function to the output
301 16800 : if (_layer + 1 == _n_layer)
302 : break;
303 :
304 : // apply activation function
305 12000 : applyLayerActivation();
306 :
307 : // next layer
308 12000 : ++_layer;
309 12000 : }
310 4800 : }
|
