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).
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, ... )
| FUN | The function to be visualized. There are two possibilities to provide it, a quoted character string, or a function object.
When executing |
|---|---|
| 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
|
| ylim | The y limits of the plot (y1, y2) to be used in the same way as
|
| invertFlip | If |
| res | Desired resolution of the plot in dots per inch (dpi). Default is
150 dpi. All other things being equal, |
| 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 |
|
| 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 |
| 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
|
| pi2Div | Angle distance for the argument reference zones added for
|
| logBase | Modulus ratio between the edges of the modulus reference zones
in |
| 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 |
| lambda | Parameter steering the shading interpolation between the higher
and the lower edges of the the modulus and argument reference zones in
|
| gamma | Parameter for adjusting the combined shading of modulus and
argument reference zones in |
| stdSaturation | Saturation value for unshaded hues which applies to the
whole plot in |
| hsvNaN |
|
| asp | Aspect ratio y/x as defined in |
| 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
|
| 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 |
| verbose | if TRUE (default), |
| ... | All parameters accepted by the |
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.
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.
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.
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.
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.
vapplyYou 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.
Wegert E (2012). Visual Complex Functions. An Introduction with Phase Portraits. Springer, Basel Heidelberg New York Dordrecht London. ISBN 978-3-0348-0179-9.
# 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 # \donttest{ # 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 # \donttest{ # 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 # \donttest{ 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. # \donttest{ 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()). # \donttest{ # 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. # \donttest{ # 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 # \donttest{ 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 # \donttest{ 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) # \donttest{ # 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) # \donttest{ # 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) # }