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 17591 days, 17 hours, 45 minutes ago.