Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : MIT License
3 :
4 : Copyright (c) 2019 Wei Chen and Andrew Ferguson
5 :
6 : Permission is hereby granted, free of charge, to any person obtaining a copy
7 : of this software and associated documentation files (the "Software"), to deal
8 : in the Software without restriction, including without limitation the rights
9 : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10 : copies of the Software, and to permit persons to whom the Software is
11 : furnished to do so, subject to the following conditions:
12 :
13 : The above copyright notice and this permission notice shall be included in all
14 : copies or substantial portions of the Software.
15 :
16 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19 : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21 : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22 : SOFTWARE.
23 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
24 : #include "function/Function.h"
25 : #include "core/ActionRegister.h"
26 : #include "cassert"
27 :
28 : #include <string>
29 : #include <cmath>
30 : #include <iostream>
31 :
32 : using namespace std;
33 :
34 :
35 : namespace PLMD {
36 : namespace function {
37 : namespace annfunc {
38 :
39 : //+PLUMEDOC ANNMOD_Function ANN
40 : /*
41 : Calculates the ANN-function.
42 :
43 : This module implements ANN class, which is a subclass of Function class.
44 : ANN class takes multi-dimensional arrays as inputs for a fully-connected feedforward neural network
45 : with specified neural network weights and generates corresponding outputs.
46 : The ANN outputs can be used as collective variables, inputs for other collective variables,
47 : or inputs for data analysis tools.
48 :
49 : ## Examples
50 :
51 : Assume we have an ANN with numbers of nodes being [2, 3, 1], and weights connecting layer 0 and 1 are
52 :
53 : [[1,2], [3,4], [5,6]]
54 :
55 : weights connecting layer 1 and 2 are
56 :
57 : [[7,8,9]]
58 :
59 : Bias for layer 1 and 2 are [10, 11, 12] and [13], respectively.
60 :
61 : All activation functions are Tanh.
62 :
63 : Then if input variables are l_0_out_0, l_0_out_1, the corresponding ANN function object can be defined using
64 : following plumed script:
65 :
66 : ```plumed
67 : ann: ANN ...
68 : ARG=l_0_out_0,l_0_out_1
69 : NUM_LAYERS=3
70 : NUM_NODES=2,3,1
71 : ACTIVATIONS=Tanh,Tanh
72 : WEIGHTS0=1,2,3,4,5,6
73 : WEIGHTS1=7,8,9
74 : BIASES0=10,11,12
75 : BIASES1=13
76 : ...
77 : ```
78 :
79 : To access its components, we use "ann.node-0", "ann.node-1", ..., which represents the components of neural network outputs.
80 :
81 :
82 : */
83 : //+ENDPLUMEDOC
84 :
85 : class ANN : public Function {
86 : private:
87 : int num_layers;
88 : vector<int> num_nodes;
89 : vector<string> activations; // activation functions
90 : vector<vector<double> > weights; // flattened weight arrays
91 : vector<vector<double> > biases;
92 : vector<vector<double> > output_of_each_layer;
93 : vector<vector<double> > input_of_each_layer;
94 : vector<double** > coeff; // weight matrix arrays, reshaped from "weights"
95 :
96 : public:
97 : static void registerKeywords( Keywords& keys );
98 : explicit ANN(const ActionOptions&);
99 : virtual void calculate();
100 : void calculate_output_of_each_layer(const vector<double>& input);
101 : void back_prop(vector<vector<double> >& derivatives_of_each_layer, int index_of_output_component);
102 : };
103 :
104 : PLUMED_REGISTER_ACTION(ANN,"ANN")
105 :
106 9 : void ANN::registerKeywords( Keywords& keys ) {
107 9 : Function::registerKeywords(keys);
108 9 : keys.use("PERIODIC");
109 9 : keys.add("compulsory", "NUM_LAYERS", "number of layers of the neural network");
110 9 : keys.add("compulsory", "NUM_NODES", "numbers of nodes in each layer of the neural network");
111 9 : keys.add("compulsory", "ACTIVATIONS", "activation functions for the neural network");
112 9 : keys.add("numbered", "WEIGHTS", "flattened weight arrays connecting adjacent layers, "
113 : "WEIGHTS0 represents flattened weight array connecting layer 0 and layer 1, "
114 : "WEIGHTS1 represents flattened weight array connecting layer 1 and layer 2, ...");
115 9 : keys.add("numbered", "BIASES", "bias array for each layer of the neural network, "
116 : "BIASES0 represents bias array for layer 1, BIASES1 represents bias array for layer 2, ...");
117 : // since v2.2 plumed requires all components be registered
118 18 : keys.addOutputComponent("node", "default", "scalar", "components of ANN outputs");
119 9 : }
120 :
121 5 : ANN::ANN(const ActionOptions&ao):
122 : Action(ao),
123 5 : Function(ao) {
124 5 : parse("NUM_LAYERS", num_layers);
125 5 : num_nodes = vector<int>(num_layers);
126 5 : activations = vector<string>(num_layers - 1);
127 5 : output_of_each_layer = vector<vector<double> >(num_layers);
128 5 : input_of_each_layer = vector<vector<double> >(num_layers);
129 5 : coeff = vector<double** >(num_layers - 1);
130 5 : parseVector("NUM_NODES", num_nodes);
131 5 : parseVector("ACTIVATIONS", activations);
132 5 : log.printf("\nactivations = ");
133 15 : for (const auto & ss: activations) {
134 10 : log.printf("%s, ", ss.c_str());
135 : }
136 5 : log.printf("\nnum_nodes = ");
137 20 : for (auto ss: num_nodes) {
138 15 : log.printf("%d, ", ss);
139 : }
140 5 : vector<double> temp_single_coeff, temp_single_bias;
141 10 : for (int ii = 0; ; ii ++) {
142 : // parse coeff
143 30 : if( !parseNumberedVector("WEIGHTS", ii, temp_single_coeff) ) {
144 5 : temp_single_coeff=weights[ii-1];
145 : break;
146 : }
147 10 : weights.push_back(temp_single_coeff);
148 10 : log.printf("size of temp_single_coeff = %lu\n", temp_single_coeff.size());
149 10 : log.printf("size of weights = %lu\n", weights.size());
150 : // parse bias
151 20 : if( !parseNumberedVector("BIASES", ii, temp_single_bias) ) {
152 0 : temp_single_bias=biases[ii-1];
153 : }
154 10 : biases.push_back(temp_single_bias);
155 10 : log.printf("size of temp_single_bias = %lu\n", temp_single_bias.size());
156 10 : log.printf("size of biases = %lu\n", biases.size());
157 : }
158 :
159 5 : if(static_cast<int>(getNumberOfArguments()) != num_nodes[0]) {
160 0 : error("Number of arguments is wrong");
161 : }
162 :
163 5 : auto temp_coeff = weights;
164 15 : for (int ii = 0; ii < num_layers - 1; ii ++) {
165 : int num_of_rows, num_of_cols; // num of rows/cols for the coeff matrix of this connection
166 10 : num_of_rows = num_nodes[ii + 1];
167 10 : num_of_cols = num_nodes[ii];
168 : assert (num_of_rows * num_of_cols == static_cast<int>(temp_coeff[ii].size())); // check whether the size matches
169 : // create a 2d array to hold coefficients
170 10 : coeff[ii] = new double*[num_of_rows];
171 32 : for (int kk = 0; kk < num_of_rows; kk ++) {
172 22 : coeff[ii][kk] = new double[num_of_cols];
173 : }
174 61 : for (int jj = 0; jj < static_cast<int>(temp_coeff[ii].size()); jj ++) {
175 51 : coeff[ii][jj / num_of_cols][jj % num_of_cols] = temp_coeff[ii][jj];
176 : }
177 : }
178 : // check coeff
179 15 : for (int ii = 0; ii < num_layers - 1; ii ++) {
180 10 : log.printf("coeff %d = \n", ii);
181 32 : for (int jj = 0; jj < num_nodes[ii + 1]; jj ++) {
182 73 : for (int kk = 0; kk < num_nodes[ii]; kk ++) {
183 51 : log.printf("%f ", coeff[ii][jj][kk]);
184 : }
185 22 : log.printf("\n");
186 : }
187 : }
188 : // check bias
189 15 : for (int ii = 0; ii < num_layers - 1; ii ++) {
190 10 : log.printf("bias %d = \n", ii);
191 32 : for (int jj = 0; jj < num_nodes[ii + 1]; jj ++) {
192 22 : log.printf("%f ", biases[ii][jj]);
193 : }
194 10 : log.printf("\n");
195 : }
196 5 : log.printf("initialization ended\n");
197 : // create components
198 12 : for (int ii = 0; ii < num_nodes[num_layers - 1]; ii ++) {
199 14 : string name_of_this_component = "node-" + to_string(ii);
200 7 : addComponentWithDerivatives(name_of_this_component);
201 7 : componentIsNotPeriodic(name_of_this_component);
202 : }
203 5 : checkRead();
204 5 : }
205 :
206 1638 : void ANN::calculate_output_of_each_layer(const vector<double>& input) {
207 : // first layer
208 1638 : output_of_each_layer[0] = input;
209 : // following layers
210 4914 : for(int ii = 1; ii < static_cast<int>(num_nodes.size()); ii ++) {
211 3276 : output_of_each_layer[ii].resize(num_nodes[ii]);
212 3276 : input_of_each_layer[ii].resize(num_nodes[ii]);
213 : // first calculate input
214 10374 : for (int jj = 0; jj < num_nodes[ii]; jj ++) {
215 7098 : input_of_each_layer[ii][jj] = biases[ii - 1][jj]; // add bias term
216 23478 : for (int kk = 0; kk < num_nodes[ii - 1]; kk ++) {
217 16380 : input_of_each_layer[ii][jj] += coeff[ii - 1][jj][kk] * output_of_each_layer[ii - 1][kk];
218 : }
219 : }
220 : // then get output
221 3276 : if (activations[ii - 1] == string("Linear")) {
222 1092 : for(int jj = 0; jj < num_nodes[ii]; jj ++) {
223 546 : output_of_each_layer[ii][jj] = input_of_each_layer[ii][jj];
224 : }
225 2730 : } else if (activations[ii - 1] == string("Tanh")) {
226 7644 : for(int jj = 0; jj < num_nodes[ii]; jj ++) {
227 5460 : output_of_each_layer[ii][jj] = tanh(input_of_each_layer[ii][jj]);
228 : }
229 546 : } else if (activations[ii - 1] == string("Circular")) {
230 : assert (num_nodes[ii] % 2 == 0);
231 1092 : for(int jj = 0; jj < num_nodes[ii] / 2; jj ++) {
232 546 : double radius = sqrt(input_of_each_layer[ii][2 * jj] * input_of_each_layer[ii][2 * jj]
233 546 : +input_of_each_layer[ii][2 * jj + 1] * input_of_each_layer[ii][2 * jj + 1]);
234 546 : output_of_each_layer[ii][2 * jj] = input_of_each_layer[ii][2 * jj] / radius;
235 546 : output_of_each_layer[ii][2 * jj + 1] = input_of_each_layer[ii][2 * jj + 1] / radius;
236 :
237 : }
238 : } else {
239 : printf("layer type not found!\n\n");
240 0 : return;
241 : }
242 : }
243 : #ifdef DEBUG_2
244 : // print out the result for debugging
245 : printf("output_of_each_layer = \n");
246 : for (int ii = 0; ii < num_layers; ii ++) {
247 : printf("layer[%d]: ", ii);
248 : if (ii != 0) {
249 : cout << activations[ii - 1] << "\t";
250 : } else {
251 : cout << "input \t" ;
252 : }
253 : for (int jj = 0; jj < num_nodes[ii]; jj ++) {
254 : printf("%lf\t", output_of_each_layer[ii][jj]);
255 : }
256 : printf("\n");
257 : }
258 : printf("\n");
259 : #endif
260 : return;
261 : }
262 :
263 2184 : void ANN::back_prop(vector<vector<double> >& derivatives_of_each_layer, int index_of_output_component) {
264 2184 : derivatives_of_each_layer = output_of_each_layer; // the data structure and size should be the same, so I simply deep copy it
265 : // first calculate derivatives for bottleneck layer
266 5460 : for (int ii = 0; ii < num_nodes[num_nodes.size() - 1]; ii ++ ) {
267 3276 : if (ii == index_of_output_component) {
268 2184 : derivatives_of_each_layer[num_nodes.size() - 1][ii] = 1;
269 : } else {
270 1092 : derivatives_of_each_layer[num_nodes.size() - 1][ii] = 0;
271 : }
272 : }
273 : // the use back propagation to calculate derivatives for previous layers
274 6552 : for (int jj = num_nodes.size() - 2; jj >= 0; jj --) {
275 4368 : if (activations[jj] == string("Circular")) {
276 1092 : vector<double> temp_derivative_of_input_for_this_layer;
277 1092 : temp_derivative_of_input_for_this_layer.resize(num_nodes[jj + 1]);
278 : #ifdef DEBUG
279 : assert (num_nodes[jj + 1] % 2 == 0);
280 : #endif
281 : // first calculate the derivative of input from derivative of output of this circular layer
282 2184 : for(int ii = 0; ii < num_nodes[jj + 1] / 2; ii ++) {
283 1092 : double x_p = input_of_each_layer[jj + 1][2 * ii];
284 1092 : double x_q = input_of_each_layer[jj + 1][2 * ii + 1];
285 1092 : double radius = sqrt(x_p * x_p + x_q * x_q);
286 1092 : temp_derivative_of_input_for_this_layer[2 * ii] = x_q / (radius * radius * radius)
287 1092 : * (x_q * derivatives_of_each_layer[jj + 1][2 * ii]
288 1092 : - x_p * derivatives_of_each_layer[jj + 1][2 * ii + 1]);
289 1092 : temp_derivative_of_input_for_this_layer[2 * ii + 1] = x_p / (radius * radius * radius)
290 1092 : * (x_p * derivatives_of_each_layer[jj + 1][2 * ii + 1]
291 1092 : - x_q * derivatives_of_each_layer[jj + 1][2 * ii]);
292 : }
293 : #ifdef DEBUG
294 : for (int mm = 0; mm < num_nodes[jj + 1]; mm ++) {
295 : printf("temp_derivative_of_input_for_this_layer[%d] = %lf\n", mm, temp_derivative_of_input_for_this_layer[mm]);
296 : printf("derivatives_of_each_layer[%d + 1][%d] = %lf\n", jj, mm, derivatives_of_each_layer[jj + 1][mm]);
297 : }
298 : #endif
299 : // the calculate the derivative of output of layer jj, from derivative of input of layer (jj + 1)
300 4368 : for (int mm = 0; mm < num_nodes[jj]; mm ++) {
301 3276 : derivatives_of_each_layer[jj][mm] = 0;
302 9828 : for (int kk = 0; kk < num_nodes[jj + 1]; kk ++) {
303 6552 : derivatives_of_each_layer[jj][mm] += coeff[jj][kk][mm] \
304 6552 : * temp_derivative_of_input_for_this_layer[kk];
305 : #ifdef DEBUG
306 : printf("derivatives_of_each_layer[%d][%d] = %lf\n", jj, mm, derivatives_of_each_layer[jj][mm]);
307 : printf("coeff[%d][%d][%d] = %lf\n", jj, kk, mm, coeff[jj][kk][mm]);
308 : #endif
309 : }
310 : }
311 : // TODO: should be fine, pass all tests, although there seems to be some problems here previously
312 1092 : } else {
313 10920 : for (int mm = 0; mm < num_nodes[jj]; mm ++) {
314 7644 : derivatives_of_each_layer[jj][mm] = 0;
315 24024 : for (int kk = 0; kk < num_nodes[jj + 1]; kk ++) {
316 16380 : if (activations[jj] == string("Tanh")) {
317 14742 : derivatives_of_each_layer[jj][mm] += derivatives_of_each_layer[jj + 1][kk] \
318 14742 : * coeff[jj][kk][mm] \
319 14742 : * (1 - output_of_each_layer[jj + 1][kk] * output_of_each_layer[jj + 1][kk]);
320 1638 : } else if (activations[jj] == string("Linear")) {
321 1638 : derivatives_of_each_layer[jj][mm] += derivatives_of_each_layer[jj + 1][kk] \
322 1638 : * coeff[jj][kk][mm] \
323 1638 : * 1;
324 : } else {
325 : printf("layer type not found!\n\n");
326 0 : return;
327 : }
328 : }
329 : }
330 : }
331 : }
332 : #ifdef DEBUG
333 : // print out the result for debugging
334 : printf("derivatives_of_each_layer = \n");
335 : for (int ii = 0; ii < num_layers; ii ++) {
336 : printf("layer[%d]: ", ii);
337 : for (int jj = 0; jj < num_nodes[ii]; jj ++) {
338 : printf("%lf\t", derivatives_of_each_layer[ii][jj]);
339 : }
340 : printf("\n");
341 : }
342 : printf("\n");
343 : #endif
344 : return;
345 : }
346 :
347 1638 : void ANN::calculate() {
348 :
349 1638 : vector<double> input_layer_data(num_nodes[0]);
350 4914 : for (int ii = 0; ii < num_nodes[0]; ii ++) {
351 3276 : input_layer_data[ii] = getArgument(ii);
352 : }
353 :
354 1638 : calculate_output_of_each_layer(input_layer_data);
355 1638 : vector<vector<double> > derivatives_of_each_layer;
356 :
357 3822 : for (int ii = 0; ii < num_nodes[num_layers - 1]; ii ++) {
358 2184 : back_prop(derivatives_of_each_layer, ii);
359 2184 : string name_of_this_component = "node-" + to_string(ii);
360 2184 : Value* value_new=getPntrToComponent(name_of_this_component);
361 2184 : value_new -> set(output_of_each_layer[num_layers - 1][ii]);
362 6552 : for (int jj = 0; jj < num_nodes[0]; jj ++) {
363 4368 : value_new -> setDerivative(jj, derivatives_of_each_layer[0][jj]); // TODO: setDerivative or addDerivative?
364 : }
365 : #ifdef DEBUG_3
366 : printf("derivatives = ");
367 : for (int jj = 0; jj < num_nodes[0]; jj ++) {
368 : printf("%f ", value_new -> getDerivative(jj));
369 : }
370 : printf("\n");
371 : #endif
372 : }
373 :
374 1638 : }
375 :
376 : }
377 : }
378 : }
|
