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://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 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.
5. 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.
6. 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.
7. 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