~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathNum_simplexmin.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA 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  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_simplexmin.c,v $
 27  * Date    :  $Date: 2009/03/19 18:36:26 $
 28  * Version :  $Revision: 1.11 $
 29  * CVS Id  :  $Id: mathNum_simplexmin.c,v 1.11 2009/03/19 18:36:26 paul Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes :
 34  *
 35  *********
 36 */
 37 /** 
 38  *  @file
 39  *  @brief  Simplex Functions.  
 40  *
 41  *  Method based upon the algorithm published by Nelder and Mead.
 42  *  See Numerical Recipes (3rd Ed) for descriptions of the equivalent functions;
 43  *
 44  *  - amoeba() and amotry() Ch 10.4 pp 408-412.
 45  *  - Simplex (general) Ch 10.8 pp 430-443.
 46 */
 47 
 48 #include "mathNum_simplexmin.h"
 49 
 50 #if HAVE_CONFIG_H
 51 #include <config.h>
 52 #endif
 53 
 54 #include <math.h>
 55 #include <stdio.h>
 56 #include <tina/sys/sysDef.h>
 57 #include <tina/sys/sysPro.h>
 58 #include <tina/math/mathNum_simplexmin.h>
 59 
 60 #define NMAX 1000
 61 #define ALPHA 1.0
 62 #define BETA 0.5
 63 #define GAMMA 2.0
 64 
 65 /* smc 2/2/93 */
 66 #ifdef Print
 67 #undef Print
 68 #endif
 69 /* smc 2/2/93 */
 70 
 71 #define Print(x)  { if(out_text!=NULL) out_text(x); }
 72 
 73 static void     (*out_text) () = NULL;
 74 
 75 double          simplexmin(int n, double *x, double lambda, double (*funk) ( /* ??? */ ), void *data, double ftol, void (*text) ( /* ??? */ ))
 76 {
 77         int             i, j, lowest = 0;
 78         double        **p;
 79         double         *y;
 80         int             nfunk;
 81         double          f;
 82         void            amoeba(double **p, double *y, int ndim, double ftol, double (*funk) ( /* ??? */ ), void*data, int *nfunk);
 83         char            temp[256];
 84 
 85         out_text = text;
 86         y = (double *) ralloc((unsigned) ((n + 1) * sizeof(double)));
 87         if ((p = darray_alloc(0, 0, n + 1, n)) == NULL)
 88                 return (0.0);
 89 
 90         for (i = 0; i < n + 1; i++)
 91         {
 92                 for (j = 0; j < n; j++)
 93                 {
 94                         p[i][j] = x[j];
 95                 }
 96                 if (i > 0)
 97                         p[i][i - 1] = x[i - 1] * (1.0 + lambda);
 98                 y[i] = (*funk) (n, p[i], data);
 99                 if (i > 0)
100                 {
101                         p[i][i - 1] = x[i - 1] * (1.0 - lambda);
102                         if ((f = (*funk) (n, p[i], data)) > y[i])
103                                 p[i][i - 1] = x[i - 1] * (1.0 + lambda);
104                         else
105                                 y[i] = f;
106                 }
107         }
108 
109         sprintf(temp, " initial value %10.4e \n", y[0]);
110         Print(temp);
111         amoeba(p, y, n, ftol, funk, data, &nfunk);
112         for (j = 0; j < n + 1; j++)
113                 if (y[j] < y[lowest])
114                         lowest = j;
115         for (i = 0; i < n; i++)
116                 x[i] = p[lowest][i];
117         f = y[lowest];
118         darray_free(p, 0, 0, n + 1, n);
119         rfree((void *) y);
120         sprintf(temp, " number of function evaluations %d \n", nfunk);
121         Print(temp);
122         sprintf(temp, " final value %10.4e \n", f);
123         Print(temp);
124         return (f);
125 }
126 
127 double          simplexmin2(int n, double *x, double *lambda, double (*funk) ( /* ??? */ ), void *data, double ftol, void (*text) ( /* ??? */ ))
128 {
129         int             i, j, lowest = 0;
130         double        **p;
131         double         *y;
132         int             nfunk;
133         double          f;
134         void            amoeba(double **p, double *y, int ndim, double ftol, double (*funk) ( /* ??? */ ), void*data, int *nfunk);
135         char            temp[256];
136 
137         out_text = text;
138         y = (double *) ralloc((unsigned) ((n + 1) * sizeof(double)));
139         if ((p = darray_alloc(0, 0, n + 1, n)) == NULL)
140                 return (0.0);
141 
142         for (i = 0; i < n + 1; i++)
143         {
144                 for (j = 0; j < n; j++)
145                 {
146                         p[i][j] = x[j];
147                 }
148                 if (i > 0)
149                         p[i][i - 1] = x[i - 1] + lambda[i - 1];
150                 y[i] = (*funk) (n, p[i], data);
151                 if (i > 0)
152                 {
153                         p[i][i - 1] = x[i - 1] - lambda[i - 1];
154                         if ((f = (*funk) (n, p[i], data)) > y[i])
155                                 p[i][i - 1] = x[i - 1] + lambda[i - 1];
156                         else
157                                 y[i] = f;
158                 }
159         }
160 
161         sprintf(temp, " initial value %10.4e \n", y[0]);
162         Print(temp);
163         amoeba(p, y, n, ftol, funk, data, &nfunk);
164         for (j = 0; j < n + 1; j++)
165                 if (y[j] < y[lowest])
166                         lowest = j;
167         for (i = 0; i < n; i++)
168                 x[i] = p[lowest][i];
169         f = y[lowest];
170         darray_free(p, 0, 0, n + 1, n);
171         rfree((void *) y);
172         sprintf(temp, " number of function evaluations %d \n", nfunk);
173         Print(temp);
174         sprintf(temp, " final value %10.4e \n", f);
175         Print(temp);
176         return (f);
177 }
178 
179 void            amoeba(double **p, double *y, int ndim, double ftol, double (*funk) ( /* ??? */ ), void *data, int *nfunk)
180 {
181         int             i, j, ilo, ihi, inhi, mpts = ndim + 1;
182         double          ytry, ysave, sum, rtol, amotry(double **p, double *y, double *psum, int ndim, double (*funk) ( /* ??? */ ), void *data, int ihi, int *nfunk, double fac), *psum;
183 
184         psum = (double *) ralloc((unsigned) ndim * sizeof(double));
185         *nfunk = 0;
186         for (j = 0; j < ndim; j++)
187         {
188                 for (i = 0, sum = 0.0; i < mpts; i++)
189                         sum += p[i][j];
190                 psum[j] = sum;
191         }
192         for (;;)
193         {
194                 ilo = 0;
195                 ihi = y[0] > y[1] ? (inhi = 1, 0) : (inhi = 0, 1);
196                 for (i = 0; i < mpts; i++)
197                 {
198                         if (y[i] < y[ilo])
199                                 ilo = i;
200                         if (y[i] > y[ihi])
201                         {
202                                 inhi = ihi;
203                                 ihi = i;
204                         } else if (y[i] > y[inhi])
205                                 if (i != ihi)
206                                         inhi = i;
207                 }
208                 rtol = 2.0 * fabs(y[ihi] - y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]));
209                 if ((rtol > ftol && *nfunk <= NMAX)||(*nfunk < 25));
210                 else
211                         break;
212 
213                 ytry = amotry(p, y, psum, ndim, funk, data, ihi, nfunk, -ALPHA);
214 
215                 if (ytry <= y[ilo])
216                         ytry = amotry(p, y, psum, ndim, funk, data, ihi, nfunk, GAMMA);
217                 else if (ytry >= y[inhi])
218                 {
219                         ysave = y[ihi];
220                         ytry = amotry(p, y, psum, ndim, funk, data, ihi, nfunk, BETA);
221                         if (ytry >= ysave)
222                         {
223                                 for (i = 0; i < mpts; i++)
224                                 {
225                                         if (i != ilo)
226                                         {
227                                                 for (j = 0; j < ndim; j++)
228                                                 {
229                                                         psum[j] = 0.5 * (p[i][j] + p[ilo][j]);
230                                                         p[i][j] = psum[j];
231                                                 }
232                                                 y[i] = (*funk) (ndim, psum, data);
233                                         }
234                                 }
235                                 *nfunk += ndim;
236                                 for (j = 0; j < ndim; j++)
237                                 {
238                                         for (i = 0, sum = 0.0; i < mpts; i++)
239                                                 sum += p[i][j];
240                                         psum[j] = sum;
241                                 }
242                         }
243                 }
244         }
245         rfree((void *) psum);
246 }
247 
248 double          amotry(double **p, double *y, double *psum, int ndim, double (*funk) ( /* ??? */ ), void *data, int ihi, int *nfunk, double fac)
249 {
250         int             j;
251         double          fac1, fac2, ytry, *ptry;
252 
253         ptry = (double *) ralloc((unsigned) ndim * sizeof(double));;
254         fac1 = (1.0 - fac) / ndim;
255         fac2 = fac1 - fac;
256         for (j = 0; j < ndim; j++)
257                 ptry[j] = psum[j] * fac1 - p[ihi][j] * fac2;
258         ytry = (*funk) (ndim, ptry, data);
259         ++(*nfunk);
260         if (ytry < y[ihi])
261         {
262                 y[ihi] = ytry;
263                 for (j = 0; j < ndim; j++)
264                 {
265                         psum[j] += ptry[j] - p[ihi][j];
266                         p[ihi][j] = ptry[j];
267                 }
268         }
269         rfree((void *) ptry);
270         return ytry;
271 }
272 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.