Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/mlpbase.cs @ 3580

Last change on this file since 3580 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 126.6 KB
Line 
1/*************************************************************************
2Copyright (c) 2007-2008, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class mlpbase
26    {
27        public struct multilayerperceptron
28        {
29            public int[] structinfo;
30            public double[] weights;
31            public double[] columnmeans;
32            public double[] columnsigmas;
33            public double[] neurons;
34            public double[] dfdnet;
35            public double[] derror;
36            public double[] x;
37            public double[] y;
38            public double[,] chunks;
39            public double[] nwbuf;
40        };
41
42
43
44
45        public const int mlpvnum = 7;
46        public const int nfieldwidth = 4;
47        public const int chunksize = 32;
48
49
50        /*************************************************************************
51        Creates  neural  network  with  NIn  inputs,  NOut outputs, without hidden
52        layers, with linear output layer. Network weights are  filled  with  small
53        random values.
54
55          -- ALGLIB --
56             Copyright 04.11.2007 by Bochkanov Sergey
57        *************************************************************************/
58        public static void mlpcreate0(int nin,
59            int nout,
60            ref multilayerperceptron network)
61        {
62            int[] lsizes = new int[0];
63            int[] ltypes = new int[0];
64            int[] lconnfirst = new int[0];
65            int[] lconnlast = new int[0];
66            int layerscount = 0;
67            int lastproc = 0;
68
69            layerscount = 1+2;
70           
71            //
72            // Allocate arrays
73            //
74            lsizes = new int[layerscount-1+1];
75            ltypes = new int[layerscount-1+1];
76            lconnfirst = new int[layerscount-1+1];
77            lconnlast = new int[layerscount-1+1];
78           
79            //
80            // Layers
81            //
82            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
83            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
84           
85            //
86            // Create
87            //
88            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
89        }
90
91
92        /*************************************************************************
93        Same  as  MLPCreate0,  but  with  one  hidden  layer  (NHid  neurons) with
94        non-linear activation function. Output layer is linear.
95
96          -- ALGLIB --
97             Copyright 04.11.2007 by Bochkanov Sergey
98        *************************************************************************/
99        public static void mlpcreate1(int nin,
100            int nhid,
101            int nout,
102            ref multilayerperceptron network)
103        {
104            int[] lsizes = new int[0];
105            int[] ltypes = new int[0];
106            int[] lconnfirst = new int[0];
107            int[] lconnlast = new int[0];
108            int layerscount = 0;
109            int lastproc = 0;
110
111            layerscount = 1+3+2;
112           
113            //
114            // Allocate arrays
115            //
116            lsizes = new int[layerscount-1+1];
117            ltypes = new int[layerscount-1+1];
118            lconnfirst = new int[layerscount-1+1];
119            lconnlast = new int[layerscount-1+1];
120           
121            //
122            // Layers
123            //
124            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
125            addbiasedsummatorlayer(nhid, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
126            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
127            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
128           
129            //
130            // Create
131            //
132            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
133        }
134
135
136        /*************************************************************************
137        Same as MLPCreate0, but with two hidden layers (NHid1 and  NHid2  neurons)
138        with non-linear activation function. Output layer is linear.
139         $ALL
140
141          -- ALGLIB --
142             Copyright 04.11.2007 by Bochkanov Sergey
143        *************************************************************************/
144        public static void mlpcreate2(int nin,
145            int nhid1,
146            int nhid2,
147            int nout,
148            ref multilayerperceptron network)
149        {
150            int[] lsizes = new int[0];
151            int[] ltypes = new int[0];
152            int[] lconnfirst = new int[0];
153            int[] lconnlast = new int[0];
154            int layerscount = 0;
155            int lastproc = 0;
156
157            layerscount = 1+3+3+2;
158           
159            //
160            // Allocate arrays
161            //
162            lsizes = new int[layerscount-1+1];
163            ltypes = new int[layerscount-1+1];
164            lconnfirst = new int[layerscount-1+1];
165            lconnlast = new int[layerscount-1+1];
166           
167            //
168            // Layers
169            //
170            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
171            addbiasedsummatorlayer(nhid1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
172            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
173            addbiasedsummatorlayer(nhid2, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
174            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
175            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
176           
177            //
178            // Create
179            //
180            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
181        }
182
183
184        /*************************************************************************
185        Creates  neural  network  with  NIn  inputs,  NOut outputs, without hidden
186        layers with non-linear output layer. Network weights are filled with small
187        random values.
188
189        Activation function of the output layer takes values:
190
191            (B, +INF), if D>=0
192
193        or
194
195            (-INF, B), if D<0.
196
197
198          -- ALGLIB --
199             Copyright 30.03.2008 by Bochkanov Sergey
200        *************************************************************************/
201        public static void mlpcreateb0(int nin,
202            int nout,
203            double b,
204            double d,
205            ref multilayerperceptron network)
206        {
207            int[] lsizes = new int[0];
208            int[] ltypes = new int[0];
209            int[] lconnfirst = new int[0];
210            int[] lconnlast = new int[0];
211            int layerscount = 0;
212            int lastproc = 0;
213            int i = 0;
214
215            layerscount = 1+3;
216            if( (double)(d)>=(double)(0) )
217            {
218                d = 1;
219            }
220            else
221            {
222                d = -1;
223            }
224           
225            //
226            // Allocate arrays
227            //
228            lsizes = new int[layerscount-1+1];
229            ltypes = new int[layerscount-1+1];
230            lconnfirst = new int[layerscount-1+1];
231            lconnlast = new int[layerscount-1+1];
232           
233            //
234            // Layers
235            //
236            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
237            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
238            addactivationlayer(3, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
239           
240            //
241            // Create
242            //
243            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
244           
245            //
246            // Turn on ouputs shift/scaling.
247            //
248            for(i=nin; i<=nin+nout-1; i++)
249            {
250                network.columnmeans[i] = b;
251                network.columnsigmas[i] = d;
252            }
253        }
254
255
256        /*************************************************************************
257        Same as MLPCreateB0 but with non-linear hidden layer.
258
259          -- ALGLIB --
260             Copyright 30.03.2008 by Bochkanov Sergey
261        *************************************************************************/
262        public static void mlpcreateb1(int nin,
263            int nhid,
264            int nout,
265            double b,
266            double d,
267            ref multilayerperceptron network)
268        {
269            int[] lsizes = new int[0];
270            int[] ltypes = new int[0];
271            int[] lconnfirst = new int[0];
272            int[] lconnlast = new int[0];
273            int layerscount = 0;
274            int lastproc = 0;
275            int i = 0;
276
277            layerscount = 1+3+3;
278            if( (double)(d)>=(double)(0) )
279            {
280                d = 1;
281            }
282            else
283            {
284                d = -1;
285            }
286           
287            //
288            // Allocate arrays
289            //
290            lsizes = new int[layerscount-1+1];
291            ltypes = new int[layerscount-1+1];
292            lconnfirst = new int[layerscount-1+1];
293            lconnlast = new int[layerscount-1+1];
294           
295            //
296            // Layers
297            //
298            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
299            addbiasedsummatorlayer(nhid, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
300            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
301            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
302            addactivationlayer(3, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
303           
304            //
305            // Create
306            //
307            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
308           
309            //
310            // Turn on ouputs shift/scaling.
311            //
312            for(i=nin; i<=nin+nout-1; i++)
313            {
314                network.columnmeans[i] = b;
315                network.columnsigmas[i] = d;
316            }
317        }
318
319
320        /*************************************************************************
321        Same as MLPCreateB0 but with two non-linear hidden layers.
322
323          -- ALGLIB --
324             Copyright 30.03.2008 by Bochkanov Sergey
325        *************************************************************************/
326        public static void mlpcreateb2(int nin,
327            int nhid1,
328            int nhid2,
329            int nout,
330            double b,
331            double d,
332            ref multilayerperceptron network)
333        {
334            int[] lsizes = new int[0];
335            int[] ltypes = new int[0];
336            int[] lconnfirst = new int[0];
337            int[] lconnlast = new int[0];
338            int layerscount = 0;
339            int lastproc = 0;
340            int i = 0;
341
342            layerscount = 1+3+3+3;
343            if( (double)(d)>=(double)(0) )
344            {
345                d = 1;
346            }
347            else
348            {
349                d = -1;
350            }
351           
352            //
353            // Allocate arrays
354            //
355            lsizes = new int[layerscount-1+1];
356            ltypes = new int[layerscount-1+1];
357            lconnfirst = new int[layerscount-1+1];
358            lconnlast = new int[layerscount-1+1];
359           
360            //
361            // Layers
362            //
363            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
364            addbiasedsummatorlayer(nhid1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
365            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
366            addbiasedsummatorlayer(nhid2, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
367            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
368            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
369            addactivationlayer(3, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
370           
371            //
372            // Create
373            //
374            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
375           
376            //
377            // Turn on ouputs shift/scaling.
378            //
379            for(i=nin; i<=nin+nout-1; i++)
380            {
381                network.columnmeans[i] = b;
382                network.columnsigmas[i] = d;
383            }
384        }
385
386
387        /*************************************************************************
388        Creates  neural  network  with  NIn  inputs,  NOut outputs, without hidden
389        layers with non-linear output layer. Network weights are filled with small
390        random values. Activation function of the output layer takes values [A,B].
391
392          -- ALGLIB --
393             Copyright 30.03.2008 by Bochkanov Sergey
394        *************************************************************************/
395        public static void mlpcreater0(int nin,
396            int nout,
397            double a,
398            double b,
399            ref multilayerperceptron network)
400        {
401            int[] lsizes = new int[0];
402            int[] ltypes = new int[0];
403            int[] lconnfirst = new int[0];
404            int[] lconnlast = new int[0];
405            int layerscount = 0;
406            int lastproc = 0;
407            int i = 0;
408
409            layerscount = 1+3;
410           
411            //
412            // Allocate arrays
413            //
414            lsizes = new int[layerscount-1+1];
415            ltypes = new int[layerscount-1+1];
416            lconnfirst = new int[layerscount-1+1];
417            lconnlast = new int[layerscount-1+1];
418           
419            //
420            // Layers
421            //
422            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
423            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
424            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
425           
426            //
427            // Create
428            //
429            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
430           
431            //
432            // Turn on outputs shift/scaling.
433            //
434            for(i=nin; i<=nin+nout-1; i++)
435            {
436                network.columnmeans[i] = 0.5*(a+b);
437                network.columnsigmas[i] = 0.5*(a-b);
438            }
439        }
440
441
442        /*************************************************************************
443        Same as MLPCreateR0, but with non-linear hidden layer.
444
445          -- ALGLIB --
446             Copyright 30.03.2008 by Bochkanov Sergey
447        *************************************************************************/
448        public static void mlpcreater1(int nin,
449            int nhid,
450            int nout,
451            double a,
452            double b,
453            ref multilayerperceptron network)
454        {
455            int[] lsizes = new int[0];
456            int[] ltypes = new int[0];
457            int[] lconnfirst = new int[0];
458            int[] lconnlast = new int[0];
459            int layerscount = 0;
460            int lastproc = 0;
461            int i = 0;
462
463            layerscount = 1+3+3;
464           
465            //
466            // Allocate arrays
467            //
468            lsizes = new int[layerscount-1+1];
469            ltypes = new int[layerscount-1+1];
470            lconnfirst = new int[layerscount-1+1];
471            lconnlast = new int[layerscount-1+1];
472           
473            //
474            // Layers
475            //
476            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
477            addbiasedsummatorlayer(nhid, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
478            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
479            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
480            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
481           
482            //
483            // Create
484            //
485            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
486           
487            //
488            // Turn on outputs shift/scaling.
489            //
490            for(i=nin; i<=nin+nout-1; i++)
491            {
492                network.columnmeans[i] = 0.5*(a+b);
493                network.columnsigmas[i] = 0.5*(a-b);
494            }
495        }
496
497
498        /*************************************************************************
499        Same as MLPCreateR0, but with two non-linear hidden layers.
500
501          -- ALGLIB --
502             Copyright 30.03.2008 by Bochkanov Sergey
503        *************************************************************************/
504        public static void mlpcreater2(int nin,
505            int nhid1,
506            int nhid2,
507            int nout,
508            double a,
509            double b,
510            ref multilayerperceptron network)
511        {
512            int[] lsizes = new int[0];
513            int[] ltypes = new int[0];
514            int[] lconnfirst = new int[0];
515            int[] lconnlast = new int[0];
516            int layerscount = 0;
517            int lastproc = 0;
518            int i = 0;
519
520            layerscount = 1+3+3+3;
521           
522            //
523            // Allocate arrays
524            //
525            lsizes = new int[layerscount-1+1];
526            ltypes = new int[layerscount-1+1];
527            lconnfirst = new int[layerscount-1+1];
528            lconnlast = new int[layerscount-1+1];
529           
530            //
531            // Layers
532            //
533            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
534            addbiasedsummatorlayer(nhid1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
535            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
536            addbiasedsummatorlayer(nhid2, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
537            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
538            addbiasedsummatorlayer(nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
539            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
540           
541            //
542            // Create
543            //
544            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, false, ref network);
545           
546            //
547            // Turn on outputs shift/scaling.
548            //
549            for(i=nin; i<=nin+nout-1; i++)
550            {
551                network.columnmeans[i] = 0.5*(a+b);
552                network.columnsigmas[i] = 0.5*(a-b);
553            }
554        }
555
556
557        /*************************************************************************
558        Creates classifier network with NIn  inputs  and  NOut  possible  classes.
559        Network contains no hidden layers and linear output  layer  with  SOFTMAX-
560        normalization  (so  outputs  sums  up  to  1.0  and  converge to posterior
561        probabilities).
562
563          -- ALGLIB --
564             Copyright 04.11.2007 by Bochkanov Sergey
565        *************************************************************************/
566        public static void mlpcreatec0(int nin,
567            int nout,
568            ref multilayerperceptron network)
569        {
570            int[] lsizes = new int[0];
571            int[] ltypes = new int[0];
572            int[] lconnfirst = new int[0];
573            int[] lconnlast = new int[0];
574            int layerscount = 0;
575            int lastproc = 0;
576
577            System.Diagnostics.Debug.Assert(nout>=2, "MLPCreateC0: NOut<2!");
578            layerscount = 1+2+1;
579           
580            //
581            // Allocate arrays
582            //
583            lsizes = new int[layerscount-1+1];
584            ltypes = new int[layerscount-1+1];
585            lconnfirst = new int[layerscount-1+1];
586            lconnlast = new int[layerscount-1+1];
587           
588            //
589            // Layers
590            //
591            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
592            addbiasedsummatorlayer(nout-1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
593            addzerolayer(ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
594           
595            //
596            // Create
597            //
598            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, true, ref network);
599        }
600
601
602        /*************************************************************************
603        Same as MLPCreateC0, but with one non-linear hidden layer.
604
605          -- ALGLIB --
606             Copyright 04.11.2007 by Bochkanov Sergey
607        *************************************************************************/
608        public static void mlpcreatec1(int nin,
609            int nhid,
610            int nout,
611            ref multilayerperceptron network)
612        {
613            int[] lsizes = new int[0];
614            int[] ltypes = new int[0];
615            int[] lconnfirst = new int[0];
616            int[] lconnlast = new int[0];
617            int layerscount = 0;
618            int lastproc = 0;
619
620            System.Diagnostics.Debug.Assert(nout>=2, "MLPCreateC1: NOut<2!");
621            layerscount = 1+3+2+1;
622           
623            //
624            // Allocate arrays
625            //
626            lsizes = new int[layerscount-1+1];
627            ltypes = new int[layerscount-1+1];
628            lconnfirst = new int[layerscount-1+1];
629            lconnlast = new int[layerscount-1+1];
630           
631            //
632            // Layers
633            //
634            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
635            addbiasedsummatorlayer(nhid, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
636            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
637            addbiasedsummatorlayer(nout-1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
638            addzerolayer(ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
639           
640            //
641            // Create
642            //
643            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, true, ref network);
644        }
645
646
647        /*************************************************************************
648        Same as MLPCreateC0, but with two non-linear hidden layers.
649
650          -- ALGLIB --
651             Copyright 04.11.2007 by Bochkanov Sergey
652        *************************************************************************/
653        public static void mlpcreatec2(int nin,
654            int nhid1,
655            int nhid2,
656            int nout,
657            ref multilayerperceptron network)
658        {
659            int[] lsizes = new int[0];
660            int[] ltypes = new int[0];
661            int[] lconnfirst = new int[0];
662            int[] lconnlast = new int[0];
663            int layerscount = 0;
664            int lastproc = 0;
665
666            System.Diagnostics.Debug.Assert(nout>=2, "MLPCreateC2: NOut<2!");
667            layerscount = 1+3+3+2+1;
668           
669            //
670            // Allocate arrays
671            //
672            lsizes = new int[layerscount-1+1];
673            ltypes = new int[layerscount-1+1];
674            lconnfirst = new int[layerscount-1+1];
675            lconnlast = new int[layerscount-1+1];
676           
677            //
678            // Layers
679            //
680            addinputlayer(nin, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
681            addbiasedsummatorlayer(nhid1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
682            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
683            addbiasedsummatorlayer(nhid2, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
684            addactivationlayer(1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
685            addbiasedsummatorlayer(nout-1, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
686            addzerolayer(ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, ref lastproc);
687           
688            //
689            // Create
690            //
691            mlpcreate(nin, nout, ref lsizes, ref ltypes, ref lconnfirst, ref lconnlast, layerscount, true, ref network);
692        }
693
694
695        /*************************************************************************
696        Copying of neural network
697
698        INPUT PARAMETERS:
699            Network1 -   original
700
701        OUTPUT PARAMETERS:
702            Network2 -   copy
703
704          -- ALGLIB --
705             Copyright 04.11.2007 by Bochkanov Sergey
706        *************************************************************************/
707        public static void mlpcopy(ref multilayerperceptron network1,
708            ref multilayerperceptron network2)
709        {
710            int i = 0;
711            int ssize = 0;
712            int ntotal = 0;
713            int nin = 0;
714            int nout = 0;
715            int wcount = 0;
716            int i_ = 0;
717
718           
719            //
720            // Unload info
721            //
722            ssize = network1.structinfo[0];
723            nin = network1.structinfo[1];
724            nout = network1.structinfo[2];
725            ntotal = network1.structinfo[3];
726            wcount = network1.structinfo[4];
727           
728            //
729            // Allocate space
730            //
731            network2.structinfo = new int[ssize-1+1];
732            network2.weights = new double[wcount-1+1];
733            if( mlpissoftmax(ref network1) )
734            {
735                network2.columnmeans = new double[nin-1+1];
736                network2.columnsigmas = new double[nin-1+1];
737            }
738            else
739            {
740                network2.columnmeans = new double[nin+nout-1+1];
741                network2.columnsigmas = new double[nin+nout-1+1];
742            }
743            network2.neurons = new double[ntotal-1+1];
744            network2.chunks = new double[3*ntotal+1, chunksize-1+1];
745            network2.nwbuf = new double[Math.Max(wcount, 2*nout)-1+1];
746            network2.dfdnet = new double[ntotal-1+1];
747            network2.x = new double[nin-1+1];
748            network2.y = new double[nout-1+1];
749            network2.derror = new double[ntotal-1+1];
750           
751            //
752            // Copy
753            //
754            for(i=0; i<=ssize-1; i++)
755            {
756                network2.structinfo[i] = network1.structinfo[i];
757            }
758            for(i_=0; i_<=wcount-1;i_++)
759            {
760                network2.weights[i_] = network1.weights[i_];
761            }
762            if( mlpissoftmax(ref network1) )
763            {
764                for(i_=0; i_<=nin-1;i_++)
765                {
766                    network2.columnmeans[i_] = network1.columnmeans[i_];
767                }
768                for(i_=0; i_<=nin-1;i_++)
769                {
770                    network2.columnsigmas[i_] = network1.columnsigmas[i_];
771                }
772            }
773            else
774            {
775                for(i_=0; i_<=nin+nout-1;i_++)
776                {
777                    network2.columnmeans[i_] = network1.columnmeans[i_];
778                }
779                for(i_=0; i_<=nin+nout-1;i_++)
780                {
781                    network2.columnsigmas[i_] = network1.columnsigmas[i_];
782                }
783            }
784            for(i_=0; i_<=ntotal-1;i_++)
785            {
786                network2.neurons[i_] = network1.neurons[i_];
787            }
788            for(i_=0; i_<=ntotal-1;i_++)
789            {
790                network2.dfdnet[i_] = network1.dfdnet[i_];
791            }
792            for(i_=0; i_<=nin-1;i_++)
793            {
794                network2.x[i_] = network1.x[i_];
795            }
796            for(i_=0; i_<=nout-1;i_++)
797            {
798                network2.y[i_] = network1.y[i_];
799            }
800            for(i_=0; i_<=ntotal-1;i_++)
801            {
802                network2.derror[i_] = network1.derror[i_];
803            }
804        }
805
806
807        /*************************************************************************
808        Serialization of MultiLayerPerceptron strucure
809
810        INPUT PARAMETERS:
811            Network -   original
812
813        OUTPUT PARAMETERS:
814            RA      -   array of real numbers which stores network,
815                        array[0..RLen-1]
816            RLen    -   RA lenght
817
818          -- ALGLIB --
819             Copyright 29.03.2008 by Bochkanov Sergey
820        *************************************************************************/
821        public static void mlpserialize(ref multilayerperceptron network,
822            ref double[] ra,
823            ref int rlen)
824        {
825            int i = 0;
826            int ssize = 0;
827            int ntotal = 0;
828            int nin = 0;
829            int nout = 0;
830            int wcount = 0;
831            int sigmalen = 0;
832            int offs = 0;
833            int i_ = 0;
834            int i1_ = 0;
835
836           
837            //
838            // Unload info
839            //
840            ssize = network.structinfo[0];
841            nin = network.structinfo[1];
842            nout = network.structinfo[2];
843            ntotal = network.structinfo[3];
844            wcount = network.structinfo[4];
845            if( mlpissoftmax(ref network) )
846            {
847                sigmalen = nin;
848            }
849            else
850            {
851                sigmalen = nin+nout;
852            }
853           
854            //
855            //  RA format:
856            //      LEN         DESRC.
857            //      1           RLen
858            //      1           version (MLPVNum)
859            //      1           StructInfo size
860            //      SSize       StructInfo
861            //      WCount      Weights
862            //      SigmaLen    ColumnMeans
863            //      SigmaLen    ColumnSigmas
864            //
865            rlen = 3+ssize+wcount+2*sigmalen;
866            ra = new double[rlen-1+1];
867            ra[0] = rlen;
868            ra[1] = mlpvnum;
869            ra[2] = ssize;
870            offs = 3;
871            for(i=0; i<=ssize-1; i++)
872            {
873                ra[offs+i] = network.structinfo[i];
874            }
875            offs = offs+ssize;
876            i1_ = (0) - (offs);
877            for(i_=offs; i_<=offs+wcount-1;i_++)
878            {
879                ra[i_] = network.weights[i_+i1_];
880            }
881            offs = offs+wcount;
882            i1_ = (0) - (offs);
883            for(i_=offs; i_<=offs+sigmalen-1;i_++)
884            {
885                ra[i_] = network.columnmeans[i_+i1_];
886            }
887            offs = offs+sigmalen;
888            i1_ = (0) - (offs);
889            for(i_=offs; i_<=offs+sigmalen-1;i_++)
890            {
891                ra[i_] = network.columnsigmas[i_+i1_];
892            }
893            offs = offs+sigmalen;
894        }
895
896
897        /*************************************************************************
898        Unserialization of MultiLayerPerceptron strucure
899
900        INPUT PARAMETERS:
901            RA      -   real array which stores network
902
903        OUTPUT PARAMETERS:
904            Network -   restored network
905
906          -- ALGLIB --
907             Copyright 29.03.2008 by Bochkanov Sergey
908        *************************************************************************/
909        public static void mlpunserialize(ref double[] ra,
910            ref multilayerperceptron network)
911        {
912            int i = 0;
913            int ssize = 0;
914            int ntotal = 0;
915            int nin = 0;
916            int nout = 0;
917            int wcount = 0;
918            int sigmalen = 0;
919            int offs = 0;
920            int i_ = 0;
921            int i1_ = 0;
922
923            System.Diagnostics.Debug.Assert((int)Math.Round(ra[1])==mlpvnum, "MLPUnserialize: incorrect array!");
924           
925            //
926            // Unload StructInfo from IA
927            //
928            offs = 3;
929            ssize = (int)Math.Round(ra[2]);
930            network.structinfo = new int[ssize-1+1];
931            for(i=0; i<=ssize-1; i++)
932            {
933                network.structinfo[i] = (int)Math.Round(ra[offs+i]);
934            }
935            offs = offs+ssize;
936           
937            //
938            // Unload info from StructInfo
939            //
940            ssize = network.structinfo[0];
941            nin = network.structinfo[1];
942            nout = network.structinfo[2];
943            ntotal = network.structinfo[3];
944            wcount = network.structinfo[4];
945            if( network.structinfo[6]==0 )
946            {
947                sigmalen = nin+nout;
948            }
949            else
950            {
951                sigmalen = nin;
952            }
953           
954            //
955            // Allocate space for other fields
956            //
957            network.weights = new double[wcount-1+1];
958            network.columnmeans = new double[sigmalen-1+1];
959            network.columnsigmas = new double[sigmalen-1+1];
960            network.neurons = new double[ntotal-1+1];
961            network.chunks = new double[3*ntotal+1, chunksize-1+1];
962            network.nwbuf = new double[Math.Max(wcount, 2*nout)-1+1];
963            network.dfdnet = new double[ntotal-1+1];
964            network.x = new double[nin-1+1];
965            network.y = new double[nout-1+1];
966            network.derror = new double[ntotal-1+1];
967           
968            //
969            // Copy parameters from RA
970            //
971            i1_ = (offs) - (0);
972            for(i_=0; i_<=wcount-1;i_++)
973            {
974                network.weights[i_] = ra[i_+i1_];
975            }
976            offs = offs+wcount;
977            i1_ = (offs) - (0);
978            for(i_=0; i_<=sigmalen-1;i_++)
979            {
980                network.columnmeans[i_] = ra[i_+i1_];
981            }
982            offs = offs+sigmalen;
983            i1_ = (offs) - (0);
984            for(i_=0; i_<=sigmalen-1;i_++)
985            {
986                network.columnsigmas[i_] = ra[i_+i1_];
987            }
988            offs = offs+sigmalen;
989        }
990
991
992        /*************************************************************************
993        Randomization of neural network weights
994
995          -- ALGLIB --
996             Copyright 06.11.2007 by Bochkanov Sergey
997        *************************************************************************/
998        public static void mlprandomize(ref multilayerperceptron network)
999        {
1000            int i = 0;
1001            int nin = 0;
1002            int nout = 0;
1003            int wcount = 0;
1004
1005            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1006            for(i=0; i<=wcount-1; i++)
1007            {
1008                network.weights[i] = AP.Math.RandomReal()-0.5;
1009            }
1010        }
1011
1012
1013        /*************************************************************************
1014        Randomization of neural network weights and standartisator
1015
1016          -- ALGLIB --
1017             Copyright 10.03.2008 by Bochkanov Sergey
1018        *************************************************************************/
1019        public static void mlprandomizefull(ref multilayerperceptron network)
1020        {
1021            int i = 0;
1022            int nin = 0;
1023            int nout = 0;
1024            int wcount = 0;
1025            int ntotal = 0;
1026            int istart = 0;
1027            int offs = 0;
1028            int ntype = 0;
1029
1030            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1031            ntotal = network.structinfo[3];
1032            istart = network.structinfo[5];
1033           
1034            //
1035            // Process network
1036            //
1037            for(i=0; i<=wcount-1; i++)
1038            {
1039                network.weights[i] = AP.Math.RandomReal()-0.5;
1040            }
1041            for(i=0; i<=nin-1; i++)
1042            {
1043                network.columnmeans[i] = 2*AP.Math.RandomReal()-1;
1044                network.columnsigmas[i] = 1.5*AP.Math.RandomReal()+0.5;
1045            }
1046            if( !mlpissoftmax(ref network) )
1047            {
1048                for(i=0; i<=nout-1; i++)
1049                {
1050                    offs = istart+(ntotal-nout+i)*nfieldwidth;
1051                    ntype = network.structinfo[offs+0];
1052                    if( ntype==0 )
1053                    {
1054                       
1055                        //
1056                        // Shifts are changed only for linear outputs neurons
1057                        //
1058                        network.columnmeans[nin+i] = 2*AP.Math.RandomReal()-1;
1059                    }
1060                    if( ntype==0 | ntype==3 )
1061                    {
1062                       
1063                        //
1064                        // Scales are changed only for linear or bounded outputs neurons.
1065                        // Note that scale randomization preserves sign.
1066                        //
1067                        network.columnsigmas[nin+i] = Math.Sign(network.columnsigmas[nin+i])*(1.5*AP.Math.RandomReal()+0.5);
1068                    }
1069                }
1070            }
1071        }
1072
1073
1074        /*************************************************************************
1075        Internal subroutine.
1076
1077          -- ALGLIB --
1078             Copyright 30.03.2008 by Bochkanov Sergey
1079        *************************************************************************/
1080        public static void mlpinitpreprocessor(ref multilayerperceptron network,
1081            ref double[,] xy,
1082            int ssize)
1083        {
1084            int i = 0;
1085            int j = 0;
1086            int jmax = 0;
1087            int nin = 0;
1088            int nout = 0;
1089            int wcount = 0;
1090            int ntotal = 0;
1091            int istart = 0;
1092            int offs = 0;
1093            int ntype = 0;
1094            double[] means = new double[0];
1095            double[] sigmas = new double[0];
1096            double s = 0;
1097
1098            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1099            ntotal = network.structinfo[3];
1100            istart = network.structinfo[5];
1101           
1102            //
1103            // Means/Sigmas
1104            //
1105            if( mlpissoftmax(ref network) )
1106            {
1107                jmax = nin-1;
1108            }
1109            else
1110            {
1111                jmax = nin+nout-1;
1112            }
1113            means = new double[jmax+1];
1114            sigmas = new double[jmax+1];
1115            for(j=0; j<=jmax; j++)
1116            {
1117                means[j] = 0;
1118                for(i=0; i<=ssize-1; i++)
1119                {
1120                    means[j] = means[j]+xy[i,j];
1121                }
1122                means[j] = means[j]/ssize;
1123                sigmas[j] = 0;
1124                for(i=0; i<=ssize-1; i++)
1125                {
1126                    sigmas[j] = sigmas[j]+AP.Math.Sqr(xy[i,j]-means[j]);
1127                }
1128                sigmas[j] = Math.Sqrt(sigmas[j]/ssize);
1129            }
1130           
1131            //
1132            // Inputs
1133            //
1134            for(i=0; i<=nin-1; i++)
1135            {
1136                network.columnmeans[i] = means[i];
1137                network.columnsigmas[i] = sigmas[i];
1138                if( (double)(network.columnsigmas[i])==(double)(0) )
1139                {
1140                    network.columnsigmas[i] = 1;
1141                }
1142            }
1143           
1144            //
1145            // Outputs
1146            //
1147            if( !mlpissoftmax(ref network) )
1148            {
1149                for(i=0; i<=nout-1; i++)
1150                {
1151                    offs = istart+(ntotal-nout+i)*nfieldwidth;
1152                    ntype = network.structinfo[offs+0];
1153                   
1154                    //
1155                    // Linear outputs
1156                    //
1157                    if( ntype==0 )
1158                    {
1159                        network.columnmeans[nin+i] = means[nin+i];
1160                        network.columnsigmas[nin+i] = sigmas[nin+i];
1161                        if( (double)(network.columnsigmas[nin+i])==(double)(0) )
1162                        {
1163                            network.columnsigmas[nin+i] = 1;
1164                        }
1165                    }
1166                   
1167                    //
1168                    // Bounded outputs (half-interval)
1169                    //
1170                    if( ntype==3 )
1171                    {
1172                        s = means[nin+i]-network.columnmeans[nin+i];
1173                        if( (double)(s)==(double)(0) )
1174                        {
1175                            s = Math.Sign(network.columnsigmas[nin+i]);
1176                        }
1177                        if( (double)(s)==(double)(0) )
1178                        {
1179                            s = 1.0;
1180                        }
1181                        network.columnsigmas[nin+i] = Math.Sign(network.columnsigmas[nin+i])*Math.Abs(s);
1182                        if( (double)(network.columnsigmas[nin+i])==(double)(0) )
1183                        {
1184                            network.columnsigmas[nin+i] = 1;
1185                        }
1186                    }
1187                }
1188            }
1189        }
1190
1191
1192        /*************************************************************************
1193        Returns information about initialized network: number of inputs, outputs,
1194        weights.
1195
1196          -- ALGLIB --
1197             Copyright 04.11.2007 by Bochkanov Sergey
1198        *************************************************************************/
1199        public static void mlpproperties(ref multilayerperceptron network,
1200            ref int nin,
1201            ref int nout,
1202            ref int wcount)
1203        {
1204            nin = network.structinfo[1];
1205            nout = network.structinfo[2];
1206            wcount = network.structinfo[4];
1207        }
1208
1209
1210        /*************************************************************************
1211        Tells whether network is SOFTMAX-normalized (i.e. classifier) or not.
1212
1213          -- ALGLIB --
1214             Copyright 04.11.2007 by Bochkanov Sergey
1215        *************************************************************************/
1216        public static bool mlpissoftmax(ref multilayerperceptron network)
1217        {
1218            bool result = new bool();
1219
1220            result = network.structinfo[6]==1;
1221            return result;
1222        }
1223
1224
1225        /*************************************************************************
1226        Procesing
1227
1228        INPUT PARAMETERS:
1229            Network -   neural network
1230            X       -   input vector,  array[0..NIn-1].
1231
1232        OUTPUT PARAMETERS:
1233            Y       -   result. Regression estimate when solving regression  task,
1234                        vector of posterior probabilities for classification task.
1235                        Subroutine does not allocate memory for this vector, it is
1236                        responsibility of a caller to allocate it. Array  must  be
1237                        at least [0..NOut-1].
1238
1239          -- ALGLIB --
1240             Copyright 04.11.2007 by Bochkanov Sergey
1241        *************************************************************************/
1242        public static void mlpprocess(ref multilayerperceptron network,
1243            ref double[] x,
1244            ref double[] y)
1245        {
1246            mlpinternalprocessvector(ref network.structinfo, ref network.weights, ref network.columnmeans, ref network.columnsigmas, ref network.neurons, ref network.dfdnet, ref x, ref y);
1247        }
1248
1249
1250        /*************************************************************************
1251        Error function for neural network, internal subroutine.
1252
1253          -- ALGLIB --
1254             Copyright 04.11.2007 by Bochkanov Sergey
1255        *************************************************************************/
1256        public static double mlperror(ref multilayerperceptron network,
1257            ref double[,] xy,
1258            int ssize)
1259        {
1260            double result = 0;
1261            int i = 0;
1262            int k = 0;
1263            int nin = 0;
1264            int nout = 0;
1265            int wcount = 0;
1266            double e = 0;
1267            int i_ = 0;
1268            int i1_ = 0;
1269
1270            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1271            result = 0;
1272            for(i=0; i<=ssize-1; i++)
1273            {
1274                for(i_=0; i_<=nin-1;i_++)
1275                {
1276                    network.x[i_] = xy[i,i_];
1277                }
1278                mlpprocess(ref network, ref network.x, ref network.y);
1279                if( mlpissoftmax(ref network) )
1280                {
1281                   
1282                    //
1283                    // class labels outputs
1284                    //
1285                    k = (int)Math.Round(xy[i,nin]);
1286                    if( k>=0 & k<nout )
1287                    {
1288                        network.y[k] = network.y[k]-1;
1289                    }
1290                }
1291                else
1292                {
1293                   
1294                    //
1295                    // real outputs
1296                    //
1297                    i1_ = (nin) - (0);
1298                    for(i_=0; i_<=nout-1;i_++)
1299                    {
1300                        network.y[i_] = network.y[i_] - xy[i,i_+i1_];
1301                    }
1302                }
1303                e = 0.0;
1304                for(i_=0; i_<=nout-1;i_++)
1305                {
1306                    e += network.y[i_]*network.y[i_];
1307                }
1308                result = result+e/2;
1309            }
1310            return result;
1311        }
1312
1313
1314        /*************************************************************************
1315        Natural error function for neural network, internal subroutine.
1316
1317          -- ALGLIB --
1318             Copyright 04.11.2007 by Bochkanov Sergey
1319        *************************************************************************/
1320        public static double mlperrorn(ref multilayerperceptron network,
1321            ref double[,] xy,
1322            int ssize)
1323        {
1324            double result = 0;
1325            int i = 0;
1326            int j = 0;
1327            int k = 0;
1328            int nin = 0;
1329            int nout = 0;
1330            int wcount = 0;
1331            double e = 0;
1332            double t = 0;
1333            double z = 0;
1334            int i_ = 0;
1335            int i1_ = 0;
1336
1337            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1338            result = 0;
1339            for(i=0; i<=ssize-1; i++)
1340            {
1341               
1342                //
1343                // Process vector
1344                //
1345                for(i_=0; i_<=nin-1;i_++)
1346                {
1347                    network.x[i_] = xy[i,i_];
1348                }
1349                mlpprocess(ref network, ref network.x, ref network.y);
1350               
1351                //
1352                // Update error function
1353                //
1354                if( network.structinfo[6]==0 )
1355                {
1356                   
1357                    //
1358                    // Least squares error function
1359                    //
1360                    i1_ = (nin) - (0);
1361                    for(i_=0; i_<=nout-1;i_++)
1362                    {
1363                        network.y[i_] = network.y[i_] - xy[i,i_+i1_];
1364                    }
1365                    e = 0.0;
1366                    for(i_=0; i_<=nout-1;i_++)
1367                    {
1368                        e += network.y[i_]*network.y[i_];
1369                    }
1370                    result = result+e/2;
1371                }
1372                else
1373                {
1374                   
1375                    //
1376                    // Cross-entropy error function
1377                    //
1378                    k = (int)Math.Round(xy[i,nin]);
1379                    if( k>=0 & k<nout )
1380                    {
1381                        result = result+safecrossentropy(1, network.y[k]);
1382                    }
1383                }
1384            }
1385            return result;
1386        }
1387
1388
1389        /*************************************************************************
1390        Classification error
1391
1392          -- ALGLIB --
1393             Copyright 04.11.2007 by Bochkanov Sergey
1394        *************************************************************************/
1395        public static int mlpclserror(ref multilayerperceptron network,
1396            ref double[,] xy,
1397            int ssize)
1398        {
1399            int result = 0;
1400            int i = 0;
1401            int j = 0;
1402            int nin = 0;
1403            int nout = 0;
1404            int wcount = 0;
1405            double[] workx = new double[0];
1406            double[] worky = new double[0];
1407            int nn = 0;
1408            int ns = 0;
1409            int nmax = 0;
1410            int i_ = 0;
1411
1412            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1413            workx = new double[nin-1+1];
1414            worky = new double[nout-1+1];
1415            result = 0;
1416            for(i=0; i<=ssize-1; i++)
1417            {
1418               
1419                //
1420                // Process
1421                //
1422                for(i_=0; i_<=nin-1;i_++)
1423                {
1424                    workx[i_] = xy[i,i_];
1425                }
1426                mlpprocess(ref network, ref workx, ref worky);
1427               
1428                //
1429                // Network version of the answer
1430                //
1431                nmax = 0;
1432                for(j=0; j<=nout-1; j++)
1433                {
1434                    if( (double)(worky[j])>(double)(worky[nmax]) )
1435                    {
1436                        nmax = j;
1437                    }
1438                }
1439                nn = nmax;
1440               
1441                //
1442                // Right answer
1443                //
1444                if( mlpissoftmax(ref network) )
1445                {
1446                    ns = (int)Math.Round(xy[i,nin]);
1447                }
1448                else
1449                {
1450                    nmax = 0;
1451                    for(j=0; j<=nout-1; j++)
1452                    {
1453                        if( (double)(xy[i,nin+j])>(double)(xy[i,nin+nmax]) )
1454                        {
1455                            nmax = j;
1456                        }
1457                    }
1458                    ns = nmax;
1459                }
1460               
1461                //
1462                // compare
1463                //
1464                if( nn!=ns )
1465                {
1466                    result = result+1;
1467                }
1468            }
1469            return result;
1470        }
1471
1472
1473        /*************************************************************************
1474        Relative classification error on the test set
1475
1476        INPUT PARAMETERS:
1477            Network -   network
1478            XY      -   test set
1479            NPoints -   test set size
1480
1481        RESULT:
1482            percent of incorrectly classified cases. Works both for
1483            classifier networks and general purpose networks used as
1484            classifiers.
1485
1486          -- ALGLIB --
1487             Copyright 25.12.2008 by Bochkanov Sergey
1488        *************************************************************************/
1489        public static double mlprelclserror(ref multilayerperceptron network,
1490            ref double[,] xy,
1491            int npoints)
1492        {
1493            double result = 0;
1494
1495            result = (double)(mlpclserror(ref network, ref xy, npoints))/(double)(npoints);
1496            return result;
1497        }
1498
1499
1500        /*************************************************************************
1501        Average cross-entropy (in bits per element) on the test set
1502
1503        INPUT PARAMETERS:
1504            Network -   neural network
1505            XY      -   test set
1506            NPoints -   test set size
1507
1508        RESULT:
1509            CrossEntropy/(NPoints*LN(2)).
1510            Zero if network solves regression task.
1511
1512          -- ALGLIB --
1513             Copyright 08.01.2009 by Bochkanov Sergey
1514        *************************************************************************/
1515        public static double mlpavgce(ref multilayerperceptron network,
1516            ref double[,] xy,
1517            int npoints)
1518        {
1519            double result = 0;
1520            int nin = 0;
1521            int nout = 0;
1522            int wcount = 0;
1523
1524            if( mlpissoftmax(ref network) )
1525            {
1526                mlpproperties(ref network, ref nin, ref nout, ref wcount);
1527                result = mlperrorn(ref network, ref xy, npoints)/(npoints*Math.Log(2));
1528            }
1529            else
1530            {
1531                result = 0;
1532            }
1533            return result;
1534        }
1535
1536
1537        /*************************************************************************
1538        RMS error on the test set
1539
1540        INPUT PARAMETERS:
1541            Network -   neural network
1542            XY      -   test set
1543            NPoints -   test set size
1544
1545        RESULT:
1546            root mean square error.
1547            Its meaning for regression task is obvious. As for
1548            classification task, RMS error means error when estimating posterior
1549            probabilities.
1550
1551          -- ALGLIB --
1552             Copyright 04.11.2007 by Bochkanov Sergey
1553        *************************************************************************/
1554        public static double mlprmserror(ref multilayerperceptron network,
1555            ref double[,] xy,
1556            int npoints)
1557        {
1558            double result = 0;
1559            int nin = 0;
1560            int nout = 0;
1561            int wcount = 0;
1562
1563            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1564            result = Math.Sqrt(2*mlperror(ref network, ref xy, npoints)/(npoints*nout));
1565            return result;
1566        }
1567
1568
1569        /*************************************************************************
1570        Average error on the test set
1571
1572        INPUT PARAMETERS:
1573            Network -   neural network
1574            XY      -   test set
1575            NPoints -   test set size
1576
1577        RESULT:
1578            Its meaning for regression task is obvious. As for
1579            classification task, it means average error when estimating posterior
1580            probabilities.
1581
1582          -- ALGLIB --
1583             Copyright 11.03.2008 by Bochkanov Sergey
1584        *************************************************************************/
1585        public static double mlpavgerror(ref multilayerperceptron network,
1586            ref double[,] xy,
1587            int npoints)
1588        {
1589            double result = 0;
1590            int i = 0;
1591            int j = 0;
1592            int k = 0;
1593            int nin = 0;
1594            int nout = 0;
1595            int wcount = 0;
1596            int i_ = 0;
1597
1598            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1599            result = 0;
1600            for(i=0; i<=npoints-1; i++)
1601            {
1602                for(i_=0; i_<=nin-1;i_++)
1603                {
1604                    network.x[i_] = xy[i,i_];
1605                }
1606                mlpprocess(ref network, ref network.x, ref network.y);
1607                if( mlpissoftmax(ref network) )
1608                {
1609                   
1610                    //
1611                    // class labels
1612                    //
1613                    k = (int)Math.Round(xy[i,nin]);
1614                    for(j=0; j<=nout-1; j++)
1615                    {
1616                        if( j==k )
1617                        {
1618                            result = result+Math.Abs(1-network.y[j]);
1619                        }
1620                        else
1621                        {
1622                            result = result+Math.Abs(network.y[j]);
1623                        }
1624                    }
1625                }
1626                else
1627                {
1628                   
1629                    //
1630                    // real outputs
1631                    //
1632                    for(j=0; j<=nout-1; j++)
1633                    {
1634                        result = result+Math.Abs(xy[i,nin+j]-network.y[j]);
1635                    }
1636                }
1637            }
1638            result = result/(npoints*nout);
1639            return result;
1640        }
1641
1642
1643        /*************************************************************************
1644        Average relative error on the test set
1645
1646        INPUT PARAMETERS:
1647            Network -   neural network
1648            XY      -   test set
1649            NPoints -   test set size
1650
1651        RESULT:
1652            Its meaning for regression task is obvious. As for
1653            classification task, it means average relative error when estimating
1654            posterior probability of belonging to the correct class.
1655
1656          -- ALGLIB --
1657             Copyright 11.03.2008 by Bochkanov Sergey
1658        *************************************************************************/
1659        public static double mlpavgrelerror(ref multilayerperceptron network,
1660            ref double[,] xy,
1661            int npoints)
1662        {
1663            double result = 0;
1664            int i = 0;
1665            int j = 0;
1666            int k = 0;
1667            int lk = 0;
1668            int nin = 0;
1669            int nout = 0;
1670            int wcount = 0;
1671            int i_ = 0;
1672
1673            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1674            result = 0;
1675            k = 0;
1676            for(i=0; i<=npoints-1; i++)
1677            {
1678                for(i_=0; i_<=nin-1;i_++)
1679                {
1680                    network.x[i_] = xy[i,i_];
1681                }
1682                mlpprocess(ref network, ref network.x, ref network.y);
1683                if( mlpissoftmax(ref network) )
1684                {
1685                   
1686                    //
1687                    // class labels
1688                    //
1689                    lk = (int)Math.Round(xy[i,nin]);
1690                    for(j=0; j<=nout-1; j++)
1691                    {
1692                        if( j==lk )
1693                        {
1694                            result = result+Math.Abs(1-network.y[j]);
1695                            k = k+1;
1696                        }
1697                    }
1698                }
1699                else
1700                {
1701                   
1702                    //
1703                    // real outputs
1704                    //
1705                    for(j=0; j<=nout-1; j++)
1706                    {
1707                        if( (double)(xy[i,nin+j])!=(double)(0) )
1708                        {
1709                            result = result+Math.Abs(xy[i,nin+j]-network.y[j])/Math.Abs(xy[i,nin+j]);
1710                            k = k+1;
1711                        }
1712                    }
1713                }
1714            }
1715            if( k!=0 )
1716            {
1717                result = result/k;
1718            }
1719            return result;
1720        }
1721
1722
1723        /*************************************************************************
1724        Gradient calculation. Internal subroutine.
1725
1726          -- ALGLIB --
1727             Copyright 04.11.2007 by Bochkanov Sergey
1728        *************************************************************************/
1729        public static void mlpgrad(ref multilayerperceptron network,
1730            ref double[] x,
1731            ref double[] desiredy,
1732            ref double e,
1733            ref double[] grad)
1734        {
1735            int i = 0;
1736            int nout = 0;
1737            int ntotal = 0;
1738
1739           
1740            //
1741            // Prepare dError/dOut, internal structures
1742            //
1743            mlpprocess(ref network, ref x, ref network.y);
1744            nout = network.structinfo[2];
1745            ntotal = network.structinfo[3];
1746            e = 0;
1747            for(i=0; i<=ntotal-1; i++)
1748            {
1749                network.derror[i] = 0;
1750            }
1751            for(i=0; i<=nout-1; i++)
1752            {
1753                network.derror[ntotal-nout+i] = network.y[i]-desiredy[i];
1754                e = e+AP.Math.Sqr(network.y[i]-desiredy[i])/2;
1755            }
1756           
1757            //
1758            // gradient
1759            //
1760            mlpinternalcalculategradient(ref network, ref network.neurons, ref network.weights, ref network.derror, ref grad, false);
1761        }
1762
1763
1764        /*************************************************************************
1765        Gradient calculation (natural error function). Internal subroutine.
1766
1767          -- ALGLIB --
1768             Copyright 04.11.2007 by Bochkanov Sergey
1769        *************************************************************************/
1770        public static void mlpgradn(ref multilayerperceptron network,
1771            ref double[] x,
1772            ref double[] desiredy,
1773            ref double e,
1774            ref double[] grad)
1775        {
1776            double s = 0;
1777            int i = 0;
1778            int nout = 0;
1779            int ntotal = 0;
1780
1781           
1782            //
1783            // Prepare dError/dOut, internal structures
1784            //
1785            mlpprocess(ref network, ref x, ref network.y);
1786            nout = network.structinfo[2];
1787            ntotal = network.structinfo[3];
1788            for(i=0; i<=ntotal-1; i++)
1789            {
1790                network.derror[i] = 0;
1791            }
1792            e = 0;
1793            if( network.structinfo[6]==0 )
1794            {
1795               
1796                //
1797                // Regression network, least squares
1798                //
1799                for(i=0; i<=nout-1; i++)
1800                {
1801                    network.derror[ntotal-nout+i] = network.y[i]-desiredy[i];
1802                    e = e+AP.Math.Sqr(network.y[i]-desiredy[i])/2;
1803                }
1804            }
1805            else
1806            {
1807               
1808                //
1809                // Classification network, cross-entropy
1810                //
1811                s = 0;
1812                for(i=0; i<=nout-1; i++)
1813                {
1814                    s = s+desiredy[i];
1815                }
1816                for(i=0; i<=nout-1; i++)
1817                {
1818                    network.derror[ntotal-nout+i] = s*network.y[i]-desiredy[i];
1819                    e = e+safecrossentropy(desiredy[i], network.y[i]);
1820                }
1821            }
1822           
1823            //
1824            // gradient
1825            //
1826            mlpinternalcalculategradient(ref network, ref network.neurons, ref network.weights, ref network.derror, ref grad, true);
1827        }
1828
1829
1830        /*************************************************************************
1831        Batch gradient calculation. Internal subroutine.
1832
1833          -- ALGLIB --
1834             Copyright 04.11.2007 by Bochkanov Sergey
1835        *************************************************************************/
1836        public static void mlpgradbatch(ref multilayerperceptron network,
1837            ref double[,] xy,
1838            int ssize,
1839            ref double e,
1840            ref double[] grad)
1841        {
1842            int i = 0;
1843            int nin = 0;
1844            int nout = 0;
1845            int wcount = 0;
1846
1847            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1848            for(i=0; i<=wcount-1; i++)
1849            {
1850                grad[i] = 0;
1851            }
1852            e = 0;
1853            i = 0;
1854            while( i<=ssize-1 )
1855            {
1856                mlpchunkedgradient(ref network, ref xy, i, Math.Min(ssize, i+chunksize)-i, ref e, ref grad, false);
1857                i = i+chunksize;
1858            }
1859        }
1860
1861
1862        /*************************************************************************
1863        Batch gradient calculation (natural error function). Internal subroutine.
1864
1865          -- ALGLIB --
1866             Copyright 04.11.2007 by Bochkanov Sergey
1867        *************************************************************************/
1868        public static void mlpgradnbatch(ref multilayerperceptron network,
1869            ref double[,] xy,
1870            int ssize,
1871            ref double e,
1872            ref double[] grad)
1873        {
1874            int i = 0;
1875            int nin = 0;
1876            int nout = 0;
1877            int wcount = 0;
1878
1879            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1880            for(i=0; i<=wcount-1; i++)
1881            {
1882                grad[i] = 0;
1883            }
1884            e = 0;
1885            i = 0;
1886            while( i<=ssize-1 )
1887            {
1888                mlpchunkedgradient(ref network, ref xy, i, Math.Min(ssize, i+chunksize)-i, ref e, ref grad, true);
1889                i = i+chunksize;
1890            }
1891        }
1892
1893
1894        /*************************************************************************
1895        Batch Hessian calculation (natural error function) using R-algorithm.
1896        Internal subroutine.
1897
1898          -- ALGLIB --
1899             Copyright 26.01.2008 by Bochkanov Sergey.
1900             
1901             Hessian calculation based on R-algorithm described in
1902             "Fast Exact Multiplication by the Hessian",
1903             B. A. Pearlmutter,
1904             Neural Computation, 1994.
1905        *************************************************************************/
1906        public static void mlphessiannbatch(ref multilayerperceptron network,
1907            ref double[,] xy,
1908            int ssize,
1909            ref double e,
1910            ref double[] grad,
1911            ref double[,] h)
1912        {
1913            mlphessianbatchinternal(ref network, ref xy, ssize, true, ref e, ref grad, ref h);
1914        }
1915
1916
1917        /*************************************************************************
1918        Batch Hessian calculation using R-algorithm.
1919        Internal subroutine.
1920
1921          -- ALGLIB --
1922             Copyright 26.01.2008 by Bochkanov Sergey.
1923
1924             Hessian calculation based on R-algorithm described in
1925             "Fast Exact Multiplication by the Hessian",
1926             B. A. Pearlmutter,
1927             Neural Computation, 1994.
1928        *************************************************************************/
1929        public static void mlphessianbatch(ref multilayerperceptron network,
1930            ref double[,] xy,
1931            int ssize,
1932            ref double e,
1933            ref double[] grad,
1934            ref double[,] h)
1935        {
1936            mlphessianbatchinternal(ref network, ref xy, ssize, false, ref e, ref grad, ref h);
1937        }
1938
1939
1940        /*************************************************************************
1941        Internal subroutine, shouldn't be called by user.
1942        *************************************************************************/
1943        public static void mlpinternalprocessvector(ref int[] structinfo,
1944            ref double[] weights,
1945            ref double[] columnmeans,
1946            ref double[] columnsigmas,
1947            ref double[] neurons,
1948            ref double[] dfdnet,
1949            ref double[] x,
1950            ref double[] y)
1951        {
1952            int i = 0;
1953            int j = 0;
1954            int n1 = 0;
1955            int n2 = 0;
1956            int w1 = 0;
1957            int w2 = 0;
1958            int ntotal = 0;
1959            int nin = 0;
1960            int nout = 0;
1961            int istart = 0;
1962            int offs = 0;
1963            double net = 0;
1964            double e = 0;
1965            double f = 0;
1966            double df = 0;
1967            double d2f = 0;
1968            double mx = 0;
1969            bool perr = new bool();
1970            int i_ = 0;
1971            int i1_ = 0;
1972
1973           
1974            //
1975            // Read network geometry
1976            //
1977            nin = structinfo[1];
1978            nout = structinfo[2];
1979            ntotal = structinfo[3];
1980            istart = structinfo[5];
1981           
1982            //
1983            // Inputs standartisation and putting in the network
1984            //
1985            for(i=0; i<=nin-1; i++)
1986            {
1987                if( (double)(columnsigmas[i])!=(double)(0) )
1988                {
1989                    neurons[i] = (x[i]-columnmeans[i])/columnsigmas[i];
1990                }
1991                else
1992                {
1993                    neurons[i] = x[i]-columnmeans[i];
1994                }
1995            }
1996           
1997            //
1998            // Process network
1999            //
2000            for(i=0; i<=ntotal-1; i++)
2001            {
2002                offs = istart+i*nfieldwidth;
2003                if( structinfo[offs+0]>0 )
2004                {
2005                   
2006                    //
2007                    // Activation function
2008                    //
2009                    mlpactivationfunction(neurons[structinfo[offs+2]], structinfo[offs+0], ref f, ref df, ref d2f);
2010                    neurons[i] = f;
2011                    dfdnet[i] = df;
2012                }
2013                if( structinfo[offs+0]==0 )
2014                {
2015                   
2016                    //
2017                    // Adaptive summator
2018                    //
2019                    n1 = structinfo[offs+2];
2020                    n2 = n1+structinfo[offs+1]-1;
2021                    w1 = structinfo[offs+3];
2022                    w2 = w1+structinfo[offs+1]-1;
2023                    i1_ = (n1)-(w1);
2024                    net = 0.0;
2025                    for(i_=w1; i_<=w2;i_++)
2026                    {
2027                        net += weights[i_]*neurons[i_+i1_];
2028                    }
2029                    neurons[i] = net;
2030                    dfdnet[i] = 1.0;
2031                }
2032                if( structinfo[offs+0]<0 )
2033                {
2034                    perr = true;
2035                    if( structinfo[offs+0]==-2 )
2036                    {
2037                       
2038                        //
2039                        // input neuron, left unchanged
2040                        //
2041                        perr = false;
2042                    }
2043                    if( structinfo[offs+0]==-3 )
2044                    {
2045                       
2046                        //
2047                        // "-1" neuron
2048                        //
2049                        neurons[i] = -1;
2050                        perr = false;
2051                    }
2052                    if( structinfo[offs+0]==-4 )
2053                    {
2054                       
2055                        //
2056                        // "0" neuron
2057                        //
2058                        neurons[i] = 0;
2059                        perr = false;
2060                    }
2061                    System.Diagnostics.Debug.Assert(!perr, "MLPInternalProcessVector: internal error - unknown neuron type!");
2062                }
2063            }
2064           
2065            //
2066            // Extract result
2067            //
2068            i1_ = (ntotal-nout) - (0);
2069            for(i_=0; i_<=nout-1;i_++)
2070            {
2071                y[i_] = neurons[i_+i1_];
2072            }
2073           
2074            //
2075            // Softmax post-processing or standardisation if needed
2076            //
2077            System.Diagnostics.Debug.Assert(structinfo[6]==0 | structinfo[6]==1, "MLPInternalProcessVector: unknown normalization type!");
2078            if( structinfo[6]==1 )
2079            {
2080               
2081                //
2082                // Softmax
2083                //
2084                mx = y[0];
2085                for(i=1; i<=nout-1; i++)
2086                {
2087                    mx = Math.Max(mx, y[i]);
2088                }
2089                net = 0;
2090                for(i=0; i<=nout-1; i++)
2091                {
2092                    y[i] = Math.Exp(y[i]-mx);
2093                    net = net+y[i];
2094                }
2095                for(i=0; i<=nout-1; i++)
2096                {
2097                    y[i] = y[i]/net;
2098                }
2099            }
2100            else
2101            {
2102               
2103                //
2104                // Standardisation
2105                //
2106                for(i=0; i<=nout-1; i++)
2107                {
2108                    y[i] = y[i]*columnsigmas[nin+i]+columnmeans[nin+i];
2109                }
2110            }
2111        }
2112
2113
2114        /*************************************************************************
2115        Internal subroutine: adding new input layer to network
2116        *************************************************************************/
2117        private static void addinputlayer(int ncount,
2118            ref int[] lsizes,
2119            ref int[] ltypes,
2120            ref int[] lconnfirst,
2121            ref int[] lconnlast,
2122            ref int lastproc)
2123        {
2124            lsizes[0] = ncount;
2125            ltypes[0] = -2;
2126            lconnfirst[0] = 0;
2127            lconnlast[0] = 0;
2128            lastproc = 0;
2129        }
2130
2131
2132        /*************************************************************************
2133        Internal subroutine: adding new summator layer to network
2134        *************************************************************************/
2135        private static void addbiasedsummatorlayer(int ncount,
2136            ref int[] lsizes,
2137            ref int[] ltypes,
2138            ref int[] lconnfirst,
2139            ref int[] lconnlast,
2140            ref int lastproc)
2141        {
2142            lsizes[lastproc+1] = 1;
2143            ltypes[lastproc+1] = -3;
2144            lconnfirst[lastproc+1] = 0;
2145            lconnlast[lastproc+1] = 0;
2146            lsizes[lastproc+2] = ncount;
2147            ltypes[lastproc+2] = 0;
2148            lconnfirst[lastproc+2] = lastproc;
2149            lconnlast[lastproc+2] = lastproc+1;
2150            lastproc = lastproc+2;
2151        }
2152
2153
2154        /*************************************************************************
2155        Internal subroutine: adding new summator layer to network
2156        *************************************************************************/
2157        private static void addactivationlayer(int functype,
2158            ref int[] lsizes,
2159            ref int[] ltypes,
2160            ref int[] lconnfirst,
2161            ref int[] lconnlast,
2162            ref int lastproc)
2163        {
2164            System.Diagnostics.Debug.Assert(functype>0, "AddActivationLayer: incorrect function type");
2165            lsizes[lastproc+1] = lsizes[lastproc];
2166            ltypes[lastproc+1] = functype;
2167            lconnfirst[lastproc+1] = lastproc;
2168            lconnlast[lastproc+1] = lastproc;
2169            lastproc = lastproc+1;
2170        }
2171
2172
2173        /*************************************************************************
2174        Internal subroutine: adding new zero layer to network
2175        *************************************************************************/
2176        private static void addzerolayer(ref int[] lsizes,
2177            ref int[] ltypes,
2178            ref int[] lconnfirst,
2179            ref int[] lconnlast,
2180            ref int lastproc)
2181        {
2182            lsizes[lastproc+1] = 1;
2183            ltypes[lastproc+1] = -4;
2184            lconnfirst[lastproc+1] = 0;
2185            lconnlast[lastproc+1] = 0;
2186            lastproc = lastproc+1;
2187        }
2188
2189
2190        /*************************************************************************
2191        Internal subroutine.
2192
2193          -- ALGLIB --
2194             Copyright 04.11.2007 by Bochkanov Sergey
2195        *************************************************************************/
2196        private static void mlpcreate(int nin,
2197            int nout,
2198            ref int[] lsizes,
2199            ref int[] ltypes,
2200            ref int[] lconnfirst,
2201            ref int[] lconnlast,
2202            int layerscount,
2203            bool isclsnet,
2204            ref multilayerperceptron network)
2205        {
2206            int i = 0;
2207            int j = 0;
2208            int ssize = 0;
2209            int ntotal = 0;
2210            int wcount = 0;
2211            int offs = 0;
2212            int nprocessed = 0;
2213            int wallocated = 0;
2214            int[] localtemp = new int[0];
2215            int[] lnfirst = new int[0];
2216            int[] lnsyn = new int[0];
2217
2218           
2219            //
2220            // Check
2221            //
2222            System.Diagnostics.Debug.Assert(layerscount>0, "MLPCreate: wrong parameters!");
2223            System.Diagnostics.Debug.Assert(ltypes[0]==-2, "MLPCreate: wrong LTypes[0] (must be -2)!");
2224            for(i=0; i<=layerscount-1; i++)
2225            {
2226                System.Diagnostics.Debug.Assert(lsizes[i]>0, "MLPCreate: wrong LSizes!");
2227                System.Diagnostics.Debug.Assert(lconnfirst[i]>=0 & (lconnfirst[i]<i | i==0), "MLPCreate: wrong LConnFirst!");
2228                System.Diagnostics.Debug.Assert(lconnlast[i]>=lconnfirst[i] & (lconnlast[i]<i | i==0), "MLPCreate: wrong LConnLast!");
2229            }
2230           
2231            //
2232            // Build network geometry
2233            //
2234            lnfirst = new int[layerscount-1+1];
2235            lnsyn = new int[layerscount-1+1];
2236            ntotal = 0;
2237            wcount = 0;
2238            for(i=0; i<=layerscount-1; i++)
2239            {
2240               
2241                //
2242                // Analyze connections.
2243                // This code must throw an assertion in case of unknown LTypes[I]
2244                //
2245                lnsyn[i] = -1;
2246                if( ltypes[i]>=0 )
2247                {
2248                    lnsyn[i] = 0;
2249                    for(j=lconnfirst[i]; j<=lconnlast[i]; j++)
2250                    {
2251                        lnsyn[i] = lnsyn[i]+lsizes[j];
2252                    }
2253                }
2254                else
2255                {
2256                    if( ltypes[i]==-2 | ltypes[i]==-3 | ltypes[i]==-4 )
2257                    {
2258                        lnsyn[i] = 0;
2259                    }
2260                }
2261                System.Diagnostics.Debug.Assert(lnsyn[i]>=0, "MLPCreate: internal error #0!");
2262               
2263                //
2264                // Other info
2265                //
2266                lnfirst[i] = ntotal;
2267                ntotal = ntotal+lsizes[i];
2268                if( ltypes[i]==0 )
2269                {
2270                    wcount = wcount+lnsyn[i]*lsizes[i];
2271                }
2272            }
2273            ssize = 7+ntotal*nfieldwidth;
2274           
2275            //
2276            // Allocate
2277            //
2278            network.structinfo = new int[ssize-1+1];
2279            network.weights = new double[wcount-1+1];
2280            if( isclsnet )
2281            {
2282                network.columnmeans = new double[nin-1+1];
2283                network.columnsigmas = new double[nin-1+1];
2284            }
2285            else
2286            {
2287                network.columnmeans = new double[nin+nout-1+1];
2288                network.columnsigmas = new double[nin+nout-1+1];
2289            }
2290            network.neurons = new double[ntotal-1+1];
2291            network.chunks = new double[3*ntotal+1, chunksize-1+1];
2292            network.nwbuf = new double[Math.Max(wcount, 2*nout)-1+1];
2293            network.dfdnet = new double[ntotal-1+1];
2294            network.x = new double[nin-1+1];
2295            network.y = new double[nout-1+1];
2296            network.derror = new double[ntotal-1+1];
2297           
2298            //
2299            // Fill structure: global info
2300            //
2301            network.structinfo[0] = ssize;
2302            network.structinfo[1] = nin;
2303            network.structinfo[2] = nout;
2304            network.structinfo[3] = ntotal;
2305            network.structinfo[4] = wcount;
2306            network.structinfo[5] = 7;
2307            if( isclsnet )
2308            {
2309                network.structinfo[6] = 1;
2310            }
2311            else
2312            {
2313                network.structinfo[6] = 0;
2314            }
2315           
2316            //
2317            // Fill structure: neuron connections
2318            //
2319            nprocessed = 0;
2320            wallocated = 0;
2321            for(i=0; i<=layerscount-1; i++)
2322            {
2323                for(j=0; j<=lsizes[i]-1; j++)
2324                {
2325                    offs = network.structinfo[5]+nprocessed*nfieldwidth;
2326                    network.structinfo[offs+0] = ltypes[i];
2327                    if( ltypes[i]==0 )
2328                    {
2329                       
2330                        //
2331                        // Adaptive summator:
2332                        // * connections with weights to previous neurons
2333                        //
2334                        network.structinfo[offs+1] = lnsyn[i];
2335                        network.structinfo[offs+2] = lnfirst[lconnfirst[i]];
2336                        network.structinfo[offs+3] = wallocated;
2337                        wallocated = wallocated+lnsyn[i];
2338                        nprocessed = nprocessed+1;
2339                    }
2340                    if( ltypes[i]>0 )
2341                    {
2342                       
2343                        //
2344                        // Activation layer:
2345                        // * each neuron connected to one (only one) of previous neurons.
2346                        // * no weights
2347                        //
2348                        network.structinfo[offs+1] = 1;
2349                        network.structinfo[offs+2] = lnfirst[lconnfirst[i]]+j;
2350                        network.structinfo[offs+3] = -1;
2351                        nprocessed = nprocessed+1;
2352                    }
2353                    if( ltypes[i]==-2 | ltypes[i]==-3 | ltypes[i]==-4 )
2354                    {
2355                        nprocessed = nprocessed+1;
2356                    }
2357                }
2358            }
2359            System.Diagnostics.Debug.Assert(wallocated==wcount, "MLPCreate: internal error #1!");
2360            System.Diagnostics.Debug.Assert(nprocessed==ntotal, "MLPCreate: internal error #2!");
2361           
2362            //
2363            // Fill weights by small random values
2364            // Initialize means and sigmas
2365            //
2366            for(i=0; i<=wcount-1; i++)
2367            {
2368                network.weights[i] = AP.Math.RandomReal()-0.5;
2369            }
2370            for(i=0; i<=nin-1; i++)
2371            {
2372                network.columnmeans[i] = 0;
2373                network.columnsigmas[i] = 1;
2374            }
2375            if( !isclsnet )
2376            {
2377                for(i=0; i<=nout-1; i++)
2378                {
2379                    network.columnmeans[nin+i] = 0;
2380                    network.columnsigmas[nin+i] = 1;
2381                }
2382            }
2383        }
2384
2385
2386        /*************************************************************************
2387        Internal subroutine
2388
2389          -- ALGLIB --
2390             Copyright 04.11.2007 by Bochkanov Sergey
2391        *************************************************************************/
2392        private static void mlpactivationfunction(double net,
2393            int k,
2394            ref double f,
2395            ref double df,
2396            ref double d2f)
2397        {
2398            double net2 = 0;
2399            double arg = 0;
2400            double root = 0;
2401            double r = 0;
2402
2403            f = 0;
2404            df = 0;
2405            if( k==1 )
2406            {
2407               
2408                //
2409                // TanH activation function
2410                //
2411                if( (double)(Math.Abs(net))<(double)(100) )
2412                {
2413                    f = Math.Tanh(net);
2414                }
2415                else
2416                {
2417                    f = Math.Sign(net);
2418                }
2419                df = 1-AP.Math.Sqr(f);
2420                d2f = -(2*f*df);
2421                return;
2422            }
2423            if( k==3 )
2424            {
2425               
2426                //
2427                // EX activation function
2428                //
2429                if( (double)(net)>=(double)(0) )
2430                {
2431                    net2 = net*net;
2432                    arg = net2+1;
2433                    root = Math.Sqrt(arg);
2434                    f = net+root;
2435                    r = net/root;
2436                    df = 1+r;
2437                    d2f = (root-net*r)/arg;
2438                }
2439                else
2440                {
2441                    f = Math.Exp(net);
2442                    df = f;
2443                    d2f = f;
2444                }
2445                return;
2446            }
2447            if( k==2 )
2448            {
2449                f = Math.Exp(-AP.Math.Sqr(net));
2450                df = -(2*net*f);
2451                d2f = -(2*(f+df*net));
2452                return;
2453            }
2454        }
2455
2456
2457        /*************************************************************************
2458        Internal subroutine for Hessian calculation.
2459
2460        WARNING!!! Unspeakable math far beyong human capabilities :)
2461        *************************************************************************/
2462        private static void mlphessianbatchinternal(ref multilayerperceptron network,
2463            ref double[,] xy,
2464            int ssize,
2465            bool naturalerr,
2466            ref double e,
2467            ref double[] grad,
2468            ref double[,] h)
2469        {
2470            int nin = 0;
2471            int nout = 0;
2472            int wcount = 0;
2473            int ntotal = 0;
2474            int istart = 0;
2475            int i = 0;
2476            int j = 0;
2477            int k = 0;
2478            int kl = 0;
2479            int offs = 0;
2480            int n1 = 0;
2481            int n2 = 0;
2482            int w1 = 0;
2483            int w2 = 0;
2484            double s = 0;
2485            double t = 0;
2486            double v = 0;
2487            double et = 0;
2488            bool bflag = new bool();
2489            double f = 0;
2490            double df = 0;
2491            double d2f = 0;
2492            double deidyj = 0;
2493            double mx = 0;
2494            double net = 0;
2495            double q = 0;
2496            double z = 0;
2497            double s2 = 0;
2498            double expi = 0;
2499            double expj = 0;
2500            double[] x = new double[0];
2501            double[] desiredy = new double[0];
2502            double[] gt = new double[0];
2503            double[] zeros = new double[0];
2504            double[,] rx = new double[0,0];
2505            double[,] ry = new double[0,0];
2506            double[,] rdx = new double[0,0];
2507            double[,] rdy = new double[0,0];
2508            int i_ = 0;
2509            int i1_ = 0;
2510
2511            mlpproperties(ref network, ref nin, ref nout, ref wcount);
2512            ntotal = network.structinfo[3];
2513            istart = network.structinfo[5];
2514           
2515            //
2516            // Prepare
2517            //
2518            x = new double[nin-1+1];
2519            desiredy = new double[nout-1+1];
2520            zeros = new double[wcount-1+1];
2521            gt = new double[wcount-1+1];
2522            rx = new double[ntotal+nout-1+1, wcount-1+1];
2523            ry = new double[ntotal+nout-1+1, wcount-1+1];
2524            rdx = new double[ntotal+nout-1+1, wcount-1+1];
2525            rdy = new double[ntotal+nout-1+1, wcount-1+1];
2526            e = 0;
2527            for(i=0; i<=wcount-1; i++)
2528            {
2529                zeros[i] = 0;
2530            }
2531            for(i_=0; i_<=wcount-1;i_++)
2532            {
2533                grad[i_] = zeros[i_];
2534            }
2535            for(i=0; i<=wcount-1; i++)
2536            {
2537                for(i_=0; i_<=wcount-1;i_++)
2538                {
2539                    h[i,i_] = zeros[i_];
2540                }
2541            }
2542           
2543            //
2544            // Process
2545            //
2546            for(k=0; k<=ssize-1; k++)
2547            {
2548               
2549                //
2550                // Process vector with MLPGradN.
2551                // Now Neurons, DFDNET and DError contains results of the last run.
2552                //
2553                for(i_=0; i_<=nin-1;i_++)
2554                {
2555                    x[i_] = xy[k,i_];
2556                }
2557                if( mlpissoftmax(ref network) )
2558                {
2559                   
2560                    //
2561                    // class labels outputs
2562                    //
2563                    kl = (int)Math.Round(xy[k,nin]);
2564                    for(i=0; i<=nout-1; i++)
2565                    {
2566                        if( i==kl )
2567                        {
2568                            desiredy[i] = 1;
2569                        }
2570                        else
2571                        {
2572                            desiredy[i] = 0;
2573                        }
2574                    }
2575                }
2576                else
2577                {
2578                   
2579                    //
2580                    // real outputs
2581                    //
2582                    i1_ = (nin) - (0);
2583                    for(i_=0; i_<=nout-1;i_++)
2584                    {
2585                        desiredy[i_] = xy[k,i_+i1_];
2586                    }
2587                }
2588                if( naturalerr )
2589                {
2590                    mlpgradn(ref network, ref x, ref desiredy, ref et, ref gt);
2591                }
2592                else
2593                {
2594                    mlpgrad(ref network, ref x, ref desiredy, ref et, ref gt);
2595                }
2596               
2597                //
2598                // grad, error
2599                //
2600                e = e+et;
2601                for(i_=0; i_<=wcount-1;i_++)
2602                {
2603                    grad[i_] = grad[i_] + gt[i_];
2604                }
2605               
2606                //
2607                // Hessian.
2608                // Forward pass of the R-algorithm
2609                //
2610                for(i=0; i<=ntotal-1; i++)
2611                {
2612                    offs = istart+i*nfieldwidth;
2613                    for(i_=0; i_<=wcount-1;i_++)
2614                    {
2615                        rx[i,i_] = zeros[i_];
2616                    }
2617                    for(i_=0; i_<=wcount-1;i_++)
2618                    {
2619                        ry[i,i_] = zeros[i_];
2620                    }
2621                    if( network.structinfo[offs+0]>0 )
2622                    {
2623                       
2624                        //
2625                        // Activation function
2626                        //
2627                        n1 = network.structinfo[offs+2];
2628                        for(i_=0; i_<=wcount-1;i_++)
2629                        {
2630                            rx[i,i_] = ry[n1,i_];
2631                        }
2632                        v = network.dfdnet[i];
2633                        for(i_=0; i_<=wcount-1;i_++)
2634                        {
2635                            ry[i,i_] = v*rx[i,i_];
2636                        }
2637                    }
2638                    if( network.structinfo[offs+0]==0 )
2639                    {
2640                       
2641                        //
2642                        // Adaptive summator
2643                        //
2644                        n1 = network.structinfo[offs+2];
2645                        n2 = n1+network.structinfo[offs+1]-1;
2646                        w1 = network.structinfo[offs+3];
2647                        w2 = w1+network.structinfo[offs+1]-1;
2648                        for(j=n1; j<=n2; j++)
2649                        {
2650                            v = network.weights[w1+j-n1];
2651                            for(i_=0; i_<=wcount-1;i_++)
2652                            {
2653                                rx[i,i_] = rx[i,i_] + v*ry[j,i_];
2654                            }
2655                            rx[i,w1+j-n1] = rx[i,w1+j-n1]+network.neurons[j];
2656                        }
2657                        for(i_=0; i_<=wcount-1;i_++)
2658                        {
2659                            ry[i,i_] = rx[i,i_];
2660                        }
2661                    }
2662                    if( network.structinfo[offs+0]<0 )
2663                    {
2664                        bflag = true;
2665                        if( network.structinfo[offs+0]==-2 )
2666                        {
2667                           
2668                            //
2669                            // input neuron, left unchanged
2670                            //
2671                            bflag = false;
2672                        }
2673                        if( network.structinfo[offs+0]==-3 )
2674                        {
2675                           
2676                            //
2677                            // "-1" neuron, left unchanged
2678                            //
2679                            bflag = false;
2680                        }
2681                        if( network.structinfo[offs+0]==-4 )
2682                        {
2683                           
2684                            //
2685                            // "0" neuron, left unchanged
2686                            //
2687                            bflag = false;
2688                        }
2689                        System.Diagnostics.Debug.Assert(!bflag, "MLPHessianNBatch: internal error - unknown neuron type!");
2690                    }
2691                }
2692               
2693                //
2694                // Hessian. Backward pass of the R-algorithm.
2695                //
2696                // Stage 1. Initialize RDY
2697                //
2698                for(i=0; i<=ntotal+nout-1; i++)
2699                {
2700                    for(i_=0; i_<=wcount-1;i_++)
2701                    {
2702                        rdy[i,i_] = zeros[i_];
2703                    }
2704                }
2705                if( network.structinfo[6]==0 )
2706                {
2707                   
2708                    //
2709                    // Standardisation.
2710                    //
2711                    // In context of the Hessian calculation standardisation
2712                    // is considered as additional layer with weightless
2713                    // activation function:
2714                    //
2715                    // F(NET) := Sigma*NET
2716                    //
2717                    // So we add one more layer to forward pass, and
2718                    // make forward/backward pass through this layer.
2719                    //
2720                    for(i=0; i<=nout-1; i++)
2721                    {
2722                        n1 = ntotal-nout+i;
2723                        n2 = ntotal+i;
2724                       
2725                        //
2726                        // Forward pass from N1 to N2
2727                        //
2728                        for(i_=0; i_<=wcount-1;i_++)
2729                        {
2730                            rx[n2,i_] = ry[n1,i_];
2731                        }
2732                        v = network.columnsigmas[nin+i];
2733                        for(i_=0; i_<=wcount-1;i_++)
2734                        {
2735                            ry[n2,i_] = v*rx[n2,i_];
2736                        }
2737                       
2738                        //
2739                        // Initialization of RDY
2740                        //
2741                        for(i_=0; i_<=wcount-1;i_++)
2742                        {
2743                            rdy[n2,i_] = ry[n2,i_];
2744                        }
2745                       
2746                        //
2747                        // Backward pass from N2 to N1:
2748                        // 1. Calculate R(dE/dX).
2749                        // 2. No R(dE/dWij) is needed since weight of activation neuron
2750                        //    is fixed to 1. So we can update R(dE/dY) for
2751                        //    the connected neuron (note that Vij=0, Wij=1)
2752                        //
2753                        df = network.columnsigmas[nin+i];
2754                        for(i_=0; i_<=wcount-1;i_++)
2755                        {
2756                            rdx[n2,i_] = df*rdy[n2,i_];
2757                        }
2758                        for(i_=0; i_<=wcount-1;i_++)
2759                        {
2760                            rdy[n1,i_] = rdy[n1,i_] + rdx[n2,i_];
2761                        }
2762                    }
2763                }
2764                else
2765                {
2766                   
2767                    //
2768                    // Softmax.
2769                    //
2770                    // Initialize RDY using generalized expression for ei'(yi)
2771                    // (see expression (9) from p. 5 of "Fast Exact Multiplication by the Hessian").
2772                    //
2773                    // When we are working with softmax network, generalized
2774                    // expression for ei'(yi) is used because softmax
2775                    // normalization leads to ei, which depends on all y's
2776                    //
2777                    if( naturalerr )
2778                    {
2779                       
2780                        //
2781                        // softmax + cross-entropy.
2782                        // We have:
2783                        //
2784                        // S = sum(exp(yk)),
2785                        // ei = sum(trn)*exp(yi)/S-trn_i
2786                        //
2787                        // j=i:   d(ei)/d(yj) = T*exp(yi)*(S-exp(yi))/S^2
2788                        // j<>i:  d(ei)/d(yj) = -T*exp(yi)*exp(yj)/S^2
2789                        //
2790                        t = 0;
2791                        for(i=0; i<=nout-1; i++)
2792                        {
2793                            t = t+desiredy[i];
2794                        }
2795                        mx = network.neurons[ntotal-nout];
2796                        for(i=0; i<=nout-1; i++)
2797                        {
2798                            mx = Math.Max(mx, network.neurons[ntotal-nout+i]);
2799                        }
2800                        s = 0;
2801                        for(i=0; i<=nout-1; i++)
2802                        {
2803                            network.nwbuf[i] = Math.Exp(network.neurons[ntotal-nout+i]-mx);
2804                            s = s+network.nwbuf[i];
2805                        }
2806                        for(i=0; i<=nout-1; i++)
2807                        {
2808                            for(j=0; j<=nout-1; j++)
2809                            {
2810                                if( j==i )
2811                                {
2812                                    deidyj = t*network.nwbuf[i]*(s-network.nwbuf[i])/AP.Math.Sqr(s);
2813                                    for(i_=0; i_<=wcount-1;i_++)
2814                                    {
2815                                        rdy[ntotal-nout+i,i_] = rdy[ntotal-nout+i,i_] + deidyj*ry[ntotal-nout+i,i_];
2816                                    }
2817                                }
2818                                else
2819                                {
2820                                    deidyj = -(t*network.nwbuf[i]*network.nwbuf[j]/AP.Math.Sqr(s));
2821                                    for(i_=0; i_<=wcount-1;i_++)
2822                                    {
2823                                        rdy[ntotal-nout+i,i_] = rdy[ntotal-nout+i,i_] + deidyj*ry[ntotal-nout+j,i_];
2824                                    }
2825                                }
2826                            }
2827                        }
2828                    }
2829                    else
2830                    {
2831                       
2832                        //
2833                        // For a softmax + squared error we have expression
2834                        // far beyond human imagination so we dont even try
2835                        // to comment on it. Just enjoy the code...
2836                        //
2837                        // P.S. That's why "natural error" is called "natural" -
2838                        // compact beatiful expressions, fast code....
2839                        //
2840                        mx = network.neurons[ntotal-nout];
2841                        for(i=0; i<=nout-1; i++)
2842                        {
2843                            mx = Math.Max(mx, network.neurons[ntotal-nout+i]);
2844                        }
2845                        s = 0;
2846                        s2 = 0;
2847                        for(i=0; i<=nout-1; i++)
2848                        {
2849                            network.nwbuf[i] = Math.Exp(network.neurons[ntotal-nout+i]-mx);
2850                            s = s+network.nwbuf[i];
2851                            s2 = s2+AP.Math.Sqr(network.nwbuf[i]);
2852                        }
2853                        q = 0;
2854                        for(i=0; i<=nout-1; i++)
2855                        {
2856                            q = q+(network.y[i]-desiredy[i])*network.nwbuf[i];
2857                        }
2858                        for(i=0; i<=nout-1; i++)
2859                        {
2860                            z = -q+(network.y[i]-desiredy[i])*s;
2861                            expi = network.nwbuf[i];
2862                            for(j=0; j<=nout-1; j++)
2863                            {
2864                                expj = network.nwbuf[j];
2865                                if( j==i )
2866                                {
2867                                    deidyj = expi/AP.Math.Sqr(s)*((z+expi)*(s-2*expi)/s+expi*s2/AP.Math.Sqr(s));
2868                                }
2869                                else
2870                                {
2871                                    deidyj = expi*expj/AP.Math.Sqr(s)*(s2/AP.Math.Sqr(s)-2*z/s-(expi+expj)/s+(network.y[i]-desiredy[i])-(network.y[j]-desiredy[j]));
2872                                }
2873                                for(i_=0; i_<=wcount-1;i_++)
2874                                {
2875                                    rdy[ntotal-nout+i,i_] = rdy[ntotal-nout+i,i_] + deidyj*ry[ntotal-nout+j,i_];
2876                                }
2877                            }
2878                        }
2879                    }
2880                }
2881               
2882                //
2883                // Hessian. Backward pass of the R-algorithm
2884                //
2885                // Stage 2. Process.
2886                //
2887                for(i=ntotal-1; i>=0; i--)
2888                {
2889                   
2890                    //
2891                    // Possible variants:
2892                    // 1. Activation function
2893                    // 2. Adaptive summator
2894                    // 3. Special neuron
2895                    //
2896                    offs = istart+i*nfieldwidth;
2897                    if( network.structinfo[offs+0]>0 )
2898                    {
2899                        n1 = network.structinfo[offs+2];
2900                       
2901                        //
2902                        // First, calculate R(dE/dX).
2903                        //
2904                        mlpactivationfunction(network.neurons[n1], network.structinfo[offs+0], ref f, ref df, ref d2f);
2905                        v = d2f*network.derror[i];
2906                        for(i_=0; i_<=wcount-1;i_++)
2907                        {
2908                            rdx[i,i_] = df*rdy[i,i_];
2909                        }
2910                        for(i_=0; i_<=wcount-1;i_++)
2911                        {
2912                            rdx[i,i_] = rdx[i,i_] + v*rx[i,i_];
2913                        }
2914                       
2915                        //
2916                        // No R(dE/dWij) is needed since weight of activation neuron
2917                        // is fixed to 1.
2918                        //
2919                        // So we can update R(dE/dY) for the connected neuron.
2920                        // (note that Vij=0, Wij=1)
2921                        //
2922                        for(i_=0; i_<=wcount-1;i_++)
2923                        {
2924                            rdy[n1,i_] = rdy[n1,i_] + rdx[i,i_];
2925                        }
2926                    }
2927                    if( network.structinfo[offs+0]==0 )
2928                    {
2929                       
2930                        //
2931                        // Adaptive summator
2932                        //
2933                        n1 = network.structinfo[offs+2];
2934                        n2 = n1+network.structinfo[offs+1]-1;
2935                        w1 = network.structinfo[offs+3];
2936                        w2 = w1+network.structinfo[offs+1]-1;
2937                       
2938                        //
2939                        // First, calculate R(dE/dX).
2940                        //
2941                        for(i_=0; i_<=wcount-1;i_++)
2942                        {
2943                            rdx[i,i_] = rdy[i,i_];
2944                        }
2945                       
2946                        //
2947                        // Then, calculate R(dE/dWij)
2948                        //
2949                        for(j=w1; j<=w2; j++)
2950                        {
2951                            v = network.neurons[n1+j-w1];
2952                            for(i_=0; i_<=wcount-1;i_++)
2953                            {
2954                                h[j,i_] = h[j,i_] + v*rdx[i,i_];
2955                            }
2956                            v = network.derror[i];
2957                            for(i_=0; i_<=wcount-1;i_++)
2958                            {
2959                                h[j,i_] = h[j,i_] + v*ry[n1+j-w1,i_];
2960                            }
2961                        }
2962                       
2963                        //
2964                        // And finally, update R(dE/dY) for connected neurons.
2965                        //
2966                        for(j=w1; j<=w2; j++)
2967                        {
2968                            v = network.weights[j];
2969                            for(i_=0; i_<=wcount-1;i_++)
2970                            {
2971                                rdy[n1+j-w1,i_] = rdy[n1+j-w1,i_] + v*rdx[i,i_];
2972                            }
2973                            rdy[n1+j-w1,j] = rdy[n1+j-w1,j]+network.derror[i];
2974                        }
2975                    }
2976                    if( network.structinfo[offs+0]<0 )
2977                    {
2978                        bflag = false;
2979                        if( network.structinfo[offs+0]==-2 | network.structinfo[offs+0]==-3 | network.structinfo[offs+0]==-4 )
2980                        {
2981                           
2982                            //
2983                            // Special neuron type, no back-propagation required
2984                            //
2985                            bflag = true;
2986                        }
2987                        System.Diagnostics.Debug.Assert(bflag, "MLPHessianNBatch: unknown neuron type!");
2988                    }
2989                }
2990            }
2991        }
2992
2993
2994        /*************************************************************************
2995        Internal subroutine
2996
2997        Network must be processed by MLPProcess on X
2998        *************************************************************************/
2999        private static void mlpinternalcalculategradient(ref multilayerperceptron network,
3000            ref double[] neurons,
3001            ref double[] weights,
3002            ref double[] derror,
3003            ref double[] grad,
3004            bool naturalerrorfunc)
3005        {
3006            int i = 0;
3007            int j = 0;
3008            int n1 = 0;
3009            int n2 = 0;
3010            int w1 = 0;
3011            int w2 = 0;
3012            int ntotal = 0;
3013            int istart = 0;
3014            int nin = 0;
3015            int nout = 0;
3016            int offs = 0;
3017            double dedf = 0;
3018            double dfdnet = 0;
3019            double v = 0;
3020            double fown = 0;
3021            double deown = 0;
3022            double net = 0;
3023            double mx = 0;
3024            bool bflag = new bool();
3025            int i_ = 0;
3026            int i1_ = 0;
3027
3028           
3029            //
3030            // Read network geometry
3031            //
3032            nin = network.structinfo[1];
3033            nout = network.structinfo[2];
3034            ntotal = network.structinfo[3];
3035            istart = network.structinfo[5];
3036           
3037            //
3038            // Pre-processing of dError/dOut:
3039            // from dError/dOut(normalized) to dError/dOut(non-normalized)
3040            //
3041            System.Diagnostics.Debug.Assert(network.structinfo[6]==0 | network.structinfo[6]==1, "MLPInternalCalculateGradient: unknown normalization type!");
3042            if( network.structinfo[6]==1 )
3043            {
3044               
3045                //
3046                // Softmax
3047                //
3048                if( !naturalerrorfunc )
3049                {
3050                    mx = network.neurons[ntotal-nout];
3051                    for(i=0; i<=nout-1; i++)
3052                    {
3053                        mx = Math.Max(mx, network.neurons[ntotal-nout+i]);
3054                    }
3055                    net = 0;
3056                    for(i=0; i<=nout-1; i++)
3057                    {
3058                        network.nwbuf[i] = Math.Exp(network.neurons[ntotal-nout+i]-mx);
3059                        net = net+network.nwbuf[i];
3060                    }
3061                    i1_ = (0)-(ntotal-nout);
3062                    v = 0.0;
3063                    for(i_=ntotal-nout; i_<=ntotal-1;i_++)
3064                    {
3065                        v += network.derror[i_]*network.nwbuf[i_+i1_];
3066                    }
3067                    for(i=0; i<=nout-1; i++)
3068                    {
3069                        fown = network.nwbuf[i];
3070                        deown = network.derror[ntotal-nout+i];
3071                        network.nwbuf[nout+i] = (-v+deown*fown+deown*(net-fown))*fown/AP.Math.Sqr(net);
3072                    }
3073                    for(i=0; i<=nout-1; i++)
3074                    {
3075                        network.derror[ntotal-nout+i] = network.nwbuf[nout+i];
3076                    }
3077                }
3078            }
3079            else
3080            {
3081               
3082                //
3083                // Un-standardisation
3084                //
3085                for(i=0; i<=nout-1; i++)
3086                {
3087                    network.derror[ntotal-nout+i] = network.derror[ntotal-nout+i]*network.columnsigmas[nin+i];
3088                }
3089            }
3090           
3091            //
3092            // Backpropagation
3093            //
3094            for(i=ntotal-1; i>=0; i--)
3095            {
3096               
3097                //
3098                // Extract info
3099                //
3100                offs = istart+i*nfieldwidth;
3101                if( network.structinfo[offs+0]>0 )
3102                {
3103                   
3104                    //
3105                    // Activation function
3106                    //
3107                    dedf = network.derror[i];
3108                    dfdnet = network.dfdnet[i];
3109                    derror[network.structinfo[offs+2]] = derror[network.structinfo[offs+2]]+dedf*dfdnet;
3110                }
3111                if( network.structinfo[offs+0]==0 )
3112                {
3113                   
3114                    //
3115                    // Adaptive summator
3116                    //
3117                    n1 = network.structinfo[offs+2];
3118                    n2 = n1+network.structinfo[offs+1]-1;
3119                    w1 = network.structinfo[offs+3];
3120                    w2 = w1+network.structinfo[offs+1]-1;
3121                    dedf = network.derror[i];
3122                    dfdnet = 1.0;
3123                    v = dedf*dfdnet;
3124                    i1_ = (n1) - (w1);
3125                    for(i_=w1; i_<=w2;i_++)
3126                    {
3127                        grad[i_] = v*neurons[i_+i1_];
3128                    }
3129                    i1_ = (w1) - (n1);
3130                    for(i_=n1; i_<=n2;i_++)
3131                    {
3132                        derror[i_] = derror[i_] + v*weights[i_+i1_];
3133                    }
3134                }
3135                if( network.structinfo[offs+0]<0 )
3136                {
3137                    bflag = false;
3138                    if( network.structinfo[offs+0]==-2 | network.structinfo[offs+0]==-3 | network.structinfo[offs+0]==-4 )
3139                    {
3140                       
3141                        //
3142                        // Special neuron type, no back-propagation required
3143                        //
3144                        bflag = true;
3145                    }
3146                    System.Diagnostics.Debug.Assert(bflag, "MLPInternalCalculateGradient: unknown neuron type!");
3147                }
3148            }
3149        }
3150
3151
3152        /*************************************************************************
3153        Internal subroutine, chunked gradient
3154        *************************************************************************/
3155        private static void mlpchunkedgradient(ref multilayerperceptron network,
3156            ref double[,] xy,
3157            int cstart,
3158            int csize,
3159            ref double e,
3160            ref double[] grad,
3161            bool naturalerrorfunc)
3162        {
3163            int i = 0;
3164            int j = 0;
3165            int k = 0;
3166            int kl = 0;
3167            int n1 = 0;
3168            int n2 = 0;
3169            int w1 = 0;
3170            int w2 = 0;
3171            int c1 = 0;
3172            int c2 = 0;
3173            int ntotal = 0;
3174            int nin = 0;
3175            int nout = 0;
3176            int offs = 0;
3177            double dedf = 0;
3178            double dfdnet = 0;
3179            double f = 0;
3180            double df = 0;
3181            double d2f = 0;
3182            double v = 0;
3183            double s = 0;
3184            double fown = 0;
3185            double deown = 0;
3186            double net = 0;
3187            double lnnet = 0;
3188            double mx = 0;
3189            bool bflag = new bool();
3190            int istart = 0;
3191            int ineurons = 0;
3192            int idfdnet = 0;
3193            int iderror = 0;
3194            int izeros = 0;
3195            int i_ = 0;
3196            int i1_ = 0;
3197
3198           
3199            //
3200            // Read network geometry, prepare data
3201            //
3202            nin = network.structinfo[1];
3203            nout = network.structinfo[2];
3204            ntotal = network.structinfo[3];
3205            istart = network.structinfo[5];
3206            c1 = cstart;
3207            c2 = cstart+csize-1;
3208            ineurons = 0;
3209            idfdnet = ntotal;
3210            iderror = 2*ntotal;
3211            izeros = 3*ntotal;
3212            for(j=0; j<=csize-1; j++)
3213            {
3214                network.chunks[izeros,j] = 0;
3215            }
3216           
3217            //
3218            // Forward pass:
3219            // 1. Load inputs from XY to Chunks[0:NIn-1,0:CSize-1]
3220            // 2. Forward pass
3221            //
3222            for(i=0; i<=nin-1; i++)
3223            {
3224                for(j=0; j<=csize-1; j++)
3225                {
3226                    if( (double)(network.columnsigmas[i])!=(double)(0) )
3227                    {
3228                        network.chunks[i,j] = (xy[c1+j,i]-network.columnmeans[i])/network.columnsigmas[i];
3229                    }
3230                    else
3231                    {
3232                        network.chunks[i,j] = xy[c1+j,i]-network.columnmeans[i];
3233                    }
3234                }
3235            }
3236            for(i=0; i<=ntotal-1; i++)
3237            {
3238                offs = istart+i*nfieldwidth;
3239                if( network.structinfo[offs+0]>0 )
3240                {
3241                   
3242                    //
3243                    // Activation function:
3244                    // * calculate F vector, F(i) = F(NET(i))
3245                    //
3246                    n1 = network.structinfo[offs+2];
3247                    for(i_=0; i_<=csize-1;i_++)
3248                    {
3249                        network.chunks[i,i_] = network.chunks[n1,i_];
3250                    }
3251                    for(j=0; j<=csize-1; j++)
3252                    {
3253                        mlpactivationfunction(network.chunks[i,j], network.structinfo[offs+0], ref f, ref df, ref d2f);
3254                        network.chunks[i,j] = f;
3255                        network.chunks[idfdnet+i,j] = df;
3256                    }
3257                }
3258                if( network.structinfo[offs+0]==0 )
3259                {
3260                   
3261                    //
3262                    // Adaptive summator:
3263                    // * calculate NET vector, NET(i) = SUM(W(j,i)*Neurons(j),j=N1..N2)
3264                    //
3265                    n1 = network.structinfo[offs+2];
3266                    n2 = n1+network.structinfo[offs+1]-1;
3267                    w1 = network.structinfo[offs+3];
3268                    w2 = w1+network.structinfo[offs+1]-1;
3269                    for(i_=0; i_<=csize-1;i_++)
3270                    {
3271                        network.chunks[i,i_] = network.chunks[izeros,i_];
3272                    }
3273                    for(j=n1; j<=n2; j++)
3274                    {
3275                        v = network.weights[w1+j-n1];
3276                        for(i_=0; i_<=csize-1;i_++)
3277                        {
3278                            network.chunks[i,i_] = network.chunks[i,i_] + v*network.chunks[j,i_];
3279                        }
3280                    }
3281                }
3282                if( network.structinfo[offs+0]<0 )
3283                {
3284                    bflag = false;
3285                    if( network.structinfo[offs+0]==-2 )
3286                    {
3287                       
3288                        //
3289                        // input neuron, left unchanged
3290                        //
3291                        bflag = true;
3292                    }
3293                    if( network.structinfo[offs+0]==-3 )
3294                    {
3295                       
3296                        //
3297                        // "-1" neuron
3298                        //
3299                        for(k=0; k<=csize-1; k++)
3300                        {
3301                            network.chunks[i,k] = -1;
3302                        }
3303                        bflag = true;
3304                    }
3305                    if( network.structinfo[offs+0]==-4 )
3306                    {
3307                       
3308                        //
3309                        // "0" neuron
3310                        //
3311                        for(k=0; k<=csize-1; k++)
3312                        {
3313                            network.chunks[i,k] = 0;
3314                        }
3315                        bflag = true;
3316                    }
3317                    System.Diagnostics.Debug.Assert(bflag, "MLPChunkedGradient: internal error - unknown neuron type!");
3318                }
3319            }
3320           
3321            //
3322            // Post-processing, error, dError/dOut
3323            //
3324            for(i=0; i<=ntotal-1; i++)
3325            {
3326                for(i_=0; i_<=csize-1;i_++)
3327                {
3328                    network.chunks[iderror+i,i_] = network.chunks[izeros,i_];
3329                }
3330            }
3331            System.Diagnostics.Debug.Assert(network.structinfo[6]==0 | network.structinfo[6]==1, "MLPChunkedGradient: unknown normalization type!");
3332            if( network.structinfo[6]==1 )
3333            {
3334               
3335                //
3336                // Softmax output, classification network.
3337                //
3338                // For each K = 0..CSize-1 do:
3339                // 1. place exp(outputs[k]) to NWBuf[0:NOut-1]
3340                // 2. place sum(exp(..)) to NET
3341                // 3. calculate dError/dOut and place it to the second block of Chunks
3342                //
3343                for(k=0; k<=csize-1; k++)
3344                {
3345                   
3346                    //
3347                    // Normalize
3348                    //
3349                    mx = network.chunks[ntotal-nout,k];
3350                    for(i=1; i<=nout-1; i++)
3351                    {
3352                        mx = Math.Max(mx, network.chunks[ntotal-nout+i,k]);
3353                    }
3354                    net = 0;
3355                    for(i=0; i<=nout-1; i++)
3356                    {
3357                        network.nwbuf[i] = Math.Exp(network.chunks[ntotal-nout+i,k]-mx);
3358                        net = net+network.nwbuf[i];
3359                    }
3360                   
3361                    //
3362                    // Calculate error function and dError/dOut
3363                    //
3364                    if( naturalerrorfunc )
3365                    {
3366                       
3367                        //
3368                        // Natural error func.
3369                        //
3370                        //
3371                        s = 1;
3372                        lnnet = Math.Log(net);
3373                        kl = (int)Math.Round(xy[cstart+k,nin]);
3374                        for(i=0; i<=nout-1; i++)
3375                        {
3376                            if( i==kl )
3377                            {
3378                                v = 1;
3379                            }
3380                            else
3381                            {
3382                                v = 0;
3383                            }
3384                            network.chunks[iderror+ntotal-nout+i,k] = s*network.nwbuf[i]/net-v;
3385                            e = e+safecrossentropy(v, network.nwbuf[i]/net);
3386                        }
3387                    }
3388                    else
3389                    {
3390                       
3391                        //
3392                        // Least squares error func
3393                        // Error, dError/dOut(normalized)
3394                        //
3395                        kl = (int)Math.Round(xy[cstart+k,nin]);
3396                        for(i=0; i<=nout-1; i++)
3397                        {
3398                            if( i==kl )
3399                            {
3400                                v = network.nwbuf[i]/net-1;
3401                            }
3402                            else
3403                            {
3404                                v = network.nwbuf[i]/net;
3405                            }
3406                            network.nwbuf[nout+i] = v;
3407                            e = e+AP.Math.Sqr(v)/2;
3408                        }
3409                       
3410                        //
3411                        // From dError/dOut(normalized) to dError/dOut(non-normalized)
3412                        //
3413                        i1_ = (0)-(nout);
3414                        v = 0.0;
3415                        for(i_=nout; i_<=2*nout-1;i_++)
3416                        {
3417                            v += network.nwbuf[i_]*network.nwbuf[i_+i1_];
3418                        }
3419                        for(i=0; i<=nout-1; i++)
3420                        {
3421                            fown = network.nwbuf[i];
3422                            deown = network.nwbuf[nout+i];
3423                            network.chunks[iderror+ntotal-nout+i,k] = (-v+deown*fown+deown*(net-fown))*fown/AP.Math.Sqr(net);
3424                        }
3425                    }
3426                }
3427            }
3428            else
3429            {
3430               
3431                //
3432                // Normal output, regression network
3433                //
3434                // For each K = 0..CSize-1 do:
3435                // 1. calculate dError/dOut and place it to the second block of Chunks
3436                //
3437                for(i=0; i<=nout-1; i++)
3438                {
3439                    for(j=0; j<=csize-1; j++)
3440                    {
3441                        v = network.chunks[ntotal-nout+i,j]*network.columnsigmas[nin+i]+network.columnmeans[nin+i]-xy[cstart+j,nin+i];
3442                        network.chunks[iderror+ntotal-nout+i,j] = v*network.columnsigmas[nin+i];
3443                        e = e+AP.Math.Sqr(v)/2;
3444                    }
3445                }
3446            }
3447           
3448            //
3449            // Backpropagation
3450            //
3451            for(i=ntotal-1; i>=0; i--)
3452            {
3453               
3454                //
3455                // Extract info
3456                //
3457                offs = istart+i*nfieldwidth;
3458                if( network.structinfo[offs+0]>0 )
3459                {
3460                   
3461                    //
3462                    // Activation function
3463                    //
3464                    n1 = network.structinfo[offs+2];
3465                    for(k=0; k<=csize-1; k++)
3466                    {
3467                        network.chunks[iderror+i,k] = network.chunks[iderror+i,k]*network.chunks[idfdnet+i,k];
3468                    }
3469                    for(i_=0; i_<=csize-1;i_++)
3470                    {
3471                        network.chunks[iderror+n1,i_] = network.chunks[iderror+n1,i_] + network.chunks[iderror+i,i_];
3472                    }
3473                }
3474                if( network.structinfo[offs+0]==0 )
3475                {
3476                   
3477                    //
3478                    // "Normal" activation function
3479                    //
3480                    n1 = network.structinfo[offs+2];
3481                    n2 = n1+network.structinfo[offs+1]-1;
3482                    w1 = network.structinfo[offs+3];
3483                    w2 = w1+network.structinfo[offs+1]-1;
3484                    for(j=w1; j<=w2; j++)
3485                    {
3486                        v = 0.0;
3487                        for(i_=0; i_<=csize-1;i_++)
3488                        {
3489                            v += network.chunks[n1+j-w1,i_]*network.chunks[iderror+i,i_];
3490                        }
3491                        grad[j] = grad[j]+v;
3492                    }
3493                    for(j=n1; j<=n2; j++)
3494                    {
3495                        v = network.weights[w1+j-n1];
3496                        for(i_=0; i_<=csize-1;i_++)
3497                        {
3498                            network.chunks[iderror+j,i_] = network.chunks[iderror+j,i_] + v*network.chunks[iderror+i,i_];
3499                        }
3500                    }
3501                }
3502                if( network.structinfo[offs+0]<0 )
3503                {
3504                    bflag = false;
3505                    if( network.structinfo[offs+0]==-2 | network.structinfo[offs+0]==-3 | network.structinfo[offs+0]==-4 )
3506                    {
3507                       
3508                        //
3509                        // Special neuron type, no back-propagation required
3510                        //
3511                        bflag = true;
3512                    }
3513                    System.Diagnostics.Debug.Assert(bflag, "MLPInternalCalculateGradient: unknown neuron type!");
3514                }
3515            }
3516        }
3517
3518
3519        /*************************************************************************
3520        Returns T*Ln(T/Z), guarded against overflow/underflow.
3521        Internal subroutine.
3522        *************************************************************************/
3523        private static double safecrossentropy(double t,
3524            double z)
3525        {
3526            double result = 0;
3527            double r = 0;
3528
3529            if( (double)(t)==(double)(0) )
3530            {
3531                result = 0;
3532            }
3533            else
3534            {
3535                if( (double)(Math.Abs(z))>(double)(1) )
3536                {
3537                   
3538                    //
3539                    // Shouldn't be the case with softmax,
3540                    // but we just want to be sure.
3541                    //
3542                    if( (double)(t/z)==(double)(0) )
3543                    {
3544                        r = AP.Math.MinRealNumber;
3545                    }
3546                    else
3547                    {
3548                        r = t/z;
3549                    }
3550                }
3551                else
3552                {
3553                   
3554                    //
3555                    // Normal case
3556                    //
3557                    if( (double)(z)==(double)(0) | (double)(Math.Abs(t))>=(double)(AP.Math.MaxRealNumber*Math.Abs(z)) )
3558                    {
3559                        r = AP.Math.MaxRealNumber;
3560                    }
3561                    else
3562                    {
3563                        r = t/z;
3564                    }
3565                }
3566                result = t*Math.Log(r);
3567            }
3568            return result;
3569        }
3570    }
3571}
Note: See TracBrowser for help on using the repository browser.