Not enough resources: Simulation and the calculation of pair-wise distances

[tesu]

This question is a follow-up from Patrick's answer on the pair-wise calculation: https://communities.sas.com/message/239952

My current working SAS code is the following:

%let datetime_start = %sysfunc(TIME()) ;
%put START TIME: %sysfunc(datetime(),datetime14.);

data original (drop = i);
   do i = 1 to 1000000;
      lon = -78 + ranuni(i);
   lat = 37 + ranuni(i+1);
   output;
   end;
run;
proc surveyselect data = original method = SRS rep = 1000
  sampsize = 5500 out = sampledata noprint;
  id _all_;
run;

%macro geodist;

%do i=1 %to 1000;

data loc&i; set sampledata; if replicate=&i; drop replicate; recid+1; run;
data loc&i; retain recid; set loc&i; run;
data geodist&i(drop=_:);
  set loc&i nobs=nobs;
  if _n_=1 then
    do;
      if 0 then set loc&i loc&i(keep=lat lon rename=(lat=hlat lon=hlon));
      declare hash h1(dataset:"loc&i(keep=recid lat lon rename=(lat=hlat lon=hlon))");
      _rc=h1.defineKey("recid");
      _rc=h1.defineData("hlat","hlon");
      _rc=h1.defineDone();
    end;
  do _i= (recid+1) to nobs;  
    _rc=h1.find(key:_i);
    geodist=geodist(lat, lon, hlat, hlon, "m");
    output;
  end;
run;
proc kde data=geodist&i;
   univar geodist/method=srot out=kout&i plots=NONE NOPRINT;
run;
data kout&i; set kout&i; keep density; run;

%end;

%mend geodist;

%geodist;

%put END TIME: %sysfunc(datetime(),datetime14.);
%put PROCESSING TIME:  %sysfunc(putn(%sysevalf(%sysfunc(TIME())-&datetime_start.),mmss.)) (mm:ss) ;

My task is (1) to create a date set with 1 million longitude-latitude observations, (2) to sample 5500 observations 1000 times (hence, 1000 data sets), (3) to calculate the pair-wise distances for every dataset, and (4) to draw kernel densities. This is a huge calculation. The SAS system crashed about at the 80-th iteration with complaining that it did not have sufficient resources.

My question is whether or not the code can be modified to increase the calculation performance by minimize the iteration process. If there is no room for improvement, I may have to go down to a low-level programming, such as Fortran, C, and so on.

Thank you very much!

[SASKiwi]

So what is the exact SAS error message and on what step within your program did it crash?

[PGStats]

What are you expecting to get in the original dataset :

data original (drop = i);

   do i = 1 to 1000000;

      lon = -78 + ranuni(i);

      lat = 37 + ranuni(i+1);

   output;

   end;

run;

using the iteration control variable i as a random number seed for RANUNI seems awkward. It is equivalent to coding:

      lon = -78 + ranuni(1);

      lat = 37 + ranuni(2);

is that what you wanted?

What is the purpose of this large simulation? Are you really interested in the distances between every points or only trying to locate neighboring clusters?

PG

[Patrick]

Using the code you've posted I've added a few minor optimizations to it. It still takes on my laptop around 15 sec per cycle and therefore would take more than 4 hours to run the 1000 cycles. So may-be first answering Pierre's questions is a good idea.

options fullstimer msglevel=i;

%let datetime_start = %sysfunc(TIME());

%put START TIME: %sysfunc(datetime(),datetime14.);

data original (drop = i);

  do i = 1 to 1000000;

    lon = -78 + ranuni(i);

    lat = 37 + ranuni(i+1);

    output;

  end;

run;

proc surveyselect data = original method = SRS rep = 1000

  sampsize = 5500 out = sampledata noprint;

  id _all_;

run;

proc sql;

  create index replicate on sampledata;

quit;

%macro geodist;

/*  %do i=1 %to 1000;*/

  %do i=1 %to 10;

    data loc&i(drop=replicate);

      set sampledata(where=(replicate=&i) idxname=replicate);

      recid=_n_;

    run;

    data geodist&i (keep=geodist);

      set loc&i nobs=nobs;

      if _n_=1 then

        do;

          if 0 then

            set loc&i loc&i(keep=lat lon rename=(lat=hlat lon=hlon));

          declare hash h1(dataset:"loc&i(keep=recid lat lon rename=(lat=hlat lon=hlon))");

          _rc=h1.defineKey("recid");

          _rc=h1.defineData("hlat","hlon");

          _rc=h1.defineDone();

        end;

      _start=(recid+1);

      do _i= _start to nobs;

        _rc=h1.find(key:_i);

        geodist=geodist(lat, lon, hlat, hlon, "m");

        output;

      end;

    run;

    proc kde data=geodist&i;

      univar geodist/method=srot out=kout&i(keep=density) plots=NONE NOPRINT;

    run;

    proc datasets lib=work nolist nowarn;

      delete loc&i geodist&i;

      run;

    quit;

  %end;

%mend geodist;

%geodist;

%put END TIME: %sysfunc(datetime(),datetime14.);

%put PROCESSING TIME:  %sysfunc(putn(%sysevalf(%sysfunc(TIME())-&datetime_start.),mmss.)) (mm:ss);

[tesu]

Thanks again, Patrick. The 4-hour running time is quite long to me. I should use some C routines to do the work. -- Your hash table solution is a killer!

[Patrick]

I'd probably would run a 4 hour job in batch mode (writing the kout&i data sets to a permanent location).

If it's a one-off then I'd accept the 4 hours, if it's something you have to run more often then there would be other options to improve performance (eg. by running stuff in parallel). The end-to-end time is eventually much better in your environment with faster hardware.

The hash lookup could become faster if you set the a better value for hashexp (5 or 6 I'd guess; hashexp is a property of the hash defined during hash declaration). If not running stuff in parallel then doing everything in one go (with by group processing for Proc KDE) would likely improve performance as well.

[tesu]

Thanks for further guidance, Patrick. In order to use the "by" command in the KDE procedure, I believe I need to make one big "geodist" data set with an additional group variable, which goes from 1 to the number of rep=1000. I don't know how to modify the hash process for the big geodist data, which should have the (5500*(5500-1)/2)*1000 = 1.512e+10 observations. Here, (5500*(5500-1)/2) is the number of pair-wise distances for one geodist data set. I need to study the hash operation, which is new to me.

The PROCESSING TIME for [%do i=1 %to 10;] ended up with 3:48 (mm:ss).

The calculation of the pair-wise distances (15,122,250 observations) for one location sample takes up about 20 seconds. I believe I should write this tedious distance calculation job in C, and make a dll. Then, I may use it in the data step. Why I stick to SAs is because the KDE procedure is convenient. Loading a dll is a separate issue, and I may be visit again this community to ask further questions.

[gergely_batho]

Small improovement by using array, view and by group processing:

options fullstimer msglevel=i;

%let datetime_start = %sysfunc(TIME());

%put START TIME: %sysfunc(datetime(),datetime14.);

data original (drop = i);

  do i = 1 to 1000000;

    lon = -78 + ranuni(i);

    lat = 37 + ranuni(i+1);

    output;

  end;

run;

%let numrep=10;

%let sampsize=5500;

proc surveyselect data = original method = SRS rep = &numrep.

  sampsize = &sampsize. out = sampledata noprint;

  id _all_;

run;

options vbufsize=500M;

data geodist (keep=Replicate geodist) / view=geodist;

  array alat[&sampsize.] _temporary_;

  array alon[&sampsize.] _temporary_;

  do i=1 to &sampsize.;

  set sampledata;

  alat=lat;alon=lon;

  end;

      do i=1 to &sampsize.;

    lat=alat;lon=alon;

    do j=i+1 to &sampsize.;

         geodist=geodist(lat, lon, alat, alon, "m");

        output;

  end;

      end;

    run;

proc kde data=geodist(spill=no);

      univar geodist/method=srot out=kout(keep=density) plots=NONE NOPRINT;

   by Replicate;

run;

%put END TIME: %sysfunc(datetime(),datetime14.);

%put PROCESSING TIME:  %sysfunc(putn(%sysevalf(%sysfunc(TIME())-&datetime_start.),mmss.)) (mm:ss);

[Patrick · Accepted solution]

Below one way of how to run the processes in parallel. You will have to amend the code (the paths) for the code to run in your environment. Below code executed within less than 90 minutes on my Quad Core laptop with a elecromechanical hard drive. If you have more cores and a faster disk then end-to-end elapsed time could be quite a bit better.

If you run this code on a server used also by other people then make sure it executes outside of business hours and/or doesn't use up all the available resources (eg. by also using statement WAITFOR) - else you will become your companies most liked person....

ods listing;

options fullstimer msglevel=i;

%let rep=1000;

%let datetime_start = %sysfunc(TIME());

%put START TIME: %sysfunc(datetime(),datetime14.);

libname test 'c:\test';

data test.original (drop = i);

  do i = 1 to 1000000;

    lon = -78 + ranuni(i);

    lat = 37 + ranuni(i+1);

    output;

  end;

run;

proc surveyselect data = test.original method = SRS rep = &rep

  sampsize = 5500 out = test.sampledata noprint;

  id _all_;

run;

proc sql;

  create index replicate on test.sampledata;

quit;

%macro geodist;

  %do i=1 %to &rep;

    filename prg "c:\test\prog&i..sas";

    data _null_;

      file prg;

      put

        "libname test 'c:\test';"//

        " data loc&i(drop=replicate);"/

        "   set test.sampledata(where=(replicate=&i) idxname=replicate);"/

        "   recid=_n_;"/

        " run;"/

        " data _null_;"/

        "   call symputx('exp',ceil(log2(nobs)));"/

        "   stop;"/

        "   set loc&i nobs=nobs;"/

        " run;"/

        /

        " data geodist&i (keep=geodist);"/

        "   set loc&i nobs=nobs;"/

        /

        "   if _n_=1 then"/

        "     do;"/

        "       if 0 then"/

        "         set loc&i loc&i(keep=lat lon rename=(lat=hlat lon=hlon));"/

        "       declare hash h1(dataset:'loc&i(keep=recid lat lon rename=(lat=hlat lon=hlon))', " 'hashexp:&exp);'/

        "       _rc=h1.defineKey('recid');"/

        "       _rc=h1.defineData('hlat','hlon');"/

        "       _rc=h1.defineDone();"/

        "     end;"/

        /

        "   _start=(recid+1);"/

        "   do _i= _start to nobs;"/

        "     _rc=h1.find(key:_i);"/

        "     geodist=geodist(lat, lon, hlat, hlon, 'm');"/

        "     output;"/

        "   end;"/

        " run;"/

        /

        " proc kde data=geodist&i;"/

        "   univar geodist/method=srot out=test.kout&i(keep=density) plots=NONE NOPRINT;"/

        " run;"/

        ;

    run;

    %let sasbat=C:\SAS\Config\Lev1\SASApp\sas.bat;

    %let pgrm=prog&i..sas;

    %let test=c:\test;

    %let DTstamp=%sysfunc(datetime(), B8601DT.);

    %let cmd=&sasbat -sysin "&test\&pgrm";

    %let cmd=&cmd -log "&test\%scan(&pgrm,1,.)_&DTstamp..log";

    %let cmd=&cmd -print "&test\%scan(&pgrm,1,.)_&DTstamp..lst";

/*    %let cmd=&cmd -Metauser "&Metauser" -Metapass "&Metapass";*/

    %let cmd=&cmd -ICON -NOSPLASH;

    /* batch submit current program */

    systask command %unquote(%str(%')&cmd%str(%'))

      status    =S012345_&DTstamp

      mname     =taskname 

      nowait

      ;

    data _null_;

      call sleep(1,1);

      stop;

    run;

  %end;

%mend geodist;

%geodist;

%put END TIME: %sysfunc(datetime(),datetime14.);

%put PROCESSING TIME:  %sysfunc(putn(%sysevalf(%sysfunc(TIME())-&datetime_start.),mmss.)) (mm:ss);

[tesu]

Thank you very much, Patrick. The running time was around 18~19 minutes on my laptop (i3-3120M, 16GB DDR3 RAM, and a SSD). Great improvement. I will try hard to study the parallel processing that you provide, but for now I have just one question. After the job had been finished, I can't find any "kout&i" dataset, which is the output dataset from the KDE procedure. I need to collect those datasets to do something else. Each kout&i dataset has just one column the variable of which is density. I need to make a big matrix that has the 1000 density&i variables. -- Lots to learn in this world, and life is great!

Miscellaneous: Why do we need two dots as in "c:\test\prog&i..sas" and two slashes as in "libname test 'c:\test';"// ?

[Patrick]

Hmm... 18-19 minutes but then you miss the output. Looks to me rather that it didn't work for you. There is the following code bit in the program:

    data _null_;

      call sleep(1,1);

      stop;

    run;

What this does is waiting 1 second after every single iteration. So submitting 1000 jobs takes already more than 16 minutes. I would assume that the additional time spend is for the pre-steps and the program waiting to get a slot for another process.

I assume you haven't amended all the paths in the previously posted program for it to work in your environment. I've moved now this environment specific stuff at the beginning of the code. Make sure you pass values/paths fit for your environment.

When a macro variable resolves and it's ended with a dot, the dot gets consumed. So in order to get to something like "c:\test\prog10.sas" I need to have a second dot to feed to the macro variable &i "c:\test\prog&i..sas"

That's part of the PUT statement. The code below

filename prg "&dir_path\prog10.sas";

data _null_;

     file prg;

     PUT "libname test 'c:\test';"//

     ;

run;

Would write to file &dir_path\prog10.sas a string libname test 'c:\test';

This is a way of dynamically generating SAS code. You then can %include this code later on or as I do use "systask" and execute the generated code as a child process in batch.

The 2 slashes belong to the PUT statement and simply add 2 line breaks so the generated code will look nicer.

And here the newest code version

######################################

/* Directory for output (generated SAS code, logs, data) */

%let dir_path=c:\test;

/* Libref for result data sets */

%let outlibref=test;

/* number of iterations */

%let rep=2;

/*%let rep=1000;*/

/***** no changes below required ******/

options fullstimer msglevel=i;

%let datetime_start = %sysfunc(TIME());

%put START TIME: %sysfunc(datetime(),datetime14.);

libname &outlibref "&dir_path";

data &outlibref..original (drop = i);

  do i = 1 to 1000000;

    lon = -78 + ranuni(i);

    lat = 37 + ranuni(i+1);

    output;

  end;

run;

proc surveyselect data = &outlibref..original method = SRS rep = &rep

  sampsize = 5500 out = &outlibref..sampledata noprint;

  id _all_;

run;

proc sql;

  create index replicate on &outlibref..sampledata;

quit;

%macro geodist;

  %do i=1 %to &rep;

    filename prg "&dir_path\prog&i..sas";

    data _null_;

      file prg;

      put

        "libname &outlibref '&dir_path';"//

        " data loc&i(drop=replicate);"/

        "   set &outlibref..sampledata(where=(replicate=&i) idxname=replicate);"/

        "   recid=_n_;"/

        " run;"/

        " data _null_;"/

        "   call symputx('exp',ceil(log2(nobs)));"/

        "   stop;"/

        "   set loc&i nobs=nobs;"/

        " run;"/

        /

        " data geodist&i (keep=geodist);"/

        "   set loc&i nobs=nobs;"/

        /

        "   if _n_=1 then"/

        "     do;"/

        "       if 0 then"/

        "         set loc&i loc&i(keep=lat lon rename=(lat=hlat lon=hlon));"/

        "       declare hash h1(dataset:'loc&i(keep=recid lat lon rename=(lat=hlat lon=hlon))', " 'hashexp:&exp);'/

        "       _rc=h1.defineKey('recid');"/

        "       _rc=h1.defineData('hlat','hlon');"/

        "       _rc=h1.defineDone();"/

        "     end;"/

        /

        "   _start=(recid+1);"/

        "   do _i= _start to nobs;"/

        "     _rc=h1.find(key:_i);"/

        "     geodist=geodist(lat, lon, hlat, hlon, 'm');"/

        "     output;"/

        "   end;"/

        " run;"/

        /

        " proc kde data=geodist&i;"/

        "   univar geodist/method=srot out=&outlibref..kout&i(keep=density) plots=NONE NOPRINT;"/

        " run;"/

        ;

    run;

    %let sasbat=C:\SAS\Config\Lev1\SASApp\sas.bat;

    %let pgrm=prog&i..sas;

    %let DTstamp=%sysfunc(datetime(), B8601DT.);

    %let cmd=&sasbat -sysin "&dir_path\&pgrm";

    %let cmd=&cmd -log "&dir_path\%scan(&pgrm,1,.)_&DTstamp..log";

/*    %let cmd=&cmd -log "&dir_path\%scan(&pgrm,1,.)_#Y#m#d#H#M#s_#p.log" -logparm "rollover=auto";*/

    %let cmd=&cmd -print "&dir_path\%scan(&pgrm,1,.)_&DTstamp..lst";

/*    %let cmd=&cmd -Metauser "&Metauser" -Metapass "&Metapass";*/

    %let cmd=&cmd -ICON -NOSPLASH;

    /* batch submit current program */

    systask command %unquote(%str(%')&cmd%str(%'))

      status    =S012345_&DTstamp

      mname     =taskname

      nowait

      ;

    data _null_;

      call sleep(1,1);

      stop;

    run;

  %end;

%mend geodist;

%geodist;

%put END TIME: %sysfunc(datetime(),datetime14.);

%put PROCESSING TIME:  %sysfunc(putn(%sysevalf(%sysfunc(TIME())-&datetime_start.),mmss.)) (mm:ss);

[Patrick]

Now to your questions of how to combine the KOUT tables. Below code should give you what your after. You can copy this code at the end of the first program (without the definitions for &dir_path and &outlibref) or you can run it separately after the first program finished.

Make sure to only run it after all child processes finished. As systask uses "nowait" the last few child processes will still be running even though you've got control back to your parent session. So when running as part of the same program make sure to check what options WAIT & SYNC and statement WAITFOR can do for you.

/* Directory for output tables (KOUTnnn tables) */

%let dir_path=c:\test;

/* Libref for kout data sets */

%let outlibref=test;

libname &outlibref "&dir_path";

proc sql noprint;

  create table kout_tables as

    select cats(libname,'.',memname,'(keep=density rename=(density=density',compress(memname,,'kd'),'))') as string

    from dictionary.tables

    where libname="%upcase(&outlibref)" and prxmatch('/kout\d+ *$/oi',memname)>0

    ;

quit;

filename codegen temp;

data _null_;

/*  file print;*/

  file codegen;

  set kout_tables end=last;

  put string;

  if last then put ';';

run;

data &outlibref..kout_all;

  merge

    %include codegen;

run;

[tesu]

EDIT: Your code works after I changed the %let sasbat statement to %let sasbat="C:\Program Files\SASHome\SASFoundation\9.3\sas.exe";

=== [ Work history ] ===

Thank you very much again for the update, Patrick. Same as before, the code did not generate any output. I DID NOT run the second section of your SAS code that combines the kout datasets. I include the log message and how the working folder and the test library look like:

I changed the working folder from c:\test to d:\patrick.

I set the number of task iteration as two: %let rep=2;

========================================

NOTE: The file PRG is:
      Filename=d:\patrick\prog1.sas,
      RECFM=V,LRECL=256,File Size (bytes)=0,
      Last Modified=02Dec2014:15:03:10,
      Create Time=02Dec2014:14:42:40

NOTE: 37 records were written to the file PRG.
      The minimum record length was 0.
      The maximum record length was 100.
NOTE: DATA statement used (Total process time):
      real time           0.00 seconds
      user cpu time       0.00 seconds
      system cpu time     0.00 seconds
      memory              277.93k
      OS Memory           13252.00k
      Timestamp           12/02/2014 03:03:10 PM

NOTE: Task "task0" produced no LOG/Output.
NOTE: DATA statement used (Total process time):
      real time           1.00 seconds
      user cpu time       0.01 seconds
      system cpu time     0.01 seconds
      memory              232.34k
      OS Memory           13252.00k
      Timestamp           12/02/2014 03:03:11 PM

NOTE: The file PRG is:
      Filename=d:\patrick\prog2.sas,
      RECFM=V,LRECL=256,File Size (bytes)=0,
      Last Modified=02Dec2014:15:03:11,
      Create Time=02Dec2014:14:42:41

NOTE: 37 records were written to the file PRG.
      The minimum record length was 0.
      The maximum record length was 100.
NOTE: DATA statement used (Total process time):
      real time           0.00 seconds
      user cpu time       0.00 seconds
      system cpu time     0.00 seconds
      memory              231.50k
      OS Memory           13252.00k
      Timestamp           12/02/2014 03:03:11 PM

NOTE: Task "task1" produced no LOG/Output.
NOTE: DATA statement used (Total process time):
      real time           1.00 seconds
      user cpu time       0.00 seconds
      system cpu time     0.00 seconds
      memory              232.34k
      OS Memory           13252.00k
      Timestamp           12/02/2014 03:03:12 PM

=======================================

Something goes wrong in the batch mode processing? The main LOG says "NOTE: Task "task0" produced no LOG/Output." The batch processing is completely a new area to me, but I will study those material. I am googling cases that SAS batch processing produces no outputs.

Thank you very much again for all of your responses.

EDIT:

I am reading this blog: SAS Programming for Data Mining: Parallel Processing with Single-Thread SAS license, which seems to be closely relevant to my problem. It seems that systask command operation is the product of SAS/CONNECT? My SAS license DOES contain that product.

EDIT:

From your comment in this thread: https://communities.sas.com/message/220782, I changed your %let sasbat statement to %let sasbat="C:\Program Files\SASHome\SASFoundation\9.3\sas.exe";. Now, it worked. I have the kout datasets in the test library, and in the d:\patrick.

[Patrick]

Forgot to move one thing to the top of the code for you to amend:  %let sasbat=C:\SAS\Config\Lev1\SASApp\sas.bat;

This is the path to the SAS executable (sas.exe) or a .bat file which then calls internally the SAS executable. This must point to the correct location in your environment.

The article you've posted shows ways of how to run stuff in parallel without SAS/Connect licensed. This includes the usage of SYSTASK. Just make sure that your environment doesn't have option NOXCMD set. It must be XCMD.

[tesu]

We are talking to each other real time right now. Yes. I did it. And it worked. Please look at my edited previous reply. I will get back to you with the total computation time for the 1,000 iterations.