Index TableOfContents

Summary

[http://wiki.na-mic.org/Wiki/images/a/ae/NA-MIC-2005-10-30-Licencing.ppt GPL is bad] and we're going to replace GSL. See also NumericalRecipesReplacement

Stub library

Stub library located at: /usr/pubsw/packages/gsl/stub/include/ This is stripped-down GSL code containing just the header defs/structs and emptied routines of the stuff from GSL we use in Freesurfer. You can change your configure line from --with-gsl-dir=/usr/pubsw/packages/gsl/current to --with-gsl-dir=/usr/pubsw/packages/gsl/stub and freesurfer will compile, but when a binary which uses GSL is executed, the emptied routine (which is not quite empty) prints a nasty message, indicating such. This was used to identify the minimum gsl set, and is the precursor to a wrapper .c file.

Freesurfer files using GSL

randomfields.h
mri_cht2p.c
optseq2.c
fmriutils.c
fsglm.c
matrix.c -- lots of matrix functionality that we can replace with vxl
randomfields.c
volcluster.c

GSL function stack trace

gsl_matrix
- matrix.c::MatrixToGSL
  - matrix.c::MatrixSVDPseudoInverse
    - gcamorph.c::gcamComputeOptimalTargetLinearTransform
      - gcamorph.c::GCAMreinitWithLTA
        mri_ca_register.c::main
- matrix.c::MatrixFromGSL
  - matrix.c::MatrixSVDPseudoInverse
    - ...
gsl_matrix_alloc
- matrix.c::MatrixToGSL
  - ...
- matrix.c::MatrixSVDPseudoInverse
  - ...
gsl_matrix_calloc
- optseq2.c::AR1WhitenMatrix
  - optseq2.c::main
gsl_matrix_free
- optseq2.c::AR1WhitenMatrix
  - ...
- matrix.c::MatrixSVDPseudoInverse
  - ...
gsl_matrix_get
- optseq2.c::AR1WhitenMatrix
  - ...
- matrix.c::MatrixFromGSL
  - ...
gsl_matrix_set
- optseq2.c::AR1WhitenMatrix
  - ...
- matrix.c::MatrixFromGSL
  - ...
gsl_vector
- matrix.c::VectorFromGSL
  - matrix.c::MatrixSVDPseudoInverse
    - ...
- matrix.c::MatrixSVDPseudoInverse
  - ...
gsl_vector_alloc
- matrix.c::MatrixSVDPseudoInverse
  - ...
gsl_vector_free
- matrix.c::MatrixSVDPseudoInverse
  - ...
gsl_vector_get
- matrix.c::VectorFromGSL
  - matrix.c::MatrixSVDPseudoInverse
    - ...
gsl_vector_set -- never used
gsl_cdf_fdist_Q
- fmriutils.c::fMRIsigF
  - mris_glm.c::main
- fsglm.c::GLMtest
  - fsglm.c::GLManalyze -- never used
  - fsglm.c::GLMprofile
    - mri_glmfit::parse_commandline
      - mri_glmfit::main
  - fsglm.c::GLMsynth -- never used
- randomfields.c::RFstat2PVal
  - randomfields.c::RFstat2P
    - randomfields.c::RFstat2Stat -- never used
gsl_cdf_flat_Q
- randomfields.c::RFstat2PVal
  - ... -- never used
gsl_cdf_gaussian_Q
- randomfields.c::RFstat2PVal
  - ... -- never used
gsl_cdf_tdist_Q
- randomfields.c::RFstat2PVal
  - ... -- never used
gsl_cdf_chisq_Q
- randomfields.c::RFstat2PVal
  - ... -- never used
gsl_cdf_flat_Qinv
- randomfields.c::RFp2StatVal
  - randomfields.c::RFp2Stat
    - randomfields.c::RFstat2Stat -- never used
gsl_cdf_gaussian_Qinv
- randomfields.c::RFp2StatVal
  - ...
gsl_cdf_tdist_Qinv
- randomfields.c::RFp2StatVal
  - ... -- never used
gsl_cdf_chisq_Qinv
- randomfields.c::RFp2StatVal
  - ... -- never used
gsl_linalg_SV_decomp
- matrix.c::MatrixSVDPseudoInverse
  - ...
gsl_linalg_cholesky_decomp
- optseq2.c::AR1WhitenMatrix
  - ...
gsl_rng_type
- RANDOM_FIELD_SPEC, RFS
  - used everywhere in randomfied.c
- randomfields.c::RFspecInit
  - mri_glmfit.c::main
  - mri_mcsim.c::main
  - mri_volsynth.c::main
gsl_rng
- RANDOM_FIELD_SPEC, RFS
gsl_rng_alloc
- randomfields.c::RFspecInit
  - ...
gsl_rng_free
- randomfields.c::RFspecFree
  - RFspecFree -- never used
gsl_rng_set
- randomfields.c::RFspecSetSeed
  - randomfields.c::RFspecInit
    - ...
gsl_ran_flat
- randomfields.c::RFdrawVal
  - randomfields.c::RFsynth
    - mri_glmfit.c::main
    - mri_mcsim.c::main
    - mri_volsynth.c::main
gsl_ran_gaussian
- randomfields.c::RFdrawVal
  - ...
gsl_ran_fdist
- randomfields.c::RFdrawVal
  - ...
gsl_ran_tdist
- randomfields.c::RFdrawVal
  - ...
gsl_ran_chisq
- randomfields.c::RFdrawVal
  - ...
gsl_ran_binomial_pdf
- mri_cht2p.c::main
gsl_ran_binomial_pdf
- volcluster.c::CSDpvalClustSize
  - mri_surfcluster.c::main

GSL types and functions used

matrix/vector

gsl_matrix
- optseq2.c::AR1WhitenMatrix
- matrix.c::MatrixToGSL
- matrix.c::MatrixFromGSL
- matrix.c::MatrixSVDPseudoInverse
gsl_matrix_alloc
- matrix.c::MatrixToGSL
- matrix.c::MatrixSVDPseudoInverse
gsl_matrix_calloc
- optseq2.c::AR1WhitenMatrix
gsl_matrix_free
- optseq2.c::AR1WhitenMatrix
- matrix.c::MatrixSVDPseudoInverse
gsl_matrix_get
- optseq2.c::AR1WhitenMatrix
- matrix.c::MatrixFromGSL
gsl_matrix_set
- optseq2.c::AR1WhitenMatrix
- matrix.c::MatrixFromGSL
gsl_vector
- matrix.c::MatrixSVDPseudoInverse
- matrix.c::VectorFromGSL
gsl_vector_alloc
- matrix.c::MatrixSVDPseudoInverse
gsl_vector_free
- matrix.c::MatrixSVDPseudoInverse
gsl_vector_get
- matrix.c::VectorFromGSL
gsl_vector_set -- not used

random number generator

gsl_rng_type
- randomfields.h::RANDOM_FIELD_SPEC, RFS
- randomfields.h::RFspecInit
- randomfields.c::RFspecInit
gsl_rng_ranlux389
- randomfields.c::RFspecInit
gsl_rng
- randomfields.h::RANDOM_FIELD_SPEC, RFS
gsl_rng_alloc
- randomfields.c::RFspecInit
gsl_rng_free
- randomfields.c::RFspecFree
gsl_rng_set
- randomfields.c::RFspecSetSeed
gsl_ran_flat
- randomfields.c::RFdrawVal
gsl_ran_gaussian
- randomfields.c::RFdrawVal
gsl_ran_fdist
- randomfields.c::RFdrawVal
gsl_ran_tdist
- randomfields.c::RFdrawVal
gsl_ran_chisq
- randomfields.c::RFdrawVal
gsl_ran_binomial_pdf
- mri_cht2p.c::main
- volcluster.c::CSDpvalClustSize

statistics functions

gsl_cdf_fdist_Q
- fmriutils.c::fMRIsigF
- fsglm.c::GLMtest
- randomfields.c::RFstat2PVal
gsl_cdf_flat_Q
- randomfields.c::RFstat2PVal
gsl_cdf_gaussian_Q
- randomfields.c::RFstat2PVal
gsl_cdf_tdist_Q
- randomfields.c::RFstat2PVal
- sig.c::sigt -- only used if MGH_GSL is defined and it never is
gsl_cdf_chisq_Q
- randomfields.c::RFstat2PVal
gsl_cdf_flat_Qinv
- randomfields.c::RFp2StatVal
gsl_cdf_gaussian_Qinv
- randomfields.c::RFp2StatVal
gsl_cdf_tdist_Qinv
- randomfields.c::RFp2StatVal
gsl_cdf_chisq_Qinv
- randomfields.c::RFp2StatVal
gsl_linalg_SV_decomp
- matrix.c::MatrixSVDPseudoInverse
gsl_linalg_cholesky_decomp
- optseq2.c::AR1WhitenMatrix

Possible replacements

http://www.netlib.org/cephes/ and http://www.moshier.net/#Cephes -- The probability functions used by the FSL people come from the Cephes Mathematical Library
VXL/VNL for SV Decomposition
http://www.thefreecountry.com/sourcecode/mathematics.shtml -- Free Numerical Libraries, Mathematical Libraries, Statistics Libraries
http://netlib.org/ -- BLAS, CBLAS, CLAPACK, ATLAS. maybe?
- Advantages: gold standard for matrix software (used by matlab), real and complex, BSD, possible to optimize for individual platforms, parallel version available.
- Disads: poorly documented (and I mean *poorly*), routines are low-level and would require work up-front to get going, no statistics.
http://math.nist.gov/MatrixMarket/ -- no
- doesn't have the stat functions we need
http://www.oonumerics.org/oon/ -- no
- a long list of stuff we can't use
http://rac.uits.iu.edu/hpc/numerics.shtml -- no
- long list of libraries
http://osl.iu.edu/research/mtl/ -- no
- doesn't have the stat functions we need
http://math.nist.gov/tnt/ -- no
- doesn't have the stat functions we need
http://www.gnu.org/software/goose/goose.html
- LGPL
- has CDF
ITK for some of the distribution functions?
http://cran.us.r-project.org R statistics package. Probably has all the stats we need. Don't know about the matrix stuff. Oops, I think it's GPL though. Yes, it is under GPL.
Probability Density Functions: http://www.scs.fsu.edu/~burkardt/cpp_src/prob/prob.html
- has CDF functions, but not for evaluating Q alone, which is what we want...
- not much in terms of random number generators either
Cumulative Density Functions: http://www.scs.fsu.edu/~burkardt/cpp_src/dcdflib/dcdflib.html -- no

Steps for replacement

Create wrapper file (gsl_wrapper.c in /utils, and gsl_wrapper.h in /include) to wrap existing gsl functions with replacement routines (example: create new routine called 'cdf_flat_Q' which wraps 'gsl_cdf_flat_Q'). The basis for this wrapper file can be taken from the /usr/pubsw/packages/gsl/stub directory,
Create test file (utils/gsl_wrapper_test.c) containing tests of gsl function, as taken from the GSL source code (found at /usr/pubsw/packages/gsl/current/src/gsl-1.6.tar.gz, where each sub-package has its own test.c file).
Change all freesurfer code to call the wrapper functions instead of the gsl functions directly. Commit the gsl wrappers, gsl test files, and code changes to CVS dev. Essentially, no code has changed, so freesurfer should behave exactly the same. This is the baseline. Tag it as 'pre_gsl_replacement', so that its possible to revert to that code line.
The GSL code we use is basically of three types: the matrix/vector stuff, random number generation, and the statistic stuff. First, replace the matrix/vector stuff with the VNL equivalent that is already being used to replace the Numerical Recipes (NRC) stuff. Once this change-over is made, run test_recon-all to compare recon-all results before and after the change. Commit changes to CVS dev.
Now replace the rng stuff with the random number generator found in utils/utils.c. It is important that all random number usage in freesurfer goes through the rng routines found in utils/utils.c. That way, we can control the source of this variability in one place. After replacement, run system tests (test_recon-all), and commit to CVS.
Now the tough part: identify suitable equivalents to the GSL statistics stuff that we use (gsl_cdf_ etc.), and replace. Make sure the replacement code is not somehow a GPL derivative! Run units tests (the test code taken from the orignal GSL code), run system tests (test_recon-all script), and commit to CVS.
Run recon-all on the buckner40, and manually check results.