#fit-test.py# #sasio.py# #sasplot.py# #test.py# .#fit-test.py .#sasplot.py .#test.py __init__.py data.DAT fit-test.py fit-test.py~ keypress.py main.py sas.pyc sas.py~ sasfit.pyc sasfit.py~ sasio.pyc sasio.py~ sasplot.pyc sasplot.py~ sastrim.pyc sastrim.py~ sastrim2.py test.py test.py~ @@ -1,4 +1,4 @@ -from numpy import * +import numpy import unittest import scipy.optimize as opt @@ -86,62 +86,6 @@ class SasData(object): out.i[j] = self.i[j] * other return out - -class ModelParameter: - """Parameter class for model fits. - - See scipy.org cookbook entry on least squares fitting - """ - - def __init__(self, value): - self.value = value - - def set(self, value): - self.value = value - - def __call__(self): - return self.value - -def guinier(q): - """Returns the intensities according to Guinier model for a sphere - - Takes a q vector and three parameters (I(0), Rg and background) and - returns calculated values of I for each value of Q. Requires a simple - list of q values as an argument. Relies on i0, Rg, and background having - been previously defined as members of the ModelParameter class. - """ - - assert len(q) != 0, 'Zero length q vector' - - return i0()*exp((-1/3)*(Rg()**2)*(q**2)) + background() - - -def fit(function, parameters, q, i): - """Fit function adapted from the scipy cookbook. - - I don't entirely understand how this works and it is not yet fully - tested but it worked ok on a test rig so will put it here for the - moment. Parameters are taken as a list of instances of the - ModelParameter class. Calling the function therefore looks like e.g. - - >>> fit(guinier, [i0, Rg, background], q_data, i_data) - - It throws a bunch of deprecation warnings but these look to be - quite deep in the library so I'm not going to fuss about them for the - moment. - """ - - def f(params): - j = 0 - for p in parameters: - p.set(params[j]) - j += 1 - return i - function(q) - - p = [param() for param in parameters] - optimize.leastsq(f, p) - - class TestSasData(unittest.TestCase): def setUp(self): @@ -223,17 +167,19 @@ class TestSasData(unittest.TestCase): class TestAnalysis(unittest.TestCase): def test_guinier(self): - test_data = SasData((arange(0,3,0.001)), (arange(4,1,-0.001))) + test_data = SasData((numpy.arange(0,3,0.001)), + (numpy.arange(4,1,-0.001))) # testing with pre-set parameter values - test_guinier = [] test_guinier = guinier(test_data.q) self.assertEqual(len(test_guinier), len(test_data)) self.assertEqual(test_guinier[0], 7) self.assertEqual(test_guinier[-1], 2) + self.assertRaises(AssertionError, guinier, []) + if __name__ == '__main__': i0 = ModelParameter(5) sas.py @@ -4,6 +4,7 @@ from scipy import * from scipy.optimize import * from scipy.io import read_array +import unittest #definitions for specific models @@ -37,3 +38,11 @@ def fit_guinier(data): least_squares_fit = leastsq(guinier_residuals, param_0, args=(i, q)) return least_squares_fit +class TestFit(unittest.TestCase): + + def test_guinier(self): + q = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5] + parameters = [10.0, 100.0, 0.0] + test = guinier(q, parameters) + + self.assertEqual(len(test), len(q)) sasfit.py @@ -1,11 +1,10 @@ # sasio # routines for loading specific data types -from scipy import * -from scipy.io import read_array +from numpy import loadtxt -def loadi22(file): - return read_array(file, lines=(3,-1)) +def loadi22(file): + return loadtxt(file, skiprows=3) sasio.py @@ -1,52 +1,36 @@ # sastrim # routines for trimming and processing scattering curves - -from numpy import * + +import numpy as np +import copy as cp def mask_data(data_to_mask, mask): - """Take an input Q-I data array and filter out points defined by mask + q = data_to_mask[:,0] + i = data_to_mask[:,1] - Takes a 2d Q-I array and generates a new array containing only - the points defined by non-zero values in the array mask. Mask - is a 1d array which is converted internally to a 2d array with - both columns the same. The numpy routine extract is then used - to generate a 2d array with the zero values in mask removed. - """ + q_masked = np.extract(mask, q) + i_masked = np.extract(mask, i) - mask_2d = [mask[:], mask[:]] - return extract(mask_2d, data_to_mask) + masked_data = np.column_stack((q_masked,i_masked)) -def generate_mask(data, mask_ranges) - """"Generates a 1d mask of len(data) from values in mask_range + return masked_data - Takes the values in mask_range which contains values in Q to - remove from the dataset in the form of a 2xn array. Iterates - along the data array determining whether the Q values lie - within the range(s) to be masked. - """ - q_mask = data[:,0] - mask = q_mask[:] # make a slice copy of q list because we will change it - - n = 0 - for a in len(q_mask): - if q_mask[a] > mask_ranges[n,0] and q_mask[a] < mask_ranges[n,1]: - mask[a] = 0 - - if q_mask[a] => mask_ranges[n,1]: - n = n + 1 - - if n > (len(mask_ranges) - 1) break - - if q_mask[a] >mask_ranges[n,0] and q_mask[a] < mask_ranges[n,1]: - mask[a] = 0 - +def generate_mask(data_to_mask, mask_ranges): + q_mask = data_to_mask[:,0] + mask = cp.deepcopy(q_mask) # copy the q list because we will change it + + for i in range(0,len(mask)): + for low, high in mask_ranges: + if q_mask[i] >= low and q_mask[i] <= high: + mask[i] = 0 + break else: - mask[a] = 1 - - else mask[a] = 1 + mask[i] = 1 return mask + + sastrim.py c3eee46318ca5c4525da8967c764e39c8c704166 Cameron Neylon cameron.neylon@stfc.ac.uk http://github.com/cameronneylon/sas/commit/eb5b50c55ce9c56b816f5589497fc5f319c35c08 eb5b50c55ce9c56b816f5589497fc5f319c35c08 2009-03-03T11:55:53-08:00 2009-03-03T11:55:53-08:00 several rounds of trying to produce fitting modules 1baf69687eb6b23c51c89aae8a9f338001872077 Cameron Neylon cameron.neylon@stfc.ac.uk