source: sasview/_modules/sas/pr/invertor.html @ a462c6a

gh-pages
Last change on this file since a462c6a was a462c6a, checked in by ajj, 9 years ago

Rebuild to fix index and modules docs
  • Property mode set to 100644
File size: 100.4 KB
Line 
1<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
2  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
3
4
5<html xmlns="http://www.w3.org/1999/xhtml">
6  <head>
7    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
8   
9    <title>sas.pr.invertor &mdash; SasView 3.0.0 documentation</title>
10   
11    <link rel="stylesheet" href="../../../_static/default.css" type="text/css" />
12    <link rel="stylesheet" href="../../../_static/pygments.css" type="text/css" />
13   
14    <script type="text/javascript">
15      var DOCUMENTATION_OPTIONS = {
16        URL_ROOT:    '../../../',
17        VERSION:     '3.0.0',
18        COLLAPSE_INDEX: false,
19        FILE_SUFFIX: '.html',
20        HAS_SOURCE:  true
21      };
22    </script>
23    <script type="text/javascript" src="../../../_static/jquery.js"></script>
24    <script type="text/javascript" src="../../../_static/underscore.js"></script>
25    <script type="text/javascript" src="../../../_static/doctools.js"></script>
26    <script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
27    <link rel="top" title="SasView 3.0.0 documentation" href="../../../index.html" />
28    <link rel="up" title="Module code" href="../../index.html" /> 
29  </head>
30  <body>
31    <div class="related">
32      <h3>Navigation</h3>
33      <ul>
34        <li class="right" style="margin-right: 10px">
35          <a href="../../../genindex.html" title="General Index"
36             accesskey="I">index</a></li>
37        <li class="right" >
38          <a href="../../../py-modindex.html" title="Python Module Index"
39             >modules</a> |</li>
40        <li><a href="../../../index.html">SasView 3.0.0 documentation</a> &raquo;</li>
41          <li><a href="../../index.html" accesskey="U">Module code</a> &raquo;</li> 
42      </ul>
43    </div> 
44
45    <div class="document">
46      <div class="documentwrapper">
47        <div class="bodywrapper">
48          <div class="body">
49           
50  <h1>Source code for sas.pr.invertor</h1><div class="highlight"><pre>
51<span class="sd">&quot;&quot;&quot;</span>
52<span class="sd">Module to perform P(r) inversion.</span>
53<span class="sd">The module contains the Invertor class.</span>
54<span class="sd">&quot;&quot;&quot;</span>
55
56<span class="kn">import</span> <span class="nn">numpy</span>
57<span class="kn">import</span> <span class="nn">sys</span>
58<span class="kn">import</span> <span class="nn">math</span>
59<span class="kn">import</span> <span class="nn">time</span>
60<span class="kn">import</span> <span class="nn">copy</span>
61<span class="kn">import</span> <span class="nn">os</span>
62<span class="kn">import</span> <span class="nn">re</span>
63<span class="kn">from</span> <span class="nn">numpy.linalg</span> <span class="kn">import</span> <span class="n">lstsq</span>
64<span class="kn">from</span> <span class="nn">scipy</span> <span class="kn">import</span> <span class="n">optimize</span>
65<span class="kn">from</span> <span class="nn">sas.pr.core.pr_inversion</span> <span class="kn">import</span> <span class="n">Cinvertor</span>
66
67<div class="viewcode-block" id="help"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.help">[docs]</a><span class="k">def</span> <span class="nf">help</span><span class="p">():</span>
68    <span class="sd">&quot;&quot;&quot;</span>
69<span class="sd">    Provide general online help text</span>
70<span class="sd">    Future work: extend this function to allow topic selection</span>
71<span class="sd">    &quot;&quot;&quot;</span>
72    <span class="n">info_txt</span>  <span class="o">=</span> <span class="s">&quot;The inversion approach is based on Moore, J. Appl. Cryst. &quot;</span>
73    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;(1980) 13, 168-175.</span><span class="se">\n\n</span><span class="s">&quot;</span>
74    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;P(r) is set to be equal to an expansion of base functions &quot;</span>
75    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;of the type &quot;</span>
76    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;phi_n(r) = 2*r*sin(pi*n*r/D_max). The coefficient of each &quot;</span>
77    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;base functions &quot;</span>
78    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;in the expansion is found by performing a least square fit &quot;</span>
79    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;with the &quot;</span>
80    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;following fit function:</span><span class="se">\n\n</span><span class="s">&quot;</span>
81    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;chi**2 = sum_i[ I_meas(q_i) - I_th(q_i) ]**2/error**2 +&quot;</span>
82    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;Reg_term</span><span class="se">\n\n</span><span class="s">&quot;</span>
83    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;where I_meas(q) is the measured scattering intensity and &quot;</span>
84    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;I_th(q) is &quot;</span>
85    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;the prediction from the Fourier transform of the P(r) &quot;</span>
86    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;expansion. &quot;</span>
87    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;The Reg_term term is a regularization term set to the second&quot;</span>
88    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot; derivative &quot;</span>
89    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;d**2P(r)/dr**2 integrated over r. It is used to produce &quot;</span>
90    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;a smooth P(r) output.</span><span class="se">\n\n</span><span class="s">&quot;</span>
91    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;The following are user inputs:</span><span class="se">\n\n</span><span class="s">&quot;</span>
92    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;   - Number of terms: the number of base functions in the P(r)&quot;</span>
93    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot; expansion.</span><span class="se">\n\n</span><span class="s">&quot;</span>
94    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;   - Regularization constant: a multiplicative constant &quot;</span>
95    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;to set the size of &quot;</span>
96    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;the regularization term.</span><span class="se">\n\n</span><span class="s">&quot;</span>
97    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;   - Maximum distance: the maximum distance between any &quot;</span>
98    <span class="n">info_txt</span> <span class="o">+=</span> <span class="s">&quot;two points in the system.</span><span class="se">\n</span><span class="s">&quot;</span>
99     
100    <span class="k">return</span> <span class="n">info_txt</span>
101   
102</div>
103<div class="viewcode-block" id="Invertor"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor">[docs]</a><span class="k">class</span> <span class="nc">Invertor</span><span class="p">(</span><span class="n">Cinvertor</span><span class="p">):</span>
104    <span class="sd">&quot;&quot;&quot;</span>
105<span class="sd">    Invertor class to perform P(r) inversion</span>
106<span class="sd">    </span>
107<span class="sd">    The problem is solved by posing the problem as  Ax = b,</span>
108<span class="sd">    where x is the set of coefficients we are looking for.</span>
109<span class="sd">    </span>
110<span class="sd">    Npts is the number of points.</span>
111<span class="sd">    </span>
112<span class="sd">    In the following i refers to the ith base function coefficient.</span>
113<span class="sd">    The matrix has its entries j in its first Npts rows set to ::</span>
114
115<span class="sd">        A[j][i] = (Fourier transformed base function for point j)</span>
116<span class="sd">        </span>
117<span class="sd">    We them choose a number of r-points, n_r, to evaluate the second</span>
118<span class="sd">    derivative of P(r) at. This is used as our regularization term.</span>
119<span class="sd">    For a vector r of length n_r, the following n_r rows are set to ::</span>
120
121<span class="sd">        A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,</span>
122<span class="sd">        evaluated at r[j])</span>
123<span class="sd">        </span>
124<span class="sd">    The vector b has its first Npts entries set to ::</span>
125
126<span class="sd">        b[j] = (I(q) observed for point j)</span>
127<span class="sd">        </span>
128<span class="sd">    The following n_r entries are set to zero.</span>
129<span class="sd">    </span>
130<span class="sd">    The result is found by using scipy.linalg.basic.lstsq to invert</span>
131<span class="sd">    the matrix and find the coefficients x.</span>
132<span class="sd">    </span>
133<span class="sd">    Methods inherited from Cinvertor:</span>
134
135<span class="sd">    * ``get_peaks(pars)``: returns the number of P(r) peaks</span>
136<span class="sd">    * ``oscillations(pars)``: returns the oscillation parameters for the output P(r)</span>
137<span class="sd">    * ``get_positive(pars)``: returns the fraction of P(r) that is above zero</span>
138<span class="sd">    * ``get_pos_err(pars)``: returns the fraction of P(r) that is 1-sigma above zero</span>
139<span class="sd">    &quot;&quot;&quot;</span>
140    <span class="c">## Chisqr of the last computation</span>
141    <span class="n">chi2</span>  <span class="o">=</span> <span class="mi">0</span>
142    <span class="c">## Time elapsed for last computation</span>
143    <span class="n">elapsed</span> <span class="o">=</span> <span class="mi">0</span>
144    <span class="c">## Alpha to get the reg term the same size as the signal</span>
145    <span class="n">suggested_alpha</span> <span class="o">=</span> <span class="mi">0</span>
146    <span class="c">## Last number of base functions used</span>
147    <span class="n">nfunc</span> <span class="o">=</span> <span class="mi">10</span>
148    <span class="c">## Last output values</span>
149    <span class="n">out</span> <span class="o">=</span> <span class="bp">None</span>
150    <span class="c">## Last errors on output values</span>
151    <span class="n">cov</span> <span class="o">=</span> <span class="bp">None</span>
152    <span class="c">## Background value</span>
153    <span class="n">background</span> <span class="o">=</span> <span class="mi">0</span>
154    <span class="c">## Information dictionary for application use</span>
155    <span class="n">info</span> <span class="o">=</span> <span class="p">{}</span>
156   
157    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
158        <span class="n">Cinvertor</span><span class="o">.</span><span class="n">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span>
159       
160    <span class="k">def</span> <span class="nf">__setstate__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">state</span><span class="p">):</span>
161        <span class="sd">&quot;&quot;&quot;</span>
162<span class="sd">        restore the state of invertor for pickle</span>
163<span class="sd">        &quot;&quot;&quot;</span>
164        <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">__dict__</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">alpha</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span><span class="p">,</span>
165         <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span><span class="p">,</span>
166         <span class="bp">self</span><span class="o">.</span><span class="n">x</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">y</span><span class="p">,</span>
167         <span class="bp">self</span><span class="o">.</span><span class="n">err</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span><span class="p">,</span>
168         <span class="bp">self</span><span class="o">.</span><span class="n">slit_height</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">slit_width</span><span class="p">)</span> <span class="o">=</span> <span class="n">state</span>
169   
170    <span class="k">def</span> <span class="nf">__reduce_ex__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">proto</span><span class="p">):</span>
171        <span class="sd">&quot;&quot;&quot;</span>
172<span class="sd">        Overwrite the __reduce_ex__</span>
173<span class="sd">        &quot;&quot;&quot;</span>
174
175        <span class="n">state</span> <span class="o">=</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">__dict__</span><span class="p">,</span>
176                 <span class="bp">self</span><span class="o">.</span><span class="n">alpha</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span><span class="p">,</span>
177                 <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span><span class="p">,</span>
178                 <span class="bp">self</span><span class="o">.</span><span class="n">x</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">y</span><span class="p">,</span>
179                 <span class="bp">self</span><span class="o">.</span><span class="n">err</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span><span class="p">,</span>
180                 <span class="bp">self</span><span class="o">.</span><span class="n">slit_height</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">slit_width</span><span class="p">,</span>
181                 <span class="p">)</span>
182        <span class="k">return</span> <span class="p">(</span><span class="n">Invertor</span><span class="p">,</span> <span class="nb">tuple</span><span class="p">(),</span> <span class="n">state</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="bp">None</span><span class="p">)</span>
183   
184    <span class="k">def</span> <span class="nf">__setattr__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">name</span><span class="p">,</span> <span class="n">value</span><span class="p">):</span>
185        <span class="sd">&quot;&quot;&quot;</span>
186<span class="sd">        Set the value of an attribute.</span>
187<span class="sd">        Access the parent class methods for</span>
188<span class="sd">        x, y, err, d_max, q_min, q_max and alpha</span>
189<span class="sd">        &quot;&quot;&quot;</span>
190        <span class="k">if</span>   <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;x&#39;</span><span class="p">:</span>
191            <span class="k">if</span> <span class="mf">0.0</span> <span class="ow">in</span> <span class="n">value</span><span class="p">:</span>
192                <span class="n">msg</span> <span class="o">=</span> <span class="s">&quot;Invertor: one of your q-values is zero. &quot;</span>
193                <span class="n">msg</span> <span class="o">+=</span> <span class="s">&quot;Delete that entry before proceeding&quot;</span>
194                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">,</span> <span class="n">msg</span>
195            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_x</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
196        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;y&#39;</span><span class="p">:</span>
197            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_y</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
198        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;err&#39;</span><span class="p">:</span>
199            <span class="n">value2</span> <span class="o">=</span> <span class="nb">abs</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
200            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_err</span><span class="p">(</span><span class="n">value2</span><span class="p">)</span>
201        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;d_max&#39;</span><span class="p">:</span>
202            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_dmax</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
203        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;q_min&#39;</span><span class="p">:</span>
204            <span class="k">if</span> <span class="n">value</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
205                <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_qmin</span><span class="p">(</span><span class="o">-</span><span class="mf">1.0</span><span class="p">)</span>
206            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_qmin</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
207        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;q_max&#39;</span><span class="p">:</span>
208            <span class="k">if</span> <span class="n">value</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
209                <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_qmax</span><span class="p">(</span><span class="o">-</span><span class="mf">1.0</span><span class="p">)</span>
210            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_qmax</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
211        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;alpha&#39;</span><span class="p">:</span>
212            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_alpha</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
213        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;slit_height&#39;</span><span class="p">:</span>
214            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_slit_height</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
215        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;slit_width&#39;</span><span class="p">:</span>
216            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_slit_width</span><span class="p">(</span><span class="n">value</span><span class="p">)</span>
217        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;has_bck&#39;</span><span class="p">:</span>
218            <span class="k">if</span> <span class="n">value</span> <span class="o">==</span> <span class="bp">True</span><span class="p">:</span>
219                <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_has_bck</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
220            <span class="k">elif</span> <span class="n">value</span> <span class="o">==</span> <span class="bp">False</span><span class="p">:</span>
221                <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">set_has_bck</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
222            <span class="k">else</span><span class="p">:</span>
223                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">,</span> <span class="s">&quot;Invertor: has_bck can only be True or False&quot;</span>
224           
225        <span class="k">return</span> <span class="n">Cinvertor</span><span class="o">.</span><span class="n">__setattr__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">name</span><span class="p">,</span> <span class="n">value</span><span class="p">)</span>
226   
227    <span class="k">def</span> <span class="nf">__getattr__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">name</span><span class="p">):</span>
228        <span class="sd">&quot;&quot;&quot;</span>
229<span class="sd">        Return the value of an attribute</span>
230<span class="sd">        &quot;&quot;&quot;</span>
231        <span class="c">#import numpy</span>
232        <span class="k">if</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;x&#39;</span><span class="p">:</span>
233            <span class="n">out</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">get_nx</span><span class="p">())</span>
234            <span class="bp">self</span><span class="o">.</span><span class="n">get_x</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
235            <span class="k">return</span> <span class="n">out</span>
236        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;y&#39;</span><span class="p">:</span>
237            <span class="n">out</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">get_ny</span><span class="p">())</span>
238            <span class="bp">self</span><span class="o">.</span><span class="n">get_y</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
239            <span class="k">return</span> <span class="n">out</span>
240        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;err&#39;</span><span class="p">:</span>
241            <span class="n">out</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">get_nerr</span><span class="p">())</span>
242            <span class="bp">self</span><span class="o">.</span><span class="n">get_err</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
243            <span class="k">return</span> <span class="n">out</span>
244        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;d_max&#39;</span><span class="p">:</span>
245            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_dmax</span><span class="p">()</span>
246        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;q_min&#39;</span><span class="p">:</span>
247            <span class="n">qmin</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_qmin</span><span class="p">()</span>
248            <span class="k">if</span> <span class="n">qmin</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span>
249                <span class="k">return</span> <span class="bp">None</span>
250            <span class="k">return</span> <span class="n">qmin</span>
251        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;q_max&#39;</span><span class="p">:</span>
252            <span class="n">qmax</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_qmax</span><span class="p">()</span>
253            <span class="k">if</span> <span class="n">qmax</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span>
254                <span class="k">return</span> <span class="bp">None</span>
255            <span class="k">return</span> <span class="n">qmax</span>
256        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;alpha&#39;</span><span class="p">:</span>
257            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_alpha</span><span class="p">()</span>
258        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;slit_height&#39;</span><span class="p">:</span>
259            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_slit_height</span><span class="p">()</span>
260        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;slit_width&#39;</span><span class="p">:</span>
261            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_slit_width</span><span class="p">()</span>
262        <span class="k">elif</span> <span class="n">name</span> <span class="o">==</span> <span class="s">&#39;has_bck&#39;</span><span class="p">:</span>
263            <span class="n">value</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_has_bck</span><span class="p">()</span>
264            <span class="k">if</span> <span class="n">value</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
265                <span class="k">return</span> <span class="bp">True</span>
266            <span class="k">else</span><span class="p">:</span>
267                <span class="k">return</span> <span class="bp">False</span>
268        <span class="k">elif</span> <span class="n">name</span> <span class="ow">in</span> <span class="bp">self</span><span class="o">.</span><span class="n">__dict__</span><span class="p">:</span>
269            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">__dict__</span><span class="p">[</span><span class="n">name</span><span class="p">]</span>
270        <span class="k">return</span> <span class="bp">None</span>
271   
272<div class="viewcode-block" id="Invertor.clone"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.clone">[docs]</a>    <span class="k">def</span> <span class="nf">clone</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
273        <span class="sd">&quot;&quot;&quot;</span>
274<span class="sd">        Return a clone of this instance</span>
275<span class="sd">        &quot;&quot;&quot;</span>
276        <span class="c">#import copy</span>
277       
278        <span class="n">invertor</span> <span class="o">=</span> <span class="n">Invertor</span><span class="p">()</span>
279        <span class="n">invertor</span><span class="o">.</span><span class="n">chi2</span>    <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">chi2</span>
280        <span class="n">invertor</span><span class="o">.</span><span class="n">elapsed</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">elapsed</span>
281        <span class="n">invertor</span><span class="o">.</span><span class="n">nfunc</span>   <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span>
282        <span class="n">invertor</span><span class="o">.</span><span class="n">alpha</span>   <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">alpha</span>
283        <span class="n">invertor</span><span class="o">.</span><span class="n">d_max</span>   <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span>
284        <span class="n">invertor</span><span class="o">.</span><span class="n">q_min</span>   <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span>
285        <span class="n">invertor</span><span class="o">.</span><span class="n">q_max</span>   <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span>
286       
287        <span class="n">invertor</span><span class="o">.</span><span class="n">x</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">x</span>
288        <span class="n">invertor</span><span class="o">.</span><span class="n">y</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">y</span>
289        <span class="n">invertor</span><span class="o">.</span><span class="n">err</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">err</span>
290        <span class="n">invertor</span><span class="o">.</span><span class="n">has_bck</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span>
291        <span class="n">invertor</span><span class="o">.</span><span class="n">slit_height</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">slit_height</span>
292        <span class="n">invertor</span><span class="o">.</span><span class="n">slit_width</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">slit_width</span>
293       
294        <span class="n">invertor</span><span class="o">.</span><span class="n">info</span> <span class="o">=</span> <span class="n">copy</span><span class="o">.</span><span class="n">deepcopy</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">info</span><span class="p">)</span>
295       
296        <span class="k">return</span> <span class="n">invertor</span>
297    </div>
298<div class="viewcode-block" id="Invertor.invert"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.invert">[docs]</a>    <span class="k">def</span> <span class="nf">invert</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">nfunc</span><span class="o">=</span><span class="mi">10</span><span class="p">,</span> <span class="n">nr</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
299        <span class="sd">&quot;&quot;&quot;</span>
300<span class="sd">        Perform inversion to P(r)</span>
301<span class="sd">        </span>
302<span class="sd">        The problem is solved by posing the problem as  Ax = b,</span>
303<span class="sd">        where x is the set of coefficients we are looking for.</span>
304<span class="sd">        </span>
305<span class="sd">        Npts is the number of points.</span>
306<span class="sd">        </span>
307<span class="sd">        In the following i refers to the ith base function coefficient.</span>
308<span class="sd">        The matrix has its entries j in its first Npts rows set to ::</span>
309
310<span class="sd">            A[i][j] = (Fourier transformed base function for point j)</span>
311<span class="sd">            </span>
312<span class="sd">        We them choose a number of r-points, n_r, to evaluate the second</span>
313<span class="sd">        derivative of P(r) at. This is used as our regularization term.</span>
314<span class="sd">        For a vector r of length n_r, the following n_r rows are set to ::</span>
315
316<span class="sd">            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j])</span>
317<span class="sd">            </span>
318<span class="sd">        The vector b has its first Npts entries set to ::</span>
319
320<span class="sd">            b[j] = (I(q) observed for point j)</span>
321<span class="sd">            </span>
322<span class="sd">        The following n_r entries are set to zero.</span>
323<span class="sd">        </span>
324<span class="sd">        The result is found by using scipy.linalg.basic.lstsq to invert</span>
325<span class="sd">        the matrix and find the coefficients x.</span>
326<span class="sd">        </span>
327<span class="sd">        :param nfunc: number of base functions to use.</span>
328<span class="sd">        :param nr: number of r points to evaluate the 2nd derivative at for the reg. term.</span>
329<span class="sd">        :return: c_out, c_cov - the coefficients with covariance matrix</span>
330<span class="sd">        &quot;&quot;&quot;</span>
331        <span class="c"># Reset the background value before proceeding</span>
332        <span class="bp">self</span><span class="o">.</span><span class="n">background</span> <span class="o">=</span> <span class="mf">0.0</span>
333        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">lstsq</span><span class="p">(</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nr</span><span class="o">=</span><span class="n">nr</span><span class="p">)</span>
334    </div>
335<div class="viewcode-block" id="Invertor.iq"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.iq">[docs]</a>    <span class="k">def</span> <span class="nf">iq</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">out</span><span class="p">,</span> <span class="n">q</span><span class="p">):</span>
336        <span class="sd">&quot;&quot;&quot;</span>
337<span class="sd">        Function to call to evaluate the scattering intensity</span>
338<span class="sd">        </span>
339<span class="sd">        :param args: c-parameters, and q</span>
340<span class="sd">        :return: I(q)</span>
341<span class="sd">        </span>
342<span class="sd">        &quot;&quot;&quot;</span>
343        <span class="k">return</span> <span class="n">Cinvertor</span><span class="o">.</span><span class="n">iq</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">out</span><span class="p">,</span> <span class="n">q</span><span class="p">)</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">background</span>
344    </div>
345<div class="viewcode-block" id="Invertor.invert_optimize"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.invert_optimize">[docs]</a>    <span class="k">def</span> <span class="nf">invert_optimize</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">nfunc</span><span class="o">=</span><span class="mi">10</span><span class="p">,</span> <span class="n">nr</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
346        <span class="sd">&quot;&quot;&quot;</span>
347<span class="sd">        Slower version of the P(r) inversion that uses scipy.optimize.leastsq.</span>
348<span class="sd">        </span>
349<span class="sd">        This probably produce more reliable results, but is much slower.</span>
350<span class="sd">        The minimization function is set to</span>
351<span class="sd">        sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term,</span>
352<span class="sd">        where the reg_term is given by Svergun: it is the integral of</span>
353<span class="sd">        the square of the first derivative</span>
354<span class="sd">        of P(r), d(P(r))/dr, integrated over the full range of r.</span>
355<span class="sd">        </span>
356<span class="sd">        :param nfunc: number of base functions to use.</span>
357<span class="sd">        :param nr: number of r points to evaluate the 2nd derivative at</span>
358<span class="sd">            for the reg. term.</span>
359<span class="sd">        </span>
360<span class="sd">        :return: c_out, c_cov - the coefficients with covariance matrix</span>
361<span class="sd">        </span>
362<span class="sd">        &quot;&quot;&quot;</span>
363        <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span> <span class="o">=</span> <span class="n">nfunc</span>
364        <span class="c"># First, check that the current data is valid</span>
365        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">is_valid</span><span class="p">()</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">:</span>
366            <span class="n">msg</span> <span class="o">=</span> <span class="s">&quot;Invertor.invert: Data array are of different length&quot;</span>
367            <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">msg</span>
368       
369        <span class="n">p</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="n">nfunc</span><span class="p">)</span>
370        <span class="n">t_0</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
371        <span class="n">out</span><span class="p">,</span> <span class="n">cov_x</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">optimize</span><span class="o">.</span><span class="n">leastsq</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">residuals</span><span class="p">,</span>
372                                                            <span class="n">p</span><span class="p">,</span> <span class="n">full_output</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
373       
374        <span class="c"># Compute chi^2</span>
375        <span class="n">res</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">residuals</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
376        <span class="n">chisqr</span> <span class="o">=</span> <span class="mi">0</span>
377        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">res</span><span class="p">)):</span>
378            <span class="n">chisqr</span> <span class="o">+=</span> <span class="n">res</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
379       
380        <span class="bp">self</span><span class="o">.</span><span class="n">chi2</span> <span class="o">=</span> <span class="n">chisqr</span>
381
382        <span class="c"># Store computation time</span>
383        <span class="bp">self</span><span class="o">.</span><span class="n">elapsed</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span> <span class="o">-</span> <span class="n">t_0</span>
384       
385        <span class="k">if</span> <span class="n">cov_x</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
386            <span class="n">cov_x</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ones</span><span class="p">([</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nfunc</span><span class="p">])</span>
387            <span class="n">cov_x</span> <span class="o">*=</span> <span class="n">math</span><span class="o">.</span><span class="n">fabs</span><span class="p">(</span><span class="n">chisqr</span><span class="p">)</span>
388        <span class="k">return</span> <span class="n">out</span><span class="p">,</span> <span class="n">cov_x</span>
389    </div>
390<div class="viewcode-block" id="Invertor.pr_fit"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.pr_fit">[docs]</a>    <span class="k">def</span> <span class="nf">pr_fit</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">nfunc</span><span class="o">=</span><span class="mi">5</span><span class="p">):</span>
391        <span class="sd">&quot;&quot;&quot;</span>
392<span class="sd">        This is a direct fit to a given P(r). It assumes that the y data</span>
393<span class="sd">        is set to some P(r) distribution that we are trying to reproduce</span>
394<span class="sd">        with a set of base functions.</span>
395<span class="sd">        </span>
396<span class="sd">        This method is provided as a test.</span>
397<span class="sd">        &quot;&quot;&quot;</span>
398        <span class="c"># First, check that the current data is valid</span>
399        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">is_valid</span><span class="p">()</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">:</span>
400            <span class="n">msg</span> <span class="o">=</span> <span class="s">&quot;Invertor.invert: Data arrays are of different length&quot;</span>
401            <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">msg</span>
402       
403        <span class="n">p</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="n">nfunc</span><span class="p">)</span>
404        <span class="n">t_0</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
405        <span class="n">out</span><span class="p">,</span> <span class="n">cov_x</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">optimize</span><span class="o">.</span><span class="n">leastsq</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">pr_residuals</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span>
406                                                            <span class="n">full_output</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
407       
408        <span class="c"># Compute chi^2</span>
409        <span class="n">res</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">pr_residuals</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
410        <span class="n">chisqr</span> <span class="o">=</span> <span class="mi">0</span>
411        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">res</span><span class="p">)):</span>
412            <span class="n">chisqr</span> <span class="o">+=</span> <span class="n">res</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
413       
414        <span class="bp">self</span><span class="o">.</span><span class="n">chisqr</span> <span class="o">=</span> <span class="n">chisqr</span>
415       
416        <span class="c"># Store computation time</span>
417        <span class="bp">self</span><span class="o">.</span><span class="n">elapsed</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span> <span class="o">-</span> <span class="n">t_0</span>
418
419        <span class="k">return</span> <span class="n">out</span><span class="p">,</span> <span class="n">cov_x</span>
420    </div>
421<div class="viewcode-block" id="Invertor.pr_err"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.pr_err">[docs]</a>    <span class="k">def</span> <span class="nf">pr_err</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">c</span><span class="p">,</span> <span class="n">c_cov</span><span class="p">,</span> <span class="n">r</span><span class="p">):</span>
422        <span class="sd">&quot;&quot;&quot;</span>
423<span class="sd">        Returns the value of P(r) for a given r, and base function</span>
424<span class="sd">        coefficients, with error.</span>
425<span class="sd">        </span>
426<span class="sd">        :param c: base function coefficients</span>
427<span class="sd">        :param c_cov: covariance matrice of the base function coefficients</span>
428<span class="sd">        :param r: r-value to evaluate P(r) at</span>
429<span class="sd">        </span>
430<span class="sd">        :return: P(r)</span>
431<span class="sd">        </span>
432<span class="sd">        &quot;&quot;&quot;</span>
433        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_pr_err</span><span class="p">(</span><span class="n">c</span><span class="p">,</span> <span class="n">c_cov</span><span class="p">,</span> <span class="n">r</span><span class="p">)</span>
434       </div>
435    <span class="k">def</span> <span class="nf">_accept_q</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">q</span><span class="p">):</span>
436        <span class="sd">&quot;&quot;&quot;</span>
437<span class="sd">        Check q-value against user-defined range</span>
438<span class="sd">        &quot;&quot;&quot;</span>
439        <span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span> <span class="o">==</span> <span class="bp">None</span> <span class="ow">and</span> <span class="n">q</span> <span class="o">&lt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span><span class="p">:</span>
440            <span class="k">return</span> <span class="bp">False</span>
441        <span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span> <span class="o">==</span> <span class="bp">None</span> <span class="ow">and</span> <span class="n">q</span> <span class="o">&gt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span><span class="p">:</span>
442            <span class="k">return</span> <span class="bp">False</span>
443        <span class="k">return</span> <span class="bp">True</span>
444       
445<div class="viewcode-block" id="Invertor.lstsq"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.lstsq">[docs]</a>    <span class="k">def</span> <span class="nf">lstsq</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">nfunc</span><span class="o">=</span><span class="mi">5</span><span class="p">,</span> <span class="n">nr</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
446        <span class="sd">&quot;&quot;&quot;</span>
447<span class="sd">        The problem is solved by posing the problem as  Ax = b,</span>
448<span class="sd">        where x is the set of coefficients we are looking for.</span>
449<span class="sd">        </span>
450<span class="sd">        Npts is the number of points.</span>
451<span class="sd">        </span>
452<span class="sd">        In the following i refers to the ith base function coefficient.</span>
453<span class="sd">        The matrix has its entries j in its first Npts rows set to ::</span>
454
455<span class="sd">            A[i][j] = (Fourier transformed base function for point j)</span>
456<span class="sd">            </span>
457<span class="sd">        We them choose a number of r-points, n_r, to evaluate the second</span>
458<span class="sd">        derivative of P(r) at. This is used as our regularization term.</span>
459<span class="sd">        For a vector r of length n_r, the following n_r rows are set to ::</span>
460
461<span class="sd">            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,</span>
462<span class="sd">            evaluated at r[j])</span>
463<span class="sd">            </span>
464<span class="sd">        The vector b has its first Npts entries set to ::</span>
465
466<span class="sd">            b[j] = (I(q) observed for point j)</span>
467<span class="sd">            </span>
468<span class="sd">        The following n_r entries are set to zero.</span>
469<span class="sd">        </span>
470<span class="sd">        The result is found by using scipy.linalg.basic.lstsq to invert</span>
471<span class="sd">        the matrix and find the coefficients x.</span>
472<span class="sd">        </span>
473<span class="sd">        :param nfunc: number of base functions to use.</span>
474<span class="sd">        :param nr: number of r points to evaluate the 2nd derivative at for the reg. term.</span>
475
476<span class="sd">        If the result does not allow us to compute the covariance matrix,</span>
477<span class="sd">        a matrix filled with zeros will be returned.</span>
478
479<span class="sd">        &quot;&quot;&quot;</span>
480        <span class="c"># Note: To make sure an array is contiguous:</span>
481        <span class="c"># blah = numpy.ascontiguousarray(blah_original)</span>
482        <span class="c"># ... before passing it to C</span>
483       
484        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">is_valid</span><span class="p">()</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span>
485            <span class="n">msg</span> <span class="o">=</span> <span class="s">&quot;Invertor: invalid data; incompatible data lengths.&quot;</span>
486            <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">msg</span>
487       
488        <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span> <span class="o">=</span> <span class="n">nfunc</span>
489        <span class="c"># a -- An M x N matrix.</span>
490        <span class="c"># b -- An M x nrhs matrix or M vector.</span>
491        <span class="n">npts</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">x</span><span class="p">)</span>
492        <span class="n">nq</span>   <span class="o">=</span> <span class="n">nr</span>
493        <span class="n">sqrt_alpha</span> <span class="o">=</span> <span class="n">math</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">math</span><span class="o">.</span><span class="n">fabs</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">alpha</span><span class="p">))</span>
494        <span class="k">if</span> <span class="n">sqrt_alpha</span> <span class="o">&lt;</span> <span class="mf">0.0</span><span class="p">:</span>
495            <span class="n">nq</span> <span class="o">=</span> <span class="mi">0</span>
496
497        <span class="c"># If we need to fit the background, add a term</span>
498        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span> <span class="o">==</span> <span class="bp">True</span><span class="p">:</span>
499            <span class="n">nfunc_0</span> <span class="o">=</span> <span class="n">nfunc</span>
500            <span class="n">nfunc</span> <span class="o">+=</span> <span class="mi">1</span>
501
502        <span class="n">a</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="n">npts</span> <span class="o">+</span> <span class="n">nq</span><span class="p">,</span> <span class="n">nfunc</span><span class="p">])</span>
503        <span class="n">b</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">npts</span> <span class="o">+</span> <span class="n">nq</span><span class="p">)</span>
504        <span class="n">err</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nfunc</span><span class="p">])</span>
505       
506        <span class="c"># Construct the a matrix and b vector that represent the problem</span>
507        <span class="n">t_0</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
508        <span class="k">try</span><span class="p">:</span>
509            <span class="bp">self</span><span class="o">.</span><span class="n">_get_matrix</span><span class="p">(</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nq</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">)</span>
510        <span class="k">except</span><span class="p">:</span>
511            <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="s">&quot;Invertor: could not invert I(Q)</span><span class="se">\n</span><span class="s">  </span><span class="si">%s</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">sys</span><span class="o">.</span><span class="n">exc_value</span>
512             
513        <span class="c"># Perform the inversion (least square fit)</span>
514        <span class="n">c</span><span class="p">,</span> <span class="n">chi2</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">lstsq</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">)</span>
515        <span class="c"># Sanity check</span>
516        <span class="k">try</span><span class="p">:</span>
517            <span class="nb">float</span><span class="p">(</span><span class="n">chi2</span><span class="p">)</span>
518        <span class="k">except</span><span class="p">:</span>
519            <span class="n">chi2</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.0</span>
520        <span class="bp">self</span><span class="o">.</span><span class="n">chi2</span> <span class="o">=</span> <span class="n">chi2</span>
521               
522        <span class="n">inv_cov</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nfunc</span><span class="p">])</span>
523        <span class="c"># Get the covariance matrix, defined as inv_cov = a_transposed * a</span>
524        <span class="bp">self</span><span class="o">.</span><span class="n">_get_invcov_matrix</span><span class="p">(</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nr</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">inv_cov</span><span class="p">)</span>
525                   
526        <span class="c"># Compute the reg term size for the output</span>
527        <span class="n">sum_sig</span><span class="p">,</span> <span class="n">sum_reg</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">_get_reg_size</span><span class="p">(</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nr</span><span class="p">,</span> <span class="n">a</span><span class="p">)</span>
528                   
529        <span class="k">if</span> <span class="n">math</span><span class="o">.</span><span class="n">fabs</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">alpha</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
530            <span class="n">new_alpha</span> <span class="o">=</span> <span class="n">sum_sig</span> <span class="o">/</span> <span class="p">(</span><span class="n">sum_reg</span> <span class="o">/</span> <span class="bp">self</span><span class="o">.</span><span class="n">alpha</span><span class="p">)</span>
531        <span class="k">else</span><span class="p">:</span>
532            <span class="n">new_alpha</span> <span class="o">=</span> <span class="mf">0.0</span>
533        <span class="bp">self</span><span class="o">.</span><span class="n">suggested_alpha</span> <span class="o">=</span> <span class="n">new_alpha</span>
534       
535        <span class="k">try</span><span class="p">:</span>
536            <span class="n">cov</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">inv_cov</span><span class="p">)</span>
537            <span class="n">err</span> <span class="o">=</span> <span class="n">math</span><span class="o">.</span><span class="n">fabs</span><span class="p">(</span><span class="n">chi2</span> <span class="o">/</span> <span class="nb">float</span><span class="p">(</span><span class="n">npts</span> <span class="o">-</span> <span class="n">nfunc</span><span class="p">))</span> <span class="o">*</span> <span class="n">cov</span>
538        <span class="k">except</span><span class="p">:</span>
539            <span class="c"># We were not able to estimate the errors</span>
540            <span class="c"># Return an empty error matrix</span>
541            <span class="k">pass</span>
542           
543        <span class="c"># Keep a copy of the last output</span>
544        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span> <span class="o">==</span> <span class="bp">False</span><span class="p">:</span>
545            <span class="bp">self</span><span class="o">.</span><span class="n">background</span> <span class="o">=</span> <span class="mi">0</span>
546            <span class="bp">self</span><span class="o">.</span><span class="n">out</span> <span class="o">=</span> <span class="n">c</span>
547            <span class="bp">self</span><span class="o">.</span><span class="n">cov</span> <span class="o">=</span> <span class="n">err</span>
548        <span class="k">else</span><span class="p">:</span>
549            <span class="bp">self</span><span class="o">.</span><span class="n">background</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
550           
551            <span class="n">err_0</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">nfunc</span><span class="p">])</span>
552            <span class="n">c_0</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nfunc</span><span class="p">)</span>
553           
554            <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nfunc_0</span><span class="p">):</span>
555                <span class="n">c_0</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span>
556                <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nfunc_0</span><span class="p">):</span>
557                    <span class="n">err_0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">err</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span>
558                   
559            <span class="bp">self</span><span class="o">.</span><span class="n">out</span> <span class="o">=</span> <span class="n">c_0</span>
560            <span class="bp">self</span><span class="o">.</span><span class="n">cov</span> <span class="o">=</span> <span class="n">err_0</span>
561           
562        <span class="c"># Store computation time</span>
563        <span class="bp">self</span><span class="o">.</span><span class="n">elapsed</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span> <span class="o">-</span> <span class="n">t_0</span>
564       
565        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">out</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">cov</span>
566        </div>
567<div class="viewcode-block" id="Invertor.estimate_numterms"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.estimate_numterms">[docs]</a>    <span class="k">def</span> <span class="nf">estimate_numterms</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">isquit_func</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
568        <span class="sd">&quot;&quot;&quot;</span>
569<span class="sd">        Returns a reasonable guess for the</span>
570<span class="sd">        number of terms</span>
571<span class="sd">        </span>
572<span class="sd">        :param isquit_func:</span>
573<span class="sd">          reference to thread function to call to check whether the computation needs to</span>
574<span class="sd">          be stopped.</span>
575<span class="sd">        </span>
576<span class="sd">        :return: number of terms, alpha, message</span>
577<span class="sd">        </span>
578<span class="sd">        &quot;&quot;&quot;</span>
579        <span class="kn">from</span> <span class="nn">num_term</span> <span class="kn">import</span> <span class="n">Num_terms</span>
580        <span class="n">estimator</span> <span class="o">=</span> <span class="n">Num_terms</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">clone</span><span class="p">())</span>
581        <span class="k">try</span><span class="p">:</span>
582            <span class="k">return</span> <span class="n">estimator</span><span class="o">.</span><span class="n">num_terms</span><span class="p">(</span><span class="n">isquit_func</span><span class="p">)</span>
583        <span class="k">except</span><span class="p">:</span>
584            <span class="c"># If we fail, estimate alpha and return the default</span>
585            <span class="c"># number of terms</span>
586            <span class="n">best_alpha</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">estimate_alpha</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span><span class="p">)</span>
587            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span><span class="p">,</span> <span class="n">best_alpha</span><span class="p">,</span> <span class="s">&quot;Could not estimate number of terms&quot;</span>
588                    </div>
589<div class="viewcode-block" id="Invertor.estimate_alpha"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.estimate_alpha">[docs]</a>    <span class="k">def</span> <span class="nf">estimate_alpha</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">nfunc</span><span class="p">):</span>
590        <span class="sd">&quot;&quot;&quot;</span>
591<span class="sd">        Returns a reasonable guess for the</span>
592<span class="sd">        regularization constant alpha</span>
593<span class="sd">        </span>
594<span class="sd">        :param nfunc: number of terms to use in the expansion.</span>
595<span class="sd">        </span>
596<span class="sd">        :return: alpha, message, elapsed</span>
597<span class="sd">        </span>
598<span class="sd">        where alpha is the estimate for alpha,</span>
599<span class="sd">        message is a message for the user,</span>
600<span class="sd">        elapsed is the computation time</span>
601<span class="sd">        &quot;&quot;&quot;</span>
602        <span class="c">#import time</span>
603        <span class="k">try</span><span class="p">:</span>
604            <span class="n">pr</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">clone</span><span class="p">()</span>
605           
606            <span class="c"># T_0 for computation time</span>
607            <span class="n">starttime</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
608            <span class="n">elapsed</span> <span class="o">=</span> <span class="mi">0</span>
609           
610            <span class="c"># If the current alpha is zero, try</span>
611            <span class="c"># another value</span>
612            <span class="k">if</span> <span class="n">pr</span><span class="o">.</span><span class="n">alpha</span> <span class="o">&lt;=</span> <span class="mi">0</span><span class="p">:</span>
613                <span class="n">pr</span><span class="o">.</span><span class="n">alpha</span> <span class="o">=</span> <span class="mf">0.0001</span>
614                 
615            <span class="c"># Perform inversion to find the largest alpha</span>
616            <span class="n">out</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">invert</span><span class="p">(</span><span class="n">nfunc</span><span class="p">)</span>
617            <span class="n">elapsed</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span> <span class="o">-</span> <span class="n">starttime</span>
618            <span class="n">initial_alpha</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">alpha</span>
619            <span class="n">initial_peaks</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">get_peaks</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
620   
621            <span class="c"># Try the inversion with the estimated alpha</span>
622            <span class="n">pr</span><span class="o">.</span><span class="n">alpha</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">suggested_alpha</span>
623            <span class="n">out</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">invert</span><span class="p">(</span><span class="n">nfunc</span><span class="p">)</span>
624   
625            <span class="n">npeaks</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">get_peaks</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
626            <span class="c"># if more than one peak to start with</span>
627            <span class="c"># just return the estimate</span>
628            <span class="k">if</span> <span class="n">npeaks</span> <span class="o">&gt;</span> <span class="mi">1</span><span class="p">:</span>
629                <span class="c">#message = &quot;Your P(r) is not smooth,</span>
630                <span class="c">#please check your inversion parameters&quot;</span>
631                <span class="n">message</span> <span class="o">=</span> <span class="bp">None</span>
632                <span class="k">return</span> <span class="n">pr</span><span class="o">.</span><span class="n">suggested_alpha</span><span class="p">,</span> <span class="n">message</span><span class="p">,</span> <span class="n">elapsed</span>
633            <span class="k">else</span><span class="p">:</span>
634               
635                <span class="c"># Look at smaller values</span>
636                <span class="c"># We assume that for the suggested alpha, we have 1 peak</span>
637                <span class="c"># if not, send a message to change parameters</span>
638                <span class="n">alpha</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">suggested_alpha</span>
639                <span class="n">best_alpha</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">suggested_alpha</span>
640                <span class="n">found</span> <span class="o">=</span> <span class="bp">False</span>
641                <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">10</span><span class="p">):</span>
642                    <span class="n">pr</span><span class="o">.</span><span class="n">alpha</span> <span class="o">=</span> <span class="p">(</span><span class="mf">0.33</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">*</span> <span class="n">alpha</span>
643                    <span class="n">out</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">invert</span><span class="p">(</span><span class="n">nfunc</span><span class="p">)</span>
644                   
645                    <span class="n">peaks</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">get_peaks</span><span class="p">(</span><span class="n">out</span><span class="p">)</span>
646                    <span class="k">if</span> <span class="n">peaks</span> <span class="o">&gt;</span> <span class="mi">1</span><span class="p">:</span>
647                        <span class="n">found</span> <span class="o">=</span> <span class="bp">True</span>
648                        <span class="k">break</span>
649                    <span class="n">best_alpha</span> <span class="o">=</span> <span class="n">pr</span><span class="o">.</span><span class="n">alpha</span>
650                   
651                <span class="c"># If we didn&#39;t find a turning point for alpha and</span>
652                <span class="c"># the initial alpha already had only one peak,</span>
653                <span class="c"># just return that</span>
654                <span class="k">if</span> <span class="ow">not</span> <span class="n">found</span> <span class="ow">and</span> <span class="n">initial_peaks</span> <span class="o">==</span> <span class="mi">1</span> <span class="ow">and</span> \
655                    <span class="n">initial_alpha</span> <span class="o">&lt;</span> <span class="n">best_alpha</span><span class="p">:</span>
656                    <span class="n">best_alpha</span> <span class="o">=</span> <span class="n">initial_alpha</span>
657                   
658                <span class="c"># Check whether the size makes sense</span>
659                <span class="n">message</span> <span class="o">=</span> <span class="s">&#39;&#39;</span>
660               
661                <span class="k">if</span> <span class="ow">not</span> <span class="n">found</span><span class="p">:</span>
662                    <span class="n">message</span> <span class="o">=</span> <span class="bp">None</span>
663                <span class="k">elif</span> <span class="n">best_alpha</span> <span class="o">&gt;=</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="n">pr</span><span class="o">.</span><span class="n">suggested_alpha</span><span class="p">:</span>
664                    <span class="c"># best alpha is too big, return a</span>
665                    <span class="c"># reasonable value</span>
666                    <span class="n">message</span>  <span class="o">=</span> <span class="s">&quot;The estimated alpha for your system is too &quot;</span>
667                    <span class="n">message</span> <span class="o">+=</span> <span class="s">&quot;large. &quot;</span>
668                    <span class="n">message</span> <span class="o">+=</span> <span class="s">&quot;Try increasing your maximum distance.&quot;</span>
669               
670                <span class="k">return</span> <span class="n">best_alpha</span><span class="p">,</span> <span class="n">message</span><span class="p">,</span> <span class="n">elapsed</span>
671   
672        <span class="k">except</span><span class="p">:</span>
673            <span class="n">message</span> <span class="o">=</span> <span class="s">&quot;Invertor.estimate_alpha: </span><span class="si">%s</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">sys</span><span class="o">.</span><span class="n">exc_value</span>
674            <span class="k">return</span> <span class="mi">0</span><span class="p">,</span> <span class="n">message</span><span class="p">,</span> <span class="n">elapsed</span>
675    </div>
676<div class="viewcode-block" id="Invertor.to_file"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.to_file">[docs]</a>    <span class="k">def</span> <span class="nf">to_file</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">path</span><span class="p">,</span> <span class="n">npts</span><span class="o">=</span><span class="mi">100</span><span class="p">):</span>
677        <span class="sd">&quot;&quot;&quot;</span>
678<span class="sd">        Save the state to a file that will be readable</span>
679<span class="sd">        by SliceView.</span>
680<span class="sd">        </span>
681<span class="sd">        :param path: path of the file to write</span>
682<span class="sd">        :param npts: number of P(r) points to be written</span>
683<span class="sd">        </span>
684<span class="sd">        &quot;&quot;&quot;</span>
685        <span class="nb">file</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">path</span><span class="p">,</span> <span class="s">&#39;w&#39;</span><span class="p">)</span>
686        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#d_max=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span><span class="p">)</span>
687        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#nfunc=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span><span class="p">)</span>
688        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#alpha=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">alpha</span><span class="p">)</span>
689        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#chi2=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">chi2</span><span class="p">)</span>
690        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#elapsed=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">elapsed</span><span class="p">)</span>
691        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#qmin=</span><span class="si">%s</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">q_min</span><span class="p">))</span>
692        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#qmax=</span><span class="si">%s</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">q_max</span><span class="p">))</span>
693        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#slit_height=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">slit_height</span><span class="p">)</span>
694        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#slit_width=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">slit_width</span><span class="p">)</span>
695        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#background=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">background</span><span class="p">)</span>
696        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span> <span class="o">==</span> <span class="bp">True</span><span class="p">:</span>
697            <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#has_bck=1</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">)</span>
698        <span class="k">else</span><span class="p">:</span>
699            <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#has_bck=0</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">)</span>
700        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#alpha_estimate=</span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">suggested_alpha</span><span class="p">)</span>
701        <span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">out</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
702            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">out</span><span class="p">)</span> <span class="o">==</span> <span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">cov</span><span class="p">):</span>
703                <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">out</span><span class="p">)):</span>
704                    <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;#C_</span><span class="si">%i</span><span class="s">=</span><span class="si">%s</span><span class="s">+-</span><span class="si">%s</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="nb">str</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">out</span><span class="p">[</span><span class="n">i</span><span class="p">]),</span>
705                                                    <span class="nb">str</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">cov</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">])))</span>
706        <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;&lt;r&gt;  &lt;Pr&gt;  &lt;dPr&gt;</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">)</span>
707        <span class="n">r</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span><span class="o">/</span><span class="n">npts</span><span class="p">)</span>
708       
709        <span class="k">for</span> <span class="n">r_i</span> <span class="ow">in</span> <span class="n">r</span><span class="p">:</span>
710            <span class="p">(</span><span class="n">value</span><span class="p">,</span> <span class="n">err</span><span class="p">)</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">pr_err</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">out</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">cov</span><span class="p">,</span> <span class="n">r_i</span><span class="p">)</span>
711            <span class="nb">file</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;</span><span class="si">%g</span><span class="s">  </span><span class="si">%g</span><span class="s">  </span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">r_i</span><span class="p">,</span> <span class="n">value</span><span class="p">,</span> <span class="n">err</span><span class="p">))</span>
712   
713        <span class="nb">file</span><span class="o">.</span><span class="n">close</span><span class="p">()</span>
714     </div>
715<div class="viewcode-block" id="Invertor.from_file"><a class="viewcode-back" href="../../../dev/api/sas.pr.html#sas.pr.invertor.Invertor.from_file">[docs]</a>    <span class="k">def</span> <span class="nf">from_file</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">path</span><span class="p">):</span>
716        <span class="sd">&quot;&quot;&quot;</span>
717<span class="sd">        Load the state of the Invertor from a file,</span>
718<span class="sd">        to be able to generate P(r) from a set of</span>
719<span class="sd">        parameters.</span>
720<span class="sd">        </span>
721<span class="sd">        :param path: path of the file to load</span>
722<span class="sd">        </span>
723<span class="sd">        &quot;&quot;&quot;</span>
724        <span class="c">#import os</span>
725        <span class="c">#import re</span>
726        <span class="k">if</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">isfile</span><span class="p">(</span><span class="n">path</span><span class="p">):</span>
727            <span class="k">try</span><span class="p">:</span>
728                <span class="n">fd</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">path</span><span class="p">,</span> <span class="s">&#39;r&#39;</span><span class="p">)</span>
729               
730                <span class="n">buff</span> <span class="o">=</span> <span class="n">fd</span><span class="o">.</span><span class="n">read</span><span class="p">()</span>
731                <span class="n">lines</span> <span class="o">=</span> <span class="n">buff</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;</span><span class="se">\n</span><span class="s">&#39;</span><span class="p">)</span>
732                <span class="k">for</span> <span class="n">line</span> <span class="ow">in</span> <span class="n">lines</span><span class="p">:</span>
733                    <span class="k">if</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#d_max=&#39;</span><span class="p">):</span>
734                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
735                        <span class="bp">self</span><span class="o">.</span><span class="n">d_max</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
736                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#nfunc=&#39;</span><span class="p">):</span>
737                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
738                        <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
739                        <span class="bp">self</span><span class="o">.</span><span class="n">out</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span><span class="p">)</span>
740                        <span class="bp">self</span><span class="o">.</span><span class="n">cov</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nfunc</span><span class="p">])</span>
741                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#alpha=&#39;</span><span class="p">):</span>
742                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
743                        <span class="bp">self</span><span class="o">.</span><span class="n">alpha</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
744                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#chi2=&#39;</span><span class="p">):</span>
745                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
746                        <span class="bp">self</span><span class="o">.</span><span class="n">chi2</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
747                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#elapsed=&#39;</span><span class="p">):</span>
748                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
749                        <span class="bp">self</span><span class="o">.</span><span class="n">elapsed</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
750                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#alpha_estimate=&#39;</span><span class="p">):</span>
751                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
752                        <span class="bp">self</span><span class="o">.</span><span class="n">suggested_alpha</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
753                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#qmin=&#39;</span><span class="p">):</span>
754                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
755                        <span class="k">try</span><span class="p">:</span>
756                            <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
757                        <span class="k">except</span><span class="p">:</span>
758                            <span class="bp">self</span><span class="o">.</span><span class="n">q_min</span> <span class="o">=</span> <span class="bp">None</span>
759                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#qmax=&#39;</span><span class="p">):</span>
760                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
761                        <span class="k">try</span><span class="p">:</span>
762                            <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
763                        <span class="k">except</span><span class="p">:</span>
764                            <span class="bp">self</span><span class="o">.</span><span class="n">q_max</span> <span class="o">=</span> <span class="bp">None</span>
765                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#slit_height=&#39;</span><span class="p">):</span>
766                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
767                        <span class="bp">self</span><span class="o">.</span><span class="n">slit_height</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
768                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#slit_width=&#39;</span><span class="p">):</span>
769                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
770                        <span class="bp">self</span><span class="o">.</span><span class="n">slit_width</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
771                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#background=&#39;</span><span class="p">):</span>
772                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
773                        <span class="bp">self</span><span class="o">.</span><span class="n">background</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
774                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#has_bck=&#39;</span><span class="p">):</span>
775                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
776                        <span class="k">if</span> <span class="nb">int</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
777                            <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span> <span class="o">=</span> <span class="bp">True</span>
778                        <span class="k">else</span><span class="p">:</span>
779                            <span class="bp">self</span><span class="o">.</span><span class="n">has_bck</span> <span class="o">=</span> <span class="bp">False</span>
780           
781                    <span class="c"># Now read in the parameters</span>
782                    <span class="k">elif</span> <span class="n">line</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#C_&#39;</span><span class="p">):</span>
783                        <span class="n">toks</span> <span class="o">=</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;=&#39;</span><span class="p">)</span>
784                        <span class="n">p</span> <span class="o">=</span> <span class="n">re</span><span class="o">.</span><span class="n">compile</span><span class="p">(</span><span class="s">&#39;#C_([0-9]+)&#39;</span><span class="p">)</span>
785                        <span class="n">m</span> <span class="o">=</span> <span class="n">p</span><span class="o">.</span><span class="n">search</span><span class="p">(</span><span class="n">toks</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
786                        <span class="n">toks2</span> <span class="o">=</span> <span class="n">toks</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&#39;+-&#39;</span><span class="p">)</span>
787                        <span class="n">i</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">m</span><span class="o">.</span><span class="n">group</span><span class="p">(</span><span class="mi">1</span><span class="p">))</span>
788                        <span class="bp">self</span><span class="o">.</span><span class="n">out</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks2</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
789                       
790                        <span class="bp">self</span><span class="o">.</span><span class="n">cov</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">toks2</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
791           
792            <span class="k">except</span><span class="p">:</span>
793                <span class="n">msg</span> <span class="o">=</span> <span class="s">&quot;Invertor.from_file: corrupted file</span><span class="se">\n</span><span class="si">%s</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">sys</span><span class="o">.</span><span class="n">exc_value</span>
794                <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">msg</span>
795        <span class="k">else</span><span class="p">:</span>
796            <span class="n">msg</span> <span class="o">=</span> <span class="s">&quot;Invertor.from_file: &#39;</span><span class="si">%s</span><span class="s">&#39; is not a file&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="n">path</span><span class="p">)</span>
797            <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">msg</span></div></div>
798</pre></div>
799
800          </div>
801        </div>
802      </div>
803      <div class="sphinxsidebar">
804        <div class="sphinxsidebarwrapper">
805<div id="searchbox" style="display: none">
806  <h3>Quick search</h3>
807    <form class="search" action="../../../search.html" method="get">
808      <input type="text" name="q" />
809      <input type="submit" value="Go" />
810      <input type="hidden" name="check_keywords" value="yes" />
811      <input type="hidden" name="area" value="default" />
812    </form>
813    <p class="searchtip" style="font-size: 90%">
814    Enter search terms or a module, class or function name.
815    </p>
816</div>
817<script type="text/javascript">$('#searchbox').show(0);</script>
818        </div>
819      </div>
820      <div class="clearer"></div>
821    </div>
822    <div class="related">
823      <h3>Navigation</h3>
824      <ul>
825        <li class="right" style="margin-right: 10px">
826          <a href="../../../genindex.html" title="General Index"
827             >index</a></li>
828        <li class="right" >
829          <a href="../../../py-modindex.html" title="Python Module Index"
830             >modules</a> |</li>
831        <li><a href="../../../index.html">SasView 3.0.0 documentation</a> &raquo;</li>
832          <li><a href="../../index.html" >Module code</a> &raquo;</li> 
833      </ul>
834    </div>
835    <div class="footer">
836        &copy; Copyright 2013, The SasView Project.
837      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.2.3.
838    </div>
839  </body>
840</html>
Note: See TracBrowser for help on using the repository browser.