rmd_node__define.pro

该程序是基于IDL平台下开法的应用随机搜索树进行搜索的优化算法实例！

                           _Extra = extra

if n_elements(extra) ne 0 then self.pextra = ptr_new(extra,/no_copy) else $
   self.pextra = ptr_new(/allocate_heap)
if n_params() eq 0 then return,0
if n_elements(decay_width) ne 0 then self.decay_width = decay_width else $
   self.decay_width = 5.0
if n_elements(parent_node) ne 0 then self.parent_node = parent_node
if n_elements(func) ne 0 then self.func = func else $
   self.func = 'rmd_node_func'
self.id = id
self.prange_ptr = ptr_new(/allocate_heap)
if n_elements(level) ne 0 then self.level = level
self.nsubdiv = 2
if n_elements(prange) ne 0 then *self.prange_ptr = prange
; The contents of PRANGE are expected to have the
; dimensions of 2 x NP where NP is the number of
; parameters.
;
; Determine the node value by evaluating the function
ret = self->evaluate_function()
self.branch_value = 1L
self.increment_value = 1L
if n_elements(increment_value) ne 0 then $
   self.increment_value = increment_value
return,1
end
; *********************************** ;
pro rmd_node__define
def = {  rmd_node,                     $
         inherits idl_container,       $
         id:'',                        $  ;  string identifier for node
         level:0L,                     $  ;  level (depth) of node
         nsubdiv:0L,                   $  ;  number of subdivisions per node (2 for now)
         visits:0L,                    $  ;  number of visits to this node
         decay_width:0.,               $  ;  std deviation for probability decay
         parent_node:obj_new(),        $  ;  obj reference for parent of this node
         node_value:0D,                $  ;  value (func evaluation) for this node
         branch_value:0L,              $  ;  branch value for branch *leading* to this node
         increment_value:0L,           $  ;  amount by which to increment the branch upon
                                          ;  creating a new *best* node
         pextra:ptr_new(),             $  ;  keyword inheritance mechanism used for passing
                                          ;  information along to the function to be minimized
         func:'',                      $  ;  string name for the function to be minimized
         prange_ptr:ptr_new()          }  ;  Parameter ranges
end
; *********************************** ;
; *********************************** ;
; Begin example usage of RMD_NODE__DEFINE
; *********************************** ;
function test_node_func,p,_EXTRA = extra
y = (p-1.25)^2
return,y
end
; *********************************** ;
function test_2d_node_func,p,_EXTRA = extra
x = p[0] & y = p[1]
z = 3.*(1.-x)^2 * exp(-x^2 - (y+1.)^2) - $
    10.*(x/5 - x^3 - y^5)*exp(-x^2-y^2) - $
    (1./3)*exp(-(x+1)^2 - y^2)
return,z
end
; *********************************** ;
function node_2d_func,x1,y1
ux = 1+bytarr(n_elements(x1))
uy = 1+bytarr(n_elements(y1))
x = x1#uy & y = ux#y1
z = 3.*(1.-x)^2 * exp(-x^2 - (y+1.)^2) - $
    10.*(x/5 - x^3 - y^5)*exp(-x^2-y^2) - $
    (1./3)*exp(-(x+1)^2 - y^2)
return,z
end
; *********************************** ;
pro test_node_min
prange = [-3.0,10.0]
func = 'test_node_func'
itmax = 15  ; max number of search iterations
; Display the function
xlo = prange[0] & xhi = prange[1] & nx = 500 & dx = (xhi-xlo)/(nx-1.)
x = xlo+dx*findgen(nx)
y = call_function(func,x)
plot,x,y,psym = 0

; Create the first node object
ostart = obj_new('rmd_node',  '0-0',                     $
                              prange = prange,           $
                              increment_value = 10,      $
                              decay_width = 3.0,         $
                              func = func,               $
                              level = 0                  )

; Set up the main iteration loop
iter = 0L
obest = obj_new()
best_eval = 1e12
ocurrent = ostart
best_iter = 0L
num_nodes = 0L

while iter le itmax do begin
   ; Is the current node a terminal node?
   if ocurrent->terminated() then begin   ; Yes
      ocurrent->get_property,parms = parms
      ; Create new nodes.  Note that the node values
      ; are actually calculated when the new node objects
      ; are instantiated using the add_nodes method.
      next_index = num_nodes+1

      ret = ocurrent->add_nodes(next_index)
      num_nodes = num_nodes+ocurrent->count()
      ; For each new node, get the best node value.
      nchildren = ocurrent->count()
      ochild = ocurrent->get(/all)
      level_values = fltarr(nchildren)
      for i = 0,nchildren-1 do begin
         ochild[i]->get_property,node_value = node_value,level = level,parms = parms
         ocurrent->get_property,parms = parent_parms
         level_values[i] = node_value
         ; Traverse up the tree and decrement each branch according
         ; to how many visits each node has had.
         ret = ochild[i]->go_to_the_top(/decrease)
      endfor
      min_level_value = min(level_values,imin)
      ; How does this compare with the current best?
      if min_level_value lt best_eval then begin
         ; Replace the current best node with this one
         obest = ochild[imin]
         best_eval = min_level_value
         best_iter = iter
         ; Traverse back up the tree and increment each branch as we go.
         ret = obest->go_to_the_top(/increase)
      endif
      iter++
      ocurrent = ostart

   endif else begin  ; No, not a terminal node...keep going

      ; Calculate the relative probabilities of all branches to the
      ; available nodes.
      nchildren = ocurrent->count()
      ochild = ocurrent->get(/all)
      branch_prob = fltarr(nchildren)
      for i = 0,nchildren-1 do begin
         ochild[i]->get_property,node_value = node_value, $
                                 branch_value = branch_value
         branch_prob[i] = branch_value
      endfor
      branch_prob = branch_prob/total(branch_prob)
      regions = [0.0,total(branch_prob,/cum)]
      select_index = value_locate(regions,randomu(s,1))
      ocurrent = ochild[select_index]
      ret = ocurrent->add_visit()

   endelse

endwhile

; If a good solution has been found then plot it on the
; function that was to be minimized.
if obj_valid(obest) then begin
   obest->get_property,parms = parms,node_value = node_value,func = func
   plots,[parms[0]],[call_function(func,parms[0])],psym = 4,symsize = 2.0, $
      thick = 2.0
endif

obj_destroy,ostart
end

; *********************************** ;
pro test_node_2d_min
device,decomposed = 0
loadct,3,/silent

prange = [[-3.0,3.0],[-3.,3.]]
func = 'test_2d_node_func'
itmax = 45  ; max number of search iterations
; Display the function
xlo = prange[0,0] & xhi = prange[1,0] & nx = 500 & dx = (xhi-xlo)/(nx-1.)
x = xlo+dx*findgen(nx)
ylo = prange[0,1] & yhi = prange[1,1] & ny = 500 & dy = (yhi-ylo)/(ny-1.)
y = ylo+dy*findgen(ny)

z = call_function('node_2d_func',x,y)
plotimage,bytscl(z),imgxrange = [xlo,xhi],imgyrange = [ylo,yhi]

; Create the first node object
ostart = obj_new('rmd_node',  '0-0',                     $
                              prange = prange,           $
                              increment_value = 50,      $
                              decay_width = 3.0,         $
                              func = func,               $
                              level = 0                  )

; Set up the main iteration loop
iter = 0L
obest = obj_new()
best_eval = 1e12
ocurrent = ostart
best_iter = 0L
num_nodes = 0L

while iter le itmax do begin
   ; Is the current node a terminal node?
   if ocurrent->terminated() then begin   ; Yes
      ocurrent->get_property,parms = parms
      ; Create new nodes.  Note that the node values
      ; are actually calculated when the new node objects
      ; are instantiated using the add_nodes method.
      next_index = num_nodes+1

      ret = ocurrent->add_nodes(next_index)
      num_nodes = num_nodes+ocurrent->count()
      ; For each new node, get the best node value.
      nchildren = ocurrent->count()
      ochild = ocurrent->get(/all)
      level_values = fltarr(nchildren)
      for i = 0,nchildren-1 do begin
         ochild[i]->get_property,node_value = node_value,level = level,parms = parms
         ocurrent->get_property,parms = parent_parms
         level_values[i] = node_value
         ; Traverse up the tree and decrement each branch according
         ; to how many visits each node has had.
         ret = ochild[i]->go_to_the_top(/decrease)
      endfor
      min_level_value = min(level_values,imin)
      ; How does this compare with the current best?
      if min_level_value lt best_eval then begin
         ; Replace the current best node with this one
         obest = ochild[imin]
         best_eval = min_level_value
         best_iter = iter
         ; Traverse back up the tree and increment each branch as we go.
         ret = obest->go_to_the_top(/increase)
      endif
      iter++
      ocurrent = ostart

   endif else begin  ; No, not a terminal node...keep going

      ; Calculate the relative probabilities of all branches to the
      ; available nodes.
      nchildren = ocurrent->count()
      ochild = ocurrent->get(/all)
      branch_prob = fltarr(nchildren)
      for i = 0,nchildren-1 do begin
         ochild[i]->get_property,node_value = node_value, $
                                 branch_value = branch_value
         branch_prob[i] = branch_value
      endfor
      branch_prob = branch_prob/total(branch_prob)
      regions = [0.0,total(branch_prob,/cum)]
      select_index = value_locate(regions,randomu(s,1))
      ocurrent = ochild[select_index]
      ret = ocurrent->add_visit()

   endelse

endwhile

; If a good solution has been found then plot it on the
; function that was to be minimized.
if obj_valid(obest) then begin
   obest->get_property,parms = parms,node_value = node_value,func = func
   print,parms
   plots,[parms[0]],[parms[1]],/data,symsize = 2.0,color = 255,psym = 4
;   plots,[parms[0]],[call_function(func,parms[0])],psym = 4,symsize = 2.0, $
;      thick = 2.0
endif
zmin = min(z,imin)
ux = 1+bytarr(nx)
uy = 1+bytarr(ny)
xm = x#uy & ym = ux#y
print,xm[imin],ym[imin]
plots,!x.crange,[ym[imin],ym[imin]],psym = 0,linestyle = 2,color = 255
plots,[xm[imin],xm[imin]],!y.crange,psym = 0,linestyle = 2,color = 255
obj_destroy,ostart
end