9
$\begingroup$

I need to evaluate a function func on a range of integers (NOT recursively, just sequentially) until the result is stable enough. I can write a not very clever step-by-step program using While:

Module[{variation,list={},j=1},

  AppendTo[list,func[1]];
  variation = Abs@First@list;

  While[variation>threshold,
    j++; AppendTo[list,func[j]]; variation = Abs[list[[j]]-list[[j-1]]];];

list    
]

but I would prefer using proper functional programming.

Ideally, I would like a function that I can map on a Range[nMax], which stops on its own if the result has become stable (say, to within some absolute threshold) and returns the list of numbers. In this way I could also control the maximum number nMax of evaluations:

stableList[func,threshold]/@Range[nMax]

I suspect I have to use some clever trick (memoization?) to allow Map to compare previous values. Any suggestions?

$\endgroup$
  • 2
    $\begingroup$ There's a FixedPoint function. $\endgroup$ – faysou Mar 16 '15 at 11:53
  • 2
    $\begingroup$ I don't need recursion though... $\endgroup$ – Ziofil Mar 16 '15 at 17:36
4
$\begingroup$

A possible implementation of stableList that evaluates fun sequentially:

stableList[func_, threshold_, nMax_: Infinity] := 
  NestWhileList[{func[#], #} &[Last@# + 1] &, {func[1], 1}, 
    Abs[First@#1 - First@#2] > threshold &, 2, nMax - 1][[All, 1]]

Usage example:

fun[x_] := 1/x

stableList[fun, 0.01]

{1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11}

stableList[fun, 0.00001, 7]

{1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7}

An alternative implementation using FixedPointList instead of NestWhileList:

stableList2[func_, threshold_, nMax_: Infinity] := 
  FixedPointList[{func[#], #} &[Last@# + 1] &, {func[1], 1}, nMax - 1, 
    SameTest -> (Abs[First@#1 - First@#2] < threshold &)][[All, 1]]

And one that is based on using While, but is more efficient than yours

stableList3[func_, threshold_, nMax_: Infinity] := Reap[
  Module[{n = 1, difference = Infinity, newResult},
    newResult = Sow[func[n++]];
    While[n <= nMax && difference > threshold,
      difference = Abs[newResult - (newResult = Sow[func[n++]])];
    ];];][[2, 1]]
$\endgroup$
  • $\begingroup$ Overloading the stableList function for the case without a nMax will be faster than defining nMax as an optional argument with default value Infinity as checking n <= nMax isn't needed if there is no nMax. $\endgroup$ – Karsten 7. Mar 17 '15 at 17:54
8
$\begingroup$

FixedPointList[f, expr, n, SameTest -> g[e1,e2]] evaluates f[expr] recursively for at most n times. It also stops if g[e1,e2] returns true. e1 and e2 being the most recent values. It is very similar to NestWhileList. There is also a version that returns only the last value, called FixedPoint.

Map is not meant to be used sequentially, and so you can't break out of it the way you can break out of a loop. This makes it unsuitable. You can build your own function using plain recursion like this:

f[e1_, e2_, e3_] := If[
  Abs[e1 - e2] > 0.001,
  f[(e1 + 2/e1)/2, e1, {(e1 + 2/e1)/2, e3}],
  Flatten@{e1, e3}
  ]

f[2.0, 10.0, {}]

(* Out: {1.41421, 1.41421, 1.41422, 1.41667, 1.5} *)

This particular function approximates the square root of 2.


EDIT. This is a comment on Jacob's answer that is to large to be posted as a comment. I would write his first function like this:

stableList[threshold_, f_] := Reap@Module[{prev, next, j},
   For[j = 1, j < 3 || Abs[next - prev] > threshold, j++,
    prev = next;
    next = Sow@f@j;
    ]
   ]

It is more readable in my opinion, and also equivalent

thres = 0.00001;
func = 1/# &;
stableList[thres, func][[2, 1]] == intermitThresDif[thres]
(* Out: True *)

I'm not convinced that a procedural approach is advantageous to a functional approach for this problem, especially not in Mathematica. Functionally this code could for example be implemented as

func2[{x_, y_}] := {1/y, y + 1}
stableList2[threshold_, f_] := First /@ NestWhileList[f, {1, 2}, Abs[First@# - First@#2] > threshold &, 2]

which I think is even more succinct, and equivalent

thres = 0.00001;
func = 1/# &;
stableList2[thres, func2] == intermitThresDif[thres]
(* True *)

If you use a procedural approach you should try to use Compile, which can indeed be much faster than a functional approach (without Compile it should be slower; Jacob tested and it turns out they are equally fast). [Deleted text re. AppendTo here]. If you compile your code you can use Internal`Bag instead of Append and AppendTo. But it's better (if you intend to compile the code) to create a list equal in length to a max number of attempts and then use Part to update elements of that list (this is not good if you are not compiling the code.) See Jacob's answer for more information.

$\endgroup$
  • $\begingroup$ Hey Pickett, thanks for your response. You make some good points. For example I like how you managed to put everything in the body of the For loop and avoided nonsense like variation = threshold+1. I suppose I should have done this but I had a weird preference for While in cases like this. I was about to implement a better functional solution when I saw your edit and I like your use of NestWhileList here. I guess my procedural code should be a bit slower because I have more intermediate definitions and I suppose NestWhileList avoids the need to ... evaluate a bunch of Sows. But... $\endgroup$ – Jacob Akkerboom Mar 16 '15 at 15:51
  • $\begingroup$ That difference in speed should be small compared to compiled vs uncompiled code. Also my advice about Append(To) extends to the compiled case, I suppose I will post code to show what I mean. I guess I was silly to say programming functionally is not too efficient in cases like these. Really what I want to say is that procedural code like this is not invalid, especially because of the possibility to compile. $\endgroup$ – Jacob Akkerboom Mar 16 '15 at 15:51
  • $\begingroup$ On my machine, your two functions in your edit are about equally fast. $\endgroup$ – Jacob Akkerboom Mar 16 '15 at 16:09
  • $\begingroup$ @JacobAkkerboom You are right, I don't know why I wrote the thing I wrote about Append, I wasn't thinking straight. Thanks for testing the speed of the solutions. If we take a broader perspective on things and focus not just on this particular problem, the advantage of recommending a functional style to other users is that Map, Table etc. have auto-compilation, which the procedural style doesn't. Simply put, functional programming allows for more optimizations which is why I like to use it whenever I can unless I know I will be compiling the function myself. $\endgroup$ – C. E. Mar 16 '15 at 16:40
  • $\begingroup$ I think your point about auto compilation is very good and I guess this even works for Table, as I suppose the first argument of Table can only be compiled if it has no side effects, which is its functional use case. I don't fully agree with everything you say, but I guess it is largely a matter of preference for procedural vs functional programming (in MMA). I would like to dive into Haskell one day, where (auto) compilation is much more advanced, giving functionally represented solutions much more speed. $\endgroup$ – Jacob Akkerboom Mar 19 '15 at 16:10
3
$\begingroup$

I think your code is not too bad. In my opinion, the main thing to learn is not to use Append(To) in a loop. Here is an improved version

func = 1/# &;
intermitThresDif =
  Function[
   {threshold},
   Reap[
     Module[
      {variation, prev, next, j = 1},
      prev = func@1;
      Sow@prev;
      variation = threshold + 1;
      While[
       variation > threshold,
       j++;
       next = func@j;
       Sow@next;
       variation = Abs[next - prev];
       prev = next
       ]
      ]
     ][[2, 1]]
   ];

Note that is usually easier to convert procedural code to a to compiled function, which is an advantage over a functional approach.

Edit: I used to have some functional code here, but that was a bit silly


Here is some code that shows that even in the case where you Compile, using Append in a loop is bad. I also show an alternative, which makes use of undocumented functions, but is still "pretty standard".

cfu1 =
 Compile[
  {{nn, _Integer}}
  ,
  Block[{bag},
   bag = Internal`Bag@ Most@{1};
   Do[Internal`StuffBag[bag, 0], {nn}];
   Internal`BagPart[bag, All]
   ]
  ,
  CompilationTarget -> "C"
  ];
cfu2 =
 Compile[
  {{nn, _Integer}}
  ,
  Block[{list},
   list = Most@{0};
   Do[AppendTo[list, 0], {nn}];
   list
   ]
  ,
  CompilationTarget -> "C"
  ];

cfu1[10^6] // Timing // First
cfu2[10^4] // Timing // First
0.023427
0.053426
$\endgroup$
  • $\begingroup$ I added some comments at the end of my answer in response to some of the things you wrote here. $\endgroup$ – C. E. Mar 16 '15 at 14:38
2
$\begingroup$

This example from the documentation may be instructive for approximation of $\sqrt{2}$:

NestWhileList[(# + 2/# )/2 &, 1, Abs[#1 - #2] > 0.001 &, 2]

->{1, 3/2, 17/12, 577/408, 665857/470832}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.