Context Navigation

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/logit.cs @ 5229

Visit:

Last change on this file since 5229 was 3839, checked in by mkommend, 14 years ago
implemented first version of LR (ticket #1012)
File size: 60.8 KB

Line
1	/*************************************************************************
2	Copyright (c) 2008, Sergey Bochkanov (ALGLIB project).
3
4	>>> SOURCE LICENSE >>>
5	This program is free software; you can redistribute it and/or modify
6	it under the terms of the GNU General Public License as published by
7	the Free Software Foundation (www.fsf.org); either version 2 of the
8	License, or (at your option) any later version.
9
10	This program is distributed in the hope that it will be useful,
11	but WITHOUT ANY WARRANTY; without even the implied warranty of
12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13	GNU General Public License for more details.
14
15	A copy of the GNU General Public License is available at
16	http://www.fsf.org/licensing/licenses
17
18	>>> END OF LICENSE >>>
19	*************************************************************************/
20
21	using System;
22
23	namespace alglib
24	{
25	public class logit
26	{
27	public struct logitmodel
28	{
29	public double[] w;
30	};
31
32
33	public struct logitmcstate
34	{
35	public bool brackt;
36	public bool stage1;
37	public int infoc;
38	public double dg;
39	public double dgm;
40	public double dginit;
41	public double dgtest;
42	public double dgx;
43	public double dgxm;
44	public double dgy;
45	public double dgym;
46	public double finit;
47	public double ftest1;
48	public double fm;
49	public double fx;
50	public double fxm;
51	public double fy;
52	public double fym;
53	public double stx;
54	public double sty;
55	public double stmin;
56	public double stmax;
57	public double width;
58	public double width1;
59	public double xtrapf;
60	};
61
62
63	/*************************************************************************
64	MNLReport structure contains information about training process:
65	* NGrad - number of gradient calculations
66	* NHess - number of Hessian calculations
67	*************************************************************************/
68	public struct mnlreport
69	{
70	public int ngrad;
71	public int nhess;
72	};
73
74
75
76
77	public const double xtol = 100*AP.Math.MachineEpsilon;
78	public const double ftol = 0.0001;
79	public const double gtol = 0.3;
80	public const int maxfev = 20;
81	public const double stpmin = 1.0E-2;
82	public const double stpmax = 1.0E5;
83	public const int logitvnum = 6;
84
85
86	/*************************************************************************
87	This subroutine trains logit model.
88
89	INPUT PARAMETERS:
90	XY - training set, array[0..NPoints-1,0..NVars]
91	First NVars columns store values of independent
92	variables, next column stores number of class (from 0
93	to NClasses-1) which dataset element belongs to. Fractional
94	values are rounded to nearest integer.
95	NPoints - training set size, NPoints>=1
96	NVars - number of independent variables, NVars>=1
97	NClasses - number of classes, NClasses>=2
98
99	OUTPUT PARAMETERS:
100	Info - return code:
101	* -2, if there is a point with class number
102	outside of [0..NClasses-1].
103	* -1, if incorrect parameters was passed
104	(NPoints<NVars+2, NVars<1, NClasses<2).
105	* 1, if task has been solved
106	LM - model built
107	Rep - training report
108
109	-- ALGLIB --
110	Copyright 10.09.2008 by Bochkanov Sergey
111	*************************************************************************/
112	public static void mnltrainh(ref double[,] xy,
113	int npoints,
114	int nvars,
115	int nclasses,
116	ref int info,
117	ref logitmodel lm,
118	ref mnlreport rep)
119	{
120	int i = 0;
121	int j = 0;
122	int k = 0;
123	int ssize = 0;
124	bool allsame = new bool();
125	int offs = 0;
126	double threshold = 0;
127	double wminstep = 0;
128	double decay = 0;
129	int wdim = 0;
130	int expoffs = 0;
131	double v = 0;
132	double s = 0;
133	mlpbase.multilayerperceptron network = new mlpbase.multilayerperceptron();
134	int nin = 0;
135	int nout = 0;
136	int wcount = 0;
137	double e = 0;
138	double[] g = new double[0];
139	double[,] h = new double[0,0];
140	bool spd = new bool();
141	double[] x = new double[0];
142	double[] y = new double[0];
143	double[] wbase = new double[0];
144	double wstep = 0;
145	double[] wdir = new double[0];
146	double[] work = new double[0];
147	int mcstage = 0;
148	logitmcstate mcstate = new logitmcstate();
149	int mcinfo = 0;
150	int mcnfev = 0;
151	int solverinfo = 0;
152	densesolver.densesolverreport solverrep = new densesolver.densesolverreport();
153	int i_ = 0;
154	int i1_ = 0;
155
156	threshold = 1000*AP.Math.MachineEpsilon;
157	wminstep = 0.001;
158	decay = 0.001;
159
160	//
161	// Test for inputs
162	//
163	if( npoints<nvars+2 \| nvars<1 \| nclasses<2 )
164	{
165	info = -1;
166	return;
167	}
168	for(i=0; i<=npoints-1; i++)
169	{
170	if( (int)Math.Round(xy[i,nvars])<0 \| (int)Math.Round(xy[i,nvars])>=nclasses )
171	{
172	info = -2;
173	return;
174	}
175	}
176	info = 1;
177
178	//
179	// Initialize data
180	//
181	rep.ngrad = 0;
182	rep.nhess = 0;
183
184	//
185	// Allocate array
186	//
187	wdim = (nvars+1)*(nclasses-1);
188	offs = 5;
189	expoffs = offs+wdim;
190	ssize = 5+(nvars+1)*(nclasses-1)+nclasses;
191	lm.w = new double[ssize-1+1];
192	lm.w[0] = ssize;
193	lm.w[1] = logitvnum;
194	lm.w[2] = nvars;
195	lm.w[3] = nclasses;
196	lm.w[4] = offs;
197
198	//
199	// Degenerate case: all outputs are equal
200	//
201	allsame = true;
202	for(i=1; i<=npoints-1; i++)
203	{
204	if( (int)Math.Round(xy[i,nvars])!=(int)Math.Round(xy[i-1,nvars]) )
205	{
206	allsame = false;
207	}
208	}
209	if( allsame )
210	{
211	for(i=0; i<=(nvars+1)*(nclasses-1)-1; i++)
212	{
213	lm.w[offs+i] = 0;
214	}
215	v = -(2*Math.Log(AP.Math.MinRealNumber));
216	k = (int)Math.Round(xy[0,nvars]);
217	if( k==nclasses-1 )
218	{
219	for(i=0; i<=nclasses-2; i++)
220	{
221	lm.w[offs+i*(nvars+1)+nvars] = -v;
222	}
223	}
224	else
225	{
226	for(i=0; i<=nclasses-2; i++)
227	{
228	if( i==k )
229	{
230	lm.w[offs+i*(nvars+1)+nvars] = +v;
231	}
232	else
233	{
234	lm.w[offs+i*(nvars+1)+nvars] = 0;
235	}
236	}
237	}
238	return;
239	}
240
241	//
242	// General case.
243	// Prepare task and network. Allocate space.
244	//
245	mlpbase.mlpcreatec0(nvars, nclasses, ref network);
246	mlpbase.mlpinitpreprocessor(ref network, ref xy, npoints);
247	mlpbase.mlpproperties(ref network, ref nin, ref nout, ref wcount);
248	for(i=0; i<=wcount-1; i++)
249	{
250	network.weights[i] = (2*AP.Math.RandomReal()-1)/nvars;
251	}
252	g = new double[wcount-1+1];
253	h = new double[wcount-1+1, wcount-1+1];
254	wbase = new double[wcount-1+1];
255	wdir = new double[wcount-1+1];
256	work = new double[wcount-1+1];
257
258	//
259	// First stage: optimize in gradient direction.
260	//
261	for(k=0; k<=wcount/3+10; k++)
262	{
263
264	//
265	// Calculate gradient in starting point
266	//
267	mlpbase.mlpgradnbatch(ref network, ref xy, npoints, ref e, ref g);
268	v = 0.0;
269	for(i_=0; i_<=wcount-1;i_++)
270	{
271	v += network.weights[i_]*network.weights[i_];
272	}
273	e = e+0.5decayv;
274	for(i_=0; i_<=wcount-1;i_++)
275	{
276	g[i_] = g[i_] + decay*network.weights[i_];
277	}
278	rep.ngrad = rep.ngrad+1;
279
280	//
281	// Setup optimization scheme
282	//
283	for(i_=0; i_<=wcount-1;i_++)
284	{
285	wdir[i_] = -g[i_];
286	}
287	v = 0.0;
288	for(i_=0; i_<=wcount-1;i_++)
289	{
290	v += wdir[i_]*wdir[i_];
291	}
292	wstep = Math.Sqrt(v);
293	v = 1/Math.Sqrt(v);
294	for(i_=0; i_<=wcount-1;i_++)
295	{
296	wdir[i_] = v*wdir[i_];
297	}
298	mcstage = 0;
299	mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
300	while( mcstage!=0 )
301	{
302	mlpbase.mlpgradnbatch(ref network, ref xy, npoints, ref e, ref g);
303	v = 0.0;
304	for(i_=0; i_<=wcount-1;i_++)
305	{
306	v += network.weights[i_]*network.weights[i_];
307	}
308	e = e+0.5decayv;
309	for(i_=0; i_<=wcount-1;i_++)
310	{
311	g[i_] = g[i_] + decay*network.weights[i_];
312	}
313	rep.ngrad = rep.ngrad+1;
314	mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
315	}
316	}
317
318	//
319	// Second stage: use Hessian when we are close to the minimum
320	//
321	while( true )
322	{
323
324	//
325	// Calculate and update E/G/H
326	//
327	mlpbase.mlphessiannbatch(ref network, ref xy, npoints, ref e, ref g, ref h);
328	v = 0.0;
329	for(i_=0; i_<=wcount-1;i_++)
330	{
331	v += network.weights[i_]*network.weights[i_];
332	}
333	e = e+0.5decayv;
334	for(i_=0; i_<=wcount-1;i_++)
335	{
336	g[i_] = g[i_] + decay*network.weights[i_];
337	}
338	for(k=0; k<=wcount-1; k++)
339	{
340	h[k,k] = h[k,k]+decay;
341	}
342	rep.nhess = rep.nhess+1;
343
344	//
345	// Select step direction
346	// NOTE: it is important to use lower-triangle Cholesky
347	// factorization since it is much faster than higher-triangle version.
348	//
349	spd = trfac.spdmatrixcholesky(ref h, wcount, false);
350	densesolver.spdmatrixcholeskysolve(ref h, wcount, false, ref g, ref solverinfo, ref solverrep, ref wdir);
351	spd = solverinfo>0;
352	if( spd )
353	{
354
355	//
356	// H is positive definite.
357	// Step in Newton direction.
358	//
359	for(i_=0; i_<=wcount-1;i_++)
360	{
361	wdir[i_] = -1*wdir[i_];
362	}
363	spd = true;
364	}
365	else
366	{
367
368	//
369	// H is indefinite.
370	// Step in gradient direction.
371	//
372	for(i_=0; i_<=wcount-1;i_++)
373	{
374	wdir[i_] = -g[i_];
375	}
376	spd = false;
377	}
378
379	//
380	// Optimize in WDir direction
381	//
382	v = 0.0;
383	for(i_=0; i_<=wcount-1;i_++)
384	{
385	v += wdir[i_]*wdir[i_];
386	}
387	wstep = Math.Sqrt(v);
388	v = 1/Math.Sqrt(v);
389	for(i_=0; i_<=wcount-1;i_++)
390	{
391	wdir[i_] = v*wdir[i_];
392	}
393	mcstage = 0;
394	mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
395	while( mcstage!=0 )
396	{
397	mlpbase.mlpgradnbatch(ref network, ref xy, npoints, ref e, ref g);
398	v = 0.0;
399	for(i_=0; i_<=wcount-1;i_++)
400	{
401	v += network.weights[i_]*network.weights[i_];
402	}
403	e = e+0.5decayv;
404	for(i_=0; i_<=wcount-1;i_++)
405	{
406	g[i_] = g[i_] + decay*network.weights[i_];
407	}
408	rep.ngrad = rep.ngrad+1;
409	mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
410	}
411	if( spd & (mcinfo==2 \| mcinfo==4 \| mcinfo==6) )
412	{
413	break;
414	}
415	}
416
417	//
418	// Convert from NN format to MNL format
419	//
420	i1_ = (0) - (offs);
421	for(i_=offs; i_<=offs+wcount-1;i_++)
422	{
423	lm.w[i_] = network.weights[i_+i1_];
424	}
425	for(k=0; k<=nvars-1; k++)
426	{
427	for(i=0; i<=nclasses-2; i++)
428	{
429	s = network.columnsigmas[k];
430	if( (double)(s)==(double)(0) )
431	{
432	s = 1;
433	}
434	j = offs+(nvars+1)*i;
435	v = lm.w[j+k];
436	lm.w[j+k] = v/s;
437	lm.w[j+nvars] = lm.w[j+nvars]+v*network.columnmeans[k]/s;
438	}
439	}
440	for(k=0; k<=nclasses-2; k++)
441	{
442	lm.w[offs+(nvars+1)k+nvars] = -lm.w[offs+(nvars+1)k+nvars];
443	}
444	}
445
446
447	/*************************************************************************
448	Procesing
449
450	INPUT PARAMETERS:
451	LM - logit model, passed by non-constant reference
452	(some fields of structure are used as temporaries
453	when calculating model output).
454	X - input vector, array[0..NVars-1].
455
456	OUTPUT PARAMETERS:
457	Y - result, array[0..NClasses-1]
458	Vector of posterior probabilities for classification task.
459	Subroutine does not allocate memory for this vector, it is
460	responsibility of a caller to allocate it. Array must be
461	at least [0..NClasses-1].
462
463	-- ALGLIB --
464	Copyright 10.09.2008 by Bochkanov Sergey
465	*************************************************************************/
466	public static void mnlprocess(ref logitmodel lm,
467	ref double[] x,
468	ref double[] y)
469	{
470	int nvars = 0;
471	int nclasses = 0;
472	int offs = 0;
473	int i = 0;
474	int i1 = 0;
475	double s = 0;
476
477	System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLProcess: unexpected model version");
478	nvars = (int)Math.Round(lm.w[2]);
479	nclasses = (int)Math.Round(lm.w[3]);
480	offs = (int)Math.Round(lm.w[4]);
481	mnliexp(ref lm.w, ref x);
482	s = 0;
483	i1 = offs+(nvars+1)*(nclasses-1);
484	for(i=i1; i<=i1+nclasses-1; i++)
485	{
486	s = s+lm.w[i];
487	}
488	for(i=0; i<=nclasses-1; i++)
489	{
490	y[i] = lm.w[i1+i]/s;
491	}
492	}
493
494
495	/*************************************************************************
496	Unpacks coefficients of logit model. Logit model have form:
497
498	P(class=i) = S(i) / (S(0) + S(1) + ... +S(M-1))
499	S(i) = Exp(A[i,0]X[0] + ... + A[i,N-1]X[N-1] + A[i,N]), when i<M-1
500	S(M-1) = 1
501
502	INPUT PARAMETERS:
503	LM - logit model in ALGLIB format
504
505	OUTPUT PARAMETERS:
506	V - coefficients, array[0..NClasses-2,0..NVars]
507	NVars - number of independent variables
508	NClasses - number of classes
509
510	-- ALGLIB --
511	Copyright 10.09.2008 by Bochkanov Sergey
512	*************************************************************************/
513	public static void mnlunpack(ref logitmodel lm,
514	ref double[,] a,
515	ref int nvars,
516	ref int nclasses)
517	{
518	int offs = 0;
519	int i = 0;
520	int i_ = 0;
521	int i1_ = 0;
522
523	System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLUnpack: unexpected model version");
524	nvars = (int)Math.Round(lm.w[2]);
525	nclasses = (int)Math.Round(lm.w[3]);
526	offs = (int)Math.Round(lm.w[4]);
527	a = new double[nclasses-2+1, nvars+1];
528	for(i=0; i<=nclasses-2; i++)
529	{
530	i1_ = (offs+i*(nvars+1)) - (0);
531	for(i_=0; i_<=nvars;i_++)
532	{
533	a[i,i_] = lm.w[i_+i1_];
534	}
535	}
536	}
537
538
539	/*************************************************************************
540	"Packs" coefficients and creates logit model in ALGLIB format (MNLUnpack
541	reversed).
542
543	INPUT PARAMETERS:
544	A - model (see MNLUnpack)
545	NVars - number of independent variables
546	NClasses - number of classes
547
548	OUTPUT PARAMETERS:
549	LM - logit model.
550
551	-- ALGLIB --
552	Copyright 10.09.2008 by Bochkanov Sergey
553	*************************************************************************/
554	public static void mnlpack(ref double[,] a,
555	int nvars,
556	int nclasses,
557	ref logitmodel lm)
558	{
559	int offs = 0;
560	int i = 0;
561	int wdim = 0;
562	int ssize = 0;
563	int i_ = 0;
564	int i1_ = 0;
565
566	wdim = (nvars+1)*(nclasses-1);
567	offs = 5;
568	ssize = 5+(nvars+1)*(nclasses-1)+nclasses;
569	lm.w = new double[ssize-1+1];
570	lm.w[0] = ssize;
571	lm.w[1] = logitvnum;
572	lm.w[2] = nvars;
573	lm.w[3] = nclasses;
574	lm.w[4] = offs;
575	for(i=0; i<=nclasses-2; i++)
576	{
577	i1_ = (0) - (offs+i*(nvars+1));
578	for(i_=offs+i(nvars+1); i_<=offs+i(nvars+1)+nvars;i_++)
579	{
580	lm.w[i_] = a[i,i_+i1_];
581	}
582	}
583	}
584
585
586	/*************************************************************************
587	Copying of LogitModel strucure
588
589	INPUT PARAMETERS:
590	LM1 - original
591
592	OUTPUT PARAMETERS:
593	LM2 - copy
594
595	-- ALGLIB --
596	Copyright 15.03.2009 by Bochkanov Sergey
597	*************************************************************************/
598	public static void mnlcopy(ref logitmodel lm1,
599	ref logitmodel lm2)
600	{
601	int k = 0;
602	int i_ = 0;
603
604	k = (int)Math.Round(lm1.w[0]);
605	lm2.w = new double[k-1+1];
606	for(i_=0; i_<=k-1;i_++)
607	{
608	lm2.w[i_] = lm1.w[i_];
609	}
610	}
611
612
613	/*************************************************************************
614	Serialization of LogitModel strucure
615
616	INPUT PARAMETERS:
617	LM - original
618
619	OUTPUT PARAMETERS:
620	RA - array of real numbers which stores model,
621	array[0..RLen-1]
622	RLen - RA lenght
623
624	-- ALGLIB --
625	Copyright 15.03.2009 by Bochkanov Sergey
626	*************************************************************************/
627	public static void mnlserialize(ref logitmodel lm,
628	ref double[] ra,
629	ref int rlen)
630	{
631	int i_ = 0;
632	int i1_ = 0;
633
634	rlen = (int)Math.Round(lm.w[0])+1;
635	ra = new double[rlen-1+1];
636	ra[0] = logitvnum;
637	i1_ = (0) - (1);
638	for(i_=1; i_<=rlen-1;i_++)
639	{
640	ra[i_] = lm.w[i_+i1_];
641	}
642	}
643
644
645	/*************************************************************************
646	Unserialization of LogitModel strucure
647
648	INPUT PARAMETERS:
649	RA - real array which stores model
650
651	OUTPUT PARAMETERS:
652	LM - restored model
653
654	-- ALGLIB --
655	Copyright 15.03.2009 by Bochkanov Sergey
656	*************************************************************************/
657	public static void mnlunserialize(ref double[] ra,
658	ref logitmodel lm)
659	{
660	int i_ = 0;
661	int i1_ = 0;
662
663	System.Diagnostics.Debug.Assert((int)Math.Round(ra[0])==logitvnum, "MNLUnserialize: incorrect array!");
664	lm.w = new double[(int)Math.Round(ra[1])-1+1];
665	i1_ = (1) - (0);
666	for(i_=0; i_<=(int)Math.Round(ra[1])-1;i_++)
667	{
668	lm.w[i_] = ra[i_+i1_];
669	}
670	}
671
672
673	/*************************************************************************
674	Average cross-entropy (in bits per element) on the test set
675
676	INPUT PARAMETERS:
677	LM - logit model
678	XY - test set
679	NPoints - test set size
680
681	RESULT:
682	CrossEntropy/(NPoints*ln(2)).
683
684	-- ALGLIB --
685	Copyright 10.09.2008 by Bochkanov Sergey
686	*************************************************************************/
687	public static double mnlavgce(ref logitmodel lm,
688	ref double[,] xy,
689	int npoints)
690	{
691	double result = 0;
692	int nvars = 0;
693	int nclasses = 0;
694	int i = 0;
695	double[] workx = new double[0];
696	double[] worky = new double[0];
697	int i_ = 0;
698
699	System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLClsError: unexpected model version");
700	nvars = (int)Math.Round(lm.w[2]);
701	nclasses = (int)Math.Round(lm.w[3]);
702	workx = new double[nvars-1+1];
703	worky = new double[nclasses-1+1];
704	result = 0;
705	for(i=0; i<=npoints-1; i++)
706	{
707	System.Diagnostics.Debug.Assert((int)Math.Round(xy[i,nvars])>=0 & (int)Math.Round(xy[i,nvars])<nclasses, "MNLAvgCE: incorrect class number!");
708
709	//
710	// Process
711	//
712	for(i_=0; i_<=nvars-1;i_++)
713	{
714	workx[i_] = xy[i,i_];
715	}
716	mnlprocess(ref lm, ref workx, ref worky);
717	if( (double)(worky[(int)Math.Round(xy[i,nvars])])>(double)(0) )
718	{
719	result = result-Math.Log(worky[(int)Math.Round(xy[i,nvars])]);
720	}
721	else
722	{
723	result = result-Math.Log(AP.Math.MinRealNumber);
724	}
725	}
726	result = result/(npoints*Math.Log(2));
727	return result;
728	}
729
730
731	/*************************************************************************
732	Relative classification error on the test set
733
734	INPUT PARAMETERS:
735	LM - logit model
736	XY - test set
737	NPoints - test set size
738
739	RESULT:
740	percent of incorrectly classified cases.
741
742	-- ALGLIB --
743	Copyright 10.09.2008 by Bochkanov Sergey
744	*************************************************************************/
745	public static double mnlrelclserror(ref logitmodel lm,
746	ref double[,] xy,
747	int npoints)
748	{
749	double result = 0;
750
751	result = (double)(mnlclserror(ref lm, ref xy, npoints))/(double)(npoints);
752	return result;
753	}
754
755
756	/*************************************************************************
757	RMS error on the test set
758
759	INPUT PARAMETERS:
760	LM - logit model
761	XY - test set
762	NPoints - test set size
763
764	RESULT:
765	root mean square error (error when estimating posterior probabilities).
766
767	-- ALGLIB --
768	Copyright 30.08.2008 by Bochkanov Sergey
769	*************************************************************************/
770	public static double mnlrmserror(ref logitmodel lm,
771	ref double[,] xy,
772	int npoints)
773	{
774	double result = 0;
775	double relcls = 0;
776	double avgce = 0;
777	double rms = 0;
778	double avg = 0;
779	double avgrel = 0;
780
781	System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNLRMSError: Incorrect MNL version!");
782	mnlallerrors(ref lm, ref xy, npoints, ref relcls, ref avgce, ref rms, ref avg, ref avgrel);
783	result = rms;
784	return result;
785	}
786
787
788	/*************************************************************************
789	Average error on the test set
790
791	INPUT PARAMETERS:
792	LM - logit model
793	XY - test set
794	NPoints - test set size
795
796	RESULT:
797	average error (error when estimating posterior probabilities).
798
799	-- ALGLIB --
800	Copyright 30.08.2008 by Bochkanov Sergey
801	*************************************************************************/
802	public static double mnlavgerror(ref logitmodel lm,
803	ref double[,] xy,
804	int npoints)
805	{
806	double result = 0;
807	double relcls = 0;
808	double avgce = 0;
809	double rms = 0;
810	double avg = 0;
811	double avgrel = 0;
812
813	System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNLRMSError: Incorrect MNL version!");
814	mnlallerrors(ref lm, ref xy, npoints, ref relcls, ref avgce, ref rms, ref avg, ref avgrel);
815	result = avg;
816	return result;
817	}
818
819
820	/*************************************************************************
821	Average relative error on the test set
822
823	INPUT PARAMETERS:
824	LM - logit model
825	XY - test set
826	NPoints - test set size
827
828	RESULT:
829	average relative error (error when estimating posterior probabilities).
830
831	-- ALGLIB --
832	Copyright 30.08.2008 by Bochkanov Sergey
833	*************************************************************************/
834	public static double mnlavgrelerror(ref logitmodel lm,
835	ref double[,] xy,
836	int ssize)
837	{
838	double result = 0;
839	double relcls = 0;
840	double avgce = 0;
841	double rms = 0;
842	double avg = 0;
843	double avgrel = 0;
844
845	System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNLRMSError: Incorrect MNL version!");
846	mnlallerrors(ref lm, ref xy, ssize, ref relcls, ref avgce, ref rms, ref avg, ref avgrel);
847	result = avgrel;
848	return result;
849	}
850
851
852	/*************************************************************************
853	Classification error on test set = MNLRelClsError*NPoints
854
855	-- ALGLIB --
856	Copyright 10.09.2008 by Bochkanov Sergey
857	*************************************************************************/
858	public static int mnlclserror(ref logitmodel lm,
859	ref double[,] xy,
860	int npoints)
861	{
862	int result = 0;
863	int nvars = 0;
864	int nclasses = 0;
865	int i = 0;
866	int j = 0;
867	double[] workx = new double[0];
868	double[] worky = new double[0];
869	int nmax = 0;
870	int i_ = 0;
871
872	System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLClsError: unexpected model version");
873	nvars = (int)Math.Round(lm.w[2]);
874	nclasses = (int)Math.Round(lm.w[3]);
875	workx = new double[nvars-1+1];
876	worky = new double[nclasses-1+1];
877	result = 0;
878	for(i=0; i<=npoints-1; i++)
879	{
880
881	//
882	// Process
883	//
884	for(i_=0; i_<=nvars-1;i_++)
885	{
886	workx[i_] = xy[i,i_];
887	}
888	mnlprocess(ref lm, ref workx, ref worky);
889
890	//
891	// Logit version of the answer
892	//
893	nmax = 0;
894	for(j=0; j<=nclasses-1; j++)
895	{
896	if( (double)(worky[j])>(double)(worky[nmax]) )
897	{
898	nmax = j;
899	}
900	}
901
902	//
903	// compare
904	//
905	if( nmax!=(int)Math.Round(xy[i,nvars]) )
906	{
907	result = result+1;
908	}
909	}
910	return result;
911	}
912
913
914	/*************************************************************************
915	Internal subroutine. Places exponents of the anti-overflow shifted
916	internal linear outputs into the service part of the W array.
917	*************************************************************************/
918	private static void mnliexp(ref double[] w,
919	ref double[] x)
920	{
921	int nvars = 0;
922	int nclasses = 0;
923	int offs = 0;
924	int i = 0;
925	int i1 = 0;
926	double v = 0;
927	double mx = 0;
928	int i_ = 0;
929	int i1_ = 0;
930
931	System.Diagnostics.Debug.Assert((double)(w[1])==(double)(logitvnum), "LOGIT: unexpected model version");
932	nvars = (int)Math.Round(w[2]);
933	nclasses = (int)Math.Round(w[3]);
934	offs = (int)Math.Round(w[4]);
935	i1 = offs+(nvars+1)*(nclasses-1);
936	for(i=0; i<=nclasses-2; i++)
937	{
938	i1_ = (0)-(offs+i*(nvars+1));
939	v = 0.0;
940	for(i_=offs+i(nvars+1); i_<=offs+i(nvars+1)+nvars-1;i_++)
941	{
942	v += w[i_]*x[i_+i1_];
943	}
944	w[i1+i] = v+w[offs+i*(nvars+1)+nvars];
945	}
946	w[i1+nclasses-1] = 0;
947	mx = 0;
948	for(i=i1; i<=i1+nclasses-1; i++)
949	{
950	mx = Math.Max(mx, w[i]);
951	}
952	for(i=i1; i<=i1+nclasses-1; i++)
953	{
954	w[i] = Math.Exp(w[i]-mx);
955	}
956	}
957
958
959	/*************************************************************************
960	Calculation of all types of errors
961
962	-- ALGLIB --
963	Copyright 30.08.2008 by Bochkanov Sergey
964	*************************************************************************/
965	private static void mnlallerrors(ref logitmodel lm,
966	ref double[,] xy,
967	int npoints,
968	ref double relcls,
969	ref double avgce,
970	ref double rms,
971	ref double avg,
972	ref double avgrel)
973	{
974	int nvars = 0;
975	int nclasses = 0;
976	int i = 0;
977	double[] buf = new double[0];
978	double[] workx = new double[0];
979	double[] y = new double[0];
980	double[] dy = new double[0];
981	int i_ = 0;
982
983	System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNL unit: Incorrect MNL version!");
984	nvars = (int)Math.Round(lm.w[2]);
985	nclasses = (int)Math.Round(lm.w[3]);
986	workx = new double[nvars-1+1];
987	y = new double[nclasses-1+1];
988	dy = new double[0+1];
989	bdss.dserrallocate(nclasses, ref buf);
990	for(i=0; i<=npoints-1; i++)
991	{
992	for(i_=0; i_<=nvars-1;i_++)
993	{
994	workx[i_] = xy[i,i_];
995	}
996	mnlprocess(ref lm, ref workx, ref y);
997	dy[0] = xy[i,nvars];
998	bdss.dserraccumulate(ref buf, ref y, ref dy);
999	}
1000	bdss.dserrfinish(ref buf);
1001	relcls = buf[0];
1002	avgce = buf[1];
1003	rms = buf[2];
1004	avg = buf[3];
1005	avgrel = buf[4];
1006	}
1007
1008
1009	/*************************************************************************
1010	THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES A SUFFICIENT
1011	DECREASE CONDITION AND A CURVATURE CONDITION.
1012
1013	AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF UNCERTAINTY WITH
1014	ENDPOINTS STX AND STY. THE INTERVAL OF UNCERTAINTY IS INITIALLY CHOSEN
1015	SO THAT IT CONTAINS A MINIMIZER OF THE MODIFIED FUNCTION
1016
1017	F(X+STPS) - F(X) - FTOLSTP*(GRADF(X)'S).
1018
1019	IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE
1020	FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, THEN THE INTERVAL OF
1021	UNCERTAINTY IS CHOSEN SO THAT IT CONTAINS A MINIMIZER OF F(X+STP*S).
1022
1023	THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES THE SUFFICIENT
1024	DECREASE CONDITION
1025
1026	F(X+STPS) .LE. F(X) + FTOLSTP*(GRADF(X)'S),
1027
1028	AND THE CURVATURE CONDITION
1029
1030	ABS(GRADF(X+STPS)'S)) .LE. GTOLABS(GRADF(X)'S).
1031
1032	IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION IS BOUNDED
1033	BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES BOTH CONDITIONS.
1034	IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH CONDITIONS, THEN THE
1035	ALGORITHM USUALLY STOPS WHEN ROUNDING ERRORS PREVENT FURTHER PROGRESS.
1036	IN THIS CASE STP ONLY SATISFIES THE SUFFICIENT DECREASE CONDITION.
1037
1038	PARAMETERS DESCRIPRION
1039
1040	N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER OF VARIABLES.
1041
1042	X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE BASE POINT FOR
1043	THE LINE SEARCH. ON OUTPUT IT CONTAINS X+STP*S.
1044
1045	F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F AT X. ON OUTPUT
1046	IT CONTAINS THE VALUE OF F AT X + STP*S.
1047
1048	G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE GRADIENT OF F AT X.
1049	ON OUTPUT IT CONTAINS THE GRADIENT OF F AT X + STP*S.
1050
1051	S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE SEARCH DIRECTION.
1052
1053	STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN INITIAL ESTIMATE
1054	OF A SATISFACTORY STEP. ON OUTPUT STP CONTAINS THE FINAL ESTIMATE.
1055
1056	FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION OCCURS WHEN THE
1057	SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
1058	SATISFIED.
1059
1060	XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS WHEN THE RELATIVE
1061	WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL.
1062
1063	STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH SPECIFY LOWER AND
1064	UPPER BOUNDS FOR THE STEP.
1065
1066	MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION OCCURS WHEN THE
1067	NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV BY THE END OF AN ITERATION.
1068
1069	INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
1070	INFO = 0 IMPROPER INPUT PARAMETERS.
1071
1072	INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE
1073	DIRECTIONAL DERIVATIVE CONDITION HOLD.
1074
1075	INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
1076	IS AT MOST XTOL.
1077
1078	INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
1079
1080	INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN.
1081
1082	INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX.
1083
1084	INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
1085	THERE MAY NOT BE A STEP WHICH SATISFIES THE
1086	SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
1087	TOLERANCES MAY BE TOO SMALL.
1088
1089	NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF CALLS TO FCN.
1090
1091	WA IS A WORK ARRAY OF LENGTH N.
1092
1093	ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
1094	JORGE J. MORE', DAVID J. THUENTE
1095	*************************************************************************/
1096	private static void mnlmcsrch(int n,
1097	ref double[] x,
1098	ref double f,
1099	ref double[] g,
1100	ref double[] s,
1101	ref double stp,
1102	ref int info,
1103	ref int nfev,
1104	ref double[] wa,
1105	ref logitmcstate state,
1106	ref int stage)
1107	{
1108	double v = 0;
1109	double p5 = 0;
1110	double p66 = 0;
1111	double zero = 0;
1112	int i_ = 0;
1113
1114
1115	//
1116	// init
1117	//
1118	p5 = 0.5;
1119	p66 = 0.66;
1120	state.xtrapf = 4.0;
1121	zero = 0;
1122
1123	//
1124	// Main cycle
1125	//
1126	while( true )
1127	{
1128	if( stage==0 )
1129	{
1130
1131	//
1132	// NEXT
1133	//
1134	stage = 2;
1135	continue;
1136	}
1137	if( stage==2 )
1138	{
1139	state.infoc = 1;
1140	info = 0;
1141
1142	//
1143	// CHECK THE INPUT PARAMETERS FOR ERRORS.
1144	//
1145	if( n<=0 \| (double)(stp)<=(double)(0) \| (double)(ftol)<(double)(0) \| (double)(gtol)<(double)(zero) \| (double)(xtol)<(double)(zero) \| (double)(stpmin)<(double)(zero) \| (double)(stpmax)<(double)(stpmin) \| maxfev<=0 )
1146	{
1147	stage = 0;
1148	return;
1149	}
1150
1151	//
1152	// COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
1153	// AND CHECK THAT S IS A DESCENT DIRECTION.
1154	//
1155	v = 0.0;
1156	for(i_=0; i_<=n-1;i_++)
1157	{
1158	v += g[i_]*s[i_];
1159	}
1160	state.dginit = v;
1161	if( (double)(state.dginit)>=(double)(0) )
1162	{
1163	stage = 0;
1164	return;
1165	}
1166
1167	//
1168	// INITIALIZE LOCAL VARIABLES.
1169	//
1170	state.brackt = false;
1171	state.stage1 = true;
1172	nfev = 0;
1173	state.finit = f;
1174	state.dgtest = ftol*state.dginit;
1175	state.width = stpmax-stpmin;
1176	state.width1 = state.width/p5;
1177	for(i_=0; i_<=n-1;i_++)
1178	{
1179	wa[i_] = x[i_];
1180	}
1181
1182	//
1183	// THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
1184	// FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
1185	// THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
1186	// FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
1187	// THE INTERVAL OF UNCERTAINTY.
1188	// THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
1189	// FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
1190	//
1191	state.stx = 0;
1192	state.fx = state.finit;
1193	state.dgx = state.dginit;
1194	state.sty = 0;
1195	state.fy = state.finit;
1196	state.dgy = state.dginit;
1197
1198	//
1199	// NEXT
1200	//
1201	stage = 3;
1202	continue;
1203	}
1204	if( stage==3 )
1205	{
1206
1207	//
1208	// START OF ITERATION.
1209	//
1210	// SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
1211	// TO THE PRESENT INTERVAL OF UNCERTAINTY.
1212	//
1213	if( state.brackt )
1214	{
1215	if( (double)(state.stx)<(double)(state.sty) )
1216	{
1217	state.stmin = state.stx;
1218	state.stmax = state.sty;
1219	}
1220	else
1221	{
1222	state.stmin = state.sty;
1223	state.stmax = state.stx;
1224	}
1225	}
1226	else
1227	{
1228	state.stmin = state.stx;
1229	state.stmax = stp+state.xtrapf*(stp-state.stx);
1230	}
1231
1232	//
1233	// FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
1234	//
1235	if( (double)(stp)>(double)(stpmax) )
1236	{
1237	stp = stpmax;
1238	}
1239	if( (double)(stp)<(double)(stpmin) )
1240	{
1241	stp = stpmin;
1242	}
1243
1244	//
1245	// IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
1246	// STP BE THE LOWEST POINT OBTAINED SO FAR.
1247	//
1248	if( state.brackt & ((double)(stp)<=(double)(state.stmin) \| (double)(stp)>=(double)(state.stmax)) \| nfev>=maxfev-1 \| state.infoc==0 \| state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
1249	{
1250	stp = state.stx;
1251	}
1252
1253	//
1254	// EVALUATE THE FUNCTION AND GRADIENT AT STP
1255	// AND COMPUTE THE DIRECTIONAL DERIVATIVE.
1256	//
1257	for(i_=0; i_<=n-1;i_++)
1258	{
1259	x[i_] = wa[i_];
1260	}
1261	for(i_=0; i_<=n-1;i_++)
1262	{
1263	x[i_] = x[i_] + stp*s[i_];
1264	}
1265
1266	//
1267	// NEXT
1268	//
1269	stage = 4;
1270	return;
1271	}
1272	if( stage==4 )
1273	{
1274	info = 0;
1275	nfev = nfev+1;
1276	v = 0.0;
1277	for(i_=0; i_<=n-1;i_++)
1278	{
1279	v += g[i_]*s[i_];
1280	}
1281	state.dg = v;
1282	state.ftest1 = state.finit+stp*state.dgtest;
1283
1284	//
1285	// TEST FOR CONVERGENCE.
1286	//
1287	if( state.brackt & ((double)(stp)<=(double)(state.stmin) \| (double)(stp)>=(double)(state.stmax)) \| state.infoc==0 )
1288	{
1289	info = 6;
1290	}
1291	if( (double)(stp)==(double)(stpmax) & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)<=(double)(state.dgtest) )
1292	{
1293	info = 5;
1294	}
1295	if( (double)(stp)==(double)(stpmin) & ((double)(f)>(double)(state.ftest1) \| (double)(state.dg)>=(double)(state.dgtest)) )
1296	{
1297	info = 4;
1298	}
1299	if( nfev>=maxfev )
1300	{
1301	info = 3;
1302	}
1303	if( state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
1304	{
1305	info = 2;
1306	}
1307	if( (double)(f)<=(double)(state.ftest1) & (double)(Math.Abs(state.dg))<=(double)(-(gtol*state.dginit)) )
1308	{
1309	info = 1;
1310	}
1311
1312	//
1313	// CHECK FOR TERMINATION.
1314	//
1315	if( info!=0 )
1316	{
1317	stage = 0;
1318	return;
1319	}
1320
1321	//
1322	// IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
1323	// FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
1324	//
1325	if( state.stage1 & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)>=(double)(Math.Min(ftol, gtol)*state.dginit) )
1326	{
1327	state.stage1 = false;
1328	}
1329
1330	//
1331	// A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
1332	// WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
1333	// FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
1334	// DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
1335	// OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
1336	//
1337	if( state.stage1 & (double)(f)<=(double)(state.fx) & (double)(f)>(double)(state.ftest1) )
1338	{
1339
1340	//
1341	// DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
1342	//
1343	state.fm = f-stp*state.dgtest;
1344	state.fxm = state.fx-state.stx*state.dgtest;
1345	state.fym = state.fy-state.sty*state.dgtest;
1346	state.dgm = state.dg-state.dgtest;
1347	state.dgxm = state.dgx-state.dgtest;
1348	state.dgym = state.dgy-state.dgtest;
1349
1350	//
1351	// CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
1352	// AND TO COMPUTE THE NEW STEP.
1353	//
1354	mnlmcstep(ref state.stx, ref state.fxm, ref state.dgxm, ref state.sty, ref state.fym, ref state.dgym, ref stp, state.fm, state.dgm, ref state.brackt, state.stmin, state.stmax, ref state.infoc);
1355
1356	//
1357	// RESET THE FUNCTION AND GRADIENT VALUES FOR F.
1358	//
1359	state.fx = state.fxm+state.stx*state.dgtest;
1360	state.fy = state.fym+state.sty*state.dgtest;
1361	state.dgx = state.dgxm+state.dgtest;
1362	state.dgy = state.dgym+state.dgtest;
1363	}
1364	else
1365	{
1366
1367	//
1368	// CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
1369	// AND TO COMPUTE THE NEW STEP.
1370	//
1371	mnlmcstep(ref state.stx, ref state.fx, ref state.dgx, ref state.sty, ref state.fy, ref state.dgy, ref stp, f, state.dg, ref state.brackt, state.stmin, state.stmax, ref state.infoc);
1372	}
1373
1374	//
1375	// FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
1376	// INTERVAL OF UNCERTAINTY.
1377	//
1378	if( state.brackt )
1379	{
1380	if( (double)(Math.Abs(state.sty-state.stx))>=(double)(p66*state.width1) )
1381	{
1382	stp = state.stx+p5*(state.sty-state.stx);
1383	}
1384	state.width1 = state.width;
1385	state.width = Math.Abs(state.sty-state.stx);
1386	}
1387
1388	//
1389	// NEXT.
1390	//
1391	stage = 3;
1392	continue;
1393	}
1394	}
1395	}
1396
1397
1398	private static void mnlmcstep(ref double stx,
1399	ref double fx,
1400	ref double dx,
1401	ref double sty,
1402	ref double fy,
1403	ref double dy,
1404	ref double stp,
1405	double fp,
1406	double dp,
1407	ref bool brackt,
1408	double stmin,
1409	double stmax,
1410	ref int info)
1411	{
1412	bool bound = new bool();
1413	double gamma = 0;
1414	double p = 0;
1415	double q = 0;
1416	double r = 0;
1417	double s = 0;
1418	double sgnd = 0;
1419	double stpc = 0;
1420	double stpf = 0;
1421	double stpq = 0;
1422	double theta = 0;
1423
1424	info = 0;
1425
1426	//
1427	// CHECK THE INPUT PARAMETERS FOR ERRORS.
1428	//
1429	if( brackt & ((double)(stp)<=(double)(Math.Min(stx, sty)) \| (double)(stp)>=(double)(Math.Max(stx, sty))) \| (double)(dx*(stp-stx))>=(double)(0) \| (double)(stmax)<(double)(stmin) )
1430	{
1431	return;
1432	}
1433
1434	//
1435	// DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
1436	//
1437	sgnd = dp*(dx/Math.Abs(dx));
1438
1439	//
1440	// FIRST CASE. A HIGHER FUNCTION VALUE.
1441	// THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
1442	// TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
1443	// ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
1444	//
1445	if( (double)(fp)>(double)(fx) )
1446	{
1447	info = 1;
1448	bound = true;
1449	theta = 3*(fx-fp)/(stp-stx)+dx+dp;
1450	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
1451	gamma = sMath.Sqrt(AP.Math.Sqr(theta/s)-dx/s(dp/s));
1452	if( (double)(stp)<(double)(stx) )
1453	{
1454	gamma = -gamma;
1455	}
1456	p = gamma-dx+theta;
1457	q = gamma-dx+gamma+dp;
1458	r = p/q;
1459	stpc = stx+r*(stp-stx);
1460	stpq = stx+dx/((fx-fp)/(stp-stx)+dx)/2*(stp-stx);
1461	if( (double)(Math.Abs(stpc-stx))<(double)(Math.Abs(stpq-stx)) )
1462	{
1463	stpf = stpc;
1464	}
1465	else
1466	{
1467	stpf = stpc+(stpq-stpc)/2;
1468	}
1469	brackt = true;
1470	}
1471	else
1472	{
1473	if( (double)(sgnd)<(double)(0) )
1474	{
1475
1476	//
1477	// SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
1478	// OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
1479	// STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
1480	// THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
1481	//
1482	info = 2;
1483	bound = false;
1484	theta = 3*(fx-fp)/(stp-stx)+dx+dp;
1485	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
1486	gamma = sMath.Sqrt(AP.Math.Sqr(theta/s)-dx/s(dp/s));
1487	if( (double)(stp)>(double)(stx) )
1488	{
1489	gamma = -gamma;
1490	}
1491	p = gamma-dp+theta;
1492	q = gamma-dp+gamma+dx;
1493	r = p/q;
1494	stpc = stp+r*(stx-stp);
1495	stpq = stp+dp/(dp-dx)*(stx-stp);
1496	if( (double)(Math.Abs(stpc-stp))>(double)(Math.Abs(stpq-stp)) )
1497	{
1498	stpf = stpc;
1499	}
1500	else
1501	{
1502	stpf = stpq;
1503	}
1504	brackt = true;
1505	}
1506	else
1507	{
1508	if( (double)(Math.Abs(dp))<(double)(Math.Abs(dx)) )
1509	{
1510
1511	//
1512	// THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
1513	// SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
1514	// THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
1515	// IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
1516	// IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
1517	// EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
1518	// COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
1519	// CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
1520	//
1521	info = 3;
1522	bound = true;
1523	theta = 3*(fx-fp)/(stp-stx)+dx+dp;
1524	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
1525
1526	//
1527	// THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
1528	// TO INFINITY IN THE DIRECTION OF THE STEP.
1529	//
1530	gamma = sMath.Sqrt(Math.Max(0, AP.Math.Sqr(theta/s)-dx/s(dp/s)));
1531	if( (double)(stp)>(double)(stx) )
1532	{
1533	gamma = -gamma;
1534	}
1535	p = gamma-dp+theta;
1536	q = gamma+(dx-dp)+gamma;
1537	r = p/q;
1538	if( (double)(r)<(double)(0) & (double)(gamma)!=(double)(0) )
1539	{
1540	stpc = stp+r*(stx-stp);
1541	}
1542	else
1543	{
1544	if( (double)(stp)>(double)(stx) )
1545	{
1546	stpc = stmax;
1547	}
1548	else
1549	{
1550	stpc = stmin;
1551	}
1552	}
1553	stpq = stp+dp/(dp-dx)*(stx-stp);
1554	if( brackt )
1555	{
1556	if( (double)(Math.Abs(stp-stpc))<(double)(Math.Abs(stp-stpq)) )
1557	{
1558	stpf = stpc;
1559	}
1560	else
1561	{
1562	stpf = stpq;
1563	}
1564	}
1565	else
1566	{
1567	if( (double)(Math.Abs(stp-stpc))>(double)(Math.Abs(stp-stpq)) )
1568	{
1569	stpf = stpc;
1570	}
1571	else
1572	{
1573	stpf = stpq;
1574	}
1575	}
1576	}
1577	else
1578	{
1579
1580	//
1581	// FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
1582	// SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
1583	// NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
1584	// IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
1585	//
1586	info = 4;
1587	bound = false;
1588	if( brackt )
1589	{
1590	theta = 3*(fp-fy)/(sty-stp)+dy+dp;
1591	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dy), Math.Abs(dp)));
1592	gamma = sMath.Sqrt(AP.Math.Sqr(theta/s)-dy/s(dp/s));
1593	if( (double)(stp)>(double)(sty) )
1594	{
1595	gamma = -gamma;
1596	}
1597	p = gamma-dp+theta;
1598	q = gamma-dp+gamma+dy;
1599	r = p/q;
1600	stpc = stp+r*(sty-stp);
1601	stpf = stpc;
1602	}
1603	else
1604	{
1605	if( (double)(stp)>(double)(stx) )
1606	{
1607	stpf = stmax;
1608	}
1609	else
1610	{
1611	stpf = stmin;
1612	}
1613	}
1614	}
1615	}
1616	}
1617
1618	//
1619	// UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
1620	// DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
1621	//
1622	if( (double)(fp)>(double)(fx) )
1623	{
1624	sty = stp;
1625	fy = fp;
1626	dy = dp;
1627	}
1628	else
1629	{
1630	if( (double)(sgnd)<(double)(0.0) )
1631	{
1632	sty = stx;
1633	fy = fx;
1634	dy = dx;
1635	}
1636	stx = stp;
1637	fx = fp;
1638	dx = dp;
1639	}
1640
1641	//
1642	// COMPUTE THE NEW STEP AND SAFEGUARD IT.
1643	//
1644	stpf = Math.Min(stmax, stpf);
1645	stpf = Math.Max(stmin, stpf);
1646	stp = stpf;
1647	if( brackt & bound )
1648	{
1649	if( (double)(sty)>(double)(stx) )
1650	{
1651	stp = Math.Min(stx+0.66*(sty-stx), stp);
1652	}
1653	else
1654	{
1655	stp = Math.Max(stx+0.66*(sty-stx), stp);
1656	}
1657	}
1658	}
1659	}
1660	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences