"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:[email protected]] On Behalf Of Stéphane
Mottelet
Sent: 01 February 2018 17:10
To: [email protected]
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
> [email protected]
> 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
[email protected]
http://lists.scilab.org/mailman/listinfo/users
_______________________________________________
users mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/users
