Rexpokit

Update, June 2017: fixed many new CRAN check issues

(link to this section)

Updates to CRAN's code checks, and updates to operating systems and their FORTRAN-compilation defaults, created new errors/warnings/notes in R CMD check, and installation problems (apparently) on some Windows systems. It looks to me like the main issue was that one of the Windows compilers no longer allows a "-W" FORTRAN compiler flag, thus F77 FORTRAN code is interpreted as F90 code. The solution was to update the F77 source code to F90.

Doing this update was quite a process, as few people in my circles know FORTRAN, let alone the details of standards changes between FORTRAN versions. I am posting my notes/discoveries, below, for google-ability.

Version 0.26 has been submitted to CRAN. Until it gets on CRAN, you can install it from the the GitHub archive via:

library(devtools)
install_github("rexpokit", username="nmatzke")

Update, June 2014: faster rexpokit

(link to this section)

Drew Schmidt, UTK Math Dept., has been working on improvements to rexpokit. These make rexpokit about 40% faster. This will be version 0.25. Until it gets on CRAN, you can install it via:

library(devtools)
install_github("rexpokit", username="wrathematics")

Introduction

(link to this section)

"rexpokit" — an R package by Nicholas J. Matzke (coauthor is Roger B. Sidje, who wrote EXPOKIT; Drew Schmidt has been added as a coauthor as well) which bundles some of the routines from the FORTRAN EXPOKIT library for rapid exponentiation of small dense matrices, and large sparse matrices.

This package is required for use of BioGeoBEARS, along with cladoRcpp.

This package wraps some of the matrix exponentiation utilities from EXPOKIT (http://www.maths.uq.edu.au/expokit/), a FORTRAN library that is widely recommended for matrix exponentiation (Sidje RB, 1998. "Expokit: A Software Package for Computing Matrix Exponentials." ACM Trans. Math. Softw. 24(1): 130-156). EXPOKIT includes functions for exponentiating both small, dense matrices, and large, sparse matrices (in sparse matrices, most of the cells have value 0). Rapid matrix exponentiation is useful in phylogenetics when we have a large number of states (as we do when we are inferring the history of transitions between the possible geographic ranges of a species), but is probably useful in other ways as well.

This package is a part of the work described in:

  • Matzke, Nicholas J. (2014). "Model Selection in Historical Biogeography Reveals that Founder-event Speciation is a Crucial Process in Island Clades." Systematic Biology, 63(6), 951-970. http://dx.doi.org/10.1093/sysbio/syu056
  • Matzke, Nicholas J. (2013). Probabilistic historical biogeography: new models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. Frontiers of Biogeography, 5(4), 242-248. http://escholarship.org/uc/item/44j7n141
  • Matzke N (2012). "Founder-event speciation in BioGeoBEARS package dramatically improves likelihoods and alters parameter inference in Dispersal-Extinction-Cladogenesis (DEC) analyses." Frontiers of Biogeography, 4(suppl. 1), pp. 210. ISSN 1948-6596, Poster abstract published in the Conference Program and Abstracts of the International Biogeography Society 6th Biannual Meeting, Miami, Florida. Poster Session P10: Historical and Paleo-Biogeography. Poster 129B. January 11, 2013. http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster.

Download

(link to this section)

The package has been accepted by CRAN. It is online here: http://cran.r-project.org/web/packages/rexpokit/index.html. Downloads of the package and PDF/example files are also linked to this page. See "Files", below.

License

(link to this section)

The license for rexpokit and all updates / additions / example scripts is the same as for the package on CRAN:

License: GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]

The GPL-3 license can be found at: http://cran.r-project.org/web/licenses/GPL-3

FORTRAN update notes: August 2017

(link to this section)

My notes on the various changes I had to make to the FORTRAN code to get it to pass R CMD check are below. They are rough, but mostly are useful for google-ability in case other R-package maintainers have similar issues.

Change requested by Kurt: Registering the C++ and FORTRAN calls:
https://stat.ethz.ch/pipermail/r-devel/2017-February/073755.html

library(tools)
package_native_routine_registration_skeleton(dir="/GitHub/rexpokit")

#######################################################
# To install by hand, from a source directory:
#######################################################
remove.packages("rexpokit")
install.packages(pkgs="/GitHub/rexpokit", repos=NULL, type="source")
library(rexpokit)
maxent(3.5, 1:6)
remove.packages("rexpokit")

rm rexpokit/src/*.o; rm rexpokit/src/*.so; rm rexpokit/src/lapack/*.o
R CMD check rexpokit

Windows:
https://win-builder.r-project.org/

cd /GitHub
rm rexpokit/src/*.o; rm rexpokit/src/*.so; rm rexpokit/src/lapack/*.o
R CMD check --as-cran rexpokit

rm rexpokit/src/*.o; rm rexpokit/src/*.so; rm rexpokit/src/lapack/*.o
R CMD build rexpokit

#######################################################
# Notes:
#
# I did this install on a Mac OS X 10.11.  I had to:
# 
# 1. Make sure Xcode's xtools were installed (xcode-select, IIRC)
#    xcode-select --install
# 
# 2. Make sure I had an up-to-date install of gcc and gfortran:
#    https://solarianprogrammer.com/2017/05/21/compiling-gcc-macos/
# 
#######################################################

void wrapdgpadm_(int *ideg, int *m, double *t, double *H, int *ldh,
    double *wsp, int *lwsp, int *ipiv, int *iexph, int *ns, int *iflag);

wrapdgpadm(ideg,m,t,H,ldh,wsp,lwsp,ipiv,iexph,
     .                        ns,iflag )
      implicit none

      integer,intent(inout) :: ideg,m,ldh,lwsp,iexph,ns,iflag,ipiv(m)
      double precision,intent(inout) :: t,H(ldh,m),wsp(lwsp)
      integer i,j

void itscale5_(double *SXT, int *ngroups, int *ntraits, double *const, 
  double *prior, double *prob, double *entropy, int *niter, double *tol, double *denom)

subroutine itscale5(SXT,ngroups,ntraits,const,
     & prior,prob,entropy,niter,tol,denom) 

          double precision SXT(ngroups,ntraits),const(ntraits)
      double precision prob(ngroups),prob2(ngroups),prior(ngroups)
      double precision gamma1(ntraits),total,test1,tol
      double precision Csums(ntraits),denom(ntraits),unstand(ngroups)
      double precision entropy

HINTS:

Modernizing Old Fortran
http://fortranwiki.org/fortran/show/Modernizing+Old+Fortran 

COMPLEX(KIND=8) (DOUBLE COMPLEX) Representation
http://kiwi.atmos.colostate.edu/rr/old/tidbits/intel/macintel/doc_files/source/extfile/bldaps_for/pg11fltc.htm
COMPLEX(8) (same as COMPLEX(KIND=8) and COMPLEX*16) data is 16 contiguous bytes containing a pair of REAL(8) values stored in IEEE T_floating format. 

BIG HELP FOR "Computed GO TO" STATEMENTS:
##############################################################      
On-Line Fortran F77 - F90 Converter
https://www.fortran.uk/plusfortonline.php

This free service, offered by Polyhedron Solutions, allows you to submit your Fortran 77 source code (up to 500 lines), and to view the result of processing it using SPAG - part of Polyhedron's plusFORT toolset. A very small subset of the customization options available in the full plusFORT package is available using the drop-down boxes. 

Instructions
Cut and paste your code into the F77 code window, or, if you don't have something ready, click "Load Our Example Code". Then, choose your conversion options and click SUBMIT. The converted code will appear in the F90 code window and a log of the conversion process can be viewed by clicking the 'show output log' button.

For more information about plusFORT, please visit our website, www.fortran.uk or contact us on sales@fortran.uk.

F77 code (must be old style fixed source form)
##############################################################

CHANGES, random notes:

       o Numerous changes were made to FORTRAN source to solve CRAN check
         warnings or installation issues on e.g. Windows machines
         (https://win-builder.r-project.org). Typically these were obsolescence
      issues with legacy F77 code. Typically the solution was to convert to
      F90 standard. Summary below.
       o Fix to a CRAN-identified NOTE: "Packages in Depends field not imported from:
        'Rcpp'  'SparseM' 'methods'."  Solution: removed methods and SparseM from 
        Depends, moved Rcpp to Imports in NAMESPACE file
    o Commented out -W flag in Makevars, in order to remove "Non-portable 
      flags in variable 'PKG_FFLAGS'" warning
    o Summary of modernization edits to FORTRAN:
    o complex*16 -> complex(kind=8)
    o CHARACTER*6 -> CHARACTER(LEN=6)
    o double complex -> complex(kind=8)
    o absx = dabs(dreal(zx(i))) -> absx = DABS(REALPART(Zx(i)))
    o absx = dabs(dimag(zx(i))) -> absx = DABS(IMAGPART(Zx(i)))
    o Removed all "statement functions", replaced with basic math in-line
    o Variable declarations for dreal, dimag, dzum, dzumi, dzumr now unnecessary and
      commented out
    o Replaced "Computed GO TO" statements loops re-coded using On-Line 
      Fortran F77 - F90 Converter: https://www.fortran.uk/plusfortonline.php
    o Problems like "Index '2' of dimension 1 of array 'dx' above upper bound of 1"
      solved with changing dx(1) -> dx(n), dy(1) -> dy(n)
    o cabs1 "statement function" replaced with several in-line calculations
    o All example matrix exponentiations checked after FORTRAN edits
    o Generated "init.c" using package_native_routine_registration_skeleton() in 
      order to register native routines. 
      (Following instructions at: 
      https://stat.ethz.ch/pipermail/r-devel/2017-February/073755.html)
    o init.c content was integrated into pre-existing "wrappers.cpp", adding
      "R_CMethodDef CEntries[]" and "R_init_rexpokit" entries

c     Putting some if statements, so that these are not dummy variables 
      if (SRNAME .EQ. 'DGEMV ') then
         continue
c        print*,"DGEMV problem: no INFO, see XERBLA in blas_mod.f"
c        print*,"Attempting to print INFO below:"
c        print*,INFO
      end if

****************************************************
* COMMON ERRORS THAT HAVE STUPID REASONS
****************************************************
* 1.
* 
* pta = REALPART(a)
* 1
* Error: Non-numeric character in statement label at (1)
*
* THIS MEANS: THE CODE MUST START AT COLUMN 7!!
* 
* 2.
*  if ((dabs(REALPART(a(1,1)))+dabs(IMAGPART(a(1,1)))).eq.0.0d0) &
* Warning: Line truncated at (1) [-Wline-truncation]
* 
* THIS MEANS: Lines have to end at about column 70
* (Because that's how long ticker-tape was in 1965,
*  or something like that)
*
****************************************************

*     NOTE
*     Fixing compiling errors noted by CRAN check in some 
*     Windows machines
*     
*     PROBLEM:
*     my_expokit.f:816:16:
*     complex*16       H(ldh,m), wsp(lwsp)
*     Warning: GNU Extension: Nonstandard type declaration COMPLEX*16 at (1)
*     
*     FIX:
*     complex*16 --> complex(kind=8)
*     
*     
*     
*     PROBLEM:
*     CHARACTER*6        SRNAME
*     Warning: Obsolescent feature: Old-style character length at (1)
*     
*     
*     FIX:
*     CHARACTER*6 --> CHARACTER(LEN=6)
*     
*     
*     
*     PROBLEM:
*     lapack/blas_mod.f:2512:20:
*     double complex zx(1),zy(1),ztemp
*     Warning: GNU Extension: DOUBLE COMPLEX at (1)
*     
*     
*     
*     FIX:
*     double complex --> complex(kind=8)
*     
*     
*     
*     PROBLEM:
*      !absx = dabs(dreal(zx(i)))
*     
*     
*     
*     FIX:
*      absx = DABS(REALPART(Zx(i)))
*     
*     
!     ERROR:
!     lapack/blas_mod.f:1884:26:
!     absx = dabs(dimag(zx(i)))
!     Error: Syntax error in argument list at (1)
!     FIX:
!     NO:   absx = dabs(dimag(zx(i)))
!     NO:   absx = dabs((0.0d0,-1.0d0)*zx(i))
!     YES:  Comment out dimag "statement function",
!           Just use IMAGPART
*             absx = DABS(IMAGPART(Zx(i)))
* 
* 
*     PROBLEM:
*     
*     
*     
*     
*     FIX:
*     
*     
*     
*     
*     PROBLEM:
*     
*     
*     
*     
*     FIX:
*     
*     
*     
*     

*     
*     NOTE -- MODIFIED by Nick Matzke to fix these warnings when
*     compiling with g77:
*     
*     2013-02-15
*     
*     gfortran   -fno-underscoring   -O3  -mtune=core2 -c blas.f -o blas.o
*     blas.f:404.72:
*     
*        10 assign 30 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:409.19:
*     
*        20    go to next,(30, 50, 70, 110)                                   
*                        1
*     Warning: Deleted feature: Assigned GOTO statement at (1)
*     blas.f:411.72:
*     
*           assign 50 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:420.72:
*     
*           assign 70 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:427.72:
*     
*           assign 110 to next                                                
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:1621.72:
*     
*        10 assign 30 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:1628.19:
*     
*              go to next,(30, 50, 70, 90, 110)                               
*                        1
*     Warning: Deleted feature: Assigned GOTO statement at (1)
*     blas.f:1630.72:
*     
*           assign 50 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:1639.72:
*     
*           assign 70 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:1644.72:
*     
*       100 assign 110 to next                                                
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:1671.72:
*     
*        85 assign 90 to next                                                 
*                                                                             1
*     Warning: Deleted feature: ASSIGN statement at (1)
*     blas.f:1689.16:
*     
*           go to next,(  50, 70, 90, 110 )                                   
*                     1
*     Warning: Deleted feature: Assigned GOTO statement at (1)
*     
*     THE FIX IS PROVIDED HERE:
*     http://ubuntuforums.org/showthread.php?t=1578045
*     
*     As you've discovered, assigned goto is deprecated in f90 (and deleted
*     from f95 I believe - not quite sure on that). To fix the code to
*     avoid using deprecated assigned gotos, you can perform the following
*     steps:
*     
*     (a) change each of the ASSIGN N TO NEXT (where N is some number)
*     statement to a simple NEXT = N statement (where N is the same
*     number), and
*     
*     (b) replace each of the GO TO NEXT, (x, y, z, ...) statement with the
*     following computed goto statement
*     
*     GO TO (1,2,3,4,5,6,7,8,9,10,11) NEXT
*     
*     H
*     ===============

Fix for common statement function "cabs1"
(absolute value of a complex number)

c      complex(kind=8) zdum
      double precision cabs1
      double precision pta, ptb

c      double precision dreal,dimag
c      complex(kind=8) zdumr,zdumi
c     dreal(zdumr) = zdumr
c     dimag(zdumi) = (0.0d0,-1.0d0)*zdumi

c     Statement function:
c      cabs1(zdum) = dabs(REALPART(zdum)) + dabs(IMAGPART(zdum))
c     FIX:
c      double precision cabs1
c      double precision pta, ptb
c      pta = REALPART(zdum)
c      ptb = IMAGPART(zdum)
c      ((dabs(pta)+dabs(ptb))

At line 525 of file lapack/blas_mod.f
Fortran runtime error: Index '2' of dimension 1 of array 'dx' above upper bound of 1
** running examples for arch 'x64' ... ERROR
Running examples in 'rexpokit-Ex.R' failed
The error most likely occurred in:

CHANGING dy(1) to dy(n)

c     FIX:
c     double precision dx(1),dy(1),da
      double precision dx(n),dy(n),da

OLD, SINCE R-FORGE SEEMS TO GET STUCK A LOT: R-Forge Project

(link to this section)

The link to the R-Forge Project is here: http://r-forge.r-project.org/projects/rexpokit/.

The latest subversion source code is here: https://r-forge.r-project.org/scm/viewvc.php/?root=rexpokit.

However, R-Forge's R packages page doesn't seem to be updating (it is stuck on an old version). The up-to-date, approximately final version of the package is in "Files", below.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License