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