[a9d5684] | 1 | """ |
---|
| 2 | """ |
---|
| 3 | from scipy import optimize |
---|
| 4 | |
---|
| 5 | |
---|
| 6 | class Parameter: |
---|
| 7 | """ |
---|
| 8 | Class to handle model parameters |
---|
| 9 | """ |
---|
| 10 | def __init__(self, model, name, value=None): |
---|
| 11 | self.model = model |
---|
| 12 | self.name = name |
---|
| 13 | if not value == None: |
---|
| 14 | self.model.setParam(self.name, value) |
---|
| 15 | |
---|
| 16 | def set(self, value): |
---|
| 17 | """ |
---|
| 18 | Set the value of the parameter |
---|
| 19 | """ |
---|
| 20 | self.model.setParam(self.name, value) |
---|
| 21 | |
---|
| 22 | def __call__(self): |
---|
| 23 | """ |
---|
| 24 | Return the current value of the parameter |
---|
| 25 | """ |
---|
| 26 | return self.model.getParam(self.name) |
---|
| 27 | |
---|
| 28 | |
---|
[b9a5f0e] | 29 | def sasfit(model, pars, x, y, err_y , qmin=None, qmax=None): |
---|
[a9d5684] | 30 | """ |
---|
| 31 | Fit function |
---|
| 32 | |
---|
[79492222] | 33 | :param model: sas model object |
---|
[a9d5684] | 34 | :param pars: list of parameters |
---|
| 35 | :param x: vector of x data |
---|
| 36 | :param y: vector of y data |
---|
| 37 | :param err_y: vector of y errors |
---|
| 38 | """ |
---|
| 39 | def f(params): |
---|
| 40 | """ |
---|
| 41 | Calculates the vector of residuals for each point |
---|
| 42 | in y for a given set of input parameters. |
---|
| 43 | |
---|
| 44 | :param params: list of parameter values |
---|
| 45 | :return: vector of residuals |
---|
| 46 | """ |
---|
| 47 | i = 0 |
---|
| 48 | for p in pars: |
---|
| 49 | p.set(params[i]) |
---|
| 50 | i += 1 |
---|
| 51 | |
---|
| 52 | residuals = [] |
---|
| 53 | for j in range(len(x)): |
---|
| 54 | if x[j] >= qmin and x[j] <= qmax: |
---|
| 55 | residuals.append((y[j] - model.runXY(x[j])) / err_y[j]) |
---|
| 56 | return residuals |
---|
| 57 | |
---|
| 58 | def chi2(params): |
---|
| 59 | """ |
---|
| 60 | Calculates chi^2 |
---|
| 61 | |
---|
| 62 | :param params: list of parameter values |
---|
| 63 | |
---|
| 64 | :return: chi^2 |
---|
| 65 | |
---|
| 66 | """ |
---|
| 67 | sum = 0 |
---|
| 68 | res = f(params) |
---|
| 69 | for item in res: |
---|
| 70 | sum += item * item |
---|
| 71 | return sum |
---|
| 72 | |
---|
| 73 | p = [param() for param in pars] |
---|
| 74 | out, cov_x, info, mesg, success = optimize.leastsq(f, p, full_output=1) |
---|
| 75 | # Calculate chi squared |
---|
| 76 | if len(pars) > 1: |
---|
| 77 | chisqr = chi2(out) |
---|
| 78 | elif len(pars) == 1: |
---|
| 79 | chisqr = chi2([out]) |
---|
| 80 | |
---|
| 81 | return chisqr, out, cov_x |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | def calcCommandline(event): |
---|
| 85 | #Testing implementation |
---|
| 86 | # Fit a Line model |
---|
| 87 | from LineModel import LineModel |
---|
| 88 | line = LineModel() |
---|
| 89 | cstA = Parameter(line, 'A', event.cstA) |
---|
| 90 | cstB = Parameter(line, 'B', event.cstB) |
---|
| 91 | y = line.run() |
---|
[b9a5f0e] | 92 | chisqr, out, cov = sasfit(line, [cstA, cstB], event.x, y, 0) |
---|
[a9d5684] | 93 | # print "Output parameters:", out |
---|
| 94 | print "The right answer is [70.0, 1.0]" |
---|
| 95 | print chisqr, out, cov |
---|