[SciPy-User] Crashing

[Andres Gomez-Lievano]

Daπid <davidmenhur <at> gmail.com> writes:

> 
> 
> 
> On 23 April 2015 at 23:08, Andres Gomez-Lievano <agomez26 <at> 
asu.edu> wrote:
> 
> # Initializing the parameters and data of the problem
> Xc = np.array([80.0, 230.0, 132.0])
> Xp = np.array([74.0, 200.0, 91.0, 77.0])
> size_c = len(Xc)
> size_p = len(Xp)
> 
> 
> Your problem is impossible. You have a 4 by 3 matrix of elements 
between 0 and 1, and are trying to get sums of more than 200. Try first 
with a problem with existing solution:_mat = np.random.random((4,3))Xc = 
_mat.sum(axis=0)Xp = _mat.sum(axis=1)
> 
> Also, Numpy offers functions like reshape, that can replace your 
vec2mat.
> 
> Hard constraints are something that minimisation algorithms don't 
like. It is usually better to add them as a penalty to your target 
function. See here:import numpy as npfrom scipy.optimize import 
minimize# Defining functionsdef entropy(p_vec):    # Log missbehaves if 
there are negative values, so just take them away.    p_vec = 
np.abs(p_vec)    return -1.0 * np.dot(p_vec, np.log(p_vec))def 
const_rows(p_vec, rows, shape):    mat = p_vec.reshape(*shape)    diff = 
mat.sum(axis=0) - rows    return diff.dot(diff)def const_cols(p_vec, 
cols,  shape):    mat = p_vec.reshape(*shape)    diff = mat.sum(axis=1) 
- cols    return diff.dot(diff)def function(x):   # Adjustable 
parameters   lambda_1 = 100.   lambda_2 = 100.   shape = (4, 3)   e = 
entropy(x)   rows_violation = const_rows(x, Xc, shape)   cols_violation 
= const_cols(x, Xp, shape)   return -e + lambda_1 * rows_violation + 
lambda_2 * cols_violation# Initializing the parameters and data of the 
problem_mat = np.random.random((4,3))Xc = _mat.sum(axis=0)Xp = 
_mat.sum(axis=1)size_p = len(Xp)size_c = len(Xc)# Generating the bounds 
for each probability. NB: we are not using them right now!zeroonebounds 
= tuple((0,1) for _ in xrange(size_p * size_c))# Generating the initial 
vector to start the 'minimize' algorithminit_vec = 
np.random.random(len(Xc) * len(Xp))# Finally, running the 
algorithmresult = minimize(function, init_vec)# Display resultsres_mat = 
result.x.reshape(size_p, size_c)print np.abs(res_mat - _mat).sum()# 
Violations of the boundsprint res_mat.sum(axis=1) - Xpprint 
res_mat.sum(axis=0) - Xc
> 
> /David.
> 
>  
> 
> 
> 
> 
> 
> 
> _______________________________________________
> SciPy-User mailing list
> SciPy-User <at> scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
> 


David, thank you for your help.

My problem was well-posed (i.e., is not impossible!). The rows or 
columns of the probability matrix do not have to add up to the vectors 
Xc and Xp that are externally provided.

But in any case, your suggestions helped me solve the problem. As you 
did, I included the constraints into the objective function, and the 
program is now working. THANK YOU!

... Now, other problems have emerged... but those will be the topic of 
another post.

Best,
Andres