ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/search.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (14 years ago) by tjod
Original Path: trunk/smi23d/src/mengine/search.c
File size: 6557 byte(s)
Log Message:
test

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