tests/testthat/test-accuracy.R

In baptiste/rcppfaddeeva: 'Rcpp' Bindings for the 'Faddeeva' Package

library(RcppFaddeeva)
library(testthat)
  
## parse stuff from c code
# library(dplyr)
# tmp  =  readLines("tmp.txt")
# tmp  =  tmp %>% gsub("C","c", .) %>% gsub("),","i),", .) %>% 
#   gsub("NaNi", "NaN", .) %>% gsub("Infi", "Inf", .)
# eval(parse(text=paste0("dput(sapply(list(",paste(tmp, collapse='', sep=""), "), sum))")))


context("Checking the accuracy")

  
tests <- list(
ref_Faddeeva_w = list(
z  =  c(624.2-0.26123i, -0.4+3i, 0.6+2i, -1+1i, -1-9i, -1+9i, -2.34545e-08+1.1234e+00i, 
       -3+5.1i, -53+30.1i, 0+0.12345i, 11+1i, -22-2i, 9-28i, 21-33i, 
       1e+05+1e+05i, 1e+14+1e+14i, -3001-1000i, 1e+160-1e+159i, -6.01+0.01i, 
       -0.7-0.7i, 26.1178+4540.90961097249i, 8e+06+3e+06i, -20-19.8081i, 
       1e-16-1.1e-16i, 2.3e-08+1.3e-08i, 6.3e+00-1e-13i, 6.3+0i, 0+6.3i, 
       0+16.3i, 9+0i, 6.01+0.11i, 8.01e+00+1.01e-10i, 28.01+0i, 10.01+0i, 
       10.01-0i, 1.001e+01+9.9e-11i, 1.001e+01-9.9e-11i, 0+7.01i, -1+7.01i, 
       5.99+7.01i, 1+0i, 55+0i, -0.1+0i, 1e-20+0i, 0+5e-14i, 0+51i, 
       Inf+0i, -Inf+0i, Inf+0i, -Inf+0i, Inf+0i, complex(real=NaN, imaginary=0), 
       complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
       complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
       complex(real=NaN, imaginary=0)),
value  =  c(-3.78270245519e-07+9.03861276433172e-04i, 0.176490622700482-0.021465505394685i, 
  0.241025071577269+0.060875796634281i, 0.304744205256913-0.208218938202832i, 
  7.31713106897238e+34+8.3218734997144e+34i, 0.0615698507236324-0.0067600578371657i, 
  3.96079300769987e-01-5.593152e-09i, 0.0821719922673945-0.0470129108764361i, 
  0.00457246000350282-0.00804900791411692i, 0.874634285960805+0i, 
  0.0046819016496544+0.0510735563901306i, -0.0023193175200188-0.0254600547397316i, 
  9.11463368405637e+304+3.97101807145263e+305i, -4.49272078577156e+281-2.80195912134231e+281i, 
  2.8209479178093e-06+2.82094791766826e-06i, 2.82094791773878e-15+2.82094791773878e-15i, 
  -5.6385128969624e-05-1.69211755126812e-04i, -5.5860354806709e-162+5.58603548067085e-161i, 
  0.0001631832513714-0.0952324565730093i, 0.69504753678407-1.89164111711036i, 
  1.24241826965328e-04+7.14597582632e-07i, 2.31858732964835e-08+6.18289954572886e-08i, 
  -0.0133426877243506-0.0148087097143221i, 1+0i, 9.9999998533107e-01+2.595272e-08i, 
  -1.5e-15+9.07276596841274e-02i, 0+0.0907276596841274i, 0.0884658993528522+0i, 
  0.034548084541919+0i, 0+0.0630820900592583i, 0.001794352332087+0.0951983814805271i, 
  9.098e-13+7.09979210725139e-02i, 0+0.0201552956479527i, 0+0.0566481651760675i, 
  0+0.0566481651760675i, 5.66e-13+5.66481651760675e-02i, -5.66e-13+5.66481651760675e-02i, 
  0.0796884251721652+0i, 0.0781719582124736-0.0109391367010358i, 
  0.0467003298099045+0.0394403896193353i, 0.367879441171442+0.607157705841394i, 
  0+0.0102596888055368i, 0.990049833749168-0.112088664364495i, 
  1+0i, 0.999999999999944+0i, 0.0110604154853277+0i, 0+0i, 0+0i, 
  0+0i, Inf+0i, 0+0i, complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
  complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
  complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0))),

ref_erf = list(
z  =  c(1e+160-1e+159i, 5.1e-03+1e-08i, -0.0049+0.00495i, 0.0049+0.5i, 
0.00049-5i, -0.000049-50i, 0.0051+0.5i, 0.00051-5i, -0.000051-50i, 
1e-06+2e-06i, 0+2e-06i, 0+2i, 0+20i, 0+200i, Inf+0i, -Inf+0i, 
Inf+0i, -Inf+0i, Inf+0i, complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
complex(real=NaN, imaginary=0), 0.07+0.07i, 7e-02-7e-04i, -9e-02+7e-04i, 
-0.09+0.09i, -7e-04+9e-02i, 0.07+0.009i, 0.081+0i),
value =  c(-0.53664356577856-5.04914370344703i, 0.53664356577856-5.04914370344703i, 
         -0.53664356577856+5.04914370344703i, 0.53664356577856+5.04914370344703i, 
         3.35947367383058e+303-1.99989613967988e+303i, 3.584459971463e+277+3.81895488525718e+279i, 
         0.999602042265715+0.000028010441169i, -1+0i, 1+0i, 5.7546838590348e-03+1.128349818e-08i, 
         -0.00552914914234182+0.00558538838786471i, 0.007099365669981+0.614934701285421i, 
         39811763.38702-8298176341.66525i, -Inf+0i, 0.007389128308257+0.614933252460166i, 
         41436719.23268-8298168216.81831i, -Inf+0i, 1.12837916709965e-06+2.25675833419178e-06i, 
         0+2.25675833419403e-06i, 0+18.5648024145756i, 0+1.47479753962879e+172i, 
         Inf+0i, 1+0i, -1+0i, Inf+0i, -Inf+0i, complex(real=NaN, imaginary=0), 
         complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
         complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
         complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
         complex(real=NaN, imaginary=0), 0.0792438040461578+0.0787277621804668i, 
         0.0788577582851228-0.0007860046704118i, -0.10128064327472+0.000783493474702i, 
         -0.10209984187981+0.101003077889231i, -0.000796289176315+0.101828938593628i, 
         0.0788640866647048+0.0101060428878087i, 0.0912192242781352+0i)
),

ref_erfi = list(
z  =  1.234+0.5678i,
value  =  1.081032284405373149432716643834106923212+1.926775520840916645838949402886591180834i),

ref_erfcx = list(
z  =  1.234+0.5678i,
value  =  0.3382187479799972294747793561190487832579-0.1116077470811648467464927471872945833154i),


ref_erfc = list(
z  =  c(1+2i, -1+2i, 1-2i, -1-2i, 9-28i, 21-33i, 1000+1000i, -3001-1000i, 
        1e+160-1e+159i, 5.1e-03+1e-08i, 0+2e-06i, 0+2i, 0+20i, 0+200i, 
        2e-06+0i, 2+0i, 20+0i, 200+0i, Inf+0i, -Inf+0i, Inf+0i, -Inf+0i, 
        Inf+0i, complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
        complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
        complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
        88+0i),
value  =  c(1.53664356577857+5.04914370344703i, 0.46335643422143+5.04914370344703i, 
            1.53664356577857-5.04914370344703i, 0.46335643422143-5.04914370344703i, 
            -3.35947367383058e+303+1.99989613967988e+303i, -3.584459971463e+277-3.81895488525718e+279i, 
            3.97957734285136e-04-2.8010441169082e-05i, 2+0i, 0+0i, 9.94245316140965e-01-1.1283498e-08i, 
            1-0.00000225675833i, 1-18.5648024145756i, 0e+00-1.47479753962879e+172i, 
            -Inf+0i, 0.999997743241666+0i, 0.00467773498104727+0i, 5.3958656116079e-176+0i, 
            0+0i, 0+0i, 2+0i, -Inf+0i, Inf+0i, complex(real=NaN, imaginary=0), 
            complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
            complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
            complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
            0+0i)
),

ref_Dawson = list(
z  =  c(2+1i, -2+1i, 2-1i, -2-1i, -28+9i, 33-21i, 1000+1000i, -1000-3001i, 
       1e-08+5.1e-03i, 0.00495-0.0049i, 0.0051+0.0051i, 0.5+0.0049i, 
       -5+0.00049i, -50-0.000049i, 5e+02+4.9e-06i, 0.5+0.0051i, -5+0.00051i, 
       -50-0.000051i, 1e-06+2e-06i, 2e-06+0i, 2+0i, 20+0i, 200+0i, 0+0.0049i, 
       0-0.0051i, 0+2e-06i, 0-2i, 0+20i, 0-200i, Inf+0i, -Inf+0i, Inf+0i, 
       -Inf+0i, Inf+0i, complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
       complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
       complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
       39+0.000064i, 41+0.0000609i, 49000000+0i, 51000000+0i, 1e+09+0e+00i, 
       1e+11+0e+00i, 1e+13+0e+00i, 1e+300+0i),
value  =  c(0.163539409434536-0.153124575537123i, -0.163539409434536-0.153124575537123i, 
0.163539409434536+0.153124575537123i, -0.163539409434536+0.153124575537123i, 
-0.016190822566816-0.0052102242033591i, 0.010783770809781+0.0068668887834338i, 
-0.580861681919674+0.668859390550556i, complex(real=NaN, imaginary=0), 
1.000052021e-08+5.10008843492007e-03i, 0.00495015683758159-0.00489983830515523i, 
0.00510017686431968+0.00509982312831979i, 0.424453484087183+0.002820278933187i, 
-0.102134073327105-0.000010456964561i, -1.00020012011921e-02+9.8058859e-09i, 
1.000002000012e-03-9.80006e-12i, 0.424454908562851+0.002935393851312i, 
-0.102134073235712-0.00001088377943i, -1.00020012011913e-02+1.02061261e-08i, 
1.00000000000733e-06+2.00000000000133e-06i, 1.99999999999467e-06+0i, 
0.301340388923792+0i, 0.0250313679264037+0i, 0.00250003125117195+0i, 
0+0.00490007843341994i, 0-0.00510008843492007i, 0+2.00000000000533e-06i, 
0-48.1600121142912i, 0+4.62740702950444e+173i, -Inf+0i, 0+0i, 
0+0i, Inf+0i, -Inf+0i, complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
complex(real=NaN, imaginary=0), complex(real=NaN, imaginary=0), 
complex(real=NaN, imaginary=0), 1.28247314848943e-02-2.10595728e-08i, 
1.21987525342363e-02-1.81304056e-08i, 1.02040816326531e-08-0e+00i, 
9.80392156862745e-09-0e+00i, 5e-10-0e+00i, 5e-12-0e+00i, 5e-14-0e+00i, 
5e-301+0i)
)
)

  for(test in c("Faddeeva_w", "erfc", "erfcx", "erfi", "Dawson")){
    z <- tests[[paste0("ref_", test)]][["z"]]
    ref <- tests[[paste0("ref_", test)]][["value"]]
    calc <- do.call(test, list(z))
    rel.error <- pmax(abs(Re(ref - calc)) / Mod(ref + calc), 
                      abs(Im(ref - calc)) / Mod(ref + calc))
    
    test_that(sprintf("Standard tests for %s are accurate to 1e-13", test), {
    expect_equal(rel.error, 0*rel.error, tolerance = 1e-13)
    })
  }