Changes in sasmodels/direct_model.py [5024a56:2c4a190] in sasmodels

File:

• sasmodels/direct_model.py (modified) (11 diffs)

sasmodels/direct_model.py

@@ -32 +32 @@
 from .details import make_kernel_args, dispersion_mesh
 from .modelinfo import DEFAULT_BACKGROUND
-from .product import RADIUS_MODE_ID
 
 # pylint: disable=unused-import
@@ -65 +64 @@
     return calculator(call_details, values, cutoff, is_magnetic)
 
-def call_Fq(calculator, pars, cutoff=0., mono=False):
-    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray
-    """
-    Like :func:`call_kernel`, but returning F, F^2, R_eff, V_shell, V_form/V_shell.
-
-    For solid objects V_shell is equal to V_form and the volume ratio is 1.
-
-    Use parameter *radius_effective_mode* to select the effective radius
-    calculation.
-    """
-    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0))
-    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono)
-    #print("pars", list(zip(*mesh))[0])
-    call_details, values, is_magnetic = make_kernel_args(calculator, mesh)
-    #print("values:", values)
-    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type)
-
-def call_profile(model_info, pars=None):
-    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]]
+def call_ER(model_info, pars):
+    # type: (ModelInfo, ParameterSet) -> float
+    """
+    Call the model ER function using *values*.
+
+    *model_info* is either *model.info* if you have a loaded model,
+    or *kernel.info* if you have a model kernel prepared for evaluation.
+    """
+    if model_info.ER is None:
+        return 1.0
+    elif not model_info.parameters.form_volume_parameters:
+        # handle the case where ER is provided but model is not polydisperse
+        return model_info.ER()
+    else:
+        value, weight = _vol_pars(model_info, pars)
+        individual_radii = model_info.ER(*value)
+        return np.sum(weight*individual_radii) / np.sum(weight)
+
+
+def call_VR(model_info, pars):
+    # type: (ModelInfo, ParameterSet) -> float
+    """
+    Call the model VR function using *pars*.
+
+    *model_info* is either *model.info* if you have a loaded model,
+    or *kernel.info* if you have a model kernel prepared for evaluation.
+    """
+    if model_info.VR is None:
+        return 1.0
+    elif not model_info.parameters.form_volume_parameters:
+        # handle the case where ER is provided but model is not polydisperse
+        return model_info.VR()
+    else:
+        value, weight = _vol_pars(model_info, pars)
+        whole, part = model_info.VR(*value)
+        return np.sum(weight*part)/np.sum(weight*whole)
+
+
+def call_profile(model_info, **pars):
+    # type: (ModelInfo, ...) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]]
     """
     Returns the profile *x, y, (xlabel, ylabel)* representing the model.
     """
-    if pars is None:
-        pars = {}
     args = {}
     for p in model_info.parameters.kernel_parameters:
@@ -224 +242 @@
             else:
                 Iq, dIq = None, None
-            #self._theory = np.zeros_like(q)
-            q_vectors = [res.q_calc]
         elif self.data_type == 'Iqxy':
             #if not model.info.parameters.has_2d:
@@ -242 +258 @@
             res = resolution2d.Pinhole2D(data=data, index=index,
                                          nsigma=3.0, accuracy=accuracy)
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = res.q_calc
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -268 +282 @@
             else:
                 res = resolution.Perfect1D(data.x[index])
-
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = [res.q_calc]
         elif self.data_type == 'Iq-oriented':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -286 +297 @@
                                       qx_width=data.dxw[index],
                                       qy_width=data.dxl[index])
-            q_vectors = res.q_calc
         else:
             raise ValueError("Unknown data type") # never gets here
@@ -292 +302 @@
         # Remember function inputs so we can delay loading the function and
         # so we can save/restore state
-        self._kernel_inputs = q_vectors
         self._kernel = None
         self.Iq, self.dIq, self.index = Iq, dIq, index
@@ -329 +338 @@
         # type: (ParameterSet, float) -> np.ndarray
         if self._kernel is None:
-            self._kernel = self._model.make_kernel(self._kernel_inputs)
+            # TODO: change interfaces so that resolution returns kernel inputs
+            # Maybe have resolution always return a tuple, or maybe have
+            # make_kernel accept either an ndarray or a pair of ndarrays.
+            kernel_inputs = self.resolution.q_calc
+            if isinstance(kernel_inputs, np.ndarray):
+                kernel_inputs = (kernel_inputs,)
+            self._kernel = self._model.make_kernel(kernel_inputs)
 
         # Need to pull background out of resolution for multiple scattering
@@ -385 +400 @@
         Generate a plottable profile.
         """
-        return call_profile(self.model.info, pars)
+        return call_profile(self.model.info, **pars)
 
 def main():
@@ -401 +416 @@
     model_name = sys.argv[1]
     call = sys.argv[2].upper()
-    try:
-        values = [float(v) for v in call.split(',')]
-    except ValueError:
-        values = []
-    if len(values) == 1:
-        q, = values
-        data = empty_data1D([q])
-    elif len(values) == 2:
-        qx, qy = values
-        data = empty_data2D([qx], [qy])
-    else:
-        print("use q or qx,qy")
-        sys.exit(1)
+    if call != "ER_VR":
+        try:
+            values = [float(v) for v in call.split(',')]
+        except ValueError:
+            values = []
+        if len(values) == 1:
+            q, = values
+            data = empty_data1D([q])
+        elif len(values) == 2:
+            qx, qy = values
+            data = empty_data2D([qx], [qy])
+        else:
+            print("use q or qx,qy or ER or VR")
+            sys.exit(1)
+    else:
+        data = empty_data1D([0.001])  # Data not used in ER/VR
 
     model_info = load_model_info(model_name)
@@ -421 +439 @@
                 for pair in sys.argv[3:]
                 for k, v in [pair.split('=')])
-    Iq = calculator(**pars)
-    print(Iq[0])
+    if call == "ER_VR":
+        ER = call_ER(model_info, pars)
+        VR = call_VR(model_info, pars)
+        print(ER, VR)
+    else:
+        Iq = calculator(**pars)
+        print(Iq[0])
 
 if __name__ == "__main__":