Most Yorick statements look like algebraic formulas. A variable name is a string like Var_1 -- upper or lower case characters (case matters), digits, or underscores in any combination except that the first character may not be a digit. Expressions consist of the usual arithmetic operations + - * /, with parentheses to indicate the order of operations (when that order is different than or unclear from the ordinary rules of precedence in algebra). Elementary mathematical functions such as exp(x), cos(x), or atan(x) look just like that.
Usually, a Yorick variable is a parametric representation of a mathematical function. The variable is an array of numbers which are values of the function at a number of points; few points to represent the function coarsely, more for an accurate rendition. The parameters of the function are the indices into the array, which rarely make an explicit appearance in Yorick programs. Thus,
theta= span(0.0, 2*pi, 100)
defines a variable theta consisting of 100 evenly spaced values starting with 0.0 and ending with 2*pi.
Now that theta has been defined as a list of 100 numbers, any function of theta has a concrete representation as a list of 100 numbers -- namely the values of the function at the 100 particular values of theta. Hence, variables x and y representing coordinates of the unit circle are defined with:
x= cos(theta); y= sin(theta)
Here, cos and sin are built-in Yorick functions. Like most Yorick functions, they operate on an entire array of numbers, returning an array of like shape. Hence both x and y are now lists of 100 numbers -- the cosines and sines of the 100 numbers theta.
The semicolon marks the end of a Yorick statement, allowing several statements to share a single line. The end of a line (i.e.- a newline) can also mark the end of a Yorick statement. However, if any parentheses are open, or if a binary operator or a comma is the last token on the line, then the newline is treated like a space or a tab character and does not terminate the Yorick statement.
If a line ends with backslash, the following newline will never terminate the Yorick statement. (That is, backslash is the continuation character in Yorick.) I recommend that you never use a backslash -- end the line to be continued with a binary operator, or leave the comma separating subroutine arguments at the end of the line, or split a parenthetic expression across the line, and it will be continued automatically.
Including the span function introduced in the previous section, there are five common ways to originate arrays -- that is, to make an array out of scalar values:
As I have said, a Yorick array often represents the values of a continuous function at a number of discrete points. In order to find the values of the function at other points, you need to know how it varies between (or beyond) the given points. In general, to interpolate (or extrapolate), you need a detailed understanding of how the function was discretized in the first place. However, when the list of values accurately represents the function, linear interpolation between the known points will suffice. A function which is linear between successive points is called "piecewise linear".
The interp function is a mechanism for converting a list of function values at discrete points into a piecewise linear function which can be evaluated at any point.
theta= span(0, pi, 100); x_circle= cos(theta); y_circle= sin(theta); x= span(-2, 2, 64); y= interp(y_circle, x_circle, x);
This code fragment produces a y array with the same number of points as x (64), with the values of the piecewise linear function defined by the points (x_circle, y_circle). Outside the range covered by x_circle, the piecewise linear function remains constant -- the simplest possible extrapolation rule.
Regarded as a function of its third argument, interp behaves just like the sin or cos function -- its first two arguments are really parameters specifying which piecewise linear function interp will evaluate.
The integ function works just like interp, except that it returns the integral of the piecewise linear function. The integration constant is chosen so that integ returns zero at the first point of the piecewise linear function. (This point will actually have the maximum value of x if the x array is decreasing.) Thus, the integral of the piecewise linear approximation to the semicircle and the exact integral of the semicircle can be computed by:
yi= integ(y_circle, x_circle, x); yi_exact= 0.5*(acos(max(min(x,1),-1)) - x*sqrt(1-min(x^2,1)));
Again, the piecewise linear function is assumed to remain constant beyond the first and last points specified. Hence, integ is a linear function when extrapolating, and piecewise parabolic when interpolating.
Use integ only when you need the indefinite integral of a piecewise linear function. Yorick has more efficient ways to compute definite integrals. Again, think of integ, like interp, as a continuous function of its third argument; the first two arguments are parameters specifying which function.
Neither interp nor integ makes sense unless its second argument is either increasing or decreasing. There is no way to decide which branch of a multi-valued function should be returned.
Internally, both interp and integ need a lookup function -- that is, a function which finds the index of the point in x_circle just beyond each of the x values. This lookup function can also be called directly; its name is digitize.
Yorick has a bewildering variety of different ways to refer to individual array elements or subsets of array elements. In order to master the language, you must learn to use them all. Nearly all of the examples later in this manual use one or more of these indexing techniques, so trust me to show you how to use them later:
An array of objects is stored in consecutive locations in memory (where each location is big enough to hold one of the objects). An array x of three numbers is stored in the order [x(1), x(2), x(3)] in three consecutive slots in memory. A three-by-two array y means nothing more than an array of two arrays of three numbers each. Thus, the six numbers are stored in two contiguous blocks of three numbers each: [[x(1,1), x(2,1), x(3,1)], [x(1,2), x(2,2), x(3,2)]].
A multi-dimensional array may be referenced using fewer indices than its number of dimensions. Hence, in the previous example, x(5) is the same as x(2,2), since the latter element is stored fifth.
Although most of Yorick's syntax follows the C language, array indexing is designed to resemble FORTRAN array indexing. In Yorick, as in FORTRAN, the first (leftmost) dimension of an array is always the index which varies fastest in memory. Furthermore, the first element along any dimension is at index 1, so that a dimension of length three can be referenced by index 1 (the first element), index 2 (the second element), or index 3 (the third element).
If this inconsistency bothers you, here is why Yorick indexing is like FORTRAN indexing: In C, an array of three numbers, for example, is a data type on the same footing as the data type of each of its three members; by this trick C sidesteps the issue of multi-dimensional arrays --- they are singly arrays of objects of an array data type. While this picture accurately reflects the way the multi-dimensional array is stored in memory, it does not reflect the way a multi-dimensional array is used in a scientific computer program.
In such a program, the fact that the array is stored with one or the other index varying fastest is irrelevent -- you are equally likely to want to consider as a "data type" a slice at a constant value of the first dimension as of the second. Furthermore, the length of every dimension varies as you vary the resolution of the calculation in the corresponding physical direction.
You can refer to several consecutive array elements by an index range: x(3:6) means the four element subarray [x(3), x(4), x(5), x(6)].
Occasionally, you also want to refer to a sparse subset of an array; you can add an increment to an index range by means of a second colon: x(3:7:2) means the three element subarray [x(3), x(5), x(7)].
A negative increment reverses the order of the elements: x(7:3:-2) represents the same three elements, but in the opposite order [x(7), x(5), x(3)]. The second element mentioned in the index range may not actually be present in the resulting subset, for example, x(7:2:-2) is the same as x(7:3:-2), and x(3:6:2) represents the two element subarray [x(3), x(5)].
Just as the increment defaults to 1 if it is omitted, the start and stop elements of an index range also have default values, namely the first and last possible index values. Hence, if x is a one-dimensional array with 10 elements, x(8:) is the same as x(8:10). With a negative increment, the defaults are reversed, so that x(:8:-1) is the same as x(10:8:-1).
A useful special case of the index range default rules is x(::-1), which represents the array x in reverse order.
Beware of a minor subtlety: x(3:3) is not the same thing as x(3). An index range always represents an array of values, while a scalar index represents a single value. Hence, x(3:3) is an array with a single element, [x(3)].
When you index a multi-dimensional array, very often you want to let one or more dimensions be "spectators". In Yorick, you accomplish this by leaving the corresponding index blank:
x(3,) x(,5:7) y(,::-1,) x(,)
In these examples, x is a 2-D array, and y is a 3-D array. The first example, x(3,), represents the 1-D array of all the elements of x with first index 3. The second represents a 2-D array of all of the elements of x whose second indices are 5, 6, or 7. In the third example, y(,::-1,) is the y array with the elements in reverse order along its middle index. The fourth expression, x(,), means the entire 2-D array x, unchanged.
An index list is an array of index values. Use an index list to specify an arbitrary subset of an array: x([5,1,2,1]) means the four element array [x(5), x(1), x(2), x(1)]. The where function returns an index list:
list= where(x > 3.5) y= x(list)
These lines define the array y to be the subset of the x array, consisting of the elements greater than 3.5.
Like the result of an index range, the result of an index list is itself an array. However, the index list follows a more general rule: The dimensionality of the result is the same as the dimensionality of the index list. Hence, x([[5, 1], [2, 1]]) refers to the two dimensional array [[x(5), x(1)], [x(2), x(1)]]. The general rule for index lists is:
Dimensions from the dimensions of the index list; values from the array being indexed.
Note that the scalar index value is a special case of an index list according to this rule.
The rule applies to multi-dimensional arrays as well: If x is a five-by-nine array, then x(, [[5, 1], [2, 1]]) is a five-by-two-by-two array. And x([[5, 1], [2, 1]], 3:6) is a two-by-two-by-four array.
A binary operation applied between two arrays of numbers yields an array of results in Yorick -- the operation is performed once for each corresponding pair of elements of the operands. Hence, if rho and vol are each four-by-three arrays, the expression rho*vol will be a four-by-three array of products, starting with rho(1,1)*vol(1,1).
This extension of binary operations to array operands is not always appropriate. Instead of operating on the corresponding elements of arrays of the same shape, you may want to perform the operation between all pairs of elements. The most common example is the outer product of two vectors x and y:
outer= x*y(-,);
Here, if x were a four element vector, and y were a three element vector, outer would be the four-by-three array [[x(1)*y(1), x(2)*y(1), x(3)*y(1), x(4)*y(1)], [x(1)*y(2), x(2)*y(2), x(3)*y(2), x(4)*y(2)], [x(1)*y(3), x(2)*y(3), x(3)*y(3), x(4)*y(3)]].
In Yorick, this type of multiplication is still commutative. That is, x*y(-,) is the same as y(-,)*x. To produce the three-by-four transpose of the array outer, you would write x(-,)*y.
I call the - sign, when used as an index, a pseudo-index because it actually inserts an additional dimension into the result array which was not present in the array being indexed. By itself, the expression y(-,) is a one-by-three array (with the same three values as the three element vector y). You may insert as many pseudo-indices into a list of subscripts as you like, at any location you like relative to the actual dimensions of the array you are indexing. Hence, outer(-,-,,-,) would be a one-by-one-by-four-by-one-by-three array.
By default, a pseudo-index produces a result dimension of length one. However, by appending an index range to the - sign, separated by a colon, you can produce a new dimension of any convenient length:
x= span(-10, 10, 100)(,-:1:50); y= span(-5, 5, 50)(-:1:100,); gauss2d= exp(-0.5*(x^2+y^2))/(2.0*pi);
computes a normalized 2-D Gaussian function on a 100-by-50 rectangular grid of points in the xy-plane. The pseudo-index -:1:50 has 50 elements, and -:1:100 has 100. For a pseudo-index of non-unit length, the values along the actual dimensions are simply copied, so that span(-10, 10, 100)(,-:1:50) is a 100-by-50 array consisting of 50 copies of span(-10, 10, 100).
If only the result gauss2d were required, a single default pseudo-index would have sufficed:
gauss2d= exp(-0.5*( span(-10,10,100)^2 +
span(-5,5,50)(-,)^2 )) / (2.0*pi);
However, the rectangular grid of points (x,y) is often required -- as input to plotting routines for example.
Array index values are subtly asymmetric: An index of 1 represents the first element, 2 represents the second element, 3 the third, and so on. In order to refer to the last, or next to last, or any element relative to the final element, you apparently need to find out the length of the dimension.
In order to remedy this asymmetry, Yorick interprets numbers less than 1 relative to the final element of an array. Hence, x(1) and x(2) are the first and second elements of x, while x(0) and x(-1) are the last and next to last elements, and so on.
With this convention for negative indices, many Yorick programs can be written without the need to determine the length of a dimension:
deriv= (f(3:0)-f(1:-2)) / (x(3:0)-x(1:-2));
computes a point-centered estimate of the derivative of a function f with values known at points x. (A better way to compute this derivative is to use the pcen and dif range functions described below. See section Rank preserving (finite difference) range functions.)
In this example, the extra effort required to compute the array length would be slight:
n= numberof(f); deriv= (f(3:n)-f(1:n-2)) / (x(3:n)-x(1:n-2));
However, using the negative index convention produces faster code, and generalizes to multi-dimensional cases in an obvious way.
The negative index convention works for scalar index values and for the start or stop field of an index range (as in the example). Dealing with negative indices in an index list would slow the code down too much, so the values in an index list may not be zero or negative.
Many Yorick functions must work on arrays with an unknown number of dimensions. Consider a filter response function, which takes as input a spectrum (brightness as a function of photon energy), and returns the response of a detector. Such a function could be passed a spectrum and return a scalar value. Or, it might be passed a two dimensional array of spectra for each of a list of rays, and be expected to return the corresponding list of responses. Or a three dimensional array of spectra at each pixel of a two dimensional image, returning a two dimensional array of response values.
In order to write such a function, you need a way to say, "and all other indices this array might have". Yorick's rubber-index, .., stands for zero or more actual indices of the array being indexed. Any indices preceding a rubber-index are "left justified", and any following it are "right justified". Using this syntax, you can easily index an array, as long as you know that the dimension (or dimensions) you are interested in will always be first, or last -- even if you don't know how many spectator dimensions might be present in addition to the one your routine processes.
Thus, as long as the spectral dimension is always the final dimension of the input brightness array,
brightness(..,i)
will always place the i index in the spectral dimension, whether brightness itself is a 1-D, 2-D, or 3-D array.
Similarly, x(i,..) selects a value of the first index of x, leaving intact all following dimensions, if any. Constructions such as x(i,j,..,k,l) are also legal, albeit rarely necessary.
A second form of rubber-index collapses zero or more dimensions into a single index. The length of the collapsed dimension will be the product of the lengths of all of the dimensions it replaces (or 1, if it replaces zero dimensions). The symbol for this type of rubber index is an asterisk *. For example, if x were a five-by-three-by-four-by-two array, then x(*) would be a 1-D array of 120 elements, while x(,*,) would be a five-by-twelve-by-two array.
If the last actual index in a subscripted array is nil, and if this index does not correspond to the final actual dimension of the array, Yorick will append a .. rubber-index to the end of the subscript list. Hence, in the previous example, x() is equivalent to x(,..), which is equivalent to x(..), which is equivalent to simply x. This rule is the only rogue in Yorick's array subscripting stable, and I am mightily tempted to remove it on grounds of linguistic purity. When you mean, "and any other dimensions which might be present," use the .. rubber-index, not a nil index. Use a trailing nil index only when you mean, "and the single remaining dimension (which I know is present)."
Yorick has a special syntax for matrix multiplication, or, more generally, inner product:
A(,+)*B(+,) B(+,)*A(,+) x(+)*y(+) P(,+,,)*Q(,,+)
In the first example, A would be an L-by-M array, B would be an M-by-N array, and the result would be the L-by-N matrix product. In the second example, the result would be the N-by-L transpose of the first result. The general rule is that all of the "spectator" dimensions of the left operand precede the spectator dimensions of the right operand in the result.
The third example shows how to form the inner product of two vectors x and y of equal length. The fourth example shows how to contract the second dimension of a 4-D array P with the third dimension of the 3-D array Q. If P were 2-by-3-by-4-by-5 and Q were 6-by-7-by-3, the result array would be 2-by-4-by-5-by-6-by-7.
Unlike all of the other special subscript symbols (nil, -, .., and * so far), the + sign marking an index for use in an inner product is actually treated specially by the Yorick parser. The + subscript is a parse error unless the array (or expression) being subscripted is the left or right operand of a binary * operator, which is then parsed as matrix multiplication instead of Yorick's usual element-by-element multiplication. A parse error will also result if only one of the operands has a dimension marked by +. Both operands must have exactly one marked dimension, and the marked dimensions must turn out to be of equal length at run time.
The beauty of Yorick's matrix multiplication syntax is that you "point" to the dimension which is to be contracted by placing the + marker in the corresponding subscript. In this section and the following section, I introduce Yorick's range functions, which share the "this dimension right here" syntax with matrix multiplication. The topic in this section is the statistical range functions. These functions reduce the rank of an array, as if they were a simple scalar index, but instead of selecting a particular element along the dimension, a statistical range function selects a value based on an examination of all of the elements along the selected dimension. The statistical function is repeated separately for each value of any spectator dimensions. The available functions are:
The min, max, sum, and avg functions may also be applied using ordinary function syntax, which is preferred if you want the function to be applied across all the dimensions of an array to yield a single scalar result.
Given the brightness array representing the spectrum incident on a detector or set of detectors, the mxx function can be used to find the photon energy at which the incident light is brightest. Assume that the final dimension of brightness is always the spectral dimension, and that the 1-D array gav of photon energies (with the same length as the final dimension of brightness) is also given:
max_index_list= brightness(.., mxx); gav_at_max= gav(max_index_list);
Note that gav_at_max would be a scalar if brightness were a 1-D spectrum for a single detector, a 2-D array if brightness were a 3-D array of spectra for each point of an image, and so on.
An arbitrary index range (start:stop or start:stop:step) may be specified for any range function, by separating the function name from the range by another colon. For example, to select only a relative maximum of brightness for photon energies above 1.0, ignoring possible larger values at smaller energies, you could use:
i= min(where(gav > 1.0)); max_index_list= brightness(.., mxx:i:0); gav_at_max= gav(max_index_list);
Note the use of min invoked as an ordinary function in the first line of this example. (Recall that where returns a list of indices where some conditional expression is true.) In the second line, mxx:i:0 is equivalent to mxx:i:. Because of the details of Yorick's current implementation, the former executes slightly faster.
More than one range function may appear in a single subscript list. If so, they are computed from left to right. In order to execute them in another order, you must explicitly subscript the expression resulting from the first application:
x= [[1, 3, 2], [8, 0, 9]]; max_min= x(max, min); min_max= x(, min)(max);
The value of max_min is 3; the value of min_max is 2.
Because Yorick arrays almost invariably represent function values, Yorick provides numercal equivalents to the common operations of differential and integral calculus. In order to handle functions of several variables in a straightforward manner, these operators are implemented as range functions. Unlike the statistical range functions, which return a scalar result, the finite difference range functions do not reduce the rank of the subscripted array. Instead, they preserve rank, in the same way that a simple index range start:stop preserves rank. The available finite difference functions are:
The derivative dy/dx of a function y(x), where y and x are represented by 1-D arrays y and x of equal length is:
deriv= y(dif)/x(dif);
Note that deriv has one fewer element than either y or x. The derivative is computed as if y(x) were the piecewise linear function passing through the given points (x,y); there is one fewer line segment (slope) than point.
The values x and y can be called "point-centered", while the values deriv can be called "zone-centered". The zcen and pcen functions provide a simple mechanism for moving back and forth between point-centered and zone-centered quantities. Usually, there will be several reasonable ways to point-center zone centered data and vice-versa. For example:
deriv_pc1= deriv(pcen); deriv_pc2= y(dif)(pcen)/x(dif)(pcen); deriv_pc3= y(pcen)(dif)/x(pcen)(dif);
For a well-resolved function, the differences among these three arrays will be negligible. That is, the differences are second order in x(dif), which is often the order of the errors in the calculation that produced x and y in the first place. If x and y represent y(x) more accurately that this, then you must know a better model of the shape of y(x) than the simple piecewise linear model, and you should use that model to select deriv_pc1, deriv_pc2, and deriv_pc3, or some other expression.
An indefinite integral may be estimated using the trapezoidal rule:
indef_integ= (y(zcen)*x(dif))(psum);
Once again, indef_integ has one fewer point than x or y, because there is one fewer trapezoid than point. This time, however, indef_integ is not zone centered. Instead, indef_integ represents values at the upper (or lower) boundaries x(2:) (or x(:-1)). Often, you want to think of the integral of y(x) as a point centered array of definite integrals from x(1) up to each of the x(i). In this case (which actually arises more frequently), use the cum function instead of psum in order to produce a result def_integs with the same number of points as the x and y arrays:
def_integs= (y(zcen)*x(dif))(cum);
For single definite integrals, the matrix multiply syntax can be used in conjunction with the dif range function. For example, suppose you know the transmission fraction of the filter, ff, at several photon energies ef. That is, ff and ef are 1-D arrays of equal length, specifying a filter transmission function as the piecewise linear function connecting the given points. The final dimension of an array brightness represents an incident spectrum (any other dimensions represent different rays, say one ray per pixel of an imaging detector). The 1-D array gb represents the group boundary energies -- that is, the photon energies at the boundaries of the spectral groups represented in brightness. The following Yorick statements compute the detector response:
filter= integ(ff, ef, gb)(dif); response= brightness(..,+)*filter(+);
For this application, the correct interpolation routine is integ. The integrated transmission function is evaluated at the boundaries gb of the groups; the effective transmission fraction for each group is the difference between the integral at the upper bin boundary and the lower bin boundary. The dif range function computes these pairwise differences. Since the gb array naturally has one more element than the final dimension of brightness (there is one more group boundary energy than group), and since dif reduces the length of a dimension by one, the filter array has the same length as the final dimension of brightness.
Note that the cum function cannot be used to integrate ff(ef), because the points gb at which the integral must be evaluated are not the same as the points ef at which the integrand is known. Whenever the points at which the integral is required are the same (or a subset of) the points at which the integrand is known, you should perform the integral using the zcen, dif, and cum index functions, instead of the more general integ interpolator.
The function sort returns an index list:
list= sort(x);
The list has the same length as the input vector x, and list(1) is the index of the smallest element of x, list(2) the index of the next smallest, and so on. Thus, x(list) will be x sorted into ascending order.
By returning an index list instead of the sorted array, sort simplifies the co-sorting of other arrays. For example, consider a series of elaborate experiments. On each experiment, thermonuclear yield and laser input energy are measured, as well as fuel mass. These might reasonably be stored in three vectors yield, input, and mass, each of which is a 1D array with as many elements as experiments performed. The order of the elements in the arrays maight naturally be the order in which the experiments were performed, so that yield(3), input(3), and mass(3) represent the third experiment. The experiments can be sorted into order of increasing gain (yield per unit input energy) as follows:
list= sort(yield/input); yield= yield(list); input= input(list); mass= mass(list);
The inverse list, which will return them to their orginal order, is:
invlist= list; /* faster than array(0, dimsof(list)) */ invlist(list)= indgen(numberof(list));
The sort function actually sorts along only one index of a multi-dimensional array. The dimensions of the returned list are the same as the dimensions of the input array; the values are indices relative to the beginning of the entire array, not indices for the dimension being sorted. Thus, x(sort(x)) is the array x sorted so that the elements along its first dimension are in ascending order. In order to sort along another dimension, pass sort a second argument -- x(sort(x,2)) will be x sorted into increasing order along its second dimension, x(sort(x,3)) along its third dimension, and so on.
A related function is median. It takes one or two arguments, just like sort, but its result -- which is, of course, the median values along the dimension being sorted -- has one fewer dimension than its input.
In a carefully designed Yorick program, array dimensions usually wind up where you need them. However, you may occasionally wind up with a five-by-seven array, instead of the seven-by-five array you wanted. Yorick has a very general transpose function:
x623451= transpose(x123456); x561234= transpose(x123456, 3); x345612= transpose(x123456, 5); x153426= transpose(x123456, [2,5]); x145326= transpose(x123456, [2,5,3,4]); x653124= transpose(x123456, [2,5], [1,4,6]);
Here, x123456 represents a six-dimensional array (hopefully these will be rare). The same array with its first and last dimensions transposed is called x623451; this is the default result of transpose. The transpose function can take any number of additional arguments to describe an arbitrary permutation of the indices of its first argument -- the array to be transposed.
A scalar integer value, as in the second and third lines, represents a cyclic permutaion of all the dimensions; the first dimension of the input becomes the Nth dimension of the result, for an argument value of N.
An array of integer values represents a cyclic permutation of the specified dimensions. Hence, in the fifth example, [2,5,3,4] means that the second dimension becomes the fifth, the fifth becomes the third, the third becomes the fourth, and the fourth becomes the second. An arbitrary permutation may be built up of a number of cyclic permutations, as shown in the sixth example.
One additional problem can arise in Yorick -- you may not know how many dimensions the array to be transposed has. In order to deal with this possibility, either a scalar argument or any of the numbers in a cyclic permutation list may be zero or negative to count from the last dimension. That is, 0 represents the last dimension, -1 the next to last, -2 the one before that, and so on.
Typical transposing tasks are: (1) Move the first dimension to be last, and all the others back one cyclically. (2) Move the last dimension first, and all the others forward one cyclically. (3) Transpose the first two dimensions, leaving all others fixed. (4) Transpose the final two dimensions, leaving all others fixed. These would be accomplished, in order, by the following four lines:
x234561= transpose(x123456, 0); x612345= transpose(x123456, 2); x213456= transpose(x123456, [1,2]); x123465= transpose(x123456, [0,-1]);
Two arrays need not have identical shape in order for a binary operation between them to make perfect sense:
y= a*x^3 + b*x^2 + c*x + d;
The obvious intent of this assignment statement is that y should have the same shape as x, each value of y being the value of the polynomial at the corresponding y.
Alternatively, array valued coefficients a, b, ..., represent an array of several polynomials -- perhaps Legendre polynomials. Then, if x is a scalar value, the meaning is again obvious; y should have the same shape as the coefficient arrays, each y being the value of the corresponding polynomial.
A binary operation is performed once for each element, so that a+b, say, means
a(i,j,k,l) + b(i,j,k,l)
for every i, j, k, l (taking a and b to be four dimensional). The lengths of corresponding indices must match in order for this procedure to make sense; Yorick signals a "conformability error" when the shapes of binary operands do not match.
However, if a had only three dimensions, a+b still makes sense as:
a(i,j,k) + b(i,j,k,l)
This sense extends to two dimensional, one dimensional, and finally to scalar a:
a + b(i,j,k,l)
which is how Yorick interpreted the monomial x^3 in the first example of this section. The shapes of a and b are conformable as long as the dimensions which they have in common all have the same lengths; the shorter operand simply repeats its values for every index of a dimension it doesn't have. This repitition is called "broadcasting".
Broadcasting is the key to Yorick's array syntax. In practical situations, it is just as likely for the a array to be missing the second (j) dimension of the b array as its last (l) dimension. To handle this case, Yorick will broadcast any unit-length dimension in addition to a missing final dimension. Hence, if the a array has a second dimension of length one, a+b means:
a(i,1,k,l) + b(i,j,k,l)
for every i, j, k, l. The pseudo-index can be used to generate such unit length indices when necessary (see section Creating a pseudo-index).
You should strive to write Yorick programs in such a way that you never need to refer to the lengths of array dimensions. Array dimensions are usually of no direct significance in a calculation; your programs will tend to be clearer if only the arrays themselves appear.
In practice, unfortunately, you can't always get by without mentioning dimension lengths. Two functions are important:
The array function (see section Creating Arrays), and the more arcane functions add_variable, add_member, and reshape all have parameter lists ending with one or more "dimension list" parameters. Each parameter in a dimension list can be either a scalar integer value, representing the length of a single dimension, or a list of integers in the format returned by dimsof to represent zero or more dimensions. Several arguments can be used to build up a complicated dimension list:
x1= array(0.0, 9, 2, 6); /* 9-by-2-by-6 array of 0.0 */ x2= array(0.0, [3,9,2,6]); /* another 9-by-2-by-6 array */ x3= array(0.0, 9, [0], [2,2,6]); /* ...and yet another */ flux= array(0.0, 3, dimsof(z), numberof(groups));
In the final example, the flux might represent the three components of a flux vector, at each of a number of positions z, and for each of a number of photon energies groups. The first dimension of flux has length three, corresponding to the three components of each flux vector. The last dimension has the same length as the groups array. In between are zero or more dimensions -- whatever the dimensions of the array of positions z.
By using the rubber index syntax (see section Using a rubber index), you can extract meaningful slices of the flux array without ever needing to know how many dimensions z had, let alone their lengths.