Wednesday, March 31, 2010

Efficient Bulk Processing: The Command Line

The most efficient way to do bulk processing (such as doing the same thing to lots of files) of the kind we need for simulations is at the command prompt. This post describes the basics of the command prompt.

What is a "Shell"?
A shell is a wrapper around the low level of the actual operating system. Linus/Unix users are generally familiar with shells. They even have a choice-the c-shell csh, the bourne shell (propriety), or it's more famous cousin, the "bourne again shell" bash (free), the tschell tsh, etc. In this post, I'll restrict myself to bash, as it is arguably the the most commonly used.

In the Windows world, it is not as well-known that there are also several shells available. The original, dinosaur is the MS DOS shell command.com. Windows XP introduced another shell, cmd.exe. Most recently, Microsoft has developed PowerShell, which is a substantially more powerful animal. I will discuss PowerShell in another post.

How Do I Get a Shell in Unix/Linux?
In the *nix world, this probably isn't an issue. Most users spend a good deal of time in the shell environment. If you are working in a desktop enviornment such as KDE or GNOME, you are looking for a program called called "terminal", "konsole", or "xterm".

How Do I Get a Shell in Windows?
In the ancient, pre-windows days, this simply wasn't a question. The computer booted and you were faces with the DOS command prompt.

In the Windows world, there are two solutions: one is to click "Start", then click "run", and type "command" (or "command.com" and hit return. A black box will open on the desktop.
The second way is for click "Start", then click "All Programs" "Accessories" and find the black-box icon called "Command Prompt". Clicking it will bring up the same black box.
Because We'll do a lot of work at the command prompt, it is worth having the command Prompt on your Windows desktop. Click on the icon and drag it onto the desktop. This will create a short-cut that you can double click.

Both bash and command are usually configured to display the current directory. You will see the current directory followed by the prompt symbol: "$" in bash, ">" in windows.

To move from place to place, we use the "cd" command for both operating systems. The only difference is that *nix uses a forward slash "/" to separate directories, while Windows uses a backslash "\". Sometimes Windows will accept a forward slash, but the rules for when it will and won't are more trouble than they are worth. When working at the Windows command line, just use the backslash and move on with your day.

To list the contents of the current directory, use "dir" on Windows and "ls" on *nix.
Navigate to the directory with your output, using the "cd" command. Say your output is a folder on your Windows desktop. The Windows command prompt opens in your home directory: usually
XP--typically c:\documents and settings\xxx
Vista and Windows 7--typically c:\users\xxx
where "xxx" is your user name (typically "owner" if you haven't changed it on a home machine, your login ID on a work machine.

To move around, use the "cd" command. this works on both *nix shells like bash and windows. To go down one directory, simply type "cd" , a space, and the name of the directory. Make sure to include the space, or you'll get an error.

NOTE: the windows command shell IS NOT case sensitive. The Unix/Linux shell IS case sensitive. Memorize that. Get used to it. And move on with your life.

Under XP you may see some rubbishy stuff before the ">" such as "C:\DOCUME~1\JDONOG~1>" Due to an old limitation of 8 characters for a name in the windows command prompt, the command prompt abbreviates long directory names with the first 6 characters, a tilde ("~") and a number. The number is generally 1, unless you have two or more folders that begin with the same name. In that case, you'll see "~2", "~3", etc., where the number indicates the order in which the folders were created. If you simply want to move to a directory with spaces in the name, like "Documents and Settings", "cd Documents and settings" will work. However, if you want to move down 2 directories, to "Documents and Settings\desktop", typing "cd Documents and settings\desktop" will yield an error.

This can be a pain. There are 3 simple solutions:
  • Move one directory at a time
  • Put the name of the directory with spaces in double quotes: eg. 'cd "Documents and settings"\desktop'
  • Figure out which is which is to use the "dir" command (which lists the names of the files in the directory). Say I have two folders, "Documents and settings" and "Documents to give to my cat to shred". I can type "dir DOCUME~1*" and the system will return "Documents and Settings", or type "dir DOCUME~2*" and get back "Documents to give to my cat to shred" Notice the aterisk at the end of the commands. Asterisk is a "wildcard" that will match any collection of characters. Using it is especially important for names that have periods in the ("jython-2.5.2"). It is important to include it.
We have covered moving down one directory ("cd child") and moving directly to a specific directory "cd \xxx\yyy". To move up on directory, we can make use of a special names, "." and ".." Whever you've type "dir" or "ls", at the top of the list are two entries with no names, just a period (".") and two periods (".."). These indicate special references. "." refers to the current directory, and ".." refers to the parent of the current directory. So, the command "cd .." will take you up one directory.

One final note for this post: If "Documents and Settings" is the only directory that begins with the characters "Doc", I can save myself some typing by just using "cd Doc*" and Windows will figure it out.

Next time we'll begin using this knowledge of the command line to extract specific information from an output file.

Sunday, February 14, 2010

Seedy Random Number Generation

One of the advantages of the multiple small runs strategy (described last time) was that it makes it easy to replicate any individual run. That gives us the capability to examine the details of a specific run, such as one in which the iterations did not converge, or that aborted.

There are two ways to do this; keep all of the individual data sets, or be able to re-generate each and every data set whenever we want.

Whether we are using a strategy of generating all of the data, or of sampling from a large dataset, or re-sampling (for a jackknife- or bootstrap-based study), we will need to generate a series of (pseudo-) random numbers. All of the pseudo-random number generators make use of a seed. We input some value, and then a series of numbers is generated. Once the seed is provided, the series is completely determined. Most implementations provide two possibilities for the seed:
  1. explicitly provide a seed value
  2. implicitly have the seed generated for us
Most generators make it simpler to one or the other of these, but pretty much every generator I have used allows both options. A common strategy for implicit seeding is to use the system's clock.

The key to being able to replicate runs is to capture the seed for every run. We can do this by using an implicit seed, and recording that seed (say by appending it to a file). I have found this approach to be extremely error prone, and don't advise it. Instead, my approach is to generate a seed file by using a random number generator (often with an implicit seed) to generate a file containing a large number of seed values., one for each replication. Then I copy the seed file to a backup (a separate disk or flash drive). In general, I will have one seed file per condition of the study. That file should contain at least the number of iterations I want. In general, I have some overage (at least 10% if I'm optimistic, 100% if I'm feeling pessimistic). That way if a specific run hiccups, my simulations can continue without running out of seeds.

Say that I call the file something clever, like "seed.dat". Each replication:
  1. creates a backup copy of the file--"seed.bak",
  2. reads the N values in the seed file,
  3. removes the top value, and
  4. writes N-1 values (without the top value) back to same file (seed.dat) for the next replication's use.
The purpose of the backup is to be able to kill the current run. I often want to pause the simulation to do some interim data analyses. Using this strategy, I can copy seed.bak to seed.dat, and resume the simulation.

This approach gives me two advantages:
  1. I can don't have to keep the individual data sets around, reducing the possibility that I'l run into storage problems
  2. I can re-generate any specific data set by simply using the corresponding seed value
  3. Finally, when I go to archive the files associated with the study, I only have to save the seed files, secure in the knowledge that I can recreate anything (from a single run to the entire study) at a moment's notice.

Saturday, February 13, 2010

Doing the Replications: (Again, and again and again...)

We'll consider how to get the data for the replications in another installment. But first, I want to talk about the general approach to doing a large simulation study.

There are two schools of though there:
  1. One BIG RUN. It goes for 4 weeks, performs all 50,000 replications, and then prints out a 7. This kind of thing is routinely done in the physical sciences; I've read stories about work on subatomic particles that do this one-long-run approach.
  2. 50,000 small runs. Then we pull the results from all of the files, and create a big data set.

I fall firmly into the second camp. There a several reasons.
  1. I have had the experience of the BIG RUN go for 3 weeks, 6 days, and facilities turns off the power for some needed maintenance. This always results in me saying a long spew of bad words, that my mother would not approve of.
  2. Somewhere during the process, a few runs hiccup. It may be that iterations fail to converge; that the dataset is defective; something. Now I have 49,936 replications, and no idea what happened to the other 64 runs. And actually, it may be that the 64 runs are trying to tell me something (maybe more interesting that what I thought I was looking for). However, I don't know what dataset produced the problem, why, or anything else. In short, I'm SOL [1].
  3. There may be something else about the data, like the skewness of the statistic, or the standard error of the estimated standard errors; (it's always weird) that a journal reviewer wants to know about before approving the manuscript. If you've done the BIG RUN, you have little choice but to redo it--more bad words. If I have the dataset of individual runs, I can re-analyze it. Sometimes, the reviewer is on to something, or they may be full of hot air. Either way, I can go back to the data. With the BIG RUN, it's not an option.
  4. It is easier to monitor the progress. If the conference paper is due in 3 days and I'm only on replication 17, it is time to re-group.
  5. (and this is my favorite), if I structure the replications so that each of the conditions is executed in round-robin fashion, I can do some preliminary analysis after 10 or 50 replications. Based on these, I can do things like power analyses for the effects I'm looking for and realize that I need more (or fewer!) replications. Or, things may look different than I expected, and I can shut it down, call it a pilot, and go off in another direction.
In short, I've been burned a bunch of different ways by the BIG RUN. In the physical sciences like subatomic physics, they may understand the domain to a depth that this approach makes sense. In statistical simulations, there are some domains where that is true, but in general we are working with newer, lesser known quantities, and circumspection is well-advised.

There are ways to reduce the risks/downsides of the BIG RUN. But, overall, I find that doing a bunch of little runs gives me many more degrees of freedom to adapt to conditions that arise. I've tried it both ways, and lots of small runs is better.

There are downsides to doing small runs:
  1. Disk space can become a problem,
  2. You often can't use simulators that are baked into the software,
  3. Figuring out how to run a commercial program (especially one that is based on a GUI) can pose a major problem,
  4. With so many tiny pieces, it is easy to get lost in the detail.
  5. Pulling the results together from lots of individual files can be a major undertaking.
In the next post, I'll talk about my general approach. Then we'll consider, how to solve some of the trickier problems.


[1] American slang meaning "Sh*t Out of Luck."

Sunday, January 24, 2010

Fortran Unit Testing Framework (part 6): Removing Some Duplication

As we look back at our progress, there is a substantial amount of duplication between the three test methods in our suite of tests for fxxUnit. The routines are:
TSTOUT, TASRTT, and TASRTF. Now we're going to refactor to clean up that duplication.
All 3 SUBROUTINEs perform the same set of actions:
  1. Open a scratch file for our results
  2. Change the output unit to a scratch file,
  3. Call the SUBROUTINE that we are testing with the argument .TRUE.
  4. Call the SUBROUTINE that we are testing with the argument .FALSE.
  5. Rewind the scratch file
  6. Read back the results from the scratch file,
  7. Compare what was written to our expectation, outputting the result character "." or "F",
  8. Clean up by closing the scratch file and re-setting the output
To remind ourselves of the structure, here is the TSTOUT SUBROUTINE,


C
C TEST CASE FOR THE ROUTINE THAT PRINTS OUT THE
C CHARACTER BASED ON THE VALUE OF THE TEST
SUBROUTINE TSTOUT()
COMMON /OUTP/KOUT
INTEGER KOUT
INTEGER I
INTEGER OUTPUT
INTEGER EXPECT
INTEGER RESULT
INTEGER KTEMP
DIMENSION RESULT(2)
DIMENSION EXPECT(2)
DATA (EXPECT(I), I = 1, 2)/1H., 1HF/,OUTPUT/10/
KTEMP = KOUT
KOUT = OUTPUT
CALL RESOUT(.TRUE.)
CALL RESOUT(.FALSE.)
REWIND(OUTPUT)
READ(OUTPUT, 100) (RESULT(I), I = 1, 2)
100 FORMAT(2A1)
CLOSE(OUTPUT)
OPEN(UNIT = OUTPUT, STATUS = 'SCRATCH')
DO 600 I = 1, 2
IF (EXPECT(I) .EQ. RESULT(I)) GO TO 300
WRITE(KTEMP, 200)
200 FORMAT(1HF, $)
GO TO 500
300 CONTINUE
WRITE(KTEMP, 400) IPASS
400 FORMAT(1H., $)
500 CONTINUE
600 CONTINUE
KOUT = KTEMP
RETURN
END


First, we pull out step 7 (the loop) into its own SUBROUTINE. Taking a look at the actual code, we see that it is comparing two INTEGER arrays. So we name it IAASRT:

C
C COMPARES THE RESULT OF TWO INTEGER ARRAYS, WRITING THE RESULT TO
C *KTEMP*. ONE RESULT IS WRITTEN PER ELEMENT OF THE ARRAY.
SUBROUTINE IAEQUL(EXPECT, RESULT, N, KTEMP)
INTEGER EXPECT
INTEGER RESULt
INTEGER N
INTEGER KTEMP
DIMENSION EXPECT(1)
DIMENSION RESULT(1)
DO 500 I = 1, N
IF (EXPECT(I) .EQ. RESULT(I)) GO TO 200
WRITE(KTEMP, 100)
100 FORMAT(1HF, $)
GO TO 400
200 CONTINUE
WRITE(KTEMP, 300)
300 FORMAT(1H., $)
400 CONTINUE
500 CONTINUE
RETURN
END


We replace the original logic with calls to the SUBROUTINE in TSTOUT, TASRTT, and TASRTF. Then we run the tests to verify that we haven't broken anything.

Next, we pull the remaining steps into a second SUBROUTINE, called COMPTF, because it will call a SUBROUTINE that we provide it, with .TRUE., and then with .FALSE. and compare the results written out to our expectations.


C
C COMPARES THE CHARACTERS WRITTEN OUT BY THE SUBROUTINE *FUNCT*
C TO THOSE EXPECTED.
SUBROUTINE COMPTF(FUNCT, PRINTS)
COMMON /OUTP/KOUT
INTEGER KOUT
INTEGER PRINTS
INTEGER KTEMP
INTEGER OUTPUT
INTEGER RESULT
INTEGER I
DIMENSION PRINTS(2)
DIMENSION RESULT(2)
DATA OUTPUT/10/
KTEMP = KOUT
KOUT = OUTPUT
OPEN(UNIT = OUTPUT, STATUS = 'SCRATCH')
CALL FUNCT(.TRUE.)
CALL FUNCT(.FALSE.)
REWIND(OUTPUT)
READ(OUTPUT, 100) (RESULT(I), I = 1, 2)
CLOSE(OUTPUT)
100 FORMAT(2A1)
CALL IAEQUL(PRINTS, RESULT, 2, KTEMP)
KOUT = KTEMP
RETURN
END


The unusual piece here is the first argument, FUNCT. FORTRAN 66 spec allows us to pass a symbolic reference to a SUBROUTINE or FUNCTION. It is similar to a function pointer in C or C++. We give the routine a symbolic name FUNCT in this case, and call it from within our COMPTF SUBROUTINE. The only restriction is that we have to specify the name of the routine being passed as EXTERNAL within the calling code.

With these two routines pulled out, the code to TSTOUT, TASRTT, and TASRTF becomes quite simple:


C
C TEST CASE FOR THE ROUTINE THAT PRINTS OUT THE
C CHARACTER BASED ON THE VALUE OF THE TEST
SUBROUTINE TSTOUT()
EXTERNAL RESOUT
INTEGER EXPECT
DIMENSION EXPECT(2)
INTEGER I
DATA (EXPECT(I), I = 1, 2)/1H., 1HF/
CALL COMPTF(RESOUT, EXPECT)
RETURN
END
C
C TEST CASE FOR ASRTT, THE ASSERTION THAT THE ARGUMENT IS TRUE
SUBROUTINE TASRTT()
EXTERNAL ASRTT
INTEGER EXPECT
DIMENSION EXPECT(2)
INTEGER I
DATA (EXPECT(I), I = 1, 2)/1H., 1HF/
CALL COMPTF(ASRTT, EXPECT)
RETURN
END
C
C TEST CASE FOR ASRTF, THE ASSERTION THAT THE ARGUMENT IS FALSE
SUBROUTINE TASRTF()
EXTERNAL ASRTF
INTEGER EXPECT
DIMENSION EXPECT(2)
DIMENSION ARGS(2)
INTEGER I
DATA (EXPECT(I), I = 1, 2)/1HF, 1H./
CALL COMPTF(ASRTF, EXPECT)
RETURN
END


Summary:

In this installment we refactored our code that tests our testing framework to remove some duplicate logic.
  • We pulled out the comparison into one SUBROUTINE, and, realizing that its logic was more general than the original context, gave it a name that reveled the general function
  • We then moved the rest of the duplicate code into a second SUBROUTINE.
  • We made use of the EXTERNAL capacity of FORTRAN 66 to simplify the code even further.

There are couple more things that we could change. Two possibilities that immediately come to mind are
Pulling the code to swap output units, and
Replacing the output characters "." and "F" with meaningful names that we use throughout.

The first only occurs in the COMPTF SUBROUTINE, so it doesn't bother me. The second is worth considering for another day.

One final, important point:
The discussion above is quite sanitized compared to the actual process I used in these refactorings. I omitted all of the intermediate steps to save space and increase the pace to something that would not be too annoying. Every small change (creating new code, calling it, removing the original code, etc.) was followed by running the tests. Several times (at least 5) I discovered that I had introduced a error, usually by forgetting to declare the type of a variable in the SUBROUTINE, or by omitting a variable in the argument list of the new SUBROUTINE. FORTRAN is a major pain here, because the compiler creates new variables and sets their types according the the implicit conventions.

In each case, it introduced a bug that I would have missed if I hadn't had a suite of tests. In the microcosm of this little example, I had the value of unit testing reinforced to me again and again.

Without a good test suite, writing correct code is hideously hard. In the context of a simulation, in which we have so many moving parts, recognizing errors based upon the final output is nearly impossible. Unit testing is a vital part of conducting high quality, professional simulations.

Wednesday, January 6, 2010

Fortran Unit Testing Framework (part 5): Writing Assertions

In this installment, we finally get around to writing the assertions. Recall that from part 2 we determined that the assertions would be SUBROUTINEs, not FUNCTIONs. We also decided to adopt the Fortran convention that the first letter will indicate what type of argument is expected. The middle 4 letters will be ASRT (for assert) and the final letter will contain information about the assertion ("T" for .TRUE., "F" for .FALSE., or "E" for .EQ.).

Before we begin, we notice that RESOUT requires the output unit number. So, we are going to have to write
ASRTT(RESULT, OUTPUT)
That gets intrusive pretty quickly. Let's pull the output unit out, into a COMMON block named OUTP:


C
C
PROGRAM MAIN
COMMON /OUTP/KOUT
INTEGER KOUT
...
STOP
END
...
C
C SUBROUTINE RESOUT
C WRITES OUT THE RESULT OF THE TEST
C '.' IF THE ARGUMENT IS .TRUE.
C 'F' OTHERWISE
C PARAM - LOGICAL *RESLT* WHETHER OR NOT THE TEST HAS PASSED
SUBROUTINE RESOUT(RESLT)
COMMON /OUTP/KOUT
INTEGER KOUT
LOGICAL RESLT
INTEGER RESCHR
INTEGER CHARA
CHARA = RESCHR(RESLT)
WRITE(KOUT, 100) CHARA
100 FORMAT(A1, $)
RETURN
END


We only include the block in the routines what need access to it, effectively limiting its scope. Now we can write assertions along the following lines:


    LOGICAL TEST
...
TEST = .TRUE.
...
CALL ASRTT(TEST)
...


which is much more natural.
With the refactoring done, we begin with the test for assert true (
ASRTT). The test code is pretty simple:


C
C TEST CASE FOR ASRTT, THE ASSERTION THAT THE ARGUMENT IS TRUE
SUBROUTINE TASRTT()
COMMON /OUTP/KOUT
EXTERNAL ASRTT
INTEGER KOUT
INTEGER I
INTEGER OUTPUT
INTEGER EXPECT
INTEGER RESULT
INTEGER KTEMP
DIMENSION RESULT(2)
DIMENSION EXPECT(2)
DATA (EXPECT(I), I= 1, 2)/1H., 1HF/, OUTPUT/10/
KTEMP = KOUT
KOUT = OUTPUT
OPEN(UNIT = OUTPUT, STATUS = 'SCRATCH')
CALL ASRTT(.TRUE.)
CALL ASRTT(.FALSE.)
REWIND(OUTPUT)
READ(OUTPUT, 100) (RESULT(I), I = 1, 2)
CLOSE(OUTPUT)
100 FORMAT(2A1)
DO 600 I = 1, 2
IF (EXPECT(I) .EQ. result(I)) GO TO 300
WRITE(KTEMP, 200)
200 FORMAT(1HF, $)
GO TO 500
300 CONTINUE
WRITE(KTEMP, 400) IPASS
400 FORMAT(1H., $)
500 CONTINUE
600 CONTINUE
KOUT = KTEMP
RETURN
END


This is essentially a repeat of the test we developed for RESOUT. Making the test pass simply requires a call to RESOUT.


C
C BASIC ASSERT_TRUE
C
C PARAM: LOGICAL TSTVAL THE VALUE TO BE TESTED
C PARAM: INTEGER OUTPUT THE UNIT TO WRITE THE RESULT TO
SUBROUTINE ASRTT(TSTVAL)
LOGICAL TSTVAL
CALL RESOUT(TSTVAL, OUTPUT)
RETURN
END


The test for ASRTF is the same, with the sole exception that the elements of EXPECT array are reversed.


C
C TEST CASE FOR ASRTF, THE ASSERTION THAT THE ARGUMENT IS TRUE
SUBROUTINE TASRTF()
COMMON /OUTP/KOUT
EXTERNAL ASRTT
INTEGER KOUT
INTEGER I
INTEGER OUTPUT
INTEGER EXPECT
INTEGER RESULT
INTEGER KTEMP
DIMENSION RESULT(2)
DIMENSION EXPECT(2)
DATA (EXPECT(I), I= 1, 2)/1H., 1HF/, OUTPUT/10/
KTEMP = KOUT
KOUT = OUTPUT
OPEN(UNIT = OUTPUT, STATUS = 'SCRATCH')
CALL ASRTF(.TRUE.)
CALL ASRTF(.FALSE.)
REWIND(OUTPUT)
READ(OUTPUT, 100) (RESULT(I), I = 1, 2)
CLOSE(OUTPUT)
100 FORMAT(2A1)
DO 600 I = 1, 2
IF (EXPECT(I) .EQ. RESULT(I)) GO TO 300
WRITE(KOUT, 200)
200 FORMAT(1HF, $)
GO TO 500
300 CONTINUE
WRITE(KOUT, 400) IPASS
400 FORMAT(1H., $)
500 CONTINUE
600 CONTINUE
KOUT = KTEMP
RETURN
END


Making the test pass involves simply negating the argument and passing it to ASRTT.


C
C BASIC ASSERT_FALSE
C
C PARAM: LOGICAL TSTVAL THE VALUE TO BE TESTED
C PARAM: INTEGER OUTPUT THE UNIT TO WRITE THE RESULT TO
      SUBROUTINE ASRTF(TSTVAL)
      LOGICAL TSTVAL
      CALL ASRTT(.NOT. TSTVAL)
      RETURN
      END


At this point we have repeated the code to set the output unit, set our expectations, and read the results back in 3 times. This is a major violation of the DRY principle, and so we really should refactor the code to remove the duplication. We will make a note to fix that next time.

As a second point, we really would like to drop the writing and reading back in part of the tests. We tested that as part of RESOUT; it seems foolish to keep testing it. Instead, let's factor out the common functionality into the the test for RESOUT, and just call that for each of our assertions. Our main logic will reside in a set of LOGICAL functions, that can be tested separately. The basis of each function name will be EQUAL, with the prefix letter that indicates the type of the argument:
  • IEQUAL - compares 2 INTEGER values
  • LEQUAL - compares 2 LOGICAL values
  • AEQUAL - compares 2 REAL values
  • DEQUAL - compares 2 DOUBLE PRECISION values
  • CEQUAL - compares 2 COMPLEX values


C
C TESTS THE FUNCTIONING OF THE INTEGER COMPARISON,
C LOGICAL FUNCTION IEQUAL
      SUBROUTINE TESTI()
      LOGICAL IEQUAL
      CALL ASRTT(IEQUAL(1, 1))
      CALL ASRTF(IEQUAL(1, 2))
      RETURN
      END



The tests verify that 2 integers are equal return .TRUE., and two that are not equal return .FALSE. We then make use of the functionality that we have created already in RESCHR and RESOUT. The code to make these pass is simple:


C
C DETERMINE IF TWO INTEGERS ARE EQUAL:
C
C INPUT EXPECT (INTEGER) - THE EXPECTED VALUE
C INPUT ACTUAL (INTEGER) - THE VALUE TO BE COMPARED TO EXPECT
C OUTPUT EXPECT .EQ. ACTUAL
      LOGICAL FUNCTION IEQUAL(EXPECT, ACTUAL)
      INTEGER EXPECT
      INTEGER ACTUAL
      IEQUAL = EXPECT .EQ. ACTUAL
      END


The test for LOGICAL variables is similarly simple:


C
C TESTS THE FUNCTIONING OF THE LOGICAL COMPARISON,
C LOGICAL FUNCTION LEQUAL
      SUBROUTINE TESTL()
      LOGICAL LEQUAL
      CALL ASRTT(LEQUAL(.TRUE., .TRUE.))
      CALL ASRTT(LEQUAL(.FALSE., .FALSE.))
      CALL ASRTF(LEQUAL(.TRUE., .FALSE.))
      CALL ASRTF(LEQUAL(.FALSE., .TRUE.))
      RETURN
      END


The assertion code has s similarly simple implementation:


C
C DETERMINE IF TWO LOGICALS ARE EQUIVALENT:
C
C INPUT EXPECT (LOGICAL) - THE EXPECTED VALUE
C INPUT ACTUAL (LOGICAL) - THE VALUE TO BE COMPARED TO EXPECT
C OUTPUT EXPECT .EQ. ACTUAL
C
C NOTE: BECAUSE F77 WANTS ".EQV." INSTEAD OF
C THE F66 ".EQ." WHEN COMPARING LOGICALS, THE
C COMPARISON IS WRITTEN IN A SOMEWHAT CONVOLUTED WAY
      LOGICAL FUNCTION LEQUAL(EXPECT, ACTUAL)
      LOGICAL EXPECT
      LOGICAL ACTUAL
      LEQUAL = (EXPECT .AND. ACTUAL) .OR.
     & (.NOT. EXPECT .AND. .NOT. ACTUAL)
      END


As the comment indicates, there is one slight complication. FORTRAN 77 prefers the use of the .EQV. relation for comparing two LOGICAL variables. To silence the compiler warnings (from g77), the routine is written in the current, somewhat unintuitive way.
For the comparison of REALs,there is a slight twist. Floating point arithmetic is imprecise, and so the exact comparison of 2 floating point numbers (REAL, DOUBLE PRCISION, or COMPLEX) is a bad idea; we may get two comparisons that should be equal that are not. For example:

C
C DEMONSTRATION OF FLOATING POINT
      PROGRAM COMP
      REAL ONE
      REAL THIRD
      REAL THREE
      LOGICAL RESULT
      ONE = 1.0
      THIRD = 1.0 / 3.0
      WRITE (6, 100) ONE, THIRD
  100 FORMAT(6HONE = , F15.10, 1X, 13H 1.0 / 3.0 = , F15.10)
      THREE = THIRD * 3.0
      WRITE (6, 200) THIRD
  200 FORMAT(20H 3.0 * (1.0 / 3.0) = , F15.10)
      RESULT = ONE .EQ. THIRD
      WRITE (6, 300) RESULT
  300 FORMAT(32H 1.0 .EQ. (3.0 * (1.0 / 3.0)) = ,L1)
      STOP
      END


yields the output


ONE =    1.0000000000  1.0 / 3.0 =    0.3333333433
3.0 * (1.0 / 3.0) = 0.3333333433
1.0 .EQ. (3.0 * (1.0 / 3.0)) = F


To account for this, xUnit framework assertions for floating point values add a third parameter, the positive value EPSILN. If the two values are within EPSILN of each others, the assertion returns true. Otherwise, it is false. This gives us the test cases for AEQUAL and DEQUAL:

C
C TESTS THE FUNCTIONING OF THE COMPARISON OF REALS,
C LOGICAL FUNCTION AEQUAL
      SUBROUTINE TESTA()
      LOGICAL AEQUAL
      REAL EXPECT
      REAL ACTUL1
      REAL ACTUL2
      REAL ACTUL3
      REAL EPS1
      REAL EPS2
      DATA EXPECT/1.0/, ACTUL1/1.0/, ACTUL2/1.00001/, ACTUL3/0.99999/
      DATA EPS1/1.E-1/, EPS2/1.E-7/
      CALL ASRTT(AEQUAL(EXPECT, ACTUL1, EPS1))
      CALL ASRTT(AEQUAL(EXPECT, ACTUL1, EPS2))
      CALL ASRTT(AEQUAL(EXPECT, ACTUL2, EPS1))
      CALL ASRTF(AEQUAL(EXPECT, ACTUL2, EPS2))
      CALL ASRTT(AEQUAL(EXPECT, ACTUL3, EPS1))
      CALL ASRTF(AEQUAL(EXPECT, ACTUL3, EPS2))
      RETURN
      END
C
C TESTS THE FUNCTIONING OF THE COMPARISON OF DOUBLES,
C LOGICAL FUNCTION EQUALD
      SUBROUTINE TESTD()
      LOGICAL DEQUAL
      DOUBLE PRECISION EXPECT
      DOUBLE PRECISION ACTUL1
      DOUBLE PRECISION ACTUL2
      DOUBLE PRECISION ACTUL3
      DOUBLE PRECISION EPS1
      DOUBLE PRECISION EPS2
      DATA EXPECT/1.0D0/, ACTUL1/1.0D0/, ACTUL2/1.00000001D0/,
     & ACTUL3/0.99999999D0/
      DATA EPS1/1.D-1/
C EPS2 SEEMS LIKE IT SHOULD BE 1.D-8, BUT DOUBLE PRECISION
C REPRESENTATION OF *ACTUL3* IS SLIGHTLY LESS THAN THE DECIMAL
C REPRESENTATION, SO THE TEST WILL FAIL IF 1.D-8 IS USED
      DATA EPS2/5.D-9/
      CALL ASRTT(DEQUAL(EXPECT, ACTUL1, EPS1))
      CALL ASRTT(DEQUAL(EXPECT, ACTUL1, EPS2))
      CALL ASRTT(DEQUAL(EXPECT, ACTUL2, EPS1))
      CALL ASRTF(DEQUAL(EXPECT, ACTUL2, EPS2))
      CALL ASRTT(DEQUAL(EXPECT, ACTUL3, EPS1))
      CALL ASRTF(DEQUAL(EXPECT, ACTUL3, EPS2))
      RETURN
      END



with the obvious implementations:


C
C DETERMINE IF TWO REALS ARE "ESSENTIALLY" EQUAL:
C ABS(EXPECT - ACTUAL) <= EPSILN
C
C INPUT EXPECT (REAL) - THE EXPECTED VALUE
C INPUT ACTUAL (REAL) - THE VALUE TO BE COMPARED TO EXPECT
C INPUT EPSILN (REAL) - THE TOLERANCE
C OUTPUT RESULT
      LOGICAL FUNCTION AEQUAL(EXPECT, ACTUAL, EPSILN)
      REAL EXPECT
      REAL ACTUAL
      REAL EPSILN
      AEQUAL = ABS(EXPECT - ACTUAL) .LE. EPSILN
      END
C
C DETERMINE IF TWO DOUBLES ARE "ESSENTIALLY" EQUAL:
C ABS(EXPECT - ACTUAL) <= EPSILN
C
C INPUT EXPECT (DOUBLE) - THE EXPECTED VALUE
C INPUT ACTUAL (DOUBLE) - THE VALUE TO BE COMPARED TO EXPECT
C INPUT EPSILN (DOUBLE) - THE TOLERANCE
C OUTPUT RESULT
      LOGICAL FUNCTION DEQUAL(EXPECT, ACTUAL, EPSILN)
      DOUBLE PRECISION EXPECT
      DOUBLE PRECISION ACTUAL
      DOUBLE PRECISION EPSILN
      DEQUAL = DABS(EXPECT - ACTUAL) .LE. EPSILN
      END



Finally, the COMPLEX data type has a similar implementation. The only difference is that EPSILN will also be of type COMPLEX. This allows the user to specify separate epsilon values for the real and imaginary parts of the comparison. The FORTRAN 66 spec specifies that both the real and imaginary parts of a COMPLEX are REALs. So, we come up with the test:

C
C TESTS THE FUNCTIONING OF THE COMPARISON OF COMPLEX VALUES,
C LOGICAL FUNCTION CEQUAL
C THIS IS SOMEWHAT INELEGANT, BUT COVERS ALL OF THE POSSIBLITIES
      SUBROUTINE TESTC
      LOGICAL CEQUAL
      COMPLEX EXPECT
      COMPLEX ACTUAL
      COMPLEX EPS1
      COMPLEX EPS2
      COMPLEX EPS3
      COMPLEX EPS4
      REAL ONE
      REAL ONEPLS
      REAL ONEMNS
      REAL EPSLN1
      REAL EPSLN2
      DATA ONE/1.0/, ONEPLS/1.00001/, ONEMNS/0.99999/
      DATA EPSLN1/1.E-1/, EPSLN2/1.E-7/
C
      EXPECT = CMPLX(ONE, ONE)
      ACTUAL = CMPLX(ONE, ONE)
      EPS1 = CMPLX(EPSLN1, EPSLN1)
      EPS2 = CMPLX(EPSLN2, EPSLN2)
      EPS4 = CMPLX(EPSLN1, EPSLN2)
      EPS3 = CMPLX(EPSLN2, EPSLN1)
C
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONEPLS, ONE)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONEMNS, ONE)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONE, ONEPLS)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONE, ONEMNS)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONEPLS, ONEPLS)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONEPLS, ONEMNS)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONEMNS, ONEPLS)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS4))
C
      ACTUAL = CMPLX(ONEMNS, ONEMNS)
      CALL ASRTT(CEQUAL(EXPECT, ACTUAL, EPS1))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS2))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS3))
      CALL ASRTF(CEQUAL(EXPECT, ACTUAL, EPS4))
      RETURN
      END


We could implement this inline (with calls to the ABS function). Instead, we again apply the DRY principle and call AEQUAL.


C
C DETERMINE IF TWO COMPLEX VALUES ARE "ESSENTIALLY" EQUAL:
C ABS(EXPECT - ACTUAL) <= EPSILN FOR BOTH THE REAL
C AND IMAGINARY PARTS
C
C INPUT EXPECT (COMPLEX) - THE EXPECTED VALUE
C INPUT ACTUAL (COMPLEX) - THE VALUE TO BE COMPARED TO EXPECT
C INPUT EPSILN (COMPLEX) - THE TOLERANCE
C OUTPUT RESULT
      LOGICAL FUNCTION
      CEQUAL(EXPECT, ACTUAL, EPSILN)
LOGICAL AEQUAL
      COMPLEX EXPECT
      COMPLEX ACTUAL
      COMPLEX EPSILN
      CEQUAL = AEQUAL(REAL(EXPECT), REAL(ACTUAL), REAL(EPSILN))
     & .AND. AEQUAL(AIMAG(EXPECT), AIMAG(ACTUAL), AIMAG(EPSILN))
      END


The next step is to write our xASRTE SUBROUTINEs in terms of xEQUAL and ASRTT:


C
C IASRTE
C ASSERTS THAT 2 INTEGERS ARE EQUAL
      SUBROUTINE IASRTE(EXPECT, ACTUAL)
      INTEGER EXPECT
      INTEGER ACTUAL
      LOGICAL RESULT
      RESULT = IEQUAL(EXPECT, ACTUAL)
      CALL ASRTT(RESULT)
      RETURN
      END
C
C LASRTE
C ASSERTS THAT 2 LOGICAL VARIABLES ARE EQUAL
      SUBROUTINE LASRTE(EXPECT, ACTUAL)
      LOGICAL EXPECT
      LOGICAL ACTUAL
      LOGICAL RESULT
      RESULT = LEQUAL(EXPECT, ACTUAL)
      CALL ASRTT(RESULT)
      RETURN
      END
C
C AASRTE
C ASSERTS THAT 2 REALS ARE ESSENTIALLY EQUAL
C ABS(EXPECT - ACTUAL) .LE. EPSILN
      SUBROUTINE AASRTE(EXPECT, ACTUAL, EPSILN)
      REAL EXPECT
      REAL ACTUAL
      REAL EPSILN
      LOGICAL RESULT
      RESULT = AEQUAL(EXPECT, ACTUAL)
      CALL ASRTT(RESULT)
      RETURN
      END
C
C DASRTE
C ASSERTS THAT 2 DOUBLE PRECISION VARIABLES ARE ESSENTIALLY EQUAL
C DABS(EXPECT - ACTUAL) .LE. EPSILN
      SUBROUTINE IASRTE(EXPECT, ACTUAL, EPSILN)
      DOUBLE PRECISION EXPECT
      DOUBLE PRECISION ACTUAL
      DOUBLE PRECISION EPSILN
      LOGICAL RESULT
      RESULT = DEQUAL(EXPECT, ACTUAL)
      CALL ASRTT(RESULT)
      RETURN
      END
C
C CASRTE
C ASSERTS THAT 2 DOUBLE PRECISION VARIABLES ARE ESSENTIALLY EQUAL
C ABS(REAL(EXPECT) - REAL(ACTUAL)) .LE. REAL(EPSILN)
C ABS(AIMAG(EXPECT) - AIMAG(ACTUAL)) .LE. AIMAG(EPSILN)
      SUBROUTINE CASRTE(EXPECT, ACTUAL, EPSILN)
      COMPLEX EXPECT
      COMPLEX ACTUAL
      COMPLEX EPSILN
      LOGICAL RESULT
      RESULT = CEQUAL(EXPECT, ACTUAL, EPSILN)
      CALL ASRTT(RESULT)
      RETURN
      END



Finally, we have a main program to run all of our tests:


      PROGRAM TCASES
      COMMON /OUTP/KOUT
      INTEGER KOUT
      KOUT = 6
      CALL TSTCHR()
      CALL TSTOUT()
      CALL TASRTT()
      CALL TASRTF()
      CALL TESTI()
      CALL TESTL()
      CALL TESTA()
      CALL TESTD()
      CALL TESTC()
      STOP
      END

Summary


In this installment, we've covered a lot of ground. We completed the basic version of our xUnit framework for FORTRAN 66.
We:
  • Wrote ASRTT by using RESOUT

  • Wrote ASRTF by using ASRTT

  • Created xEQUAL classes for each of the types
    • INTEGER
    • LOGICAL
    • REAL (accounting for the imprecision in floating point arithmetic)
    • DOUBLE PRECISION (again accounting for the imprecision in floating point arithmetic)
    • COMPLEX
    • (by using AEQUAL)

  • Wrote our xASRTE classes using xEQUAL and ASRTT

  • Did all of it using Test-driven development (TDD), so that we have a complete set of test cases to support refactoring.



Next time we will tighten up the code, cleaning up a a few remaining bits of duplication, and enhance the framework to include some useful feedback.

Friday, January 1, 2010

Creating A Fortran Unit Testing Framework (part 4): Outputting the Test Result

In this installment, we are going to implement the next set of functionality, wiring out the results of the tests. Recall that last time we wrote RESCHR, the routine that determined what to print out ("." or "F") depending upon the result of the test. Our next piece of functionality will be to write that character out, a routine we call RESOUT.
As with RESCHR, we will write RESOUT using test-driven development TDD). The first problem is how to figure out whether something is printing correctly. The simple answer is to just print it to the screen and look at it (and frankly is the very first approach I use). The problem with this approach is that it requires us to look at the output for every test every time we run them. This gets old in a hurry. To be really effective, unit tests must be quick and easy to run. The more effort required to run and evaluate them, the less often they will be used. So, we want the testing framework to automatically figure out whether or not the right thing is being printed out.
To make this automatic, we need to print to a place other than the default output and then compare what was written to our expectation. Most modern languages have some sort of internal representation of a string or file. In FORTRAN 66, we could use Hollerith fields encoded into integer arrays to achieve that effect, but it seems way too much effort of this particular job (we'll tackle that another day). Instead, I'm going to take the easy way out. I'll create a scratch file, write out to it, read the results back in, and compare them to my expectation. After a couple of false starts, here is the test that I came up with:



C
C TEST CASE FOR THE ROUTINE THAT PRINTS OUT THE
C CHARACTER BASED ON THE VALUE OF THE TEST
      SUBROUTINE TSTOUT()
      INTEGER I
      INTEGER PASS
      INTEGER FAIL
      INTEGER RESULT
      DIMENSION RESULT(2)
      DATA PASS/1H./, FAIL/1HF/
      OPEN(UNIT=10, STATUS='SCRATCH')
      CALL RESOUT(.TRUE.)
      CALL RESOUT(.FALSE.)
      REWIND(10)
      READ(10, 3000) (RESULT(I), I = 1, 2)
      IF (RESULT(1) .EQ. PASS) GO TO 100
          WRITE(6, 1000)
          GO TO 200
  100 CONTINUE
          WRITE(6, 2000)
  200 CONTINUE
      IF (RESULT(2) .EQ. FAIL) GO TO 300
          WRITE(6, 1000)
          GO TO 400
  300 CONTINUE
          WRITE(6, 2000)
  400 CONTINUE
 1000 FORMAT(1HF, $)
 2000 FORMAT(1H., $)
 3000 FORMAT(2A1)
      RETURN
      END



The next TDD step is to write our minimal code that will compile, and cause the test to fail:



C
C WRITES OUT '.' IF THE ARGUMENT IS .TRUE.
C 'F' OTHERWISE
      SUBROUTINE RESOUT(A)
      LOGICAL A
      WRITE(10, 100)
  100 FORMAT(2HXX)
      RETURN
      END



We run it, and see the expected output:


FF



Next, we remove the hard-coded unit (10), and take it an argument to the routine, and revise the code to make the tests pass (I actually do these in separate steps, but from the last post you get the idea of TDD, so I'm picking up the pace a bit):



C
C WRITES OUT '.' IF THE ARGUMENT IS .TRUE.
C 'F' OTHERWISE
      SUBROUTINE RESOUT(A, OUTPUT)
      LOGICAL A
      INTEGER OUTPUT
      IF (A) GO TO 200
          WRITE(OUTPUT, 100)
  100     FORMAT(1HF, $)
          GO TO 400
  200 CONTINUE
          WRITE(OUTPUT, 300)
  300     FORMAT(1H., $)
  400 CONTINUE
      RETURN
      END


And the output becomes:


..


indicating that both of the tests now pass.

Now, we can leave this version as it is, and some purists would say we should. However, those purists would be wrong. As it stands, RESOUT replicates the logic of RESCHR (choosing which character to print, based on the value of the LOGICAL argument A. The replication is a violation of the "Don't repeat yourself" rule (know as DRY, from The Pragmatic Programmer). There is no reason to have that code twice, so we modify RESOUT to make use of RESCHR routine that we developed last time. So, our final version of RESOUT becomes:



C
C WRITES OUT '.' IF THE ARGUMENT IS .TRUE.
C 'F' OTHERWISE
C
C LOGICAL A: WHETHER OR NOT THE TEST PASSED
C INTEGER OUTPUT: THE UNIT TO WRITE THE RESULT TO
C
      SUBROUTINE RESOUT(A, OUTPUT)
      LOGICAL A
      INTEGER OUTPUT
      INTEGER RESCHR
      INTEGER OUTCHR
      OUTCHR = RESCHR(A)
      WRITE(OUTPUT, 100) OUTCHR
  100 FORMAT(A1, $)
      RETURN
      END


We run this, and again the tests pass.
Finally, we tighten up the test code:



C
C TEST CASE FOR THE ROUTINE THAT PRINTS OUT THE
C CHARACTER BASED ON THE VALUE OF THE TEST
      SUBROUTINE TSTOUT()
      INTEGER I
      INTEGER OUTPUT
      INTEGER EXPECT
      INTEGER RESULT
      DIMENSION RESULT(2)
      DATA (EXPECT(I = 1 ,2))/1H., 1HF/, OUTPUT/10/
      OPEN(UNIT=OUTPUT, STATUS='SCRATCH')
      CALL RESOUT(.TRUE., OUTPUT)
      CALL RESOUT(.FALSE., OUTPUT)
      REWIND(10)
      READ(10, 100) (RESULT(I), I = 1, 2)
  100 FORMAT(2A1)
      DO 600 I = 1, 2
          IF (EXPECT(I) .EQ. RESULT(I)) go to 300
              WRITE(6, 200)
  200         FORMAT(1HF, $)
              GO TO 500
  300     CONTINUE
              WRITE(6, 400)
  400         FORMAT(1H., $)
  500     CONTINUE
  600 CONTINUE
      RETURN
      END