phasePortrait {viscomplexr}

Create phase portraits of complex functions

Description

phasePortrait makes phase portraits of functions in the complex number plane. It uses a technique often (but not quite correctly) called domain coloring (https://en.wikipedia.org/wiki/Domain_coloring). While many varieties of this technique exist, this book relates closely to the standards proposed by E. Wegert in his book Visual Complex Functions (Wegert 2012). In a nutshell, the argument (Arg) of any complex function value is displayed as a color from the chromatic circle. The fundamental colors red, green, and blue relate to the arguments (angles) of 0, 2/3pi, and 4/3pi (with smooth color transitions in between), respectively. Options for displaying the modulus (Mod) of the complex values and additional reference lines for the argument are available. This function is designed for being used inside the framework of R base graphics. It makes use of parallel computing, and depending on the desired resolution it may create extensive sets of large temporary files (see Details and Examples).

Usage

phasePortrait(
  FUN,
  moreArgs = NULL,
  xlim,
  ylim,
  invertFlip = FALSE,
  res = 150,
  blockSizePx = 2250000,
  tempDir = NULL,
  nCores = max(1, parallel::detectCores() - 1),
  pType = "pma",
  pi2Div = 9,
  logBase = exp(2 * pi/pi2Div),
  argOffset = 0,
  darkestShade = 0.1,
  lambda = 7,
  gamma = 0.9,
  stdSaturation = 0.8,
  hsvNaN = c(0, 0, 0.5),
  asp = 1,
  deleteTempFiles = TRUE,
  noScreenDevice = FALSE,
  autoDereg = FALSE,
  verbose = TRUE,
  ...
)

Arguments

FUN

The function to be visualized. There are two possibilities to provide it, a quoted character string, or a function object.

Quoted character string

A function can be provided as a quoted character string containing an expression R can interpret as a function of a complex number z. Examples: "sin(z)", "(z^2 - 1i)/(tan(z))", "1/4*z^2 - 10*z/(z^4+4)". Names of functions known in your R session can be used in a standalone way, without mentioning z, e.g. "sin", "tan", "sqrt". Obviously, this also works for functions you defined yourself, e.g. "myIncredibleFunction" would be valid if you coded a function with this name before. This is especially useful for functions which require additional parameters beside the complex number they are supposed to calculate with. Such arguments can be provided via the parameter moreArgs. One-liner expressions provided as strings are also compatible with moreArgs (see examples).

While it is not the way we recommend for most purposes, you can even define more complicated functions of your own as character strings. In this case, you need to use vapply as a wrapper for your actual function (see Details, and Examples). Such constructions allow to provide additional input variables as a part of the character string by using the vapply-mechanism (see Details and Examples). The helper function vector2String) can be useful for that matter. However, the parameter moreArgs is not applicable in this context. Probably, the most useful application of the function-as-string concept is when the user defined function, possibly including values for additional arguments, is to be pasted together at runtime.

Function object

It is also possible to directly provide function objects to FUN. This can be any function known to R in the current session. Simply put, for functions like sin, tan, cos, and sqrt you do not even have to quote their names when passing them to phasePortrait. Same applies to functions you defined yourself. It is also possible to hand over an anonymous function to FUN when calling phasePortrait. In all these cases, the parameter moreArgs can be used for providing additional arguments to FUN. In general, providing a function as an object, and using moreArgs in case additional arguments are required, is what we recommend for user-defined functions.

When executing phasePortrait, FUN is first evaluated with match.fun. If this is not successful, an attempt to interpret FUN as an expression will be made. If this fails, phasePortrait terminates with an error.

moreArgs

A named list of other arguments to FUN. The names must match the names of the arguments in FUN's definition.

xlim

The x limits (x1, x2) of the plot as a two-element numeric vector. Follows exactly the same definition as in plot.default. Here, xlim has to be interpreted as the plot limits on the real axis.

ylim

The y limits of the plot (y1, y2) to be used in the same way as xlim. Evidently, ylim indicates the plot limits on the imaginary axis.

invertFlip

If TRUE, the function is mapped over a z plane, which has been transformed to 1/z * exp(1i*pi). This is the projection required to plot the north Riemann hemisphere in the way proposed by Wegert (2012), p. 41. Defaults to FALSE. If this option is chosen, the numbers at the axis ticks have another meaning than in the normal case. Along the real axis, they represent the real part of 1/z, and along the imaginary axis, they represent the imaginary part of 1/z. Thus, if you want annotation, you should choose appropriate axis labels like xlab = Re(1/z), and ylab = Im(1/z).

res

Desired resolution of the plot in dots per inch (dpi). Default is 150 dpi. All other things being equal, res has a strong influence on computing times (double res means fourfold number of pixels to compute). A good approach could be to make a plot with low resolution (e.g. the default 150 dpi) first, adjust whatever required, and plot into a graphics file with high resolution after that.

blockSizePx

Number of pixels and corresponding complex values to be processed at the same time (see Details). Default is 2250000. This value gave good performance on older systems as well as on a high-end gaming machine, but some tweaking for your individual system might even improve things.

tempDir

NULL or a character string, specifying the name of the directory where the temporary files written by phasePortrait are stored. Default is NULL, which makes phasePortrait use the current R session's temporary directory. Note that if you specify another directory, it will be created if it does not exist already. Even though the temporary files are deleted after completing a phase portrait (unless the user specifies deleteTempFiles = FALSE, see below), the directory will remain alive even if has been created by phasePortrait.

nCores

Number of processor cores to be used in the parallel computing tasks. Defaults to the maximum number of cores available minus 1. Any number between 1 (serial computation) and the maximum number of cores available as indicated by parallel::detectCores() is accepted. If nCores is set to a value greater than the available number of cores, the function will use one core less than available.

pType

One of the four options for plotting, "p", "pa", "pm", and "pma" as a character string. Defaults to "pma". Option "p" produces a mere phase plot, i.e. contains only colors for the complex numbers' arguments, but no reference lines at all. the option "pa" introduces shading zones that emphasize the arguments. These zones each cover an angle defined by 2*pi/pi2Div, where p2Div is another parameter of this function (see there). These zones are shaded darkest at the lowest angle (counter clockwise). Option "pm" displays the modulus by indicating zones, where the moduli at the higher edge of each zone are in a constant ratio with the moduli at the lower edge of the zone. Default is a ratio of almost exactly 2 (see parameter logBase) for details. At the lower edge, color saturation is lowest and highest at the higher edge (see parameters darkestShade, and stdSaturation). Option "pma" (default) includes both shading schemes.

pi2Div

Angle distance for the argument reference zones added for pType = "pma" or pType = "pa". The value has to be given as an integer (reasonably) fraction of 2*pi (i.e. 360 degrees). The default is 9; thus, reference zones are delineated by default in distances of 2*pi/9, i.e. (40 degrees), starting with 0, i.e. the color red if not defined otherwise with the parameter argOffset. In contrast to the borders delimiting the modulus zones, the borders of the reference zones for the argument always follow the same color (by definition).

logBase

Modulus ratio between the edges of the modulus reference zones in pType "pm" and "pma". As recommended by Wegert (2012), the default setting is logBase = exp(2*pi/pi2Div). This relation between the parameters logBase and pi2Div ensures an analogue scaling of the modulus and argument reference zones (see Details). Conveniently, for the default pi2Div = 9, we obtain logBase == 2.0099..., which is very close to 2. Thus, the modulus at the higher edge of a given zone is almost exactly two times the value at the lower edge.

argOffset

The (complex number) argument in radians counterclockwise, at which the argument reference zones are fixed. Default is 0, i.e. all argument reference zones align to the center of the red area.

darkestShade

Darkest possible shading of modulus and angle reference zones for pType "pm" and "pma". It corresponds to the value "v" in the hsv color model. darkestShade = 0 means no brightness at all, i.e. black, while darkestShade = 1 indicates maximum brightness. Defaults to 0.1, i.e. very dark, but hue still discernible.

lambda

Parameter steering the shading interpolation between the higher and the lower edges of the the modulus and argument reference zones in pType "pm" and "pm". Should be > 0, default and reference is lambda = 7. Values < 7 increase the contrast at the zone borders, values > 7 weaken the contrast.

gamma

Parameter for adjusting the combined shading of modulus and argument reference zones in pType "pma". Should be in the interval [0, 1]. Default is 0.9. The higher the value, the more the smaller of both shading values will dominate the outcome and vice versa.

stdSaturation

Saturation value for unshaded hues which applies to the whole plot in pType "p" and to the (almost) unshaded zones in pType "pm" and "p". This corresponds to the value "s" in the hsv color model. Must be between 0 and 1, where 1 indicates full saturation and 0 indicates a neutral grey. Defaults to 0.8.

hsvNaN

hsv coded color for being used in areas where the function to be plotted is not defined. Must be given as a numeric vector with containing the values h, s, and v in this order. Defaults to c(0, 0, 0.5) which is a neutral grey.

asp

Aspect ratio y/x as defined in plot.window. Default is 1, ensuring an accurate representation of distances between points on the screen.

deleteTempFiles

If TRUE (default), all temporary files are deleted after the plot is completed. Set it on FALSE only, if you know exactly what you are doing - the temporary files can occupy large amounts of hard disk space (see details).

noScreenDevice

Suppresses any graphical output if TRUE. This is only intended for test purposes and makes probably only sense together with deleteTempFiles == FALSE. For dimensioning purposes, phasePortrait will use a 1 x 1 inch pseudo graphics device in this case. The default for this parameter is FALSE, and you should change it only if you really know what you are doing.

autoDereg

if TRUE, automatically register sequential backend after the phase portrait is completed. Default is FALSE, because registering a parallel backend can be time consuming. Thus, if you want to make several phase portraits in succession, you should set autoDereg only for the last one, or simply type foreach::registerDoSEQ after you are done. In any case, you don't want to have an unused parallel backend lying about.

verbose

if TRUE (default), phasePortrait will continuously write progress messages to the console. This is convenient for normal purposes, as calculating larger phase portraits in higher resolution may take several minutes. The setting verbose = FALSE, will suppress any output to the console.

...

All parameters accepted by the plot.default function.

Details

This function is intended to be used inside the framework of R base graphics. It plots into the active open graphics device where it will display the phase plot of a user defined function as a raster image. If no graphics device is open when called, the function will plot into the default graphics device. This principle allows to utilize the full functionality of R base graphics. All graphics parameters (par) can be freely set and the function phasePortrait accepts all parameters that can be passed to the plot.default function. This allows all kinds of plots - from scientific representations with annotated axes and auxiliary lines, notation, etc. to poster-like artistic pictures.

Mode of operation

After being called, phasePortrait gets the size in inch of the plot region of the graphics device it is plotting into. With the parameter res which is the desired plot resolution in dpi, the horizontal and vertical number of pixels is known. As xlim and ylim are provided by the user, each pixel can be attributed a complex number z from the complex plane. In that way a two-dimensional array is built, where each cell represents a point of the complex plane, containing the corresponding complex number z. This array is set up in horizontal strips (i.e. split along the imaginary axis), each strip containing approximately blockSizePx pixels. In a parallel computing loop, the strips are constructed, saved as temporary files and immediately deleted from the RAM in order to avoid memory overflow. After that, the strips are sequentially loaded and subdivided into a number of chunks that corresponds to the number of registered parallel workers (parameter nCores). By parallely processing each chunk, the function f(z) defined by the user in the argument FUN is applied to each cell of the strip. This results in an array of function values that has exactly the same size as the original strip. The new array is saved as a temporary file, the RAM is cleared, and the next strip is loaded. This continues until all strips are processed. In a similar way, all strips containing the function values are loaded sequentially, and in a parallel process the complex values are translated into colors which are stored in a raster object. While the strips are deleted from the RAM after processing, the color values obtained from each new strip are appended to the color raster. After all strips are processed, the raster is plotted into the plot region of the graphics device. If not explicitly defined otherwise by the user, all temporary files are deleted after that.

Temporary file system

By default, the above-mentioned temporary files are deleted after use. This will not happen, if the parameter deleteTempFiles is set to FALSE or if phasePortrait does not terminate properly. In both cases, you will find the files in the directory specified with the parameter tempDir. These files are .RData files, each one contains a two-dimensional array of complex numbers. The file names follow a strict convention, see the following examples:

0001zmat2238046385.RData
0001wmat2238046385.RData

Both names begin with '0001', indicating that the array's top line is the first line of the whole clipping of the complex number plane where the phase portrait relates to. The array which follows below can e.g. begin with a number like '0470', indicating that its first line is line number 470 of the whole clipping. The number of digits for these line numbers is not fixed. It is determined by the greatest number required. Numbers with less digits are zero-padded. The second part of the file name is either zmat or wmat. The former indicates an array whose cells contain untransformed numbers of the complex number plane. The latter contains the values obtained from applying the function of interest to the first array. Thus, cells at the same position in both arrays exactly relate to each other. The third part of the file names is a ten-digit integer. This is a random number which all temporary files stemming from the same call of phasePortrait have in common. This guarantees that no temporary files will be confounded by the function, even if undeleted temporary files from previous runs are still present.

HSV color model

For color-coding the argument of a complex number, phasePortrait uses the hsv (hue, saturation, value) color model. Hereby, the argument is mapped to a position on the chromatic circle, where the fundamental colors red, green, and blue relate to the arguments (angles) of 0, 2/3pi, and 4/3pi, respectively. This affects only the hue component of the color model. The value component is used for shading modulus and/or argument zones. The saturation component for all colors can be defined with the parameter stdSaturation.

Zone definitions and shading

In addition to displaying colors for the arguments of complex numbers, zones for the modulus and/or the argument are shaded for pType other than "p". The modulus zones are defined in a way that each zone covers moduli whose logarithms to the base logBase have the same integer part. Thus, from the lower edge of one modulus zone to its upper edge, the modulus multiplies with the value of logBase. The shading of a modulus zone depends on the fractional parts x of the above-mentioned logarithms, which cover the interval [0, 1[. This translates into the value component v of the hsv color model as follows:

v = darkestShade + (1 - darkestShade) * x^(1/lambda)

where darkestShade and lambda are parameters that can be defined by the user. Modifying the parameters lambda and darkestShade is useful for adjusting contrasts in the phase portraits. The argument zone definition is somewhat simpler: Each zone covers an angle domain of 2*pi / pi2Div, the "zero reference" for all zones being argOffset. The angle domain of one zone is linearly mapped to a value x from the range [0, 1[. The value component of the color to be displayed is calculated as a function of x with the same equation as shown above. In case the user has chosen pType = "pma", x-values xMod and xArg are calculated separately for the modulus and the argument, respectively. They are transformed into preliminary v-values as follows:

vMod = xMod^(1/lambda) and vArg = xArg^(1/lambda)

From these, the final v value results as

v = darkestShade + (1-darkestShade) * (gamma * vMod * vArg + (1-gamma) * (1 - (1-vMod) * (1-vArg)))

The parameter gamma (range [0, 1]) determines they way how vMod and vArg are combined. The closer gamma is to one, the more the smaller of both values will dominate the outcome and vice versa.

Defining more complicated functions as strings with vapply

You might want to write and use functions which require more code than a single statement like (z-3)^2+1i*z. As mentioned in the description of the parameter FUN, we recommend to define such functions as separate objects and hand them over as such. There might be, however, cases, where it is more convenient, to define a function as a single long string, and pass this string to FUN. In order to make this work, vapply should be be used for wrapping the actual code of the function. This is probably not the use of vapply intended by its developers, but it works nicely and performs well. The character string has to have the following structure "vapply(z, function(z, other arguments if required) {define function code in here}, define other arguments here, FUN.VALUE = complex(1))". See examples.

References

Wegert E (2012). Visual Complex Functions. An Introduction with Phase Portraits. Springer, Basel Heidelberg New York Dordrecht London. ISBN 978-3-0348-0179-9.

Examples

# Map the complex plane on itself

# x11(width = 8, height = 8)   # Screen device commented out
                               # due to CRAN test requirements.
                               # Use it when trying this example
phasePortrait("z", xlim = c(-2, 2), ylim = c(-2, 2),
              xlab = "real", ylab = "imaginary",
              verbose = FALSE, # Suppress progress messages
              nCores = 2)      # Max. two cores allowed on CRAN
                               # not a limit for your own use



# A rational function

# x11(width = 10, height = 8) # Screen device commented out
                              # due to CRAN test requirements.
                              # Use it when trying this example
phasePortrait("(2-z)^2*(-1i+z)^3*(4-3i-z)/((2+2i+z)^4)",
              xlim = c(-8, 8), ylim = c(-6.3, 4.3),
              xlab = "real", ylab = "imaginary",
              nCores = 2)     # Max. two cores allowed on CRAN
                              # not a limit for your own use
  



# Different pType options by example of the sine function.
# Note the different equivalent definitions of the sine
# function in the calls to phasePortrait

# x11(width = 9, height = 9) # Screen device commented out
                             # due to CRAN test requirements.
                             # Use it when trying this example
op <- par(mfrow = c(2, 2), mar = c(2.1, 2.1, 2.1, 2.1))
phasePortrait("sin(z)", xlim = c(-pi, pi), ylim = c(-pi, pi),
              pType = "p",   main = "pType = 'p'",   axes = FALSE,
              nCores = 2) # Max. two cores on CRAN, not a limit for your use
  
phasePortrait("sin(z)", xlim = c(-pi, pi), ylim = c(-pi, pi),
              pType = "pm",  main = "pType = 'pm'",  axes = FALSE,
              nCores = 2)
  
phasePortrait("sin",    xlim = c(-pi, pi), ylim = c(-pi, pi),
              pType = "pa",  main = "pType = 'pa'",  axes = FALSE,
              nCores = 2)
  
phasePortrait(sin,      xlim = c(-pi, pi), ylim = c(-pi, pi),
              pType = "pma", main = "pType = 'pma'", axes = FALSE,
              nCores = 2)
  
par(op)


# I called this one 'nuclear fusion'

# x11(width = 16/9*8, height = 8) # Screen device commented out
                                  # due to CRAN test requirements.
                                  # Use it when trying this example

op <- par(mar = c(0, 0, 0, 0), omi = c(0.2, 0.2, 0.2, 0.2), bg = "black")
phasePortrait("cos((z + 1/z)/(1i/2 * (z-1)^10))",
              xlim = 16/9*c(-2, 2), ylim = c(-2, 2),
              axes = FALSE, xaxs = "i", yaxs = "i",
              nCores = 2) # Max. two cores allowed on CRAN
                          # not a limit for your own use
  
par(op)


# Passing function objects to phasePortrait:
# Two mathematical celebrities - Riemann's zeta function
# and the gamma function, both from the package pracma.
# R's built-in gamma is not useful, as it does not work
# with complex input values.

if(requireNamespace("pracma", quietly = TRUE)) {
# x11(width = 16, height = 8) # Screen device commented out
                              # due to CRAN test requirements.
                              # Use it when trying this example
op <- par(mfrow = c(1, 2))
phasePortrait(pracma::zeta,  xlim = c(-35, 15), ylim = c(-25, 25),
              xlab = "real", ylab = "imaginary",
              main = expression(zeta(z)), cex.main = 2,
              nCores = 2) # Max. two cores on CRAN, not a limit for your use
  
phasePortrait(pracma::gammaz, xlim = c(-10, 10), ylim = c(-10, 10),
              xlab = "real", ylab = "imaginary",
              main = expression(Gamma(z)), cex.main = 2,
              nCores = 2) # Max. two cores allowed on CRAN
                          # not a limit for your own use
  
}



# Using vapply for defining a whole function as a string.
# This is a Blaschke product with a sequence a of twenty numbers.
# See the documentation of the function vector2String for a more
# convenient space-saving definition of a.
# But note that a C++ version of the Blaschke product is available
# in this package (function blaschkeProd()).

# x11(width = 10, height = 8) # Screen device commented out
                              # due to CRAN test requirements.
                              # Use it when trying this example
phasePortrait("vapply(z, function(z, a) {
                fct <- ifelse(abs(a) != 0,
                  abs(a)/a * (a-z)/(1-Conj(a)*z), z)
                return(prod(fct))
              },
              a = c(0.12152611+0.06171533i,  0.53730315+0.32797530i,
                    0.35269601-0.53259644i, -0.57862039+0.33328986i,
                   -0.94623221+0.06869166i, -0.02392968-0.21993132i,
                    0.04060671+0.05644165i,  0.15534449-0.14559097i,
                    0.32884452-0.19524764i,  0.58631745+0.05218419i,
                    0.02562213+0.36822933i, -0.80418478+0.58621875i,
                   -0.15296208-0.94175193i, -0.02942663+0.38039250i,
                   -0.35184130-0.24438324i, -0.09048155+0.18131963i,
                    0.63791697+0.47284679i,  0.25651928-0.46341192i,
                    0.04353117-0.73472528i, -0.04606189+0.76068461i),
              FUN.VALUE = complex(1))",
              pType = "p",
              xlim = c(-4, 2), ylim = c(-2, 2),
              xlab = "real", ylab = "imaginary",
              nCores = 2) # Max. two cores allowed on CRAN
                          # not a limit for your own use
  



# Much more elegant: Define the function outside.
# Here comes a Blaschke product with 200 random points.

# define function for calculating blaschke products, even
# possible as a one-liner
blaschke <- function(z, a) {
  return(prod(ifelse(abs(a) != 0, abs(a)/a * (a-z)/(1-Conj(a)*z), z)))
}
# define 200 random numbers inside the unit circle
n <- 200
a <- complex(modulus = runif(n), argument = runif(n)*2*pi)
# Plot it
# x11(width = 10, height = 8) # Screen device commented out
                              # due to CRAN test requirements.
                              # Use it when trying this example
phasePortrait(blaschke,
  moreArgs = list(a = a),
  pType = "p",
  xlim = c(-2.5, 2.5), ylim = c(-1.7, 1.7),
  xlab = "real", ylab = "imaginary",
  nCores = 2) # Max. two cores allowed on CRAN
              # not a limit for your own use
  



# A hybrid solution: A one-liner expression given as a character string
# can be provided additional arguments with moreArgs

n <- 73
a <- complex(modulus = runif(n), argument = runif(n)*2*pi)
# x11(width = 10, height = 8) # Screen device commented out
                              # due to CRAN test requirements.
                              # Use it when trying this example
phasePortrait("prod(ifelse(abs(a) != 0,
  abs(a)/a * (a-z)/(1-Conj(a)*z), z))",
  moreArgs = list(a = a),
  pType = "p",
  xlim = c(-2.5, 2.5), ylim = c(-1.7, 1.7),
  xlab = "real", ylab = "imaginary",
  nCores = 1) # Max. two cores allowed on CRAN
              # not a limit for your own use
  



# Note the difference in performance when using the C++ defined
# function blaschkeProd() provided in this package

n <- 73
a <- complex(modulus = runif(n), argument = runif(n)*2*pi)
# Plot it
# x11(width = 10, height = 8) # Screen device commented out
                              # due to CRAN test requirements.
                              # Use it when trying this example
phasePortrait(blaschkeProd,
  moreArgs = list(a = a),
  pType = "p",
  xlim = c(-2.5, 2.5), ylim = c(-1.7, 1.7),
  xlab = "real", ylab = "imaginary",
  nCores = 1) # Max. two cores allowed on CRAN
              # not a limit for your own use
  



# Interesting reunion with Benoit Mandelbrot.
# The function mandelbrot() is part of this package (defined
# in C++ for performance)

# x11(width = 11.7, height = 9/16*11.7) # Screen device commented out
                                        # due to CRAN test requirements.
                                        # Use it when trying this example
op <- par(mar = c(0, 0, 0, 0), bg = "black")
phasePortrait(mandelbrot,
              moreArgs = list(itDepth = 100),
              xlim = c(-0.847, -0.403), ylim = c(0.25, 0.50),
              axes = TRUE, pType = "pma",
              hsvNaN = c(0, 0, 0), xaxs = "i", yaxs = "i",
              nCores = 1) # Max. two cores allowed on CRAN
                          # not a limit for your own use
par(op)
  



# Here comes a Julia set.
# The function juliaNormal() is part of this package (defined
# in C++ for performance)

# x11(width = 11.7, height = 9/16*11.7) # Screen device commented out
                                        # due to CRAN test requirements.
                                        # Use it when trying this example
op <- par(mar = c(0, 0, 0, 0), bg = "black")
phasePortrait(juliaNormal,
  moreArgs = list(c = -0.09 - 0.649i, R_esc = 2),
  xlim = c(-2, 2),
  ylim = c(-1.3, 1.3),
  hsvNaN = c(0, 0, 0),
  nCores = 1) # Max. two cores allowed on CRAN
              # not a limit for your own use
par(op)