[Scilab-users] need a more efficient and faster code: suggestions welcome

classic Classic list List threaded Threaded
20 messages Options
Heinz Nabielek Heinz Nabielek
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] need a more efficient and faster code: suggestions welcome

Dear friends,

problem [below] is solved, but I need a more efficient and faster code: this is a Monte-Carlo simulation of n random points in a 3d spherical volume of radius r. A typical case is n=20,000 and r=23 mm. And from this, I need the resulting probability distribution of nearest neighbours. 

First attempt with doubly nested for-loops took 6 hours execution time. Version below with more vecctor operations is down to 6 minutes, but suggestions for further improvements would be welcome. 

Greetings

Heinz

 

n=20000;
r=23;

radius = r*grand(n,1,'def').^(1/3);

phi = 2*%pi*grand(n,1, 'def');

costheta = 1 - 2*grand(n,1, 'def');

radsintheta = radius.*sin(acos(costheta));

X = [radsintheta.*cos(phi),radsintheta.*sin(phi), radius.*costheta];

ONE=(ones(1:n-1))';

MinDist=[];

for i=1:n;

    this=X(i,:);XX=X;XX(i,:)=[];

    DIFF=XX-ONE*this;

    MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];

end;

jj=-0.05:0.1:3.05;nj=length(jj);

[0.05+jj(1:nj-1)' (histc(jj,MinDist))']

RESULTING OUTPUT
DIST.    PROBABILITY
[mm]
   0.    0.0003 
   0.1   0.0044 
   0.2   0.0202 
   0.3   0.0417 
   0.4   0.06735
   0.5   0.0984 
   0.6   0.12555
   0.7   0.13335
   0.8   0.1357 
   0.9   0.1155 
   1.    0.09735
   1.1   0.0684 
   1.2   0.0438 
   1.3   0.02495
   1.4   0.01385
   1.5   0.00535
   1.6   0.0023 
   1.7   0.00095
   1.8   0.0005 
   1.9   0.00005
   2.    0.00005
   2.1   0. 
   ...  ....


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

Hello,

The following suggestions will probably not have a drastic influence
(I don't see how it could be more vectorised)
but his a little thing I see:

> De : users [mailto:[hidden email]] De la part de Heinz Nabielek
> Envoyé : mercredi 31 janvier 2018 00:13
>
>    MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];

Maybe you could concatenate the squares of the distance
and then compute the square root of the whole vector in the end:

sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];



end



MinDist = sqrt(sqMinDist)

Hope this helps,

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

Thanks Christophe- I should have thought of it myself...
I have difficulties handling the SciLab timer: how do you do this exactly?
Heinz

PS: Had you noticed that the resulting probability had been predicted by Prof Paul Hertz in Heidelberg, Mathematische Annalen, Vol 67 (1909), 387ff as early as 20 December 1908? My 'real' materials problem is much more complex, but this part shown here would determine the time consuming part.

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Dang Ngoc Chan, Christophe
Sent: 31 January 2018 09:36
To: Users mailing list for Scilab <[hidden email]>
Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code: suggestions welcome

Hello,

The following suggestions will probably not have a drastic influence (I don't see how it could be more vectorised) but his a little thing I see:

> De : users [mailto:[hidden email]] De la part de Heinz
> Nabielek Envoyé : mercredi 31 janvier 2018 00:13
>
>    MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];

Maybe you could concatenate the squares of the distance and then compute the square root of the whole vector in the end:

sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];



end



MinDist = sqrt(sqMinDist)

Hope this helps,

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by Christophe Dang Ngoc Chan
Replacing

     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];

by

     MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];

will be at least twice faster. Crunching elapsed time could be done by
using parallel_run (with 5.5.2 version) if you have a multi-core processor.

S.

Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :

> Hello,
>
> The following suggestions will probably not have a drastic influence
> (I don't see how it could be more vectorised)
> but his a little thing I see:
>
>> De : users [mailto:[hidden email]] De la part de Heinz Nabielek
>> Envoyé : mercredi 31 janvier 2018 00:13
>>
>>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
> Maybe you could concatenate the squares of the distance
> and then compute the square root of the whole vector in the end:
>
> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
>
> …
>
> end
>
> …
>
> MinDist = sqrt(sqMinDist)
>
> Hope this helps,
>
> Regards
>
> --
> Christophe Dang Ngoc Chan
> Mechanical calculation engineer
> This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users


--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

moreover,

     MinDist=[MinDist sqrt(min((DIFF.*DIFF)*[1;1;1]))];

is even faster.

S.

Le 31/01/2018 à 10:53, Stéphane Mottelet a écrit :

> Replacing
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>
> by
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
>
> will be at least twice faster. Crunching elapsed time could be done by
> using parallel_run (with 5.5.2 version) if you have a multi-core
> processor.
>
> S.
>
> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
>> Hello,
>>
>> The following suggestions will probably not have a drastic influence
>> (I don't see how it could be more vectorised)
>> but his a little thing I see:
>>
>>> De : users [mailto:[hidden email]] De la part de
>>> Heinz Nabielek
>>> Envoyé : mercredi 31 janvier 2018 00:13
>>>
>>>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>> Maybe you could concatenate the squares of the distance
>> and then compute the square root of the whole vector in the end:
>>
>> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
>>
>> …
>>
>> end
>>
>> …
>>
>> MinDist = sqrt(sqMinDist)
>>
>> Hope this helps,
>>
>> Regards
>>
>> --
>> Christophe Dang Ngoc Chan
>> Mechanical calculation engineer
>> This e-mail may contain confidential and/or privileged information.
>> If you are not the intended recipient (or have received this e-mail
>> in error), please notify the sender immediately and destroy this
>> e-mail. Any unauthorized copying, disclosure or distribution of the
>> material in this e-mail is strictly forbidden.
>> _______________________________________________
>> users mailing list
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users
>
>

--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
amonmayr amonmayr
Reply | Threaded
Open this post in threaded view
|

Re: ?==?utf-8?q? {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by mottelet
Hello Stéphane,

Sorry to hijack the discussion but I didn't know that there was such a difference between A.*A and A.^2.
Could you tell us more about it?
Why is is twice faster to use the A.*A form?
Is this documented somewhere?

Cheers,

Antoine
 
 
Le Mercredi, Janvier 31, 2018 10:53 CET, Stéphane Mottelet <[hidden email]> a écrit:
 

> Replacing
>
>      MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>
> by
>
>      MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
>
> will be at least twice faster. Crunching elapsed time could be done by
> using parallel_run (with 5.5.2 version) if you have a multi-core processor.
>
> S.
>
> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
> > Hello,
> >
> > The following suggestions will probably not have a drastic influence
> > (I don't see how it could be more vectorised)
> > but his a little thing I see:
> >
> >> De : users [mailto:[hidden email]] De la part de Heinz Nabielek
> >> Envoyé : mercredi 31 janvier 2018 00:13
> >>
> >>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
> > Maybe you could concatenate the squares of the distance
> > and then compute the square root of the whole vector in the end:
> >
> > sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
> >
> > …
> >
> > end
> >
> > …
> >
> > MinDist = sqrt(sqMinDist)
> >
> > Hope this helps,
> >
> > Regards
> >
> > --
> > Christophe Dang Ngoc Chan
> > Mechanical calculation engineer
> > This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
> > _______________________________________________
> > users mailing list
> > [hidden email]
> > http://lists.scilab.org/mailman/listinfo/users
>
>
> --
> Stéphane Mottelet
> Ingénieur de recherche
> EA 4297 Transformations Intégrées de la Matière Renouvelable
> Département Génie des Procédés Industriels
> Sorbonne Universités - Université de Technologie de Compiègne
> CS 60319, 60203 Compiègne cedex
> Tel : +33(0)344234688
> http://www.utc.fr/~mottelet
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users
>

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: ?==?utf-8?q? {EXT} need a more efficient and faster code: suggestions welcome

I am an old school guy and have learned scientific computing with
Fortran... Using mutiplication instead of power elevation is an old
trick which should not be necessary with a more clever interpreter
(which should detect that 2 is actually a (small) integer and use
multiplication instead)

S.

Le 31/01/2018 à 11:22, Antoine Monmayrant a écrit :

> Hello Stéphane,
>
> Sorry to hijack the discussion but I didn't know that there was such a difference between A.*A and A.^2.
> Could you tell us more about it?
> Why is is twice faster to use the A.*A form?
> Is this documented somewhere?
>
> Cheers,
>
> Antoine
>  
>  
> Le Mercredi, Janvier 31, 2018 10:53 CET, Stéphane Mottelet <[hidden email]> a écrit:
>  
>> Replacing
>>
>>       MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>>
>> by
>>
>>       MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
>>
>> will be at least twice faster. Crunching elapsed time could be done by
>> using parallel_run (with 5.5.2 version) if you have a multi-core processor.
>>
>> S.
>>
>> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
>>> Hello,
>>>
>>> The following suggestions will probably not have a drastic influence
>>> (I don't see how it could be more vectorised)
>>> but his a little thing I see:
>>>
>>>> De : users [mailto:[hidden email]] De la part de Heinz Nabielek
>>>> Envoyé : mercredi 31 janvier 2018 00:13
>>>>
>>>>      MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>>> Maybe you could concatenate the squares of the distance
>>> and then compute the square root of the whole vector in the end:
>>>
>>> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
>>>
>>> …
>>>
>>> end
>>>
>>> …
>>>
>>> MinDist = sqrt(sqMinDist)
>>>
>>> Hope this helps,
>>>
>>> Regards
>>>
>>> --
>>> Christophe Dang Ngoc Chan
>>> Mechanical calculation engineer
>>> This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
>>> _______________________________________________
>>> users mailing list
>>> [hidden email]
>>> http://lists.scilab.org/mailman/listinfo/users
>>
>> --
>> Stéphane Mottelet
>> Ingénieur de recherche
>> EA 4297 Transformations Intégrées de la Matière Renouvelable
>> Département Génie des Procédés Industriels
>> Sorbonne Universités - Université de Technologie de Compiègne
>> CS 60319, 60203 Compiègne cedex
>> Tel : +33(0)344234688
>> http://www.utc.fr/~mottelet
>>
>> _______________________________________________
>> users mailing list
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users
>>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users


--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
amonmayr amonmayr
Reply | Threaded
Open this post in threaded view
|

Re: ?==?utf-8?q? ?==?utf-8?q? ?= {EXT} need a more efficient and faster code: suggestions welcom

Argh, OK, I get it: scilab treats A.^2 exactly like A.^%pi and not like "square it".
Makes sense, thank you for the info.
I wonder how julia is performing with respect to A.^2 compared to A.*A...

Antoine
 
 
Le Mercredi, Janvier 31, 2018 11:30 CET, Stéphane Mottelet <[hidden email]> a écrit:
 

> I am an old school guy and have learned scientific computing with
> Fortran... Using mutiplication instead of power elevation is an old
> trick which should not be necessary with a more clever interpreter
> (which should detect that 2 is actually a (small) integer and use
> multiplication instead)
>
> S.
>
> Le 31/01/2018 à 11:22, Antoine Monmayrant a écrit :
> > Hello Stéphane,
> >
> > Sorry to hijack the discussion but I didn't know that there was such a difference between A.*A and A.^2.
> > Could you tell us more about it?
> > Why is is twice faster to use the A.*A form?
> > Is this documented somewhere?
> >
> > Cheers,
> >
> > Antoine
> >  
> >  
> > Le Mercredi, Janvier 31, 2018 10:53 CET, Stéphane Mottelet <[hidden email]> a écrit:
> >  
> >> Replacing
> >>
> >>       MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
> >>
> >> by
> >>
> >>       MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
> >>
> >> will be at least twice faster. Crunching elapsed time could be done by
> >> using parallel_run (with 5.5.2 version) if you have a multi-core processor.
> >>
> >> S.
> >>
> >> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
> >>> Hello,
> >>>
> >>> The following suggestions will probably not have a drastic influence
> >>> (I don't see how it could be more vectorised)
> >>> but his a little thing I see:
> >>>
> >>>> De : users [mailto:[hidden email]] De la part de Heinz Nabielek
> >>>> Envoyé : mercredi 31 janvier 2018 00:13
> >>>>
> >>>>      MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
> >>> Maybe you could concatenate the squares of the distance
> >>> and then compute the square root of the whole vector in the end:
> >>>
> >>> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
> >>>
> >>> …
> >>>
> >>> end
> >>>
> >>> …
> >>>
> >>> MinDist = sqrt(sqMinDist)
> >>>
> >>> Hope this helps,
> >>>
> >>> Regards
> >>>
> >>> --
> >>> Christophe Dang Ngoc Chan
> >>> Mechanical calculation engineer
> >>> This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
> >>> _______________________________________________
> >>> users mailing list
> >>> [hidden email]
> >>> http://lists.scilab.org/mailman/listinfo/users
> >>
> >> --
> >> Stéphane Mottelet
> >> Ingénieur de recherche
> >> EA 4297 Transformations Intégrées de la Matière Renouvelable
> >> Département Génie des Procédés Industriels
> >> Sorbonne Universités - Université de Technologie de Compiègne
> >> CS 60319, 60203 Compiègne cedex
> >> Tel : +33(0)344234688
> >> http://www.utc.fr/~mottelet
> >>
> >> _______________________________________________
> >> users mailing list
> >> [hidden email]
> >> http://lists.scilab.org/mailman/listinfo/users
> >>
> > _______________________________________________
> > users mailing list
> > [hidden email]
> > http://lists.scilab.org/mailman/listinfo/users
>
>
> --
> Stéphane Mottelet
> Ingénieur de recherche
> EA 4297 Transformations Intégrées de la Matière Renouvelable
> Département Génie des Procédés Industriels
> Sorbonne Universités - Université de Technologie de Compiègne
> CS 60319, 60203 Compiègne cedex
> Tel : +33(0)344234688
> http://www.utc.fr/~mottelet
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users
>

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by mottelet
Dear Heinz,

Here is what can be done in Scilab-5.5.2 to accelerate your computations
on a multi-core architecture under Linux :

function out=distances(X)
     function out=distance(k)
         this=X(k,:);XX=X;XX(k,:)=[];
         DIFF=XX-ONE*this;
         out=sqrt(min((DIFF.*DIFF)*[1;1;1]));
     endfunction
     n=size(X,1);
     ONE=ones(n-1,1);
     out=parallel_run(1:n,distance);
endfunction

n=40000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1, 'def');
costheta = 1 - 2*grand(n,1, 'def');
radsintheta = radius.*sin(acos(costheta));
X = [radsintheta.*cos(phi),radsintheta.*sin(phi), radius.*costheta];

On a modest server (two 10-core/20-threads Intel(R) Xeon(R) CPU E5-2660
v2 @ 2.20GHz = 40 threads),

-->tic;distances(X);toc
  ans  =

     13.063

On a smarter server (sixteen 6-core/6-threads Intel(R) Xeon(R) CPU X7542
@ 2.67GHz = 96 threads)

-->tic;d=distances(X);toc
  ans  =

     2.498

I hope that parallel_run will be available again in the next version of
scilab.

S.

Le 31/01/2018 à 10:53, Stéphane Mottelet a écrit :

> Replacing
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>
> by
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
>
> will be at least twice faster. Crunching elapsed time could be done by
> using parallel_run (with 5.5.2 version) if you have a multi-core
> processor.
>
> S.
>
> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
>> Hello,
>>
>> The following suggestions will probably not have a drastic influence
>> (I don't see how it could be more vectorised)
>> but his a little thing I see:
>>
>>> De : users [mailto:[hidden email]] De la part de
>>> Heinz Nabielek
>>> Envoyé : mercredi 31 janvier 2018 00:13
>>>
>>>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>> Maybe you could concatenate the squares of the distance
>> and then compute the square root of the whole vector in the end:
>>
>> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
>>
>> …
>>
>> end
>>
>> …
>>
>> MinDist = sqrt(sqMinDist)
>>
>> Hope this helps,
>>
>> Regards
>>
>> --
>> Christophe Dang Ngoc Chan
>> Mechanical calculation engineer
>> This e-mail may contain confidential and/or privileged information.
>> If you are not the intended recipient (or have received this e-mail
>> in error), please notify the sender immediately and destroy this
>> e-mail. Any unauthorized copying, disclosure or distribution of the
>> material in this e-mail is strictly forbidden.
>> _______________________________________________
>> users mailing list
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users
>
>

--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

My latest tictoc is 366: would that be 6.1 minutes?
I am away from my 2 iMacs at home and running Scilab on a lowly Intel Pentium CPU N3540 @ 2.16GHz Win10 laptop. Your "2.498" for 40,000 points are phenomenal regarding the fact that the running time increases with the square of the number of points.

The later iMac has a Quad-Core-i7 at 3.2 GHz: can I run your code version there under the older Scilab? Not that I understand how parallel processing works...
Heinz

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Stéphane Mottelet
Sent: 31 January 2018 16:12
To: [hidden email]
Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code: suggestions welcome

Dear Heinz,

Here is what can be done in Scilab-5.5.2 to accelerate your computations on a multi-core architecture under Linux :

function out=distances(X)
     function out=distance(k)
         this=X(k,:);XX=X;XX(k,:)=[];
         DIFF=XX-ONE*this;
         out=sqrt(min((DIFF.*DIFF)*[1;1;1]));
     endfunction
     n=size(X,1);
     ONE=ones(n-1,1);
     out=parallel_run(1:n,distance);
endfunction

n=40000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1, 'def');
costheta = 1 - 2*grand(n,1, 'def');
radsintheta = radius.*sin(acos(costheta)); X = [radsintheta.*cos(phi),radsintheta.*sin(phi), radius.*costheta];

On a modest server (two 10-core/20-threads Intel(R) Xeon(R) CPU E5-2660
v2 @ 2.20GHz = 40 threads),

-->tic;distances(X);toc
  ans  =

     13.063

On a smarter server (sixteen 6-core/6-threads Intel(R) Xeon(R) CPU X7542 @ 2.67GHz = 96 threads)

-->tic;d=distances(X);toc
  ans  =

     2.498

I hope that parallel_run will be available again in the next version of scilab.

S.

Le 31/01/2018 à 10:53, Stéphane Mottelet a écrit :

> Replacing
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>
> by
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
>
> will be at least twice faster. Crunching elapsed time could be done by
> using parallel_run (with 5.5.2 version) if you have a multi-core
> processor.
>
> S.
>
> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
>> Hello,
>>
>> The following suggestions will probably not have a drastic influence
>> (I don't see how it could be more vectorised) but his a little thing
>> I see:
>>
>>> De : users [mailto:[hidden email]] De la part de
>>> Heinz Nabielek Envoyé : mercredi 31 janvier 2018 00:13
>>>
>>>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>> Maybe you could concatenate the squares of the distance and then
>> compute the square root of the whole vector in the end:
>>
>> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
>>
>> …
>>
>> end
>>
>> …
>>
>> MinDist = sqrt(sqMinDist)
>>
>> Hope this helps,
>>
>> Regards
>>
>> --
>> Christophe Dang Ngoc Chan
>> Mechanical calculation engineer
>> This e-mail may contain confidential and/or privileged information.
>> If you are not the intended recipient (or have received this e-mail
>> in error), please notify the sender immediately and destroy this
>> e-mail. Any unauthorized copying, disclosure or distribution of the
>> material in this e-mail is strictly forbidden.
>> _______________________________________________
>> users mailing list
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users
>
>

--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

Hello,

Don't even think about using parallel_run under OSX. It used to work (with some tweaking) with scilab-5.5.1 under Mavericks, but since 5.5.2 version, it has become completely unstable.

S.

Heinz <[hidden email]> a écrit :

My latest tictoc is 366: would that be 6.1 minutes?
I am away from my 2 iMacs at home and running Scilab on a lowly Intel Pentium CPU N3540 @ 2.16GHz Win10 laptop. Your "2.498" for 40,000 points are phenomenal regarding the fact that the running time increases with the square of the number of points.

The later iMac has a Quad-Core-i7 at 3.2 GHz: can I run your code version there under the older Scilab? Not that I understand how parallel processing works...
Heinz

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Stéphane Mottelet
Sent: 31 January 2018 16:12
To: [hidden email]
Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code: suggestions welcome

Dear Heinz,

Here is what can be done in Scilab-5.5.2 to accelerate your computations on a multi-core architecture under Linux :

function out=distances(X)
    function out=distance(k)
        this=X(k,:);XX=X;XX(k,:)=[];
        DIFF=XX-ONE*this;
        out=sqrt(min((DIFF.*DIFF)*[1;1;1]));
    endfunction
    n=size(X,1);
    ONE=ones(n-1,1);
    out=parallel_run(1:n,distance);
endfunction

n=40000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1, 'def');
costheta = 1 - 2*grand(n,1, 'def');
radsintheta = radius.*sin(acos(costheta)); X = [radsintheta.*cos(phi),radsintheta.*sin(phi), radius.*costheta];

On a modest server (two 10-core/20-threads Intel(R) Xeon(R) CPU E5-2660
v2 @ 2.20GHz = 40 threads),

-->tic;distances(X);toc
ans  =

    13.063

On a smarter server (sixteen 6-core/6-threads Intel(R) Xeon(R) CPU X7542 @ 2.67GHz = 96 threads)

-->tic;d=distances(X);toc
ans  =

    2.498

I hope that parallel_run will be available again in the next version of scilab.

S.

Le 31/01/2018 à 10:53, Stéphane Mottelet a écrit :

Replacing

    MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];

by

    MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];

will be at least twice faster. Crunching elapsed time could be done by
using parallel_run (with 5.5.2 version) if you have a multi-core
processor.

S.

Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :

Hello,

The following suggestions will probably not have a drastic influence
(I don't see how it could be more vectorised) but his a little thing
I see:

De : users [mailto:[hidden email]] De la part de
Heinz Nabielek Envoyé : mercredi 31 janvier 2018 00:13

    MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];

Maybe you could concatenate the squares of the distance and then
compute the square root of the whole vector in the end:

sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];



end



MinDist = sqrt(sqMinDist)

Hope this helps,

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information.
If you are not the intended recipient (or have received this e-mail
in error), please notify the sender immediately and destroy this
e-mail. Any unauthorized copying, disclosure or distribution of the
material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email].<a href="orghttp://lists.scilab.org/mailman/listinfo/users" target="_blank">orghttp://lists.scilab.org/mailman/listinfo/users




_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by mottelet
On my puny 250 Euro Win10 travel-laptop, I have achieved 66 tictocs for 20,000 random points with your code [not that I would understand it] in Scilab 5.5.2 and that is phenomenal and opens up new simulation possibilities.....

Thanks so much....
Heinz

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Stéphane Mottelet
Sent: 31 January 2018 16:12
To: [hidden email]
Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code: suggestions welcome

Dear Heinz,

Here is what can be done in Scilab-5.5.2 to accelerate your computations on a multi-core architecture under Linux :

function out=distances(X)
     function out=distance(k)
         this=X(k,:);XX=X;XX(k,:)=[];
         DIFF=XX-ONE*this;
         out=sqrt(min((DIFF.*DIFF)*[1;1;1]));
     endfunction
     n=size(X,1);
     ONE=ones(n-1,1);
     out=parallel_run(1:n,distance);
endfunction

n=40000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1, 'def');
costheta = 1 - 2*grand(n,1, 'def');
radsintheta = radius.*sin(acos(costheta)); X = [radsintheta.*cos(phi),radsintheta.*sin(phi), radius.*costheta];

On a modest server (two 10-core/20-threads Intel(R) Xeon(R) CPU E5-2660
v2 @ 2.20GHz = 40 threads),

-->tic;distances(X);toc
  ans  =

     13.063

On a smarter server (sixteen 6-core/6-threads Intel(R) Xeon(R) CPU X7542 @ 2.67GHz = 96 threads)

-->tic;d=distances(X);toc
  ans  =

     2.498

I hope that parallel_run will be available again in the next version of scilab.

S.

Le 31/01/2018 à 10:53, Stéphane Mottelet a écrit :

> Replacing
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>
> by
>
>     MinDist=[MinDist sqrt(min(sum(DIFF.*DIFF,2)))];
>
> will be at least twice faster. Crunching elapsed time could be done by
> using parallel_run (with 5.5.2 version) if you have a multi-core
> processor.
>
> S.
>
> Le 31/01/2018 à 09:36, Dang Ngoc Chan, Christophe a écrit :
>> Hello,
>>
>> The following suggestions will probably not have a drastic influence
>> (I don't see how it could be more vectorised) but his a little thing
>> I see:
>>
>>> De : users [mailto:[hidden email]] De la part de
>>> Heinz Nabielek Envoyé : mercredi 31 janvier 2018 00:13
>>>
>>>     MinDist=[MinDist sqrt(min(sum(DIFF.^2,2)))];
>> Maybe you could concatenate the squares of the distance and then
>> compute the square root of the whole vector in the end:
>>
>> sqMinDist=[sqMinDist min(sum(DIFF.^2,2))];
>>
>> …
>>
>> end
>>
>> …
>>
>> MinDist = sqrt(sqMinDist)
>>
>> Hope this helps,
>>
>> Regards
>>
>> --
>> Christophe Dang Ngoc Chan
>> Mechanical calculation engineer
>> This e-mail may contain confidential and/or privileged information.
>> If you are not the intended recipient (or have received this e-mail
>> in error), please notify the sender immediately and destroy this
>> e-mail. Any unauthorized copying, disclosure or distribution of the
>> material in this e-mail is strictly forbidden.
>> _______________________________________________
>> users mailing list
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users
>
>

--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

Hi Heinz,

Find herein a vectorised memory hog solution for Scilab 6, which is fast (~30 s for 20K points on Win7 laptop):

// START OF CODE (Scilab 6)
Clear
n=20000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1, 'def');
costheta = 1 - 2*grand(n,1, 'def');
radsintheta = radius.*sin(acos(costheta));
X = [radsintheta.*cos(phi), radsintheta.*sin(phi), radius.*costheta];
tic;
X2 = sum(X.*X, 'c');
X3 = X2.*.ones(n,1)';
X3 = X3 + X3';
D = abs(X3 - 2*(X*X'))';
D(1:n+1:$) = [];  // remove diagonal 0's
D = matrix(D,n-1,n); // reshape D to (n-1) x n size
MinDist1 = sqrt(min(D,'r'));
t1=toc();
disp(t1)
jj=-0.05:0.1:3.05;
H1 = histc(jj,MinDist1);
disp([0.05+jj(1:$-1)' H1']);
// END OF CODE

Regards,
Rafael

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Wozai Wozai
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

Hi,

On a Intel(R) Core(TM) i3-3220T CPU @ 2.80GHz (Linux), the following code
needs ~8s and
it doesn't allocate that much memory.

n=20000;
r=23;

radius = r*grand(n,1,'def').^(1/3);

phi = 2*%pi*grand(n,1, 'def');

costheta = 1 - 2*grand(n,1, 'def');

radsintheta = radius.*sin(acos(costheta));

X = [radsintheta.*cos(phi),radsintheta.*sin(phi), radius.*costheta];

ONE=ones(n,1);
NAN=[%nan,%nan,%nan];
SUM=[1;1;1];

MinDist=zeros(1,n);

for i=1:n;

    XX = X;
    XX(i,:) = NAN;
    DIFF = XX - ONE*X(i,:);
    MinDist(i) = sqrt(min((DIFF.*DIFF)*SUM));

end;

Cheers,

wozai




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by Rafael Guerra
Hi all,

I hope that the following will stop the (useless) competition : if speed
is a real bottleneck for this particular computation, using plain C is
the only option. The self-contained script (n=20000) given at the end of
this message gives the following timings on various hardware (only one
processor core is used) :

2,8 GHz Quad-Core Intel Xeon : 1.08 s (OSX)
2.2 GHz Xeon(R) CPU E5-2660 v2 @  : 1.84 s (Linux)

Don't use plain Scilab if speed is a crucial issue....

//  start of code
source=['#include <math.h>'
'#define SQR(x) ((x)*(x))'
'#define MIN(x, y) (((x) < (y)) ? (x) : (y))'
'void mindist(double d[],double x[],int *n) {'
'  int i,j,k;'
'  double dist;'
'  for (i=0;i<*n;i++) { '
'    d[i]=INFINITY;'
'    for (j=0;j<*n;j++) if (i!=j) { '
'      dist=0;'
'      for (k=0;k<3;k++) dist+=SQR(x[j*3+k]-x[i*3+k]);'
'      d[i]=MIN(d[i],dist);'
'    }'
'  }'
'}']
mputl(source,'mindist.c')
ilib_for_link('mindist','mindist.c',[],"c");
exec loader.sce

n=20000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1,'def');
costheta = 1 - 2*grand(n,1,'def');
radsintheta = radius.*sin(acos(costheta));

X=[radsintheta.*cos(phi) radsintheta.*sin(phi) radius.*costheta];

n=size(X,1);
tic;
dist=sqrt(call("mindist",X',2,"d",n,3,"i","out",[1,n],1,"d"));
disp(toc())
// end of code

Le 31/01/2018 à 23:31, Rafael Guerra a écrit :

> Hi Heinz,
>
> Find herein a vectorised memory hog solution for Scilab 6, which is fast (~30 s for 20K points on Win7 laptop):
>
> // START OF CODE (Scilab 6)
> Clear
> n=20000;
> r=23;
> radius = r*grand(n,1,'def').^(1/3);
> phi = 2*%pi*grand(n,1, 'def');
> costheta = 1 - 2*grand(n,1, 'def');
> radsintheta = radius.*sin(acos(costheta));
> X = [radsintheta.*cos(phi), radsintheta.*sin(phi), radius.*costheta];
> tic;
> X2 = sum(X.*X, 'c');
> X3 = X2.*.ones(n,1)';
> X3 = X3 + X3';
> D = abs(X3 - 2*(X*X'))';
> D(1:n+1:$) = [];  // remove diagonal 0's
> D = matrix(D,n-1,n); // reshape D to (n-1) x n size
> MinDist1 = sqrt(min(D,'r'));
> t1=toc();
> disp(t1)
> jj=-0.05:0.1:3.05;
> H1 = histc(jj,MinDist1);
> disp([0.05+jj(1:$-1)' H1']);
> // END OF CODE
>
> Regards,
> Rafael
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users


--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

For me it's fun to see thevarious attempts at improving speed. Thanks to Stephane for the x*x tip being better than x^2. 

Best regards
Claus 

On Feb 1, 2018 17:11, "Stéphane Mottelet" <[hidden email]> wrote:
Hi all,

I hope that the following will stop the (useless) competition : if speed is a real bottleneck for this particular computation, using plain C is the only option. The self-contained script (n=20000) given at the end of this message gives the following timings on various hardware (only one processor core is used) :

2,8 GHz Quad-Core Intel Xeon : 1.08 s (OSX)
2.2 GHz Xeon(R) CPU E5-2660 v2 @  : 1.84 s (Linux)

Don't use plain Scilab if speed is a crucial issue....

//  start of code
source=['#include <math.h>'
'#define SQR(x) ((x)*(x))'
'#define MIN(x, y) (((x) < (y)) ? (x) : (y))'
'void mindist(double d[],double x[],int *n) {'
'  int i,j,k;'
'  double dist;'
'  for (i=0;i<*n;i++) { '
'    d[i]=INFINITY;'
'    for (j=0;j<*n;j++) if (i!=j) { '
'      dist=0;'
'      for (k=0;k<3;k++) dist+=SQR(x[j*3+k]-x[i*3+k]);'
'      d[i]=MIN(d[i],dist);'
'    }'
'  }'
'}']
mputl(source,'mindist.c')
ilib_for_link('mindist','mindist.c',[],"c");
exec loader.sce

n=20000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1,'def');
costheta = 1 - 2*grand(n,1,'def');
radsintheta = radius.*sin(acos(costheta));

X=[radsintheta.*cos(phi) radsintheta.*sin(phi) radius.*costheta];

n=size(X,1);
tic;
dist=sqrt(call("mindist",X',2,"d",n,3,"i","out",[1,n],1,"d"));
disp(toc())
// end of code

Le 31/01/2018 à 23:31, Rafael Guerra a écrit :
Hi Heinz,

Find herein a vectorised memory hog solution for Scilab 6, which is fast (~30 s for 20K points on Win7 laptop):

// START OF CODE (Scilab 6)
Clear
n=20000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1, 'def');
costheta = 1 - 2*grand(n,1, 'def');
radsintheta = radius.*sin(acos(costheta));
X = [radsintheta.*cos(phi), radsintheta.*sin(phi), radius.*costheta];
tic;
X2 = sum(X.*X, 'c');
X3 = X2.*.ones(n,1)';
X3 = X3 + X3';
D = abs(X3 - 2*(X*X'))';
D(1:n+1:$) = [];  // remove diagonal 0's
D = matrix(D,n-1,n); // reshape D to (n-1) x n size
MinDist1 = sqrt(min(D,'r'));
t1=toc();
disp(t1)
jj=-0.05:0.1:3.05;
H1 = histc(jj,MinDist1);
disp([0.05+jj(1:$-1)' H1']);
// END OF CODE

Regards,
Rafael

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : <a href="tel:%2B33%280%29344234688" value="+33344234688" target="_blank">+33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by mottelet
"plain C is the only option": absolutely correct. As a matter of fact, I had programmed the problem in C in the first place. Idea was to probe the possibilty of Scilab, because that would make me more flexible with quick data handling and quick diagram sketches.

On my lowly Win10 laptop, the first rewrite from C to Scilab with doubly nested loop took 6h, change to vectorization of the inner loop 6min and further obvious smoothings 4min.

I am immensely greatful how to do the SciLab->C handover, all my previous attempts had failed and everything was extremely fast: tictoc=3.02

And the coding results in the correct solution as can be shown with
//  PLOT COMPARISON TO THEORY
dd=gsort(dist,'g','i');yy=(1:n)/(n+1);
scf;ii=0:.05:2.2;plot(ii,1-exp(-n*((ii/r).^3)),'r+');plot(dd,yy);
xtitle("Nearest Neighbour Distribution in 3d","distance (mm)","cumulative probability");
legend('prediction Prof Paul Hertz 1909','Monte-Carlo simulation of 20,000 points',pos=4);
//   END OF COMPARISON
I can send the diagram if somebody is interested.

Thanks a great lot, Stéphane...
Heinz

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Stéphane Mottelet
Sent: 01 February 2018 17:10
To: [hidden email]
Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code: suggestions welcome

Hi all,

I hope that the following will stop the (useless) competition : if speed is a real bottleneck for this particular computation, using plain C is the only option. The self-contained script (n=20000) given at the end of this message gives the following timings on various hardware (only one processor core is used) :

2,8 GHz Quad-Core Intel Xeon : 1.08 s (OSX)
2.2 GHz Xeon(R) CPU E5-2660 v2 @  : 1.84 s (Linux)

Don't use plain Scilab if speed is a crucial issue....

//  start of code
source=['#include <math.h>'
'#define SQR(x) ((x)*(x))'
'#define MIN(x, y) (((x) < (y)) ? (x) : (y))'
'void mindist(double d[],double x[],int *n) {'
'  int i,j,k;'
'  double dist;'
'  for (i=0;i<*n;i++) { '
'    d[i]=INFINITY;'
'    for (j=0;j<*n;j++) if (i!=j) { '
'      dist=0;'
'      for (k=0;k<3;k++) dist+=SQR(x[j*3+k]-x[i*3+k]);'
'      d[i]=MIN(d[i],dist);'
'    }'
'  }'
'}']
mputl(source,'mindist.c')
ilib_for_link('mindist','mindist.c',[],"c");
exec loader.sce

n=20000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1,'def');
costheta = 1 - 2*grand(n,1,'def');
radsintheta = radius.*sin(acos(costheta));

X=[radsintheta.*cos(phi) radsintheta.*sin(phi) radius.*costheta];

n=size(X,1);
tic;
dist=sqrt(call("mindist",X',2,"d",n,3,"i","out",[1,n],1,"d"));
disp(toc())
// end of code

Le 31/01/2018 à 23:31, Rafael Guerra a écrit :

> Hi Heinz,
>
> Find herein a vectorised memory hog solution for Scilab 6, which is fast (~30 s for 20K points on Win7 laptop):
>
> // START OF CODE (Scilab 6)
> Clear
> n=20000;
> r=23;
> radius = r*grand(n,1,'def').^(1/3);
> phi = 2*%pi*grand(n,1, 'def');
> costheta = 1 - 2*grand(n,1, 'def');
> radsintheta = radius.*sin(acos(costheta)); X = [radsintheta.*cos(phi),
> radsintheta.*sin(phi), radius.*costheta]; tic;
> X2 = sum(X.*X, 'c');
> X3 = X2.*.ones(n,1)';
> X3 = X3 + X3';
> D = abs(X3 - 2*(X*X'))';
> D(1:n+1:$) = [];  // remove diagonal 0's D = matrix(D,n-1,n); //
> reshape D to (n-1) x n size
> MinDist1 = sqrt(min(D,'r'));
> t1=toc();
> disp(t1)
> jj=-0.05:0.1:3.05;
> H1 = histc(jj,MinDist1);
> disp([0.05+jj(1:$-1)' H1']);
> // END OF CODE
>
> Regards,
> Rafael
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users


--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

In reply to this post by Wozai
Hello,

Heinz's interesting problem allows testing the Scilab limits and to learn
more about Scilab code optimization for very large datasets.
In the figure here below the execution times are plotted for the efficient
single loop code provided by Wozai and the modified vectorized code provided
further below.
<http://mailinglists.scilab.org/file/t495698/NN_problem_loop_and_vectorized.png>

The execution times seem to increase with N^2 or worse... which does not
look very good. Maybe faster algorithms exist but lower level languages seem
definitely required for very large number of points. Thanks Stephane for the
example in C.

Scilab 6.0.0 complains about arrays with more than 2147483647 elements and
memory allocation problems are observed in Win7 32GB laptop when doing
computations with smaller matrices, far from that limit. A dummy example:
--> n=4e4; A = ones(n,1)*ones(1,n) + ones(n,1)*ones(1,n);
    Can not allocate 12800.00 MB memory.
However, >21GB of RAM are available.
Are these matrix size constraints going to be lifted in future Scilab 6
releases?


// START OF CODE (Scilab 6)
clear;
n= 30000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1,'def');
costheta = 1 - 2*grand(n,1,'def');
radsintheta = radius.*sin(acos(costheta));
X = [radsintheta.*cos(phi), radsintheta.*sin(phi), radius.*costheta];
Xs = sum(X.*X,'c');
XX = Xs.*.ones(n,1)' + ones(n,1).*.Xs';  // faster than using transpose
XX = XX - 2*X*X';
XX(1:n+1:$) = %nan;  // remove diagonal 0's
MinDist1 = min(XX,'r').^0.5;   //min taken along rows is faster
clear  Xs  XX;
//END OF CODE

Thanks and regards,
Rafael




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

Nearest neighbour search of n points requires n.(n+1)/2 comparisons,
therefore execution times proportional to n^2.
Heinz

Same as number of matches between n soccer teams........

-----Original Message-----
From: users [mailto:[hidden email]] On Behalf Of Rafael
Guerra
Sent: 03 February 2018 12:42
To: [hidden email]
Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code:
suggestions welcome

Hello,

Heinz's interesting problem allows testing the Scilab limits and to learn
more about Scilab code optimization for very large datasets.
In the figure here below the execution times are plotted for the efficient
single loop code provided by Wozai and the modified vectorized code provided
further below.
<http://mailinglists.scilab.org/file/t495698/NN_problem_loop_and_vectorized.
png>

The execution times seem to increase with N^2 or worse... which does not
look very good. Maybe faster algorithms exist but lower level languages seem
definitely required for very large number of points. Thanks Stephane for the
example in C.

Scilab 6.0.0 complains about arrays with more than 2147483647 elements and
memory allocation problems are observed in Win7 32GB laptop when doing
computations with smaller matrices, far from that limit. A dummy example:
--> n=4e4; A = ones(n,1)*ones(1,n) + ones(n,1)*ones(1,n);
    Can not allocate 12800.00 MB memory.
However, >21GB of RAM are available.
Are these matrix size constraints going to be lifted in future Scilab 6
releases?


// START OF CODE (Scilab 6)
clear;
n= 30000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1,'def');
costheta = 1 - 2*grand(n,1,'def');
radsintheta = radius.*sin(acos(costheta)); X = [radsintheta.*cos(phi),
radsintheta.*sin(phi), radius.*costheta]; Xs = sum(X.*X,'c'); XX =
Xs.*.ones(n,1)' + ones(n,1).*.Xs';  // faster than using transpose XX = XX -
2*X*X';
XX(1:n+1:$) = %nan;  // remove diagonal 0's
MinDist1 = min(XX,'r').^0.5;   //min taken along rows is faster
clear  Xs  XX;
//END OF CODE

Thanks and regards,
Rafael




--
Sent from:
http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.
html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} need a more efficient and faster code: suggestions welcome

SORRY: n.(n-1)/2 but stll n^2

Sent from Heinz Nabielek

> On 03 Feb 2018, at 15:57, Heinz <[hidden email]> wrote:
>
> Nearest neighbour search of n points requires n.(n+1)/2 comparisons,
> therefore execution times proportional to n^2.
> Heinz
>
> Same as number of matches between n soccer teams........
>
> -----Original Message-----
> From: users [mailto:[hidden email]] On Behalf Of Rafael
> Guerra
> Sent: 03 February 2018 12:42
> To: [hidden email]
> Subject: Re: [Scilab-users] {EXT} need a more efficient and faster code:
> suggestions welcome
>
> Hello,
>
> Heinz's interesting problem allows testing the Scilab limits and to learn
> more about Scilab code optimization for very large datasets.
> In the figure here below the execution times are plotted for the efficient
> single loop code provided by Wozai and the modified vectorized code provided
> further below.
> <http://mailinglists.scilab.org/file/t495698/NN_problem_loop_and_vectorized.
> png>
>
> The execution times seem to increase with N^2 or worse... which does not
> look very good. Maybe faster algorithms exist but lower level languages seem
> definitely required for very large number of points. Thanks Stephane for the
> example in C.
>
> Scilab 6.0.0 complains about arrays with more than 2147483647 elements and
> memory allocation problems are observed in Win7 32GB laptop when doing
> computations with smaller matrices, far from that limit. A dummy example:
> --> n=4e4; A = ones(n,1)*ones(1,n) + ones(n,1)*ones(1,n);
>    Can not allocate 12800.00 MB memory.
> However, >21GB of RAM are available.
> Are these matrix size constraints going to be lifted in future Scilab 6
> releases?
>
>
> // START OF CODE (Scilab 6)
> clear;
> n= 30000;
> r=23;
> radius = r*grand(n,1,'def').^(1/3);
> phi = 2*%pi*grand(n,1,'def');
> costheta = 1 - 2*grand(n,1,'def');
> radsintheta = radius.*sin(acos(costheta)); X = [radsintheta.*cos(phi),
> radsintheta.*sin(phi), radius.*costheta]; Xs = sum(X.*X,'c'); XX =
> Xs.*.ones(n,1)' + ones(n,1).*.Xs';  // faster than using transpose XX = XX -
> 2*X*X';
> XX(1:n+1:$) = %nan;  // remove diagonal 0's
> MinDist1 = min(XX,'r').^0.5;   //min taken along rows is faster
> clear  Xs  XX;
> //END OF CODE
>
> Thanks and regards,
> Rafael
>
>
>
>
> --
> Sent from:
> http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.
> html
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users
>
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users