Carl Love

Carl Love

28070 Reputation

25 Badges

13 years, 28 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

Your test function has a discontinuous noise component: a normal random variable of fixed variance. I don't think that that's what Calvin's algorithm was intended for. Here is an example with a continuous noise component, a Wiener process:

Fi:= Finance:  CF:= CurveFitting:
npts:= 10^4:
WS:= Fi:-SamplePath(
    Fi:-WienerProcess(<<.1>>)(t), t= 0..1, timesteps= npts-1
)[1,1]:
T:= rtable(1..npts, k-> (k-1)/(npts-1), datatype= float[8]):
WF:= t-> CF:-ArrayInterpolation(T, WS, t, uniform):
W:= t-> (t-0.5)^2 + WF(t):
CodeTools:-Usage(Calvin(W, 350));
memory used=457.37MiB, alloc change=0 bytes, cpu time=7.61s, 
real time=7.47s, gc time=437.50ms

            [0.702660743159389, -0.527027281446182]

plot(W, 0..1);

The computations take much longer because I discretized the Wiener process to a very large number of points, 10,000.

@ ArrayTools:-Insert was introduced in Maple 2016. Here's a retrofit for the above code that doesn't use it:

Calvin:= proc(
    X::procedure,
    n::posint, 
    gamma::procedure:= (n-> sqrt(n+1)/ln(n+2))
)
description "Minimize continuous function with random noise on 0..1";
option
    `Reference: `,
        `https://www.sciencedirect.com/science/article/pii/S0885064X01905746 `,
        `James M. Calvin, "A One-Dimensional Optimization Algorithm and Its `,
        `Convergence Rate under the Wiener Measure", _Journal of Complexity_ `,
        `17, 306-344 (2001)`,
    `Maple author: Carl Love <carl.j.love@gmail.com> 2020-Aug-27`
;     
local 
    T:= rtable(1..n+2, [0, 1], datatype= float[8], subtype= Vector[row]),
    XT:= X~(T), Y:= copy(XT),
    t, tmin, M, z:= evalf(gamma(0)), tau:= 1, k, i,
    Insert:= proc(V::Vector, i::posint, k::posint, x)
        V[i+1..k+2]:= V[i..k+1];
        V[i]:= x
    end proc    
;
    (tmin,M):= `if`(XT[1] < XT[2], [T[1],XT[1]], [T[2],XT[2]])[];
    for k to n do
        Y[..k+1]:= XT[..k+1] -~ (M - z);
        i:= 1 + max[index]( (T[2..k+1] - T[..k]) /~ Y[2..k+1] /~ Y[..k]);
        t:= T[i-1] + (T[i]-T[i-1])*Y[i-1]/(Y[i-1]+Y[i]);
        tau:= min(tau, t - T[i-1], T[i] - t);
        z:= min(z, evalf(gamma(k))*sqrt(tau));
        Insert~([T, XT], i, k, [t, evalf(X(t))]); 
        if XT[i] < M then (tmin,M):= (t, XT[i]) fi;
    od;
    [tmin, M]
end proc
:
N:= Statistics:-Sample(Statistics:-RandomVariable(Normal(0, 0.1))):
S:= Array(1..1, datatype= float[8]):
W:= t-> 1 + (t - 0.5)^2 + N(S)[1]:
CodeTools:-Usage(Calvin(W, 350)); 
memory used=20.76MiB, alloc change=0 bytes, cpu time=141.00ms, 
real time=146.00ms, gc time=0ns

             [0.495821549079493, 0.710534027290144]

Regarding the precision, or Digits setting: It's ridiculous to use any precision higher than hardware floats for this algorithm or for an objective function with this much noise.

@dharr In this case, I think that Nu stands for the Nusselt number. There are so many named parameters in boundary-layer fluid dynamics that they usually go by two-letter abbreviations.

@nguyenhuyenag The embedded assignment was introduced in Maple 2018. By using it, I was able to avoid having the same statement, vars:= indets(f) minus kv, appear both before the loop and inside the loop. That's just a stylistic preference, and has little effect on efficiency.

@ In your code, you've used a sequence theta(n) for what Calvin calls gamma(n). You've defined theta(n) as 1/(1+n), but this does not have the properties listed at the bottom of p 308: "an increasing sequence of positive numbers such that gamma(n) => infinity, gamma(n)^2*ln(gamma(n))/n => 0, and gamma(n)/gamma(n+1) => 1." So, this is your error, or perhaps one of your errors.

@ Okay, I can read the paper now. But I don't see a formal listing of an algorithm or any pseudocode. What part of the paper have you translated into Maple?

I noticed that you added a try...catch statement. What type of errors has it been catching?

@ Your Maple procedure WP declares several variables that aren't used in the code. That's a sign that you missed something. The variable a is initialized for each iteration of the main loop, but it isn't otherwise used.

Do the simplifications depend on anything that's not explicitly in the expressions being simplified? I'm primarily thinking of assume commands.

What's the value of nops(dummy_UU1)?

What's the average CPU time per simplification for the serial code? In other words, what's the number returned by

CodeTools:-Usage(map(simplify, dummy_UU1), output= cputime, quiet)/nops(dummy_UU1);

Have you tried Grid:-Map without the tasksize option?

You didn't reply to my extensive Answer to one of your other Questions: https://www.mapleprimes.com/questions/230277-What-Is-The-Conception-Of-process

Nor did you respond to the two emails that I sent you.

I spent several hours researching and writing that Answer.

@ The body of the Elsevier / Science Direct article that you linked to is not accessible without paying a fee.

Does your Maple version produce the wrong results? Or does it simply take too long? If it's the latter, I can fix it without needing to see the original. 

If this is an attempt by you to find the longest path, it won't work.

@vs140580 From the Wikipedia article "Hypercube graph", I learned that for distinct vertices u and v in HyperCubeGraph(n) there is a Hamiltonian path between them iff their Hamming distance is odd. Thus, if the distance is odd, then the length of the longest path is 2^n - 1. If the Hamming distance is even, I don't know, but I suspect there's an equally simple formula.

@vs140580 If you'd want to restrict attention to certain highly symmetric special cases, such as hypercubes, then a fast answer is possible. The vertices of HyperCubeGraph(n) can be represented as n-lists of 0s and 1s. Because of the symmetry, the relationship between any two vertices is (up to isomorphism) entirely determined by the number of positions that are different in their 0-1 representations (known as their Hamming distance). Thus, there are only n possible relationships between pairs of distinct vertices. Thus there are only n possible values for the length of the longest path.

@sand15 The issue is that Sample uses an hfloat array by default. If you map those integerizing functions over an hfloat array, the result is still an hfloat array, and thus can only contain hfloats. This allows those commands to work under evalhf. But if you use an sfloat array, as you did, then mapping those functions changes the datatype of the array, because now there's no worry about evalhf.

But if you first convert to a list, as I did, before mapping the integerizing function, it works fine. As my code shows, this is many times faster than using an sfloat array.

You asked:

  • Couldn't  it be interesting that these three functions also operate on hfloats?

Yes, it is, and exploiting that fact can lead to remarkable speed gains (20X - 30X). Hardware floats can exactly represent and compute with integers under 2^50 (or so). This is one reason why the LinearAlgebra:-Modular package is so fast. For the functions that can work with evalhf, see ?evalhf,fcnlist.

@vs140580 The longest-path problem for general graphs and digraphs is known to be NP-hard, meaning that there's provably no fast algorithm for it. You need to check every path.

First 168 169 170 171 172 173 174 Last Page 170 of 709