Implement feature: adjustable maximum displacement in Monte-Carlo #52

g-bauer · 2016-11-01T20:15:57Z

This implementation addresses issue #23 .

There is still work to do to make the changes work with the lumol input crate.

I'd appreciate feedback/critique/suggestions and questions.

Luthaf · 2016-11-01T21:30:34Z

src/core/src/sim/mc/moves/mod.rs

+    /// Current sample range
+    pub delta_range: Range<f64>,
+    /// Number of move calls
+    pub ncalls: u64, 


What is the difference between ncall and nattempted?

nattempted (as well as naccepted) will be set to zero after an "successful" update. ncalls continues to count. Currently, ncalls is only used to check if the displacement should be updated (and the information about the total numbers of calls to the move). We could also use it in conjunction with a nequilibration variable which automatically switches off updates after a certain number of trials.

We could eliminate ncalls and compare the frequency directly to nattempts in update.

Luthaf · 2016-11-01T21:30:55Z

src/core/src/sim/mc/moves/mod.rs

+    /// Current maximum displacement
+    delta: f64,
+    /// Current sample range
+    pub delta_range: Range<f64>,


this could just be called range

Luthaf · 2016-11-01T21:33:43Z

src/core/src/sim/mc/moves/mod.rs

+        }
+    }
+
+    /// Change `update_frequency` to zero, which will prevent further updates.


We should also have a way to change the update frequency (pub fn update_frequency(&mut self, frequency: usize)) here.

I would also prefer the documentation not to refer internal implementation details. How about:

Disable future updates of the move amplitude, until a call to update_frequency is done with a non zero value

I thought about that as well. The question is, where should the user call the function? After adding the move to MonteCarlo and while running the simulation, the user currently has no control over the states of the moves.

I will change the documentation as you requested.

It could be called before adding the move to the Monte-Carlo simulation. After adding it, you are right that there is no more control over the moves. See also how MonteCarlo::add send an error if a move is added after the simulation is started.

Luthaf · 2016-11-01T21:35:46Z

src/core/src/sim/mc/moves/mod.rs

+        let quotient = current_acceptance / self.target_acceptance;
+        match quotient {
+            _ if quotient > 1.5 => self.delta *= 1.5,
+            _ if quotient < 0.5 => self.delta *= 0.5,


50% seems to be a very big change. I have seen 10% in other software, do you have a specific reason to pick 50% ?

Not really, I use it in my code. It only limits the change though. We can easily run tests on what would be the most efficient number here (which would heavily depend on the system I guess).

With a large value, initially there will be very large changes which can lead to faster convergence to the target acceptance.

Luthaf · 2016-11-01T21:36:09Z

src/core/src/sim/mc/moves/rotate.rs

@@ -5,9 +5,11 @@ use rand::Rng;

 use std::usize;
 use std::f64;
+use std::f64::consts::{PI};


style: no need for brackets here

Luthaf · 2016-11-01T21:39:11Z

src/core/src/sim/mc/moves/rotate.rs

 }

 impl Rotate {
    /// Create a new `Rotate` move, with maximum angular displacement of `theta`,
    /// rotating all the molecules in the system.
-    pub fn new(theta: f64) -> Rotate {
-        Rotate::create(theta, None)
+    pub fn new(theta: f64, update_freq: u64, target_acceptance: f64) -> Rotate {


I would not add all the parameters to the constructor, but rather use the default and have setters when deviating from the default.

Something like this:

pub fn new(theta: f64) -> Rotate { // same as before, using Displacement::default for the initial value of displacement } pub fn set_target_acceptance(&mut self, target: f64) { /* ... */ } pub fn set_update_frequency(&mut self, frequency: f64) { /* ... */ }

I see. You would then create a new move and set everything explicitly before adding it to the propagator?

Yes, this feels more clean to me.

Luthaf · 2016-11-01T21:41:29Z

src/core/src/sim/mc/moves/rotate.rs

+/// Returns a random unit vector,
+/// i.e. a random point on an unit sphere.
+/// @TODO: Maybe move this to types::Vector3D?
+fn random_unit_vector3d (rng: &mut Box<Rng>) -> Vector3D {


The algorithm changed here. Why not use the nice 3D normal distribution => uniform distribution on the R3 sphere feature?

Luthaf · 2016-11-01T22:09:49Z

From the design point of view, I still prefer to have the MonteCarlo struct
manage the moves instead of all the moves having a Displacement field. Three
reason for that:

update is just a delegation to displacement.update() in the current code;
you forgot to update the displacement counter in Rotate, so other will make the same mistake and forget this too;
for Hybrid MC (Implement Hybrid Monte-Carlo #16), the adjustable parameter is not the amplitude of the move, but the number of steps of the MD run (an integer and not a float).

I would go for something like the following.

PS: This was going to be small, but turned into a huge pile of code. Ask me if you
do not get what I wrote, and tell me if you disagree!

So we start with a counter for moves frequency:

struct MoveCounter {
    pub ncalls: u64,
    pub naccepted: u64,
    pub nattempted: u64,
    target_acceptance: f64,
}

Then the MCMove win a new method to update itself in a given direction

pub enum UpdateDirection {
    /// We need the acceptance to raise by the given ratio
    More(f64),
    /// We need the acceptance to lower by the given ratio
    Less(f64),
    /// The acceptance does not need to change
    Same,
}

pub trait MCMove {
    // ...

    fn update_acceptance(&mut self, direction: UpdateDirection);
}

Method that can be implemented like this

impl MCMove for Translate {
    // ...
    fn update_acceptance(&mut self, direction: UpdateDirection) {
        match direction {
            UpdateDirection::Same => {return}
            UpdateDirection::More(mut ratio) => {
                debug_assert!(ratio > 1.0);
                if ratio > 1.5 {ratio = 1.5}
                self.delta *= ratio;
            }
            UpdateDirection::Less(mut ratio) => {
                debug_assert!(ratio < 1.0);
                assert!(ratio > 0.0);
                if ratio < 0.5 {ratio = 0.5}
                self.delta *= ratio;
            }
        }
        self.range = Range(-self.delta, self.delta);
    }
}

And finally we can use this in the MonteCarlo struct

pub struct MonteCarlo {
    // ...
    /// List of possible Monte-Carlo moves and associated counters
    moves: Vec<(Box<MCMove>, MoveCounter)>,
    /// Use the same update frequency for all the moves
    update_frequency: u64,
    /// Number of moves counter
    nmoves: u64,
}

impl Propagator for MonteCarlo {
    fn propagate(&mut self, system: &mut System) {
        let mcmove = {
            let probability = self.rng.next_f64();
            // Get the index of the first move with frequency >= probability.
            let (i, _) = self.frequencies.iter()
                                         .enumerate()
                                         .find(|&(_, f)| probability <= *f)
                                         .expect("Could not find a move in MonteCarlo moves list");
            &mut self.moves[i]
        };

        // Update the counters
        mcmove.1.nattempted += 1;
        self.nmoves += 1;

        trace!("Selected move is '{}'", mcmove.0.describe());

        if !mcmove.0.prepare(system, &mut self.rng) {
            trace!("    --> Can not perform the move");
            return;
        }

        let cost = mcmove.0.cost(system, self.beta, &mut self.cache);
        trace!("    --> Move cost is {}", cost);

        let accepted = cost <= 0.0 || self.rng.next_f64() < f64::exp(-cost);

        if accepted {
            // Update the counters
            mcmove.1.naccepted += 1;
            trace!("    --> Move was accepted");
            mcmove.0.apply(system);
            self.cache.update(system);
        } else {
            trace!("    --> Move was rejected");
            mcmove.0.restore(system);
        }

        // Do the adjustments in moves amplitude as needed
        if self.nmoves == self.update_frequency {
            for mcmove in &mut self.moves {
                let direction = mcmove.1.direction();
                mcmove.0.update_acceptance(direction);
            }
            self.nmoves = 0;
        }
    }
}

g-bauer · 2016-11-01T22:48:32Z

All you wrote makes sense. Shifting counting into propagate will certainly save some troubles.

for Hybrid MC (#16), the adjustable parameter is not the amplitude of the move, but the number of steps of the MD run (an integer and not a float).

That is true, future moves could have all kinds of adjustable parameters (number of configurational bias moves, number of MD steps,...). But aren't you shifting the problem towards UpdateDirection where More and Less will encapsulate floats?

How do we progress? I could implement your suggestions and commit the updates.

Luthaf · 2016-11-01T22:56:46Z

That is true, future moves could have all kinds of adjustable parameters (number of configurational bias moves, number of MD steps,...). But aren't you shifting the problem towards UpdateDirection where More and Less will encapsulate floats?

Yeah maybe. I think we need to pass the acceptance ratio, because some moves will need it. if other moves do not need the ratio, they always can ignore the value. But the fact that I needed to write debug_assert!(ratio < 1.0); assert!(ratio > 0.0); in the match do not please me. I'll think about a solution here.

How do we progress? I could implement your suggestions and commit the updates.

Go for it! We always can rebase this branch later to clean the git history. I'll also have some time on Tuesday evening if you need some help.

g-bauer · 2016-11-15T00:08:29Z

Still WIP. Maybe you can have a look at the changes. I decided to compute the scaling factor for the updates as method of MoveCounter. Still not very good for the general case (i.e. decide the number of MD steps or configurational bias moves according to acceptance). That way, it is completely decoupled from the move implementation and so I put it in monte_carlo.rs.

Other option: hand over the acceptance ratio and the target acceptance ...

TODO:

Better documentation
How to set update_frequency in MonteCarlo?
Communication with input
write tests

Luthaf

I like this version better!

I just have some comments about naming improvement, tell me what you think about them.

Luthaf · 2016-11-15T13:24:33Z

src/core/src/sim/mc/monte_carlo.rs

        self.frequencies.push(frequency);
    }

+    pub fn add_with_acceptance(&mut self, mcmove: Box<MCMove>, frequency: f64, target: Option<f64>) {


I feel like the name of this function should start with add_move, but I can not find a nice one. I though about add_move_with_acceptance, add_move_auto_update, ...

What do you think?

Luthaf · 2016-11-15T13:25:13Z

src/core/src/sim/mc/monte_carlo.rs

        self.frequencies.push(frequency);
    }

+    pub fn add_with_acceptance(&mut self, mcmove: Box<MCMove>, frequency: f64, target: Option<f64>) {


The target parameter should be a f64 here. there is already a function for passing None.

Luthaf · 2016-11-15T13:30:50Z

src/core/src/sim/mc/moves/mod.rs

    /// `prepare`.
    fn restore(&mut self, system: &mut System);
+
+    /// Update the displacement
+    fn change_range(&mut self, scaling_factor: Option<f64>);


This could have a better name, how about something like update_amplitude(&mut self, factor: Option<f64>).

Also, what does the Option means here? If it is None, no changes should be performed?
EDIT: looks like this. Can you confirm?

Yes you are right.

Luthaf · 2016-11-15T13:31:14Z

src/core/src/sim/mc/moves/rotate.rs

    /// Range distribution, for generation of the angle
-    angle_rng: Range<f64>,
+    theta_range: Range<f64>,


This could just be range

Luthaf · 2016-11-15T13:31:58Z

src/core/src/sim/mc/moves/rotate.rs

@@ -153,4 +163,4 @@ fn rotation_matrix(axis: &Vector3D, angle: f64) -> Matrix3 {
    );

    return rotation;
-}
+}


Please leave an empty new line at the end of the files.

Luthaf · 2016-11-15T13:32:38Z

src/core/src/sim/mc/moves/translate.rs

-    /// Translation vector
-    delta: Vector3D,
+    /// Maximum displacement value
+    delta: f64,
    /// Translation range for random number generation
    delta_range: Range<f64>,


This could just be range too

Luthaf · 2016-11-15T13:33:15Z

src/core/src/sim/mc/moves/translate.rs

    pub fn with_moltype(dr: f64, moltype: u64) -> Translate {
        Translate::create(dr, Some(moltype))
    }

    /// Factorizing the constructors
    fn create(dr: f64, moltype: Option<u64>) -> Translate {
        assert!(dr > 0.0, "dr must be positive in Translate move");
-        let dr = dr / f64::sqrt(3.0);


I think you need to keep this division by sqrt(3) to ensure that the vector norm will be lower than the given dr.

Luthaf · 2016-11-15T13:33:51Z

src/core/src/sim/mc/moves/translate.rs

@@ -72,15 +71,15 @@ impl MCMove for Translate {
            return false;
        }

-        self.delta = Vector3D::new(
+        let displace_vec = Vector3D::new(


naming: displacement? delta?

Luthaf · 2016-11-15T13:37:36Z

How to set update_frequency in MonteCarlo?

What is wrong with a simple MonteCarlo::update_frequency(&mut self, frequency: usize) for that?

Communication with input

This can be made in a separated PR. That was part of the idea for separating the crates, you do not always need to update all of them at the same time.

Luthaf · 2016-11-15T13:41:01Z

You will also need to rebase on top of master (as there are conflicts) and squash related commits together. Here is a small guide about it, just use rebase -i master instead of HEAD~3. I also found a longer explanation here.

g-bauer · 2016-11-25T09:42:46Z

I need to further work on the documentation. Skimming through the API doc, there is some confusing stuff. Also, I should add some examples.

Luthaf · 2016-11-25T11:33:35Z

So Travis is failing on the lumol-core tests, here. This is missing the new update_amplitude function.

You can run the tests locally by doing cd src/core; cargo test. This will be easier to do when rust-lang/cargo#3221 lands.

…_carlo.rs

Luthaf · 2016-11-26T14:47:57Z

I need to further work on the documentation. Skimming through the API doc, there is some confusing stuff. Also, I should add some examples.

You will also need to add a few unit tests for the codecov check to pass. For example, you can test that compute_scaling_factor does return None when called too early, and return the right value of the ratio for the given values of self.naccepted and self.nattempted. And same for the other functions.

… tests.

Luthaf · 2016-11-27T15:44:56Z

The codecov diff hit check is good for me, missing path where not previously checked. Do you want to add more doc/examples or should I merge this?

g-bauer · 2016-11-27T15:58:27Z

I'd be happy if you merge if you think the coverage is ok.

Luthaf · 2016-11-27T18:02:09Z

Merged 🎉 !

Luthaf reviewed Nov 1, 2016

View reviewed changes

Luthaf requested changes Nov 15, 2016

View reviewed changes

Implementation of adjustments of MC moves ranges.

f5d8e52

Luthaf force-pushed the adjust_displacement branch from 773c8d1 to f5d8e52 Compare November 25, 2016 10:51

Luthaf approved these changes Nov 25, 2016

View reviewed changes

Added missing implementation of update_amplitude in testcase of monte…

3cc219c

…_carlo.rs

Moved error handling for invalid acceptances to set_acceptance. Added…

a951e22

… tests.

Luthaf merged commit b6db02e into lumol-org:master Nov 27, 2016

g-bauer deleted the adjust_displacement branch November 28, 2016 15:43

Implement feature: adjustable maximum displacement in Monte-Carlo #52

Implement feature: adjustable maximum displacement in Monte-Carlo #52

Conversation

g-bauer commented Nov 1, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Luthaf commented Nov 1, 2016

g-bauer commented Nov 1, 2016

Luthaf commented Nov 1, 2016

g-bauer commented Nov 15, 2016 • edited Loading

Luthaf left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Luthaf Nov 15, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Luthaf commented Nov 15, 2016

Luthaf commented Nov 15, 2016

g-bauer commented Nov 25, 2016

Luthaf commented Nov 25, 2016

Luthaf commented Nov 26, 2016

Luthaf commented Nov 27, 2016

g-bauer commented Nov 27, 2016

Luthaf commented Nov 27, 2016

g-bauer commented Nov 15, 2016 •

edited

Loading

Luthaf Nov 15, 2016 •

edited

Loading