1 | """ |
---|
2 | """ |
---|
3 | from scipy import optimize |
---|
4 | |
---|
5 | |
---|
6 | class Parameter(object): |
---|
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 | |
---|
29 | def sasfit(model, pars, x, y, err_y, qmin=None, qmax=None): |
---|
30 | """ |
---|
31 | Fit function |
---|
32 | |
---|
33 | :param model: sas model object |
---|
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() |
---|
92 | chisqr, out, cov = sasfit(line, [cstA, cstB], event.x, y, 0) |
---|
93 | # print "Output parameters:", out |
---|
94 | print "The right answer is [70.0, 1.0]" |
---|
95 | print chisqr, out, cov |
---|