Gaussian HMC - integration within cooling gaussians and complexity improvements #301
Conversation
|
Thanks @vgnecula! This is a useful PR. Could you please change the title of the PR to represent the actual changes that this PR is carrying? |
vgnecula
changed the title
| @@ -457,7 +455,7 @@ double volume_cooling_gaussians(Polytope& Pin, | |||
|
|
|||
| template | |||
| < | |||
| typename WalkTypePolicy = GaussianCDHRWalk, | |||
There was a problem hiding this comment.
Please don't change the default random walk in volume_cooling_gaussians<>.
You can create a new example in folder examples that it will call this function templated with GaussianHamiltonianMonteCarloExactWalk. You can write 3-4 examples that compute volumes of polytopes from 10 to 50 dimensions. You can use the volesti polytope generators.
| endif () | ||
|
|
||
|
|
||
| add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard |
| typedef BoostRandomNumberGenerator<boost::mt19937, NT, 3> RNGType; | ||
| typedef HPolytope<Point> HPOLYTOPE; | ||
|
|
||
| void calculateAndVerifyVolume(HPOLYTOPE& polytope, const std::string& description, NT expectedVolume, NT epsilon, NT L) { |
There was a problem hiding this comment.
We are not very consistent but we use lowercase and "_" for function names.
|
|
||
| // birkhoff polytope dimension 3 | ||
| HPOLYTOPE birkhoffPolytope = generate_birkhoff<HPOLYTOPE>(3); | ||
| calculateAndVerifyVolume(birkhoffPolytope, "Birkhoff polytope (dim 3)", -1, epsilon, 8); // Theoretical volume not easily available |
There was a problem hiding this comment.
Note that the Birkhoff polytope B(3) is not 3-dimensional, but 4-dimensional. In general B(n) is (n-1)^2 dimensional, see https://en.wikipedia.org/wiki/Birkhoff_polytope
Also (theoretical) volumes for n<=10 are available here https://oeis.org/A037302 (could be useful for experiments).
| p.set_coord(i, C * std::cos(omega * T + Phi)); | ||
| v.set_coord(i, -C * omega * std::sin(omega * T + Phi)); | ||
| } | ||
| inline void update_position(Point& p, Point& v, const NT& T, const NT& omega) { |
There was a problem hiding this comment.
I think I prefer the constant reference of type T to be written as T const& is there any reason for rewriting as const T&?
| @@ -487,4 +498,4 @@ double volume_cooling_gaussians(Polytope &Pin, | |||
| return volume_cooling_gaussians<WalkTypePolicy>(Pin, rng, error, walk_length); | |||
| } | |||
|
|
|||
| #endif // VOLUME_COOLING_GAUSSIANS_HPP | |||
| #endif // VOLUME_COOLING_GAUSSIANS_HPP | |||
There was a problem hiding this comment.
please keep the empty ending line
|
|
||
| update_delta<WalkType> |
There was a problem hiding this comment.
why are you removing update_delta?
| #endif | ||
|
|
||
| NT a_stop = 0.0; | ||
| const NT tol = 0.001; | ||
| unsigned int it = 0; | ||
| unsigned int n = P.dimension(); | ||
| const unsigned int totalSteps = ((int)150/((1.0 - frac) * error))+1; | ||
|
|
||
| if (a_vals[0]<a_stop) a_vals[0] = a_stop; |
There was a problem hiding this comment.
Why are you removing this check? If the computed (first) gaussian has negative variance then we set it to zero. @TolisChal is this check actually needed?
| curr_fn += eval_exp(p, next_a) / eval_exp(p, a_vals[it]); | ||
| for (unsigned int j = 0; j < totalSteps; ++j) { | ||
| walk.template apply(P, p, a_vals[it], walk_length, rng); | ||
| curr_its++; |
There was a problem hiding this comment.
this not an integer, although I do not know the reason (@TolisChal ?)
| // Set the last a_i equal to zero | ||
| if (next_a>0 && curr_fn/curr_its>(1.0+tol)) | ||
| { | ||
| // pick the last gaussian |
There was a problem hiding this comment.
This is the opposite of the last comment, why? It seems that the code indeed remove the last Gaussian (found by that ratio tolerance test) by setting its variance equals to zero.
There was a problem hiding this comment.
I don't see any reason to delete this check. @vgnecula could you please comment on that?
There was a problem hiding this comment.
Hmm, I don't know why I modified the comment, as the code seems to do the exact same thing. It probably remained like that from when I was trying to check a bug. Sorry!
…elta update on each step
Integrate the Exact Gaussian HMC in cooling gaussian