[9e8dc22] | 1 | from sans.pr.core.pr_inversion import Cinvertor |
---|
| 2 | import numpy |
---|
| 3 | |
---|
| 4 | class Invertor(Cinvertor): |
---|
| 5 | |
---|
[eca05c8] | 6 | ## Chisqr of the last computation |
---|
| 7 | chisqr = 0 |
---|
| 8 | |
---|
[9e8dc22] | 9 | def __init__(self): |
---|
| 10 | Cinvertor.__init__(self) |
---|
| 11 | |
---|
| 12 | def __setattr__(self, name, value): |
---|
| 13 | """ |
---|
| 14 | Set the value of an attribute. |
---|
| 15 | Access the parent class methods for |
---|
| 16 | x, y, err and d_max. |
---|
| 17 | """ |
---|
| 18 | if name=='x': |
---|
[eca05c8] | 19 | if 0.0 in value: |
---|
| 20 | raise ValueError, "Invertor: one of your q-values is zero. Delete that entry before proceeding" |
---|
[9e8dc22] | 21 | return self.set_x(value) |
---|
| 22 | elif name=='y': |
---|
| 23 | return self.set_y(value) |
---|
| 24 | elif name=='err': |
---|
| 25 | return self.set_err(value) |
---|
| 26 | elif name=='d_max': |
---|
| 27 | return self.set_dmax(value) |
---|
[eca05c8] | 28 | elif name=='alpha': |
---|
| 29 | return self.set_alpha(value) |
---|
[9e8dc22] | 30 | |
---|
| 31 | return Cinvertor.__setattr__(self, name, value) |
---|
| 32 | |
---|
| 33 | def __getattr__(self, name): |
---|
| 34 | """ |
---|
| 35 | Return the value of an attribute |
---|
| 36 | For the moment x, y, err and d_max are write-only |
---|
| 37 | TODO: change that! |
---|
| 38 | """ |
---|
| 39 | import numpy |
---|
| 40 | if name=='x': |
---|
| 41 | out = numpy.ones(self.get_nx()) |
---|
| 42 | self.get_x(out) |
---|
| 43 | return out |
---|
| 44 | elif name=='y': |
---|
| 45 | out = numpy.ones(self.get_ny()) |
---|
| 46 | self.get_y(out) |
---|
| 47 | return out |
---|
| 48 | elif name=='err': |
---|
| 49 | out = numpy.ones(self.get_nerr()) |
---|
| 50 | self.get_err(out) |
---|
| 51 | return out |
---|
| 52 | elif name=='d_max': |
---|
| 53 | return self.get_dmax() |
---|
[eca05c8] | 54 | elif name=='alpha': |
---|
| 55 | return self.get_alpha() |
---|
[9e8dc22] | 56 | elif name in self.__dict__: |
---|
| 57 | return self.__dict__[name] |
---|
| 58 | return None |
---|
| 59 | |
---|
| 60 | def invert(self, nfunc=5): |
---|
| 61 | """ |
---|
| 62 | Perform inversion to P(r) |
---|
| 63 | """ |
---|
| 64 | from scipy import optimize |
---|
| 65 | |
---|
| 66 | # First, check that the current data is valid |
---|
| 67 | if self.is_valid()<=0: |
---|
| 68 | raise RuntimeError, "Invertor.invert: Data array are of different length" |
---|
| 69 | |
---|
| 70 | p = numpy.ones(nfunc) |
---|
| 71 | out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True) |
---|
| 72 | |
---|
[eca05c8] | 73 | # Compute chi^2 |
---|
| 74 | res = self.residuals(out) |
---|
| 75 | chisqr = 0 |
---|
| 76 | for i in range(len(res)): |
---|
| 77 | chisqr += res[i] |
---|
| 78 | |
---|
| 79 | self.chi2 = chisqr |
---|
| 80 | |
---|
| 81 | return out, cov_x |
---|
| 82 | |
---|
| 83 | def pr_fit(self, nfunc=5): |
---|
| 84 | """ |
---|
| 85 | Perform inversion to P(r) |
---|
| 86 | """ |
---|
| 87 | from scipy import optimize |
---|
| 88 | |
---|
| 89 | # First, check that the current data is valid |
---|
| 90 | if self.is_valid()<=0: |
---|
| 91 | raise RuntimeError, "Invertor.invert: Data arrays are of different length" |
---|
| 92 | |
---|
| 93 | p = numpy.ones(nfunc) |
---|
| 94 | out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True) |
---|
| 95 | |
---|
| 96 | # Compute chi^2 |
---|
| 97 | res = self.pr_residuals(out) |
---|
| 98 | chisqr = 0 |
---|
| 99 | for i in range(len(res)): |
---|
| 100 | chisqr += res[i] |
---|
| 101 | |
---|
| 102 | self.chisqr = chisqr |
---|
| 103 | |
---|
[9e8dc22] | 104 | return out, cov_x |
---|
| 105 | |
---|
[eca05c8] | 106 | def pr_err(self, c, c_cov, r): |
---|
| 107 | import math |
---|
| 108 | c_err = numpy.zeros(len(c)) |
---|
| 109 | for i in range(len(c)): |
---|
| 110 | try: |
---|
| 111 | c_err[i] = math.sqrt(math.fabs(c_cov[i][i])) |
---|
| 112 | except: |
---|
| 113 | import sys |
---|
| 114 | print sys.exc_value |
---|
| 115 | print "oups", c_cov[i][i] |
---|
| 116 | c_err[i] = c[i] |
---|
| 117 | |
---|
| 118 | return self.get_pr_err(c, c_err, r) |
---|
| 119 | |
---|
| 120 | |
---|
[9e8dc22] | 121 | if __name__ == "__main__": |
---|
| 122 | o = Invertor() |
---|
| 123 | |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | |
---|