Changeset 560 for branches


Ignore:
Timestamp:
09/14/12 16:42:47 (4 years ago)
Author:
mmckerns
Message:

update Lan's collapse code to as used in SURF final report;
add Lan's semilog plotter script

Location:
branches/collapse
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • branches/collapse/TEST_OUQ_surrogate_diam.py

    r542 r560  
    11#!/usr/bin/env python 
    22 
    3 debug = False 
     3debug = True #False 
    44MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    55####################################################################### 
     
    88####################################################################### 
    99npop = 40 
    10 maxiter = 1000 
     10maxiter = 6000 
    1111maxfun = 1e+6 
    12 convergence_tol = 1e-6; ngen = 40 
    1312crossover = 0.9 
    1413percent_change = 0.9 
     14convergence_tol = 1e-6 
     15tolW = 0.05; tolP_range_fr = 0.05 
     16ngen = 200; ngcol = 200 
    1517 
    1618 
     
    1820# the model function 
    1921####################################################################### 
    20 from surrogate import marc_surr as model 
    21  
     22from surrogate import marc_surr as model0 
    2223 
    2324####################################################################### 
     
    2627def optimize(cost,_bounds,_constraints): 
    2728  from mystic.solvers import DifferentialEvolutionSolver2 
    28  #from mystic.termination import ChangeOverGeneration as COG 
    29   from mystic.termination import Or 
    30   from collapse_termination import CLPS_Weight as WCLPS 
    31   from collapse_termination import CLPS_Location as LCLPS 
    3229  from mystic.strategy import Best1Exp 
    3330  from mystic.monitors import VerboseMonitor, Monitor 
    3431  from mystic.tools import random_seed 
    3532 
    36  #random_seed(123) 
    37  
    38   stepmon = VerboseMonitor(2) 
     33  random_seed(122) 
     34 
     35  stepmon = VerboseMonitor(10,10) 
    3936 #stepmon = Monitor() 
    4037  evalmon = Monitor() 
     
    5047  solver.SetGenerationMonitor(stepmon) 
    5148  solver.SetConstraints(_constraints) 
    52  
     49  
     50  ############################################################################## 
     51   
     52                 ## THE FOLLOWING TERMINATION AND CONSTRAINTS ## 
     53                 ## CAN BE MOVED TO A NEW MODULE AND IMPORTED ## 
     54   
     55  ######################### TERMINATION CONDITIONS ############################# 
     56  from mystic.termination import Or   
     57  from mystic.termination import ChangeOverGeneration as COG 
     58  from collapse_termination_v2 import CLPS_Weight as WCLPS 
     59  from collapse_termination_v2 import CLPS_Position as PCLPS   
     60  from mystic.math.measures import _nested_split   
     61   
    5362  tol = convergence_tol 
    54   _WCLPS = WCLPS(npts, tolG=tol, generations=ngen) 
    55   _LCLPS = LCLPS(npts, tolL=[1,1,0.03], tolG=tol, generations=ngen) #FIXME 
    56   CLPS = Or( _WCLPS, _LCLPS ) 
    57   solver.Solve(cost,termination=CLPS,strategy=Best1Exp, \ 
    58                CrossProbability=crossover,ScalingFactor=percent_change) 
    59  
     63  _COG = COG(tol,ngen) 
     64   
     65  w_lb, x_lb = _nested_split(lb, npts) 
     66  w_ub, x_ub = _nested_split(ub, npts) 
     67  range_x = [] 
     68  for i in range(len(x_lb)): 
     69      range_x.append(abs(x_lb[i][0] - x_ub[i][0]))   
     70  tolP = [ele*tolP_range_fr for ele in range_x]  
     71  if debug: print 'tolP = %s' %tolP 
     72   
     73  indexW_ids = [range(nx) for nx in npts]  # list of indices for CLPS_Weight termination check 
     74  indexP_ids = []  # list of indices for CLPS_Position termination check   
     75   
     76  x = [] 
     77  _WCLPS = [] 
     78  _PCLPS = [] 
     79   
     80  for m in range(len(indexW_ids)): 
     81      for idx in indexW_ids[m]: 
     82          _WCLPS.append(WCLPS(npts, tol_weight = tolW, generations=ngcol, axisW = m, indexW = idx)) 
     83           
     84  for m in range(len(indexP_ids)): 
     85      for pair in indexP_ids[m]: 
     86          _PCLPS.append( PCLPS(npts, tol_position=tolP, generations=ngcol, \ 
     87                                axisP = m, index1 = pair[0], index2  = pair[1]))  
     88  x.extend(_WCLPS) 
     89  x.extend(_PCLPS) 
     90  x.append(_COG) 
     91  trm = Or(*x) 
     92   
     93  ######################## NEW CONSTRAINTS AFTER COLLAPSE ###################### 
     94   
     95  list_vanished_weights = [] # list of [measure, index] lists 
     96  list_collapsed_pairs = [] # list of [measure, pair] lists, where pair = [index1, index2]   
     97   
     98  for i in range(len(npts)): 
     99    c = [] 
     100    for idx in range(npts[i]): 
     101        for jdx in range(idx,npts[i]): 
     102            if idx !=jdx: 
     103                c.append([idx, jdx]) 
     104    indexP_ids.append(c)   
     105   
     106  def _set_weight_to_zero(rv, npts, axisW, indexW): 
     107    from mystic.math.dirac_measure import product_measure 
     108    c = product_measure() 
     109    c.load(rv, npts) 
     110    c[axisW][indexW].weight = 0 
     111    return c.flatten() 
     112   
     113  def _set_same_positions(rv, npts, axisP, index1, index2): 
     114    from mystic.math.dirac_measure import product_measure 
     115    c = product_measure() 
     116    c.load(rv, npts) 
     117    c[axisP][index1].weight += c[axisP][index2].weight 
     118    c[axisP][index2].weight = 0 
     119    c[axisP][index2].position = c[axisP][index1].position 
     120    return c.flatten()  
     121   
     122  def new_constraints(rv, list_collapsed_pairs, list_vanished_weights): 
     123    for ele in list_vanished_weights: 
     124      rv = _set_weight_to_zero(rv, npts, axisW = ele[0], indexW = ele[1]) 
     125    for mem in list_collapsed_pairs: 
     126      rv = _set_same_positions(rv, npts, axisP = mem[0], \ 
     127              index1 = mem[1][0], index2 = mem[1][1])             
     128    return _constraints(rv)    
     129   
     130  ############################################################################## 
     131   
     132   
     133  ############################ RUN THE SOLVER ################################## 
     134  collapse_data = {} 
     135  tot_evaluations = 0 # total number of function eveluations 
     136  while _COG(solver) != True:          
     137    solver.Solve(cost,termination=trm,strategy=Best1Exp, \ 
     138               CrossProbability=crossover,ScalingFactor=percent_change)  
     139    tot_evaluations += solver.evaluations     
     140    collapse_data[stepmon._step] = trm(solver, True) 
     141     
     142    if debug: 
     143      print 'last_itr_evaluations = %s' %solver.evaluations 
     144      print "Generation  %d " %stepmon._step  
     145     
     146    int_trm = trm(solver, info = 'self')   
     147    for ele in range(len(int_trm)): 
     148      cond = int_trm[ele] 
     149      string = cond(solver, info = True) 
     150      string1 = string.replace(',', ' ') 
     151      string2 = string1.replace('&', ' ') 
     152      string3 = string2.replace(':', ' ') 
     153      lst = string3.split() 
     154       
     155      if lst[0] == 'CLPS_Weight': 
     156          for i in range(len(lst)): 
     157              if lst[i] == 'measure': 
     158                  axisW = int(lst[i+1]) 
     159                  if lst[i+2] == 'index': indexW = int(lst[i+3]) 
     160                   
     161          if debug: 
     162            print 'axisW = %s' %axisW 
     163            print 'indexW = %s' %indexW 
     164           
     165          list_vanished_weights.append([axisW, indexW])  
     166          indexW_ids[axisW].remove(indexW) 
     167                  
     168      if lst[0] == 'CLPS_Position': 
     169          for i in range(len(lst)):   
     170              if lst[i] == 'measure':         
     171                  axisP = int(lst[i+1]) 
     172                  if lst[i+2] == 'indices':  
     173                      index1 = int(lst[i+3].lstrip('(')) 
     174                      index2 = int(lst[i+4].rstrip(')')) 
     175          if debug: 
     176            print "axisP = %s" %axisP 
     177            print "index1 = %s" %index1 
     178            print "index2 = %s" %index2 
     179           
     180          list_collapsed_pairs.append([axisP, [index1, index2]]) 
     181          indexP_ids[axisP].remove([index1,index2])           
     182       
     183      ##### AFTER COLLAPSE OCCURRED CHANGE THE TERMINATION AND CONSTRAINTS ##### 
     184                     
     185      def _new_constraints(rv): 
     186        return new_constraints(rv, list_collapsed_pairs, list_vanished_weights)       
     187      solver.SetConstraints(_new_constraints) 
     188       
     189      x = [] 
     190      _WCLPS = [] 
     191      _PCLPS = [] 
     192      for m in range(len(indexW_ids)): 
     193          for idx in indexW_ids[m]: 
     194              _WCLPS.append(WCLPS(npts, generations=ngcol, axisW = m, indexW = idx)) 
     195               
     196      for m in range(len(indexP_ids)): 
     197          for pair in indexP_ids[m]: 
     198              _PCLPS.append( PCLPS(npts, tol_position=tolP, generations=ngcol, \ 
     199                                    axisP = m, index1 = pair[0], index2  = pair[1]))  
     200      x.extend(_WCLPS) 
     201      x.extend(_PCLPS) 
     202      x.append(_COG) 
     203      trm = Or(*x) 
     204       
     205      if debug: 
     206        print 'list_vanished_weights = ' 
     207        print  list_vanished_weights 
     208        print 'list_collapsed_pairs = ' 
     209        print  list_collapsed_pairs                          
     210   
    60211  solved = solver.bestSolution 
    61212 #print "solved: %s" % solver.Solution() 
     
    65216  from mystic.munge import write_support_file 
    66217  write_support_file(stepmon) 
    67   return solved, func_max, func_evals 
     218  ''' 
     219  if collapse_data: 
     220    f = open('collapse_data.txt','w') 
     221    f.write('COLLAPSE DATA = %s' %collapse_data) 
     222    f.close()''' 
     223  return solved, func_max, tot_evaluations, collapse_data 
    68224 
    69225 
     
    90246 
    91247  # generate primary constraints function 
    92   def constraints(rv): 
     248  def constraints(rv, cnstr_x=None): 
    93249    c = product_measure() 
    94250    c.load(rv, npts) 
     
    103259    if not (E <= float(target[0] + error[0])) \ 
    104260    or not (float(target[0] - error[0]) <= E): 
    105       c.set_expect((target[0],error[0]), model, (x_lb,x_ub)) 
     261      c.set_expect((target[0],error[0]), model, (x_lb,x_ub), cnstr_x) 
    106262    ###################### end function-specific ###################### 
    107263    # extract weights and positions 
    108264    return c.flatten() 
    109265 
    110   # generate maximizing function 
     266  # generate maximizing functio 
    111267  def cost(rv): 
    112268    c = product_measure() 
     
    119275 
    120276  # maximize 
    121   solved, func_max, func_evals = optimize(cost,(lb,ub,npts),constraints) 
     277  solved, func_max, tot_evaluations, collapse_data = optimize(cost,(lb,ub,npts),constraints) 
    122278 
    123279  if MINMAX == 1: 
     
    125281  else: 
    126282    print "func_maximum: %s" % func_max  # sup 
    127   print "func_evals: %s" % func_evals 
    128  
    129   return solved, func_max 
     283  print "tot_evaluations: %s" % tot_evaluations 
     284 
     285  return solved, func_max, collapse_data 
    130286 
    131287 
     
    134290####################################################################### 
    135291if __name__ == '__main__': 
    136   function_name = model.__name__ 
     292  from time import clock 
     293  start = clock() 
     294   
     295  from mystic.tools import wrap_function 
     296  from mystic.monitors import Null 
     297  fmon = Null() 
     298  model_evaluations, model = wrap_function(model0, [], fmon)   
     299  function_name = model0.__name__ 
    137300 
    138301  H_mean = 6.5    #NOTE: SET THE 'mean' HERE! 
    139302  H_range = 1.0   #NOTE: SET THE 'range' HERE! 
    140   nx = 2  #NOTE: SET THE NUMBER OF 'h' POINTS HERE! 
    141   ny = 2  #NOTE: SET THE NUMBER OF 'a' POINTS HERE! 
    142   nz = 2  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
     303  nx = 5  #NOTE: SET THE NUMBER OF 'h' POINTS HERE! 
     304  ny = 5  #NOTE: SET THE NUMBER OF 'a' POINTS HERE! 
     305  nz = 5  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    143306  target = (H_mean,) 
    144307  error = (H_range,) 
     
    158321  print "...SETTINGS..." 
    159322  print "npop = %s" % npop 
     323  print "ngen = %s" %ngen 
     324  print "ngcol = %s" %ngcol 
    160325  print "maxiter = %s" % maxiter 
    161326  print "maxfun = %s" % maxfun 
     
    192357  npts = (nx,ny,nz) 
    193358  bounds = (lower_bounds,upper_bounds) 
    194   solved, diameter = maximize(pars,npts,bounds) 
     359  solved, diameter, collapse_data = maximize(pars,npts,bounds) 
    195360 
    196361  from numpy import array 
     
    203368 
    204369  print "expect: %s" % str( c.get_expect(model) ) 
    205  
     370   
     371  elapsed = (clock() - start) 
     372   
     373  print "########################" 
     374  print 'RUNTIME = %s' %elapsed 
     375  print 'COLLAPSE DATA= %s' %collapse_data 
     376  print 'MODEL Evaluations = %s' %model_evaluations 
    206377# EOF 
  • branches/collapse/collapse_termination.py

    r541 r560  
    1 debug = True #False 
    2  
    3 def CLPS_Weight(npts, tolW=0.05, tolG=10e-6, generations=70): 
     1debug = False #True 
     2 
     3def CLPS_Weight(npts, tol_weight=0.05, generations=200, axisW = '', \ 
     4                indexW = ''): 
    45    """ Termination condition for vanishing weights. 
    56    if info = False: 
     
    78    if info = True:  
    89       prints out a message saying for which points the weight vanished 
     10        
     11    In this code, weight collapse is the event where the weight of a support  
     12    point (its maximum accros last given number of generations) has been  
     13    less than or equal to the tollerance level, tol_weight.  
     14     
     15    If axisW and index are not declared, the function returns whenever weight  
     16    has vanished for any point. If declared (axisW and index are tuple or  
     17    integer objects) then the function returns only if collapse had occured f 
     18    or points with given indices. Note that neither or both axisW and indexW  
     19    must be declared, otherwise an error is printed. 
    920     
    1021    npts = (nx, ny, ...) is the tuple of the support point numbers for each 
    1122    input.""" 
     23    # from mystic.tools import isNull 
     24    from mystic.monitors import Null, Monitor, VerboseMonitor 
     25    monitor = Monitor() 
     26    if debug: monitor = VerboseMonitor(2, 2) 
     27       
     28 
     29    def _CLPS_Weight(inst, info=False): 
     30        ###FIXME: this should probably be in Solver ####### 
     31        # if isNull(inst._stepmon): inst._stepmon = monitor 
     32        ################################################### 
     33         
     34        ######################## TO BE DELETED################################## 
     35        #print "HERE I AM and the last 10 best solutions are: " 
     36        #print inst._stepmon.x[-10:] 
     37        ########################################################################           
     38         
     39        ############## PARAMETERS (length, and number of emasures) ############# 
     40        hist_lg = len(inst._stepmon.x) 
     41        if hist_lg <= generations: return False 
     42         
     43        from mystic.math.measures import _nested_split    
     44        w_bestSolution, x_bestSolution = \ 
     45                _nested_split(inst.bestSolution, npts) 
     46        ######################################################################## 
     47        collapseWeight_msg = '' 
     48                 
     49        if isinstance(axisW, int) and isinstance(indexW, int): 
     50            max_weight = 0 
     51            for n in range(-generations, 0): 
     52                weights, positions = _nested_split(inst._stepmon.x[n], npts) 
     53                if max_weight < weights[axisW][indexW]: 
     54                    max_weight = weights[axisW][indexW] 
     55            if max_weight <= tol_weight: 
     56                collapseWeight_msg += 'measure: %s & index: %s, ' \ 
     57                %(axisW,indexW)                 
     58             
     59        elif isinstance(axisW, tuple) and isinstance(indexW, int): 
     60            max_weight = [0]*len(axisW) 
     61             
     62            for n in range(-generations, 0): 
     63                weights, positions = _nested_split(inst._stepmon.x[n], npts) 
     64                for i in range(len(max_weight)): 
     65                    if max_weight[i] < weights[axisW[i]][indexW]: 
     66                        max_weight[i] = weights[axisW[i]][indexW] 
     67            for k in range(len(max_weight)): 
     68                if max_weight[k] <= tol_weight: 
     69                    collapseWeight_msg += 'measure: %s & index: %s, ' \ 
     70                    %(axisW[k],indexW)    
     71                     
     72        elif isinstance(axisW, int) and isinstance(indexW, tuple): 
     73            max_weight = [0]*len(indexW) 
     74                         
     75            for n in range(-generations, 0): 
     76                weights, positions = _nested_split(inst._stepmon.x[n], npts) 
     77                for i in range(len(max_weight)): 
     78                    if max_weight[i] < weights[axisW][indexW[i]]: 
     79                        max_weight[i] = weights[axisW][indexW[i]] 
     80            for k in range(len(max_weight)): 
     81                if max_weight[k] <= tol_weight: 
     82                    collapseWeight_msg += 'measure: %s & index: %s, ' \ 
     83                    %(axisW,indexW[k]) 
     84                     
     85        elif (isinstance(axisW, tuple) and isinstance(indexW, tuple)) or \ 
     86            (axisW == '' and indexW == ''): #Note that can't change to not axisW 
     87                                            # and not indexW since those can be zero indices 
     88            if (axisW == '' and indexW == ''):  
     89                max_weight = [[0] * nx for nx in npts] 
     90            else: 
     91                max_weight = [[0]*len(axisW)]*len(indexW) 
     92                 
     93            for n in range(-generations, 0): 
     94                weights, positions = _nested_split(inst._stepmon.x[n], npts) 
     95                 
     96                for i in range(len(max_weight)):    
     97                    for j in range(len(max_weight[i])): 
     98                        if (axisW =='' and indexW == ''): 
     99                            if  max_weight[i][j] < weights[i][j]: 
     100                                max_weight[i][j] = weights[i][j] 
     101                        else: 
     102                            if  max_weight[i][j] < weights[axisW[i]][index[j]]: 
     103                                max_weight[i][j] = weights[axisW[i]][index[j]] 
     104                                             
     105            for k in range(len(max_weight)): 
     106                for l in range(len(max_weight[k])): 
     107                    if max_weight[k][l] <= tol_weight: 
     108                        if (axisW == '' and indexW == ''): 
     109                            collapseWeight_msg += 'measure: %s & index: %s, ' \ 
     110                            %(k,l) 
     111                        else:                            
     112                            collapseWeight_msg += 'measure: %s & index: %s, ' \ 
     113                            %(axisW[k],indexW[l]) 
     114                             
     115        else:  
     116            print 'Error, wrong arguments of CLPS_Weight' 
     117            return False 
     118                                                
     119        if info: 
     120            if collapseWeight_msg == '': msg = '' 
     121            else: 
     122                msg0 = 'CLPS_Weight with: ' 
     123                msg = msg0 + collapseWeight_msg 
     124                msg2 = msg0 + '\n npts = %s, tol_weight = %s, generations = %s: \n' \ 
     125                % (npts, tol_weight, generations)  
     126                msg2 += 'vanishing weights of support points: ' 
     127                msg2 += collapseWeight_msg 
     128                if debug: print msg2 
     129            return msg 
     130 
     131        if collapseWeight_msg == '': return False 
     132        else: return True 
     133     
     134    return _CLPS_Weight 
     135             
     136     
     137                         
     138def CLPS_Position(npts, tol_position, generations=200, axisP = '', \ 
     139                  index1 = '', index2 = ''): 
     140     
     141    """ Termination condition for collapsing positions of support points. 
     142    if info = False: 
     143       returns true if the positions collapse for any 2 support  
     144       points in a measure 
     145    if info = True:  
     146       prints out a message saying for which point-pairs the location collapsed. 
     147     
     148    In this code, collapse is the event where the distance between 2 support  
     149    points (its maximum accros last given number of generations) has been  
     150    less than the tollerance level, tol_position.  
     151     
     152    If axisP and index1 and index2 are not declared, the function returns  
     153    whenever any 2 points collapse in position. If declared (axisP is a tuple  
     154    or and integer, index1 and index2 are integers) then the function returns  
     155    only if collapse had occured for points with given indices. 
     156     
     157    npts = (nx, ny, ...) is the tuple of the support point numbers for each 
     158    input.""" 
     159     
    12160    from mystic.tools import isNull 
    13161    from mystic.monitors import Null, Monitor, VerboseMonitor 
    14162    monitor = Monitor() 
    15     if debug: monitor = VerboseMonitor(2) 
    16  
    17     def _CLPS_Weight(inst, info=False): 
     163    if debug: monitor = VerboseMonitor(2, 2) 
     164     
     165    # CHECK THAT index1, and index2 are integers if declared, and either none or  
     166    # both are empty 
     167     
     168    def _CLPS_Position(inst, info=False): 
    18169        ###FIXME: this should probably be in Solver ####### 
    19170        if isNull(inst._stepmon): inst._stepmon = monitor 
    20171        ################################################### 
    21172 
    22         lg = len(inst._stepmon) 
    23         if lg <= generations: return False 
    24          
    25         from numpy import array, ones, transpose 
    26         from mystic.math.dirac_measure import product_measure 
    27          
    28         Last_bestSolution = product_measure() 
    29         Last_bestSolution.load(inst._stepmon.x[-1], npts) 
    30              
    31         Gen_bestSolution = product_measure() 
    32         Gen_bestSolution.load(inst._stepmon.x[-generations], npts) 
    33          
    34         dim = len(Last_bestSolution) # Number of measures 
    35         collapseWeight  = ''       
    36  
     173        ############## PARAMETERS (length, and number of emasures) ############# 
     174        hist_lg = len(inst._stepmon) 
     175        if hist_lg <= generations: return False 
     176         
     177        from mystic.math.measures import _nested_split   
     178        w_bestSolution, x_bestSolution = \ 
     179                _nested_split(inst.bestSolution, npts) 
     180         
     181        if axisP == '': 
     182            dim = len(w_bestSolution) # Number of measures 
     183            axis = range(dim) 
     184        else: 
     185            axis = axisP             
     186            if isinstance( axisP, int ): dim  = 1 
     187            else: dim = len(list(axisP)) 
     188         
     189        ######################################################################## 
     190        collapsePairs_msg  = ''       
     191         
     192        from numpy import array, zeros, ones, transpose             
    37193        for i in range(dim): 
    38             Last_collapseWeight = array(Last_bestSolution[i].weights) <= tolW 
    39             Gen_collapseWeight = array(Gen_bestSolution[i].weights) <= tolW 
    40             for k in range(len(Last_collapseWeight)): 
    41                 if all([Last_collapseWeight[k], Gen_collapseWeight[k]]): 
    42                        collapseWeight += ('measure %s: index: %s; ' %(i,k)) 
     194            if isinstance(axis, int): ax_idx = axis 
     195            else: ax_idx = axis[i] 
     196                         
     197            if isinstance(index1, int): 
     198                if isinstance(index2, int): 
     199                    nx = 1 
     200                    max_measure_diffLoc = 0 
     201                else: 
     202                    print "Error argument declaration index2 is \ 
     203                    empty while index1 is not" 
     204                    return False   #FIXME 
     205            else: 
     206                if isinstance(index2, int): 
     207                    print "Error argument declaration index1 is \ 
     208                    empty while index2 is not" 
     209                    return False   #FIXME 
     210                else: 
     211                    nx = npts[ax_idx] 
     212                    max_measure_diffLoc = zeros((nx, nx))      
     213                             
     214            for n in range(-generations, 0): 
     215                weights, positions = _nested_split(inst._stepmon.x[n], npts) 
     216                if nx == 1:                        
     217                    diffLoc = abs(positions[ax_idx][index1] - \ 
     218                                  positions[ax_idx][index2]) 
     219                    if max_measure_diffLoc < diffLoc: 
     220                                            max_measure_diffLoc= diffLoc                       
     221                else: 
     222                    diffLoc = abs(transpose(positions[ax_idx]*ones((nx, nx))) - \ 
     223                                      positions[ax_idx]*ones((nx, nx)))                                       
     224                    for k in range(nx): 
     225                        for l in range(nx): 
     226                            if max_measure_diffLoc[(k, l)] < diffLoc[(k, l)]: 
     227                                max_measure_diffLoc[(k, l)] = diffLoc[(k, l)] 
     228 
     229            if nx ==1 and max_measure_diffLoc <= tol_position[ax_idx]:  
     230                    collapsePairs_msg += ('measure: %s & indices: (%s,%s),'\ 
     231                                          %(ax_idx,index1,index2)) 
     232                
     233            else:             
     234                for k in range(nx): 
     235                    for l in range(k,nx): 
     236                        if k != l and max_measure_diffLoc[(k,l)] <= \ 
     237                           tol_position[ax_idx]: 
     238                            collapsePairs_msg += ('measure: %s & indices: (%s,%s),' \ 
     239                            %(ax_idx,k,l)) 
     240     
    43241        if info: 
    44             if collapseWeight == '': msg = '' 
    45             else: 
    46                 msg = 'CLPS_Weight:\n' 
    47                 msg += 'npts = %s, tolW = %s, tolG = %s, generations = %s;\n' \ 
    48                 % (npts, tolW, tolG, generations)  
    49                 msg += 'vanishing weights of support points: ' 
    50                 msg += collapseWeight 
    51                 #if debug: print msg 
     242            if collapsePairs_msg == '': msg = '' 
     243            else: 
     244                msg0 = 'CLPS_Position with: ' 
     245                msg = msg0 + collapsePairs_msg 
     246                msg2 = msg0 + '\n npts = %s, tol_position = %s, generations = %s:\n' \ 
     247                % (npts, tol_position, generations)  
     248                msg2 += 'collapse of support points: ' 
     249                msg2 += collapsePairs_msg 
     250                if debug: print msg2 
    52251            return msg 
    53  
    54         if collapseWeight == '': return False 
     252         
     253        if collapsePairs_msg == '': return False 
    55254        else: return True 
    56255     
    57     return _CLPS_Weight 
    58              
    59      
    60              
    61              
    62 def CLPS_Location(npts, tolL, tolG=10e-6, generations=30): 
    63      
    64     """ Termination condition for collapsing locations of support points. 
    65     if info = False: 
    66        returns true if the locations collapse for any 2 support  
    67        points in a measure 
    68     if info = True:  
    69        prints out a message saying for which point-pairs the location collapsed. 
    70      
    71     npts = (nx, ny, ...) is the tuple of the support point numbers for each 
    72     input.""" 
    73     from mystic.tools import isNull 
    74     from mystic.monitors import Null, Monitor, VerboseMonitor 
    75     monitor = Monitor() 
    76     if debug: monitor = VerboseMonitor(2) 
    77      
    78     def _CLPS_Location(inst, info=False): 
    79         ###FIXME: this should probably be in Solver ####### 
    80         if isNull(inst._stepmon): inst._stepmon = monitor 
    81         ################################################### 
    82  
    83         lg = len(inst._stepmon) 
    84         if lg <= generations: return False 
    85          
    86         from numpy import array, ones, transpose 
    87         from mystic.math.dirac_measure import product_measure 
    88          
    89         Last_bestSolution = product_measure() 
    90         Last_bestSolution.load(inst._stepmon.x[-1], npts) 
    91          
    92         Gen_bestSolution = product_measure() 
    93         Gen_bestSolution.load(inst._stepmon.x[-generations], npts) 
    94          
    95         dim = len(Last_bestSolution) # Number of measures 
    96         collapsePairs  = ''       
    97    
    98         for i in range(dim): 
    99             x = array(Last_bestSolution[i].coords)    
    100             Last_DiffLoc = transpose(x*ones((len(x), len(x)))) - \ 
    101                     x*ones((len(x), len(x))) 
    102                      
    103             y = array(Gen_bestSolution[i].coords)      
    104                      
    105             Gen_DiffLoc = transpose(y*ones((len(y), len(y)))) - \ 
    106                             y*ones((len(y), len(y)))  
    107                                        
    108             collapseLoc1 = (abs(Last_DiffLoc) <= tolL[i])  
    109             collapseLoc2 = (abs(Gen_DiffLoc) <= tolL[i]) 
    110              
    111             for k in range(len(x)): 
    112                 for l in range(len(x)): 
    113                     if k != l and all([collapseLoc1[(k,l)], \ 
    114                                        collapseLoc2[(k, l)]]): 
    115                         collapsePairs += ('measure %s: indices (%s,%s); ' \ 
    116                                           %(i,k, l)) 
    117          
    118         if info: 
    119             if collapsePairs == '': msg = '' 
    120             else: 
    121                 msg = 'CLPS_Location:\n' 
    122                 msg += 'npts = %s, tolL = %s, tolG = %s, generations = %s;\n' \ 
    123                 % (npts, tolL, tolG, generations)  
    124                 msg += 'collapse of support points: ' 
    125                 msg += collapsePairs 
    126                 #if debug: print msg 
    127             return msg 
    128  
    129         if collapsePairs == '': return False 
    130         else: return True 
    131      
    132     return _CLPS_Location 
    133          
    134  
     256    return _CLPS_Position 
     257         
     258 
Note: See TracChangeset for help on using the changeset viewer.