In a previous post, I described ways to create SAS/IML vectors that contain uniformly spaced values. The methods did not involve writing any loops.

This post describes how to perform a similar operation: creating evenly spaced values on a two-dimensional grid. The DATA step solution is simple, but an efficient SAS/IML solution requires thinking about the problem in terms of vectors and matrices.

The DATA Step Solution

A grid of values is useful in many contexts, such as scoring and visualizing a statistical model. In the DATA step, you can use loops to create a grid of values. For example, to generate a 3 x 2 grid of values on the rectangle [0,2] x [0,1] you can use the following program:

data grid;
keep x y;
/** set parameters for the general problem **/
XMin = 0; XMax = 2; Nx = 3;
YMin = 0; YMax = 1; Ny = 2;
 
dx = (XMax-XMin)/(Nx-1);
dy = (YMax-YMin)/(Ny-1);
do y = YMin to YMax by dy;
   do x = XMin to XMax by dx;
      output;
   end;
end;
run;
 
proc print; var x y; run;

This solution uses a double loop to compute the values on a uniform grid. There are Nx points in the X direction, evenly spaced on the interval [XMin, XMax], and Ny points in the Y direction, evenly spaced on the interval [YMin, YMax].

As I've written before, the SAS/IML language is a vector-matrix language, and an efficient SAS/IML program is one that avoids writing unnecessary loops. How can you write a similar program in PROC IML that avoids one or both loops?

A SAS/IML Language Solution

UPDATE 2014: In SAS/IML 12.3 (released with SAS 9.4), you can use the EXPANDGRID function to generate a grid in multiple dimensions! The following program is not needed if you have a recent version of SAS.

There are three simple SAS/IML functions that you can use to create a module that computes a grid of evenly spaced values:

• The DO function for generating sequential values, which is described in a previous post.
• The REPEAT function, which repeats a pattern a specified number of times.
• The SHAPE function, which reshapes a vector into a matrix. I also wrote a previous post on reshaping matrices.

The first step is to use the parameters in the problem to determine the points in the X and Y directions that will form the grid:

proc iml;
XMin = 0; XMax = 2; Nx = 3;
YMin = 0; YMax = 1; Ny = 2;
XDiv = do(XMin, XMax, (XMax-XMin)/(Nx-1));
YDiv = do(YMin, YMax, (YMax-YMin)/(Ny-1));
print XDiv, YDiv;

These are the coordinates of the grid points, but you need to arrange them into a matrix of ordered pairs that contains Nx * Ny rows and 2 columns.

The X coordinates are easy to construct. To form the X coordinates you can repeat the XDiv vector Ny times, and then reshape those values into a column vector:

X = repeat(XDiv, Ny);
X = shape(X, Nx*Ny);  /** or use COLVEC module in IMLMLIB **/
print X;

The Y coordinates are a little more difficult to construct. Recall that SAS/IML matrices are stored row-wise. This means that if you use the REPEAT function to form a set of points from the YDiv vector, the points will be in the wrong order:

Y = repeat(YDiv, Nx);
print Y; /** wrong order: elements stored row-wise **/

What you really want is to traverse the Y matrix down the columns instead of across the rows. No problem! Use the T function to transpose the values, and then reshape those transposed values into a column:

Y = shape(T(Y), Nx*Ny);  /** or use COLVEC module in IMLMLIB **/

A SAS/IML Module That Defines Ordered Pairs on a Grid

You can put these statements into a module, which enables you to conveniently generate a grid of points:

/** Module to define ordered pairs in a uniform grid of points.
   Input: XMin = minimum value of the X coordinate of the grid
          XMax = maximum value of the X coordinate of the grid
          Nx = number of grid lines along the X coordinate
          YMin = minimum value of the Y coordinate of the grid
          YMax = maximum value of the Y coordinate of the grid
          Ny = number of grid lines along the Y coordinate
 
    Output: Ordered pairs are returned in a matrix with
         (Nx x Ny) rows and two columns.
**/
start Define2DGrid(XMin, XMax, Nx, YMin, YMax, Ny);
   XDiv = do(XMin, XMax, (XMax-XMin)/(Nx-1));
   YDiv = do(YMin, YMax, (YMax-YMin)/(Ny-1));
   X = repeat(XDiv, Ny);
   X = shape(X, Nx*Ny);
   Y = repeat(YDiv, Nx);
   Y = shape(T(Y), Nx*Ny);
   return ( X || Y );
finish;
 
/** Example: Use the module to define a grid of points and 
      then evaluate a function f(x,y) on the grid. **/
p = Define2DGrid(-2, 2, 23,  /** 23 points in [-2,2] **/
                 -1, 1, 11); /** 11 points in [-1,1] **/
x = p[,1];  y = p[,2];
z = x##2 - 2*y##2;  /** evaluate a function z = f(x,y) **/
 
/** the following SAS/IML Studio statement creates a rotating plot **/
RotatingPlot.Create("Quadratic Surface", x, y, z);