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

[Heinz]

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