Skip to content

Commit

Permalink
Removed constraint array and refactored APIS.
Browse files Browse the repository at this point in the history
  • Loading branch information
@frankpermenter
frankpermenter committed Dec 24, 2020
1 parent 06c0f58 commit 4334d35
Show file tree
Hide file tree
Showing 9 changed files with 172 additions and 132 deletions.
34 changes: 9 additions & 25 deletions conex/cone_program.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,11 +68,10 @@ void GetMuSelectionParameters(ConstraintManager<Container>* constraints,
}
}

template <typename T>
int Rank(const std::vector<T*>& c) {
int Rank(const ConstraintManager<Container>& kkt) {
int rank = 0;
for (const auto& ci : c) {
rank += Rank(*ci);
for (const auto& ci : kkt.eqs) {
rank += Rank(ci.constraint);
}
return rank;
}
Expand All @@ -92,14 +91,14 @@ void Initialize(Program& prog, const SolverConfiguration& config) {
auto& solver = prog.solver;
auto& kkt = prog.kkt;
prog.InitializeWorkspace();
SetIdentity(&prog.constraints);

solver = std::make_unique<Solver>(prog.kkt_system_manager_.cliques,
prog.kkt_system_manager_.dual_vars);

kkt.clear();
for (auto& c : prog.kkt_system_manager_.eqs) {
c.kkt_assembler.Reset();
SetIdentity(&c.constraint);
}

for (auto& c : prog.kkt_system_manager_.eqs) {
Expand Down Expand Up @@ -145,10 +144,7 @@ bool Solve(const DenseMatrix& b, Program& prog,
std::cout << "CONEX: MKL Enabled";
#endif

auto& constraints = prog.constraints;
auto& sys = prog.sys;
auto& solver = prog.solver;

bool solved = 1;

std::cout.precision(2);
Expand Down Expand Up @@ -179,7 +175,7 @@ bool Solve(const DenseMatrix& b, Program& prog,
newton_step_parameters.inv_sqrt_mu = 0;
newton_step_parameters.affine = false;

int rankK = Rank(constraints);
int rankK = Rank(prog.kkt_system_manager_);
int centering_steps = 0;

Eigen::VectorXd AW(prog.kkt_system_manager_.SizeOfKKTSystem());
Expand Down Expand Up @@ -256,22 +252,10 @@ bool Solve(const DenseMatrix& b, Program& prog,

if (config.prepare_dual_variables) {
DenseMatrix y2;
bool iter_refine = false;
if (iter_refine) {
ConstructSchurComplementSystem(&constraints, true, &sys);
Eigen::LLT<Eigen::Ref<DenseMatrix>> llt(sys.G);
DenseMatrix L = llt.matrixL();
DenseMatrix bres = newton_step_parameters.inv_sqrt_mu * b - 1 * sys.AW;
y2 = bres * 0;
for (int i = 0; i < 1; i++) {
y2 += llt.solve(bres - L * L.transpose() * y2);
}
} else {
solver->Assemble(&AW, &AQc);
solver->Factor();
DenseMatrix bres = newton_step_parameters.inv_sqrt_mu * b - 1 * AW;
y2 = solver->Solve(bres);
}
solver->Assemble(&AW, &AQc);
solver->Factor();
DenseMatrix bres = newton_step_parameters.inv_sqrt_mu * b - 1 * AW;
y2 = solver->Solve(bres);

newton_step_parameters.affine = true;
newton_step_parameters.e_weight = 0;
Expand Down
65 changes: 56 additions & 9 deletions conex/cone_program.h
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
#pragma once
#include "conex/error_checking_macros.h"
#include "conex/constraint.h"
#include "conex/constraint_manager.h"
#include "conex/equality_constraint.h"
#include "conex/kkt_assembler.h"
#include "conex/kkt_solver.h"
#include "constraint.h"
#include "workspace.h"

namespace conex {
Expand Down Expand Up @@ -44,6 +46,54 @@ class Program {

int GetNumberOfVariables() { return sys.m_; }


template<typename T>
void GetDualVariable(int i, T* xi) {
int cnt = 0;
for (auto& ci : kkt_system_manager_.eqs) {
if (cnt == i) {
ci.constraint.get_dual_variable(xi->data());
xi->array() /= stats.sqrt_inv_mu[stats.num_iter - 1];
return;
}
cnt++;
}
}

int GetDualVariableSize(int i) {
int cnt = 0;
for (auto& ci : kkt_system_manager_.eqs) {
if (cnt == i) {
return ci.constraint.dual_variable_size();
}
cnt++;
}
CONEX_DEMAND(false, "Invalid Constraint");
}

int UpdateLinearOperatorOfConstraint(int i, double value, int variable,
int row, int col, int hyper_complex_dim) {
int cnt = 0;
for (auto& ci : kkt_system_manager_.eqs) {
if (cnt == i) {
return UpdateLinearOperator(&ci.constraint, value, variable, row, col, hyper_complex_dim);
}
cnt++;
}
CONEX_DEMAND(false, "Invalid Constraint");
}

int UpdateAffineTermOfConstraint(int i, double value, int row, int col, int hyper_complex_dim) {
int cnt = 0;
for (auto& ci : kkt_system_manager_.eqs) {
if (cnt == i) {
return UpdateAffineTerm(&ci.constraint, value, row, col, hyper_complex_dim);
}
cnt++;
}
CONEX_DEMAND(false, "Invalid Constraint");
}

int SetNumberOfConstraints() const { return kkt_system_manager_.eqs.size(); }

void InitializeWorkspace() {
Expand All @@ -55,30 +105,27 @@ class Program {
memory.resize(SizeOf(workspaces));
Initialize(&workspaces, &memory[0]);

constraints.clear();
for (auto& ci : kkt_system_manager_.eqs) {
constraints.push_back(&ci.constraint);
}

is_initialized = true;
}

template <typename T>
void AddConstraint(T&& d) {
kkt_system_manager_.AddConstraint(d);
constraints.push_back(&kkt_system_manager_.eqs.back().constraint);
}

template <typename T>
void AddConstraint(T&& d, const std::vector<int>& variables) {
kkt_system_manager_.AddConstraint(d, variables);
constraints.push_back(&kkt_system_manager_.eqs.back().constraint);
}

// void AddConstraint(EqualityConstraints&& d, const std::vector<int>&
// variables) {
// kkt_system_manager_.AddConstraint(d, variables);
// }

int NumberOfConstraints() { return kkt_system_manager_.eqs.size(); }

ConstraintManager<Container> kkt_system_manager_;
std::vector<Constraint*> constraints;
SchurComplementSystem sys;
WorkspaceStats stats;
std::vector<Workspace> workspaces;
Expand Down
21 changes: 9 additions & 12 deletions conex/test/test_lp.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,9 +146,9 @@ Eigen::VectorXd SolveSparseHelper(bool sparse) {
VectorXd slack = C.at(i) - A.at(i) * Vars(y, variables.at(i));
EXPECT_TRUE((slack).minCoeff() >= -eps);

VectorXd xi(A.at(i).rows());
prog.constraints.at(i)->get_dual_variable(xi.data());
xi.array() /= prog.stats.sqrt_inv_mu[prog.stats.num_iter - 1];
MatrixXd xi(A.at(i).rows(), 1);
prog.GetDualVariable(i, &xi);

VectorXd temp = A.at(i).transpose() * xi;
int k = 0;
for (auto vk : variables.at(i)) {
Expand All @@ -165,9 +165,8 @@ Eigen::VectorXd SolveSparseHelper(bool sparse) {
auto L = Combine(A, C, variables);
VectorXd slack = L.constraint_affine_ - L.constraint_matrix_ * y;

VectorXd xi(L.constraint_affine_.rows());
prog.constraints.at(0)->get_dual_variable(xi.data());
xi.array() /= prog.stats.sqrt_inv_mu[prog.stats.num_iter - 1];
MatrixXd xi(L.constraint_affine_.rows(), 1);
prog.GetDualVariable(0, &xi);
EXPECT_NEAR((b - L.constraint_matrix_.transpose() * xi).norm(), 0, 1e-8);
}

Expand Down Expand Up @@ -234,9 +233,8 @@ Eigen::VectorXd SolveFillIn(bool sparse) {
VectorXd slack = C.at(i) - A.at(i) * Vars(y, variables.at(i));
EXPECT_TRUE((slack).minCoeff() >= -eps);

VectorXd xi(A.at(i).rows());
prog.constraints.at(i)->get_dual_variable(xi.data());
xi.array() /= prog.stats.sqrt_inv_mu[prog.stats.num_iter - 1];
MatrixXd xi(A.at(i).rows(), 1);
prog.GetDualVariable(i, &xi);
VectorXd temp = A.at(i).transpose() * xi;
int k = 0;
for (auto vk : variables.at(i)) {
Expand All @@ -253,9 +251,8 @@ Eigen::VectorXd SolveFillIn(bool sparse) {
auto L = Combine(A, C, variables);
VectorXd slack = L.constraint_affine_ - L.constraint_matrix_ * y;

VectorXd xi(L.constraint_affine_.rows());
prog.constraints.at(0)->get_dual_variable(xi.data());
xi.array() /= prog.stats.sqrt_inv_mu[prog.stats.num_iter - 1];
MatrixXd xi(L.constraint_affine_.rows(), 1);
prog.GetDualVariable(0, &xi);
EXPECT_NEAR((b - L.constraint_matrix_.transpose() * xi).norm(), 0, 1e-8);
}

Expand Down
4 changes: 2 additions & 2 deletions examples/performance_plots.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -308,9 +308,9 @@ def SolveRandomSDP(num_variables, n, config, w0 = []):
#PlotMuUpdate(1, "Real", show_plot)
#PlotMuUpdate(2, "Complex", show_plot)
#PlotMuUpdate(4, "Quaternion", show_plot)
#PlotMuUpdate(-1, "special", show_plot)
PlotMuUpdate(-1, "special", show_plot)

PlotMuUpdate(-2, "sdp", show_plot)
#lotMuUpdate(-2, "sdp", show_plot)
#PlotMuUpdate(8, "exceptional", show_plot)
#PlotMuUpdateVsDivBound()

48 changes: 22 additions & 26 deletions interfaces/conex.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,19 +16,19 @@
#include "conex/error_checking_macros.h"
// TODO(FrankPermenter): check for null pointers.

#define SAFER_CAST_TO_Program(x, prog) \
CONEX_DEMAND(x, "Program pointer is null."); \
prog = static_cast<Program*>(x); \
if (prog->is_initialized) { \
if (prog->NumberOfConstraints() + 2 != \
static_cast<int>(prog->workspaces.size())) { \
CONEX_DEMAND(false, "Program corrupted or invalid pointer."); \
} \
} else { \
if (prog->workspaces.size() != 0) { \
CONEX_DEMAND(false, "Program corrupted or invalid pointer."); \
} \
} \
#define SAFER_CAST_TO_Program(x, prog) \
CONEX_DEMAND(x, "Program pointer is null."); \
prog = static_cast<Program*>(x); \
if (prog->is_initialized) { \
if (prog->NumberOfConstraints() + 2 != \
static_cast<int>(prog->workspaces.size())) { \
CONEX_DEMAND(false, "Program corrupted or invalid pointer."); \
} \
} else { \
if (prog->workspaces.size() != 0 || prog->NumberOfConstraints() < 0) { \
CONEX_DEMAND(false, "Program corrupted or invalid pointer."); \
} \
} \
CONEX_DEMAND(prog, "Program corrupted or invalid pointer.");

using DenseMatrix = Eigen::MatrixXd;
Expand All @@ -38,7 +38,7 @@ using conex::Program;
using conex::SolverConfiguration;

int CONEX_Solve(void* prog_ptr, const double* b, int br,
const CONEX_SolverConfiguration* config, double* y, int yr) {
const CONEX_SolverConfiguration* config, double* y, int yr) {
using InputMatrix = Eigen::Map<const DenseMatrix>;
InputMatrix bmap(b, br, 1);
DenseMatrix blinear = bmap;
Expand All @@ -61,19 +61,17 @@ int CONEX_Solve(void* prog_ptr, const double* b, int br,

void CONEX_GetDualVariable(void* prog_ptr, int i, double* x, int xr, int xc) {
Program& prog = *reinterpret_cast<Program*>(prog_ptr);
assert(prog.constraints.at(i)->dual_variable_size() == xr * xc);
assert(prog.GetDualVariableSize(i) == xr * xc);

using InputMatrix = Eigen::Map<DenseMatrix>;

prog.constraints.at(i)->get_dual_variable(x);

InputMatrix xmap(x, xr, xc);
xmap.array() /= prog.stats.sqrt_inv_mu[prog.stats.num_iter - 1];
prog.GetDualVariable(i, &xmap);
}

int CONEX_GetDualVariableSize(void* prog_ptr, int i) {
Program& prog = *reinterpret_cast<Program*>(prog_ptr);
return prog.constraints.at(i)->dual_variable_size();
return prog.GetDualVariableSize(i);
}

void* CONEX_CreateConeProgram() {
Expand Down Expand Up @@ -197,7 +195,6 @@ void CONEX_GetIterationStats(void* prog, CONEX_IterationStats* stats,
stats->iteration_number = iter_num;
}


CONEX_STATUS CONEX_NewLinearMatrixInequality(void* p, int order,
int hyper_complex_dim,
int* constraint_id) {
Expand Down Expand Up @@ -235,18 +232,17 @@ CONEX_STATUS CONEX_UpdateLinearOperator(void* p, int constraint, double value,
Program* prg;
SAFER_CAST_TO_Program(p, prg);
CONEX_DEMAND(constraint < prg->NumberOfConstraints(), "Invalid Constraint.");
return UpdateLinearOperator(prg->constraints.at(constraint), value, variable,
row, col, hyper_complex_dim);
return prg->UpdateLinearOperatorOfConstraint(constraint, value, variable, row,
col, hyper_complex_dim);
}

CONEX_STATUS CONEX_UpdateAffineTerm(void* p, int constraint, double value,
int row, int col, int hyper_complex_dim) {
Program* prg;
SAFER_CAST_TO_Program(p, prg);
CONEX_DEMAND(constraint < static_cast<int>(prg->constraints.size()),
"Invalid Constraint.");
return UpdateAffineTerm(prg->constraints.at(constraint), value, row, col,
hyper_complex_dim);
CONEX_DEMAND(constraint < prg->NumberOfConstraints(), "Invalid Constraint.");
return prg->UpdateAffineTermOfConstraint(constraint, value, row, col,
hyper_complex_dim);
}

CONEX_STATUS CONEX_NewLorentzConeConstraint(void* p, int order,
Expand Down
2 changes: 1 addition & 1 deletion interfaces/conex.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ int CONEX_AddSparseLMIConstraint(void* prog, const double* Aarray, int Aarrayr,
int cc, const long* vars, int vars_c);

int CONEX_Solve(void* prog, const double* b, int br,
const CONEX_SolverConfiguration* config, double* y, int yr);
const CONEX_SolverConfiguration* config, double* y, int yr);

void CONEX_GetDualVariable(void* prog, int i, double* x, int xr, int xc);

Expand Down
Loading

0 comments on commit 4334d35

Please sign in to comment.