"Fossies" - the Fresh Open Source Software Archive

Member "darktable-3.6.1/src/iop/ashift_nmsimplex.c" (10 Sep 2021, 10997 Bytes) of package /linux/misc/darktable-3.6.1.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "ashift_nmsimplex.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 3.4.1.1_vs_3.6.0.

    1 /*
    2   This file is part of darktable,
    3   Copyright (C) 2016-2020 darktable developers.
    4 
    5   darktable 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, either version 3 of the License, or
    8   (at your option) any later version.
    9 
   10   darktable 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 darktable.  If not, see <http://www.gnu.org/licenses/>.
   17 */
   18 
   19 /* For parameter optimization we are using the Nelder-Mead simplex method
   20  * implemented by Michael F. Hutt.
   21  * Changes versus the original code:
   22  *      do not include "nmsimplex.h" (not needed)
   23  *      renamed configuration variables to NMS_*
   24  *      add additional argument to objfun for arbitrary parameters
   25  *      simplex() returns number of used iterations instead of min value
   26  *      maximum number of iterations as function parameter
   27  *      make interface function simplex() static
   28  *      initialize i and j to avoid compiler warnings
   29  *      comment out printing of status inormation
   30  *      reformat according to darktable's clang standards
   31  */
   32 
   33 /*==================================================================================
   34  * begin nmsimplex code downloaded from http://www.mikehutt.com/neldermead.html
   35  * on February 6, 2016
   36  *==================================================================================*/
   37 /*
   38  * Program: nmsimplex.c
   39  * Author : Michael F. Hutt
   40  * http://www.mikehutt.com
   41  * 11/3/97
   42  *
   43  * An implementation of the Nelder-Mead simplex method.
   44  *
   45  * Copyright (c) 1997-2011 <Michael F. Hutt>
   46  *
   47  * Permission is hereby granted, free of charge, to any person obtaining
   48  * a copy of this software and associated documentation files (the
   49  * "Software"), to deal in the Software without restriction, including
   50  * without limitation the rights to use, copy, modify, merge, publish,
   51  * distribute, sublicense, and/or sell copies of the Software, and to
   52  * permit persons to whom the Software is furnished to do so, subject to
   53  * the following conditions:
   54  *
   55  * The above copyright notice and this permission notice shall be
   56  * included in all copies or substantial portions of the Software.
   57  *
   58  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   59  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   60  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   61  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
   62  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
   63  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
   64  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
   65  *
   66  *
   67  * Jan. 6, 1999
   68  * Modified to conform to the algorithm presented
   69  * in Margaret H. Wright's paper on Direct Search Methods.
   70  *
   71  * Jul. 23, 2007
   72  * Fixed memory leak.
   73  *
   74  * Mar. 1, 2011
   75  * Added constraints.
   76  */
   77 
   78 //#include "nmsimplex.h"
   79 
   80 static int simplex(double (*objfunc)(double[], void *params), double start[], int n, double EPSILON, double scale,
   81                    int maxiter, void (*constrain)(double[], int n), void *params)
   82 {
   83 
   84   int vs; /* vertex with smallest value */
   85   int vh; /* vertex with next smallest value */
   86   int vg; /* vertex with largest value */
   87 
   88   int i = 0, j = 0, m, row;
   89   int k;   /* track the number of function evaluations */
   90   int itr; /* track the number of iterations */
   91 
   92   double **v;    /* holds vertices of simplex */
   93   double pn, qn; /* values used to create initial simplex */
   94   double *f;     /* value of function at each vertex */
   95   double fr;     /* value of function at reflection point */
   96   double fe;     /* value of function at expansion point */
   97   double fc;     /* value of function at contraction point */
   98   double *vr;    /* reflection - coordinates */
   99   double *ve;    /* expansion - coordinates */
  100   double *vc;    /* contraction - coordinates */
  101   double *vm;    /* centroid - coordinates */
  102   //double min;
  103 
  104   double fsum, favg, s, cent;
  105 
  106   /* dynamically allocate arrays */
  107 
  108   /* allocate the rows of the arrays */
  109   v = (double **)malloc(sizeof(double *) * (n + 1));
  110   f = (double *)malloc(sizeof(double) * (n + 1));
  111   vr = (double *)malloc(sizeof(double) * n);
  112   ve = (double *)malloc(sizeof(double) * n);
  113   vc = (double *)malloc(sizeof(double) * n);
  114   vm = (double *)malloc(sizeof(double) * n);
  115 
  116   /* allocate the columns of the arrays */
  117   for(i = 0; i <= n; i++)
  118   {
  119     v[i] = (double *)malloc(sizeof(double) * n);
  120   }
  121 
  122   /* create the initial simplex */
  123   /* assume one of the vertices is 0,0 */
  124 
  125   pn = scale * (sqrt(n + 1) - 1 + n) / (n * sqrt(2));
  126   qn = scale * (sqrt(n + 1) - 1) / (n * sqrt(2));
  127 
  128   for(i = 0; i < n; i++)
  129   {
  130     v[0][i] = start[i];
  131   }
  132 
  133   for(i = 1; i <= n; i++)
  134   {
  135     for(j = 0; j < n; j++)
  136     {
  137       if(i - 1 == j)
  138       {
  139         v[i][j] = pn + start[j];
  140       }
  141       else
  142       {
  143         v[i][j] = qn + start[j];
  144       }
  145     }
  146   }
  147 
  148   if(constrain != NULL)
  149   {
  150     constrain(v[j], n);
  151   }
  152   /* find the initial function values */
  153   for(j = 0; j <= n; j++)
  154   {
  155     f[j] = objfunc(v[j], params);
  156   }
  157 
  158   k = n + 1;
  159 #if 0
  160   /* print out the initial values */
  161   printf("Initial Values\n");
  162   for(j = 0; j <= n; j++)
  163   {
  164     for(i = 0; i < n; i++)
  165     {
  166       printf("%f %f\n", v[j][i], f[j]);
  167     }
  168   }
  169 #endif
  170 
  171   /* begin the main loop of the minimization */
  172   for(itr = 1; itr <= maxiter; itr++)
  173   {
  174     /* find the index of the largest value */
  175     vg = 0;
  176     for(j = 0; j <= n; j++)
  177     {
  178       if(f[j] > f[vg])
  179       {
  180         vg = j;
  181       }
  182     }
  183 
  184     /* find the index of the smallest value */
  185     vs = 0;
  186     for(j = 0; j <= n; j++)
  187     {
  188       if(f[j] < f[vs])
  189       {
  190         vs = j;
  191       }
  192     }
  193 
  194     /* find the index of the second largest value */
  195     vh = vs;
  196     for(j = 0; j <= n; j++)
  197     {
  198       if(f[j] > f[vh] && f[j] < f[vg])
  199       {
  200         vh = j;
  201       }
  202     }
  203 
  204     /* calculate the centroid */
  205     for(j = 0; j <= n - 1; j++)
  206     {
  207       cent = 0.0;
  208       for(m = 0; m <= n; m++)
  209       {
  210         if(m != vg)
  211         {
  212           cent += v[m][j];
  213         }
  214       }
  215       vm[j] = cent / n;
  216     }
  217 
  218     /* reflect vg to new vertex vr */
  219     for(j = 0; j <= n - 1; j++)
  220     {
  221       /*vr[j] = (1+NMS_ALPHA)*vm[j] - NMS_ALPHA*v[vg][j];*/
  222       vr[j] = vm[j] + NMS_ALPHA * (vm[j] - v[vg][j]);
  223     }
  224     if(constrain != NULL)
  225     {
  226       constrain(vr, n);
  227     }
  228     fr = objfunc(vr, params);
  229     k++;
  230 
  231     if(fr < f[vh] && fr >= f[vs])
  232     {
  233       for(j = 0; j <= n - 1; j++)
  234       {
  235         v[vg][j] = vr[j];
  236       }
  237       f[vg] = fr;
  238     }
  239 
  240     /* investigate a step further in this direction */
  241     if(fr < f[vs])
  242     {
  243       for(j = 0; j <= n - 1; j++)
  244       {
  245         /*ve[j] = NMS_GAMMA*vr[j] + (1-NMS_GAMMA)*vm[j];*/
  246         ve[j] = vm[j] + NMS_GAMMA * (vr[j] - vm[j]);
  247       }
  248       if(constrain != NULL)
  249       {
  250         constrain(ve, n);
  251       }
  252       fe = objfunc(ve, params);
  253       k++;
  254 
  255       /* by making fe < fr as opposed to fe < f[vs],
  256          Rosenbrocks function takes 63 iterations as opposed
  257          to 64 when using double variables. */
  258 
  259       if(fe < fr)
  260       {
  261         for(j = 0; j <= n - 1; j++)
  262         {
  263           v[vg][j] = ve[j];
  264         }
  265         f[vg] = fe;
  266       }
  267       else
  268       {
  269         for(j = 0; j <= n - 1; j++)
  270         {
  271           v[vg][j] = vr[j];
  272         }
  273         f[vg] = fr;
  274       }
  275     }
  276 
  277     /* check to see if a contraction is necessary */
  278     if(fr >= f[vh])
  279     {
  280       if(fr < f[vg] && fr >= f[vh])
  281       {
  282         /* perform outside contraction */
  283         for(j = 0; j <= n - 1; j++)
  284         {
  285           /*vc[j] = NMS_BETA*v[vg][j] + (1-NMS_BETA)*vm[j];*/
  286           vc[j] = vm[j] + NMS_BETA * (vr[j] - vm[j]);
  287         }
  288         if(constrain != NULL)
  289         {
  290           constrain(vc, n);
  291         }
  292         fc = objfunc(vc, params);
  293         k++;
  294       }
  295       else
  296       {
  297         /* perform inside contraction */
  298         for(j = 0; j <= n - 1; j++)
  299         {
  300           /*vc[j] = NMS_BETA*v[vg][j] + (1-NMS_BETA)*vm[j];*/
  301           vc[j] = vm[j] - NMS_BETA * (vm[j] - v[vg][j]);
  302         }
  303         if(constrain != NULL)
  304         {
  305           constrain(vc, n);
  306         }
  307         fc = objfunc(vc, params);
  308         k++;
  309       }
  310 
  311 
  312       if(fc < f[vg])
  313       {
  314         for(j = 0; j <= n - 1; j++)
  315         {
  316           v[vg][j] = vc[j];
  317         }
  318         f[vg] = fc;
  319       }
  320       /* at this point the contraction is not successful,
  321          we must halve the distance from vs to all the
  322          vertices of the simplex and then continue.
  323          10/31/97 - modified to account for ALL vertices.
  324       */
  325       else
  326       {
  327         for(row = 0; row <= n; row++)
  328         {
  329           if(row != vs)
  330           {
  331             for(j = 0; j <= n - 1; j++)
  332             {
  333               v[row][j] = v[vs][j] + (v[row][j] - v[vs][j]) / 2.0;
  334             }
  335           }
  336         }
  337         if(constrain != NULL)
  338         {
  339           constrain(v[vg], n);
  340         }
  341         f[vg] = objfunc(v[vg], params);
  342         k++;
  343         if(constrain != NULL)
  344         {
  345           constrain(v[vh], n);
  346         }
  347         f[vh] = objfunc(v[vh], params);
  348         k++;
  349       }
  350     }
  351 #if 0
  352     /* print out the value at each iteration */
  353     printf("Iteration %d\n", itr);
  354     for(j = 0; j <= n; j++)
  355     {
  356       for(i = 0; i < n; i++)
  357       {
  358         printf("%f %f\n", v[j][i], f[j]);
  359       }
  360     }
  361 #endif
  362     /* test for convergence */
  363     fsum = 0.0;
  364     for(j = 0; j <= n; j++)
  365     {
  366       fsum += f[j];
  367     }
  368     favg = fsum / (n + 1);
  369     s = 0.0;
  370     for(j = 0; j <= n; j++)
  371     {
  372       s += pow((f[j] - favg), 2.0) / (n);
  373     }
  374     s = sqrt(s);
  375     if(s < EPSILON) break;
  376   }
  377   /* end main loop of the minimization */
  378 
  379   /* find the index of the smallest value */
  380   vs = 0;
  381   for(j = 0; j <= n; j++)
  382   {
  383     if(f[j] < f[vs])
  384     {
  385       vs = j;
  386     }
  387   }
  388 #if 0
  389   printf("The minimum was found at\n");
  390   for(j = 0; j < n; j++)
  391   {
  392     printf("%e\n", v[vs][j]);
  393     start[j] = v[vs][j];
  394   }
  395   double min = objfunc(v[vs], params);
  396   k++;
  397   printf("The minimum value is %f\n", min);
  398   printf("%d Function Evaluations\n", k);
  399   printf("%d Iterations through program\n", itr);
  400 #else
  401   for(j = 0; j < n; j++)
  402   {
  403     start[j] = v[vs][j];
  404   }
  405 #endif
  406   free(f);
  407   free(vr);
  408   free(ve);
  409   free(vc);
  410   free(vm);
  411   for(i = 0; i <= n; i++)
  412   {
  413     free(v[i]);
  414   }
  415   free(v);
  416   return itr;
  417 }
  418 
  419 /*==================================================================================
  420  * end of nmsimplex code
  421  *==================================================================================*/
  422 
  423 // modelines: These editor modelines have been set for all relevant files by tools/update_modelines.sh
  424 // vim: shiftwidth=2 expandtab tabstop=2 cindent
  425 // kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;