Context Navigation

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/linmin.cs @ 5192

Visit:

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

Line
1	/*************************************************************************
2	ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
3	JORGE J. MORE', DAVID J. THUENTE
4
5	>>> SOURCE LICENSE >>>
6	This program is free software; you can redistribute it and/or modify
7	it under the terms of the GNU General Public License as published by
8	the Free Software Foundation (www.fsf.org); either version 2 of the
9	License, or (at your option) any later version.
10
11	This program is distributed in the hope that it will be useful,
12	but WITHOUT ANY WARRANTY; without even the implied warranty of
13	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14	GNU General Public License for more details.
15
16	A copy of the GNU General Public License is available at
17	http://www.fsf.org/licensing/licenses
18
19	>>> END OF LICENSE >>>
20	*************************************************************************/
21
22	using System;
23
24	namespace alglib
25	{
26	public class linmin
27	{
28	public struct linminstate
29	{
30	public bool brackt;
31	public bool stage1;
32	public int infoc;
33	public double dg;
34	public double dgm;
35	public double dginit;
36	public double dgtest;
37	public double dgx;
38	public double dgxm;
39	public double dgy;
40	public double dgym;
41	public double finit;
42	public double ftest1;
43	public double fm;
44	public double fx;
45	public double fxm;
46	public double fy;
47	public double fym;
48	public double stx;
49	public double sty;
50	public double stmin;
51	public double stmax;
52	public double width;
53	public double width1;
54	public double xtrapf;
55	};
56
57
58
59
60	public const double ftol = 0.001;
61	public const double xtol = 100*AP.Math.MachineEpsilon;
62	public const double gtol = 0.3;
63	public const int maxfev = 20;
64	public const double stpmin = 1.0E-50;
65	public const double defstpmax = 1.0E+50;
66
67
68	/*************************************************************************
69	Normalizes direction/step pair: makes \|D\|=1, scales Stp.
70	If \|D\|=0, it returns, leavind D/Stp unchanged.
71
72	-- ALGLIB --
73	Copyright 01.04.2010 by Bochkanov Sergey
74	*************************************************************************/
75	public static void linminnormalized(ref double[] d,
76	ref double stp,
77	int n)
78	{
79	double mx = 0;
80	double s = 0;
81	int i = 0;
82	int i_ = 0;
83
84
85	//
86	// first, scale D to avoid underflow/overflow durng squaring
87	//
88	mx = 0;
89	for(i=0; i<=n-1; i++)
90	{
91	mx = Math.Max(mx, Math.Abs(d[i]));
92	}
93	if( (double)(mx)==(double)(0) )
94	{
95	return;
96	}
97	s = 1/mx;
98	for(i_=0; i_<=n-1;i_++)
99	{
100	d[i_] = s*d[i_];
101	}
102	stp = stp/s;
103
104	//
105	// normalize D
106	//
107	s = 0.0;
108	for(i_=0; i_<=n-1;i_++)
109	{
110	s += d[i_]*d[i_];
111	}
112	s = 1/Math.Sqrt(s);
113	for(i_=0; i_<=n-1;i_++)
114	{
115	d[i_] = s*d[i_];
116	}
117	stp = stp/s;
118	}
119
120
121	/*************************************************************************
122	THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES A SUFFICIENT
123	DECREASE CONDITION AND A CURVATURE CONDITION.
124
125	AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF UNCERTAINTY WITH
126	ENDPOINTS STX AND STY. THE INTERVAL OF UNCERTAINTY IS INITIALLY CHOSEN
127	SO THAT IT CONTAINS A MINIMIZER OF THE MODIFIED FUNCTION
128
129	F(X+STPS) - F(X) - FTOLSTP*(GRADF(X)'S).
130
131	IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE
132	FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, THEN THE INTERVAL OF
133	UNCERTAINTY IS CHOSEN SO THAT IT CONTAINS A MINIMIZER OF F(X+STP*S).
134
135	THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES THE SUFFICIENT
136	DECREASE CONDITION
137
138	F(X+STPS) .LE. F(X) + FTOLSTP*(GRADF(X)'S),
139
140	AND THE CURVATURE CONDITION
141
142	ABS(GRADF(X+STPS)'S)) .LE. GTOLABS(GRADF(X)'S).
143
144	IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION IS BOUNDED
145	BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES BOTH CONDITIONS.
146	IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH CONDITIONS, THEN THE
147	ALGORITHM USUALLY STOPS WHEN ROUNDING ERRORS PREVENT FURTHER PROGRESS.
148	IN THIS CASE STP ONLY SATISFIES THE SUFFICIENT DECREASE CONDITION.
149
150	PARAMETERS DESCRIPRION
151
152	N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER OF VARIABLES.
153
154	X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE BASE POINT FOR
155	THE LINE SEARCH. ON OUTPUT IT CONTAINS X+STP*S.
156
157	F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F AT X. ON OUTPUT
158	IT CONTAINS THE VALUE OF F AT X + STP*S.
159
160	G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE GRADIENT OF F AT X.
161	ON OUTPUT IT CONTAINS THE GRADIENT OF F AT X + STP*S.
162
163	S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE SEARCH DIRECTION.
164
165	STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN INITIAL ESTIMATE
166	OF A SATISFACTORY STEP. ON OUTPUT STP CONTAINS THE FINAL ESTIMATE.
167
168	FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION OCCURS WHEN THE
169	SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
170	SATISFIED.
171
172	XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS WHEN THE RELATIVE
173	WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL.
174
175	STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH SPECIFY LOWER AND
176	UPPER BOUNDS FOR THE STEP.
177
178	MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION OCCURS WHEN THE
179	NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV BY THE END OF AN ITERATION.
180
181	INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
182	INFO = 0 IMPROPER INPUT PARAMETERS.
183
184	INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE
185	DIRECTIONAL DERIVATIVE CONDITION HOLD.
186
187	INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
188	IS AT MOST XTOL.
189
190	INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
191
192	INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN.
193
194	INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX.
195
196	INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
197	THERE MAY NOT BE A STEP WHICH SATISFIES THE
198	SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
199	TOLERANCES MAY BE TOO SMALL.
200
201	NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF CALLS TO FCN.
202
203	WA IS A WORK ARRAY OF LENGTH N.
204
205	ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
206	JORGE J. MORE', DAVID J. THUENTE
207	*************************************************************************/
208	public static void mcsrch(int n,
209	ref double[] x,
210	ref double f,
211	ref double[] g,
212	ref double[] s,
213	ref double stp,
214	double stpmax,
215	ref int info,
216	ref int nfev,
217	ref double[] wa,
218	ref linminstate state,
219	ref int stage)
220	{
221	double v = 0;
222	double p5 = 0;
223	double p66 = 0;
224	double zero = 0;
225	int i_ = 0;
226
227
228	//
229	// init
230	//
231	p5 = 0.5;
232	p66 = 0.66;
233	state.xtrapf = 4.0;
234	zero = 0;
235	if( (double)(stpmax)==(double)(0) )
236	{
237	stpmax = defstpmax;
238	}
239	if( (double)(stp)<(double)(stpmin) )
240	{
241	stp = stpmin;
242	}
243	if( (double)(stp)>(double)(stpmax) )
244	{
245	stp = stpmax;
246	}
247
248	//
249	// Main cycle
250	//
251	while( true )
252	{
253	if( stage==0 )
254	{
255
256	//
257	// NEXT
258	//
259	stage = 2;
260	continue;
261	}
262	if( stage==2 )
263	{
264	state.infoc = 1;
265	info = 0;
266
267	//
268	// CHECK THE INPUT PARAMETERS FOR ERRORS.
269	//
270	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 )
271	{
272	stage = 0;
273	return;
274	}
275
276	//
277	// COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
278	// AND CHECK THAT S IS A DESCENT DIRECTION.
279	//
280	v = 0.0;
281	for(i_=0; i_<=n-1;i_++)
282	{
283	v += g[i_]*s[i_];
284	}
285	state.dginit = v;
286	if( (double)(state.dginit)>=(double)(0) )
287	{
288	stage = 0;
289	return;
290	}
291
292	//
293	// INITIALIZE LOCAL VARIABLES.
294	//
295	state.brackt = false;
296	state.stage1 = true;
297	nfev = 0;
298	state.finit = f;
299	state.dgtest = ftol*state.dginit;
300	state.width = stpmax-stpmin;
301	state.width1 = state.width/p5;
302	for(i_=0; i_<=n-1;i_++)
303	{
304	wa[i_] = x[i_];
305	}
306
307	//
308	// THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
309	// FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
310	// THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
311	// FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
312	// THE INTERVAL OF UNCERTAINTY.
313	// THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
314	// FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
315	//
316	state.stx = 0;
317	state.fx = state.finit;
318	state.dgx = state.dginit;
319	state.sty = 0;
320	state.fy = state.finit;
321	state.dgy = state.dginit;
322
323	//
324	// NEXT
325	//
326	stage = 3;
327	continue;
328	}
329	if( stage==3 )
330	{
331
332	//
333	// START OF ITERATION.
334	//
335	// SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
336	// TO THE PRESENT INTERVAL OF UNCERTAINTY.
337	//
338	if( state.brackt )
339	{
340	if( (double)(state.stx)<(double)(state.sty) )
341	{
342	state.stmin = state.stx;
343	state.stmax = state.sty;
344	}
345	else
346	{
347	state.stmin = state.sty;
348	state.stmax = state.stx;
349	}
350	}
351	else
352	{
353	state.stmin = state.stx;
354	state.stmax = stp+state.xtrapf*(stp-state.stx);
355	}
356
357	//
358	// FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
359	//
360	if( (double)(stp)>(double)(stpmax) )
361	{
362	stp = stpmax;
363	}
364	if( (double)(stp)<(double)(stpmin) )
365	{
366	stp = stpmin;
367	}
368
369	//
370	// IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
371	// STP BE THE LOWEST POINT OBTAINED SO FAR.
372	//
373	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) )
374	{
375	stp = state.stx;
376	}
377
378	//
379	// EVALUATE THE FUNCTION AND GRADIENT AT STP
380	// AND COMPUTE THE DIRECTIONAL DERIVATIVE.
381	//
382	for(i_=0; i_<=n-1;i_++)
383	{
384	x[i_] = wa[i_];
385	}
386	for(i_=0; i_<=n-1;i_++)
387	{
388	x[i_] = x[i_] + stp*s[i_];
389	}
390
391	//
392	// NEXT
393	//
394	stage = 4;
395	return;
396	}
397	if( stage==4 )
398	{
399	info = 0;
400	nfev = nfev+1;
401	v = 0.0;
402	for(i_=0; i_<=n-1;i_++)
403	{
404	v += g[i_]*s[i_];
405	}
406	state.dg = v;
407	state.ftest1 = state.finit+stp*state.dgtest;
408
409	//
410	// TEST FOR CONVERGENCE.
411	//
412	if( state.brackt & ((double)(stp)<=(double)(state.stmin) \| (double)(stp)>=(double)(state.stmax)) \| state.infoc==0 )
413	{
414	info = 6;
415	}
416	if( (double)(stp)==(double)(stpmax) & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)<=(double)(state.dgtest) )
417	{
418	info = 5;
419	}
420	if( (double)(stp)==(double)(stpmin) & ((double)(f)>(double)(state.ftest1) \| (double)(state.dg)>=(double)(state.dgtest)) )
421	{
422	info = 4;
423	}
424	if( nfev>=maxfev )
425	{
426	info = 3;
427	}
428	if( state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
429	{
430	info = 2;
431	}
432	if( (double)(f)<=(double)(state.ftest1) & (double)(Math.Abs(state.dg))<=(double)(-(gtol*state.dginit)) )
433	{
434	info = 1;
435	}
436
437	//
438	// CHECK FOR TERMINATION.
439	//
440	if( info!=0 )
441	{
442	stage = 0;
443	return;
444	}
445
446	//
447	// IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
448	// FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
449	//
450	if( state.stage1 & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)>=(double)(Math.Min(ftol, gtol)*state.dginit) )
451	{
452	state.stage1 = false;
453	}
454
455	//
456	// A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
457	// WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
458	// FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
459	// DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
460	// OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
461	//
462	if( state.stage1 & (double)(f)<=(double)(state.fx) & (double)(f)>(double)(state.ftest1) )
463	{
464
465	//
466	// DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
467	//
468	state.fm = f-stp*state.dgtest;
469	state.fxm = state.fx-state.stx*state.dgtest;
470	state.fym = state.fy-state.sty*state.dgtest;
471	state.dgm = state.dg-state.dgtest;
472	state.dgxm = state.dgx-state.dgtest;
473	state.dgym = state.dgy-state.dgtest;
474
475	//
476	// CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
477	// AND TO COMPUTE THE NEW STEP.
478	//
479	mcstep(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);
480
481	//
482	// RESET THE FUNCTION AND GRADIENT VALUES FOR F.
483	//
484	state.fx = state.fxm+state.stx*state.dgtest;
485	state.fy = state.fym+state.sty*state.dgtest;
486	state.dgx = state.dgxm+state.dgtest;
487	state.dgy = state.dgym+state.dgtest;
488	}
489	else
490	{
491
492	//
493	// CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
494	// AND TO COMPUTE THE NEW STEP.
495	//
496	mcstep(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);
497	}
498
499	//
500	// FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
501	// INTERVAL OF UNCERTAINTY.
502	//
503	if( state.brackt )
504	{
505	if( (double)(Math.Abs(state.sty-state.stx))>=(double)(p66*state.width1) )
506	{
507	stp = state.stx+p5*(state.sty-state.stx);
508	}
509	state.width1 = state.width;
510	state.width = Math.Abs(state.sty-state.stx);
511	}
512
513	//
514	// NEXT.
515	//
516	stage = 3;
517	continue;
518	}
519	}
520	}
521
522
523	private static void mcstep(ref double stx,
524	ref double fx,
525	ref double dx,
526	ref double sty,
527	ref double fy,
528	ref double dy,
529	ref double stp,
530	double fp,
531	double dp,
532	ref bool brackt,
533	double stmin,
534	double stmax,
535	ref int info)
536	{
537	bool bound = new bool();
538	double gamma = 0;
539	double p = 0;
540	double q = 0;
541	double r = 0;
542	double s = 0;
543	double sgnd = 0;
544	double stpc = 0;
545	double stpf = 0;
546	double stpq = 0;
547	double theta = 0;
548
549	info = 0;
550
551	//
552	// CHECK THE INPUT PARAMETERS FOR ERRORS.
553	//
554	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) )
555	{
556	return;
557	}
558
559	//
560	// DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
561	//
562	sgnd = dp*(dx/Math.Abs(dx));
563
564	//
565	// FIRST CASE. A HIGHER FUNCTION VALUE.
566	// THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
567	// TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
568	// ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
569	//
570	if( (double)(fp)>(double)(fx) )
571	{
572	info = 1;
573	bound = true;
574	theta = 3*(fx-fp)/(stp-stx)+dx+dp;
575	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
576	gamma = sMath.Sqrt(AP.Math.Sqr(theta/s)-dx/s(dp/s));
577	if( (double)(stp)<(double)(stx) )
578	{
579	gamma = -gamma;
580	}
581	p = gamma-dx+theta;
582	q = gamma-dx+gamma+dp;
583	r = p/q;
584	stpc = stx+r*(stp-stx);
585	stpq = stx+dx/((fx-fp)/(stp-stx)+dx)/2*(stp-stx);
586	if( (double)(Math.Abs(stpc-stx))<(double)(Math.Abs(stpq-stx)) )
587	{
588	stpf = stpc;
589	}
590	else
591	{
592	stpf = stpc+(stpq-stpc)/2;
593	}
594	brackt = true;
595	}
596	else
597	{
598	if( (double)(sgnd)<(double)(0) )
599	{
600
601	//
602	// SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
603	// OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
604	// STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
605	// THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
606	//
607	info = 2;
608	bound = false;
609	theta = 3*(fx-fp)/(stp-stx)+dx+dp;
610	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
611	gamma = sMath.Sqrt(AP.Math.Sqr(theta/s)-dx/s(dp/s));
612	if( (double)(stp)>(double)(stx) )
613	{
614	gamma = -gamma;
615	}
616	p = gamma-dp+theta;
617	q = gamma-dp+gamma+dx;
618	r = p/q;
619	stpc = stp+r*(stx-stp);
620	stpq = stp+dp/(dp-dx)*(stx-stp);
621	if( (double)(Math.Abs(stpc-stp))>(double)(Math.Abs(stpq-stp)) )
622	{
623	stpf = stpc;
624	}
625	else
626	{
627	stpf = stpq;
628	}
629	brackt = true;
630	}
631	else
632	{
633	if( (double)(Math.Abs(dp))<(double)(Math.Abs(dx)) )
634	{
635
636	//
637	// THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
638	// SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
639	// THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
640	// IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
641	// IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
642	// EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
643	// COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
644	// CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
645	//
646	info = 3;
647	bound = true;
648	theta = 3*(fx-fp)/(stp-stx)+dx+dp;
649	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
650
651	//
652	// THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
653	// TO INFINITY IN THE DIRECTION OF THE STEP.
654	//
655	gamma = sMath.Sqrt(Math.Max(0, AP.Math.Sqr(theta/s)-dx/s(dp/s)));
656	if( (double)(stp)>(double)(stx) )
657	{
658	gamma = -gamma;
659	}
660	p = gamma-dp+theta;
661	q = gamma+(dx-dp)+gamma;
662	r = p/q;
663	if( (double)(r)<(double)(0) & (double)(gamma)!=(double)(0) )
664	{
665	stpc = stp+r*(stx-stp);
666	}
667	else
668	{
669	if( (double)(stp)>(double)(stx) )
670	{
671	stpc = stmax;
672	}
673	else
674	{
675	stpc = stmin;
676	}
677	}
678	stpq = stp+dp/(dp-dx)*(stx-stp);
679	if( brackt )
680	{
681	if( (double)(Math.Abs(stp-stpc))<(double)(Math.Abs(stp-stpq)) )
682	{
683	stpf = stpc;
684	}
685	else
686	{
687	stpf = stpq;
688	}
689	}
690	else
691	{
692	if( (double)(Math.Abs(stp-stpc))>(double)(Math.Abs(stp-stpq)) )
693	{
694	stpf = stpc;
695	}
696	else
697	{
698	stpf = stpq;
699	}
700	}
701	}
702	else
703	{
704
705	//
706	// FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
707	// SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
708	// NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
709	// IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
710	//
711	info = 4;
712	bound = false;
713	if( brackt )
714	{
715	theta = 3*(fp-fy)/(sty-stp)+dy+dp;
716	s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dy), Math.Abs(dp)));
717	gamma = sMath.Sqrt(AP.Math.Sqr(theta/s)-dy/s(dp/s));
718	if( (double)(stp)>(double)(sty) )
719	{
720	gamma = -gamma;
721	}
722	p = gamma-dp+theta;
723	q = gamma-dp+gamma+dy;
724	r = p/q;
725	stpc = stp+r*(sty-stp);
726	stpf = stpc;
727	}
728	else
729	{
730	if( (double)(stp)>(double)(stx) )
731	{
732	stpf = stmax;
733	}
734	else
735	{
736	stpf = stmin;
737	}
738	}
739	}
740	}
741	}
742
743	//
744	// UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
745	// DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
746	//
747	if( (double)(fp)>(double)(fx) )
748	{
749	sty = stp;
750	fy = fp;
751	dy = dp;
752	}
753	else
754	{
755	if( (double)(sgnd)<(double)(0.0) )
756	{
757	sty = stx;
758	fy = fx;
759	dy = dx;
760	}
761	stx = stp;
762	fx = fp;
763	dx = dp;
764	}
765
766	//
767	// COMPUTE THE NEW STEP AND SAFEGUARD IT.
768	//
769	stpf = Math.Min(stmax, stpf);
770	stpf = Math.Max(stmin, stpf);
771	stp = stpf;
772	if( brackt & bound )
773	{
774	if( (double)(sty)>(double)(stx) )
775	{
776	stp = Math.Min(stx+0.66*(sty-stx), stp);
777	}
778	else
779	{
780	stp = Math.Max(stx+0.66*(sty-stx), stp);
781	}
782	}
783	}
784	}
785	}

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

Download in other formats:

Update cookies preferences