1/*2 This file is part of darktable,3 Copyright (C) 2016-2020 darktable developers.45 darktable is free software: you can redistribute it and/or modify6 it under the terms of the GNU General Public License as published by7 the Free Software Foundation, either version 3 of the License, or8 (at your option) any later version.910 darktable is distributed in the hope that it will be useful,11 but WITHOUT ANY WARRANTY; without even the implied warranty of12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the13 GNU General Public License for more details.1415 You should have received a copy of the GNU General Public License16 along with darktable. If not, see <http://www.gnu.org/licenses/>.17*/1819/* For parameter optimization we are using the Nelder-Mead simplex method20 * 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 parameters25 * simplex() returns number of used iterations instead of min value26 * maximum number of iterations as function parameter27 * make interface function simplex() static28 * initialize i and j to avoid compiler warnings29 * comment out printing of status inormation30 * reformat according to darktable's clang standards31 */3233/*==================================================================================34 * begin nmsimplex code downloaded from http://www.mikehutt.com/neldermead.html35 * on February 6, 201636 *==================================================================================*/37/*38 * Program: nmsimplex.c39 * Author : Michael F. Hutt40 * http://www.mikehutt.com41 * 11/3/9742 *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 obtaining48 * a copy of this software and associated documentation files (the49 * "Software"), to deal in the Software without restriction, including50 * without limitation the rights to use, copy, modify, merge, publish,51 * distribute, sublicense, and/or sell copies of the Software, and to52 * permit persons to whom the Software is furnished to do so, subject to53 * the following conditions:54 *55 * The above copyright notice and this permission notice shall be56 * 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 OF60 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND61 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE62 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION63 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION64 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.65 *66 *67 * Jan. 6, 199968 * Modified to conform to the algorithm presented69 * in Margaret H. Wright's paper on Direct Search Methods.70 *71 * Jul. 23, 200772 * Fixed memory leak.73 *74 * Mar. 1, 201175 * Added constraints.76 */7778//#include "nmsimplex.h"7980static 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{8384 int vs; /* vertex with smallest value */85 int vh; /* vertex with next smallest value */86 int vg; /* vertex with largest value */8788 int i = 0, j = 0, m, row;89 int itr; /* track the number of iterations */9091 double **v; /* holds vertices of simplex */92 double pn, qn; /* values used to create initial simplex */93 double *f; /* value of function at each vertex */94 double fr; /* value of function at reflection point */95 double fe; /* value of function at expansion point */96 double fc; /* value of function at contraction point */97 double *vr; /* reflection - coordinates */98 double *ve; /* expansion - coordinates */99 double *vc; /* contraction - coordinates */100 double *vm; /* centroid - coordinates */101 //double min;102103 double fsum, favg, s, cent;104105 /* dynamically allocate arrays */106107 /* allocate the rows of the arrays */108 v = (double **)malloc(sizeof(double *) * (n + 1));109 f = (double *)malloc(sizeof(double) * (n + 1));110 vr = (double *)malloc(sizeof(double) * n);111 ve = (double *)malloc(sizeof(double) * n);112 vc = (double *)malloc(sizeof(double) * n);113 vm = (double *)malloc(sizeof(double) * n);114115 /* allocate the columns of the arrays */116 for(i = 0; i <= n; i++)117 {118 v[i] = (double *)malloc(sizeof(double) * n);119 }120121 /* create the initial simplex */122 /* assume one of the vertices is 0,0 */123124 pn = scale * (sqrt(n + 1) - 1 + n) / (n * sqrt(2));125 qn = scale * (sqrt(n + 1) - 1) / (n * sqrt(2));126127 for(i = 0; i < n; i++)128 {129 v[0][i] = start[i];130 }131132 for(i = 1; i <= n; i++)133 {134 for(j = 0; j < n; j++)135 {136 if(i - 1 == j)137 {138 v[i][j] = pn + start[j];139 }140 else141 {142 v[i][j] = qn + start[j];143 }144 }145 }146147 if(constrain != NULL)148 {149 constrain(v[j], n);150 }151 /* find the initial function values */152 for(j = 0; j <= n; j++)153 {154 f[j] = objfunc(v[j], params);155 }156157158 /* print out the initial values */159 printf("Initial Values\n");160 for(j = 0; j <= n; j++)161 {162 for(i = 0; i < n; i++)163 {164 printf("%f %f\n", v[j][i], f[j]);165 }166 }167168169 /* begin the main loop of the minimization */170 for(itr = 1; itr <= maxiter; itr++)171 {172 /* find the index of the largest value */173 vg = 0;174 for(j = 0; j <= n; j++)175 {176 if(f[j] > f[vg])177 {178 vg = j;179 }180 }181182 /* find the index of the smallest value */183 vs = 0;184 for(j = 0; j <= n; j++)185 {186 if(f[j] < f[vs])187 {188 vs = j;189 }190 }191192 /* find the index of the second largest value */193 vh = vs;194 for(j = 0; j <= n; j++)195 {196 if(f[j] > f[vh] && f[j] < f[vg])197 {198 vh = j;199 }200 }201202 /* calculate the centroid */203 for(j = 0; j <= n - 1; j++)204 {205 cent = 0.0;206 for(m = 0; m <= n; m++)207 {208 if(m != vg)209 {210 cent += v[m][j];211 }212 }213 vm[j] = cent / n;214 }215216 /* reflect vg to new vertex vr */217 for(j = 0; j <= n - 1; j++)218 {219 /*vr[j] = (1+NMS_ALPHA)*vm[j] - NMS_ALPHA*v[vg][j];*/220 vr[j] = vm[j] + NMS_ALPHA * (vm[j] - v[vg][j]);221 }222 if(constrain != NULL)223 {224 constrain(vr, n);225 }226 fr = objfunc(vr, params);227228 if(fr < f[vh] && fr >= f[vs])229 {230 for(j = 0; j <= n - 1; j++)231 {232 v[vg][j] = vr[j];233 }234 f[vg] = fr;235 }236237 /* investigate a step further in this direction */238 if(fr < f[vs])239 {240 for(j = 0; j <= n - 1; j++)241 {242 /*ve[j] = NMS_GAMMA*vr[j] + (1-NMS_GAMMA)*vm[j];*/243 ve[j] = vm[j] + NMS_GAMMA * (vr[j] - vm[j]);244 }245 if(constrain != NULL)246 {247 constrain(ve, n);248 }249 fe = objfunc(ve, params);250251 /* by making fe < fr as opposed to fe < f[vs],252 Rosenbrocks function takes 63 iterations as opposed253 to 64 when using double variables. */254255 if(fe < fr)256 {257 for(j = 0; j <= n - 1; j++)258 {259 v[vg][j] = ve[j];260 }261 f[vg] = fe;262 }263 else264 {265 for(j = 0; j <= n - 1; j++)266 {267 v[vg][j] = vr[j];268 }269 f[vg] = fr;270 }271 }272273 /* check to see if a contraction is necessary */274 if(fr >= f[vh])275 {276 if(fr < f[vg] && fr >= f[vh])277 {278 /* perform outside contraction */279 for(j = 0; j <= n - 1; j++)280 {281 /*vc[j] = NMS_BETA*v[vg][j] + (1-NMS_BETA)*vm[j];*/282 vc[j] = vm[j] + NMS_BETA * (vr[j] - vm[j]);283 }284 if(constrain != NULL)285 {286 constrain(vc, n);287 }288 fc = objfunc(vc, params);289 }290 else291 {292 /* perform inside contraction */293 for(j = 0; j <= n - 1; j++)294 {295 /*vc[j] = NMS_BETA*v[vg][j] + (1-NMS_BETA)*vm[j];*/296 vc[j] = vm[j] - NMS_BETA * (vm[j] - v[vg][j]);297 }298 if(constrain != NULL)299 {300 constrain(vc, n);301 }302 fc = objfunc(vc, params);303 }304305306 if(fc < f[vg])307 {308 for(j = 0; j <= n - 1; j++)309 {310 v[vg][j] = vc[j];311 }312 f[vg] = fc;313 }314 /* at this point the contraction is not successful,315 we must halve the distance from vs to all the316 vertices of the simplex and then continue.317 10/31/97 - modified to account for ALL vertices.318 */319 else320 {321 for(row = 0; row <= n; row++)322 {323 if(row != vs)324 {325 for(j = 0; j <= n - 1; j++)326 {327 v[row][j] = v[vs][j] + (v[row][j] - v[vs][j]) / 2.0;328 }329 }330 }331 if(constrain != NULL)332 {333 constrain(v[vg], n);334 }335 f[vg] = objfunc(v[vg], params);336 if(constrain != NULL)337 {338 constrain(v[vh], n);339 }340 f[vh] = objfunc(v[vh], params);341 }342 }343344 /* print out the value at each iteration */345 printf("Iteration %d\n", itr);346 for(j = 0; j <= n; j++)347 {348 for(i = 0; i < n; i++)349 {350 printf("%f %f\n", v[j][i], f[j]);351 }352 }353354 /* test for convergence */355 fsum = 0.0;356 for(j = 0; j <= n; j++)357 {358 fsum += f[j];359 }360 favg = fsum / (n + 1);361 s = 0.0;362 for(j = 0; j <= n; j++)363 {364 s += pow((f[j] - favg), 2.0) / (n);365 }366 s = sqrt(s);367 if(s < EPSILON) break;368 }369 /* end main loop of the minimization */370371 /* find the index of the smallest value */372 vs = 0;373 for(j = 0; j <= n; j++)374 {375 if(f[j] < f[vs])376 {377 vs = j;378 }379 }380381 printf("The minimum was found at\n");382 for(j = 0; j < n; j++)383 {384 printf("%e\n", v[vs][j]);385 start[j] = v[vs][j];386 }387 double min = objfunc(v[vs], params);388 k++;389 printf("The minimum value is %f\n", min);390 printf("%d Function Evaluations\n", k);391 printf("%d Iterations through program\n", itr);392393 for(j = 0; j < n; j++)394 {395 start[j] = v[vs][j];396 }397398 free(f);399 free(vr);400 free(ve);401 free(vc);402 free(vm);403 for(i = 0; i <= n; i++)404 {405 free(v[i]);406 }407 free(v);408 return itr;409}410411/*==================================================================================412 * end of nmsimplex code413 *==================================================================================*/414415// clang-format off416// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py417// vim: shiftwidth=2 expandtab tabstop=2 cindent418// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified;419// clang-format on420