maximize.frink

Download or view maximize.frink in plain text format


use allTransforms.frink
use functionUtils.frink

/** This file contains routines to minimize or maximize a function over a
    certain range using interval arithmetic techniques.  And, maybe, even
    better, using symbolic differentiation techniques! */


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

minimize[f, minX, maxX, minXDiff] := intervalOptimize[f, minX, maxX, minXDiff, -1]

intervalOptimize[f, minX, maxX, minXDiff, factor=1] :=
{
   widthStep = maxX - minX
   fullX = new interval[minX, minX + widthStep/2, maxX]
   
   // solutions is an array of arrays.  Each array contained is a contiguous
   // range of intervals of still-valid solutions.
   solutions = [[minX, fullX, maxX]]
//   solutions.push[[fullX]]

   fullY = factor * f[fullX]
   lowerBound = infimum[fullY]
   hardBound = lowerBound

//   println["lowerBound is $lowerBound"]

   // Try the endpoints because most functions are monotonic
   minY = factor f[minX]
   ly = infimum[minY]
   if ly < hardBound
      hardBound = ly
   
   maxY = factor f[maxX]
   ly = infimum[maxY]
   if ly < hardBound
      hardBound = ly

   while (widthStep > minXDiff)
   {
      widthStep = widthStep / 2
//      println["input solutions are: $solutions"]
      newSolutions = new array

      for solutionBlock = solutions
      {
         newSolutionBlock = undef
         SOLUTION:
         for solution = solutionBlock
         {
            y = factor f[solution]
            ly = infimum[y]
//            println["ly is $ly"]
            if ly > lowerBound
            {
               lowerBound = ly
//               println["new lower bound is $lowerBound"]
            }
            
            // If y is possibly greater or equal to the lower bound, then this
            // interval may contain a solution.
            if y PGE lowerBound
            {
               //      println["width is $width"]
               
               // Bifurcate the x interval and add the 2 new intervals to the
               // newSolutionBlock
               lower = infimum[solution]
               upper = supremum[solution]
               width = (upper-lower)/2
               halfWidth = width/2
               
               mid = lower + width
               midLower = lower + halfWidth
               nx = new interval[lower, midLower, mid]
               ny = factor f[nx]
//               println["nx is $nx, ny is $ny"]
               my = mainValue[ny]
               if my == undef
                  println["No mid for $nx, $ny"]
               if my != undef && my > hardBound
               {
//                  println["hardBound was $hardBound"]
//                  println["hardBound getting set by $my"]
                  hardBound = my
                  lowerBound = hardBound
//                  println["new hard bound is $hardBound for $nx (1)"]
               }
               
               if ny PGE hardBound and ny PGE lowerBound
               {
                  ly = infimum[ny]
                  if ly > lowerBound
                  {
                     lowerBound = ly
//                     println["new lower bound is $lowerBound (1)"]
                  }

//                  println["Pushing $nx (1)"]
                  // Construct a new solution block if it doesn't exist
                  if newSolutionBlock == undef
                     newSolutionBlock = new array
                  
                  newSolutionBlock.push[nx]
               } else
               {
                  if newSolutionBlock != undef
                  {
                     newSolutions.push[newSolutionBlock]
                     newSolutionBlock = undef
                  }
               }

               if lower == upper
                  next SOLUTION
               
               midUpper = mid + halfWidth
               nx = new interval[mid, midUpper, upper]
               ny = factor f[nx]
               my = mainValue[ny]
               if my != undef && my > hardBound
               {
//                  println["hardBound was $hardBound"]
//                  println["hardBound getting set by $my"]
                  hardBound = my
                  lowerBound = hardBound
//                  println["new hard bound is $hardBound for $nx (2)"]
               }
               
               if ny PGE hardBound and ny PGE lowerBound
               {
                  ly = infimum[ny]
                  if ly > lowerBound
                  {
                     lowerBound = ly
//                     println["new lower bound is $lowerBound (2)"]
                  }
//                  println["Pushing $nx (2)"]
                  // Construct a new solution block if it doesn't exist
                  if newSolutionBlock == undef
                     newSolutionBlock = new array

                  newSolutionBlock.push[nx]
               } else
               {
                  if newSolutionBlock != undef
                  {
                     newSolutions.push[newSolutionBlock]
                     newSolutionBlock = undef
                  }
               }
            }
         }
         
         if newSolutionBlock != undef
            newSolutions.push[newSolutionBlock]
      }
      solutions = newSolutions
      println["lowerBound is $lowerBound"]
      println["There are " + length[flatten[solutions]] + " solutions"]
//      println["Solutions: $solutions"]
      println["width is $widthStep"]
   }

//   println[solutions]
   s = new array
   for sols = solutions
   {
      maxMain = undef
      xMax = minX
   
      s1 = undef
      if (length[sols] > 0)
      {
         for i = 0 to length[sols] - 1
         {
            nx = sols@i
            ny = factor f[nx]
//            println["Prospective solution " + sols@i + " = " + ny]
            
            my = mainValue[ny]
            if my != undef && (maxMain == undef || my > maxMain)
            {
               maxMain = my
               if mainValue[nx] != undef
               {
                  xMax = mainValue[nx]
//                  println["New xMax = $xMax at $nx"]
               }
            }
            
            if ny PGE lowerBound
               if s1 == undef
                  s1 = sols@i
               else
                  s1 = union[s1, sols@i]
         }
      }

      if xMax != undef
         s1 = new interval[infimum[s1], xMax, supremum[s1]]
      
      if s1 != undef
         s.push[[s1, f[s1]]]
   }
   
   return s
}


symbolicOptimize[f, minX, maxX, factor=1] :=
{
   args = functionArgumentsAsSymbols[f]
   if length[args] != 1
   {
      println["Wrong number of arguments"]
      return undef
   }

   println["Arguments are $args"]
   derivativeFunc2 = makeDerivative[f]
   println["Before solve:   $derivativeFunc2"]
   derivative2 = transformExpression[derivativeFunc2]
   println["Derivative is:  $derivative2"]

   solve = makeSolve[derivative2, 0, args@0]
   println["Before solve: $solve"]
   solved = transformExpression[solve]
   println["After solve: $solved"]
   type = type[solved]
   println["Type is " + type]
   if (type == "Solve")  // Is this a solved expression with ===
   {
      right = getChild[solved, 1]
      println["Right is " + right]
      solution = eval[right]
      println["Solution is " + solution]
      return solution
   } else
      return undef
}


//setPrecision[15]
//println[maximize[{|x| -x^2+2x-3}, 0, 4 pi]]
//println[maximize[{|x| x}, 0, 10, 1e-20]]
//println[]
//println[join["\n",maximize[{|x| x}, 0, 10, 1e-20]]]
//println[]
//println[join["\n",maximize[{|x| x sin[x]}, -15, 15, 1e-8]]]
//println[]
//println[join["\n",maximize[{|x| x cos[x]}, -23, 23, 1e-7]]]
//println[]
//println[join["\n",minimize[{|x| x^x}, 0.01, 4 pi, 1e-8]]]
//println[]
//println[join["\n",maximize[{|x| x^2}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",minimize[{|x| -x^2}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",maximize[{|x| -x^2}, -4., 4., 1e-20]]]
//println[]
//println[join["\n",maximize[{|x| -abs[x^2]}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",maximize[{|x| abs[x^2]}, -4., 4., 1e-10]]]
//println[]
//println[join["\n",maximize[{|x| x cos[x] + sin[x]}, -10, 10, 1e-8]]]

/*
use sun.frink">sun.frink
println[]
println[maximize[{|date| sunMoonAngularSeparation[date, 40 degrees North, 105 degrees West]}, #2019-01-20 7:30 PM #, #2019-01-20 11:59 PM#, .1 s]]
*/


/*
showApproximations[false]
symbolicMode[true]
a = symbolicOptimize[{|x| sin[2 x]}, -pi, pi]
*/


Download or view 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 20136 days, 5 hours, 22 minutes ago.