maximize.frink

View or download maximize.frink in plain text format


/** This file contains routines to minimize or maximize a function over a
    certain range using interval arithmetic techniques. */


maximize[f, minX, maxX, minXDiff] := doOptimization[f, minX, maxX, minXDiff, 1]

minimize[f, minX, maxX, minXDiff] :=
{
   result = doOptimization[f, minX, maxX, minXDiff, -1]
   retval = new array
   for [x,y] = result
      retval.push[[x, -y]]

   return retval
}

doOptimization[f, minX, maxX, minXDiff, factor] :=
{
   // Create an orderedList that sorts by the second item in a list.
   orderingFunction = {|a,b| infimum[a@1] <=> infimum[b@1]}
   boxes          = new OrderedList[orderingFunction]
   foundSolutions = new OrderedList[orderingFunction]
   
   fullInterval = new interval[minX, maxX]
   y = factor * f[fullInterval]
   boxes.insert[[fullInterval, y]]
   lowerBound = infimum[y]
   do
   {
      count = length[boxes]
      newBoxes = new OrderedList[orderingFunction]
      for i=0 to count - 1
      {
         [box,val] = boxes@i
         if (val PGE lowerBound)   // Make sure previous (larger bound) is
         {                         // larger than the lower bound.
            y = factor * f[box]
            println["f[$box]=$y, bound=$lowerBound"]
            if (y PGE lowerBound)
            {
               width = supremum[box]-infimum[box]
               if width*2 > minXDiff
               {
                  mid = infimum[box] + width/2

                  nx = new interval[infimum[box], mid]
                  ny = factor * f[nx]
                  if infimum[ny] > lowerBound
                     lowerBound = infimum[ny]
               
                  if ny PGE lowerBound
                     newBoxes.insert[[nx,ny]]

                  nx = new interval[mid, supremum[box]]
                  ny = factor * f[nx]
                  if infimum[ny] > lowerBound
                     lowerBound = infimum[ny]
               
                  if ny PGE lowerBound
                     newBoxes.insert[[nx,ny]]
               
                  if infimum[y] > lowerBound
                     lowerBound = infimum[y]
               } else
                  foundSolutions.insert[[box,val]]
               //println[boxes]
            }
         }
      }

      // Now coalesce "close enough" solutions into unions of their intervals.
      //boxes = coalesce[newBoxes, foundSolutions, maxXIntervalWidth]
      boxes = newBoxes
   } while length[boxes] > 0  && length[boxes] < 1000

   retval = new array
   
   count = length[boxes]
   for i=0 to count-1
   {
      [box, y] = boxes@i
      if y PGE lowerBound
         retval.push[[box,y]]
   }

//   println["lower bound is $lowerBound"]
   count = length[foundSolutions]
   for i=0 to count-1
   {
      [box, y] = foundSolutions@i
      if y PGE lowerBound
         retval.push[[box,y]]
   }
//   println["Found solutions: " + foundSolutions]

   return coalesce[retval]
}

/** Coalesce "close enough" solutions into solutions. */
// TODO:  Most "close" values are 
coalesce[retval] :=
{
   last = length[retval]-1
   i = 0
   while i<=last-1
   {
      [lastX, lastY] = retval@i
      j = i+1
      while j<=last
      {
         [jx, jy] = retval@j
         int = intersection[lastX, jx]
         if int == undef   // No intersection, keep going
            j=j+1
         else
         {
            lastX = union[lastX, jx]
            lastY = union[lastY, jy]
//            println["$lastX intersects $jx, union is $lastX"]
            retval.remove[j]
            retval@i = [lastX, lastY]
            last = last - 1
         }
      }
      i = i + 1
   }
   return retval
}


//println[maximize[{|x| -x^2+2x-3}, 0, 4 pi]]
//setPrecision[15]
//println[join["\n",maximize[{|x| sin[x]}, 0, 4 pi, 1e-15]]]
println[]
println[join["\n",minimize[{|x| sin[x]}, 0, 4 pi, 1e-15]]]
println[]
//println[join["\n",maximize[{|x| -(x^2)}, -4, 4, 1e-30]]]

use sun.frink
println[]
//println[join["\n",minimize[{|date| sunMoonAngularSeparation[date, 40 degrees North, 105 degrees West]}, #2012-05-20 7:30 PM #, #2012-05-20 8:00 PM#, 1s]]]


View or download maximize.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17646 days, 9 hours, 35 minutes ago.