Changeset c094758 in sasmodels

Timestamp:

Mar 18, 2016 12:14:07 AM (9 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

3a45c2c, a5af4e1

Parents:

cbd37a7

Message:

fix autogen 2d plots (mask passes 0 not 1)

Files:

• doc/genmodel.py (modified) (2 diffs)
• sasmodels/data.py (modified) (1 diff)
• sasmodels/kernelcl.py (modified) (1 diff)

doc/genmodel.py

@@ -10 +10 @@
 # Convert ../sasmodels/models/name.py to name
 model_name = os.path.basename(sys.argv[1])[:-3]
+model_info = core.load_model_info(model_name)
+model = core.build_model(model_info)
 
 # Load the doc string from the module definition file and store it in rst
-docstr = generate.make_doc(core.load_model_info(model_name))
-    
-# Generate automatically plot of the model and add it to rst documentation
+docstr = generate.make_doc(model_info)
 
-info = core.load_model_info(model_name)
 
 # Calculate 1D curve for default parameters
-pars = dict((p[0], p[2]) for p in info['parameters'])
+pars = dict((p.name, p.default) for p in model_info['parameters'])
 
 # Plotting ranges and options
 opts = {
-        'xscale'    : 'log',
-        'yscale'    : 'log' if not info['structure_factor'] else 'linear',
-        'qmin'      : 0.001,
-        'qmax'      : 1.0,
-        'nq'        : 1000,
-        'nq2d'      : 100,
+    'xscale'    : 'log',
+    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear',
+    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear',
+    'q_min'     : 0.001,
+    'q_max'     : 1.0,
+    'nq'        : 1000,
+    'nq2d'      : 400,
+    'vmin'      : 1e-3,  # floor for the 2D data results
+    'qx_max'    : 0.5,
 }
 
-qmin, qmax, nq = opts['qmin'], opts['qmax'], opts['nq']
-qmin = math.log10(qmin)
-qmax = math.log10(qmax)
-q = np.logspace(qmin, qmax, nq)
-data = empty_data1D(q)
-model = core.load_model(model_name)
-calculator = DirectModel(data, model)
-Iq1D = calculator()
 
-# TO DO: Generation of 2D plots 
-# Problem in sasmodels.direct_model._calc_theory
-# There self._kernel.q_input.nq  gets a value of 0 in the 2D case
-# and returns a 0 numpy array (it does not call the C code)
+def plot_1d(model, opts, ax):
+    q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq']
+    q_min = math.log10(q_min)
+    q_max = math.log10(q_max)
+    q = np.logspace(q_min, q_max, nq)
+    data = empty_data1D(q)
+    calculator = DirectModel(data, model)
+    Iq1D = calculator()
 
-# If 2D model, compute 2D image
-#if info['has_2d'] != []:
-#    qmax, nq2d = opts['qmax'], opts['nq2d']
-#    data2d = empty_data2D(np.linspace(-qmax, qmax, nq2d), resolution=0.0) 
-#    #model = core.load_model(model_name)
-#    calculator = DirectModel(data2d, model)
-#    Iq2D = calculator()
+    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name'])
+    ax.set_xlabel(r'$Q \/(\AA^{-1})$')
+    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
+    ax.set_xscale(opts['xscale'])
+    ax.set_yscale(opts['yscale'])
+    #ax.legend(loc='best')
+
+def plot_2d(model, opts, ax):
+    qx_max, nq2d = opts['qx_max'], opts['nq2d']
+    q = np.linspace(-qx_max, qx_max, nq2d)
+    data2d = empty_data2D(q, resolution=0.0)
+    calculator = DirectModel(data2d, model)
+    Iq2D = calculator() #background=0)
+    Iq2D = Iq2D.reshape(nq2d, nq2d)
+    if opts['zscale'] == 'log':
+        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf))
+    h = ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='upper',
+           extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=ice_cm())
+    # , vmin=vmin, vmax=vmax)
+    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$')
+    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$')
+
+def ice_cm():
+    from matplotlib._cm import _Blues_data
+    from matplotlib import colors
+    from matplotlib import rcParams
+    def from_white(segments):
+        scale = 1.0/segments[0][1]
+        return [(k, v*scale, w*scale) for k, v, w in segments]
+    ice_data = dict((k,from_white(v)) for k,v in _Blues_data.items())
+    ice = colors.LinearSegmentedColormap("ice", ice_data, rcParams['image.lut'])
+    return ice
+
 
 # Generate image (comment IF for 1D/2D for the moment) and generate only 1D
-#if info['has_2d'] == []:
-#    fig = plt.figure()
-#    ax = fig.add_subplot(1,1,1)
-#    ax.plot(q, Iq1D, color='blue', lw=2, label=model_name)
-#    ax.set_xlabel(r'$Q \/(\AA^{-1})$')
-#    ax.set_xscale(opts['xscale'])   
-#    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
-#    ax.set_yscale(opts['yscale'])  
-#    ax.legend()
-#else:
-#    # need figure with 1D + 2D
-#    pass
-fig = plt.figure()
-ax = fig.add_subplot(1,1,1)
-ax.plot(q, Iq1D, color='blue', lw=2, label=info['name'])
-ax.set_xlabel(r'$Q \/(\AA^{-1})$')
-ax.set_xscale(opts['xscale'])   
-ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
-ax.set_yscale(opts['yscale'])  
-ax.legend()
- 
+fig_height = 3.0 # in
+fig_left = 0.6 # in
+fig_right = 0.5 # in
+fig_top = 0.6*0.25 # in
+fig_bottom = 0.6*0.75
+if model_info['has_2d']:
+    plot_height = fig_height - (fig_top+fig_bottom)
+    plot_width = plot_height
+    fig_width = 2*(plot_width + fig_left + fig_right)
+    aspect = (fig_width, fig_height)
+    ratio = aspect[0]/aspect[1]
+    ax_left = fig_left/fig_width
+    ax_bottom = fig_bottom/fig_height
+    ax_height = plot_height/fig_height
+    ax_width = ax_height/ratio # square axes
+    fig = plt.figure(figsize=aspect)
+    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height])
+    plot_2d(model, opts, ax2d)
+    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+    #ax.set_aspect('square')
+else:
+    plot_height = fig_height - (fig_top+fig_bottom)
+    plot_width = (1+np.sqrt(5))/2*fig_height
+    fig_width = plot_width + fig_left + fig_right
+    ax_left = fig_left/fig_width
+    ax_bottom = fig_bottom/fig_height
+    ax_width = plot_width/fig_width
+    ax_height = plot_height/fig_height
+    aspect = (fig_width, fig_height)
+    fig = plt.figure(figsize=aspect)
+    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+plot_1d(model, opts, ax1d)
 
 # Save image in model/img
 figname = model_name + '_autogenfig.png'
 filename = os.path.join('model', 'img', figname)
-plt.savefig(filename)
+plt.savefig(filename, bbox_inches='tight')
+#print "figure saved in",filename
 
 # Auto caption for figure
 captionstr = '\n'
-captionstr += '.. figure:: img/' + model_name + '_autogenfig.png\n'
+captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n'
 captionstr += '\n'
-#if info['has_2d'] == []:
-#    captionstr += '    1D plot corresponding to the default parameters of the model.\n'
-#else:
-#    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n'
 captionstr += '    1D plot corresponding to the default parameters of the model.\n'
 captionstr += '\n'
@@ -102 +134 @@
 else:
     print '------------------------------------------------------------------'
-    print 'References NOT FOUND for model: ', model_name
+    print 'References NOT FOUND for model: ', model_info['id']
     print '------------------------------------------------------------------'
     docstr = docstr + captionstr

sasmodels/data.py

@@ -178 +178 @@
         self.qy_data, self.dqy_data = y, dy
         self.data, self.err_data = z, dz
-        self.mask = (~np.isnan(z) if z is not None
-                     else np.ones_like(x) if x is not None
+        self.mask = (np.isnan(z) if z is not None
+                     else np.zeros_like(x, dtype='bool') if x is not None
                      else None)
         self.q_data = np.sqrt(x**2 + y**2)

sasmodels/kernelcl.py

@@ -367 +367 @@
         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
         context = env.get_context(self.dtype)
+        self.global_size = [self.q_vectors[0].size]
+        #print("creating inputs of size", self.global_size)
         self.q_buffers = [
             cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
             for q in self.q_vectors
         ]
-        self.global_size = [self.q_vectors[0].size]
 
     def release(self):