Skip to content

Commit

Permalink
Tried various things. Still need to make it work
Browse files Browse the repository at this point in the history
  • Loading branch information
@bstellato
bstellato committed Sep 11, 2017
1 parent d3fecb9 commit 3590088
Showing 1 changed file with 95 additions and 62 deletions.
157 changes: 95 additions & 62 deletions interfaces/python/modulepurepy/_osqp.py
View file
Expand Up @@ -401,7 +401,9 @@ def scale_data(self):
s_temp = np.ones(n + m)
c = 1.0 # Cost scaling
c_temp = 1.0
c_temp_prev = c_temp
# c_temp_prev = c_temp
# avg_norm_P_cols = 1.0
# avg_norm_P_cols_prev = 1.0

# Define reduced KKT matrix to scale
# KKT = spspa.vstack([
Expand Down Expand Up @@ -451,13 +453,15 @@ def scale_data(self):
D = D_temp.dot(D)
E = E_temp.dot(E)

# # Second Step cost normalization
# Second Step cost normalization
# avg_norm_P_cols = spla.norm(P, np.inf, axis=0).mean()
# inf_norm_q = np.linalg.norm(q, np.inf)
# # scale_cost = np.maximum(inf_norm_q, avg_norm_P_cols)
# scale_cost = np.maximum(inf_norm_q, avg_norm_P_cols)
# scale_cost = 1. / np.maximum(np.minimum(
# scale_cost, MAX_SCALING), MIN_SCALING)
# if inf_norm_q > 1e-06:
# scale_cost = np.maximum(inf_norm_q, avg_norm_P_cols)
# scale_cost = 1. / np.maximum(np.minimum(
# scale_cost, MAX_SCALING), MIN_SCALING)
# else:
# scale_cost = 1.0
#
# print("avg_norm_P_cols", avg_norm_P_cols)
# print("inf_norm_q", inf_norm_q)
Expand All @@ -471,24 +475,28 @@ def scale_data(self):
# # c_temp = 1.0
#
# # Normalize cost
# # P = c_temp * P
# # q = c_temp * q
# P = c_temp * P
# q = c_temp * q
#
# # Update scaling
# c = c_temp * c
#
# print("Total scale cost = %.2e" % c)

# if (abs(c_temp - c_temp_prev) < 1e-03):
# if (abs(avg_norm_P_cols_prev - avg_norm_P_cols) < 1e-03):
# # Stop loop if it converges
# break

# Independent cost normalization
avg_norm_P_cols = spla.norm(P, np.inf, axis=0).mean()
inf_norm_q = np.linalg.norm(q, np.inf)
scale_cost = np.maximum(inf_norm_q, avg_norm_P_cols)
if inf_norm_q > 1e-06:
scale_cost = np.maximum(inf_norm_q, avg_norm_P_cols)
else:
scale_cost = 1.0
# scale_cost = inf_norm_q + avg_norm_P_cols / 2

scale_cost = inf_norm_q if inf_norm_q > 1e-06 else 1
# scale_cost = inf_norm_q if inf_norm_q > 1e-06 else 1

scale_cost = 1. / np.maximum(np.minimum(
scale_cost, MAX_SCALING), MIN_SCALING)
Expand All @@ -503,8 +511,8 @@ def scale_data(self):
# Normalize cost
P = c * P
q = c * q

print("Total scale cost = %.2e" % c)
#
# print("Total scale cost = %.2e" % c)

# DEBUG: check distance l and u
# test_lu = u - l
Expand Down Expand Up @@ -543,67 +551,92 @@ def scale_data(self):
self.work.scaling.c = c
self.work.scaling.cinv = 1. / c

def set_rho_vec(self):
def rho_i(self, l_i, u_i):
"""
Set values of rho vector based on constraint types
Compute rho_i from l_i, u_i
"""
self.work.settings.rho = np.minimum(np.maximum(self.work.settings.rho,
RHO_MIN), RHO_MAX)
RHO_0 = 1e6
RHO_INF = 1e-1
BETA_RHO = 1e7

# Find indices of loose bounds, equality constr and one-sided constr
loose_ind = np.where(np.logical_and(
self.work.data.l < -OSQP_INFTY*MIN_SCALING,
self.work.data.u > OSQP_INFTY*MAX_SCALING))[0]
eq_ind = np.where(self.work.data.u - self.work.data.l < RHO_TOL)[0]
ineq_ind = np.setdiff1d(np.setdiff1d(np.arange(self.work.data.m),
loose_ind), eq_ind)
m = u_i - l_i

# Type of constraints
self.work.constr_type[loose_ind] = -1
self.work.constr_type[eq_ind] = 1
self.work.constr_type[ineq_ind] = 0
rho_i = (RHO_0 + RHO_INF * BETA_RHO * m) / (1 + BETA_RHO * m)

self.work.rho_vec[loose_ind] = RHO_MIN
self.work.rho_vec[eq_ind] = RHO_MAX
self.work.rho_vec[ineq_ind] = self.work.settings.rho

self.work.rho_inv_vec = np.reciprocal(self.work.rho_vec)
return rho_i

def update_rho_vec(self):
def set_rho_vec(self):
"""
Update values of rho_vec and return 1 if type of some constraints
changed.
Set values of rho vector based on constraint types
"""
# Find indices of loose bounds, equality constr and one-sided constr
loose_ind = np.where(np.logical_and(
self.work.data.l < -OSQP_INFTY*MIN_SCALING,
self.work.data.u > OSQP_INFTY*MAX_SCALING))[0]
eq_ind = np.where(self.work.data.u - self.work.data.l < RHO_TOL)[0]
ineq_ind = np.setdiff1d(np.setdiff1d(np.arange(self.work.data.m),
loose_ind), eq_ind)

# Find indices of current constraint types
old_loose_ind = np.where(self.work.constr_type == -1)
old_eq_ind = np.where(self.work.constr_type == 1)
old_ineq_ind = np.where(self.work.constr_type == 0)

# Check if type of any constraint changed
constr_type_changed = (loose_ind != old_loose_ind).any() or \
(eq_ind != old_eq_ind).any() or \
(ineq_ind != old_ineq_ind).any()

# Update type of constraints
self.work.constr_type[loose_ind] = -1
self.work.constr_type[eq_ind] = 1
self.work.constr_type[ineq_ind] = 0
# New stuff
l = self.work.data.l
u = self.work.data.u

self.work.rho_vec[loose_ind] = RHO_MIN
self.work.rho_vec[eq_ind] = RHO_MAX
self.work.rho_vec[ineq_ind] = self.work.settings.rho
# Set rho based on smooth rule depending on slacks
for i in range(self.work.data.m):
self.work.rho_vec[i] = self.rho_i(l[i], u[i])

self.work.rho_inv_vec = np.reciprocal(self.work.rho_vec)

return constr_type_changed

# self.work.settings.rho = np.minimum(np.maximum(self.work.settings.rho,
# RHO_MIN), RHO_MAX)
#
# # Find indices of loose bounds, equality constr and one-sided constr
# loose_ind = np.where(np.logical_and(
# self.work.data.l < -OSQP_INFTY*MIN_SCALING,
# self.work.data.u > OSQP_INFTY*MAX_SCALING))[0]
# eq_ind = np.where(self.work.data.u - self.work.data.l < RHO_TOL)[0]
# ineq_ind = np.setdiff1d(np.setdiff1d(np.arange(self.work.data.m),
# loose_ind), eq_ind)
#
# # Type of constraints
# self.work.constr_type[loose_ind] = -1
# self.work.constr_type[eq_ind] = 1
# self.work.constr_type[ineq_ind] = 0
#
# self.work.rho_vec[loose_ind] = RHO_MIN
# self.work.rho_vec[eq_ind] = RHO_MAX
# self.work.rho_vec[ineq_ind] = self.work.settings.rho
#
# self.work.rho_inv_vec = np.reciprocal(self.work.rho_vec)

# def update_rho_vec(self):
# """
# Update values of rho_vec and return 1 if type of some constraints
# changed.
# """
# # Find indices of loose bounds, equality constr and one-sided constr
# loose_ind = np.where(np.logical_and(
# self.work.data.l < -OSQP_INFTY*MIN_SCALING,
# self.work.data.u > OSQP_INFTY*MAX_SCALING))[0]
# eq_ind = np.where(self.work.data.u - self.work.data.l < RHO_TOL)[0]
# ineq_ind = np.setdiff1d(np.setdiff1d(np.arange(self.work.data.m),
# loose_ind), eq_ind)
#
# # Find indices of current constraint types
# old_loose_ind = np.where(self.work.constr_type == -1)
# old_eq_ind = np.where(self.work.constr_type == 1)
# old_ineq_ind = np.where(self.work.constr_type == 0)
#
# # Check if type of any constraint changed
# constr_type_changed = (loose_ind != old_loose_ind).any() or \
# (eq_ind != old_eq_ind).any() or \
# (ineq_ind != old_ineq_ind).any()
#
# # Update type of constraints
# self.work.constr_type[loose_ind] = -1
# self.work.constr_type[eq_ind] = 1
# self.work.constr_type[ineq_ind] = 0
#
# self.work.rho_vec[loose_ind] = RHO_MIN
# self.work.rho_vec[eq_ind] = RHO_MAX
# self.work.rho_vec[ineq_ind] = self.work.settings.rho
#
# self.work.rho_inv_vec = np.reciprocal(self.work.rho_vec)
#
# return constr_type_changed

def print_setup_header(self, data, settings):
"""Print solver header
Expand Down

0 comments on commit 3590088

Please sign in to comment.