(* :Name: Optimize` *) (* :Title: Expression Optimization *) (* :Author: Terry Robb *) (* :Copyright: Copyright 1993 T.D.Robb *) (* :Summary: Basic structural optimization of expressions. See ?Optimize *) (* :Keywords: Optimize, efficiency, compiling *) (* :Version: 1.2 *) (* :Discussion: See examples at the end of this file. *) BeginPackage["Optimize`"] Optimize::usage = StringJoin[ "Optimize[expr] optimizes an expression and returns a HoldBlock object. ", "An example usage is: Optimize[Sin[x] + Cos[Sin[x]]] ", "and it returns an object that can be used when defining a function, ", "for example: f[x_Real] := Optimize[Sin[x] + Cos[Sin[x]]]; ", "or it can be used as an argument inside Compile or Function as in: ", "g = Compile[x, Optimize[Sin[x] + Cos[Sin[x]]]];. ", "An optimized object can also be evaluated directly using ", "Normal or ReleaseHold as in: ", "Normal[Optimize[Sin[x] + Cos[Sin[x]]] /. x->2.2]." ]; HoldBlock::usage = StringJoin[ "HoldBlock[vars, body] represents an optimized expression. ", "The commands SetDelayed, Compile, Function, Normal and ReleaseHold ", "all know about HoldBlock and will convert it to Block as necessary." ]; Cost::usage = StringJoin[ "Cost[expr] estimates the cost of evaluating an expression by simply ", "counting the number of basic elementary operations that are used. ", "An example usage is: Cost[Sin[x] + Cos[Sin[x]]] which returns ", "Cos + Plus + 2*Sin. Cost also works with optimized objects, so if ", "you apply (Cost[#] - Cost[Optimize[#]])& to an expression you will ", "see the operations that might be saved." ]; Begin["`Private`"] Optimize[expr_, Times->True] := Block[{TIMES, a, b}, Optimize[expr //. Times[a_, b__] :> TIMES[a, Times[b]]] /. TIMES->Times]; Optimize[expr_] := Block[{COMP,dups,objs,unis,offs,vars,ruls,sets,body,i,d,v}, Attributes[COMP] = {HoldAll}; dups = Table[objs = Level[expr, {-i}, Hold (*, Heads->True *)]; unis = List @@ Map[Literal, Union[objs]]; Select[unis, (* Count[objs, #]>1 &, *) Length[Position[objs, #, {1}, 2, Heads->False]]>1 &, Length[objs] - Length[unis]], {i, 2, Depth[expr] - 1}]; dups = DeleteCases[dups, {}]; If[dups==={}, Return[expr]]; offs = Partition[FoldList[Plus, 0, Map[Length, dups]], 2, 1]; vars = GETVARS[Last[Last[offs]]]; ruls = Apply[Rule, Transpose[{Flatten[dups], vars}], {1}]; ruls = Map[Reverse[Take[ruls, # + {1,0}]]&, offs] // Reverse; sets = MapIndexed[Fold[ReplaceAll, #1, Drop[ruls, First[#2]]]&, ruls, {2}]; sets = Reverse[COMP @@ Flatten[sets]]; body = Append[sets /. Rule[_[d_], v_] :> Set[v, d], If[Length[ruls]===1, expr /. ruls[[1]], expr //. Flatten[ruls]]]; HoldBlock @@ {vars, body} /. COMP -> CompoundExpression ]; Optimize /: (op:SetDelayed|Compile|Function)[args_, o_Optimize, opts___] := op[args, Evaluate[o], opts]; Attributes[HoldBlock] = {HoldAll}; HoldBlock /: (ReleaseHold|Normal)[HoldBlock[vars_,body_]] := Block[vars,body]; HoldBlock /: (func_ := HoldBlock[vars_,body_]) := (func := Block[vars,body]); HoldBlock /: (op:Compile|Function)[args_, HoldBlock[vars_,body_], opts___] := op[args, Block[vars, body], opts]; Cost[HoldBlock[vars_, body_] | body_] := Plus @@ Map[Part[#, 1, 0]&, Map[Literal, Level[Hold[body], -2, Hold]]] /. CompoundExpression -> 0; GETVARS[n_] := ( If[n>NVARS, NVARS = Length[OPTVARS = Join[OPTVARS, Table[Unique["System`$"(*System`O, Protected*)], {n-NVARS+64}]]]]; Take[OPTVARS, n] ); If[!ListQ[OPTVARS], OPTVARS = {}; NVARS = 0; GETVARS[64]]; SetAttributes[{Optimize, HoldBlock, Cost}, {ReadProtected, Protected}]; End[] EndPackage[] (* :Examples: expr = Nest[(#+Sin[#])&, x, 10]; f1[x_Real] := Optimize[expr]; f2 = Compile[x, Optimize[expr]]; Plot[expr//Evaluate, {x, 0, 0.01}];//Timing (* 2.3 Seconds *) Plot[f1[x], {x, 0, 0.01}];//Timing (* 0.5 Seconds *) Plot[f2[x], {x, 0, 0.01}];//Timing (* 0.1 Seconds *) ?f1 expr = Expand[(Sin[x*y] + ArcTan[x] + Log[y] + Cosh[x^-y]^5)^10]; oexpr = Optimize[expr]; (expr /. {x->1.1, y->0.5})//Timing (* 1.0 Seconds *) Normal[oexpr /. {x->1.1, y->0.5}]//Timing (* 0.3 Seconds *) Cost[expr] Cost[oexpr] roots = z /. Solve[z^4 - 12*z^3 - 16*z^2 + 4*z + c == 0, z]; Optimize[roots] expr = (1+x*y)/(1-x*y); Optimize[expr] FullForm[expr] Optimize[expr //. Times[a_, b__] :> TIMES[a, Times[b]]] /. TIMES->Times Optimize[expr, Times->True] *) (* :Limitations: Expressions that contain objects with a HoldAll attribute cause problems. *) (* :Inefficiencies: The implementation used here works well if only a smallish number (say less than 500) of the optimize O$* variables are ultimately needed. The worst case is if the number of O$* variables is half the number of components in a very large expression, but such cases are rare. The code above works well for all the common cases though, and the only obvious slow parts are 1) the part where the Select[] is used, and 2) the part where the expr//.ruls is used. Part 1) could easily be sped up, for example by using something like: Clear[count]; Scan[(count[#]=0)&, unis]; Scan[(count[#]++)&, objs]; Select[unis, count[#]>1 &, Length[objs] - Length[unis]] which would asymptotically run much faster for the worst cases, but is a lot slower (even asymptotically) for the common cases. One could dynamically decide when to use that alternative though. Another way to speed up the above is to replace the lines: ------------------- unis = List @@ Map[Literal, Union[objs]]; Select[unis, (* Count[objs, #]>1 &, *) Length[Position[objs, #, {1}, 2, Heads->False]]>1 &, Length[objs] - Length[unis]], ------------------- with the lines ------------------- objs = List @@ Map[Literal, objs]; If[OddQ[Length[objs]], AppendTo[objs, Null]]; Intersection @@ Transpose @ Partition[Sort[objs], 2], ------------------- The above change speeds up the hard cases a lot, but it seems to slow down by just a little bit the common cases. Part 2) could be sped up by using Dispatch[ruls] which would hash the rules and so make the operation linear not quadratic. However there is a bug in Dispatch (in both Mma2.1 and Mma2.2) which stopped me doing that :-(. *) (* :Miscellaneous: Here's my first (but slower) version of Optimize: Optimize[expr_] := Block[{COMP,dups,objs,unis,vars,ruls,sets,body,i,d,v}, Attributes[COMP] = {HoldAll}; dups = Table[objs = Level[expr, {-i}, Hold]; unis = List @@ Map[Literal, Union[objs]]; Select[unis, Count[objs, #]>1 &, Length[objs] - Length[unis]], {i, 2, Depth[expr] - 1}] // Flatten; vars = GETVARS[Length[dups]]; If[vars==={}, Return[expr]]; ruls = Apply[Rule, Transpose[{dups, vars}], {1}] // Reverse; sets = COMP @@ MapIndexed[(#1 //. Drop[ruls, #2[[1]]])&, ruls] // Reverse; body = Append[sets /. Rule[_[d_], v_] :> Set[v, d], expr //. ruls]; HoldBlock @@ {vars, body} /. COMP -> CompoundExpression ]; *) Null;