Skip to content

Commit

Permalink
[fit] improve default parameters
Browse files Browse the repository at this point in the history
  • Loading branch information
@KernerM
KernerM committed Nov 21, 2018
1 parent 032b3dc commit 09850a2
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 17 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -287,7 +287,7 @@ Default: `3e-3`.

See `gtol` in the nonlinear least-squares fitting GSL documentation.

Default: `1e-4`.
Default: `1e-8`.

---

Expand Down
34 changes: 19 additions & 15 deletions src/fitfunctions/polysingular.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,17 @@ namespace integrators
struct PolySingularFunction
{
static const int num_parameters = 6;
const std::vector<std::vector<D>> initial_parameters = { {1.1,-0.1, 0.1,0.1, 1.0,0.0} };
const std::vector<std::vector<D>> initial_parameters = { {1.1,-0.1, 0.1,0.1, 0.3,0.3} };

D operator()(const D x, const double* p) const
{
// constraint: no singularity and singular terms have positive coefficients
if (p[0]<=static_cast<D>(1.001) or p[0]>=static_cast<D>(1.5) or p[1]>=static_cast<D>(-0.001) or p[1]<=static_cast<D>(-0.5) or p[2]<static_cast<D>(0) or p[3]<static_cast<D>(0))
return std::numeric_limits<D>::max();

D y = p[2]*(x*(p[0]-D(1)))/(p[0]-x) + p[3]*(x*(p[1]-D(1)))/(p[1]-x) + x*(p[4]+x*(p[5]+x*(D(1)-p[2]-p[3]-p[4]-p[5])));
if (p[0]<=static_cast<D>(1.001) or p[0]>=static_cast<D>(5) or p[1]>=static_cast<D>(-0.001) or p[1]<=static_cast<D>(-4) )
return D(10.); // std::numeric_limits<D>::max() will sometimes result in fit parameters being NaN

D p2 = abs(p[2]);
D p3 = abs(p[3]);
D y = p2*(x*(p[0]-D(1)))/(p[0]-x) + p3*(x*(p[1]-D(1)))/(p[1]-x) + x*(p[4]+x*(p[5]+x*(D(1)-p2-p3-p[4]-p[5])));

// constraint: transformed variable within unit hypercube
if ( y<static_cast<D>(0) || y>static_cast<D>(1) )
Expand All @@ -39,13 +41,13 @@ namespace integrators
{

if (parameter == 0) {
return p[2]*((D(1) - x)*x)/(x - p[0])/(x - p[0]);
return abs(p[2])*((D(1) - x)*x)/(x - p[0])/(x - p[0]);
} else if (parameter == 1) {
return p[3]*((D(1) - x)*x)/(x - p[1])/(x - p[1]);
return abs(p[3])*((D(1) - x)*x)/(x - p[1])/(x - p[1]);
} else if (parameter == 2) {
return (x*(p[0]-D(1)))/(p[0]-x) -x*x*x;
return ((x*(p[0]-D(1)))/(p[0]-x) -x*x*x) * ((p[2] < 0) ? D(-1) : D(1));
} else if (parameter == 3) {
return (x*(p[1]-D(1)))/(p[1]-x) -x*x*x;
return ((x*(p[1]-D(1)))/(p[1]-x) -x*x*x) * ((p[3] < 0) ? D(-1) : D(1));
} else if (parameter == 4) {
return x*(D(1)-x*x);
} else if (parameter == 5) {
Expand All @@ -60,10 +62,10 @@ namespace integrators
{
D operator()(const D x, const double* v, const double* p) const
{
D xmp0 = x-p[0];
D xmp1 = x-p[1];
return x*(D(1)-x)*D(2)* ( v[0]*(p[2]*v[0]+(x - p[0])*v[2])/xmp0/xmp0/xmp0 + v[1]*(p[3]*v[1]+(x - p[1])*v[3])/xmp1/xmp1/xmp1 );
}
D xmp0 = x-p[0];
D xmp1 = x-p[1];
return x*(D(1)-x)*D(2)*(v[0]*(abs(p[2])*v[0]+(x - p[0])*v[2])/xmp0/xmp0/xmp0 + v[1]*(abs(p[3])*v[1]+(x - p[1])*v[3])/xmp1/xmp1/xmp1);
}
};

template<typename I, typename D, U M>
Expand All @@ -85,8 +87,10 @@ namespace integrators
D wgt = 1;
for (U d = 0; d < number_of_integration_variables ; ++d)
{
wgt *= p[d][2]*p[d][0]*(p[d][0]-D(1))/(p[d][0]-x[d])/(p[d][0]-x[d]) + p[d][3]*p[d][1]*(p[d][1]-D(1))/(p[d][1]-x[d])/(p[d][1]-x[d]) + p[d][4] + x[d]*(D(2)*p[d][5]+x[d]*D(3)*(D(1)-p[d][2]-p[d][3]-p[d][4]-p[d][5]));
x[d] = p[d][2]*(x[d]*(p[d][0]-D(1)))/(p[d][0]-x[d]) + p[d][3]*(x[d]*(p[d][1]-D(1)))/(p[d][1]-x[d]) + x[d]*(p[d][4]+x[d]*(p[d][5]+x[d]*(D(1)-p[d][2]-p[d][3]-p[d][4]-p[d][5])));
D p2 = abs(p[d][2]);
D p3 = abs(p[d][3]);
wgt *= p2*p[d][0]*(p[d][0]-D(1))/(p[d][0]-x[d])/(p[d][0]-x[d]) + p3*p[d][1]*(p[d][1]-D(1))/(p[d][1]-x[d])/(p[d][1]-x[d]) + p[d][4] + x[d]*(D(2)*p[d][5]+x[d]*D(3)*(D(1)-p2-p3-p[d][4]-p[d][5]));
x[d] = p2*(x[d]*(p[d][0]-D(1)))/(p[d][0]-x[d]) + p3*(x[d]*(p[d][1]-D(1)))/(p[d][1]-x[d]) + x[d]*(p[d][4]+x[d]*(p[d][5]+x[d]*(D(1)-p2-p3-p[d][4]-p[d][5])));
if ( x[d] > D(1) || x[d] < D(0) ) return D(0);
}
return wgt * f(x);
Expand Down
2 changes: 1 addition & 1 deletion src/members.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -720,7 +720,7 @@ namespace integrators

template <typename T, typename D, U M, template<typename,typename,U> class P, template<typename,typename,U> class F, typename G, typename H>
Qmc<T,D,M,P,F,G,H>::Qmc() :
logger(std::cout), randomgenerator( G( std::random_device{}() ) ), minn(8191), minm(32), epsrel(0.01), epsabs(1e-7), maxeval(1000000), maxnperpackage(1), maxmperpackage(1024), errormode(integrators::ErrorMode::all), cputhreads(std::thread::hardware_concurrency()), cudablocks(1024), cudathreadsperblock(256), devices({-1}), generatingvectors(integrators::generatingvectors::cbcpt_dn1_100()), verbosity(0), evaluateminn(100000), fitstepsize(10), fitmaxiter(40), fitxtol(3e-3), fitgtol(1e-4), fitftol(1e-8), fitparametersgsl({})
logger(std::cout), randomgenerator( G( std::random_device{}() ) ), minn(8191), minm(32), epsrel(0.01), epsabs(1e-7), maxeval(1000000), maxnperpackage(1), maxmperpackage(1024), errormode(integrators::ErrorMode::all), cputhreads(std::thread::hardware_concurrency()), cudablocks(1024), cudathreadsperblock(256), devices({-1}), generatingvectors(integrators::generatingvectors::cbcpt_dn1_100()), verbosity(0), evaluateminn(100000), fitstepsize(10), fitmaxiter(40), fitxtol(3e-3), fitgtol(1e-8), fitftol(1e-8), fitparametersgsl({})
{
// Check U satisfies requirements of mod_mul implementation
static_assert( std::numeric_limits<U>::is_modulo, "Qmc integrator constructed with a type U that is not modulo. Please use a different unsigned integer type for U.");
Expand Down

0 comments on commit 09850a2

Please sign in to comment.