Integrate1D with MultiThread

Share icon

Hi all,

I have a function which performs multiple Integrate1D's (one for each point in a large 4D matrix). Needless to say it's very slow, taking up to 48+ hours depending on the size of the matrix. This seems like a pretty good opportunity to use MultiThread, but apparently Integrate1D isn't ThreadSafe.

Does anyone know of a workaround for this?

Any help would be appreciated.

Cheers,
Hugh
How big is the 4D matrix? Also, can you post the function that is given to integrate1D?

May 3, 2011 at 01:16 am - Permalink

Thanks for the response.

The matrix varies in size, depending on how much time I can afford to let it run for. Generally from 21^4 to 151^4 (the bigger the better).

Anyway, here's the function I integrate to:

Function/C LatScat(t)<br />
    Variable t<br />
    NVAR Gu0=Gu0<br />
    NVAR Gu1=Gu1<br />
    NVAR Gu2=Gu2<br />
    NVAR Gu3=Gu3<br />
    NVAR Gh1=Gh1<br />
    NVAR Gh2=Gh2<br />
    NVAR Gm=Gm<br />
    <br />
    Variable/C i=cmplx(0,1)<br />
    Variable U=(Gu0*cos(2*Pi*t))+(Gu1*cos(4*Pi*t))+(Gu2*cos(6*Pi*t))+(Gu3*cos(8*Pi*t))<br />
<br />
    return exp(2*Pi*i*((Gh1*U)+(Gh2*U)+(Gm*t)))<br />
End


The variables Gh1, Gh2 & Gm are constant across the matrix, wheras Gu0-Gu3 vary (they're linear with the row,column,layer,chunk index values). That is to say, my matrix plots LatScat vs Gu0 & Gu1 & Gu2 & Gu3.

Cheers,
Hugh

May 3, 2011 at 02:00 am - Permalink

Another thing springs to mind. What does Mathematica have to say about the integral?

Or, can you evaluate the function between the two limits, then use integrate (which is threadsafe)?

May 3, 2011 at 05:36 am - Permalink

andyfaff wrote:
Another thing springs to mind. What does Mathematica have to say about the integral?

Or, can you evaluate the function between the two limits, then use integrate (which is threadsafe)?


Mathematica has a hard time with that particular integral (you get the form exp[ix sin[x]]). Generally this is solved by approximating with a limited number of terms of an infinite series of Bessell functions, but Mathematica gives me a very long solution after about 5mins for a single term. And that's still a solution approximated by Bessel terms. Numerical is definately the way to go here (plus all my other work on this particular problem is on Igor anyway).

Actually, I originally tried using integrate, but found the output wave was of the same dimension as the input wave (i.e. it should be a single number, right?). Clearly I was doing something wrong.

I think having another go with integrate may be the solution here, as it means I can do away with all the global values since I can just use them when I evaluate the function to begin with, right?

Are there any differences in the quality of the integration methods used in integrate vs integrate1D that you know of?

Thanks for your help!

Cheers,
Hugh

May 3, 2011 at 06:27 am - Permalink

hugh50935][quote=andyfaff wrote:
Actually, I originally tried using integrate, but found the output wave was of the same dimension as the input wave (i.e. it should be a single number, right?). Clearly I was doing something wrong.


If I remember correctly, the output from Integrate is not a single value but rather a wave, where the value of each point is the cumulative integral up to that point. So to get the total integral you need to take the last point of the output.

Integrate and Integrate1D differ in that the first integrates only tabulated values - Integrate does not know of any user-defined function. Integrate1D, on the other hand, does. So you could try to do each integration by tabulating all the values first, and then calling Integrate. This should very easy to parallellize using a wave assignment and the MultiThread keyword. Just be careful with the global variables, and consider replacing them with function arguments instead.

Unfortunately, this also means that you will have to decide at which values to tabulate the function. Integrate1D takes care of this for you, meaning that it can potentially make fewer function calls and/or obtain a higher accuracy. Integrate1D would also be conceptually much cleaner. I hope WaveMetrics chimes in on why it's not threadsafe, and whether they can address this to help you out.

Another way to obtain a speedup is to try to rewrite the integrated function in an XOP, which might make it much faster. But again the problem here is to get the extra information in there (i.e. the global variables). This can be hacked around, but it will not be pretty.

May 3, 2011 at 06:50 am - Permalink

741 wrote:
I hope WaveMetrics chimes in on why it's not threadsafe, and whether they can address this to help you out.

AG, the author of Integrate1D is on vacation right now. I'm sure he will respond when he gets back next week.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 3, 2011 at 09:25 am - Permalink

Have you played around with the Integrate1D options?
Especially options and count look like they influence the required time.

Your user defined function is also still pretty complicated. Have you checked that it runs fast enough?
For example you could precompute some of the terms, e.g.
(Gh1*U)+(Gh2*U) = (Gh1+Gh2)*U

and you can also precompute Gh1+Gh2 therefore skipping one multiplication. But maybe this might not be the real issue here.

Or maybe there is one fitting addition theorem which lets you calculate cos(a*PI*t) less than three times.

May 3, 2011 at 02:20 pm - Permalink

Hi all,

Thanks for all your suggestions.

I rewrote my function using Integrate, and tried to cut out superfluous multiplications. I made a 1D wave, integrated it and took the last value. 100 points in the wave seems to get pretty close to Integrate1D, which is nice.

Doing this means I could make my variables local again, and reduce the total number of functions, leaving me with this:

ThreadSafe Function F(n,h1,h2,h3,m,k,u0,u1,u2,u3,mode)
    Variable n,h1,h2,h3,m,k,u0,u1,u2,u3,mode
    Variable/C i=cmplx(0,1)
    Variable U
    Make/C/O/N=(n+1) wdF=0
   
    Variable t
    for(t=0;t<=n;t+=1)
        U=(u0*sin(2*Pi*(t/n)))+(u1*sin(4*Pi*(t/n)))+(u2*sin(6*Pi*(t/n)))+(u3*sin(8*Pi*(t/n)))
        if(mode==0)
            wdF[t]=exp(2*Pi*i*(((h1+h2)*U)+(m*(t/n))))
        else
            wdF[t]=exp(2*Pi*i*(((h3+(m*k))*U)+(m*(t/n))))
        endif
    endfor
   
    Integrate wdF  
    Variable F=magsqr(wdF[n])
    KillWaves wdF
    return F
End


So, I just ran this on my desktop with a 20^4 data set, which wasn't a whole lot faster than the Integrate1D method (about 10-20% I guess). But on my less cluttered sever (4-core Xeon 5570) multithreading reduced the solving time by a factor of 4. The giveaway was that the CPU load went from 13% to 100% with multithreading. Good work, guys!

I'm not entirely happy with my integration (i.e. the bit in the FOR loop). Is there a way to streamline this using wave assignments? I feel like I want to replace 't' with 'p', but I don't know how since the Variable U is a function of 't' as well.

Cheers,
Hugh

May 3, 2011 at 06:43 pm - Permalink

hugh50935 wrote:

I'm not entirely happy with my integration (i.e. the bit in the FOR loop). Is there a way to streamline this using wave assignments? I feel like I want to replace 't' with 'p', but I don't know how since the Variable U is a function of 't' as well.


Right. You won't be able to replace it with a wave assignment unless you directly substitute the expression for 'U' itself, which is certainly fine but will lead to a more cluttered expression.

So instead of this
 
U=(u0*sin(2*Pi*(t/n)))+(u1*sin(4*Pi*(t/n)))+(u2*sin(6*Pi*(t/n)))+(u3*sin(8*Pi*(t/n)))
wdF[t]=exp(2*Pi*i*(((h1+h2)*U)+(m*(t/n))))


you write
wdF[t]=exp(2*Pi*i*(((h1+h2)*(u0*sin(2*Pi*(t/n)))+(u1*sin(4*Pi*(t/n)))+(u2*sin(6*Pi*(t/n)))+(u3*sin(8*Pi*(t/n))))+(m*(t/n))))


Then just replace 't' with 'p' and you'll be all set. Sure, it's a mess, but the compiler doesn't care.

Strictly speaking using 'p' means that you would sample your function at unit intervals, but you're addressing that by having t/n as the effective variable. I'd recommend you to use wave scaling to achieve this. Simply execute
SetScale /P x, 0, 1/n, wdf

after making the wave, and then replace any occurrence of p/n in the wave assignment with x. That will clean it up some more. But in any case, make sure that you're sampling your function at a sufficiently high density to have a reliable result.

May 3, 2011 at 07:02 pm - Permalink

Ahh brilliant. Thanks for that.

Given my new increase in calculation speed from threading, increasing the sampling density of the integration shouldn't be an issue.

Now if I can get igor to email me when the calculation is finished on the server, life will be perfect.

Cheers,
Hugh

May 3, 2011 at 08:06 pm - Permalink

hugh50935 wrote:

Now if I can get igor to email me when the calculation is finished on the server, life will be perfect.


Well, you can get the Sockit XOP available from IgorExchange, and go to http://www.igorexchange.com/node/626.

But I think it wouldn't hurt if WaveMetrics made sending email a built-in operation. The reason I know of this snippet is because I'm also having Igor send me an email when a particular data-acquisition (that takes a few hours) is done.

May 3, 2011 at 09:01 pm - Permalink

If you have access to an SMS gateway (e.g. CLickatell), then you can send yourself an SMS when the thing has finished (using easyHttp). I do this for data acquisition and it's very useful.

May 3, 2011 at 10:33 pm - Permalink

You may also be able to use the area function instead of integrate, not sure which is faster.

A.

May 3, 2011 at 10:35 pm - Permalink

Unfortunately my backwards mobile provider doesn't support an SMS gateway. An email will suffice for now though.

I looked at SOCKIT but chickened out and wrote a script in C# which I execute at the end of the program.

I agree it would be great if Igor had some socket functions and, more importantly, a convenient email notification tool.

Cheers,
Hugh

May 3, 2011 at 11:19 pm - Permalink

hugh50935 wrote:
I looked at SOCKIT but chickened out and wrote a script in C# which I execute at the end of the program.


It's not as bad as it seems. Here's the function that I use to send email:
Function SendEmail(fromAddr, toAddr, subject, message)
    string fromAddr, toAddr, subject, message
   
    String SMTPServer = // your smtp server that can be reached on port 25
   
    variable sockID
    Make /T /O /FREE bufferWave
   
    variable nToAddresses = ItemsInList(toAddr, ","), i
   
    try
        SockitOpenConnection/q sockID, SMTPServer, 25, bufferWave
        if (sockID < 1)
            Abort
        endif
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, "EHLO <my domain>\n" // replace with your domain
        if (V_Flag != 0)
            Abort
        endif
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, "MAIL FROM:" + fromAddr + "\n"
        if (V_Flag != 0)
            Abort
        endif
       
        for (i = 0; i < nToAddresses; i+=1)
            SockitSendNRecv/SMAL/TIME=0.5 sockID, "RCPT TO:" + StringFromList(i, toAddr, ",") + "\n"
            if (V_Flag != 0)
                Abort
            endif
        endfor
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, "DATA\n"
        if (V_Flag != 0)
            Abort
        endif
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, "From:" + fromAddr + "\n"
        if (V_Flag != 0)
            Abort
        endif
       
        for (i = 0; i < nToAddresses; i+=1)
            SockitSendNRecv/SMAL/TIME=0.5 sockID, "To: " + StringFromList(i, toAddr, ",") + "\n"
            if (V_flag != 0)
                Abort
            endif
        endfor
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, "Subject:" + subject + "\n\n"
        if (V_Flag != 0)
            Abort
        endif
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, message
        if (V_Flag != 0)
            Abort
        endif
       
        SockitSendNRecv/SMAL/TIME=0.5 sockID, "\r\n.\r\n"
        if (V_Flag != 0)
            Abort
        endif
    catch
        Print "Failed to send email"
    endtry
    SockitCloseConnection(sockID)
End

May 4, 2011 at 06:31 am - Permalink