Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/mlpbase.cs @ 5955

Last change on this file since 5955 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 126.4 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 k = 0;
1327            int nin = 0;
1328            int nout = 0;
1329            int wcount = 0;
1330            double e = 0;
1331            int i_ = 0;
1332            int i1_ = 0;
1333
1334            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1335            result = 0;
1336            for(i=0; i<=ssize-1; i++)
1337            {
1338               
1339                //
1340                // Process vector
1341                //
1342                for(i_=0; i_<=nin-1;i_++)
1343                {
1344                    network.x[i_] = xy[i,i_];
1345                }
1346                mlpprocess(ref network, ref network.x, ref network.y);
1347               
1348                //
1349                // Update error function
1350                //
1351                if( network.structinfo[6]==0 )
1352                {
1353                   
1354                    //
1355                    // Least squares error function
1356                    //
1357                    i1_ = (nin) - (0);
1358                    for(i_=0; i_<=nout-1;i_++)
1359                    {
1360                        network.y[i_] = network.y[i_] - xy[i,i_+i1_];
1361                    }
1362                    e = 0.0;
1363                    for(i_=0; i_<=nout-1;i_++)
1364                    {
1365                        e += network.y[i_]*network.y[i_];
1366                    }
1367                    result = result+e/2;
1368                }
1369                else
1370                {
1371                   
1372                    //
1373                    // Cross-entropy error function
1374                    //
1375                    k = (int)Math.Round(xy[i,nin]);
1376                    if( k>=0 & k<nout )
1377                    {
1378                        result = result+safecrossentropy(1, network.y[k]);
1379                    }
1380                }
1381            }
1382            return result;
1383        }
1384
1385
1386        /*************************************************************************
1387        Classification error
1388
1389          -- ALGLIB --
1390             Copyright 04.11.2007 by Bochkanov Sergey
1391        *************************************************************************/
1392        public static int mlpclserror(ref multilayerperceptron network,
1393            ref double[,] xy,
1394            int ssize)
1395        {
1396            int result = 0;
1397            int i = 0;
1398            int j = 0;
1399            int nin = 0;
1400            int nout = 0;
1401            int wcount = 0;
1402            double[] workx = new double[0];
1403            double[] worky = new double[0];
1404            int nn = 0;
1405            int ns = 0;
1406            int nmax = 0;
1407            int i_ = 0;
1408
1409            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1410            workx = new double[nin-1+1];
1411            worky = new double[nout-1+1];
1412            result = 0;
1413            for(i=0; i<=ssize-1; i++)
1414            {
1415               
1416                //
1417                // Process
1418                //
1419                for(i_=0; i_<=nin-1;i_++)
1420                {
1421                    workx[i_] = xy[i,i_];
1422                }
1423                mlpprocess(ref network, ref workx, ref worky);
1424               
1425                //
1426                // Network version of the answer
1427                //
1428                nmax = 0;
1429                for(j=0; j<=nout-1; j++)
1430                {
1431                    if( (double)(worky[j])>(double)(worky[nmax]) )
1432                    {
1433                        nmax = j;
1434                    }
1435                }
1436                nn = nmax;
1437               
1438                //
1439                // Right answer
1440                //
1441                if( mlpissoftmax(ref network) )
1442                {
1443                    ns = (int)Math.Round(xy[i,nin]);
1444                }
1445                else
1446                {
1447                    nmax = 0;
1448                    for(j=0; j<=nout-1; j++)
1449                    {
1450                        if( (double)(xy[i,nin+j])>(double)(xy[i,nin+nmax]) )
1451                        {
1452                            nmax = j;
1453                        }
1454                    }
1455                    ns = nmax;
1456                }
1457               
1458                //
1459                // compare
1460                //
1461                if( nn!=ns )
1462                {
1463                    result = result+1;
1464                }
1465            }
1466            return result;
1467        }
1468
1469
1470        /*************************************************************************
1471        Relative classification error on the test set
1472
1473        INPUT PARAMETERS:
1474            Network -   network
1475            XY      -   test set
1476            NPoints -   test set size
1477
1478        RESULT:
1479            percent of incorrectly classified cases. Works both for
1480            classifier networks and general purpose networks used as
1481            classifiers.
1482
1483          -- ALGLIB --
1484             Copyright 25.12.2008 by Bochkanov Sergey
1485        *************************************************************************/
1486        public static double mlprelclserror(ref multilayerperceptron network,
1487            ref double[,] xy,
1488            int npoints)
1489        {
1490            double result = 0;
1491
1492            result = (double)(mlpclserror(ref network, ref xy, npoints))/(double)(npoints);
1493            return result;
1494        }
1495
1496
1497        /*************************************************************************
1498        Average cross-entropy (in bits per element) on the test set
1499
1500        INPUT PARAMETERS:
1501            Network -   neural network
1502            XY      -   test set
1503            NPoints -   test set size
1504
1505        RESULT:
1506            CrossEntropy/(NPoints*LN(2)).
1507            Zero if network solves regression task.
1508
1509          -- ALGLIB --
1510             Copyright 08.01.2009 by Bochkanov Sergey
1511        *************************************************************************/
1512        public static double mlpavgce(ref multilayerperceptron network,
1513            ref double[,] xy,
1514            int npoints)
1515        {
1516            double result = 0;
1517            int nin = 0;
1518            int nout = 0;
1519            int wcount = 0;
1520
1521            if( mlpissoftmax(ref network) )
1522            {
1523                mlpproperties(ref network, ref nin, ref nout, ref wcount);
1524                result = mlperrorn(ref network, ref xy, npoints)/(npoints*Math.Log(2));
1525            }
1526            else
1527            {
1528                result = 0;
1529            }
1530            return result;
1531        }
1532
1533
1534        /*************************************************************************
1535        RMS error on the test set
1536
1537        INPUT PARAMETERS:
1538            Network -   neural network
1539            XY      -   test set
1540            NPoints -   test set size
1541
1542        RESULT:
1543            root mean square error.
1544            Its meaning for regression task is obvious. As for
1545            classification task, RMS error means error when estimating posterior
1546            probabilities.
1547
1548          -- ALGLIB --
1549             Copyright 04.11.2007 by Bochkanov Sergey
1550        *************************************************************************/
1551        public static double mlprmserror(ref multilayerperceptron network,
1552            ref double[,] xy,
1553            int npoints)
1554        {
1555            double result = 0;
1556            int nin = 0;
1557            int nout = 0;
1558            int wcount = 0;
1559
1560            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1561            result = Math.Sqrt(2*mlperror(ref network, ref xy, npoints)/(npoints*nout));
1562            return result;
1563        }
1564
1565
1566        /*************************************************************************
1567        Average error on the test set
1568
1569        INPUT PARAMETERS:
1570            Network -   neural network
1571            XY      -   test set
1572            NPoints -   test set size
1573
1574        RESULT:
1575            Its meaning for regression task is obvious. As for
1576            classification task, it means average error when estimating posterior
1577            probabilities.
1578
1579          -- ALGLIB --
1580             Copyright 11.03.2008 by Bochkanov Sergey
1581        *************************************************************************/
1582        public static double mlpavgerror(ref multilayerperceptron network,
1583            ref double[,] xy,
1584            int npoints)
1585        {
1586            double result = 0;
1587            int i = 0;
1588            int j = 0;
1589            int k = 0;
1590            int nin = 0;
1591            int nout = 0;
1592            int wcount = 0;
1593            int i_ = 0;
1594
1595            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1596            result = 0;
1597            for(i=0; i<=npoints-1; i++)
1598            {
1599                for(i_=0; i_<=nin-1;i_++)
1600                {
1601                    network.x[i_] = xy[i,i_];
1602                }
1603                mlpprocess(ref network, ref network.x, ref network.y);
1604                if( mlpissoftmax(ref network) )
1605                {
1606                   
1607                    //
1608                    // class labels
1609                    //
1610                    k = (int)Math.Round(xy[i,nin]);
1611                    for(j=0; j<=nout-1; j++)
1612                    {
1613                        if( j==k )
1614                        {
1615                            result = result+Math.Abs(1-network.y[j]);
1616                        }
1617                        else
1618                        {
1619                            result = result+Math.Abs(network.y[j]);
1620                        }
1621                    }
1622                }
1623                else
1624                {
1625                   
1626                    //
1627                    // real outputs
1628                    //
1629                    for(j=0; j<=nout-1; j++)
1630                    {
1631                        result = result+Math.Abs(xy[i,nin+j]-network.y[j]);
1632                    }
1633                }
1634            }
1635            result = result/(npoints*nout);
1636            return result;
1637        }
1638
1639
1640        /*************************************************************************
1641        Average relative error on the test set
1642
1643        INPUT PARAMETERS:
1644            Network -   neural network
1645            XY      -   test set
1646            NPoints -   test set size
1647
1648        RESULT:
1649            Its meaning for regression task is obvious. As for
1650            classification task, it means average relative error when estimating
1651            posterior probability of belonging to the correct class.
1652
1653          -- ALGLIB --
1654             Copyright 11.03.2008 by Bochkanov Sergey
1655        *************************************************************************/
1656        public static double mlpavgrelerror(ref multilayerperceptron network,
1657            ref double[,] xy,
1658            int npoints)
1659        {
1660            double result = 0;
1661            int i = 0;
1662            int j = 0;
1663            int k = 0;
1664            int lk = 0;
1665            int nin = 0;
1666            int nout = 0;
1667            int wcount = 0;
1668            int i_ = 0;
1669
1670            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1671            result = 0;
1672            k = 0;
1673            for(i=0; i<=npoints-1; i++)
1674            {
1675                for(i_=0; i_<=nin-1;i_++)
1676                {
1677                    network.x[i_] = xy[i,i_];
1678                }
1679                mlpprocess(ref network, ref network.x, ref network.y);
1680                if( mlpissoftmax(ref network) )
1681                {
1682                   
1683                    //
1684                    // class labels
1685                    //
1686                    lk = (int)Math.Round(xy[i,nin]);
1687                    for(j=0; j<=nout-1; j++)
1688                    {
1689                        if( j==lk )
1690                        {
1691                            result = result+Math.Abs(1-network.y[j]);
1692                            k = k+1;
1693                        }
1694                    }
1695                }
1696                else
1697                {
1698                   
1699                    //
1700                    // real outputs
1701                    //
1702                    for(j=0; j<=nout-1; j++)
1703                    {
1704                        if( (double)(xy[i,nin+j])!=(double)(0) )
1705                        {
1706                            result = result+Math.Abs(xy[i,nin+j]-network.y[j])/Math.Abs(xy[i,nin+j]);
1707                            k = k+1;
1708                        }
1709                    }
1710                }
1711            }
1712            if( k!=0 )
1713            {
1714                result = result/k;
1715            }
1716            return result;
1717        }
1718
1719
1720        /*************************************************************************
1721        Gradient calculation. Internal subroutine.
1722
1723          -- ALGLIB --
1724             Copyright 04.11.2007 by Bochkanov Sergey
1725        *************************************************************************/
1726        public static void mlpgrad(ref multilayerperceptron network,
1727            ref double[] x,
1728            ref double[] desiredy,
1729            ref double e,
1730            ref double[] grad)
1731        {
1732            int i = 0;
1733            int nout = 0;
1734            int ntotal = 0;
1735
1736           
1737            //
1738            // Prepare dError/dOut, internal structures
1739            //
1740            mlpprocess(ref network, ref x, ref network.y);
1741            nout = network.structinfo[2];
1742            ntotal = network.structinfo[3];
1743            e = 0;
1744            for(i=0; i<=ntotal-1; i++)
1745            {
1746                network.derror[i] = 0;
1747            }
1748            for(i=0; i<=nout-1; i++)
1749            {
1750                network.derror[ntotal-nout+i] = network.y[i]-desiredy[i];
1751                e = e+AP.Math.Sqr(network.y[i]-desiredy[i])/2;
1752            }
1753           
1754            //
1755            // gradient
1756            //
1757            mlpinternalcalculategradient(ref network, ref network.neurons, ref network.weights, ref network.derror, ref grad, false);
1758        }
1759
1760
1761        /*************************************************************************
1762        Gradient calculation (natural error function). Internal subroutine.
1763
1764          -- ALGLIB --
1765             Copyright 04.11.2007 by Bochkanov Sergey
1766        *************************************************************************/
1767        public static void mlpgradn(ref multilayerperceptron network,
1768            ref double[] x,
1769            ref double[] desiredy,
1770            ref double e,
1771            ref double[] grad)
1772        {
1773            double s = 0;
1774            int i = 0;
1775            int nout = 0;
1776            int ntotal = 0;
1777
1778           
1779            //
1780            // Prepare dError/dOut, internal structures
1781            //
1782            mlpprocess(ref network, ref x, ref network.y);
1783            nout = network.structinfo[2];
1784            ntotal = network.structinfo[3];
1785            for(i=0; i<=ntotal-1; i++)
1786            {
1787                network.derror[i] = 0;
1788            }
1789            e = 0;
1790            if( network.structinfo[6]==0 )
1791            {
1792               
1793                //
1794                // Regression network, least squares
1795                //
1796                for(i=0; i<=nout-1; i++)
1797                {
1798                    network.derror[ntotal-nout+i] = network.y[i]-desiredy[i];
1799                    e = e+AP.Math.Sqr(network.y[i]-desiredy[i])/2;
1800                }
1801            }
1802            else
1803            {
1804               
1805                //
1806                // Classification network, cross-entropy
1807                //
1808                s = 0;
1809                for(i=0; i<=nout-1; i++)
1810                {
1811                    s = s+desiredy[i];
1812                }
1813                for(i=0; i<=nout-1; i++)
1814                {
1815                    network.derror[ntotal-nout+i] = s*network.y[i]-desiredy[i];
1816                    e = e+safecrossentropy(desiredy[i], network.y[i]);
1817                }
1818            }
1819           
1820            //
1821            // gradient
1822            //
1823            mlpinternalcalculategradient(ref network, ref network.neurons, ref network.weights, ref network.derror, ref grad, true);
1824        }
1825
1826
1827        /*************************************************************************
1828        Batch gradient calculation. Internal subroutine.
1829
1830          -- ALGLIB --
1831             Copyright 04.11.2007 by Bochkanov Sergey
1832        *************************************************************************/
1833        public static void mlpgradbatch(ref multilayerperceptron network,
1834            ref double[,] xy,
1835            int ssize,
1836            ref double e,
1837            ref double[] grad)
1838        {
1839            int i = 0;
1840            int nin = 0;
1841            int nout = 0;
1842            int wcount = 0;
1843
1844            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1845            for(i=0; i<=wcount-1; i++)
1846            {
1847                grad[i] = 0;
1848            }
1849            e = 0;
1850            i = 0;
1851            while( i<=ssize-1 )
1852            {
1853                mlpchunkedgradient(ref network, ref xy, i, Math.Min(ssize, i+chunksize)-i, ref e, ref grad, false);
1854                i = i+chunksize;
1855            }
1856        }
1857
1858
1859        /*************************************************************************
1860        Batch gradient calculation (natural error function). Internal subroutine.
1861
1862          -- ALGLIB --
1863             Copyright 04.11.2007 by Bochkanov Sergey
1864        *************************************************************************/
1865        public static void mlpgradnbatch(ref multilayerperceptron network,
1866            ref double[,] xy,
1867            int ssize,
1868            ref double e,
1869            ref double[] grad)
1870        {
1871            int i = 0;
1872            int nin = 0;
1873            int nout = 0;
1874            int wcount = 0;
1875
1876            mlpproperties(ref network, ref nin, ref nout, ref wcount);
1877            for(i=0; i<=wcount-1; i++)
1878            {
1879                grad[i] = 0;
1880            }
1881            e = 0;
1882            i = 0;
1883            while( i<=ssize-1 )
1884            {
1885                mlpchunkedgradient(ref network, ref xy, i, Math.Min(ssize, i+chunksize)-i, ref e, ref grad, true);
1886                i = i+chunksize;
1887            }
1888        }
1889
1890
1891        /*************************************************************************
1892        Batch Hessian calculation (natural error function) using R-algorithm.
1893        Internal subroutine.
1894
1895          -- ALGLIB --
1896             Copyright 26.01.2008 by Bochkanov Sergey.
1897             
1898             Hessian calculation based on R-algorithm described in
1899             "Fast Exact Multiplication by the Hessian",
1900             B. A. Pearlmutter,
1901             Neural Computation, 1994.
1902        *************************************************************************/
1903        public static void mlphessiannbatch(ref multilayerperceptron network,
1904            ref double[,] xy,
1905            int ssize,
1906            ref double e,
1907            ref double[] grad,
1908            ref double[,] h)
1909        {
1910            mlphessianbatchinternal(ref network, ref xy, ssize, true, ref e, ref grad, ref h);
1911        }
1912
1913
1914        /*************************************************************************
1915        Batch Hessian calculation using R-algorithm.
1916        Internal subroutine.
1917
1918          -- ALGLIB --
1919             Copyright 26.01.2008 by Bochkanov Sergey.
1920
1921             Hessian calculation based on R-algorithm described in
1922             "Fast Exact Multiplication by the Hessian",
1923             B. A. Pearlmutter,
1924             Neural Computation, 1994.
1925        *************************************************************************/
1926        public static void mlphessianbatch(ref multilayerperceptron network,
1927            ref double[,] xy,
1928            int ssize,
1929            ref double e,
1930            ref double[] grad,
1931            ref double[,] h)
1932        {
1933            mlphessianbatchinternal(ref network, ref xy, ssize, false, ref e, ref grad, ref h);
1934        }
1935
1936
1937        /*************************************************************************
1938        Internal subroutine, shouldn't be called by user.
1939        *************************************************************************/
1940        public static void mlpinternalprocessvector(ref int[] structinfo,
1941            ref double[] weights,
1942            ref double[] columnmeans,
1943            ref double[] columnsigmas,
1944            ref double[] neurons,
1945            ref double[] dfdnet,
1946            ref double[] x,
1947            ref double[] y)
1948        {
1949            int i = 0;
1950            int n1 = 0;
1951            int n2 = 0;
1952            int w1 = 0;
1953            int w2 = 0;
1954            int ntotal = 0;
1955            int nin = 0;
1956            int nout = 0;
1957            int istart = 0;
1958            int offs = 0;
1959            double net = 0;
1960            double f = 0;
1961            double df = 0;
1962            double d2f = 0;
1963            double mx = 0;
1964            bool perr = new bool();
1965            int i_ = 0;
1966            int i1_ = 0;
1967
1968           
1969            //
1970            // Read network geometry
1971            //
1972            nin = structinfo[1];
1973            nout = structinfo[2];
1974            ntotal = structinfo[3];
1975            istart = structinfo[5];
1976           
1977            //
1978            // Inputs standartisation and putting in the network
1979            //
1980            for(i=0; i<=nin-1; i++)
1981            {
1982                if( (double)(columnsigmas[i])!=(double)(0) )
1983                {
1984                    neurons[i] = (x[i]-columnmeans[i])/columnsigmas[i];
1985                }
1986                else
1987                {
1988                    neurons[i] = x[i]-columnmeans[i];
1989                }
1990            }
1991           
1992            //
1993            // Process network
1994            //
1995            for(i=0; i<=ntotal-1; i++)
1996            {
1997                offs = istart+i*nfieldwidth;
1998                if( structinfo[offs+0]>0 )
1999                {
2000                   
2001                    //
2002                    // Activation function
2003                    //
2004                    mlpactivationfunction(neurons[structinfo[offs+2]], structinfo[offs+0], ref f, ref df, ref d2f);
2005                    neurons[i] = f;
2006                    dfdnet[i] = df;
2007                }
2008                if( structinfo[offs+0]==0 )
2009                {
2010                   
2011                    //
2012                    // Adaptive summator
2013                    //
2014                    n1 = structinfo[offs+2];
2015                    n2 = n1+structinfo[offs+1]-1;
2016                    w1 = structinfo[offs+3];
2017                    w2 = w1+structinfo[offs+1]-1;
2018                    i1_ = (n1)-(w1);
2019                    net = 0.0;
2020                    for(i_=w1; i_<=w2;i_++)
2021                    {
2022                        net += weights[i_]*neurons[i_+i1_];
2023                    }
2024                    neurons[i] = net;
2025                    dfdnet[i] = 1.0;
2026                }
2027                if( structinfo[offs+0]<0 )
2028                {
2029                    perr = true;
2030                    if( structinfo[offs+0]==-2 )
2031                    {
2032                       
2033                        //
2034                        // input neuron, left unchanged
2035                        //
2036                        perr = false;
2037                    }
2038                    if( structinfo[offs+0]==-3 )
2039                    {
2040                       
2041                        //
2042                        // "-1" neuron
2043                        //
2044                        neurons[i] = -1;
2045                        perr = false;
2046                    }
2047                    if( structinfo[offs+0]==-4 )
2048                    {
2049                       
2050                        //
2051                        // "0" neuron
2052                        //
2053                        neurons[i] = 0;
2054                        perr = false;
2055                    }
2056                    System.Diagnostics.Debug.Assert(!perr, "MLPInternalProcessVector: internal error - unknown neuron type!");
2057                }
2058            }
2059           
2060            //
2061            // Extract result
2062            //
2063            i1_ = (ntotal-nout) - (0);
2064            for(i_=0; i_<=nout-1;i_++)
2065            {
2066                y[i_] = neurons[i_+i1_];
2067            }
2068           
2069            //
2070            // Softmax post-processing or standardisation if needed
2071            //
2072            System.Diagnostics.Debug.Assert(structinfo[6]==0 | structinfo[6]==1, "MLPInternalProcessVector: unknown normalization type!");
2073            if( structinfo[6]==1 )
2074            {
2075               
2076                //
2077                // Softmax
2078                //
2079                mx = y[0];
2080                for(i=1; i<=nout-1; i++)
2081                {
2082                    mx = Math.Max(mx, y[i]);
2083                }
2084                net = 0;
2085                for(i=0; i<=nout-1; i++)
2086                {
2087                    y[i] = Math.Exp(y[i]-mx);
2088                    net = net+y[i];
2089                }
2090                for(i=0; i<=nout-1; i++)
2091                {
2092                    y[i] = y[i]/net;
2093                }
2094            }
2095            else
2096            {
2097               
2098                //
2099                // Standardisation
2100                //
2101                for(i=0; i<=nout-1; i++)
2102                {
2103                    y[i] = y[i]*columnsigmas[nin+i]+columnmeans[nin+i];
2104                }
2105            }
2106        }
2107
2108
2109        /*************************************************************************
2110        Internal subroutine: adding new input layer to network
2111        *************************************************************************/
2112        private static void addinputlayer(int ncount,
2113            ref int[] lsizes,
2114            ref int[] ltypes,
2115            ref int[] lconnfirst,
2116            ref int[] lconnlast,
2117            ref int lastproc)
2118        {
2119            lsizes[0] = ncount;
2120            ltypes[0] = -2;
2121            lconnfirst[0] = 0;
2122            lconnlast[0] = 0;
2123            lastproc = 0;
2124        }
2125
2126
2127        /*************************************************************************
2128        Internal subroutine: adding new summator layer to network
2129        *************************************************************************/
2130        private static void addbiasedsummatorlayer(int ncount,
2131            ref int[] lsizes,
2132            ref int[] ltypes,
2133            ref int[] lconnfirst,
2134            ref int[] lconnlast,
2135            ref int lastproc)
2136        {
2137            lsizes[lastproc+1] = 1;
2138            ltypes[lastproc+1] = -3;
2139            lconnfirst[lastproc+1] = 0;
2140            lconnlast[lastproc+1] = 0;
2141            lsizes[lastproc+2] = ncount;
2142            ltypes[lastproc+2] = 0;
2143            lconnfirst[lastproc+2] = lastproc;
2144            lconnlast[lastproc+2] = lastproc+1;
2145            lastproc = lastproc+2;
2146        }
2147
2148
2149        /*************************************************************************
2150        Internal subroutine: adding new summator layer to network
2151        *************************************************************************/
2152        private static void addactivationlayer(int functype,
2153            ref int[] lsizes,
2154            ref int[] ltypes,
2155            ref int[] lconnfirst,
2156            ref int[] lconnlast,
2157            ref int lastproc)
2158        {
2159            System.Diagnostics.Debug.Assert(functype>0, "AddActivationLayer: incorrect function type");
2160            lsizes[lastproc+1] = lsizes[lastproc];
2161            ltypes[lastproc+1] = functype;
2162            lconnfirst[lastproc+1] = lastproc;
2163            lconnlast[lastproc+1] = lastproc;
2164            lastproc = lastproc+1;
2165        }
2166
2167
2168        /*************************************************************************
2169        Internal subroutine: adding new zero layer to network
2170        *************************************************************************/
2171        private static void addzerolayer(ref int[] lsizes,
2172            ref int[] ltypes,
2173            ref int[] lconnfirst,
2174            ref int[] lconnlast,
2175            ref int lastproc)
2176        {
2177            lsizes[lastproc+1] = 1;
2178            ltypes[lastproc+1] = -4;
2179            lconnfirst[lastproc+1] = 0;
2180            lconnlast[lastproc+1] = 0;
2181            lastproc = lastproc+1;
2182        }
2183
2184
2185        /*************************************************************************
2186        Internal subroutine.
2187
2188          -- ALGLIB --
2189             Copyright 04.11.2007 by Bochkanov Sergey
2190        *************************************************************************/
2191        private static void mlpcreate(int nin,
2192            int nout,
2193            ref int[] lsizes,
2194            ref int[] ltypes,
2195            ref int[] lconnfirst,
2196            ref int[] lconnlast,
2197            int layerscount,
2198            bool isclsnet,
2199            ref multilayerperceptron network)
2200        {
2201            int i = 0;
2202            int j = 0;
2203            int ssize = 0;
2204            int ntotal = 0;
2205            int wcount = 0;
2206            int offs = 0;
2207            int nprocessed = 0;
2208            int wallocated = 0;
2209            int[] localtemp = new int[0];
2210            int[] lnfirst = new int[0];
2211            int[] lnsyn = new int[0];
2212
2213           
2214            //
2215            // Check
2216            //
2217            System.Diagnostics.Debug.Assert(layerscount>0, "MLPCreate: wrong parameters!");
2218            System.Diagnostics.Debug.Assert(ltypes[0]==-2, "MLPCreate: wrong LTypes[0] (must be -2)!");
2219            for(i=0; i<=layerscount-1; i++)
2220            {
2221                System.Diagnostics.Debug.Assert(lsizes[i]>0, "MLPCreate: wrong LSizes!");
2222                System.Diagnostics.Debug.Assert(lconnfirst[i]>=0 & (lconnfirst[i]<i | i==0), "MLPCreate: wrong LConnFirst!");
2223                System.Diagnostics.Debug.Assert(lconnlast[i]>=lconnfirst[i] & (lconnlast[i]<i | i==0), "MLPCreate: wrong LConnLast!");
2224            }
2225           
2226            //
2227            // Build network geometry
2228            //
2229            lnfirst = new int[layerscount-1+1];
2230            lnsyn = new int[layerscount-1+1];
2231            ntotal = 0;
2232            wcount = 0;
2233            for(i=0; i<=layerscount-1; i++)
2234            {
2235               
2236                //
2237                // Analyze connections.
2238                // This code must throw an assertion in case of unknown LTypes[I]
2239                //
2240                lnsyn[i] = -1;
2241                if( ltypes[i]>=0 )
2242                {
2243                    lnsyn[i] = 0;
2244                    for(j=lconnfirst[i]; j<=lconnlast[i]; j++)
2245                    {
2246                        lnsyn[i] = lnsyn[i]+lsizes[j];
2247                    }
2248                }
2249                else
2250                {
2251                    if( ltypes[i]==-2 | ltypes[i]==-3 | ltypes[i]==-4 )
2252                    {
2253                        lnsyn[i] = 0;
2254                    }
2255                }
2256                System.Diagnostics.Debug.Assert(lnsyn[i]>=0, "MLPCreate: internal error #0!");
2257               
2258                //
2259                // Other info
2260                //
2261                lnfirst[i] = ntotal;
2262                ntotal = ntotal+lsizes[i];
2263                if( ltypes[i]==0 )
2264                {
2265                    wcount = wcount+lnsyn[i]*lsizes[i];
2266                }
2267            }
2268            ssize = 7+ntotal*nfieldwidth;
2269           
2270            //
2271            // Allocate
2272            //
2273            network.structinfo = new int[ssize-1+1];
2274            network.weights = new double[wcount-1+1];
2275            if( isclsnet )
2276            {
2277                network.columnmeans = new double[nin-1+1];
2278                network.columnsigmas = new double[nin-1+1];
2279            }
2280            else
2281            {
2282                network.columnmeans = new double[nin+nout-1+1];
2283                network.columnsigmas = new double[nin+nout-1+1];
2284            }
2285            network.neurons = new double[ntotal-1+1];
2286            network.chunks = new double[3*ntotal+1, chunksize-1+1];
2287            network.nwbuf = new double[Math.Max(wcount, 2*nout)-1+1];
2288            network.dfdnet = new double[ntotal-1+1];
2289            network.x = new double[nin-1+1];
2290            network.y = new double[nout-1+1];
2291            network.derror = new double[ntotal-1+1];
2292           
2293            //
2294            // Fill structure: global info
2295            //
2296            network.structinfo[0] = ssize;
2297            network.structinfo[1] = nin;
2298            network.structinfo[2] = nout;
2299            network.structinfo[3] = ntotal;
2300            network.structinfo[4] = wcount;
2301            network.structinfo[5] = 7;
2302            if( isclsnet )
2303            {
2304                network.structinfo[6] = 1;
2305            }
2306            else
2307            {
2308                network.structinfo[6] = 0;
2309            }
2310           
2311            //
2312            // Fill structure: neuron connections
2313            //
2314            nprocessed = 0;
2315            wallocated = 0;
2316            for(i=0; i<=layerscount-1; i++)
2317            {
2318                for(j=0; j<=lsizes[i]-1; j++)
2319                {
2320                    offs = network.structinfo[5]+nprocessed*nfieldwidth;
2321                    network.structinfo[offs+0] = ltypes[i];
2322                    if( ltypes[i]==0 )
2323                    {
2324                       
2325                        //
2326                        // Adaptive summator:
2327                        // * connections with weights to previous neurons
2328                        //
2329                        network.structinfo[offs+1] = lnsyn[i];
2330                        network.structinfo[offs+2] = lnfirst[lconnfirst[i]];
2331                        network.structinfo[offs+3] = wallocated;
2332                        wallocated = wallocated+lnsyn[i];
2333                        nprocessed = nprocessed+1;
2334                    }
2335                    if( ltypes[i]>0 )
2336                    {
2337                       
2338                        //
2339                        // Activation layer:
2340                        // * each neuron connected to one (only one) of previous neurons.
2341                        // * no weights
2342                        //
2343                        network.structinfo[offs+1] = 1;
2344                        network.structinfo[offs+2] = lnfirst[lconnfirst[i]]+j;
2345                        network.structinfo[offs+3] = -1;
2346                        nprocessed = nprocessed+1;
2347                    }
2348                    if( ltypes[i]==-2 | ltypes[i]==-3 | ltypes[i]==-4 )
2349                    {
2350                        nprocessed = nprocessed+1;
2351                    }
2352                }
2353            }
2354            System.Diagnostics.Debug.Assert(wallocated==wcount, "MLPCreate: internal error #1!");
2355            System.Diagnostics.Debug.Assert(nprocessed==ntotal, "MLPCreate: internal error #2!");
2356           
2357            //
2358            // Fill weights by small random values
2359            // Initialize means and sigmas
2360            //
2361            for(i=0; i<=wcount-1; i++)
2362            {
2363                network.weights[i] = AP.Math.RandomReal()-0.5;
2364            }
2365            for(i=0; i<=nin-1; i++)
2366            {
2367                network.columnmeans[i] = 0;
2368                network.columnsigmas[i] = 1;
2369            }
2370            if( !isclsnet )
2371            {
2372                for(i=0; i<=nout-1; i++)
2373                {
2374                    network.columnmeans[nin+i] = 0;
2375                    network.columnsigmas[nin+i] = 1;
2376                }
2377            }
2378        }
2379
2380
2381        /*************************************************************************
2382        Internal subroutine
2383
2384          -- ALGLIB --
2385             Copyright 04.11.2007 by Bochkanov Sergey
2386        *************************************************************************/
2387        private static void mlpactivationfunction(double net,
2388            int k,
2389            ref double f,
2390            ref double df,
2391            ref double d2f)
2392        {
2393            double net2 = 0;
2394            double arg = 0;
2395            double root = 0;
2396            double r = 0;
2397
2398            f = 0;
2399            df = 0;
2400            if( k==1 )
2401            {
2402               
2403                //
2404                // TanH activation function
2405                //
2406                if( (double)(Math.Abs(net))<(double)(100) )
2407                {
2408                    f = Math.Tanh(net);
2409                }
2410                else
2411                {
2412                    f = Math.Sign(net);
2413                }
2414                df = 1-AP.Math.Sqr(f);
2415                d2f = -(2*f*df);
2416                return;
2417            }
2418            if( k==3 )
2419            {
2420               
2421                //
2422                // EX activation function
2423                //
2424                if( (double)(net)>=(double)(0) )
2425                {
2426                    net2 = net*net;
2427                    arg = net2+1;
2428                    root = Math.Sqrt(arg);
2429                    f = net+root;
2430                    r = net/root;
2431                    df = 1+r;
2432                    d2f = (root-net*r)/arg;
2433                }
2434                else
2435                {
2436                    f = Math.Exp(net);
2437                    df = f;
2438                    d2f = f;
2439                }
2440                return;
2441            }
2442            if( k==2 )
2443            {
2444                f = Math.Exp(-AP.Math.Sqr(net));
2445                df = -(2*net*f);
2446                d2f = -(2*(f+df*net));
2447                return;
2448            }
2449        }
2450
2451
2452        /*************************************************************************
2453        Internal subroutine for Hessian calculation.
2454
2455        WARNING!!! Unspeakable math far beyong human capabilities :)
2456        *************************************************************************/
2457        private static void mlphessianbatchinternal(ref multilayerperceptron network,
2458            ref double[,] xy,
2459            int ssize,
2460            bool naturalerr,
2461            ref double e,
2462            ref double[] grad,
2463            ref double[,] h)
2464        {
2465            int nin = 0;
2466            int nout = 0;
2467            int wcount = 0;
2468            int ntotal = 0;
2469            int istart = 0;
2470            int i = 0;
2471            int j = 0;
2472            int k = 0;
2473            int kl = 0;
2474            int offs = 0;
2475            int n1 = 0;
2476            int n2 = 0;
2477            int w1 = 0;
2478            int w2 = 0;
2479            double s = 0;
2480            double t = 0;
2481            double v = 0;
2482            double et = 0;
2483            bool bflag = new bool();
2484            double f = 0;
2485            double df = 0;
2486            double d2f = 0;
2487            double deidyj = 0;
2488            double mx = 0;
2489            double q = 0;
2490            double z = 0;
2491            double s2 = 0;
2492            double expi = 0;
2493            double expj = 0;
2494            double[] x = new double[0];
2495            double[] desiredy = new double[0];
2496            double[] gt = new double[0];
2497            double[] zeros = new double[0];
2498            double[,] rx = new double[0,0];
2499            double[,] ry = new double[0,0];
2500            double[,] rdx = new double[0,0];
2501            double[,] rdy = new double[0,0];
2502            int i_ = 0;
2503            int i1_ = 0;
2504
2505            mlpproperties(ref network, ref nin, ref nout, ref wcount);
2506            ntotal = network.structinfo[3];
2507            istart = network.structinfo[5];
2508           
2509            //
2510            // Prepare
2511            //
2512            x = new double[nin-1+1];
2513            desiredy = new double[nout-1+1];
2514            zeros = new double[wcount-1+1];
2515            gt = new double[wcount-1+1];
2516            rx = new double[ntotal+nout-1+1, wcount-1+1];
2517            ry = new double[ntotal+nout-1+1, wcount-1+1];
2518            rdx = new double[ntotal+nout-1+1, wcount-1+1];
2519            rdy = new double[ntotal+nout-1+1, wcount-1+1];
2520            e = 0;
2521            for(i=0; i<=wcount-1; i++)
2522            {
2523                zeros[i] = 0;
2524            }
2525            for(i_=0; i_<=wcount-1;i_++)
2526            {
2527                grad[i_] = zeros[i_];
2528            }
2529            for(i=0; i<=wcount-1; i++)
2530            {
2531                for(i_=0; i_<=wcount-1;i_++)
2532                {
2533                    h[i,i_] = zeros[i_];
2534                }
2535            }
2536           
2537            //
2538            // Process
2539            //
2540            for(k=0; k<=ssize-1; k++)
2541            {
2542               
2543                //
2544                // Process vector with MLPGradN.
2545                // Now Neurons, DFDNET and DError contains results of the last run.
2546                //
2547                for(i_=0; i_<=nin-1;i_++)
2548                {
2549                    x[i_] = xy[k,i_];
2550                }
2551                if( mlpissoftmax(ref network) )
2552                {
2553                   
2554                    //
2555                    // class labels outputs
2556                    //
2557                    kl = (int)Math.Round(xy[k,nin]);
2558                    for(i=0; i<=nout-1; i++)
2559                    {
2560                        if( i==kl )
2561                        {
2562                            desiredy[i] = 1;
2563                        }
2564                        else
2565                        {
2566                            desiredy[i] = 0;
2567                        }
2568                    }
2569                }
2570                else
2571                {
2572                   
2573                    //
2574                    // real outputs
2575                    //
2576                    i1_ = (nin) - (0);
2577                    for(i_=0; i_<=nout-1;i_++)
2578                    {
2579                        desiredy[i_] = xy[k,i_+i1_];
2580                    }
2581                }
2582                if( naturalerr )
2583                {
2584                    mlpgradn(ref network, ref x, ref desiredy, ref et, ref gt);
2585                }
2586                else
2587                {
2588                    mlpgrad(ref network, ref x, ref desiredy, ref et, ref gt);
2589                }
2590               
2591                //
2592                // grad, error
2593                //
2594                e = e+et;
2595                for(i_=0; i_<=wcount-1;i_++)
2596                {
2597                    grad[i_] = grad[i_] + gt[i_];
2598                }
2599               
2600                //
2601                // Hessian.
2602                // Forward pass of the R-algorithm
2603                //
2604                for(i=0; i<=ntotal-1; i++)
2605                {
2606                    offs = istart+i*nfieldwidth;
2607                    for(i_=0; i_<=wcount-1;i_++)
2608                    {
2609                        rx[i,i_] = zeros[i_];
2610                    }
2611                    for(i_=0; i_<=wcount-1;i_++)
2612                    {
2613                        ry[i,i_] = zeros[i_];
2614                    }
2615                    if( network.structinfo[offs+0]>0 )
2616                    {
2617                       
2618                        //
2619                        // Activation function
2620                        //
2621                        n1 = network.structinfo[offs+2];
2622                        for(i_=0; i_<=wcount-1;i_++)
2623                        {
2624                            rx[i,i_] = ry[n1,i_];
2625                        }
2626                        v = network.dfdnet[i];
2627                        for(i_=0; i_<=wcount-1;i_++)
2628                        {
2629                            ry[i,i_] = v*rx[i,i_];
2630                        }
2631                    }
2632                    if( network.structinfo[offs+0]==0 )
2633                    {
2634                       
2635                        //
2636                        // Adaptive summator
2637                        //
2638                        n1 = network.structinfo[offs+2];
2639                        n2 = n1+network.structinfo[offs+1]-1;
2640                        w1 = network.structinfo[offs+3];
2641                        w2 = w1+network.structinfo[offs+1]-1;
2642                        for(j=n1; j<=n2; j++)
2643                        {
2644                            v = network.weights[w1+j-n1];
2645                            for(i_=0; i_<=wcount-1;i_++)
2646                            {
2647                                rx[i,i_] = rx[i,i_] + v*ry[j,i_];
2648                            }
2649                            rx[i,w1+j-n1] = rx[i,w1+j-n1]+network.neurons[j];
2650                        }
2651                        for(i_=0; i_<=wcount-1;i_++)
2652                        {
2653                            ry[i,i_] = rx[i,i_];
2654                        }
2655                    }
2656                    if( network.structinfo[offs+0]<0 )
2657                    {
2658                        bflag = true;
2659                        if( network.structinfo[offs+0]==-2 )
2660                        {
2661                           
2662                            //
2663                            // input neuron, left unchanged
2664                            //
2665                            bflag = false;
2666                        }
2667                        if( network.structinfo[offs+0]==-3 )
2668                        {
2669                           
2670                            //
2671                            // "-1" neuron, left unchanged
2672                            //
2673                            bflag = false;
2674                        }
2675                        if( network.structinfo[offs+0]==-4 )
2676                        {
2677                           
2678                            //
2679                            // "0" neuron, left unchanged
2680                            //
2681                            bflag = false;
2682                        }
2683                        System.Diagnostics.Debug.Assert(!bflag, "MLPHessianNBatch: internal error - unknown neuron type!");
2684                    }
2685                }
2686               
2687                //
2688                // Hessian. Backward pass of the R-algorithm.
2689                //
2690                // Stage 1. Initialize RDY
2691                //
2692                for(i=0; i<=ntotal+nout-1; i++)
2693                {
2694                    for(i_=0; i_<=wcount-1;i_++)
2695                    {
2696                        rdy[i,i_] = zeros[i_];
2697                    }
2698                }
2699                if( network.structinfo[6]==0 )
2700                {
2701                   
2702                    //
2703                    // Standardisation.
2704                    //
2705                    // In context of the Hessian calculation standardisation
2706                    // is considered as additional layer with weightless
2707                    // activation function:
2708                    //
2709                    // F(NET) := Sigma*NET
2710                    //
2711                    // So we add one more layer to forward pass, and
2712                    // make forward/backward pass through this layer.
2713                    //
2714                    for(i=0; i<=nout-1; i++)
2715                    {
2716                        n1 = ntotal-nout+i;
2717                        n2 = ntotal+i;
2718                       
2719                        //
2720                        // Forward pass from N1 to N2
2721                        //
2722                        for(i_=0; i_<=wcount-1;i_++)
2723                        {
2724                            rx[n2,i_] = ry[n1,i_];
2725                        }
2726                        v = network.columnsigmas[nin+i];
2727                        for(i_=0; i_<=wcount-1;i_++)
2728                        {
2729                            ry[n2,i_] = v*rx[n2,i_];
2730                        }
2731                       
2732                        //
2733                        // Initialization of RDY
2734                        //
2735                        for(i_=0; i_<=wcount-1;i_++)
2736                        {
2737                            rdy[n2,i_] = ry[n2,i_];
2738                        }
2739                       
2740                        //
2741                        // Backward pass from N2 to N1:
2742                        // 1. Calculate R(dE/dX).
2743                        // 2. No R(dE/dWij) is needed since weight of activation neuron
2744                        //    is fixed to 1. So we can update R(dE/dY) for
2745                        //    the connected neuron (note that Vij=0, Wij=1)
2746                        //
2747                        df = network.columnsigmas[nin+i];
2748                        for(i_=0; i_<=wcount-1;i_++)
2749                        {
2750                            rdx[n2,i_] = df*rdy[n2,i_];
2751                        }
2752                        for(i_=0; i_<=wcount-1;i_++)
2753                        {
2754                            rdy[n1,i_] = rdy[n1,i_] + rdx[n2,i_];
2755                        }
2756                    }
2757                }
2758                else
2759                {
2760                   
2761                    //
2762                    // Softmax.
2763                    //
2764                    // Initialize RDY using generalized expression for ei'(yi)
2765                    // (see expression (9) from p. 5 of "Fast Exact Multiplication by the Hessian").
2766                    //
2767                    // When we are working with softmax network, generalized
2768                    // expression for ei'(yi) is used because softmax
2769                    // normalization leads to ei, which depends on all y's
2770                    //
2771                    if( naturalerr )
2772                    {
2773                       
2774                        //
2775                        // softmax + cross-entropy.
2776                        // We have:
2777                        //
2778                        // S = sum(exp(yk)),
2779                        // ei = sum(trn)*exp(yi)/S-trn_i
2780                        //
2781                        // j=i:   d(ei)/d(yj) = T*exp(yi)*(S-exp(yi))/S^2
2782                        // j<>i:  d(ei)/d(yj) = -T*exp(yi)*exp(yj)/S^2
2783                        //
2784                        t = 0;
2785                        for(i=0; i<=nout-1; i++)
2786                        {
2787                            t = t+desiredy[i];
2788                        }
2789                        mx = network.neurons[ntotal-nout];
2790                        for(i=0; i<=nout-1; i++)
2791                        {
2792                            mx = Math.Max(mx, network.neurons[ntotal-nout+i]);
2793                        }
2794                        s = 0;
2795                        for(i=0; i<=nout-1; i++)
2796                        {
2797                            network.nwbuf[i] = Math.Exp(network.neurons[ntotal-nout+i]-mx);
2798                            s = s+network.nwbuf[i];
2799                        }
2800                        for(i=0; i<=nout-1; i++)
2801                        {
2802                            for(j=0; j<=nout-1; j++)
2803                            {
2804                                if( j==i )
2805                                {
2806                                    deidyj = t*network.nwbuf[i]*(s-network.nwbuf[i])/AP.Math.Sqr(s);
2807                                    for(i_=0; i_<=wcount-1;i_++)
2808                                    {
2809                                        rdy[ntotal-nout+i,i_] = rdy[ntotal-nout+i,i_] + deidyj*ry[ntotal-nout+i,i_];
2810                                    }
2811                                }
2812                                else
2813                                {
2814                                    deidyj = -(t*network.nwbuf[i]*network.nwbuf[j]/AP.Math.Sqr(s));
2815                                    for(i_=0; i_<=wcount-1;i_++)
2816                                    {
2817                                        rdy[ntotal-nout+i,i_] = rdy[ntotal-nout+i,i_] + deidyj*ry[ntotal-nout+j,i_];
2818                                    }
2819                                }
2820                            }
2821                        }
2822                    }
2823                    else
2824                    {
2825                       
2826                        //
2827                        // For a softmax + squared error we have expression
2828                        // far beyond human imagination so we dont even try
2829                        // to comment on it. Just enjoy the code...
2830                        //
2831                        // P.S. That's why "natural error" is called "natural" -
2832                        // compact beatiful expressions, fast code....
2833                        //
2834                        mx = network.neurons[ntotal-nout];
2835                        for(i=0; i<=nout-1; i++)
2836                        {
2837                            mx = Math.Max(mx, network.neurons[ntotal-nout+i]);
2838                        }
2839                        s = 0;
2840                        s2 = 0;
2841                        for(i=0; i<=nout-1; i++)
2842                        {
2843                            network.nwbuf[i] = Math.Exp(network.neurons[ntotal-nout+i]-mx);
2844                            s = s+network.nwbuf[i];
2845                            s2 = s2+AP.Math.Sqr(network.nwbuf[i]);
2846                        }
2847                        q = 0;
2848                        for(i=0; i<=nout-1; i++)
2849                        {
2850                            q = q+(network.y[i]-desiredy[i])*network.nwbuf[i];
2851                        }
2852                        for(i=0; i<=nout-1; i++)
2853                        {
2854                            z = -q+(network.y[i]-desiredy[i])*s;
2855                            expi = network.nwbuf[i];
2856                            for(j=0; j<=nout-1; j++)
2857                            {
2858                                expj = network.nwbuf[j];
2859                                if( j==i )
2860                                {
2861                                    deidyj = expi/AP.Math.Sqr(s)*((z+expi)*(s-2*expi)/s+expi*s2/AP.Math.Sqr(s));
2862                                }
2863                                else
2864                                {
2865                                    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]));
2866                                }
2867                                for(i_=0; i_<=wcount-1;i_++)
2868                                {
2869                                    rdy[ntotal-nout+i,i_] = rdy[ntotal-nout+i,i_] + deidyj*ry[ntotal-nout+j,i_];
2870                                }
2871                            }
2872                        }
2873                    }
2874                }
2875               
2876                //
2877                // Hessian. Backward pass of the R-algorithm
2878                //
2879                // Stage 2. Process.
2880                //
2881                for(i=ntotal-1; i>=0; i--)
2882                {
2883                   
2884                    //
2885                    // Possible variants:
2886                    // 1. Activation function
2887                    // 2. Adaptive summator
2888                    // 3. Special neuron
2889                    //
2890                    offs = istart+i*nfieldwidth;
2891                    if( network.structinfo[offs+0]>0 )
2892                    {
2893                        n1 = network.structinfo[offs+2];
2894                       
2895                        //
2896                        // First, calculate R(dE/dX).
2897                        //
2898                        mlpactivationfunction(network.neurons[n1], network.structinfo[offs+0], ref f, ref df, ref d2f);
2899                        v = d2f*network.derror[i];
2900                        for(i_=0; i_<=wcount-1;i_++)
2901                        {
2902                            rdx[i,i_] = df*rdy[i,i_];
2903                        }
2904                        for(i_=0; i_<=wcount-1;i_++)
2905                        {
2906                            rdx[i,i_] = rdx[i,i_] + v*rx[i,i_];
2907                        }
2908                       
2909                        //
2910                        // No R(dE/dWij) is needed since weight of activation neuron
2911                        // is fixed to 1.
2912                        //
2913                        // So we can update R(dE/dY) for the connected neuron.
2914                        // (note that Vij=0, Wij=1)
2915                        //
2916                        for(i_=0; i_<=wcount-1;i_++)
2917                        {
2918                            rdy[n1,i_] = rdy[n1,i_] + rdx[i,i_];
2919                        }
2920                    }
2921                    if( network.structinfo[offs+0]==0 )
2922                    {
2923                       
2924                        //
2925                        // Adaptive summator
2926                        //
2927                        n1 = network.structinfo[offs+2];
2928                        n2 = n1+network.structinfo[offs+1]-1;
2929                        w1 = network.structinfo[offs+3];
2930                        w2 = w1+network.structinfo[offs+1]-1;
2931                       
2932                        //
2933                        // First, calculate R(dE/dX).
2934                        //
2935                        for(i_=0; i_<=wcount-1;i_++)
2936                        {
2937                            rdx[i,i_] = rdy[i,i_];
2938                        }
2939                       
2940                        //
2941                        // Then, calculate R(dE/dWij)
2942                        //
2943                        for(j=w1; j<=w2; j++)
2944                        {
2945                            v = network.neurons[n1+j-w1];
2946                            for(i_=0; i_<=wcount-1;i_++)
2947                            {
2948                                h[j,i_] = h[j,i_] + v*rdx[i,i_];
2949                            }
2950                            v = network.derror[i];
2951                            for(i_=0; i_<=wcount-1;i_++)
2952                            {
2953                                h[j,i_] = h[j,i_] + v*ry[n1+j-w1,i_];
2954                            }
2955                        }
2956                       
2957                        //
2958                        // And finally, update R(dE/dY) for connected neurons.
2959                        //
2960                        for(j=w1; j<=w2; j++)
2961                        {
2962                            v = network.weights[j];
2963                            for(i_=0; i_<=wcount-1;i_++)
2964                            {
2965                                rdy[n1+j-w1,i_] = rdy[n1+j-w1,i_] + v*rdx[i,i_];
2966                            }
2967                            rdy[n1+j-w1,j] = rdy[n1+j-w1,j]+network.derror[i];
2968                        }
2969                    }
2970                    if( network.structinfo[offs+0]<0 )
2971                    {
2972                        bflag = false;
2973                        if( network.structinfo[offs+0]==-2 | network.structinfo[offs+0]==-3 | network.structinfo[offs+0]==-4 )
2974                        {
2975                           
2976                            //
2977                            // Special neuron type, no back-propagation required
2978                            //
2979                            bflag = true;
2980                        }
2981                        System.Diagnostics.Debug.Assert(bflag, "MLPHessianNBatch: unknown neuron type!");
2982                    }
2983                }
2984            }
2985        }
2986
2987
2988        /*************************************************************************
2989        Internal subroutine
2990
2991        Network must be processed by MLPProcess on X
2992        *************************************************************************/
2993        private static void mlpinternalcalculategradient(ref multilayerperceptron network,
2994            ref double[] neurons,
2995            ref double[] weights,
2996            ref double[] derror,
2997            ref double[] grad,
2998            bool naturalerrorfunc)
2999        {
3000            int i = 0;
3001            int n1 = 0;
3002            int n2 = 0;
3003            int w1 = 0;
3004            int w2 = 0;
3005            int ntotal = 0;
3006            int istart = 0;
3007            int nin = 0;
3008            int nout = 0;
3009            int offs = 0;
3010            double dedf = 0;
3011            double dfdnet = 0;
3012            double v = 0;
3013            double fown = 0;
3014            double deown = 0;
3015            double net = 0;
3016            double mx = 0;
3017            bool bflag = new bool();
3018            int i_ = 0;
3019            int i1_ = 0;
3020
3021           
3022            //
3023            // Read network geometry
3024            //
3025            nin = network.structinfo[1];
3026            nout = network.structinfo[2];
3027            ntotal = network.structinfo[3];
3028            istart = network.structinfo[5];
3029           
3030            //
3031            // Pre-processing of dError/dOut:
3032            // from dError/dOut(normalized) to dError/dOut(non-normalized)
3033            //
3034            System.Diagnostics.Debug.Assert(network.structinfo[6]==0 | network.structinfo[6]==1, "MLPInternalCalculateGradient: unknown normalization type!");
3035            if( network.structinfo[6]==1 )
3036            {
3037               
3038                //
3039                // Softmax
3040                //
3041                if( !naturalerrorfunc )
3042                {
3043                    mx = network.neurons[ntotal-nout];
3044                    for(i=0; i<=nout-1; i++)
3045                    {
3046                        mx = Math.Max(mx, network.neurons[ntotal-nout+i]);
3047                    }
3048                    net = 0;
3049                    for(i=0; i<=nout-1; i++)
3050                    {
3051                        network.nwbuf[i] = Math.Exp(network.neurons[ntotal-nout+i]-mx);
3052                        net = net+network.nwbuf[i];
3053                    }
3054                    i1_ = (0)-(ntotal-nout);
3055                    v = 0.0;
3056                    for(i_=ntotal-nout; i_<=ntotal-1;i_++)
3057                    {
3058                        v += network.derror[i_]*network.nwbuf[i_+i1_];
3059                    }
3060                    for(i=0; i<=nout-1; i++)
3061                    {
3062                        fown = network.nwbuf[i];
3063                        deown = network.derror[ntotal-nout+i];
3064                        network.nwbuf[nout+i] = (-v+deown*fown+deown*(net-fown))*fown/AP.Math.Sqr(net);
3065                    }
3066                    for(i=0; i<=nout-1; i++)
3067                    {
3068                        network.derror[ntotal-nout+i] = network.nwbuf[nout+i];
3069                    }
3070                }
3071            }
3072            else
3073            {
3074               
3075                //
3076                // Un-standardisation
3077                //
3078                for(i=0; i<=nout-1; i++)
3079                {
3080                    network.derror[ntotal-nout+i] = network.derror[ntotal-nout+i]*network.columnsigmas[nin+i];
3081                }
3082            }
3083           
3084            //
3085            // Backpropagation
3086            //
3087            for(i=ntotal-1; i>=0; i--)
3088            {
3089               
3090                //
3091                // Extract info
3092                //
3093                offs = istart+i*nfieldwidth;
3094                if( network.structinfo[offs+0]>0 )
3095                {
3096                   
3097                    //
3098                    // Activation function
3099                    //
3100                    dedf = network.derror[i];
3101                    dfdnet = network.dfdnet[i];
3102                    derror[network.structinfo[offs+2]] = derror[network.structinfo[offs+2]]+dedf*dfdnet;
3103                }
3104                if( network.structinfo[offs+0]==0 )
3105                {
3106                   
3107                    //
3108                    // Adaptive summator
3109                    //
3110                    n1 = network.structinfo[offs+2];
3111                    n2 = n1+network.structinfo[offs+1]-1;
3112                    w1 = network.structinfo[offs+3];
3113                    w2 = w1+network.structinfo[offs+1]-1;
3114                    dedf = network.derror[i];
3115                    dfdnet = 1.0;
3116                    v = dedf*dfdnet;
3117                    i1_ = (n1) - (w1);
3118                    for(i_=w1; i_<=w2;i_++)
3119                    {
3120                        grad[i_] = v*neurons[i_+i1_];
3121                    }
3122                    i1_ = (w1) - (n1);
3123                    for(i_=n1; i_<=n2;i_++)
3124                    {
3125                        derror[i_] = derror[i_] + v*weights[i_+i1_];
3126                    }
3127                }
3128                if( network.structinfo[offs+0]<0 )
3129                {
3130                    bflag = false;
3131                    if( network.structinfo[offs+0]==-2 | network.structinfo[offs+0]==-3 | network.structinfo[offs+0]==-4 )
3132                    {
3133                       
3134                        //
3135                        // Special neuron type, no back-propagation required
3136                        //
3137                        bflag = true;
3138                    }
3139                    System.Diagnostics.Debug.Assert(bflag, "MLPInternalCalculateGradient: unknown neuron type!");
3140                }
3141            }
3142        }
3143
3144
3145        /*************************************************************************
3146        Internal subroutine, chunked gradient
3147        *************************************************************************/
3148        private static void mlpchunkedgradient(ref multilayerperceptron network,
3149            ref double[,] xy,
3150            int cstart,
3151            int csize,
3152            ref double e,
3153            ref double[] grad,
3154            bool naturalerrorfunc)
3155        {
3156            int i = 0;
3157            int j = 0;
3158            int k = 0;
3159            int kl = 0;
3160            int n1 = 0;
3161            int n2 = 0;
3162            int w1 = 0;
3163            int w2 = 0;
3164            int c1 = 0;
3165            int c2 = 0;
3166            int ntotal = 0;
3167            int nin = 0;
3168            int nout = 0;
3169            int offs = 0;
3170            double f = 0;
3171            double df = 0;
3172            double d2f = 0;
3173            double v = 0;
3174            double s = 0;
3175            double fown = 0;
3176            double deown = 0;
3177            double net = 0;
3178            double lnnet = 0;
3179            double mx = 0;
3180            bool bflag = new bool();
3181            int istart = 0;
3182            int ineurons = 0;
3183            int idfdnet = 0;
3184            int iderror = 0;
3185            int izeros = 0;
3186            int i_ = 0;
3187            int i1_ = 0;
3188
3189           
3190            //
3191            // Read network geometry, prepare data
3192            //
3193            nin = network.structinfo[1];
3194            nout = network.structinfo[2];
3195            ntotal = network.structinfo[3];
3196            istart = network.structinfo[5];
3197            c1 = cstart;
3198            c2 = cstart+csize-1;
3199            ineurons = 0;
3200            idfdnet = ntotal;
3201            iderror = 2*ntotal;
3202            izeros = 3*ntotal;
3203            for(j=0; j<=csize-1; j++)
3204            {
3205                network.chunks[izeros,j] = 0;
3206            }
3207           
3208            //
3209            // Forward pass:
3210            // 1. Load inputs from XY to Chunks[0:NIn-1,0:CSize-1]
3211            // 2. Forward pass
3212            //
3213            for(i=0; i<=nin-1; i++)
3214            {
3215                for(j=0; j<=csize-1; j++)
3216                {
3217                    if( (double)(network.columnsigmas[i])!=(double)(0) )
3218                    {
3219                        network.chunks[i,j] = (xy[c1+j,i]-network.columnmeans[i])/network.columnsigmas[i];
3220                    }
3221                    else
3222                    {
3223                        network.chunks[i,j] = xy[c1+j,i]-network.columnmeans[i];
3224                    }
3225                }
3226            }
3227            for(i=0; i<=ntotal-1; i++)
3228            {
3229                offs = istart+i*nfieldwidth;
3230                if( network.structinfo[offs+0]>0 )
3231                {
3232                   
3233                    //
3234                    // Activation function:
3235                    // * calculate F vector, F(i) = F(NET(i))
3236                    //
3237                    n1 = network.structinfo[offs+2];
3238                    for(i_=0; i_<=csize-1;i_++)
3239                    {
3240                        network.chunks[i,i_] = network.chunks[n1,i_];
3241                    }
3242                    for(j=0; j<=csize-1; j++)
3243                    {
3244                        mlpactivationfunction(network.chunks[i,j], network.structinfo[offs+0], ref f, ref df, ref d2f);
3245                        network.chunks[i,j] = f;
3246                        network.chunks[idfdnet+i,j] = df;
3247                    }
3248                }
3249                if( network.structinfo[offs+0]==0 )
3250                {
3251                   
3252                    //
3253                    // Adaptive summator:
3254                    // * calculate NET vector, NET(i) = SUM(W(j,i)*Neurons(j),j=N1..N2)
3255                    //
3256                    n1 = network.structinfo[offs+2];
3257                    n2 = n1+network.structinfo[offs+1]-1;
3258                    w1 = network.structinfo[offs+3];
3259                    w2 = w1+network.structinfo[offs+1]-1;
3260                    for(i_=0; i_<=csize-1;i_++)
3261                    {
3262                        network.chunks[i,i_] = network.chunks[izeros,i_];
3263                    }
3264                    for(j=n1; j<=n2; j++)
3265                    {
3266                        v = network.weights[w1+j-n1];
3267                        for(i_=0; i_<=csize-1;i_++)
3268                        {
3269                            network.chunks[i,i_] = network.chunks[i,i_] + v*network.chunks[j,i_];
3270                        }
3271                    }
3272                }
3273                if( network.structinfo[offs+0]<0 )
3274                {
3275                    bflag = false;
3276                    if( network.structinfo[offs+0]==-2 )
3277                    {
3278                       
3279                        //
3280                        // input neuron, left unchanged
3281                        //
3282                        bflag = true;
3283                    }
3284                    if( network.structinfo[offs+0]==-3 )
3285                    {
3286                       
3287                        //
3288                        // "-1" neuron
3289                        //
3290                        for(k=0; k<=csize-1; k++)
3291                        {
3292                            network.chunks[i,k] = -1;
3293                        }
3294                        bflag = true;
3295                    }
3296                    if( network.structinfo[offs+0]==-4 )
3297                    {
3298                       
3299                        //
3300                        // "0" neuron
3301                        //
3302                        for(k=0; k<=csize-1; k++)
3303                        {
3304                            network.chunks[i,k] = 0;
3305                        }
3306                        bflag = true;
3307                    }
3308                    System.Diagnostics.Debug.Assert(bflag, "MLPChunkedGradient: internal error - unknown neuron type!");
3309                }
3310            }
3311           
3312            //
3313            // Post-processing, error, dError/dOut
3314            //
3315            for(i=0; i<=ntotal-1; i++)
3316            {
3317                for(i_=0; i_<=csize-1;i_++)
3318                {
3319                    network.chunks[iderror+i,i_] = network.chunks[izeros,i_];
3320                }
3321            }
3322            System.Diagnostics.Debug.Assert(network.structinfo[6]==0 | network.structinfo[6]==1, "MLPChunkedGradient: unknown normalization type!");
3323            if( network.structinfo[6]==1 )
3324            {
3325               
3326                //
3327                // Softmax output, classification network.
3328                //
3329                // For each K = 0..CSize-1 do:
3330                // 1. place exp(outputs[k]) to NWBuf[0:NOut-1]
3331                // 2. place sum(exp(..)) to NET
3332                // 3. calculate dError/dOut and place it to the second block of Chunks
3333                //
3334                for(k=0; k<=csize-1; k++)
3335                {
3336                   
3337                    //
3338                    // Normalize
3339                    //
3340                    mx = network.chunks[ntotal-nout,k];
3341                    for(i=1; i<=nout-1; i++)
3342                    {
3343                        mx = Math.Max(mx, network.chunks[ntotal-nout+i,k]);
3344                    }
3345                    net = 0;
3346                    for(i=0; i<=nout-1; i++)
3347                    {
3348                        network.nwbuf[i] = Math.Exp(network.chunks[ntotal-nout+i,k]-mx);
3349                        net = net+network.nwbuf[i];
3350                    }
3351                   
3352                    //
3353                    // Calculate error function and dError/dOut
3354                    //
3355                    if( naturalerrorfunc )
3356                    {
3357                       
3358                        //
3359                        // Natural error func.
3360                        //
3361                        //
3362                        s = 1;
3363                        lnnet = Math.Log(net);
3364                        kl = (int)Math.Round(xy[cstart+k,nin]);
3365                        for(i=0; i<=nout-1; i++)
3366                        {
3367                            if( i==kl )
3368                            {
3369                                v = 1;
3370                            }
3371                            else
3372                            {
3373                                v = 0;
3374                            }
3375                            network.chunks[iderror+ntotal-nout+i,k] = s*network.nwbuf[i]/net-v;
3376                            e = e+safecrossentropy(v, network.nwbuf[i]/net);
3377                        }
3378                    }
3379                    else
3380                    {
3381                       
3382                        //
3383                        // Least squares error func
3384                        // Error, dError/dOut(normalized)
3385                        //
3386                        kl = (int)Math.Round(xy[cstart+k,nin]);
3387                        for(i=0; i<=nout-1; i++)
3388                        {
3389                            if( i==kl )
3390                            {
3391                                v = network.nwbuf[i]/net-1;
3392                            }
3393                            else
3394                            {
3395                                v = network.nwbuf[i]/net;
3396                            }
3397                            network.nwbuf[nout+i] = v;
3398                            e = e+AP.Math.Sqr(v)/2;
3399                        }
3400                       
3401                        //
3402                        // From dError/dOut(normalized) to dError/dOut(non-normalized)
3403                        //
3404                        i1_ = (0)-(nout);
3405                        v = 0.0;
3406                        for(i_=nout; i_<=2*nout-1;i_++)
3407                        {
3408                            v += network.nwbuf[i_]*network.nwbuf[i_+i1_];
3409                        }
3410                        for(i=0; i<=nout-1; i++)
3411                        {
3412                            fown = network.nwbuf[i];
3413                            deown = network.nwbuf[nout+i];
3414                            network.chunks[iderror+ntotal-nout+i,k] = (-v+deown*fown+deown*(net-fown))*fown/AP.Math.Sqr(net);
3415                        }
3416                    }
3417                }
3418            }
3419            else
3420            {
3421               
3422                //
3423                // Normal output, regression network
3424                //
3425                // For each K = 0..CSize-1 do:
3426                // 1. calculate dError/dOut and place it to the second block of Chunks
3427                //
3428                for(i=0; i<=nout-1; i++)
3429                {
3430                    for(j=0; j<=csize-1; j++)
3431                    {
3432                        v = network.chunks[ntotal-nout+i,j]*network.columnsigmas[nin+i]+network.columnmeans[nin+i]-xy[cstart+j,nin+i];
3433                        network.chunks[iderror+ntotal-nout+i,j] = v*network.columnsigmas[nin+i];
3434                        e = e+AP.Math.Sqr(v)/2;
3435                    }
3436                }
3437            }
3438           
3439            //
3440            // Backpropagation
3441            //
3442            for(i=ntotal-1; i>=0; i--)
3443            {
3444               
3445                //
3446                // Extract info
3447                //
3448                offs = istart+i*nfieldwidth;
3449                if( network.structinfo[offs+0]>0 )
3450                {
3451                   
3452                    //
3453                    // Activation function
3454                    //
3455                    n1 = network.structinfo[offs+2];
3456                    for(k=0; k<=csize-1; k++)
3457                    {
3458                        network.chunks[iderror+i,k] = network.chunks[iderror+i,k]*network.chunks[idfdnet+i,k];
3459                    }
3460                    for(i_=0; i_<=csize-1;i_++)
3461                    {
3462                        network.chunks[iderror+n1,i_] = network.chunks[iderror+n1,i_] + network.chunks[iderror+i,i_];
3463                    }
3464                }
3465                if( network.structinfo[offs+0]==0 )
3466                {
3467                   
3468                    //
3469                    // "Normal" activation function
3470                    //
3471                    n1 = network.structinfo[offs+2];
3472                    n2 = n1+network.structinfo[offs+1]-1;
3473                    w1 = network.structinfo[offs+3];
3474                    w2 = w1+network.structinfo[offs+1]-1;
3475                    for(j=w1; j<=w2; j++)
3476                    {
3477                        v = 0.0;
3478                        for(i_=0; i_<=csize-1;i_++)
3479                        {
3480                            v += network.chunks[n1+j-w1,i_]*network.chunks[iderror+i,i_];
3481                        }
3482                        grad[j] = grad[j]+v;
3483                    }
3484                    for(j=n1; j<=n2; j++)
3485                    {
3486                        v = network.weights[w1+j-n1];
3487                        for(i_=0; i_<=csize-1;i_++)
3488                        {
3489                            network.chunks[iderror+j,i_] = network.chunks[iderror+j,i_] + v*network.chunks[iderror+i,i_];
3490                        }
3491                    }
3492                }
3493                if( network.structinfo[offs+0]<0 )
3494                {
3495                    bflag = false;
3496                    if( network.structinfo[offs+0]==-2 | network.structinfo[offs+0]==-3 | network.structinfo[offs+0]==-4 )
3497                    {
3498                       
3499                        //
3500                        // Special neuron type, no back-propagation required
3501                        //
3502                        bflag = true;
3503                    }
3504                    System.Diagnostics.Debug.Assert(bflag, "MLPInternalCalculateGradient: unknown neuron type!");
3505                }
3506            }
3507        }
3508
3509
3510        /*************************************************************************
3511        Returns T*Ln(T/Z), guarded against overflow/underflow.
3512        Internal subroutine.
3513        *************************************************************************/
3514        private static double safecrossentropy(double t,
3515            double z)
3516        {
3517            double result = 0;
3518            double r = 0;
3519
3520            if( (double)(t)==(double)(0) )
3521            {
3522                result = 0;
3523            }
3524            else
3525            {
3526                if( (double)(Math.Abs(z))>(double)(1) )
3527                {
3528                   
3529                    //
3530                    // Shouldn't be the case with softmax,
3531                    // but we just want to be sure.
3532                    //
3533                    if( (double)(t/z)==(double)(0) )
3534                    {
3535                        r = AP.Math.MinRealNumber;
3536                    }
3537                    else
3538                    {
3539                        r = t/z;
3540                    }
3541                }
3542                else
3543                {
3544                   
3545                    //
3546                    // Normal case
3547                    //
3548                    if( (double)(z)==(double)(0) | (double)(Math.Abs(t))>=(double)(AP.Math.MaxRealNumber*Math.Abs(z)) )
3549                    {
3550                        r = AP.Math.MaxRealNumber;
3551                    }
3552                    else
3553                    {
3554                        r = t/z;
3555                    }
3556                }
3557                result = t*Math.Log(r);
3558            }
3559            return result;
3560        }
3561    }
3562}
Note: See TracBrowser for help on using the repository browser.