Author: bugman Date: Sun Aug 17 19:06:16 2008 New Revision: 7204 URL: http://svn.gna.org/viewcvs/relax?rev=7204&view=rev Log: PCS errors in a file separate from the data can now be read. This matches the changes for the RDC at r7202. Modified: branches/rdc_analysis/generic_fns/pcs.py branches/rdc_analysis/maths_fns/n_state_model.py branches/rdc_analysis/specific_fns/n_state_model.py Modified: branches/rdc_analysis/generic_fns/pcs.py URL: http://svn.gna.org/viewcvs/relax/branches/rdc_analysis/generic_fns/pcs.py?rev=7204&r1=7203&r2=7204&view=diff ============================================================================== --- branches/rdc_analysis/generic_fns/pcs.py (original) +++ branches/rdc_analysis/generic_fns/pcs.py Sun Aug 17 19:06:16 2008 @@ -596,7 +596,7 @@ @type spin_name_col: int or None @param spin_num_col: The column containing the spin number information. @type spin_num_col: int or None - @param sep: The column seperator which, if None, defaults to whitespace. + @param sep: The column separator which, if None, defaults to whitespace. @type sep: str or None """ @@ -611,8 +611,12 @@ cdp = ds[ds.current_pipe] # Test if PCS data corresponding to 'id' already exists. - if hasattr(cdp, 'pcs_ids') and id in cdp.pcs_ids: + if data_col != None and hasattr(cdp, 'pcs_ids') and id in cdp.pcs_ids: raise RelaxPCSError, id + + # Either the data or error column must be supplied. + if data_col == None and error_col == None: + raise RelaxError, "One of either the data or error column must be supplied." # Minimum number of columns. min_col_num = max(mol_name_col, res_num_col, res_name_col, spin_num_col, spin_name_col, data_col, error_col) @@ -626,7 +630,10 @@ header_lines = 0 for i in xrange(len(file_data)): try: - float(file_data[i][data_col]) + if data_col != None: + float(file_data[i][data_col]) + else: + float(file_data[i][error_col]) except: header_lines = header_lines + 1 else: @@ -643,9 +650,9 @@ # Skip missing data. if len(file_data[i]) <= min_col_num: continue - elif file_data[i][data_col] == 'None': + elif data_col != None and file_data[i][data_col] == 'None': continue - elif error_col and file_data[i][error_col] == 'None': + elif error_col != None and file_data[i][error_col] == 'None': continue # Test that the data are numbers. @@ -654,8 +661,9 @@ int(file_data[i][res_num_col]) if spin_num_col != None: int(file_data[i][spin_num_col]) - float(file_data[i][data_col]) - if error_col: + if data_col != None: + float(file_data[i][data_col]) + if error_col != None: float(file_data[i][error_col]) except ValueError: raise RelaxError, "The PCS data in the line " + `file_data[i]` + " is invalid." @@ -685,9 +693,16 @@ id = generate_spin_id_data_array(data=file_data[i], mol_name_col=mol_name_col, res_num_col=res_num_col, res_name_col=res_name_col, spin_num_col=spin_num_col, spin_name_col=spin_name_col) # Convert the data. - value = eval(file_data[i][data_col]) - if error_col: + value = None + if data_col != None: + value = eval(file_data[i][data_col]) + error = None + if error_col != None: error = eval(file_data[i][error_col]) + + # Test the error value (cannot be 0.0). + if error == 0.0: + raise RelaxError, "An invalid error value of zero has been encountered." # Get the corresponding spin container. spin = return_spin(id) @@ -695,12 +710,21 @@ raise RelaxNoSpinError, id # Add the data. - if not hasattr(spin, 'pcs'): - spin.pcs = [] - spin.pcs.append(value) - if error_col: + if data_col != None: + # Initialise. + if not hasattr(spin, 'pcs'): + spin.pcs = [] + + # Append the value. + spin.pcs.append(value) + + # Add the error. + if error_col != None: + # Initialise. if not hasattr(spin, 'pcs_err'): spin.pcs_err = [] + + # Append the error. spin.pcs_err.append(error) Modified: branches/rdc_analysis/maths_fns/n_state_model.py URL: http://svn.gna.org/viewcvs/relax/branches/rdc_analysis/maths_fns/n_state_model.py?rev=7204&r1=7203&r2=7204&view=diff ============================================================================== --- branches/rdc_analysis/maths_fns/n_state_model.py (original) +++ branches/rdc_analysis/maths_fns/n_state_model.py Sun Aug 17 19:06:16 2008 @@ -231,18 +231,30 @@ self.num_align_params = self.num_align * 5 # PCS errors. - if pcs_errors == None: - # Missing errors (the values need to be small, close to ppm units, so the chi-squared value is comparable to the RDC). - self.pcs_sigma_ij = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64) - else: - self.pcs_sigma_ij = pcs_errors + if self.pcs_flag: + err = False + for i in xrange(len(pcs_errors)): + for j in xrange(len(pcs_errors[i])): + if not isNaN(pcs_errors[i, j]): + err = True + if err: + self.pcs_sigma_ij = pcs_errors + else: + # Missing errors (the values need to be small, close to ppm units, so the chi-squared value is comparable to the RDC). + self.pcs_sigma_ij = 0.03 * 1e-6 * ones((self.num_align, self.num_spins), float64) # RDC errors. - if rdc_errors == None: - # Missing errors. - self.rdc_sigma_ij = ones((self.num_align, self.num_spins), float64) - else: - self.rdc_sigma_ij = rdc_errors + if self.rdc_flag: + err = False + for i in xrange(len(rdc_errors)): + for j in xrange(len(rdc_errors[i])): + if not isNaN(rdc_errors[i, j]): + err = True + if err: + self.rdc_sigma_ij = rdc_errors + else: + # Missing errors. + self.rdc_sigma_ij = ones((self.num_align, self.num_spins), float64) # Alignment tensor function and gradient matrices. self.A = zeros((self.num_align, 3, 3), float64) @@ -267,6 +279,9 @@ # Change the NaN to zero. self.Dij[i, j] = 0.0 + # Change the error to one (to avoid zero division). + self.rdc_sigma_ij[i, j] = 1.0 + # Missing data matrices (PCS). if self.pcs_flag: self.missing_deltaij = zeros((self.num_align, self.num_spins), float64) @@ -278,6 +293,9 @@ # Change the NaN to zero. self.deltaij[i, j] = 0.0 + + # Change the error to one (to avoid zero division). + self.pcs_sigma_ij[i, j] = 1.0 # PCS function, gradient, and Hessian matrices. self.deltaij_theta = zeros((self.num_align, self.num_spins), float64) Modified: branches/rdc_analysis/specific_fns/n_state_model.py URL: http://svn.gna.org/viewcvs/relax/branches/rdc_analysis/specific_fns/n_state_model.py?rev=7204&r1=7203&r2=7204&view=diff ============================================================================== --- branches/rdc_analysis/specific_fns/n_state_model.py (original) +++ branches/rdc_analysis/specific_fns/n_state_model.py Sun Aug 17 19:06:16 2008 @@ -402,6 +402,7 @@ # Initialise. pcs = [] + pcs_err = [] unit_vect = [] r = [] pcs_const = [] @@ -418,6 +419,12 @@ # Append the PCSs to the list. pcs.append(spin.pcs) + + # Append the PCS errors (or a list of None). + if hasattr(spin, 'pcs_err'): + pcs_err.append(spin.pcs_err) + else: + pcs_err.append([None]*len(spin.pcs)) # Add empty lists to the r and unit_vector lists. unit_vect.append([]) @@ -473,14 +480,16 @@ # Initialise the numpy objects (the pcs matrix is transposed!). pcs_numpy = zeros((len(pcs[0]), len(pcs)), float64) + pcs_err_numpy = zeros((len(pcs[0]), len(pcs)), float64) unit_vect_numpy = zeros((len(unit_vect), len(unit_vect[0]), 3), float64) # Loop over the spins. for spin_index in xrange(len(pcs)): # Loop over the alignments. for align_index in xrange(len(pcs[spin_index])): - # Transpose and store the PCS value. + # Transpose and store the PCS value and error. pcs_numpy[align_index, spin_index] = pcs[spin_index][align_index] + pcs_err_numpy[align_index, spin_index] = pcs_err[spin_index][align_index] # Loop over the N states. for state_index in xrange(len(unit_vect[spin_index])): @@ -489,9 +498,10 @@ # Convert the PCS from ppm to no units. pcs_numpy = pcs_numpy * 1e-6 + pcs_err_numpy = pcs_err_numpy * 1e-6 # Return the data structures. - return pcs_numpy, unit_vect_numpy, array(pcs_const) + return pcs_numpy, pcs_err_numpy, unit_vect_numpy, array(pcs_const) def __minimise_setup_rdcs(self, param_vector=None, scaling_matrix=None): @@ -510,6 +520,7 @@ # Initialise. rdcs = [] + rdc_err = [] xh_vectors = [] dj = [] @@ -524,6 +535,7 @@ # Add rows of None if other data exists. if hasattr(spin, 'pcs'): rdcs.append([None]*len(cdp.align_tensors)) + rdc_err.append([None]*len(cdp.align_tensors)) xh_vectors.append([None]*3) dj.append(None) @@ -538,6 +550,7 @@ # Add rows of None if other data exists. if hasattr(spin, 'pcs'): rdcs.append([None]*len(cdp.align_tensors)) + rdc_err.append([None]*len(cdp.align_tensors)) xh_vectors.append([None]*3) dj.append(None) @@ -548,6 +561,12 @@ rdcs.append(spin.rdc) xh_vectors.append(spin.xh_vect) + # Append the PCS errors (or a list of None). + if hasattr(spin, 'rdc_err'): + rdc_err.append(spin.rdc_err) + else: + rdc_err.append([None]*len(cdp.align_tensors)) + # Gyromagnetic ratios. gx = return_gyromagnetic_ratio(spin.heteronuc_type) gh = return_gyromagnetic_ratio(spin.proton_type) @@ -557,14 +576,16 @@ # Initialise the numpy objects (the rdc matrix is transposed!). rdcs_numpy = zeros((len(rdcs[0]), len(rdcs)), float64) + rdc_err_numpy = zeros((len(rdcs[0]), len(rdcs)), float64) xh_vect_numpy = zeros((len(xh_vectors), len(xh_vectors[0]), 3), float64) # Loop over the spins. for spin_index in xrange(len(rdcs)): # Loop over the alignments. for align_index in xrange(len(rdcs[spin_index])): - # Transpose and store the RDC value. + # Transpose and store the RDC value and error. rdcs_numpy[align_index, spin_index] = rdcs[spin_index][align_index] + rdc_err_numpy[align_index, spin_index] = rdc_err[spin_index][align_index] # Loop over the N states. for state_index in xrange(len(xh_vectors[spin_index])): @@ -572,7 +593,7 @@ xh_vect_numpy[spin_index, state_index] = xh_vectors[spin_index][state_index] # Return the data structures. - return rdcs_numpy, xh_vect_numpy, array(dj, float64) + return rdcs_numpy, rdc_err_numpy, xh_vect_numpy, array(dj, float64) def __minimise_setup_tensors(self): @@ -1081,17 +1102,17 @@ full_tensors, red_tensor_elem, red_tensor_err, full_in_ref_frame = self.__minimise_setup_tensors() # Get the data structures for optimisation using PCSs as base data sets. - pcs, pcs_vect, pcs_dj = None, None, None + pcs, pcs_err, pcs_vect, pcs_dj = None, None, None, None if 'pcs' in data_types: - pcs, pcs_vect, pcs_dj = self.__minimise_setup_pcs() + pcs, pcs_err, pcs_vect, pcs_dj = self.__minimise_setup_pcs() # Get the data structures for optimisation using RDCs as base data sets. - rdcs, xh_vect, rdc_dj = None, None, None + rdcs, rdc_err, xh_vect, rdc_dj = None, None, None, None if 'rdc' in data_types: - rdcs, xh_vect, rdc_dj = self.__minimise_setup_rdcs() + rdcs, rdc_err, xh_vect, rdc_dj = self.__minimise_setup_rdcs() # Set up the class instance containing the target function. - model = N_state_opt(model=cdp.model, N=cdp.N, init_params=param_vector, full_tensors=full_tensors, red_data=red_tensor_elem, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame, pcs=pcs, rdcs=rdcs, pcs_vect=pcs_vect, xh_vect=xh_vect, pcs_const=pcs_dj, dip_const=rdc_dj, scaling_matrix=scaling_matrix) + model = N_state_opt(model=cdp.model, N=cdp.N, init_params=param_vector, full_tensors=full_tensors, red_data=red_tensor_elem, red_errors=red_tensor_err, full_in_ref_frame=full_in_ref_frame, pcs=pcs, rdcs=rdcs, pcs_errors=pcs_err, rdc_errors=rdc_err, pcs_vect=pcs_vect, xh_vect=xh_vect, pcs_const=pcs_dj, dip_const=rdc_dj, scaling_matrix=scaling_matrix) # Minimisation. if constraints: