Changeset 560 for branches/collapse/TEST_OUQ_surrogate_diam.py

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

File:

• branches/collapse/TEST_OUQ_surrogate_diam.py (modified) (14 diffs)

branches/collapse/TEST_OUQ_surrogate_diam.py

@@ -1 +1 @@
 #!/usr/bin/env python
 
-debug = False
+debug = True #False
 MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1
 #######################################################################
@@ -8 +8 @@
 #######################################################################
 npop = 40
-maxiter = 1000
+maxiter = 6000
 maxfun = 1e+6
-convergence_tol = 1e-6; ngen = 40
 crossover = 0.9
 percent_change = 0.9
+convergence_tol = 1e-6
+tolW = 0.05; tolP_range_fr = 0.05
+ngen = 200; ngcol = 200
 
 
@@ -18 +20 @@
 # the model function
 #######################################################################
-from surrogate import marc_surr as model
-
+from surrogate import marc_surr as model0
 
 #######################################################################
@@ -26 +27 @@
 def optimize(cost,_bounds,_constraints):
   from mystic.solvers import DifferentialEvolutionSolver2
- #from mystic.termination import ChangeOverGeneration as COG
-  from mystic.termination import Or
-  from collapse_termination import CLPS_Weight as WCLPS
-  from collapse_termination import CLPS_Location as LCLPS
   from mystic.strategy import Best1Exp
   from mystic.monitors import VerboseMonitor, Monitor
   from mystic.tools import random_seed
 
- #random_seed(123)
-
-  stepmon = VerboseMonitor(2)
+  random_seed(122)
+
+  stepmon = VerboseMonitor(10,10)
  #stepmon = Monitor()
   evalmon = Monitor()
@@ -50 +47 @@
   solver.SetGenerationMonitor(stepmon)
   solver.SetConstraints(_constraints)
-
+ 
+  ##############################################################################
+  
+                 ## THE FOLLOWING TERMINATION AND CONSTRAINTS ##
+                 ## CAN BE MOVED TO A NEW MODULE AND IMPORTED ##
+  
+  ######################### TERMINATION CONDITIONS #############################
+  from mystic.termination import Or  
+  from mystic.termination import ChangeOverGeneration as COG
+  from collapse_termination_v2 import CLPS_Weight as WCLPS
+  from collapse_termination_v2 import CLPS_Position as PCLPS  
+  from mystic.math.measures import _nested_split  
+  
   tol = convergence_tol
-  _WCLPS = WCLPS(npts, tolG=tol, generations=ngen)
-  _LCLPS = LCLPS(npts, tolL=[1,1,0.03], tolG=tol, generations=ngen) #FIXME
-  CLPS = Or( _WCLPS, _LCLPS )
-  solver.Solve(cost,termination=CLPS,strategy=Best1Exp, \
-               CrossProbability=crossover,ScalingFactor=percent_change)
-
+  _COG = COG(tol,ngen)
+  
+  w_lb, x_lb = _nested_split(lb, npts)
+  w_ub, x_ub = _nested_split(ub, npts)
+  range_x = []
+  for i in range(len(x_lb)):
+      range_x.append(abs(x_lb[i][0] - x_ub[i][0]))  
+  tolP = [ele*tolP_range_fr for ele in range_x] 
+  if debug: print 'tolP = %s' %tolP
+  
+  indexW_ids = [range(nx) for nx in npts]  # list of indices for CLPS_Weight termination check
+  indexP_ids = []  # list of indices for CLPS_Position termination check  
+  
+  x = []
+  _WCLPS = []
+  _PCLPS = []
+  
+  for m in range(len(indexW_ids)):
+      for idx in indexW_ids[m]:
+          _WCLPS.append(WCLPS(npts, tol_weight = tolW, generations=ngcol, axisW = m, indexW = idx))
+          
+  for m in range(len(indexP_ids)):
+      for pair in indexP_ids[m]:
+          _PCLPS.append( PCLPS(npts, tol_position=tolP, generations=ngcol, \
+                                axisP = m, index1 = pair[0], index2  = pair[1])) 
+  x.extend(_WCLPS)
+  x.extend(_PCLPS)
+  x.append(_COG)
+  trm = Or(*x)
+  
+  ######################## NEW CONSTRAINTS AFTER COLLAPSE ######################
+  
+  list_vanished_weights = [] # list of [measure, index] lists
+  list_collapsed_pairs = [] # list of [measure, pair] lists, where pair = [index1, index2]  
+  
+  for i in range(len(npts)):
+    c = []
+    for idx in range(npts[i]):
+        for jdx in range(idx,npts[i]):
+            if idx !=jdx:
+                c.append([idx, jdx])
+    indexP_ids.append(c)  
+  
+  def _set_weight_to_zero(rv, npts, axisW, indexW):
+    from mystic.math.dirac_measure import product_measure
+    c = product_measure()
+    c.load(rv, npts)
+    c[axisW][indexW].weight = 0
+    return c.flatten()
+  
+  def _set_same_positions(rv, npts, axisP, index1, index2):
+    from mystic.math.dirac_measure import product_measure
+    c = product_measure()
+    c.load(rv, npts)
+    c[axisP][index1].weight += c[axisP][index2].weight
+    c[axisP][index2].weight = 0
+    c[axisP][index2].position = c[axisP][index1].position
+    return c.flatten() 
+  
+  def new_constraints(rv, list_collapsed_pairs, list_vanished_weights):
+    for ele in list_vanished_weights:
+      rv = _set_weight_to_zero(rv, npts, axisW = ele[0], indexW = ele[1])
+    for mem in list_collapsed_pairs:
+      rv = _set_same_positions(rv, npts, axisP = mem[0], \
+              index1 = mem[1][0], index2 = mem[1][1])            
+    return _constraints(rv)   
+  
+  ##############################################################################
+  
+  
+  ############################ RUN THE SOLVER ##################################
+  collapse_data = {}
+  tot_evaluations = 0 # total number of function eveluations
+  while _COG(solver) != True:         
+    solver.Solve(cost,termination=trm,strategy=Best1Exp, \
+               CrossProbability=crossover,ScalingFactor=percent_change) 
+    tot_evaluations += solver.evaluations    
+    collapse_data[stepmon._step] = trm(solver, True)
+    
+    if debug:
+      print 'last_itr_evaluations = %s' %solver.evaluations
+      print "Generation  %d " %stepmon._step 
+    
+    int_trm = trm(solver, info = 'self')  
+    for ele in range(len(int_trm)):
+      cond = int_trm[ele]
+      string = cond(solver, info = True)
+      string1 = string.replace(',', ' ')
+      string2 = string1.replace('&', ' ')
+      string3 = string2.replace(':', ' ')
+      lst = string3.split()
+      
+      if lst[0] == 'CLPS_Weight':
+          for i in range(len(lst)):
+              if lst[i] == 'measure':
+                  axisW = int(lst[i+1])
+                  if lst[i+2] == 'index': indexW = int(lst[i+3])
+                  
+          if debug:
+            print 'axisW = %s' %axisW
+            print 'indexW = %s' %indexW
+          
+          list_vanished_weights.append([axisW, indexW]) 
+          indexW_ids[axisW].remove(indexW)
+                 
+      if lst[0] == 'CLPS_Position':
+          for i in range(len(lst)):  
+              if lst[i] == 'measure':        
+                  axisP = int(lst[i+1])
+                  if lst[i+2] == 'indices': 
+                      index1 = int(lst[i+3].lstrip('('))
+                      index2 = int(lst[i+4].rstrip(')'))
+          if debug:
+            print "axisP = %s" %axisP
+            print "index1 = %s" %index1
+            print "index2 = %s" %index2
+          
+          list_collapsed_pairs.append([axisP, [index1, index2]])
+          indexP_ids[axisP].remove([index1,index2])          
+      
+      ##### AFTER COLLAPSE OCCURRED CHANGE THE TERMINATION AND CONSTRAINTS #####
+                    
+      def _new_constraints(rv):
+        return new_constraints(rv, list_collapsed_pairs, list_vanished_weights)      
+      solver.SetConstraints(_new_constraints)
+      
+      x = []
+      _WCLPS = []
+      _PCLPS = []
+      for m in range(len(indexW_ids)):
+          for idx in indexW_ids[m]:
+              _WCLPS.append(WCLPS(npts, generations=ngcol, axisW = m, indexW = idx))
+              
+      for m in range(len(indexP_ids)):
+          for pair in indexP_ids[m]:
+              _PCLPS.append( PCLPS(npts, tol_position=tolP, generations=ngcol, \
+                                    axisP = m, index1 = pair[0], index2  = pair[1])) 
+      x.extend(_WCLPS)
+      x.extend(_PCLPS)
+      x.append(_COG)
+      trm = Or(*x)
+      
+      if debug:
+        print 'list_vanished_weights = '
+        print  list_vanished_weights
+        print 'list_collapsed_pairs = '
+        print  list_collapsed_pairs                         
+  
   solved = solver.bestSolution
  #print "solved: %s" % solver.Solution()
@@ -65 +216 @@
   from mystic.munge import write_support_file
   write_support_file(stepmon)
-  return solved, func_max, func_evals
+  '''
+  if collapse_data:
+    f = open('collapse_data.txt','w')
+    f.write('COLLAPSE DATA = %s' %collapse_data)
+    f.close()'''
+  return solved, func_max, tot_evaluations, collapse_data
 
 
@@ -90 +246 @@
 
   # generate primary constraints function
-  def constraints(rv):
+  def constraints(rv, cnstr_x=None):
     c = product_measure()
     c.load(rv, npts)
@@ -103 +259 @@
     if not (E <= float(target[0] + error[0])) \
     or not (float(target[0] - error[0]) <= E):
-      c.set_expect((target[0],error[0]), model, (x_lb,x_ub))
+      c.set_expect((target[0],error[0]), model, (x_lb,x_ub), cnstr_x)
     ###################### end function-specific ######################
     # extract weights and positions
     return c.flatten()
 
-  # generate maximizing function
+  # generate maximizing functio
   def cost(rv):
     c = product_measure()
@@ -119 +275 @@
 
   # maximize
-  solved, func_max, func_evals = optimize(cost,(lb,ub,npts),constraints)
+  solved, func_max, tot_evaluations, collapse_data = optimize(cost,(lb,ub,npts),constraints)
 
   if MINMAX == 1:
@@ -125 +281 @@
   else:
     print "func_maximum: %s" % func_max  # sup
-  print "func_evals: %s" % func_evals
-
-  return solved, func_max
+  print "tot_evaluations: %s" % tot_evaluations
+
+  return solved, func_max, collapse_data
 
 
@@ -134 +290 @@
 #######################################################################
 if __name__ == '__main__':
-  function_name = model.__name__
+  from time import clock
+  start = clock()
+  
+  from mystic.tools import wrap_function
+  from mystic.monitors import Null
+  fmon = Null()
+  model_evaluations, model = wrap_function(model0, [], fmon)  
+  function_name = model0.__name__
 
   H_mean = 6.5    #NOTE: SET THE 'mean' HERE!
   H_range = 1.0   #NOTE: SET THE 'range' HERE!
-  nx = 2  #NOTE: SET THE NUMBER OF 'h' POINTS HERE!
-  ny = 2  #NOTE: SET THE NUMBER OF 'a' POINTS HERE!
-  nz = 2  #NOTE: SET THE NUMBER OF 'v' POINTS HERE!
+  nx = 5  #NOTE: SET THE NUMBER OF 'h' POINTS HERE!
+  ny = 5  #NOTE: SET THE NUMBER OF 'a' POINTS HERE!
+  nz = 5  #NOTE: SET THE NUMBER OF 'v' POINTS HERE!
   target = (H_mean,)
   error = (H_range,)
@@ -158 +321 @@
   print "...SETTINGS..."
   print "npop = %s" % npop
+  print "ngen = %s" %ngen
+  print "ngcol = %s" %ngcol
   print "maxiter = %s" % maxiter
   print "maxfun = %s" % maxfun
@@ -192 +357 @@
   npts = (nx,ny,nz)
   bounds = (lower_bounds,upper_bounds)
-  solved, diameter = maximize(pars,npts,bounds)
+  solved, diameter, collapse_data = maximize(pars,npts,bounds)
 
   from numpy import array
@@ -203 +368 @@
 
   print "expect: %s" % str( c.get_expect(model) )
-
+  
+  elapsed = (clock() - start)
+  
+  print "########################"
+  print 'RUNTIME = %s' %elapsed
+  print 'COLLAPSE DATA= %s' %collapse_data
+  print 'MODEL Evaluations = %s' %model_evaluations
 # EOF