use crate::sse::*;
use crate::util::allocator::Factory;
use crate::util::bondcontainer::BondContainer;
use rand::Rng;
pub trait EdgeNavigator {
fn n_bonds(&self) -> usize;
fn bonds_for_var(&self, var: usize) -> &[usize];
fn vars_for_bond(&self, bond: usize) -> (usize, usize);
fn bond_prefers_aligned(&self, bond: usize) -> bool;
fn other_var_for_bond(&self, var: usize, bond: usize) -> Option<usize> {
let (a, b) = self.vars_for_bond(bond);
if var == a {
Some(b)
} else if var == b {
Some(a)
} else {
None
}
}
}
pub trait RVBUpdater:
DiagonalSubsection
+ Factory<Vec<usize>>
+ Factory<Vec<bool>>
+ Factory<BondContainer<usize>>
+ Factory<Vec<Option<usize>>>
{
fn constant_ops_on_var(&self, var: usize, ps: &mut Vec<usize>);
fn spin_flips_on_var(&self, var: usize, ps: &mut Vec<usize>);
fn rvb_update<R: Rng, EN: EdgeNavigator, H>(
&mut self,
edges: &EN,
state: &mut [bool],
updates: usize,
diagonal_edge_hamiltonian: H,
rng: &mut R,
) -> usize
where
H: Fn(usize, bool, bool) -> f64,
{
self.rvb_update_with_ising_weight(
edges,
state,
updates,
diagonal_edge_hamiltonian,
|_| 1.0,
rng,
)
}
fn rvb_update_with_ising_weight<R: Rng, EN: EdgeNavigator, H, F>(
&mut self,
edges: &EN,
state: &mut [bool],
updates: usize,
diagonal_edge_hamiltonian: H,
ising_ratio: F,
rng: &mut R,
) -> usize
where
H: Fn(usize, bool, bool) -> f64,
F: Fn(&Self::Op) -> f64,
{
let mut var_starts: Vec<usize> = self.get_instance();
let mut var_lengths: Vec<usize> = self.get_instance();
let mut constant_ps: Vec<usize> = self.get_instance();
let mut vars_with_zero_ops: Vec<usize> = self.get_instance();
find_constants(
self,
&mut var_starts,
&mut var_lengths,
&mut constant_ps,
&mut vars_with_zero_ops,
);
let mut num_succ = 0;
for _ in 0..updates {
let choice = rng.gen_range(0, constant_ps.len() + vars_with_zero_ops.len());
let (v, flip) = if choice < constant_ps.len() {
let res = var_starts.binary_search(&choice);
let v = match res {
Err(i) => i - 1,
Ok(i) => i,
};
(v, Some(choice))
} else {
let choice = choice - constant_ps.len();
(vars_with_zero_ops[choice], None)
};
let mut cluster_vars: Vec<usize> = self.get_instance();
let mut cluster_flips: Vec<Option<usize>> = self.get_instance();
let mut boundary_vars: Vec<usize> = self.get_instance();
let mut boundary_flips_pos: Vec<Option<usize>> = self.get_instance();
let cluster_size = contiguous_bits(rng) + 1;
build_cluster(
cluster_size,
(v, flip),
(self.get_nvars(), self.get_cutoff()),
(&mut cluster_vars, &mut cluster_flips),
(&mut boundary_vars, &mut boundary_flips_pos),
(&var_starts, &var_lengths),
&constant_ps,
edges,
self,
rng,
);
let mut cluster_starting_state: Vec<bool> = self.get_instance();
let mut cluster_toggle_ps: Vec<usize> = self.get_instance();
let mut subvars: Vec<usize> = self.get_instance();
let mut var_to_subvar: Vec<Option<usize>> = self.get_instance();
var_to_subvar.resize(self.get_nvars(), None);
subvars.extend(cluster_vars.iter().chain(boundary_vars.iter()).cloned());
subvars.sort_unstable();
subvars.dedup();
subvars
.iter()
.cloned()
.enumerate()
.for_each(|(subvar, var)| var_to_subvar[var] = Some(subvar));
cluster_starting_state.resize(subvars.len(), false);
cluster_vars
.iter()
.cloned()
.zip(cluster_flips.iter().cloned())
.for_each(|(v, fi)| {
let subvar = var_to_subvar[v].unwrap();
if let Some(fi) = fi {
let vstart = var_starts[v];
let fi_rel = fi - vstart;
if fi_rel + 1 >= var_lengths[v] {
cluster_starting_state[subvar] = true;
cluster_toggle_ps.push(constant_ps[fi]);
cluster_toggle_ps.push(constant_ps[vstart]);
} else {
cluster_toggle_ps.push(constant_ps[fi]);
cluster_toggle_ps.push(constant_ps[fi + 1]);
}
} else {
cluster_starting_state[subvar] = true;
}
});
self.return_instance(cluster_flips);
self.return_instance(cluster_vars);
let mut substate: Vec<bool> = self.get_instance();
substate.extend(subvars.iter().cloned().map(|v| state[v]));
let mut subvar_boundary_tops: Vec<Option<usize>> = self.get_instance();
subvar_boundary_tops.resize(subvars.len(), None);
boundary_vars
.iter()
.zip(boundary_flips_pos.iter())
.for_each(|(bv, bfp)| {
let subvar =
var_to_subvar[*bv].expect("Boundary must be in var_to_subvar array.");
match (bfp, subvar_boundary_tops[subvar]) {
(Some(bfp), Some(bt)) => {
if *bfp < bt {
subvar_boundary_tops[subvar] = Some(constant_ps[*bfp])
}
}
(Some(_), None) | (None, None) => {
subvar_boundary_tops[subvar] = bfp.map(|bfp| constant_ps[bfp])
}
(None, Some(_)) => unreachable!(),
};
});
self.return_instance(boundary_vars);
self.return_instance(boundary_flips_pos);
cluster_toggle_ps.sort_unstable();
remove_doubles(&mut cluster_toggle_ps);
let p_to_flip = calculate_flip_prob(
self,
(state, &mut substate),
(&mut cluster_starting_state, &cluster_toggle_ps),
&subvar_boundary_tops,
(&subvars, |v| var_to_subvar[v]),
(&diagonal_edge_hamiltonian, &ising_ratio, edges),
);
let should_mutate = if p_to_flip >= 1.0 {
true
} else {
if p_to_flip < 0. {
println!("P = {}", p_to_flip);
debug_print_diagonal(self, state);
}
rng.gen_bool(p_to_flip)
};
if should_mutate {
mutate_graph(
self,
(state, &mut substate),
(&mut cluster_starting_state, &cluster_toggle_ps),
&subvar_boundary_tops,
(&subvars, |v| var_to_subvar[v]),
(&diagonal_edge_hamiltonian, edges),
rng,
);
let starting_cluster = cluster_starting_state
.iter()
.cloned()
.filter(|b| *b)
.count()
> 0;
if starting_cluster {
subvars
.iter()
.cloned()
.zip(cluster_starting_state.iter().cloned())
.for_each(|(v, c)| {
state[v] = state[v] != c;
});
}
num_succ += 1;
}
self.return_instance(var_to_subvar);
self.return_instance(subvar_boundary_tops);
self.return_instance(substate);
self.return_instance(subvars);
self.return_instance(cluster_toggle_ps);
self.return_instance(cluster_starting_state);
}
self.return_instance(vars_with_zero_ops);
self.return_instance(constant_ps);
self.return_instance(var_lengths);
self.return_instance(var_starts);
num_succ
}
}
fn mutate_graph<RVB: RVBUpdater + ?Sized, VS, EN: EdgeNavigator + ?Sized, R: Rng, H>(
rvb: &mut RVB,
(state, substate): (&[bool], &mut [bool]),
(cluster_state, cluster_flips): (&mut [bool], &[usize]),
boundary_tops: &[Option<usize>],
(vars, var_to_subvar): (&[usize], VS),
(diagonal_hamiltonian, edges): (H, &EN),
rng: &mut R,
) where
VS: Fn(usize) -> Option<usize>,
H: Fn(usize, bool, bool) -> f64,
{
let mut jump_to: Vec<usize> = rvb.get_instance();
let mut continue_until: Vec<usize> = rvb.get_instance();
let mut count = cluster_state.iter().filter(|b| **b).count();
let has_starting_cluster = if count != 0 {
jump_to.push(0);
substate
.iter_mut()
.zip(cluster_state.iter().cloned())
.for_each(|(b, c)| *b = *b != c);
true
} else {
false
};
cluster_flips.iter().cloned().for_each(|p| {
if count == 0 {
jump_to.push(p)
}
let op = rvb.get_node_ref(p).unwrap().get_op_ref();
count = op
.get_vars()
.iter()
.cloned()
.filter_map(|v| var_to_subvar(v))
.fold(count, |count, subvar| {
cluster_state[subvar] = !cluster_state[subvar];
if cluster_state[subvar] {
count + 1
} else {
count - 1
}
});
if count == 0 {
continue_until.push(p)
}
});
if count != 0 {
debug_assert!(has_starting_cluster);
continue_until.push(rvb.get_cutoff());
}
debug_assert_eq!(
jump_to.len(),
continue_until.len(),
"Should have the same number of starts and ends."
);
let mut bonds: BondContainer<usize> = rvb.get_instance();
let mut next_cluster_index = 0;
vars.iter()
.cloned()
.filter(|v| cluster_state[var_to_subvar(*v).unwrap()])
.for_each(|v| {
edges.bonds_for_var(v).iter().cloned().for_each(|b| {
let ov = edges.other_var_for_bond(v, b).unwrap();
if !cluster_state[var_to_subvar(ov).unwrap()] {
let (va, vb) = edges.vars_for_bond(b);
let subva = var_to_subvar(va).unwrap();
let subvb = var_to_subvar(vb).unwrap();
let w = diagonal_hamiltonian(b, substate[subva], substate[subvb]);
bonds.insert(b, w);
}
})
});
jump_to
.iter()
.cloned()
.zip(continue_until.iter().cloned())
.fold(
(substate, cluster_state, rng),
|(substate, cluster_state, rng), (from, until)| {
rvb.get_propagated_substate_with_hint(
from,
substate,
state,
vars,
boundary_tops.iter().cloned(),
);
substate
.iter_mut()
.zip(cluster_state.iter().cloned())
.for_each(|(b, c)| *b = *b != c);
let mut args = rvb.get_empty_args(SubvarAccess::VARLIST(&vars));
rvb.fill_args_at_p_with_hint(from, &mut args, vars, boundary_tops.iter().cloned());
let acc = (next_cluster_index, &mut bonds, substate, cluster_state, rng);
let ret = rvb.mutate_subsection_ops(
from,
until,
acc,
|_, op, p, acc| {
let (mut next_cluster_index, bonds, substate, cluster_state, mut rng) = acc;
let in_bonds = bonds.contains(&op.get_bond());
let at_next_cluster_flip = next_cluster_index < cluster_flips.len()
&& p == cluster_flips[next_cluster_index];
let newop = if in_bonds {
debug_assert!(op.is_diagonal());
let new_bond = bonds.get_random(&mut rng).unwrap().0;
let (new_a, new_b) = edges.vars_for_bond(new_bond);
let vars = RVB::Op::make_vars([new_a, new_b].iter().cloned());
let state = RVB::Op::make_substate(
[new_a, new_b]
.iter()
.cloned()
.map(|v| var_to_subvar(v).unwrap())
.map(|subvar| substate[subvar]),
);
let new_op = RVB::Op::diagonal(vars, new_bond, state, op.is_constant());
Some(Some(new_op))
} else {
let new_op = if at_next_cluster_flip {
debug_assert!(op.is_constant());
debug_assert!(op
.get_vars()
.iter()
.cloned()
.all(|v| var_to_subvar(v).is_some()));
let new_op = op.clone_and_edit_in_out(|ins, outs| {
op.get_vars()
.iter()
.cloned()
.map(|v| var_to_subvar(v).unwrap())
.zip(ins.iter_mut().zip(outs.iter_mut()))
.for_each(|(subvar, (bin, bout))| {
*bin = *bin != cluster_state[subvar];
*bout = *bout != !cluster_state[subvar];
});
});
new_op
.get_vars()
.iter()
.cloned()
.map(|v| var_to_subvar(v).unwrap())
.zip(new_op.get_outputs().iter().cloned())
.for_each(|(subvar, bout)| {
cluster_state[subvar] = !cluster_state[subvar];
substate[subvar] = bout
});
next_cluster_index += 1;
Some(Some(new_op))
} else {
debug_assert!({
let all_in = op.get_vars().iter().cloned().all(|v| {
var_to_subvar(v)
.map(|subvar| cluster_state[subvar])
.unwrap_or(false)
});
let all_out = op.get_vars().iter().cloned().all(|v| {
var_to_subvar(v)
.map(|subvar| !cluster_state[subvar])
.unwrap_or(true)
});
let succ = (all_in != all_out) && (all_in || all_out);
if !succ {
println!("subvars: {:?}", vars);
println!(
"op: {:?}\t{:?} -> {:?}",
op.get_vars(),
op.get_inputs(),
op.get_outputs()
);
println!("all_in: {}\tall_out: {}", all_in, all_out);
}
succ
});
let any_subvars = op
.get_vars()
.iter()
.cloned()
.any(|v| var_to_subvar(v).is_some());
let any_in_cluster = op
.get_vars()
.iter()
.cloned()
.filter_map(&var_to_subvar)
.any(|subvar| cluster_state[subvar]);
if !any_subvars || (!any_in_cluster && op.is_diagonal()) {
None
} else {
let new_op = if any_in_cluster {
debug_assert!(op
.get_vars()
.iter()
.cloned()
.all(|v| var_to_subvar(v).is_some()));
let new_op = op.clone_and_edit_in_out_symmetric(|state| {
state.iter_mut().for_each(|b| *b = !*b);
});
if !new_op.is_diagonal() {
new_op
.get_vars()
.iter()
.cloned()
.map(|v| var_to_subvar(v).unwrap())
.zip(new_op.get_outputs().iter().cloned())
.for_each(|(subvar, bout)| {
substate[subvar] = bout;
});
}
Some(Some(new_op))
} else {
if !op.is_diagonal() {
op.get_vars()
.iter()
.cloned()
.filter_map(|v| var_to_subvar(v))
.zip(op.get_outputs().iter().cloned())
.for_each(|(subvar, bout)| {
substate[subvar] = bout;
});
}
None
};
new_op
}
};
op.get_vars()
.iter()
.cloned()
.filter_map(|v| var_to_subvar(v).map(|subvar| (v, subvar)))
.for_each(|(v, subvar)| {
edges
.bonds_for_var(v)
.iter()
.cloned()
.filter_map(|b| {
let ov = edges.other_var_for_bond(v, b).unwrap();
var_to_subvar(ov).map(|o_subvar| (b, subvar, o_subvar))
})
.for_each(|(b, subvar, o_subvar)| {
if cluster_state[subvar] == cluster_state[o_subvar] {
if bonds.contains(&b) {
bonds.remove(&b);
}
} else {
let (va, vb) = edges.vars_for_bond(b);
let subva = var_to_subvar(va).unwrap();
let subvb = var_to_subvar(vb).unwrap();
let w = diagonal_hamiltonian(
b,
substate[subva],
substate[subvb],
);
bonds.insert(b, w);
}
})
});
new_op
};
(
newop,
(next_cluster_index, bonds, substate, cluster_state, rng),
)
},
Some(args),
);
next_cluster_index = ret.0;
let substate = ret.2;
let cluster_state = ret.3;
let rng = ret.4;
(substate, cluster_state, rng)
},
);
rvb.return_instance(bonds);
rvb.return_instance(jump_to);
rvb.return_instance(continue_until);
}
fn calculate_flip_prob<RVB: RVBUpdater + ?Sized, VS, EN: EdgeNavigator + ?Sized, H, F>(
rvb: &mut RVB,
(state, substate): (&[bool], &mut [bool]),
(cluster_state, cluster_flips): (&mut [bool], &[usize]),
boundary_tops: &[Option<usize>],
(vars, var_to_subvar): (&[usize], VS),
(diagonal_hamiltonian, ising_ratio, edges): (H, F, &EN),
) -> f64
where
VS: Fn(usize) -> Option<usize>,
H: Fn(usize, bool, bool) -> f64,
F: Fn(&RVB::Op) -> f64,
{
let mut cluster_size = cluster_state.iter().cloned().filter(|x| *x).count();
let mut psel = rvb.get_first_p();
let mut next_cluster_index = 0;
let mut mult = 1.0;
let ws_for_flip = |b: usize, subvar_to_flip: usize, substate: &[bool]| {
let (va, vb) = edges.vars_for_bond(b);
let suba = var_to_subvar(va).unwrap();
let subb = var_to_subvar(vb).unwrap();
debug_assert!(subvar_to_flip == suba || subvar_to_flip == subb);
let ba = substate[suba];
let bb = substate[subb];
let w_before = diagonal_hamiltonian(b, ba, bb);
let (ba, bb) = if subvar_to_flip == suba {
(!ba, bb)
} else {
(ba, !bb)
};
let w_after = diagonal_hamiltonian(b, ba, bb);
(w_before, w_after)
};
let mut bonds_before: BondContainer<usize> = rvb.get_instance();
let mut bonds_after: BondContainer<usize> = rvb.get_instance();
let mut n_bonds = 0;
if cluster_size != 0 {
vars.iter()
.cloned()
.map(|v| (v, var_to_subvar(v).unwrap()))
.filter(|(_, subvar)| cluster_state[*subvar])
.for_each(|(v, subvar)| {
edges.bonds_for_var(v).iter().cloned().for_each(|b| {
let ov = edges.other_var_for_bond(v, b).unwrap();
let o_subvar = var_to_subvar(ov).unwrap();
if !cluster_state[o_subvar] {
let (wbef, waft) = ws_for_flip(b, subvar, substate);
bonds_before.insert(b, wbef);
bonds_after.insert(b, waft);
}
})
});
}
while let Some(mut p) = psel {
if cluster_size == 0 {
debug_assert_eq!(bonds_before.len(), 0);
debug_assert_eq!(bonds_after.len(), 0);
debug_assert_eq!(n_bonds, 0);
if next_cluster_index < cluster_flips.len() {
p = cluster_flips[next_cluster_index];
rvb.get_propagated_substate_with_hint(
p,
substate,
state,
vars,
boundary_tops.iter().cloned(),
);
} else {
break;
}
}
let node = rvb.get_node_ref(p).unwrap();
let op = node.get_op_ref();
let near_cluster = op
.get_vars()
.iter()
.cloned()
.any(|v| var_to_subvar(v).is_some());
if near_cluster {
let is_cluster_bound =
next_cluster_index < cluster_flips.len() && p == cluster_flips[next_cluster_index];
let will_flip_spins = !op.is_diagonal();
let will_change_bonds = will_flip_spins || is_cluster_bound;
let completely_in_cluster = op.get_vars().iter().cloned().all(|v| {
var_to_subvar(v)
.map(|subvar| cluster_state[subvar])
.unwrap_or(false)
});
let b = op.get_bond();
if bonds_before.contains(&b) {
debug_assert!(!is_cluster_bound);
debug_assert!(!will_flip_spins);
debug_assert!(bonds_after.contains(&b));
n_bonds += 1;
} else {
if is_cluster_bound {
debug_assert!(op.is_constant());
debug_assert_eq!(op.get_vars().len(), 1);
let v = op.get_vars()[0];
let subvar = var_to_subvar(v).unwrap();
cluster_state[subvar] = !cluster_state[subvar];
if cluster_state[subvar] {
cluster_size += 1
} else {
cluster_size -= 1
};
next_cluster_index += 1;
}
if will_flip_spins {
op.get_vars()
.iter()
.cloned()
.zip(op.get_outputs().iter().cloned())
.filter_map(|(v, bout)| {
var_to_subvar(v)
.zip(Some(bout))
.map(|(subvar, bout)| (subvar, bout))
})
.for_each(|(subvar, bout)| {
substate[subvar] = bout;
});
}
if completely_in_cluster {
let ising_flip_weight = ising_ratio(op);
mult *= ising_flip_weight;
if mult < std::f64::EPSILON {
break;
}
}
if will_change_bonds {
mult *= calculate_mult(&bonds_before, &bonds_after, n_bonds);
n_bonds = 0;
if mult < std::f64::EPSILON {
break;
}
op.get_vars()
.iter()
.cloned()
.filter_map(|v| var_to_subvar(v).map(|subvar| (v, subvar)))
.for_each(|(v, subvar)| {
edges
.bonds_for_var(v)
.iter()
.cloned()
.filter_map(|b| {
let ov = edges.other_var_for_bond(v, b).unwrap();
var_to_subvar(ov).map(|o_subvar| (b, subvar, o_subvar))
})
.for_each(|(b, subvar, o_subvar)| {
if cluster_state[subvar] == cluster_state[o_subvar] {
if bonds_before.contains(&b) {
bonds_before.remove(&b);
bonds_after.remove(&b);
}
} else {
let subvar = if cluster_state[subvar] {
subvar
} else {
o_subvar
};
let (wbef, waft) = ws_for_flip(b, subvar, substate);
bonds_before.insert(b, wbef);
bonds_after.insert(b, waft);
}
})
});
}
}
}
psel = rvb.get_next_p(node);
}
mult *= calculate_mult(&bonds_before, &bonds_after, n_bonds);
rvb.return_instance(bonds_before);
rvb.return_instance(bonds_after);
mult
}
fn remove_doubles<T: Eq + Copy>(v: &mut Vec<T>) {
let mut ii = 0;
let mut jj = 0;
while jj + 1 < v.len() {
if v[jj] == v[jj + 1] {
jj += 2;
} else {
v[ii] = v[jj];
ii += 1;
jj += 1;
}
}
if jj < v.len() {
v[ii] = v[jj];
ii += 1;
jj += 1;
}
while jj > ii {
v.pop();
jj -= 1;
}
}
fn build_cluster<Fact: ?Sized, EN, R: Rng>(
mut cluster_size: usize,
(init_var, init_flip): (usize, Option<usize>),
(nvars, cutoff): (usize, usize),
(cluster_vars, cluster_flips): (&mut Vec<usize>, &mut Vec<Option<usize>>),
(boundary_vars, boundary_flips_pos): (&mut Vec<usize>, &mut Vec<Option<usize>>),
(var_starts, var_lengths): (&[usize], &[usize]),
constant_ps: &[usize],
edges: &EN,
fact: &mut Fact,
rng: &mut R,
) where
Fact: Factory<Vec<bool>> + Factory<Vec<usize>> + Factory<Vec<Option<usize>>>,
EN: EdgeNavigator,
{
let mut empty_var_in_cluster: Vec<bool> = fact.get_instance();
let mut op_p_in_cluster: Vec<bool> = fact.get_instance();
empty_var_in_cluster.resize(nvars, false);
op_p_in_cluster.resize(cutoff, false);
boundary_vars.push(init_var);
boundary_flips_pos.push(init_flip);
if let Some(init_flip) = init_flip {
op_p_in_cluster[constant_ps[init_flip]] = true;
} else {
empty_var_in_cluster[init_var] = true;
}
while cluster_size > 0 && !boundary_vars.is_empty() {
let to_pop = rng.gen_range(0, boundary_vars.len());
let v = pop_index(boundary_vars, to_pop).unwrap();
let flip = pop_index(boundary_flips_pos, to_pop).unwrap();
debug_assert!(flip.map(|f| f >= var_starts[v]).unwrap_or(true));
debug_assert!(flip
.map(|f| f < var_starts[v] + var_lengths[v])
.unwrap_or(var_lengths[v] == 0));
cluster_vars.push(v);
cluster_flips.push(flip);
if let Some(flip) = flip {
let relflip = flip - var_starts[v];
let flip_dec = (relflip + var_lengths[v] - 1) % var_lengths[v] + var_starts[v];
let flip_inc = (relflip + 1) % var_lengths[v] + var_starts[v];
if !op_p_in_cluster[constant_ps[flip_dec]] {
op_p_in_cluster[constant_ps[flip_dec]] = true;
boundary_vars.push(v);
boundary_flips_pos.push(Some(flip_dec));
}
if !op_p_in_cluster[constant_ps[flip_inc]] {
op_p_in_cluster[constant_ps[flip_inc]] = true;
boundary_vars.push(v);
boundary_flips_pos.push(Some(flip_inc));
}
}
edges.bonds_for_var(v).iter().for_each(|b| {
let ov = edges.other_var_for_bond(v, *b).unwrap();
if var_lengths[ov] == 0 {
if !empty_var_in_cluster[ov] {
empty_var_in_cluster[ov] = true;
boundary_vars.push(ov);
boundary_flips_pos.push(None);
}
} else {
let pis = var_starts[ov]..var_starts[ov] + var_lengths[ov];
if let Some(flip) = flip {
let relflip = flip - var_starts[v];
let flip_inc = (relflip + 1) % var_lengths[v] + var_starts[v];
let pstart = constant_ps[flip];
let pend = constant_ps[flip_inc];
find_overlapping_starts(pstart, pend, cutoff, &constant_ps[pis])
.map(|i| i + var_starts[ov])
.for_each(|flip_pos| {
if !op_p_in_cluster[constant_ps[flip_pos]] {
op_p_in_cluster[constant_ps[flip_pos]] = true;
boundary_vars.push(ov);
boundary_flips_pos.push(Some(flip_pos));
}
})
} else {
pis.for_each(|pi| {
if !op_p_in_cluster[constant_ps[pi]] {
op_p_in_cluster[constant_ps[pi]] = true;
boundary_vars.push(ov);
boundary_flips_pos.push(Some(pi));
}
})
}
}
});
cluster_size -= 1;
}
fact.return_instance(empty_var_in_cluster);
fact.return_instance(op_p_in_cluster);
}
fn pop_index<T>(v: &mut Vec<T>, index: usize) -> Option<T> {
let len = v.len();
if index < len {
v.swap(index, len - 1);
v.pop()
} else {
None
}
}
fn find_overlapping_starts(
p_start: usize,
p_end: usize,
cutoff: usize,
flip_positions: &[usize],
) -> impl Iterator<Item = usize> + '_ {
debug_assert_ne!(flip_positions.len(), 0);
let bin_found = flip_positions.binary_search(&p_start).unwrap_err();
let prev_bin_found = (bin_found + flip_positions.len() - 1) % flip_positions.len();
let lowest_ps = flip_positions[prev_bin_found];
let offset_p_start = (p_start + cutoff - lowest_ps) % cutoff;
let offset_p_end = (p_end + cutoff - lowest_ps) % cutoff;
flip_positions[prev_bin_found..]
.iter()
.cloned()
.zip(prev_bin_found..)
.chain(
flip_positions[..prev_bin_found]
.iter()
.cloned()
.zip(0..prev_bin_found),
)
.take_while(move |(p, ip)| {
let check_start = (*p + cutoff - lowest_ps) % cutoff;
let next_p = flip_positions[(*ip + 1) % flip_positions.len()];
let check_end = (next_p + cutoff - lowest_ps) % cutoff;
let has_overlap_start = check_start < offset_p_start && offset_p_start < check_end;
let has_start_within = offset_p_start < check_start && check_start < offset_p_end;
let eq = (p_start == p_end) || (check_start == check_end);
let overlap = has_overlap_start || has_start_within;
eq || overlap
})
.map(|(_, ip)| ip)
}
fn find_constants<RVB>(
rvb: &RVB,
var_starts: &mut Vec<usize>,
var_lengths: &mut Vec<usize>,
constant_ps: &mut Vec<usize>,
vars_with_zero_ops: &mut Vec<usize>,
) where
RVB: RVBUpdater + ?Sized,
{
(0..rvb.get_nvars()).for_each(|v| {
let start = constant_ps.len();
var_starts.push(start);
rvb.constant_ops_on_var(v, constant_ps);
debug_assert!(
constant_ps.iter().cloned().all(|p| {
let op = rvb.get_node_ref(p).unwrap().get_op_ref();
op.get_vars().len() == 1
}),
"RVB cluster only supports constant ops with a single variable."
);
var_lengths.push(constant_ps.len() - start);
if constant_ps.len() == *var_starts.last().unwrap() {
vars_with_zero_ops.push(v);
}
});
}
pub fn contiguous_bits<R: Rng>(r: &mut R) -> usize {
let mut acc = 0;
let mut v = r.next_u64();
loop {
if v == std::u64::MAX {
acc += 64;
v = r.next_u64();
} else {
while v & 1 == 1 {
acc += 1;
v >>= 1;
}
break acc;
}
}
}
fn calculate_mult(
bonds_before: &BondContainer<usize>,
bonds_after: &BondContainer<usize>,
n: usize,
) -> f64 {
let close = (bonds_before.get_total_weight() - bonds_after.get_total_weight()).abs()
< std::f64::EPSILON;
if n == 0 || close {
1.0
} else {
let new_mult =
(bonds_after.get_total_weight() / bonds_before.get_total_weight()).powi(n as i32);
debug_assert!({
let valid = new_mult >= 0.;
if !valid {
println!(
"Negative multiplier: {}\t{}\t{}",
bonds_after.get_total_weight(),
bonds_before.get_total_weight(),
n
);
}
valid
});
new_mult
}
}
#[cfg(test)]
mod sc_tests {
use super::*;
use crate::sse::fast_ops::*;
use smallvec::smallvec;
#[test]
fn test_overlapping_regions_simple() {
let cutoff = 10;
let flips = [0, 2, 4, 6, 8];
let p_start = 1;
let p_end = 7;
let overlaps = find_overlapping_starts(p_start, p_end, cutoff, &flips).collect::<Vec<_>>();
println!("{:?}", overlaps);
assert_eq!(overlaps, vec![0, 1, 2, 3]);
}
#[test]
fn test_overlapping_regions() {
let cutoff = 10;
let flips = [0, 2, 4, 6, 8];
let p_start = 5;
let p_end = 7;
let overlaps = find_overlapping_starts(p_start, p_end, cutoff, &flips).collect::<Vec<_>>();
println!("{:?}", overlaps);
assert_eq!(overlaps, vec![2, 3]);
}
#[test]
fn test_wrap_around() {
let cutoff = 10;
let flips = [0, 2, 4, 6, 8];
let p_start = 7;
let p_end = 1;
let overlaps = find_overlapping_starts(p_start, p_end, cutoff, &flips).collect::<Vec<_>>();
println!("{:?}", overlaps);
assert_eq!(overlaps, vec![3, 4, 0]);
}
#[test]
fn test_remove_dups() {
let mut v = vec![0, 0, 1, 2, 3, 3];
remove_doubles(&mut v);
assert_eq!(v, vec![1, 2])
}
#[test]
fn test_remove_dups_again() {
let mut v = vec![0, 0, 1, 2, 2, 3];
remove_doubles(&mut v);
assert_eq!(v, vec![1, 3])
}
struct EN {
bonds_for_var: Vec<Vec<usize>>,
bonds: Vec<((usize, usize), bool)>,
}
impl EdgeNavigator for EN {
fn n_bonds(&self) -> usize {
self.bonds.len()
}
fn bonds_for_var(&self, var: usize) -> &[usize] {
&self.bonds_for_var[var]
}
fn vars_for_bond(&self, bond: usize) -> (usize, usize) {
self.bonds[bond].0
}
fn bond_prefers_aligned(&self, bond: usize) -> bool {
self.bonds[bond].1
}
}
fn large_joined_manager() -> (FastOps, EN) {
(
FastOps::new_from_ops(
2,
vec![
(
0,
FastOp::offdiagonal(
smallvec![0],
2,
smallvec![false],
smallvec![false],
true,
),
),
(
1,
FastOp::offdiagonal(
smallvec![1],
3,
smallvec![false],
smallvec![false],
true,
),
),
(
2,
FastOp::offdiagonal(
smallvec![0],
2,
smallvec![false],
smallvec![false],
true,
),
),
(
3,
FastOp::offdiagonal(
smallvec![1],
3,
smallvec![false],
smallvec![false],
true,
),
),
(
4,
FastOp::diagonal(smallvec![0, 1], 0, smallvec![false, false], false),
),
(
5,
FastOp::offdiagonal(
smallvec![0],
2,
smallvec![false],
smallvec![false],
true,
),
),
(
6,
FastOp::offdiagonal(
smallvec![1],
3,
smallvec![false],
smallvec![false],
true,
),
),
(
7,
FastOp::offdiagonal(
smallvec![0],
2,
smallvec![false],
smallvec![false],
true,
),
),
(
8,
FastOp::offdiagonal(
smallvec![1],
3,
smallvec![false],
smallvec![false],
true,
),
),
]
.into_iter(),
),
EN {
bonds_for_var: vec![vec![0, 1], vec![0, 1]],
bonds: vec![((0, 1), true), ((0, 1), false)],
},
)
}
#[test]
fn run_two_joined_check_flip_p() {
let (mut manager, edges) = large_joined_manager();
debug_print_diagonal(&manager, &[false, false]);
let p = calculate_flip_prob(
&mut manager,
(&[false, false], &mut [false, false]),
(&mut [false, false], &[3, 6]),
&[Some(2), Some(1)],
(&[0, 1], |v| Some(v)),
(
|b, sa, sb| {
let pref = edges.bond_prefers_aligned(b);
let aligned = sa == sb;
if aligned == pref {
0.0
} else {
1.0
}
},
|_| 1.0,
&edges,
),
);
println!("{}", p);
}
}