Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

User provided proposal distribution for insert and remove moves #176

Open
wants to merge 1 commit into
base: unstable
Choose a base branch
from
Open

User provided proposal distribution for insert and remove moves #176

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 69 additions & 8 deletions c++/triqs_cthyb/moves/insert.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,21 +23,30 @@

namespace triqs_cthyb {

histogram *move_insert_c_cdag::add_histo(std::string const &name, histo_map_t *histos) {
histogram *move_insert_c_cdag::add_histo(std::string const &name, histo_map_t *histos, int nbins) {
if (!histos) return nullptr;
auto new_histo = histos->insert({name, {.0, config.beta(), 100}});
auto new_histo = histos->insert({name, {.0, config.beta(), nbins}});
return &(new_histo.first->second);
}

move_insert_c_cdag::move_insert_c_cdag(int block_index, int block_size, std::string const &block_name, qmc_data &data,
mc_tools::random_generator &rng, histo_map_t *histos)
mc_tools::random_generator &rng, histo_map_t *histos, int nbins,
std::vector<double> const *hist_insert, std::vector<double> const *hist_remove,
std::vector<time_pt> const *taus_bin, bool use_improved_sampling)
: data(data),
config(data.config),
rng(rng),
block_index(block_index),
block_size(block_size),
histo_proposed(add_histo("insert_length_proposed_" + block_name, histos)),
histo_accepted(add_histo("insert_length_accepted_" + block_name, histos)) {}
histo_proposed(add_histo("insert_length_proposed_" + block_name, histos, nbins)),
histo_accepted(add_histo("insert_length_accepted_" + block_name, histos, nbins)),
hist_insert(hist_insert ? hist_insert : nullptr),
hist_remove(hist_remove ? hist_remove : nullptr),
step_i(time_pt::Nmax / (nbins - 1)),
step_d(config.beta() / double(nbins - 1)),
use_improved_sampling(use_improved_sampling),
taus_bin(taus_bin ? taus_bin : nullptr),
t1(time_pt(1, config.beta())) {}

mc_weight_t move_insert_c_cdag::attempt() {

Expand All @@ -52,9 +61,62 @@ namespace triqs_cthyb {
op1 = op_desc{block_index, rs1, true, data.linindex[std::make_pair(block_index, rs1)]};
op2 = op_desc{block_index, rs2, false, data.linindex[std::make_pair(block_index, rs2)]};

auto &det = data.dets[block_index];
int det_size = det.size();

// Choice of times for insertion. Find the time as double and them put them on the grid.
tau1 = data.tau_seg.get_random_pt(rng);
tau2 = data.tau_seg.get_random_pt(rng);
double fac = 1.;
if (use_improved_sampling) {
// first choose the bin, each bin being weighted by the probability hist_insert[bin]*length(bin)
double ran = double(rng(time_pt::Nmax)) / double(time_pt::Nmax - 1); // random number in [0,1]
double csum = 0;
int nbins = (*hist_insert).size();
int ibin1 = 0, ibin2 = 0;
for (int i = 0; i < nbins; ++i) {
double step = step_d;
if (i==0 || i==nbins-1) step /= 2.; // first and last bins are half-sized
csum += (*hist_insert)[i] * step;
if (csum >= ran || i == (nbins - 1) ) {
ibin1 = i;
break;
}
}
// now draw a time point uniformly within this bin
tau2 = tau1 + data.tau_seg.get_random_pt(rng, (*taus_bin)[ibin1], (*taus_bin)[ibin1+1]);

// compute the probability of proposing the current config from the trial one
// this is simply hist_remove[bin(tau1-tau2)] / sum_i(hist_remove(bin(tau_i - tau2)))
// where the sum is performed over all creation operators of the current block, including the trial one
// (careful, the convention for dtau in the remove move is t_cdag - t_c, while it is t_c - t_cdag in the insert move)
time_pt dtau_r = tau1 - tau2;
// find the bin of tau1 - tau2 (different from ibin, which is the bin of tau2 - tau1)
if (dtau_r < (*taus_bin)[1]) // special treatment for first and last bin since they're not the same size
ibin2 = 0;
else if (dtau_r >= (*taus_bin)[nbins-1])
ibin2 = nbins-1;
else
ibin2 = (floor_div(dtau_r, t1) - step_i / 2 ) / step_i + 1;
int ind;
double s = (*hist_remove)[ibin2]; // normalization constant = sum_i(hist_remove(bin(tau_i - tau2)))

for (int i = 0; i < det_size; ++i) {
dtau_r = det.get_x(i).first - tau2;
if (dtau_r < (*taus_bin)[1])
ind = 0;
else if (dtau_r >= (*taus_bin)[nbins-1])
ind = nbins-1;
else
ind = (floor_div(dtau_r, t1) - step_i / 2) / step_i + 1;
s += (*hist_remove)[ind];
}

// corrective factor for t_ratio
if (std::abs(s) < 1.e-15 || (*hist_remove)[ibin2] == 0.0) return 0; // quick return
fac = double(data.dets[block_index].size() + 1) * (*hist_remove)[ibin2] / (s * config.beta() * (*hist_insert)[ibin1]);
}
else
tau2 = data.tau_seg.get_random_pt(rng);

#ifdef EXT_DEBUG
std::cerr << "* Proposing to insert:" << std::endl;
Expand All @@ -80,8 +142,6 @@ namespace triqs_cthyb {
}

// Computation of det ratio
auto &det = data.dets[block_index];
int det_size = det.size();

// Find the position for insertion in the determinant
// NB : the determinant stores the C in decreasing time order.
Expand All @@ -98,6 +158,7 @@ namespace triqs_cthyb {

// proposition probability
mc_weight_t t_ratio = std::pow(block_size * config.beta() / double(det.size() + 1), 2);
if (use_improved_sampling) t_ratio *= fac;

// For quick abandon
double random_number = rng.preview();
Expand Down
12 changes: 10 additions & 2 deletions c++/triqs_cthyb/moves/insert.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,12 +36,20 @@ namespace triqs_cthyb {
h_scalar_t new_atomic_weight, new_atomic_reweighting;
time_pt tau1, tau2;
op_desc op1, op2;
const std::vector<double> *hist_insert;
const std::vector<double> *hist_remove;
const uint64_t step_i;
const double step_d;
bool use_improved_sampling;
const std::vector<time_pt> *taus_bin;
const time_pt t1;

histogram *add_histo(std::string const &name, histo_map_t *histos);
histogram *add_histo(std::string const &name, histo_map_t *histos, int nbins);

public:
move_insert_c_cdag(int block_index, int block_size, std::string const &block_name, qmc_data &data, mc_tools::random_generator &rng,
histo_map_t *histos);
histo_map_t *histos, int nbins, std::vector<double> const *hist_insert,
std::vector<double> const *hist_remove, std::vector<time_pt> const *taus_bin, bool use_improved_sampling);

mc_weight_t attempt();
mc_weight_t accept();
Expand Down
95 changes: 81 additions & 14 deletions c++/triqs_cthyb/moves/remove.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,21 +23,31 @@

namespace triqs_cthyb {

histogram * move_remove_c_cdag::add_histo(std::string const &name, histo_map_t *histos) {
histogram * move_remove_c_cdag::add_histo(std::string const &name, histo_map_t *histos, int nbins) {
if (!histos) return nullptr;
auto new_histo = histos->insert({name, {.0, config.beta(), 100}});
auto new_histo = histos->insert({name, {.0, config.beta(), nbins}});
return &(new_histo.first->second);
}

move_remove_c_cdag::move_remove_c_cdag(int block_index, int block_size, std::string const &block_name, qmc_data &data, mc_tools::random_generator &rng,
histo_map_t *histos)
histo_map_t *histos, int nbins, std::vector<double> const *hist_insert,
std::vector<double> const *hist_remove, std::vector<time_pt> const *taus_bin, bool use_improved_sampling)
: data(data),
config(data.config),
rng(rng),
block_index(block_index),
block_size(block_size),
histo_proposed(add_histo("remove_length_proposed_" + block_name, histos)),
histo_accepted(add_histo("remove_length_accepted_" + block_name, histos)) {}
histo_proposed(add_histo("remove_length_proposed_" + block_name, histos, nbins)),
histo_accepted(add_histo("remove_length_accepted_" + block_name, histos, nbins)),
hist_insert(hist_insert ? hist_insert : nullptr),
hist_remove(hist_remove ? hist_remove : nullptr),
step_i(time_pt::Nmax / (nbins - 1)),
use_improved_sampling(use_improved_sampling),
taus_bin(taus_bin ? taus_bin : nullptr),
t1(time_pt(1, config.beta())),
Nmax(100) {
bins.reserve(Nmax);
}

mc_weight_t move_remove_c_cdag::attempt() {

Expand All @@ -63,21 +73,78 @@ namespace triqs_cthyb {

// now mark 2 nodes for deletion
tau1 = data.imp_trace.try_delete(num_c, block_index, false);
tau2 = data.imp_trace.try_delete(num_c_dag, block_index, true);

double fac = 1.;
if (use_improved_sampling) {
// choose the creation operator to remove, weighted by probability
// hist_remove[bin(tau2-tau1)] / sum_i(hist_remove(bin(tau_i - tau1)))
double s = 0; // normalization constant sum_i(hist_remove(bin(tau_i - tau1)))
int nbins = (*hist_insert).size();
if (det_size > Nmax) {
Nmax *= 2;
bins.reserve(Nmax);
}
bins.resize(det_size);
for (int i = 0; i < det_size; ++i) {
int ind;
// need to call try_delete to get binary tree index
time_pt dtau_r = data.imp_trace.try_delete(i, block_index, true) - tau1;
if (dtau_r < (*taus_bin)[1])
ind = 0;
else if (dtau_r >= (*taus_bin)[nbins-1])
ind = nbins-1;
else
ind = (floor_div(dtau_r, t1) - step_i / 2) / step_i + 1;
bins[i] = ind;
s += (*hist_remove)[ind];
}

data.imp_trace.cancel_delete();
if (std::abs(s) <= 1.e-15) return 0; // quick return

// draw a uniform variable on [0,1]
double ran = double(rng(time_pt::Nmax)) / double(time_pt::Nmax - 1);
// choose the creation operator
double csum = 0;
for (int i = 0; i < det_size; ++i) {
csum += (*hist_remove)[bins[i]] / s;
if (csum >= ran || i == (nbins-1)) {
num_c_dag = i;
break;
}
}

tau1 = data.imp_trace.try_delete(num_c, block_index, false);
tau2 = data.imp_trace.try_delete(num_c_dag, block_index, true);

// compute the probability of proposing the current config from the trial one
// this is simply hist_insert[bin(tau1-tau2)] (since for the insert move, dtau=t_c - t_cdag)
time_pt dtau_i = tau1 - tau2;
int ibin;
if (dtau_i < (*taus_bin)[1])
ibin = 0;
else if (dtau_i >= (*taus_bin)[nbins-1])
ibin = nbins-1;
else
ibin = (floor_div(dtau_i, t1) - step_i / 2) / step_i + 1;
if ((*hist_insert)[ibin] == 0.0) return 0; // quick return
fac = s * config.beta() * (*hist_insert)[ibin] / (double(det_size) * (*hist_remove)[bins[num_c_dag]]);
}
else
tau2 = data.imp_trace.try_delete(num_c_dag, block_index, true);
// record the length of the proposed removal
dtau = double(tau2 - tau1);
if (histo_proposed) *histo_proposed << dtau;

auto det_ratio = det.try_remove(num_c_dag, num_c);

// proposition probability
auto t_ratio = std::pow(block_size * config.beta() / double(det_size), 2); // Size of the det before the try_delete!
auto t_ratio = std::pow(double(det_size) / (block_size * config.beta()), 2); // Size of the det before the try_delete!
if (use_improved_sampling) t_ratio *= fac;

// For quick abandon
double random_number = rng.preview();
if (random_number == 0.0) return 0;
double p_yee = std::abs(det_ratio / t_ratio / data.atomic_weight);
double p_yee = std::abs(det_ratio * t_ratio / data.atomic_weight);

// recompute the atomic_weight
std::tie(new_atomic_weight, new_atomic_reweighting) = data.imp_trace.compute(p_yee, random_number);
Expand All @@ -98,22 +165,22 @@ namespace triqs_cthyb {
std::cerr << "Trace ratio: " << atomic_weight_ratio << '\t';
std::cerr << "Det ratio: " << det_ratio << '\t';
std::cerr << "Prefactor: " << t_ratio << '\t';
std::cerr << "Weight: " << p / t_ratio << std::endl;
std::cerr << "Weight: " << p * t_ratio << std::endl;
#endif

if (!isfinite(p)) {
std::cerr << "Remove move info\n";
std::cerr << "Trace ratio: " << atomic_weight_ratio << '\t';
std::cerr << "Det ratio: " << det_ratio << '\t';
std::cerr << "Prefactor: " << t_ratio << '\t';
std::cerr << "Weight: " << p / t_ratio << std::endl;
std::cerr << "Weight: " << p * t_ratio << std::endl;
TRIQS_RUNTIME_ERROR << "(remove) p not finite :" << p << " in config " << config.get_id();
}

if (!isfinite(p / t_ratio)){
TRIQS_RUNTIME_ERROR << "(remove) p / t_ratio not finite p : " << p << " t_ratio : " << t_ratio << " in config " << config.get_id();
if (!isfinite(p * t_ratio)){
TRIQS_RUNTIME_ERROR << "(remove) p * t_ratio not finite p : " << p << " t_ratio : " << t_ratio << " in config " << config.get_id();
}
return p / t_ratio;
return p * t_ratio;
}

mc_weight_t move_remove_c_cdag::accept() {
Expand Down
13 changes: 11 additions & 2 deletions c++/triqs_cthyb/moves/remove.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,12 +36,21 @@ namespace triqs_cthyb {
double dtau;
h_scalar_t new_atomic_weight, new_atomic_reweighting;
time_pt tau1, tau2;
const std::vector<double> *hist_insert;
const std::vector<double> *hist_remove;
const uint64_t step_i;
bool use_improved_sampling;
const std::vector<time_pt> *taus_bin;
const time_pt t1;
int Nmax;
std::vector<int> bins;

histogram *add_histo(std::string const &name, histo_map_t *histos);
histogram *add_histo(std::string const &name, histo_map_t *histos, int nbins);

public:
move_remove_c_cdag(int block_index, int block_size, std::string const &block_name, qmc_data &data, mc_tools::random_generator &rng,
histo_map_t *histos);
histo_map_t *histos, int nbins, std::vector<double> const *hist_insert,
std::vector<double> const *hist_remove, std::vector<time_pt> const *taus_bin, bool use_improved_sampling);

mc_weight_t attempt();
mc_weight_t accept();
Expand Down
8 changes: 8 additions & 0 deletions c++/triqs_cthyb/parameters.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,10 @@ namespace triqs_cthyb {
h5_write(grp, "performance_analysis", sp.performance_analysis);
h5_write(grp, "proposal_prob", sp.proposal_prob);

h5_write(grp, "nbins_histo", sp.nbins_histo);
h5_write(grp, "hist_insert", sp.hist_insert);
h5_write(grp, "hist_remove", sp.hist_remove);

//h5_write(grp, "move_global", sp.move_global);
if( sp.move_global.size() != 0 )
TRIQS_RUNTIME_ERROR << "Error serailizing: CTHYB solve_parameters, can not serialize the global moves data type.";
Expand Down Expand Up @@ -197,6 +201,10 @@ namespace triqs_cthyb {
h5_read(grp, "performance_analysis", sp.performance_analysis);
h5_read(grp, "proposal_prob", sp.proposal_prob);

h5_read(grp, "nbins_histo", sp.nbins_histo);
h5_read(grp, "hist_insert", sp.hist_insert);
h5_read(grp, "hist_remove", sp.hist_remove);

//h5_read(grp, "move_global", sp.move_global);
if( grp.has_key("move_global") )
TRIQS_RUNTIME_ERROR << "Error reading: CTHYB solve_parameters, can not de-serialize the global moves data type.";
Expand Down
9 changes: 9 additions & 0 deletions c++/triqs_cthyb/parameters.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,6 +201,15 @@ namespace triqs_cthyb {
/// default: {}
std::map<std::string, double> proposal_prob = {};

/// Number of bins for the histograms
int nbins_histo = 100;

/// Proposal distribution for insert move
std::map<std::string, std::vector<double>> hist_insert = {};

/// Proposal distribution for remove move
std::map<std::string, std::vector<double>> hist_remove = {};

/// List of global moves (with their names).
/// Each move is specified with an index substitution dictionary.
/// type: dict(str : dict(indices : indices))
Expand Down
Loading
Loading