Beyond the Manual
The TDA User’s Manual stops at chapter 8, and development did not. The spatial commands, the file readers for formats that postdate the manual, and TDA’s online help are all in the program and in none of its documentation. This is what is there.
What happens inside
There is no reimplementation. TDA’s C, context-threaded and compiled into the package, runs unchanged; a wrapper writes the command file TDA would have been given, hands it over in memory, and reads the answer back. Every object therefore carries the run it came from.
Box 1One fit, and the run underneath it
d <- tda_rrdat()
fit <- tda_rate(Surv(TFP, DES) ~ EDU, d)
cat(fit$run$commands, sep = "\n")rdataframe;
edef(
ts = Tstart,
tf = Tfin,
org = Org,
des = Des,
);
rate(
mfmt = 24.16,
tfmt = 24.16,
xa (0,1) = EDU,
) = 2;The answer arrives on two channels. The console is what TDA printed –
tda_output() shows all of it, and
tda_diagnostics() the lines that are errors, warnings or
notes, which TDA otherwise leaves in the middle of its output for the
reader to find. The exports are the doubles behind the numbers TDA
printed: coefficients, covariance matrices, tables, and a tap on every
number the console formatter was handed. A wrapper takes its numbers
from the exports and uses the text only to find its way around, so
coef() holds the double TDA computed, not a number re-read
from a print format.
Box 2Printed, and computed
grep("EDU", fit$run$output, value = TRUE)[1] " 5 EDU 3 8 24.16 c5"
[2] " 2 1 0 1 A EDU 0.0150660299026325 0.0211574987658471 0.7120893669601657 0.5236"sprintf("%.17g", coef(fit)[["EDU"]])[1] "0.015066029902632524"names(fit$run$exports) [1] "episodes.table" "episodes.excluded" "rate.logLik.null" "ml.vcov"
[5] "ml.logLik" "ml.gradient" "ml.change.f" "ml.change.p"
[9] "rate.est" "rate.est.names" "ml.logLik.pair" "print.values"
[13] "print.labels" tda_output_audit() is the check that this holds: it
classifies every line TDA printed as boilerplate, an echo of the input,
a setting, a diagnostic, or a number that is in the object or in the
exports, and reports whatever is none of those. The package’s tests run
it over every example in this manual and over the 634 command files of
TDA’s own regression suite.
Box 3Where every line went
attr(tda_output_audit(fit), "summary")cls
blank boilerplate echo settings stored
20 13 9 11 18 Spatial data
TDA holds points, lines and polygons in a spatial structure of its own, and has about twenty commands over it. None of them appears in the manual.
tda_spatial() builds the structure from a data frame.
type = 1 is points; type = 3 groups
consecutive rows sharing an id into one polygon each.
Box 4A spatial structure, and what is in it
set.seed(7)
pts <- data.frame(id = 1:12, x = round(runif(12), 3), y = round(runif(12), 3))
s <- tda_spatial(pts)
tda_sd_info(s)Call: tda_sd_info(s)
Cases: 12
points lines polygons xmin ymin xmax ymax
12 0 0 0.07 0.009 0.989 0.997tda_sd_voronoi() is the Voronoi diagram. What
sdvd writes is itself a spatial data file – points for the
vertices, lines for the edges – so the result plots with the same
machinery the input does.
Box 5The Voronoi diagram of those points
v <- tda_sd_voronoi(s, opt = 2)
plot(v, xlim = c(-0.2, 1.2), ylim = c(-0.2, 1.2))
tda_sd_select() keeps what falls inside a rectangle and
tda_sd_clip() cuts to it. Both want the rectangle given;
TDA will not guess one.
Box 6Selecting and clipping to a region
tda_sd_select(s, rec = "0,0,0.5,0.5")$table8 records, lengths 2-5
[1] 1 1 1 2 1
[2] 0.398 0.096
[3] 2 1 1 3 1
[4] 0.116 0.453
[5] 3 1 1 4 1
[6] 0.070 0.085
[7] 4 1 1 10 1
[8] 0.459 0.295 tda_sd_clip(s, c(0, 0, 0.5, 0.5))$table8 records, lengths 2-5
[1] 1 1 1 2 1
[2] 0.398 0.096
[3] 2 1 1 3 1
[4] 0.116 0.453
[5] 3 1 1 4 1
[6] 0.070 0.085
[7] 4 1 1 10 1
[8] 0.459 0.295 Polygons are the same structure with type = 3.
tda_sd_data() gives the vertices back as an ordinary data
frame, which is the form most things in R want.
Box 7Two square polygons
poly <- data.frame(id = c(1, 1, 1, 1, 2, 2, 2, 2),
x = c(0, 2, 2, 0, 3, 5, 5, 3),
y = c(0, 0, 2, 2, 0, 0, 2, 2))
sp <- tda_spatial(poly, type = 3)
tda_sd_data(sp)$table V1 V2 V3
1 1 0 0
2 1 2 0
3 1 2 2
4 1 0 2
5 2 3 0
6 2 5 0
7 2 5 2
8 2 3 2Reading a shapefile
tda_read_shapefile() reads geometry and attaches the
.dbf attribute table as $attributes. The
manual documents SPSS and Stata files and nothing else; shapefiles, E00,
DCW, GSHHS, .dbf and BIFF .xls are all
readable and all undocumented. The file is the North Carolina SIDS
shapefile that the sf package ships, read here by TDA itself; these
boxes run when sf is installed.
Box 8North Carolina counties, read through TDA
nc <- tda_read_shapefile(system.file("shape/nc.shp", package = "sf"))
tda_sd_info(nc)Call: tda_sd_info(nc)
Cases: 108
points lines polygons xmin ymin xmax ymax
0 0 108 -84.32385 33.88199 -75.45698 36.58965head(nc$attributes[, c("AREA", "PERIMETER", "NAME", "FIPS")]) AREA PERIMETER NAME FIPS
1 0.114 1.442 Ashe 37009
2 0.061 1.231 Alleghany 37005
3 0.143 1.630 Surry 37171
4 0.070 2.968 Currituck 37053
5 0.153 2.206 Northampton 37131
6 0.097 1.670 Hertford 37091The geometry is TDA’s, so TDA’s plotting draws it.
Box 9The same counties, drawn by TDA
plot(nc)
tda_sd_enclosing() gives the bounding box or the convex
hull of everything in the structure, and tda_sd_data() the
vertices.
Box 10The bounding box, and the vertices
tda_sd_enclosing(nc)$table5 records, lengths 2-3
[1] 1 3 4
[2] -84.32385 33.88199
[3] -75.45698 33.88199
[4] -75.45698 36.58965
[5] -84.32385 36.58965 head(tda_sd_data(nc)$table) V1 V2 V3
1 1 -81.47276 36.23436
2 1 -81.54084 36.27251
3 1 -81.56198 36.27359
4 1 -81.63306 36.34069
5 1 -81.74107 36.39178
6 1 -81.69828 36.47178Some of these commands compare one object type against another in the same structure. A set of polygons on its own has no lines or points for them to find, and saying so is the right answer rather than an error.
Box 11Nothing to report, correctly
tda_sd_neighbours(nc)$text[1] "Creating a network from lines. Current memory: 407281 bytes."
[2] "Number of lines: 0"
[3] "Number of segments: 0"
[4] "Tolerance: 0.0001"
[5] "Maximal number of segments: 0"
[6] "Number of nodes: 0"
[7] "Number of edges: 0"
[8] "0 records written to: out.txt" TDA’s help
The program carries its help text for every command, and
tda_help() reads it. This is the fastest way to find what a
command’s options are actually called, since the R wrappers deliberately
do not expose all of them by name.
Box 12tda_help
cat(head(tda_help("ploth"), 20), sep = "\n")ploth: The ploth command can be used to plot a histogram. Syntax is:
ploth(
x=..., intervals, required
w=..., variable for weights
dscal=..., scaling factor, def. 1.0
s=..., s=0 (default) classes are left closed
s=1 classes are left open
ns=..., ns=1, don't plot vertical lines, def. ns=0
lt=..., line type, def. 1 (solid line)
lw=..., line width, def. 0.2 (mm)
gs=..., grey scale value, def. 1 (white)
nc=..., nc=1 if no clipping, def. 0
df=..., print table to output file
fmt=..., print format for output file, def. 10.4
) = variable;
There must be exactly one variable name on the right-hand side. Also
one has to use the x parameter to specify a set of intervals (classes)
for the histogram. The syntax is
x=x1,x2,...,xn, or x=x1(d)x2,
x1 must be not greater than the minimum of the values of the variable,
and xn must not smaller than the maximum of the values of the variable.
The command requires a valid setup for plots and uses the currently
defined coordinate system. The heights of the histogram are calculated
in such a way that the areas equal the frequencies. A corresponding
frequency table can be requested with the df parameter. The output file
will then contain n - 1 lines with the following entries:
1) begin of class
2) end of class
3) absolute frequency
4) relative frequency
5) height of histogram
See also: atab
ploth: The ploth command can be used to plot a histogram. Syntax is:
ploth(
x=..., intervals, required
w=..., variable for weights
dscal=..., scaling factor, def. 1.0
s=..., s=0 (default) classes are left closed
s=1 classes are left open
ns=..., ns=1, don't plot vertical lines, def. ns=0
lt=..., line type, def. 1 (solid line)
lw=..., line width, def. 0.2 (mm)
gs=..., grey scale value, def. 1 (white)
nc=..., nc=1 if no clipping, def. 0
df=..., print table to output file
fmt=..., print format for output file, def. 10.4
) = variable;
There must be exactly one variable name on the right-hand side. Also
one has to use the x parameter to specify a set of intervals (classes)
for the histogram. The syntax isRunning TDA directly
tda_run() takes command-file text and runs it, which is
the escape hatch for anything without a wrapper.
tda_block() builds one command from named arguments so the
syntax does not have to be written by hand.
Box 13A command file, from R
d <- data.frame(X = c(1, 2, 3, 4), Y = c(2, 4, 6, 9))
res <- tda_run(c(tda_nvar(d), "dstat;"), data = d)
cat(grep("^ *(Variable|X|Y)", res$output, value = TRUE), sep = "\n")Variable Minimum Maximum Mean Std.Dev. Sum of values
X 1.0000 4.0000 2.5000 1.2910 10.0000
Y 2.0000 9.0000 5.2500 2.9861 21.0000Every wrapper returns its $run, so the command file a
call generated can always be read back. This is worth knowing when a
result looks wrong: the question is usually what TDA was asked, not what
it computed.
Box 14What a wrapper actually sent
The expression language
TDA’s operators are available without a data set through
tda_evalf(), and its distribution functions are the ones
the manual lists in 5.2.5 – nd, ndf,
td, chd, fd, mr,
poisson, negbin. They agree with R’s, which is
checked by a differential fuzzer on every run of the test suite, so
either side gives the same numbers.
Box 15The same quantity, twice
Character data
TDA holds fixed-width string variables, and a character or factor
column reaches it: tda_nvar() declares it as
str(n,m) and tda_strings() brings it back. The
four string operators – strlen, strsp,
strv, strvp – work on it.
Box 16Strings through TDA
d2 <- data.frame(ID = 1:3, Name = c("Ashe", "Alleghany", "Surry"),
Area = c(0.114, 0.061, 0.143))
res <- tda_run(c(tda_nvar(d2), "pdata() = out.txt;"), data = d2)
tda_strings(res)$Name
[1] "Ashe" "Alleghany" "Surry" Graphs
Chapter 7 of the manual covers relational data, and the 31
tda_g_* commands implement it: components, cliques, cut
points, spanning trees, cycles, flows, shortest paths, and the
social-network measures of 7.6. They are documented but easy to miss, so
one example, on a small graph.
Box 17A graph, and some of what TDA asks of it
e <- data.frame(from = c(1, 1, 2, 3, 4, 5, 5),
to = c(2, 3, 3, 4, 5, 6, 7))
g <- tda_graph(e, directed = FALSE)
tda_g_degrees(g) index node in_degree out_degree loops degree
1 1 1 2 2 0 2
2 2 2 2 2 0 2
3 3 3 3 3 0 3
4 4 4 2 2 0 2
5 5 5 3 3 0 3
6 6 6 1 1 0 1
7 7 7 1 1 0 1tda_g_components(g)Call: tda_g_components(g)
Cases: 7
component n i node_i
1 7 1 1
1 7 2 2
1 7 3 3
1 7 4 4
1 7 5 5
1 7 6 6
1 7 7 7tda_g_cutpoints(g)Call: tda_g_cutpoints(g)
Cases: 7
V1 V2 V3 V4 V5 V6
1 7 3 3 4 5Box 18The same graph, drawn by TDA
plot(g)