Changes in / [f9c8474:a0243f7] in sasmodels


Ignore:
Files:
2 added
1 deleted
7 edited

Legend:

Unmodified
Added
Removed
  • example/fit.py

    rf4b36fa r1ee5ac6  
    179179 
    180180else: 
    181     print "No parameters for %s"%name 
     181    print("No parameters for %s"%name) 
    182182    sys.exit(1) 
    183183 
     
    192192else: 
    193193   problem = FitProblem(M) 
     194 
     195if __name__ == "__main__": 
     196   problem.plot() 
     197   import pylab; pylab.show() 
  • sasmodels/bumps_model.py

    r2d81cfe r49d1f8b8  
    137137    """ 
    138138    _cache = None # type: Dict[str, np.ndarray] 
    139     def __init__(self, data, model, cutoff=1e-5, name=None): 
     139    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140140        # type: (Data, Model, float) -> None 
    141141        # remember inputs so we can inspect from outside 
     
    145145        self._interpret_data(data, model.sasmodel) 
    146146        self._cache = {} 
     147        self.extra_pars = extra_pars 
    147148 
    148149    def update(self): 
     
    166167        Return a dictionary of parameters 
    167168        """ 
    168         return self.model.parameters() 
     169        pars = self.model.parameters() 
     170        if self.extra_pars: 
     171            pars.update(self.extra_pars) 
     172        return pars 
    169173 
    170174    def theory(self): 
  • sasmodels/data.py

    r2d81cfe rf549e37  
    706706    else: 
    707707        vmin_scaled, vmax_scaled = vmin, vmax 
    708     plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
     708    nx, ny = len(data.x_bins), len(data.y_bins) 
     709    x_bins, y_bins, image = _build_matrix(data, plottable) 
     710    plt.imshow(image, 
    709711               interpolation='nearest', aspect=1, origin='lower', 
    710712               extent=[xmin, xmax, ymin, ymax], 
     
    714716    return vmin, vmax 
    715717 
     718 
     719# === The following is modified from sas.sasgui.plottools.PlotPanel 
     720def _build_matrix(self, plottable): 
     721    """ 
     722    Build a matrix for 2d plot from a vector 
     723    Returns a matrix (image) with ~ square binning 
     724    Requirement: need 1d array formats of 
     725    self.data, self.qx_data, and self.qy_data 
     726    where each one corresponds to z, x, or y axis values 
     727 
     728    """ 
     729    # No qx or qy given in a vector format 
     730    if self.qx_data is None or self.qy_data is None \ 
     731            or self.qx_data.ndim != 1 or self.qy_data.ndim != 1: 
     732        return self.x_bins, self.y_bins, plottable 
     733 
     734    # maximum # of loops to fillup_pixels 
     735    # otherwise, loop could never stop depending on data 
     736    max_loop = 1 
     737    # get the x and y_bin arrays. 
     738    x_bins, y_bins = _get_bins(self) 
     739    # set zero to None 
     740 
     741    #Note: Can not use scipy.interpolate.Rbf: 
     742    # 'cause too many data points (>10000)<=JHC. 
     743    # 1d array to use for weighting the data point averaging 
     744    #when they fall into a same bin. 
     745    weights_data = np.ones([self.data.size]) 
     746    # get histogram of ones w/len(data); this will provide 
     747    #the weights of data on each bins 
     748    weights, xedges, yedges = np.histogram2d(x=self.qy_data, 
     749                                             y=self.qx_data, 
     750                                             bins=[y_bins, x_bins], 
     751                                             weights=weights_data) 
     752    # get histogram of data, all points into a bin in a way of summing 
     753    image, xedges, yedges = np.histogram2d(x=self.qy_data, 
     754                                           y=self.qx_data, 
     755                                           bins=[y_bins, x_bins], 
     756                                           weights=plottable) 
     757    # Now, normalize the image by weights only for weights>1: 
     758    # If weight == 1, there is only one data point in the bin so 
     759    # that no normalization is required. 
     760    image[weights > 1] = image[weights > 1] / weights[weights > 1] 
     761    # Set image bins w/o a data point (weight==0) as None (was set to zero 
     762    # by histogram2d.) 
     763    image[weights == 0] = None 
     764 
     765    # Fill empty bins with 8 nearest neighbors only when at least 
     766    #one None point exists 
     767    loop = 0 
     768 
     769    # do while loop until all vacant bins are filled up up 
     770    #to loop = max_loop 
     771    while (weights == 0).any(): 
     772        if loop >= max_loop:  # this protects never-ending loop 
     773            break 
     774        image = _fillup_pixels(self, image=image, weights=weights) 
     775        loop += 1 
     776 
     777    return x_bins, y_bins, image 
     778 
     779def _get_bins(self): 
     780    """ 
     781    get bins 
     782    set x_bins and y_bins into self, 1d arrays of the index with 
     783    ~ square binning 
     784    Requirement: need 1d array formats of 
     785    self.qx_data, and self.qy_data 
     786    where each one corresponds to  x, or y axis values 
     787    """ 
     788    # find max and min values of qx and qy 
     789    xmax = self.qx_data.max() 
     790    xmin = self.qx_data.min() 
     791    ymax = self.qy_data.max() 
     792    ymin = self.qy_data.min() 
     793 
     794    # calculate the range of qx and qy: this way, it is a little 
     795    # more independent 
     796    x_size = xmax - xmin 
     797    y_size = ymax - ymin 
     798 
     799    # estimate the # of pixels on each axes 
     800    npix_y = int(np.floor(np.sqrt(len(self.qy_data)))) 
     801    npix_x = int(np.floor(len(self.qy_data) / npix_y)) 
     802 
     803    # bin size: x- & y-directions 
     804    xstep = x_size / (npix_x - 1) 
     805    ystep = y_size / (npix_y - 1) 
     806 
     807    # max and min taking account of the bin sizes 
     808    xmax = xmax + xstep / 2.0 
     809    xmin = xmin - xstep / 2.0 
     810    ymax = ymax + ystep / 2.0 
     811    ymin = ymin - ystep / 2.0 
     812 
     813    # store x and y bin centers in q space 
     814    x_bins = np.linspace(xmin, xmax, npix_x) 
     815    y_bins = np.linspace(ymin, ymax, npix_y) 
     816 
     817    return x_bins, y_bins 
     818 
     819def _fillup_pixels(self, image=None, weights=None): 
     820    """ 
     821    Fill z values of the empty cells of 2d image matrix 
     822    with the average over up-to next nearest neighbor points 
     823 
     824    :param image: (2d matrix with some zi = None) 
     825 
     826    :return: image (2d array ) 
     827 
     828    :TODO: Find better way to do for-loop below 
     829 
     830    """ 
     831    # No image matrix given 
     832    if image is None or np.ndim(image) != 2 \ 
     833            or np.isfinite(image).all() \ 
     834            or weights is None: 
     835        return image 
     836    # Get bin size in y and x directions 
     837    len_y = len(image) 
     838    len_x = len(image[1]) 
     839    temp_image = np.zeros([len_y, len_x]) 
     840    weit = np.zeros([len_y, len_x]) 
     841    # do for-loop for all pixels 
     842    for n_y in range(len(image)): 
     843        for n_x in range(len(image[1])): 
     844            # find only null pixels 
     845            if weights[n_y][n_x] > 0 or np.isfinite(image[n_y][n_x]): 
     846                continue 
     847            else: 
     848                # find 4 nearest neighbors 
     849                # check where or not it is at the corner 
     850                if n_y != 0 and np.isfinite(image[n_y - 1][n_x]): 
     851                    temp_image[n_y][n_x] += image[n_y - 1][n_x] 
     852                    weit[n_y][n_x] += 1 
     853                if n_x != 0 and np.isfinite(image[n_y][n_x - 1]): 
     854                    temp_image[n_y][n_x] += image[n_y][n_x - 1] 
     855                    weit[n_y][n_x] += 1 
     856                if n_y != len_y - 1 and np.isfinite(image[n_y + 1][n_x]): 
     857                    temp_image[n_y][n_x] += image[n_y + 1][n_x] 
     858                    weit[n_y][n_x] += 1 
     859                if n_x != len_x - 1 and np.isfinite(image[n_y][n_x + 1]): 
     860                    temp_image[n_y][n_x] += image[n_y][n_x + 1] 
     861                    weit[n_y][n_x] += 1 
     862                # go 4 next nearest neighbors when no non-zero 
     863                # neighbor exists 
     864                if n_y != 0 and n_x != 0 and \ 
     865                        np.isfinite(image[n_y - 1][n_x - 1]): 
     866                    temp_image[n_y][n_x] += image[n_y - 1][n_x - 1] 
     867                    weit[n_y][n_x] += 1 
     868                if n_y != len_y - 1 and n_x != 0 and \ 
     869                        np.isfinite(image[n_y + 1][n_x - 1]): 
     870                    temp_image[n_y][n_x] += image[n_y + 1][n_x - 1] 
     871                    weit[n_y][n_x] += 1 
     872                if n_y != len_y and n_x != len_x - 1 and \ 
     873                        np.isfinite(image[n_y - 1][n_x + 1]): 
     874                    temp_image[n_y][n_x] += image[n_y - 1][n_x + 1] 
     875                    weit[n_y][n_x] += 1 
     876                if n_y != len_y - 1 and n_x != len_x - 1 and \ 
     877                        np.isfinite(image[n_y + 1][n_x + 1]): 
     878                    temp_image[n_y][n_x] += image[n_y + 1][n_x + 1] 
     879                    weit[n_y][n_x] += 1 
     880 
     881    # get it normalized 
     882    ind = (weit > 0) 
     883    image[ind] = temp_image[ind] / weit[ind] 
     884 
     885    return image 
     886 
     887 
    716888def demo(): 
    717889    # type: () -> None 
  • sasmodels/direct_model.py

    r2d81cfe rd18d6dd  
    345345            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    346346 
     347        # Need to pull background out of resolution for multiple scattering 
     348        background = pars.get('background', 0.) 
     349        pars = pars.copy() 
     350        pars['background'] = 0. 
     351 
    347352        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    348353        # Storing the calculated Iq values so that they can be plotted. 
     
    357362                np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
    358363            ) 
    359         return result 
     364        return result + background 
    360365 
    361366 
  • sasmodels/generate.py

    rf9c8474 rf9c8474  
    215215 
    216216EXTERNAL_DIR = 'sasmodels-data' 
    217 DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_template.c') 
     217DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c') 
    218218MODEL_PATH = joinpath(DATA_PATH, 'models') 
    219219 
     
    390390    # TODO: fails DRY; templates appear two places. 
    391391    model_templates = [joinpath(DATA_PATH, filename) 
    392                        for filename in ('kernel_header.c', 'kernel_iq.cl')] 
     392                       for filename in ('kernel_header.c', 'kernel_iq.c')] 
    393393    source_files = (model_sources(model_info) 
    394394                    + model_templates 
  • sasmodels/kernelcl.py

    rb4272a2 r6cbdcd4  
    208208    Build a model to run on the gpu. 
    209209 
    210     Returns the compiled program and its type.  The returned type will 
    211     be float32 even if the desired type is float64 if any of the 
    212     devices in the context do not support the cl_khr_fp64 extension. 
     210    Returns the compiled program and its type. 
     211 
     212    Raises an error if the desired precision is not available. 
    213213    """ 
    214214    dtype = np.dtype(dtype) 
  • sasmodels/kerneldll.py

    r2d81cfe rbf94e6e  
    158158        return CC + [source, "-o", output, "-lm"] 
    159159 
    160 # Windows-specific solution 
    161 if os.name == 'nt': 
    162     # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    163     DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    164     if not os.path.exists(DLL_PATH): 
    165         os.makedirs(DLL_PATH) 
    166 else: 
    167     # Set up the default path for compiled modules. 
    168     DLL_PATH = tempfile.gettempdir() 
     160# Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     161DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    169162 
    170163ALLOW_SINGLE_PRECISION_DLLS = True 
     
    233226 
    234227    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
    235     The default is the system temporary directory. 
     228    The default is in ~/.sasmodels/compiled_models. 
    236229    """ 
    237230    if dtype == F16: 
     
    250243        need_recompile = dll_time < newest_source 
    251244    if need_recompile: 
     245        # Make sure the DLL path exists 
     246        if not os.path.exists(DLL_PATH): 
     247            os.makedirs(DLL_PATH) 
    252248        basename = splitext(os.path.basename(dll))[0] + "_" 
    253249        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
Note: See TracChangeset for help on using the changeset viewer.