added lorentz peak bench function in examples

examples/Makefile.am

@@ -1,7 +1,8 @@
-bin_PROGRAMS=sample_code sample_code_genopheno sample_code_pfunc
+bin_PROGRAMS=sample_code sample_code_genopheno sample_code_pfunc lorentzpeakbench
 sample_code_SOURCES=sample-code.cc
 sample_code_genopheno_SOURCES=sample-code-genopheno.cc
 sample_code_pfunc_SOURCES=sample-code-pfunc.cc
+lorentzpeakbench_SOURCES=lorentzpeakbench.cc
 
 AM_CPPFLAGS=-I../src/ -I$(EIGEN3_INC)
 AM_CXXFLAGS=-g -Wall -fopenmp

examples/lorentzpeakbench.cc

@@ -0,0 +1,88 @@
+/**
+ * CMA-ES, Covariance Matrix Adaptation Evolution Strategy
+ * Copyright (c) 2014 Inria
+ * Author: Emmanuel Benazera <[email protected]>
+ *
+ * This file is part of libcmaes.
+ *
+ * libcmaes is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * libcmaes is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with libcmaes.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "cmaes.h"
+#include <cmath>
+#include <fstream>
+#include <iostream>
+
+using namespace libcmaes;
+
+double background(const double *x, const double *par)
+{
+  return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
+}
+
+double lorentzianpeak(const double *x, const double *par)
+{
+  return (0.5*par[0]*par[1]/M_PI) /
+    std::max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]);
+}
+
+double blfunc(const double *x, const double *par)
+{
+  return background(x,par) + lorentzianpeak(x,&par[3]);
+}
+
+//double par[6] = {1,1,1,6,.03,1};
+double points[201];
+double values[201];
+
+// chi2
+FitFunc ff = [](const double *x, const int N)
+{
+  double sum = 0.0;
+  for (int i=0;i<201;i++)
+    {
+      //std::cout << "x=" << points[i] << " / f=" << values[i] << std::endl;
+      if (values[i] != 0.0)
+	sum += pow((values[i]-blfunc(&points[i],x))/sqrt(values[i]),2);
+    }
+  //std::cout << "sum=" << sum << std::endl;
+  return sum;
+};
+
+void loaddata(const std::string filename)
+{
+  std::ifstream fin(filename);
+  std::string line;
+  std::getline(fin,line); // skip header line.
+  int a, b, i = 0;
+  double c;
+  while(fin >> a >> b >> c)
+    {
+      points[i] = c;
+      values[i] = b;
+      //std::cout << "x=" << points[i] << " / f=" << values[i] << std::endl;
+      i++;
+    }
+}
+
+int main(int argc, char *argv[])
+{
+  int dim = 6;
+  loaddata("lorentzpeakbench.dat");
+  CMAParameters<> cmaparams(dim);
+  CMASolutions cmasols = cmaes<>(ff,cmaparams);
+  std::cout << "best solution: " << cmasols << std::endl;
+  std::cout << "optimization took " << cmasols._elapsed_time / 1000.0 << " seconds\n";
+  return cmasols._run_status;
+}

examples/lorentzpeakbench.dat

@@ -0,0 +1,203 @@
+bin val x
+0 0 -0.0075
+1 3 0.0075
+2 6 0.0225
+3 3 0.0375
+4 6 0.0525
+5 9 0.0675
+6 3 0.0825
+7 2 0.0975
+8 2 0.1125
+9 3 0.1275
+10 5 0.1425
+11 2 0.1575
+12 3 0.1725
+13 4 0.1875
+14 3 0.2025
+15 4 0.2175
+16 7 0.2325
+17 1 0.2475
+18 5 0.2625
+19 7 0.2775
+20 3 0.2925
+21 6 0.3075
+22 5 0.3225
+23 3 0.3375
+24 8 0.3525
+25 9 0.3675
+26 5 0.3825
+27 2 0.3975
+28 9 0.4125
+29 5 0.4275
+30 7 0.4425
+31 10 0.4575
+32 5 0.4725
+33 6 0.4875
+34 5 0.5025
+35 5 0.5175
+36 3 0.5325
+37 3 0.5475
+38 10 0.5625
+39 6 0.5775
+40 8 0.5925
+41 6 0.6075
+42 8 0.6225
+43 3 0.6375
+44 15 0.6525
+45 8 0.6675
+46 9 0.6825
+47 7 0.6975
+48 3 0.7125
+49 10 0.7275
+50 10 0.7425
+51 13 0.7575
+52 12 0.7725
+53 5 0.7875
+54 11 0.8025
+55 7 0.8175
+56 11 0.8325
+57 14 0.8475
+58 17 0.8625
+59 15 0.8775
+60 14 0.8925
+61 25 0.9075
+62 26 0.9225
+63 26 0.9375
+64 49 0.9525
+65 92 0.9675
+66 203 0.9825
+67 400 0.9975
+68 254 1.0125
+69 96 1.0275
+70 55 1.0425
+71 35 1.0575
+72 35 1.0725
+73 27 1.0875
+74 22 1.1025
+75 23 1.1175
+76 15 1.1325
+77 17 1.1475
+78 11 1.1625
+79 18 1.1775
+80 15 1.1925
+81 12 1.2075
+82 13 1.2225
+83 14 1.2375
+84 20 1.2525
+85 13 1.2675
+86 15 1.2825
+87 14 1.2975
+88 11 1.3125
+89 21 1.3275
+90 20 1.3425
+91 22 1.3575
+92 15 1.3725
+93 13 1.3875
+94 13 1.4025
+95 12 1.4175
+96 10 1.4325
+97 17 1.4475
+98 20 1.4625
+99 7 1.4775
+100 21 1.4925
+101 18 1.5075
+102 10 1.5225
+103 18 1.5375
+104 22 1.5525
+105 17 1.5675
+106 15 1.5825
+107 13 1.5975
+108 24 1.6125
+109 13 1.6275
+110 17 1.6425
+111 17 1.6575
+112 15 1.6725
+113 17 1.6875
+114 25 1.7025
+115 19 1.7175
+116 23 1.7325
+117 23 1.7475
+118 23 1.7625
+119 16 1.7775
+120 22 1.7925
+121 15 1.8075
+122 16 1.8225
+123 18 1.8375
+124 20 1.8525
+125 17 1.8675
+126 11 1.8825
+127 31 1.8975
+128 17 1.9125
+129 23 1.9275
+130 20 1.9425
+131 28 1.9575
+132 32 1.9725
+133 31 1.9875
+134 22 2.0025
+135 28 2.0175
+136 25 2.0325
+137 27 2.0475
+138 23 2.0625
+139 21 2.0775
+140 24 2.0925
+141 24 2.1075
+142 28 2.1225
+143 31 2.1375
+144 23 2.1525
+145 24 2.1675
+146 27 2.1825
+147 30 2.1975
+148 19 2.2125
+149 31 2.2275
+150 25 2.2425
+151 34 2.2575
+152 22 2.2725
+153 19 2.2875
+154 35 2.3025
+155 28 2.3175
+156 34 2.3325
+157 42 2.3475
+158 18 2.3625
+159 31 2.3775
+160 24 2.3925
+161 27 2.4075
+162 36 2.4225
+163 21 2.4375
+164 24 2.4525
+165 21 2.4675
+166 38 2.4825
+167 37 2.4975
+168 32 2.5125
+169 33 2.5275
+170 27 2.5425
+171 37 2.5575
+172 50 2.5725
+173 25 2.5875
+174 31 2.6025
+175 30 2.6175
+176 48 2.6325
+177 29 2.6475
+178 43 2.6625
+179 43 2.6775
+180 29 2.6925
+181 36 2.7075
+182 44 2.7225
+183 42 2.7375
+184 45 2.7525
+185 40 2.7675
+186 32 2.7825
+187 37 2.7975
+188 30 2.8125
+189 42 2.8275
+190 35 2.8425
+191 49 2.8575
+192 26 2.8725
+193 51 2.8875
+194 52 2.9025
+195 40 2.9175
+196 44 2.9325
+197 45 2.9475
+198 35 2.9625
+199 53 2.9775
+200 50 2.9925
+201 0 3.0075