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 19944 days, 21 hours, 37 minutes ago.