Index: pr_inversion/c_extensions/Cinvertor.c
===================================================================
--- pr_inversion/c_extensions/Cinvertor.c	(revision 4f63160d4aa8e41b4460f847aceb7324d62f7175)
+++ pr_inversion/c_extensions/Cinvertor.c	(revision f71287f458148547d69aae7ddc900f1afde663f6)
@@ -266,4 +266,42 @@
 	return Py_BuildValue("d", self->params.d_max);	
 }
+
+/**
+ * Sets the minimum q
+ */
+static PyObject * set_qmin(Cinvertor *self, PyObject *args) {
+	double q_min;
+  
+	if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL;
+	self->params.q_min = q_min;
+	return Py_BuildValue("d", self->params.q_min);	
+}
+
+/**
+ * Gets the minimum q
+ */
+static PyObject * get_qmin(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.q_min);	
+}
+
+
+/**
+ * Sets the maximum q
+ */
+static PyObject * set_qmax(Cinvertor *self, PyObject *args) {
+	double q_max;
+  
+	if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL;
+	self->params.q_max = q_max;
+	return Py_BuildValue("d", self->params.q_max);	
+}
+
+/**
+ * Gets the maximum q
+ */
+static PyObject * get_qmax(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.q_max);	
+}
+
 
 static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
@@ -513,4 +551,8 @@
 		   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
 		   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
+		   {"set_qmin", (PyCFunction)set_qmin, METH_VARARGS, ""},
+		   {"get_qmin", (PyCFunction)get_qmin, METH_VARARGS, ""},
+		   {"set_qmax", (PyCFunction)set_qmax, METH_VARARGS, ""},
+		   {"get_qmax", (PyCFunction)get_qmax, METH_VARARGS, ""},
 		   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
 		   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},
Index: pr_inversion/c_extensions/invertor.c
===================================================================
--- pr_inversion/c_extensions/invertor.c	(revision 4f63160d4aa8e41b4460f847aceb7324d62f7175)
+++ pr_inversion/c_extensions/invertor.c	(revision f71287f458148547d69aae7ddc900f1afde663f6)
@@ -18,4 +18,6 @@
 void invertor_init(Invertor_params *pars) {
 	pars->d_max = 180;
+	pars->q_min = -1.0;
+	pars->q_max = -1.0;
 }
 
Index: pr_inversion/c_extensions/invertor.h
===================================================================
--- pr_inversion/c_extensions/invertor.h	(revision 4f63160d4aa8e41b4460f847aceb7324d62f7175)
+++ pr_inversion/c_extensions/invertor.h	(revision f71287f458148547d69aae7ddc900f1afde663f6)
@@ -12,4 +12,6 @@
     int nerr;  
     double alpha;
+    double q_min;
+    double q_max;
 } Invertor_params; 
 
Index: pr_inversion/invertor.py
===================================================================
--- pr_inversion/invertor.py	(revision 4241e4821a36cddd1498d4438e88227a592f295f)
+++ pr_inversion/invertor.py	(revision f71287f458148547d69aae7ddc900f1afde663f6)
@@ -1,4 +1,5 @@
 from sans.pr.core.pr_inversion import Cinvertor
 import numpy
+import sys
 
 class Invertor(Cinvertor):
@@ -10,4 +11,10 @@
     ## Alpha to get the reg term the same size as the signal
     suggested_alpha = 0
+    ## Last number of base functions used
+    nfunc = 0
+    ## Last output values
+    out = None
+    ## Last errors on output values
+    cov = None
     
     def __init__(self):
@@ -30,4 +37,12 @@
         elif name=='d_max':
             return self.set_dmax(value)
+        elif name=='q_min':
+            if value==None:
+                return self.set_qmin(-1.0)
+            return self.set_qmin(value)
+        elif name=='q_max':
+            if value==None:
+                return self.set_qmax(-1.0)
+            return self.set_qmax(value)
         elif name=='alpha':
             return self.set_alpha(value)
@@ -56,4 +71,14 @@
         elif name=='d_max':
             return self.get_dmax()
+        elif name=='q_min':
+            qmin = self.get_qmin()
+            if qmin<0:
+                return None
+            return qmin
+        elif name=='q_max':
+            qmax = self.get_qmax()
+            if qmax<0:
+                return None
+            return qmax
         elif name=='alpha':
             return self.get_alpha()
@@ -71,4 +96,6 @@
         invertor.alpha   = self.alpha
         invertor.d_max   = self.d_max
+        invertor.q_min   = self.q_min
+        invertor.q_max   = self.q_max
         
         invertor.x = self.x
@@ -85,4 +112,5 @@
         import time
         
+        self.nfunc = nfunc
         # First, check that the current data is valid
         if self.is_valid()<=0:
@@ -147,8 +175,25 @@
         return self.get_pr_err(c, c_err, r)
        
+    def _accept_q(self, q):
+        """
+            Check q-value against user-defined range
+        """
+        if not self.q_min==None and q<self.q_min:
+            return False
+        if not self.q_max==None and q>self.q_max:
+            return False
+        return True
+       
     def lstsq(self, nfunc=5):
+        #TODO: do this on the C side
+        #
+        # To make sure an array is contiguous:
+        # blah = numpy.ascontiguousarray(blah_original)
+        # ... before passing it to C
+        
         import math
         from scipy.linalg.basic import lstsq
         
+        self.nfunc = nfunc
         # a -- An M x N matrix.
         # b -- An M x nrhs matrix or M vector.
@@ -163,5 +208,6 @@
         for j in range(nfunc):
             for i in range(npts):
-                a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
+                if self._accept_q(self.x[i]):
+                    a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
             for i_q in range(nq):
                 r = self.d_max/nq*i_q
@@ -170,5 +216,6 @@
         
         for i in range(npts):
-            b[i] = self.y[i]/self.err[i]
+            if self._accept_q(self.x[i]):
+                b[i] = self.y[i]/self.err[i]
             
         c, chi2, rank, n = lstsq(a, b)
@@ -181,5 +228,6 @@
                 inv_cov[i][j] = 0.0
                 for k in range(npts):
-                    inv_cov[i][j] = at[i][k]*a[k][j]
+                    if self._accept_q(self.x[i]):
+                        inv_cov[i][j] = at[i][k]*a[k][j]
                     
         # Compute the reg term size for the output
@@ -188,7 +236,8 @@
         for j in range(nfunc):
             for i in range(npts):
-                sum_sig += (a[i][j])**2
+                if self._accept_q(self.x[i]):
+                    sum_sig += (a[i][j])**2
             for i in range(nq):
-                sum_reg += (a[i_q+npts][j])**2
+                sum_reg += (a[i+npts][j])**2
                     
         if math.fabs(self.alpha)>0:
@@ -203,4 +252,7 @@
             print "Error estimating uncertainties"
             
+        # Keep a copy of the last output
+        self.out = c
+        self.cov = err
         
         return c, err
@@ -248,7 +300,170 @@
         return c, err
         
-        
-        
-        
+    def estimate_alpha(self, nfunc):
+        """
+            Returns a reasonable guess for the
+            regularization constant alpha
+            
+            @return: alpha, message, elapsed
+            
+            where alpha is the estimate for alpha,
+            message is a message for the user,
+            elapsed is the computation time
+        """
+        import time
+        try:            
+            pr = self.clone()
+            
+            # T_0 for computation time
+            starttime = time.time()
+            
+            # If the current alpha is zero, try
+            # another value
+            if pr.alpha<=0:
+                pr.alpha = 0.0001
+                 
+            # Perform inversion to find the largest alpha
+            out, cov = pr.lstsq(nfunc)
+            elapsed = time.time()-starttime
+            initial_alpha = pr.alpha
+            initial_peaks = pr.get_peaks(out)
+    
+            # Try the inversion with the estimated alpha
+            pr.alpha = pr.suggested_alpha
+            out, cov = pr.lstsq(nfunc)
+    
+            npeaks = pr.get_peaks(out)
+            # if more than one peak to start with
+            # just return the estimate
+            if npeaks>1:
+                message = "Your P(r) is not smooth, please check your inversion parameters"
+                return pr.suggested_alpha, message, elapsed
+            else:
+                
+                # Look at smaller values
+                # We assume that for the suggested alpha, we have 1 peak
+                # if not, send a message to change parameters
+                alpha = pr.suggested_alpha
+                best_alpha = pr.suggested_alpha
+                found = False
+                for i in range(10):
+                    pr.alpha = (0.33)**(i+1)*alpha
+                    out, cov = pr.lstsq(nfunc)
+                    
+                    peaks = pr.get_peaks(out)
+                    print pr.alpha, peaks
+                    if peaks>1:
+                        found = True
+                        break
+                    best_alpha = pr.alpha
+                    
+                # If we didn't find a turning point for alpha and
+                # the initial alpha already had only one peak,
+                # just return that
+                if not found and initial_peaks==1 and initial_alpha<best_alpha:
+                    best_alpha = initial_alpha
+                    
+                # Check whether the size makes sense
+                message=''
+                
+                if not found:
+                    message = "None"
+                elif best_alpha>=0.5*pr.suggested_alpha:
+                    # best alpha is too big, return a 
+                    # reasonable value
+                    message  = "The estimated alpha for your system is too large. "
+                    message += "Try increasing your maximum distance."
+                
+                return best_alpha, message, elapsed
+    
+        except:
+            message = "Invertor.estimate_alpha: %s" % sys.exc_value
+            return 0, message, elapsed
+    
+        
+    def to_file(self, path, npts=100):
+        """
+            Save the state to a file that will be readable
+            by SliceView.
+            @param path: path of the file to write
+            @param npts: number of P(r) points to be written
+        """
+        import pylab
+        
+        file = open(path, 'w')
+        file.write("#d_max=%g\n" % self.d_max)
+        file.write("#nfunc=%g\n" % self.nfunc)
+        file.write("#alpha=%g\n" % self.alpha)
+        file.write("#chi2=%g\n" % self.chi2)
+        file.write("#elapsed=%g\n" % self.elapsed)
+        file.write("#alpha_estimate=%g\n" % self.suggested_alpha)
+        if not self.out==None:
+            if len(self.out)==len(self.cov):
+                for i in range(len(self.out)):
+                    file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), str(self.cov[i][i])))
+        file.write("<r>  <Pr>  <dPr>\n")
+        r = pylab.arange(0.0, self.d_max, self.d_max/npts)
+        
+        for r_i in r:
+            (value, err) = self.pr_err(self.out, self.cov, r_i)
+            file.write("%g  %g  %g\n" % (r_i, value, err))
+    
+        file.close()
+        
+    def from_file(self, path):
+        """
+            Load the state of the Invertor from a file,
+            to be able to generate P(r) from a set of
+            parameters.
+            @param path: path of the file to load
+        """
+        import os
+        import re
+        if os.path.isfile(path):
+            try:
+                fd = open(path, 'r')
+                
+                buff    = fd.read()
+                lines   = buff.split('\n')
+                for line in lines:
+                    if line.startswith('#d_max='):
+                        toks = line.split('=')
+                        self.d_max = float(toks[1])
+                    elif line.startswith('#nfunc='):
+                        toks = line.split('=')
+                        self.nfunc = int(toks[1])
+                        self.out = numpy.zeros(self.nfunc)
+                        self.cov = numpy.zeros([self.nfunc, self.nfunc])
+                    elif line.startswith('#alpha='):
+                        toks = line.split('=')
+                        self.alpha = float(toks[1])
+                    elif line.startswith('#chi2='):
+                        toks = line.split('=')
+                        self.chi2 = float(toks[1])
+                    elif line.startswith('#elapsed='):
+                        toks = line.split('=')
+                        self.elapsed = float(toks[1])
+                    elif line.startswith('#alpha_estimate='):
+                        toks = line.split('=')
+                        self.suggested_alpha = float(toks[1])
+            
+                    # Now read in the parameters
+                    elif line.startswith('#C_'):
+                        toks = line.split('=')
+                        p = re.compile('#C_([0-9]+)')
+                        m = p.search(toks[0])
+                        toks2 = toks[1].split('+-')
+                        i = int(m.group(1))
+                        self.out[i] = float(toks2[0])
+                        
+                        self.cov[i][i] = float(toks2[1])                        
+            
+            except:
+                raise RuntimeError, "Invertor.from_file: corrupted file\n%s" % sys.exc_value
+        else:
+            raise RuntimeError, "Invertor.from_file: '%s' is not a file" % str(path) 
+        
+        
+    
     
 if __name__ == "__main__":
Index: pr_inversion/test/test_output.txt
===================================================================
--- pr_inversion/test/test_output.txt	(revision f71287f458148547d69aae7ddc900f1afde663f6)
+++ pr_inversion/test/test_output.txt	(revision f71287f458148547d69aae7ddc900f1afde663f6)
@@ -0,0 +1,117 @@
+#d_max=160
+#nfunc=10
+#alpha=0.0007
+#chi2=16654.1
+#elapsed=0
+#alpha_estimate=2.23431
+#C_0=1254.55002508+-12.0889464446
+#C_1=386.59602793+-62.4851961458
+#C_2=5.10066421813+-235.84586924
+#C_3=48.454288318+-1217.94687676
+#C_4=0.4174373785+-213577.449972
+#C_5=17.3551976035+-2673.36274352
+#C_6=-0.228001082457+-692.661759936
+#C_7=1.95934734382+-332.134754291
+#C_8=-0.885447495+-200.862052487
+#C_9=0.401866788245+-137.064680609
+<r>  <Pr>  <dPr>
+0  0  0
+1.6  236.371  235.678
+3.2  943.103  930.744
+4.8  2113.12  2049.82
+6.4  3734.86  3535.61
+8  5792.64  5311.37
+9.6  8267.09  7284.04
+11.2  11135.7  9348.1
+12.8  14373.5  11389.7
+14.4  17953.5  13291.4
+16  21847.7  14936.2
+17.6  26027.5  16212.7
+19.2  30464.2  17019.1
+20.8  35129.9  17267.8
+22.4  39997.6  16889.2
+24  45041.7  15836.1
+25.6  50238.2  14088.7
+27.2  55565.1  11662.9
+28.8  61002  8635.92
+30.4  66530.4  5281.73
+32  72133.5  3272.64
+33.6  77796  6067.26
+35.2  83503.6  10814.7
+36.8  89243.2  16028.2
+38.4  95001.9  21350.2
+40  100767  26561.2
+41.6  106527  31468
+43.2  112267  35887.3
+44.8  117975  39646.4
+46.4  123636  42586.2
+48  129233  44565.8
+49.6  134751  45467.6
+51.2  140171  45201.5
+52.8  145473  43709.8
+54.4  150638  40971.1
+56  155645  37005.5
+57.6  160473  31882.9
+59.2  165101  25742.1
+60.8  169508  18857.3
+62.4  173675  11970.1
+64  177583  8265.15
+65.6  181215  12750.9
+67.2  184556  21038.7
+68.8  187595  30104.6
+70.4  190321  39131.5
+72  192727  47705.7
+73.6  194808  55513.2
+75.2  196563  62282.7
+76.8  197991  67775.3
+78.4  199093  71784.5
+80  199873  74140.4
+81.6  200335  74714.1
+83.2  200484  73422.4
+84.8  200325  70232.5
+86.4  199862  65166
+88  199101  58304.9
+89.6  198044  49801.3
+91.2  196695  39903.5
+92.8  195056  29050.8
+94.4  193129  18346.5
+96  190915  12397.7
+97.6  188414  18724.8
+99.2  185629  30770.1
+100.8  182559  43834
+102.4  179206  56683.4
+104  175572  68726.9
+105.6  171661  79534.1
+107.2  167475  88745.8
+108.8  163019  96054.4
+110.4  158298  101203
+112  153320  103987
+113.6  148091  104262
+115.2  142620  101947
+116.8  136919  97026.9
+118.4  130998  89560.4
+120  124872  79680.9
+121.6  118558  67605.7
+123.2  112075  53655.8
+124.8  105448  38339.4
+126.4  98705.1  22821
+128  91877.8  13090.4
+129.6  85003.8  22520.4
+131.2  78125.5  39345.4
+132.8  71289.5  56946.4
+134.4  64546.5  73969.2
+136  57949.9  89741.1
+137.6  51554.5  103751
+139.2  45414.9  115564
+140.8  39583.1  124808
+142.4  34106.6  131176
+144  29025.7  134426
+145.6  24371.3  134390
+147.2  20162.8  130981
+148.8  16406  124194
+150.4  13091.6  114114
+152  10194.5  100912
+153.6  7673.45  84850.4
+155.2  5471.36  66272.7
+156.8  3516.48  45601.3
+158.4  1724.23  23326.7
Index: pr_inversion/test/utest_invertor.py
===================================================================
--- pr_inversion/test/utest_invertor.py	(revision 2d06beb109031c888260d5a73af0ae5e57628b01)
+++ pr_inversion/test/utest_invertor.py	(revision f71287f458148547d69aae7ddc900f1afde663f6)
@@ -333,4 +333,46 @@
         for i in range(len(self.x_in)):
             self.assertEqual(self.x_in[i], clone.x[i])
+        
+    def test_save(self):
+        x, y, err = load("sphere_80.txt")
+
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = .0007
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+        # Perform inversion
+        
+        out, cov = self.invertor.lstsq(10)
+        
+        # Save
+        self.invertor.to_file("test_output.txt")
+    
+    def test_load(self):
+        self.invertor.from_file("test_output.txt")
+        self.assertEqual(self.invertor.d_max, 160.0)
+        self.assertEqual(self.invertor.alpha, 0.0007)
+        self.assertEqual(self.invertor.chi2, 16654.1)
+        self.assertAlmostEqual(self.invertor.pr(self.invertor.out, 10.0), 8948.22689927, 4)
+        
+    def test_qmin(self):
+        self.invertor.q_min = 1.0
+        print self.invertor.q_min
+        self.assertEqual(self.invertor.q_min, 1.0)
+        
+        self.invertor.q_min = None
+        self.assertEqual(self.invertor.q_min, None)
+        
+                         
+    def test_qmax(self):
+        self.invertor.q_max = 1.0
+        self.assertEqual(self.invertor.q_max, 1.0)
+       
+        self.invertor.q_max = None
+        self.assertEqual(self.invertor.q_max, None)
+
         
 def pr_theory(r, R):
