Changeset b00b487 in sasview for pr_inversion


Ignore:
Timestamp:
Jun 13, 2008 1:08:35 PM (17 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
b659551
Parents:
96f13a9
Message:

Dealt with error conditions, fixed the uncertainty on the output.

Location:
pr_inversion
Files:
6 added
3 edited

Legend:

Unmodified
Added
Removed

  • pr_inversion/invertor.py

    re39640f rb00b487  
    9292            return self.set_y(value) 
    9393        elif name=='err': 
    94             return self.set_err(value) 
     94            value2 = abs(value) 
     95            return self.set_err(value2) 
    9596        elif name=='d_max': 
    9697            return self.set_dmax(value) 
    … …  
    319320            @param nfunc: number of base functions to use. 
    320321            @param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 
     322 
     323            If the result does not allow us to compute the covariance matrix, 
     324            a matrix filled with zeros will be returned. 
     325 
    321326        """ 
    322327        import math 
    323328        from scipy.linalg.basic import lstsq 
     329         
     330        if self.is_valid()<0: 
     331            raise RuntimeError, "Invertor: invalid data; incompatible data lengths." 
    324332         
    325333        self.nfunc = nfunc 
    … …  
    328336        npts = len(self.x) 
    329337        nq   = nr 
    330         sqrt_alpha = math.sqrt(self.alpha) 
     338        sqrt_alpha = math.sqrt(math.fabs(self.alpha)) 
     339        if sqrt_alpha<0.0: 
     340            nq = 0 
    331341         
    332342        a = numpy.zeros([npts+nq, nfunc]) 
    333343        b = numpy.zeros(npts+nq) 
    334         err = numpy.zeros(nfunc) 
     344        err = numpy.zeros([nfunc, nfunc]) 
    335345         
    336346        for j in range(nfunc): 
    … …  
    350360             
    351361        c, chi2, rank, n = lstsq(a, b) 
     362        # Sanity check 
     363        try: 
     364            float(chi2) 
     365        except: 
     366            chi2 = -1.0 
    352367        self.chi2 = chi2 
    353368                 
    … …  
    358373                inv_cov[i][j] = 0.0 
    359374                for k in range(npts+nr): 
    360                     if self._accept_q(self.x[i]): 
    361                         inv_cov[i][j] = at[i][k]*a[k][j] 
     375                    #if self._accept_q(self.x[i]): 
     376                    inv_cov[i][j] += at[i][k]*a[k][j] 
    362377                     
    363378        # Compute the reg term size for the output 
    … …  
    378393         
    379394        try: 
    380             err = math.fabs(chi2/(npts-nfunc))* inv_cov 
     395            cov = numpy.linalg.pinv(inv_cov) 
     396            err = math.fabs(chi2/float(npts-nfunc)) * cov 
    381397        except: 
    382             print "Error estimating uncertainties" 
     398            # We were not able to estimate the errors, 
     399            # returns an empty covariance matrix 
     400            print "lstsq:", sys.exc_value 
     401            print chi2 
     402            pass 
    383403             
    384404        # Keep a copy of the last output 

  • pr_inversion/test/test_output.txt

    rf71287f4 rb00b487  
    22#nfunc=10 
    33#alpha=0.0007 
    4 #chi2=16654.1 
     4#chi2=836.797 
    55#elapsed=0 
    6 #alpha_estimate=2.23431 
    7 #C_0=1254.55002508+-12.0889464446 
    8 #C_1=386.59602793+-62.4851961458 
    9 #C_2=5.10066421813+-235.84586924 
    10 #C_3=48.454288318+-1217.94687676 
    11 #C_4=0.4174373785+-213577.449972 
    12 #C_5=17.3551976035+-2673.36274352 
    13 #C_6=-0.228001082457+-692.661759936 
    14 #C_7=1.95934734382+-332.134754291 
    15 #C_8=-0.885447495+-200.862052487 
    16 #C_9=0.401866788245+-137.064680609 
     6#alpha_estimate=4.58991e-006 
     7#C_0=217.876478953+-0.0394946442786 
     8#C_1=0.30164617344+-0.407412922728 
     9#C_2=7.15694574584+-2.73726481828 
     10#C_3=-0.79690812244+-5.7183647083 
     11#C_4=0.948333847403+-16.2086432642 
     12#C_5=-0.751612279429+-28.2108893047 
     13#C_6=-0.0796650783501+-52.0568306543 
     14#C_7=-0.518109559543+-88.6068891309 
     15#C_8=-0.129910831984+-141.7178308 
     16#C_9=-0.258134624641+-215.767430098 
    1717<r>  <Pr>  <dPr> 
    18180  0  0 
    19 1.6  236.371  235.678 
    20 3.2  943.103  930.744 
    21 4.8  2113.12  2049.82 
    22 6.4  3734.86  3535.61 
    23 8  5792.64  5311.37 
    24 9.6  8267.09  7284.04 
    25 11.2  11135.7  9348.1 
    26 12.8  14373.5  11389.7 
    27 14.4  17953.5  13291.4 
    28 16  21847.7  14936.2 
    29 17.6  26027.5  16212.7 
    30 19.2  30464.2  17019.1 
    31 20.8  35129.9  17267.8 
    32 22.4  39997.6  16889.2 
    33 24  45041.7  15836.1 
    34 25.6  50238.2  14088.7 
    35 27.2  55565.1  11662.9 
    36 28.8  61002  8635.92 
    37 30.4  66530.4  5281.73 
    38 32  72133.5  3272.64 
    39 33.6  77796  6067.26 
    40 35.2  83503.6  10814.7 
    41 36.8  89243.2  16028.2 
    42 38.4  95001.9  21350.2 
    43 40  100767  26561.2 
    44 41.6  106527  31468 
    45 43.2  112267  35887.3 
    46 44.8  117975  39646.4 
    47 46.4  123636  42586.2 
    48 48  129233  44565.8 
    49 49.6  134751  45467.6 
    50 51.2  140171  45201.5 
    51 52.8  145473  43709.8 
    52 54.4  150638  40971.1 
    53 56  155645  37005.5 
    54 57.6  160473  31882.9 
    55 59.2  165101  25742.1 
    56 60.8  169508  18857.3 
    57 62.4  173675  11970.1 
    58 64  177583  8265.15 
    59 65.6  181215  12750.9 
    60 67.2  184556  21038.7 
    61 68.8  187595  30104.6 
    62 70.4  190321  39131.5 
    63 72  192727  47705.7 
    64 73.6  194808  55513.2 
    65 75.2  196563  62282.7 
    66 76.8  197991  67775.3 
    67 78.4  199093  71784.5 
    68 80  199873  74140.4 
    69 81.6  200335  74714.1 
    70 83.2  200484  73422.4 
    71 84.8  200325  70232.5 
    72 86.4  199862  65166 
    73 88  199101  58304.9 
    74 89.6  198044  49801.3 
    75 91.2  196695  39903.5 
    76 92.8  195056  29050.8 
    77 94.4  193129  18346.5 
    78 96  190915  12397.7 
    79 97.6  188414  18724.8 
    80 99.2  185629  30770.1 
    81 100.8  182559  43834 
    82 102.4  179206  56683.4 
    83 104  175572  68726.9 
    84 105.6  171661  79534.1 
    85 107.2  167475  88745.8 
    86 108.8  163019  96054.4 
    87 110.4  158298  101203 
    88 112  153320  103987 
    89 113.6  148091  104262 
    90 115.2  142620  101947 
    91 116.8  136919  97026.9 
    92 118.4  130998  89560.4 
    93 120  124872  79680.9 
    94 121.6  118558  67605.7 
    95 123.2  112075  53655.8 
    96 124.8  105448  38339.4 
    97 126.4  98705.1  22821 
    98 128  91877.8  13090.4 
    99 129.6  85003.8  22520.4 
    100 131.2  78125.5  39345.4 
    101 132.8  71289.5  56946.4 
    102 134.4  64546.5  73969.2 
    103 136  57949.9  89741.1 
    104 137.6  51554.5  103751 
    105 139.2  45414.9  115564 
    106 140.8  39583.1  124808 
    107 142.4  34106.6  131176 
    108 144  29025.7  134426 
    109 145.6  24371.3  134390 
    110 147.2  20162.8  130981 
    111 148.8  16406  124194 
    112 150.4  13091.6  114114 
    113 152  10194.5  100912 
    114 153.6  7673.45  84850.4 
    115 155.2  5471.36  66272.7 
    116 156.8  3516.48  45601.3 
    117 158.4  1724.23  23326.7 
     191.6  22.981  20.5137 
     203.2  91.986  78.6612 
     214.8  207.187  164.741 
     226.4  368.819  264.129 
     238  577.117  359.427 
     249.6  832.224  433.162 
     2511.2  1134.11  470.89 
     2612.8  1482.46  464.814 
     2714.4  1876.63  419.086 
     2816  2315.57  360.693 
     2917.6  2797.81  353.245 
     3019.2  3321.46  448.357 
     3120.8  3884.25  611.481 
     3222.4  4483.6  786.464 
     3324  5116.72  934.559 
     3425.6  5780.7  1030.78 
     3527.2  6472.66  1062.21 
     3628.8  7189.82  1030.28 
     3730.4  7929.59  955.681 
     3832  8689.7  883.97 
     3933.6  9468.17  878.871 
     4035.2  10263.4  977.553 
     4136.8  11073.9  1152 
     4238.4  11898.7  1342.59 
     4340  12737  1498.3 
     4441.6  13587.7  1587.08 
     4543.2  14450.2  1596.95 
     4644.8  15323.5  1537.99 
     4746.4  16206.5  1445.58 
     4848  17098  1379.39 
     4949.6  17996.3  1401.86 
     5051.2  18899.8  1532.06 
     5152.8  19806.5  1728.86 
     5254.4  20714.3  1927.37 
     5356  21620.8  2073.39 
     5457.6  22524  2135.21 
     5559.2  23421.6  2106.97 
     5660.8  24311.7  2011.28 
     5762.4  25192.3  1900.75 
     5864  26061.8  1848.67 
     5965.6  26918.7  1912.83 
     6067.2  27761.7  2089.96 
     6168.8  28589.6  2319.48 
     6270.4  29401.3  2527.71 
     6372  30195.6  2658.76 
     6473.6  30971.4  2685.17 
     6575.2  31727.2  2611.68 
     6676.8  32461.5  2478.05 
     6778.4  33172.7  2357.19 
     6880  33858.8  2333.81 
     6981.6  34517.8  2453.43 
     7083.2  35147.6  2684.61 
     7184.8  35745.9  2945.13 
     7286.4  36310.5  3152.18 
     7388  36839.5  3249.59 
     7489.6  37330.9  3217.04 
     7591.2  37783.1  3074.61 
     7692.8  38194.6  2886.09 
     7794.4  38564.3  2752.58 
     7896  38891.1  2773.03 
     7997.6  39174.1  2973.02 
     8099.2  39412.6  3281.62 
     81100.8  39605.8  3587.59 
     82102.4  39752.6  3795.95 
     83104  39851.9  3850.56 
     84105.6  39902.4  3741.32 
     85107.2  39902.5  3510.03 
     86108.8  39850.1  3255.58 
     87110.4  39743.2  3120.25 
     88112  39579.6  3218.61 
     89113.6  39356.7  3539.26 
     90115.2  39072.3  3954.92 
     91116.8  38724.1  4317.73 
     92118.4  38310.3  4517.08 
     93120  37829  4494.86 
     94121.6  37278.9  4251.44 
     95123.2  36659  3856.58 
     96124.8  35968.4  3465.78 
     97126.4  35206.6  3306.22 
     98128  34373  3535.97 
     99129.6  33466.9  4074.35 
     100131.2  32487.5  4693.64 
     101132.8  31433.4  5186.14 
     102134.4  30303.1  5411.6 
     103136  29094.1  5295.75 
     104137.6  27803.8  4830.98 
     105139.2  26428.9  4092.02 
     106140.8  24965.9  3287.78 
     107142.4  23411.1  2858.06 
     108144  21761.3  3246.41 
     109145.6  20013.4  4237.67 
     110147.2  18165.4  5345.59 
     111148.8  16216.4  6256.27 
     112150.4  14166.9  6786.26 
     113152  12019  6829.04 
     114153.6  9776.61  6338.91 
     115155.2  7445.34  5326.31 
     116156.8  5032.57  3854 
     117158.4  2547.22  2030.4 

  • pr_inversion/test/utest_invertor.py

    rffca8f2 rb00b487  
    224224        self.invertor.d_max = 160.0 
    225225        # Set a small alpha 
    226         self.invertor.alpha = .0007 
     226        self.invertor.alpha = .005 
    227227        # Set data 
    228228        self.invertor.x   = x 
    … …  
    399399        self.assertEqual(self.invertor.d_max, 160.0) 
    400400        self.assertEqual(self.invertor.alpha, 0.0007) 
    401         self.assertEqual(self.invertor.chi2, 16654.1) 
    402         self.assertAlmostEqual(self.invertor.pr(self.invertor.out, 10.0), 8948.22689927, 4) 
     401        self.assertEqual(self.invertor.chi2, 836.797) 
     402        self.assertAlmostEqual(self.invertor.pr(self.invertor.out, 10.0), 903.31577041, 4) 
    403403         
    404404    def test_qmin(self): 
    … …  
    416416        self.invertor.q_max = None 
    417417        self.assertEqual(self.invertor.q_max, None) 
     418 
     419class TestErrorConditions(unittest.TestCase): 
     420     
     421    def setUp(self): 
     422        self.invertor = Invertor() 
     423        self.invertor.d_max = 100.0 
     424         
     425        # Test array 
     426        self.ntest = 5 
     427        self.x_in = numpy.ones(self.ntest) 
     428        for i in range(self.ntest): 
     429            self.x_in[i] = 1.0*(i+1) 
     430 
     431    def test_negative_errs(self): 
     432        """ 
     433            Test an inversion for which we know the answer 
     434        """ 
     435        x, y, err = load("data_error_1.txt") 
     436 
     437        # Choose the right d_max... 
     438        self.invertor.d_max = 160.0 
     439        # Set a small alpha 
     440        self.invertor.alpha = .0007 
     441        # Set data 
     442        self.invertor.x   = x 
     443        self.invertor.y   = y 
     444        self.invertor.err = err 
     445        # Perform inversion 
     446         
     447        out, cov = self.invertor.lstsq(10) 
     448          
     449    def test_zero_errs(self): 
     450        """ 
     451            Have zero as an error should raise an exception 
     452        """ 
     453        x, y, err = load("data_error_2.txt") 
     454 
     455        # Set data 
     456        self.invertor.x   = x 
     457        self.invertor.y   = y 
     458        self.invertor.err = err 
     459        # Perform inversion 
     460        self.assertRaises(ValueError, self.invertor.invert, 10) 
     461          
     462         
     463    def test_invalid(self): 
     464        """ 
     465            Test an inversion for which we know the answer 
     466        """ 
     467        x, y, err = load("data_error_1.txt") 
     468 
     469        # Set data 
     470        self.invertor.x   = x 
     471        self.invertor.y   = y 
     472        err = numpy.zeros(len(x)-1) 
     473        self.invertor.err = err 
     474        # Perform inversion 
     475        self.assertRaises(RuntimeError, self.invertor.invert, 10) 
     476          
     477         
     478    def test_zero_q(self): 
     479        """ 
     480            One of the q-values is zero. 
     481            An exception should be raised. 
     482        """ 
     483        x, y, err = load("data_error_3.txt") 
     484 
     485        # Set data 
     486        self.assertRaises(ValueError, self.invertor.__setattr__, 'x', x) 
     487         
     488    def test_zero_Iq(self): 
     489        """ 
     490            One of the I(q) points has a value of zero 
     491            Should not complain or crash. 
     492        """ 
     493        x, y, err = load("data_error_4.txt") 
     494 
     495        # Set data 
     496        self.invertor.x   = x 
     497        self.invertor.y   = y 
     498        self.invertor.err = err 
     499        # Perform inversion 
     500        out, cov = self.invertor.lstsq(10) 
     501     
     502    def test_negative_q(self): 
     503        """ 
     504            One q value is negative. 
     505            Makes not sense, but should not complain or crash. 
     506        """ 
     507        x, y, err = load("data_error_5.txt") 
     508 
     509        # Set data 
     510        self.invertor.x   = x 
     511        self.invertor.y   = y 
     512        self.invertor.err = err 
     513        # Perform inversion 
     514        out, cov = self.invertor.lstsq(10) 
     515     
     516    def test_negative_Iq(self): 
     517        """ 
     518            One I(q) value is negative. 
     519            Makes not sense, but should not complain or crash. 
     520        """ 
     521        x, y, err = load("data_error_6.txt") 
     522 
     523        # Set data 
     524        self.invertor.x   = x 
     525        self.invertor.y   = y 
     526        self.invertor.err = err 
     527        # Perform inversion 
     528        out, cov = self.invertor.lstsq(10) 
     529         
     530    def test_nodata(self): 
     531        """ 
     532             No data was loaded. Should not complain. 
     533        """ 
     534        out, cov = self.invertor.lstsq(10) 
     535     
    418536 
    419537         
    … …  
    433551    data_y   = numpy.zeros(0) 
    434552    data_err = numpy.zeros(0) 
     553    scale    = None 
    435554    if not path == None: 
    436555        input_f = open(path,'r') 
    … …  
    442561                x = float(toks[0]) 
    443562                y = float(toks[1]) 
     563                if len(toks)>2: 
     564                    err = float(toks[2]) 
     565                else: 
     566                    if scale==None: 
     567                        scale = 0.15*math.sqrt(y) 
     568                    err = scale*math.sqrt(y) 
    444569                data_x = numpy.append(data_x, x) 
    445570                data_y = numpy.append(data_y, y) 
    446                 # Set the error of the first point to 5% 
    447                 # to make theory look like data 
    448                 scale = 0.1/math.sqrt(data_x[0]) 
    449                 data_err = numpy.append(data_err, scale*math.sqrt(y)) 
     571                data_err = numpy.append(data_err, err) 
    450572            except: 
    451573                pass 
Note: See TracChangeset for help on using the changeset viewer.