Author: bugman
Date: Tue Apr 15 14:03:51 2014
New Revision: 22766
URL: http://svn.gna.org/viewcvs/relax?rev=22766&view=rev
Log:
Implemented linear constraints for the frame order analysis.
This uses the log-barrier constraint algorithm in the minfx library
https://gna.org/projects/minfx/
to provide constraints without requiring gradients.
Modified:
trunk/specific_analyses/frame_order/api.py
trunk/specific_analyses/frame_order/parameters.py
Modified: trunk/specific_analyses/frame_order/api.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/frame_order/api.py?rev=22766&r1=22765&r2=22766&view=diff
==============================================================================
--- trunk/specific_analyses/frame_order/api.py (original)
+++ trunk/specific_analyses/frame_order/api.py Tue Apr 15 14:03:51 2014
@@ -43,7 +43,7 @@
from specific_analyses.frame_order.data import domain_moving
from specific_analyses.frame_order.optimisation import grid_row,
store_bc_data, target_fn_setup, unpack_opt_results
from specific_analyses.frame_order.parameter_object import Frame_order_params
-from specific_analyses.frame_order.parameters import assemble_param_vector,
assemble_scaling_matrix, param_num, update_model
+from specific_analyses.frame_order.parameters import assemble_param_vector,
assemble_scaling_matrix, linear_constraints, param_num, update_model
class Frame_order(API_base, API_common):
@@ -132,6 +132,17 @@
# Printout.
print("Chi2: %s" % chi2)
+
+
+ def constraint_algorithm(self):
+ """Return the 'Log barrier' optimisation constraint algorithm.
+
+ @return: The 'Log barrier' constraint algorithm.
+ @rtype: str
+ """
+
+ # The log barrier algorithm, as required by minfx.
+ return 'Log barrier'
def create_mc_data(self, data_id=None):
@@ -538,18 +549,10 @@
# Set up the target function for direct calculation.
model, param_vector, scaling_matrix =
target_fn_setup(sim_index=sim_index, verbosity=verbosity, scaling=scaling)
- # Constraints not implemented yet.
+ # Linear constraints.
+ A, b = None, None
if constraints:
- # Turn the constraints off.
- constraints = False
-
- # Remove the method of multipliers arg.
- if not search('^[Gg]rid', min_algor):
- min_algor = min_options[0]
- min_options = min_options[1:]
-
- # Throw a warning.
- warn(RelaxWarning("Constraints are as of yet not implemented -
turning this option off."))
+ A, b = linear_constraints(scaling_matrix=scaling_matrix)
# Grid search.
if search('^[Gg]rid', min_algor):
@@ -557,7 +560,7 @@
# Minimisation.
else:
- results = generic_minimise(func=model.func, args=(),
x0=param_vector, min_algor=min_algor, min_options=min_options,
func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations,
full_output=True, print_flag=verbosity)
+ results = generic_minimise(func=model.func, args=(),
x0=param_vector, min_algor=min_algor, min_options=min_options,
func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, A=A, b=b,
full_output=True, print_flag=verbosity)
# Unpack the results.
unpack_opt_results(results, scaling, scaling_matrix, sim_index)
Modified: trunk/specific_analyses/frame_order/parameters.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/frame_order/parameters.py?rev=22766&r1=22765&r2=22766&view=diff
==============================================================================
--- trunk/specific_analyses/frame_order/parameters.py (original)
+++ trunk/specific_analyses/frame_order/parameters.py Tue Apr 15 14:03:51
2014
@@ -23,7 +23,8 @@
"""Module for handling the frame order model parameters."""
# Python module imports.
-from numpy import array, float64, identity
+from math import pi
+from numpy import array, float64, identity, zeros
# relax module imports.
from specific_analyses.frame_order.data import pivot_fixed, translation_fixed
@@ -86,6 +87,106 @@
return scaling_matrix
+def linear_constraints(scaling_matrix=None):
+ """Create the linear constraint matrices A and b.
+
+ Standard notation
+ =================
+
+ The parameter constraints for the motional amplitude parameters are::
+
+ 0 <= theta <= pi,
+ 0 <= theta_x <= pi,
+ 0 <= theta_y <= pi,
+ -0.125 <= S <= 1,
+ 0 <= sigma_max <= pi,
+
+ The pivot point parameter, domain position parameters, and eigenframe
parameters are unconstrained.
+
+
+ Matrix notation
+ ===============
+
+ In the notation A.x >= b, where A is an matrix of coefficients, x is an
array of parameter values, and b is a vector of scalars, these inequality
constraints are::
+
+ | 1 0 0 0 0 | | 0 |
+ | | | |
+ |-1 0 0 0 0 | | -pi |
+ | | | |
+ | 0 1 0 0 0 | | 0 |
+ | | | theta | | |
+ | 0 -1 0 0 0 | | | | -pi |
+ | | | theta_x | | |
+ | 0 0 1 0 0 | | | | 0 |
+ | | . | theta_y | >= | |
+ | 0 0 -1 0 0 | | | | -pi |
+ | | | S | | |
+ | 0 0 0 1 0 | | | | -0.125 |
+ | | | sigma_max | | |
+ | 0 0 0 -1 0 | | -1 |
+ | | | |
+ | 0 0 0 0 1 | | 0 |
+ | | | |
+ | 0 0 0 0 -1 | | -pi |
+
+
+ @keyword scaling_matrix: The diagonal, square scaling matrix.
+ @type scaling_matrix: numpy rank-2 square matrix
+ @return: The matrices A and b.
+ @rtype: numpy rank-2 NxM array, numpy rank-1 N array
+ """
+
+ # Initialisation (0..j..m).
+ A = []
+ b = []
+ n = param_num()
+ zero_array = zeros(n, float64)
+ i = 0
+ j = 0
+
+ # Loop over the parameters of the model.
+ for i in range(n):
+ # The cone opening angles and sigma_max.
+ if cdp.params[i] in ['cone_theta', 'cone_theta_x', 'cone_theta_y',
'cone_sigma_max']:
+ # 0 <= theta <= pi.
+ A.append(zero_array * 0.0)
+ A.append(zero_array * 0.0)
+ A[j][i] = 1.0
+ A[j+1][i] = -1.0
+ b.append(0.0)
+ b.append(-pi / scaling_matrix[i, i])
+ j = j + 2
+
+ # The pseudo-ellipse restriction (theta_x >= theta_y).
+ if cdp.params[i] == 'cone_theta_x':
+ for m in range(n):
+ if cdp.params[m] == 'cone_theta_y':
+ A.append(zero_array * 0.0)
+ A[j][i] = 1.0
+ A[j][m] = -1.0
+ b.append(0.0)
+ j = j + 1
+
+
+ # The order parameter.
+ if cdp.params[i] == 'cone_s1':
+ # -0.125 <= S <= 1.
+ A.append(zero_array * 0.0)
+ A.append(zero_array * 0.0)
+ A[j][i] = 1.0
+ A[j+1][i] = -1.0
+ b.append(-0.125 / scaling_matrix[i, i])
+ b.append(-1 / scaling_matrix[i, i])
+ j = j + 2
+
+ # Convert to numpy data structures.
+ A = array(A, float64)
+ b = array(b, float64)
+
+ # Return the constraint objects.
+ return A, b
+
+
def param_num():
"""Determine the number of parameters in the model.
r22766 - in /trunk/specific_analyses/frame_order: api.py parameters.py