Dear all,
I am searching for alternatives to mpfit for a C++ program. The optimization subpackage of alglib seems to be a good candidate. Before going any further, I wrote a small example where the same problem is solved using both libraries. I managed to get the same result, but alglib is way slower (factor 8). Before jumping to conclusions, I want to make sure that I haven't missed any important things that could impact on performance. The code is compiled with cmake, and I did nothing specific in terms of compilation flag, except a release mode -O3 flag.
Could an alglib expert see what I'm doing wrong? I have also added a full example in attachement. I had to remove the alglib sources to keep the file under 1Mo, but you just have to copy/paste the sources in the alglib folder and it should work.
If the commercial edition brings a performance boost, I would definitely consider bying it.
Here is the code:
Code:
#include <algorithm>
#include <cmath>
#include <chrono>
#include <iostream>
#include <vector>
// alglib
#include "src/ap.h"
#include "src/stdafx.h"
#include "src/optimization.h"
struct vars_struct
{
double *x;
double *y;
};
void planck_alglib(const alglib::real_1d_array &p, alglib::real_1d_array &vy, void *ptr)
{
struct vars_struct *v = (struct vars_struct *) ptr;
for (int i=0; i<vy.length(); i++)
vy[i] = ((p[0] / (std::exp(p[1] / v->x[i]) - 1) + p[2]) - v->y[i]);
}
void planck_alglib_jac(const alglib::real_1d_array &p, alglib::real_1d_array &vy, alglib::real_2d_array &jac, void *ptr)
{
struct vars_struct *v = (struct vars_struct *) ptr;
double B_X;
for (int i=0; i<vy.length(); i++) {
B_X = std::exp(p[1] / v->x[i]);
vy[i] = (p[0] / (B_X - 1) + p[2]) - v->y[i];
jac[i][0] = 1 / (B_X - 1);
jac[i][1] = - p[0] * B_X / ((B_X - 1) * (B_X - 1) * v->x[i]);
jac[i][2] = 1;
}
}
int main() {
// 50 data points
std::vector<double> T = {
-5.38739, -3.7718, -2.23958, -0.78166, 0.60961, 1.94071, 3.21719, 4.44385,
5.62487, 6.76392, 7.86421, 8.9286, 9.95963, 10.95957, 11.93047, 12.87417,
13.79233, 14.68647, 15.55796, 16.40807, 17.23796, 18.04867, 18.84118, 19.61641,
20.37518, 21.11827, 21.84639, 22.56021, 23.26036, 23.94741, 24.62192, 25.2844,
25.93531, 26.57511, 27.20423, 27.82306, 28.43198, 29.03135, 29.62149, 30.20273,
30.77538, 31.3397, 31.89598, 32.44447, 32.98542, 33.51906, 34.04562, 34.5653,
35.0783, 35.58483};
std::transform(T.begin(), T.end(), T.begin(), [&](double& value){ return value + 273.15; });
std::vector<double> DV = {
3762, 3910, 4058, 4206, 4354, 4502, 4650, 4798, 4946,5094, 5242, 5390, 5538,
5686, 5834, 5982, 6130, 6278, 6426, 6574, 6722, 6870, 7018, 7166, 7314, 7462,
7610, 7758, 7906, 8054, 8202, 8350, 8498, 8646, 8794, 8942, 9090, 9238, 9386,
9534, 9682, 9830, 9978, 10126, 10274, 10422, 10570, 10718, 10866, 11014};
struct vars_struct v;
v.x = T.data();
v.y = DV.data();
alglib::real_1d_array s = "[100000, 1, 1]";
alglib::real_1d_array p0 = "[100000000, 1000, 1000]";
alglib::minlmstate state;
double epsx = 0.000001;
alglib::ae_int_t maxits = 200;
alglib::minlmreport rep;
alglib::minlmcreatevj(3, 50, p0, state);
alglib::minlmsetcond(state, epsx, maxits);
alglib::minlmsetscale(state, s);
begin = std::chrono::steady_clock::now();
// we simulate here various optimization, to check if minlmrestartfrom can help us gain some time
for (int n = 0; n < 100; n++) {
p0[0] = 100000000;
p0[1] = 1000;
p0[2] = 1000;
alglib::minlmrestartfrom(state, p0);
alglib::minlmoptimize(state, planck_alglib, planck_alglib_jac, nullptr, (void*)&v);
minlmresults(state, p0, rep);
}
end = std::chrono::steady_clock::now();
std::cout << "R: " << p0[0] << std::endl; // R: 1.476e+08
std::cout << "B: " << p0[1] << std::endl; // B: 2983.1
std::cout << "O: " << p0[2] << std::endl; // O: 1620.77
std::cout << "duration : " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << " ms" << std::endl;
return 0;
}
#include <cmath>
#include <chrono>
#include <iostream>
#include <vector>
// alglib
#include "src/ap.h"
#include "src/stdafx.h"
#include "src/optimization.h"
struct vars_struct
{
double *x;
double *y;
};
void planck_alglib(const alglib::real_1d_array &p, alglib::real_1d_array &vy, void *ptr)
{
struct vars_struct *v = (struct vars_struct *) ptr;
for (int i=0; i<vy.length(); i++)
vy[i] = ((p[0] / (std::exp(p[1] / v->x[i]) - 1) + p[2]) - v->y[i]);
}
void planck_alglib_jac(const alglib::real_1d_array &p, alglib::real_1d_array &vy, alglib::real_2d_array &jac, void *ptr)
{
struct vars_struct *v = (struct vars_struct *) ptr;
double B_X;
for (int i=0; i<vy.length(); i++) {
B_X = std::exp(p[1] / v->x[i]);
vy[i] = (p[0] / (B_X - 1) + p[2]) - v->y[i];
jac[i][0] = 1 / (B_X - 1);
jac[i][1] = - p[0] * B_X / ((B_X - 1) * (B_X - 1) * v->x[i]);
jac[i][2] = 1;
}
}
int main() {
// 50 data points
std::vector<double> T = {
-5.38739, -3.7718, -2.23958, -0.78166, 0.60961, 1.94071, 3.21719, 4.44385,
5.62487, 6.76392, 7.86421, 8.9286, 9.95963, 10.95957, 11.93047, 12.87417,
13.79233, 14.68647, 15.55796, 16.40807, 17.23796, 18.04867, 18.84118, 19.61641,
20.37518, 21.11827, 21.84639, 22.56021, 23.26036, 23.94741, 24.62192, 25.2844,
25.93531, 26.57511, 27.20423, 27.82306, 28.43198, 29.03135, 29.62149, 30.20273,
30.77538, 31.3397, 31.89598, 32.44447, 32.98542, 33.51906, 34.04562, 34.5653,
35.0783, 35.58483};
std::transform(T.begin(), T.end(), T.begin(), [&](double& value){ return value + 273.15; });
std::vector<double> DV = {
3762, 3910, 4058, 4206, 4354, 4502, 4650, 4798, 4946,5094, 5242, 5390, 5538,
5686, 5834, 5982, 6130, 6278, 6426, 6574, 6722, 6870, 7018, 7166, 7314, 7462,
7610, 7758, 7906, 8054, 8202, 8350, 8498, 8646, 8794, 8942, 9090, 9238, 9386,
9534, 9682, 9830, 9978, 10126, 10274, 10422, 10570, 10718, 10866, 11014};
struct vars_struct v;
v.x = T.data();
v.y = DV.data();
alglib::real_1d_array s = "[100000, 1, 1]";
alglib::real_1d_array p0 = "[100000000, 1000, 1000]";
alglib::minlmstate state;
double epsx = 0.000001;
alglib::ae_int_t maxits = 200;
alglib::minlmreport rep;
alglib::minlmcreatevj(3, 50, p0, state);
alglib::minlmsetcond(state, epsx, maxits);
alglib::minlmsetscale(state, s);
begin = std::chrono::steady_clock::now();
// we simulate here various optimization, to check if minlmrestartfrom can help us gain some time
for (int n = 0; n < 100; n++) {
p0[0] = 100000000;
p0[1] = 1000;
p0[2] = 1000;
alglib::minlmrestartfrom(state, p0);
alglib::minlmoptimize(state, planck_alglib, planck_alglib_jac, nullptr, (void*)&v);
minlmresults(state, p0, rep);
}
end = std::chrono::steady_clock::now();
std::cout << "R: " << p0[0] << std::endl; // R: 1.476e+08
std::cout << "B: " << p0[1] << std::endl; // B: 2983.1
std::cout << "O: " << p0[2] << std::endl; // O: 1620.77
std::cout << "duration : " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << " ms" << std::endl;
return 0;
}
I woud really appreciate if someone could help me understand what is going on. Many thanks in advance
Vincent
| Attachments: |
ful_example.zip [61.19 KiB]
Downloaded 129 times |