ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/search.c
Revision: 85
Committed: Mon Jan 12 15:05:44 2009 UTC (13 years, 5 months ago) by gilbertke
Original Path: trunk/src/mengine/src/search.c
File size: 6547 byte(s)
Log Message:
more cleanup and code removal
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     #include "energies.h"
7     #include "utility.h"
8    
9     EXTERN struct t_minvar{
10     double cappa, stpmin, stpmax, angmax;
11     int intmax;
12     } minvar;
13    
14     #define Success 0
15     #define ReSearch 1
16     #define WideAngle 2
17     #define BadIntpln 3
18     #define IntplnErr 4
19     #define blank 5
20     #define Failure 6
21    
22    
23     //#define min(a,b) (((a) < (b) ? (a) : (b))
24    
25     double minimize1(double *, double *);
26     void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *);
27    
28     void search(int nvar, double *f, double *g, double *x, double *p,
29     double f_move,double *angle,int *ncalls,double (*fgvalue)(),int *status)
30     {
31     int i, intpln;
32     int maxvar;
33     double *s; // s[maxvar];
34     double s_norm,g_norm,cosang;
35     double step,parab,cube,cubstp,sss,ttt;
36     double f_0,f_1,f_a,f_b,f_c,sg_0,sg_1,sg_a,sg_b,sg_c;
37    
38     maxvar = 3*natom;
39     s = dvector(0,maxvar);
40    
41     if (minvar.cappa == 0.0) minvar.cappa = 0.1;
42     if (minvar.stpmin == 0.0) minvar.stpmin = 1.0e-20;
43     if (minvar.stpmax == 0.0) minvar.stpmax = 2.0;
44     if (minvar.angmax == 0.0) minvar.angmax = 180.0;
45     if (minvar.intmax == 0) minvar.intmax = 5;
46    
47     for (i=0; i < nvar; i++)
48     s[i] = p[i];
49    
50     g_norm = 0.0;
51     s_norm = 0.0;
52     for (i=0; i < nvar; i++)
53     {
54     g_norm += g[i]*g[i];
55     s_norm += s[i]*s[i];
56     }
57     g_norm = sqrt(g_norm);
58     s_norm = sqrt(s_norm);
59     if (isnan(g_norm) != 0)
60     {
61 wdelano 58 fprintf(pcmlogfile,"Found nan in search\n");
62 tjod 3 *status = Failure;
63     free_dvector(s,0,maxvar);
64     return;
65     }
66     f_0 = *f;
67     sg_0 = 0.0;
68     for (i=0; i < nvar; i++)
69     {
70     s[i] = s[i] / s_norm;
71     sg_0 += s[i]*g[i];
72     }
73     cosang = -sg_0 / g_norm;
74     if (cosang < -1.00)
75     cosang = -1.0;
76     if (cosang > 1.00)
77     cosang = 1.0;
78     *angle = radian * acos(cosang);
79     if (*angle > minvar.angmax)
80     {
81     *status = WideAngle;
82     free_dvector(s,0,maxvar);
83     return;
84     }
85     step = 2.0 * fabs(f_move/sg_0);
86     step = MinFun(step,s_norm);
87     if (step > minvar.stpmax) step = minvar.stpmax;
88     if (step < minvar.stpmin) step = minvar.stpmin;
89    
90     L_10:
91     intpln = 0;
92     f_b = f_0;
93     sg_b = sg_0;
94    
95     L_20:
96     f_a = f_b;
97     sg_a = sg_b;
98     for (i=0; i < nvar; i++)
99     x[i] = x[i] + step*s[i];
100    
101     // get new function and gradient
102    
103     ncalls++;
104     f_b = fgvalue(x,g);
105     sg_b = 0.0;
106    
107     for (i=0; i < nvar; i++)
108     sg_b += s[i]*g[i];
109    
110     if (fabs(sg_b/sg_0) <= minvar.cappa && f_b < f_a)
111     {
112     *f = f_b;
113     *status = Success;
114     free_dvector(s,0,maxvar);
115     return;
116     }
117    
118     if (sg_a*sg_b < 0.0 || f_a < f_b) goto L_30;
119     if ( fabs(sg_a-sg_b) < 0.0000001)
120     {
121     *f = f_b;
122     *status = Success;
123     free_dvector(s,0,maxvar);
124     return;
125     }
126     step = 2.0 * step;
127     if (sg_b > sg_a)
128     {
129     parab = (f_a-f_b) / (sg_b-sg_a);
130     if (parab > 2.0*step) parab = 2.0 * step;
131     if (parab < 0.5*step) parab = 0.5 * step;
132     step = parab;
133     }
134     if (step > minvar.stpmax) step = minvar.stpmax;
135     goto L_20;
136    
137     L_30:
138     intpln++;
139     sss = 3.0*(f_b-f_a)/step - sg_a - sg_b;
140     ttt = sss*sss - sg_a*sg_b;
141     if (ttt < 0.0)
142     {
143     *f = f_b;
144     *status = IntplnErr;
145     free_dvector(s,0,maxvar);
146     return;
147     }
148     ttt = sqrt(ttt);
149     cube = step * (sg_b+ttt+sss)/(sg_b-sg_a+2.0*ttt);
150     if (cube < 0.0 || cube > step)
151     {
152     *f = f_b;
153     *status = IntplnErr;
154     free_dvector(s,0,maxvar);
155     return;
156     }
157     for (i=0; i < nvar; i++)
158     x[i] -= cube*s[i];
159     // get new function and gradient and test
160     ncalls++;
161     f_c = fgvalue(x,g);
162     sg_c = 0.0;
163     for (i=0; i < nvar; i++)
164     sg_c += s[i]*g[i];
165    
166     if (fabs(sg_c/sg_0) <= minvar.cappa)
167     {
168     *f = f_c;
169     *status = Success;
170     free_dvector(s,0,maxvar);
171     return;
172     }
173     if (f_c <= f_a || f_c <= f_b)
174     {
175     cubstp = MinFun(fabs(cube),fabs(step-cube));
176     if (cubstp >= minvar.stpmin && intpln < minvar.intmax)
177     {
178     if (sg_a*sg_b < 0.0)
179     {
180     if (sg_a*sg_c < 0.0)
181     {
182     f_b = f_c;
183     sg_b = sg_c;
184     step = step - cube;
185     }else
186     {
187     f_a = f_c;
188     sg_a = sg_c;
189     step = cube;
190     for(i=0; i < nvar; i++)
191     x[i] = x[i] + cube*s[i];
192    
193     }
194     } else
195     {
196     if (sg_a*sg_c < 0.0 || f_a <= f_c)
197     {
198     f_b = f_c;
199     sg_b = sg_c;
200     step = step - cube;
201     } else
202     {
203     f_a = f_c;
204     sg_a = sg_c;
205     step = cube;
206     for(i=0; i < nvar; i++)
207     x[i] = x[i] + cube*s[i];
208    
209     }
210     }
211     goto L_30;
212     }
213     }
214    
215     // failed interpolation
216     if (f_a < f_b && f_a < f_c)
217     {
218     f_1 = f_a;
219     sg_1 = sg_a;
220     for (i=0; i < nvar; i++)
221     x[i] += (cube-step)*s[i];
222     } else if (f_b < f_a && f_b < f_c)
223     {
224     f_1 = f_b;
225     sg_1 = sg_b;
226     for (i=0; i < nvar; i++)
227     x[i] += cube*s[i];
228     } else
229     {
230     f_1 = f_c;
231     sg_1 = sg_c;
232     }
233     // try to restart
234     if (f_1 > f_0)
235     {
236     ncalls++;
237     *f = fgvalue(x,g);
238     *status = IntplnErr;
239     free_dvector(s,0,maxvar);
240     return;
241     }
242     f_0 = f_1;
243     sg_0 = sg_1;
244     if (sg_1 > 0.0)
245     {
246     for (i=0; i < nvar; i++)
247     s[i] = -s[i];
248     sg_0 = -sg_1;
249     }
250     if (cube > (step-cube) )
251     step = cube/10.0;
252     else
253     step = (step-cube)/10.0;
254    
255     if (step < minvar.stpmin) step = minvar.stpmin;
256    
257     if (*status == ReSearch)
258     {
259     ncalls++;
260     *f = fgvalue(x,g);
261     *status = BadIntpln;
262     free_dvector(s,0,maxvar);
263     return;
264     } else
265     {
266     *status = ReSearch;
267     goto L_10;
268     }
269     }