GSL Replacement

GPL is bad and we're going to replace GSL.

Stub library located at: /usr/pubsw/packages/gsl/stub/include/

See also NumericalRecipesReplacement

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

• 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
• 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
• 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

Possible Replacements

• 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
• VXL/VNL for SV Decomposition
• 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

1. 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,
2. 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).
3. 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.
4. The GSL code is basically of two types: the matrix/vector stuff, 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.
5. Now the tough part: identify suitable equivalents to the GSL statistics stuff that we use (gsl_cdf_ etc.), and replace. Run units tests (the test code take from the orignal GSL code), run system tests (test_recon-all script), and commit to CVS.
6. Run recon-all on the buckner40, and manually check results.

Other Useful Links

• http://www.gnu.org/software/gsl/manual/html_node/Function-Index.html