Browse Source

Tried to fix the problem, may have made it worse.

master
Thomas Johnson 2 years ago
parent
commit
56c3434298
  1. 2
      src/clause_gen.rs
  2. 340
      src/karmarkar.rs
  3. 2
      src/main.rs
  4. 30
      src/maxsat.rs

2
src/clause_gen.rs

@ -61,7 +61,7 @@ pub fn gen_clause_list(count: usize, max_vars: usize, max_size: usize) -> Clause
cls
}
#[cfg(test)]
#[cfg(all(test, never))]
mod test {
use super::*;

340
src/karmarkar.rs

@ -327,7 +327,7 @@ impl LeqInequality {
for term in self.rhs.drain_terms() {
self.lhs += -term;
}
self.rhs += self.lhs.get_offset();
self.rhs += -self.lhs.get_offset();
self.lhs.set_offset(0.0);
}
@ -543,6 +543,22 @@ impl StandardFormLP {
self.A[(nvars, i)] -= row.sum();
}
}
#[cfg(debug_mip)]
pub fn print_paste(&self) {
print!("maximize ");
for el in self.c.iter() {
print!("{} ", el);
}
println!("");
println!("subject to");
for (i, row) in self.A.row_iter().enumerate() {
for el in row.iter() {
print!("{} ", el);
}
println!("0 {}", self.b[i]);
}
}
}
impl Display for StandardFormLP {
@ -626,7 +642,7 @@ impl HomogenousLP {
pub fn project_b_from(&mut self, b: &DVector<Fl>) {
let n1 = self.A.ncols();
self.A.column_mut(n1 - 1).copy_from(&(n1 as f64 * b));
self.A.column_mut(n1 - 1).copy_from(&(-(n1 as f64) * b));
}
pub fn project_from_sflp(sflp: StandardFormLP) -> Self {
@ -642,6 +658,26 @@ impl HomogenousLP {
self.A.copy_from(&other.A);
self.c.copy_from(&other.c);
}
#[cfg(debug_mip)]
pub fn print_paste(&self) {
print!("maximize ");
for el in self.c.iter() {
print!("{} ", el);
}
println!("");
println!("subject to");
for (i, row) in self.A.row_iter().enumerate() {
for el in row.iter() {
print!("{} ", el);
}
println!("0 0");
}
for _ in 0..self.A.ncols() {
print!("1 ");
}
println!("0 1");
}
}
impl Display for HomogenousLP {
@ -720,7 +756,9 @@ impl Karmarkar {
#[cfg(debug_karmarkar)]
{
print!("guess: {}", &self.guess);
print!("infeasiblity: {}", &problem.A * &self.guess);
print!("objective: {}", &self.problem.c.dot(&self.guess));
}
// First transform c, giving Dc
@ -736,6 +774,10 @@ impl Karmarkar {
}
self.B.row_mut(self.nconstr).fill(1.0);
}
#[cfg(debug_karmarkar)]
{
print!("B: {}", self.B);
}
self.B.transpose_to(&mut self.BT);
self.B.mul_to(&self.BT, &mut self.BBT);
if !self.BBT.try_inverse_mut() {
@ -747,16 +789,71 @@ impl Karmarkar {
self.BTBBTB[(i, i)] -= 1.0;
}
self.BTBBTB.mul_to(&self.Dc, &mut self.cp);
#[cfg(debug_karmarkar)]
{
for el in self.BTBBTB.iter_mut() {
if el.abs() < DEFAULT_ε {
*el = 0.0;
}
}
print!("BTBBTB: {}", &self.BTBBTB);
let mut check = (&self.B.transpose()
* (&self.B * &self.B.transpose()).try_inverse().unwrap()
* &self.B);
check.iter_mut().for_each(|x| {
if x.abs() < DEFAULT_ε {
*x = 0.0;
}
});
print!("BTBBTBi: {}", &check);
print!(
"ATAATAi: {}",
&problem.A.transpose()
* (&problem.A * &problem.A.transpose()).try_inverse().unwrap()
* &problem.A
);
print!("B * cp: {}", &self.B * &self.cp);
print!(
"AD * cp: {}",
&self.problem.A * DMatrix::from_diagonal(&self.guess) * &self.cp
);
print!(
"AD * a0: {}",
&self.problem.A
* DMatrix::from_diagonal(&self.guess)
* DVector::from_element(self.guess.len(), self.a0)
);
print!("D * cp: {}", DMatrix::from_diagonal(&self.guess) * &self.cp);
print!("cp {}", &self.cp);
}
// Now normalize the result and step in that direction
self.cp.normalize_mut();
self.cp *= -self.αr;
#[cfg(debug_karmarkar)]
print!("^c {}", &self.cp);
// This part shouldn't be necessary. It's a hack to avoid some of the perils of floating
// point numbers.
self.BTBBTB.mul_to(&self.cp, &mut self.Dc);
self.cp.copy_from(&self.Dc);
#[cfg(debug_karmarkar)]
{
print!("^c, revised: {}", -&self.cp);
print!(
"AD * ^c: {}",
&self.problem.A * DMatrix::from_diagonal(&self.guess) * &self.cp
);
}
self.cp *= self.αr;
self.cp.add_scalar_mut(self.a0);
#[cfg(debug_karmarkar)]
print!("A * D * b': {}", &self.problem.A * &self.cp);
// And finally, transform back into the problem space.
self.cp.component_mul_assign(&self.guess);
let scale = 1.0 / self.cp.sum();
self.cp *= scale;
if self.cp.iter().all(|x| !x.is_nan()) {
let est_dist = self.guess.metric_distance(&self.cp) * self.;
#[cfg(debug_karmarkar)]
println!("estimated remaining distance: {}", est_dist);
self.guess.copy_from(&self.cp);
Distance(est_dist)
} else {
@ -781,7 +878,12 @@ impl Karmarkar {
self.ε
};
if d < tolerance {
#[cfg(debug_karmarkar)]
println!("exiting. distance: {}", d);
break;
} else {
#[cfg(debug_karmarkar)]
println!("insufficient accuracy. distance: {}", d);
}
}
_ => {}
@ -880,11 +982,12 @@ impl BMIPSolver {
DVector::zeros(problem.A.ncols() + 2),
)
.unwrap();
let feasible_point = DVector::zeros(problem.A.ncols() + 1);
problems.push(problem);
projected_problems.push(proj_problem);
feasibility_problems.push(feas_problem);
projected_feasibility_problems.push(proj_feas_problem);
feasible_points.push(DVector::zeros(0)); // The first one should never end up being used.
feasible_points.push(feasible_point); // The first one should never end up being used.
objective_offsets.push(0.0);
for i in 0..bvars.len() {
let problem = Self::remove_var(problems[i].clone());
@ -931,6 +1034,13 @@ impl BMIPSolver {
&self.projected_feasibility_problems[nassigned].guess,
&mut self.feasible_points[nassigned],
);
#[cfg(debug_mip_verbose)]
{
print!(
"feasibility solution vector: {}",
&self.feasible_points[nassigned]
);
}
if self.feasible_points[nassigned][nvars] > DEFAULT_ε {
// The problem as given isn't feasible
return None;
@ -946,28 +1056,55 @@ impl BMIPSolver {
&self.projected_problems[nassigned].guess,
&mut self.feasible_points[nassigned].rows_mut(0, nvars),
);
Some(
self.problems[nassigned]
.c
.dot(&self.feasible_points[nassigned].rows_mut(0, nvars))
+ *self.objective_offsets.last().unwrap(),
)
let obj = self.problems[nassigned]
.c
.dot(&self.feasible_points[nassigned].rows(0, nvars))
+ *self.objective_offsets.last().unwrap();
#[cfg(debug_mip_verbose)]
{
print!(
"solution vector: {}",
&self.feasible_points[nassigned].rows(0, nvars)
);
println!("objective function: {}", obj);
}
Some(obj)
}
// `true` if we can prune
fn can_prune(&mut self, nassigned: usize) -> bool {
match self.best_solution {
None => false,
Some((best, _)) => match self.solve_at_depth(nassigned) {
None => true,
Some(value) => value < best + DEFAULT_ε,
},
match self.solve_at_depth(nassigned) {
None => {
#[cfg(debug_mip)]
println!("{: ^width$}no solution!", "", width = nassigned);
true
}
Some(value) => {
#[cfg(debug_mip)]
println!(
"{: ^width$}got {} for objective function",
"",
value,
width = nassigned
);
match self.best_solution {
None => false,
Some((best, _)) => value < best,
}
}
}
}
fn record_maximum(&mut self, value: Fl) {
#[cfg(debug_mip)]
println!("{: ^width$}record!", "", width = self.bvars.len());
println!("{: ^width$}record! {}", "", value, width = self.bvars.len());
#[cfg(debug_mip_verbose)]
{
for i in 0..self.bvars.len() {
print!("{}", self.problems[i].A);
print!("{}", self.problems[i].b);
}
}
let best_solution = match &mut self.best_solution {
x @ None => {
*x = Some((value, DVector::zeros(self.bvars.len())));
@ -993,6 +1130,9 @@ impl BMIPSolver {
}
}
}
} else {
#[cfg(debug_mip)]
println!("{: ^width$}no solution!", "", width = self.bvars.len());
}
}
@ -1005,6 +1145,19 @@ impl BMIPSolver {
.problem
.project_b_from(&newb);
self.feasibility_problems[nassigned + 1].b.copy_from(&newb);
for (i, el) in self.feasibility_problems[nassigned + 1]
.A
.column_mut(nvars - 1)
.iter_mut()
.enumerate()
{
*el = newb[i] - self.problems[nassigned + 1].A.row(i).sum();
}
self.projected_feasibility_problems[nassigned + 1]
.problem
.A
.column_mut(nvars - 1)
.copy_from(&self.feasibility_problems[nassigned + 1].A.column(nvars - 1));
self.projected_feasibility_problems[nassigned + 1]
.problem
.project_b_from(&newb);
@ -1019,19 +1172,43 @@ impl BMIPSolver {
fn search_at_depth(&mut self, nassigned: usize) {
#[cfg(debug_mip)]
println!(
"{: ^width$}searching at depth {}",
"",
nassigned,
width = nassigned
);
{
println!(
"{: ^width$}searching at depth {}",
"",
nassigned,
width = nassigned
);
print!("{: ^width$}current branch:", "", width = nassigned);
self.bvars
.iter()
.rev()
.filter_map(|x| x.get_assignment())
.for_each(|x| print!("{} ", x));
println!("");
}
#[cfg(all(debug_mip, debug_mip_verbose))]
{
println!("problem:");
self.problems[nassigned].print_paste();
//print!("{}", self.problems[nassigned]);
println!("feasibility problem:");
self.feasibility_problems[nassigned].print_paste();
//print!("{}", self.feasibility_problems[nassigned]);
println!("projected feasibility problem:");
print!("{}", self.projected_feasibility_problems[nassigned].problem);
}
if nassigned == self.bvars.len() {
self.accumulate_maximum(nassigned);
} else if !self.can_prune(nassigned) {
self.assign(nassigned, 0.0);
#[cfg(debug_mip)]
println!("{: ^width$}trying 0", "", width = nassigned);
self.search_at_depth(nassigned + 1);
self.unassign(nassigned);
self.assign(nassigned, 1.0);
#[cfg(debug_mip)]
println!("{: ^width$}trying 1", "", width = nassigned);
self.search_at_depth(nassigned + 1);
self.unassign(nassigned);
} else {
@ -1096,6 +1273,14 @@ where
if let Some(Some(weight)) = weights.get(i) {
counter += 1;
let b = clause::Variable::from_name(format!("_b{}", counter));
#[cfg(debug_mip)]
{
if let Some(true) = clause.eval() {
b.assign(false);
} else {
b.assign(true);
}
}
weightmap.insert(b.clone(), Fl::from(*weight));
clause.push(b.get_pos_literal());
}
@ -1105,6 +1290,12 @@ where
let vars = clause_list.get_vars();
for var in &vars {
let lvar = Variable::from_name(var.get_name().clone());
#[cfg(debug_mip)]
{
var.get_assignment()
.map(|x| if x { 1.0 } else { 0.0 })
.map(|a| lvar.assign(a));
}
weightmap
.get(&var)
.map(|x| objective.insert_variable(lvar.clone(), -*x));
@ -1122,21 +1313,55 @@ where
Fl: From<T>,
{
let (ineqs, objective, varmap) = clause_list_to_constraints(&mut clause_list, weights);
println!("maximize {}", objective);
println!("subject to");
for ineq in &ineqs {
println!(" {}", ineq);
#[cfg(debug_mip)]
{
println!("maximize {}", objective);
println!("subject to");
for ineq in &ineqs {
println!(" {}", ineq);
}
}
let varset = varmap.values().cloned();
let mut solver = BMIPSolver::from_ineqs_and_objective(ineqs, varset, objective);
#[cfg(debug_mip)]
{
print!("variables (forward): ");
solver
.bvars
.iter()
.rev()
.for_each(|x| print!("{} ", x.get_name()));
println!("");
let rltable = varmap
.iter()
.map(|(k, v)| (v.clone(), k.clone()))
.collect::<HashMap<_, _>>();
print!("valid assignment: ");
solver
.bvars
.iter()
.rev()
.filter_map(|x| rltable.get(&x))
.map(|x| x.get_assignment())
.for_each(|x| {
x.map(|x| if x { 1 } else { 0 })
.map(|x| print!("{} ", x))
.unwrap_or_else(|| print!("- "))
});
println!("");
solver.bvars.iter().skip(1).for_each(|x| x.unassign());
}
#[cfg(debug_mip)]
println!("solving");
solver.solve();
println!("done");
#[cfg(debug_mip)]
println!("done: {:?}", solver.best_solution.map(|x| x.0));
for (bv, lv) in varmap {
lv.get_assignment()
.map(|x| x > 0.5)
.map(|x| bv.assign(x))
.or_else(|| {
#[cfg(debug_mip)]
println!("unassigning");
bv.unassign();
None
@ -1150,7 +1375,7 @@ mod test {
// Designed to be a --nocapture test for human analysis. Pretty much just for debugging.
#[test]
fn inequality_test() {
fn inequality_test_1() {
linear_vars![x, y];
let mut constraints = Vec::<LeqInequality>::new();
constraints.push(&x * 50.0 + &y * 24.0 << 2400.0);
@ -1225,4 +1450,61 @@ mod test {
println!("vars: {:?}", vars);
print!("final x: {}", solution);
}
#[test]
fn inequality_test_2() {
let A = DMatrix::from_row_slice(
7,
8,
&[
0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0,
0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0,
],
);
let c = DVector::from_row_slice(&[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0]);
let b = DVector::from_row_slice(&[1.0, 1.0, 0.0, 2.0, 1.0, 1.0, 0.0]);
let system = StandardFormLP { A, b, c };
system.print_paste();
let proj_system = HomogenousLP::project_from_sflp(system.clone());
println!("projected system:");
proj_system.print_paste();
print!("{}", proj_system);
let fp = DVector::from_element(8, 1.0);
let proj_fp = HomogenousLP::project_point(&fp);
let mut km = Karmarkar::from_feasible_value(proj_system.clone(), proj_fp).unwrap();
km.iterate_until_failed_or_precise(true, 1000);
print!("final x': {}", km.guess);
let solution = HomogenousLP::project_point_inv(&km.guess);
print!("final x: {}", solution);
}
#[test]
fn inequality_test_3() {
let A = DMatrix::from_row_slice(
7,
8,
&[
0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0,
1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0,
0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -2.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0,
0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0,
],
);
let c = DVector::from_row_slice(&[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0]);
let b = DVector::from_row_slice(&[0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0]);
let system = StandardFormLP { A, b, c };
let proj_system = HomogenousLP::project_from_sflp(system.clone());
println!("projected system:");
proj_system.print_paste();
print!("{}", proj_system);
let fp = DVector::from_element(8, 1.0);
let proj_fp = HomogenousLP::project_point(&fp);
let mut km = Karmarkar::from_feasible_value(proj_system.clone(), proj_fp).unwrap();
km.iterate_until_failed_or_precise(true, 1000000);
print!("final x': {}", km.guess);
let solution = HomogenousLP::project_point_inv(&km.guess);
print!("final x: {}", solution);
}
}

2
src/main.rs

@ -149,7 +149,7 @@ fn main() -> io::Result<()> {
//dispatcher(16, 10, job_runner_generator, get_job, handler);
let _ = result_loop.join();
// let _ = result_loop.join();
exit(0);
}

30
src/maxsat.rs

@ -23,7 +23,35 @@ fn bfs_step(cs: &ClauseList, vs: &Vec<Variable>, n: usize) -> usize {
}
}
fn bfs_step_bounded(cs: &ClauseList, vs: &Vec<Variable>, n: usize, bound: usize) -> bool {
if n == vs.len() {
let r: usize = cs
.iter()
.filter_map(|c| c.eval())
.map(|b| if b { 1 } else { 0 })
.sum();
#[cfg(debug_bfs)]
println!("{}\t= {}", cs, r);
r == bound
} else {
let a = {
vs[n].assign(false);
bfs_step_bounded(cs, vs, n + 1, bound)
};
if a {
true
} else {
{
vs[n].assign(true);
bfs_step_bounded(cs, vs, n + 1, bound)
}
}
}
}
pub fn solve_by_bfs(cs: &ClauseList) -> usize {
let vars: Vec<Variable> = cs.get_vars().into_iter().collect();
bfs_step(cs, &vars, 0)
let bound = bfs_step(cs, &vars, 0);
bfs_step_bounded(cs, &vars, 0, bound);
bound
}
Loading…
Cancel
Save