Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# can probably ignore as handled by constraint
sub_countrol_importance[total_hh_sub_control_index] = 0
LOG_OVERFLOW = -725
log_resid_weights = np.log(np.maximum(sub_resid_weights, np.exp(LOG_OVERFLOW))).flatten('F')
assert not np.isnan(log_resid_weights).any()
log_parent_resid_weights = \
np.log(np.maximum(parent_resid_weights, np.exp(LOG_OVERFLOW))).flatten('F')
assert not np.isnan(log_parent_resid_weights).any()
# subzone and parent objective and relaxation penalties
# note: cvxpy overloads * so * in following is matrix multiplication
objective = cvx.Maximize(
cvx.sum_entries(cvx.mul_elemwise(log_resid_weights, cvx.vec(x))) +
cvx.sum_entries(cvx.mul_elemwise(log_parent_resid_weights, cvx.vec(cvx.sum_entries(x, axis=0)))) - # nopep8
cvx.sum_entries(relax_le * sub_countrol_importance) -
cvx.sum_entries(relax_ge * sub_countrol_importance) -
cvx.sum_entries(cvx.mul_elemwise(parent_countrol_importance, parent_relax_le)) -
cvx.sum_entries(cvx.mul_elemwise(parent_countrol_importance, parent_relax_ge))
)
constraints = [
(x * sub_incidence) - relax_le >= 0,
(x * sub_incidence) - relax_le <= lp_right_hand_side,
(x * sub_incidence) + relax_ge >= lp_right_hand_side,
(x * sub_incidence) + relax_ge <= hh_constraint_ge_bound,
x >= 0.0,
x <= x_max,
relax_le >= 0.0,
def Hawkes_log_lik(T, alpha_opt, lambda_opt, lambda_ti, survival, for_cvx=False):
L = 0
for i in range(len(lambda_ti)):
if for_cvx and len(lambda_ti) > 0:
L += CVX.sum_entries(CVX.log(lambda_opt + alpha_opt * lambda_ti[i]))
else:
L += np.sum(np.log(lambda_opt + alpha_opt * lambda_ti[i]))
L -= lambda_opt * T[i] + alpha_opt * survival[i]
return L
# Get residuals in new marginals from truncating to int
A_residuals = np.dot(x, hh_table) - np.dot(x_int, hh_table)
x_residuals = x - x_int
# Coefficients in objective function
x_log = np.log(x_residuals)
# Decision variables for optimization
y = cvx.Variable(n_tracts, n_samples)
# Relaxation factors
U = cvx.Variable(n_tracts, n_controls)
V = cvx.Variable(n_tracts, n_controls)
objective = cvx.Maximize(
cvx.sum_entries(
cvx.sum_entries(cvx.mul_elemwise(x_log, y), axis=1) -
(gamma) * cvx.sum_entries(U, axis=1) -
(gamma) * cvx.sum_entries(V, axis=1)
)
)
constraints = [
y * hh_table <= A_residuals + U,
y * hh_table >= A_residuals - V,
U >= 0,
V >= 0,
y >= 0,
y <= 1.0
]
prob = cvx.Problem(objective, constraints)
sample_count, control_count = incidence.shape
# - Decision variables for optimization
x = cvx.Variable(1, sample_count)
# - Create positive continuous constraint relaxation variables
relax_le = cvx.Variable(control_count)
relax_ge = cvx.Variable(control_count)
# FIXME - could ignore as handled by constraint?
control_importance_weights[total_hh_control_index] = 0
# - Set objective
objective = cvx.Maximize(
cvx.sum_entries(cvx.mul_elemwise(log_resid_weights, cvx.vec(x))) -
cvx.sum_entries(cvx.mul_elemwise(control_importance_weights, relax_le)) -
cvx.sum_entries(cvx.mul_elemwise(control_importance_weights, relax_ge))
)
total_hh_constraint = lp_right_hand_side[total_hh_control_index]
# 1.0 unless resid_weights is zero
max_x = (~(resid_weights == 0.0)).astype(float).reshape((1, -1))
constraints = [
# - inequality constraints
cvx.vec(x * incidence) - relax_le >= 0,
cvx.vec(x * incidence) - relax_le <= lp_right_hand_side,
cvx.vec(x * incidence) + relax_ge >= lp_right_hand_side,
cvx.vec(x * incidence) + relax_ge <= hh_constraint_ge_bound,
constraints = [
x >= 0,
x.T * hh_table == A,
]
prob = cvx.Problem(objective, constraints)
prob.solve(solver=cvx.SCS, verbose=verbose_solver)
return x.value
else:
# With relaxation factors
z = cvx.Variable(n_controls)
objective = cvx.Maximize(
cvx.sum_entries(cvx.entr(x) + cvx.mul_elemwise(cvx.log(w.T), x)) +
cvx.sum_entries(mu * (cvx.entr(z)))
)
constraints = [
x >= 0,
z >= 0,
x.T * hh_table == cvx.mul_elemwise(A, z.T),
]
prob = cvx.Problem(objective, constraints)
prob.solve(solver=cvx.SCS, verbose=verbose_solver)
return x.value, z.value
def solve_relaxation(prob, *args, **kwargs):
"""Solve the SDP relaxation.
"""
# lifted variables and semidefinite constraint
X = cvx.Semidef(prob.n + 1)
W = prob.f0.homogeneous_form()
rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))
rel_constr = [X[-1, -1] == 1]
for f in prob.fs:
W = f.homogeneous_form()
lhs = cvx.sum_entries(cvx.mul_elemwise(W, X))
if f.relop == '==':
rel_constr.append(lhs == 0)
else:
rel_constr.append(lhs <= 0)
rel_prob = cvx.Problem(rel_obj, rel_constr)
rel_prob.solve(*args, **kwargs)
if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
raise Exception("Relaxation problem status: %s" % rel_prob.status)
"""
cvxopt.solvers.lp Solves a pair of primal and dual LPs
minimize c'*x
subject to G*x + s = h
A*x = b
s >= 0
maximize -h'*z - b'*y
subject to G'*z + A'*y + c = 0
z >= 0.
cvxopt.solvers.lp(c, G, h[, A, b[, solver[, primalstart[, dualstart]]]])
"""
obj = cvx.Minimize(cvx.sum_entries(g_ikl * alpha_var))
prob = cvx.Problem(obj, cons)
prob.solve()
alpha_new = np.array(alpha_var.value).T[0]
""" now the problem collapse into a really simple qp problem """
# compute beta_new
for i in range(sample_num):
alpha_range = class_info.get_range(i)
alpha_piece = alpha_new[alpha_range[0]:alpha_range[1]]
c_list = c[i]
for k in range(class_num):
beta_new[k][i] = np.inner(c_list[k], alpha_piece)
"""
We need to find lambda which will make W(alpha + lambda*alpha_new) be maximum