Ignore:
Timestamp:
May 8, 2008 7:20:22 PM (17 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
f71287f4
Parents:
4f63160
Message:

Added alpha estimator and help

File:
1 edited

Legend:

Unmodified
Added
Removed

  • prview/perspectives/pr/pr_thread.py

    rf3d51f6 r32dffae4  
    7575        except KeyboardInterrupt: 
    7676            printEVT("P(r) calc interrupted") 
    77             raise KeyboardInterrupt 
     77            raise KeyboardInterrupt     
    7878         
    7979    def compute(self): 
    8080        import time 
    81         try: 
    82              
    83             print "Alpha  Oscill" 
    84              
     81        try:             
    8582            self.starttime = time.time() 
    8683            # If the current alpha is zero, try 
    … …  
    8986                self.pr.alpha = 0.0001 
    9087                  
     88            # Perform inversion to find the largest alpha 
    9189            out, cov = self.pr.lstsq(self.nfunc) 
    9290            elapsed = time.time()-self.starttime 
    93              
     91            initial_alpha = self.pr.alpha 
     92            initial_peaks = self.pr.get_peaks(out) 
    9493 
    95             # Take the default and try to find  
    96             # a better value 
    97             best_alpha = self.pr.alpha 
    98             best_osc   = self.pr.oscillations(out)  
    99              
    100             print best_alpha, best_osc 
    101              
    102             alpha = self.pr.suggested_alpha 
    103             print "initial:", alpha 
     94            # Try the inversion with the estimated alpha 
     95            self.pr.alpha = self.pr.suggested_alpha 
     96            out, cov = self.pr.lstsq(self.nfunc) 
    10497 
    105             # Look at smaller values 
    106             for i in range(5): 
    107                 self.pr.alpha = (0.1)**(i)*alpha 
    108                 out, cov = self.pr.lstsq(self.nfunc) 
    109                 osc = self.pr.oscillations(out)  
    110                 print self.pr.alpha, osc 
    111                 if  osc < best_osc: 
    112                     best_osc = osc 
    113                     best_alpha = alpha 
    114              
    115             ## Look at larger values 
    116             #for i in range(4): 
    117             #    self.pr.alpha = (10.0)**(i+1)*alpha 
    118             #    out, cov = self.pr.lstsq(self.nfunc) 
    119             #    osc = self.pr.oscillations(out)  
    120             #    print self.pr.alpha, osc 
    121             #    if  osc < best_osc: 
    122             #        best_osc = osc 
    123             #        best_alpha = alpha 
    124              
    125              
    126             self.complete(alpha=best_alpha, elapsed=elapsed) 
    127              
    128              
    129              
    130              
     98            npeaks = self.pr.get_peaks(out) 
     99            # if more than one peak to start with 
     100            # just return the estimate 
     101            if npeaks>1: 
     102                message = "Your P(r) is not smooth, please check your inversion parameters" 
     103                self.complete(alpha=self.pr.suggested_alpha, message=message, elapsed=elapsed) 
     104            else: 
     105                 
     106                # Look at smaller values 
     107                # We assume that for the suggested alpha, we have 1 peak 
     108                # if not, send a message to change parameters 
     109                alpha = self.pr.suggested_alpha 
     110                best_alpha = self.pr.suggested_alpha 
     111                found = False 
     112                for i in range(10): 
     113                    self.pr.alpha = (0.33)**(i+1)*alpha 
     114                    out, cov = self.pr.lstsq(self.nfunc) 
     115                    #osc = self.pr.oscillations(out)  
     116                    #print self.pr.alpha, osc 
     117                     
     118                    peaks = self.pr.get_peaks(out) 
     119                    print self.pr.alpha, peaks 
     120                    if peaks>1: 
     121                        found = True 
     122                        break 
     123                    best_alpha = self.pr.alpha 
     124                     
     125                # If we didn't find a turning point for alpha and 
     126                # the initial alpha already had only one peak, 
     127                # just return that 
     128                if not found and initial_peaks==1 and initial_alpha<best_alpha: 
     129                    best_alpha = initial_alpha 
     130                     
     131                # Check whether the size makes sense 
     132                message=None 
     133                 
     134                if not found: 
     135                    message = "None" 
     136                elif best_alpha>=0.5*self.pr.suggested_alpha: 
     137                    # best alpha is too big, return a  
     138                    # reasonable value 
     139                    message  = "The estimated alpha for your system is too large. " 
     140                    message += "Try increasing your maximum distance." 
     141                 
     142                self.complete(alpha=best_alpha, message=None, elapsed=elapsed) 
     143 
    131144        except: 
    132145            if not self.error_func==None: 
Note: See TracChangeset for help on using the changeset viewer.