This program implements the False Discovery Rate (FDR) algorithm for
thresholding of voxelwise statistics.

Program input consists of a functional dataset containing one (or more)
statistical sub-bricks.  Output consists of a bucket dataset with one
sub-brick for each input sub-brick.  For non-statistical input sub-bricks,
the output is a copy of the input.  However, statistical input sub-bricks
are replaced by their corresponding FDR values, as follows:

For each voxel, the minimum value of q is determined such that
                               E(FDR) <= q
leads to rejection of the null hypothesis in that voxel. Only voxels inside
the user specified mask will be considered.  These q-values are then mapped
to z-scores for compatibility with the AFNI statistical threshold display:

               stat ==> p-value ==> FDR q-value ==> FDR z-score

The reason for the final conversion from q to z is so that larger values
are more 'significant', which is how the usual thresholding procedure
in the AFNI GUI works.

    -input fname       fname = filename of input 3d functional dataset
    -input1D dname     dname = .1D file containing column of p-values

    -mask_file mname   Use mask values from file mname.
     *OR*              Note: If file mname contains more than 1 sub-brick,
    -mask mname        the mask sub-brick must be specified!
                       Default: No mask
                     ** Generally speaking, you really should use a mask
                        to avoid counting non-brain voxels.  However, with
                        the changes described below, the program will
                        automatically ignore voxels where the statistics
                        are set to 0, so if the program that created the
                        dataset used a mask, then you don't need one here.

    -mask_thr m        Only voxels whose corresponding mask value is
                       greater than or equal to m in absolute value will
                       be considered.  Default: m=1

                       Constant c(N) depends on assumption about p-values:
    -cind              c(N) = 1   p-values are independent across N voxels
    -cdep              c(N) = sum(1/i), i=1,...,N   any joint distribution
                       Default:  c(N) = 1

    -quiet             Flag to suppress screen output

    -list              Write sorted list of voxel q-values to screen

    -prefix pname      Use 'pname' for the output dataset prefix name.
    -output pname


January 2008: Changes to 3dFDR
The default mode of operation of 3dFDR has altered somewhat:

 * Voxel p-values of exactly 1 (e.g., from t=0 or F=0 or correlation=0)
     are ignored by default; in the old mode of operation, they were
     included in the count which goes into the FDR algorithm.  The old
     process tends to increase the q-values and so decrease the z-scores.

 * The array of voxel p-values are now sorted via Quicksort, rather than
     by binning, as in the old mode.  This (by itself) probably has no
     discernible effect on the results, but should be faster.

New Options:
    -old     = Use the old mode of operation (for compatibility/nostalgia)
    -new     = Use the new mode of operation [now the default]
                N.B.: '-list' does not work in the new mode!
    -pmask   = Instruct the program to ignore p=1 voxels
                [the default in the new mode, but not in the old mode]
               N.B.: voxels that were masked in 3dDeconvolve (etc.)
                     will have their statistics set to 0, which means p=1,
                     which means that such voxels are implicitly masked
                     with '-new', and so don't need to be explicitly
                     masked with the '-mask' option.
    -nopmask = Instruct the program to count p=1 voxels
                [the default in the old mode, but NOT in the new mode]
    -force   = Force the conversion of all sub-bricks, even if they
                are not marked as with a statistical code; such
                sub-bricks are treated as though they were p-values.
    -float   = Force the output of z-scores in floating point format.
    -qval    = Force the output of q-values rather than z-scores.
                N.B.: A smaller q-value is more significant!
                [-float is strongly recommended when -qval is used]

* To be clear, you can use '-new -nopmask' to have the new mode of computing
   carried out, but with p=1 voxels included (which should give results
   nearly identical to '-old').

* Or you can use '-old -pmask' to use the old mode of computing but where
   p=1 voxels are not counted (which should give results virtually
   identical to '-new').

* However, the combination of '-new', '-nopmask' and '-mask_file' does not
   work -- if you try it, '-pmask' will be turned back on and a warning
   message printed to aid your path towards elucidation and enlightenment.

Other Notes:
* '3drefit -addFDR' can be used to add FDR curves of z(q) as a function
    of threshold for all statistic sub-bricks in a dataset; in turn, these
    curves let you see the (estimated) q-value as you move the threshold
    slider in AFNI.
   - Since 3drefit doesn't have a '-mask' option, you will have to mask
     statistical sub-bricks yourself via 3dcalc (if desired):
       3dcalc -a stat+orig -b mask+orig -expr 'a*step(b)' -prefix statmm
   - '-addFDR' runs as if '-new -pmask' were given to 3dFDR, so that
     stat values == 0 are ignored in the FDR calculations.
   - most AFNI statistical programs now automatically add FDR curves to
     the output dataset header, so you can see the q-value as you adjust
     the threshold slider.

* q-values are estimates of the False Discovery Rate at a given threshold;
   that is, about 5% of all voxels with q <= 0.05 (z >= 1.96) are
   (presumably) 'false positive' detections, and the other 95% are
   (presumably) 'true positives'.  Of course, there is no way to tell
   which above-threshold voxels are 'true' detections and which are 'false'.

* Note the use of the words 'estimate' and 'about' in the above statement!
   In particular, the accuracy of the q-value calculation depends on the
   assumption that the p-values calculated from the input statistics are
   correctly distributed (e.g., that the DOF parameters are correct).

* The z-score is the conversion of the q-value to a double-sided tail
   probability of the unit Gaussian N(0,1) distribution; that is, z(q)
   is the value such that if x is a N(0,1) random variable, then
   Prob[|x|>z] = q: for example, z(0.05) = 1.95996.
  The reason for using z-scores here is simply that their range is
   highly compressed relative to the range of q-values
   (e.g., z(1e-9) = 6.10941), so z-scores are easily stored as shorts,
   whereas q-values are much better stored as floats.

* Changes above by RWCox -- 18 Jan 2008 == Cary Grant's Birthday!

26 Mar 2009 -- Yet Another Change [RWCox]
* FDR calculations in AFNI now 'adjust' the q-values downwards by
   estimating the number of true negatives [m0 in the statistics
   literature], and then reporting
     q_new = q_old * m0 / m, where m = number of voxels being tested.
   If you do NOT want this adjustment, then set environment variable
   AFNI_DONT_ADJUST_FDR to YES.  You can do this on the 3dFDR command
   line with the option '-DAFNI_DONT_ADJUST_FDR=YES'

For Further Reading and Amusement
* cf. http://en.wikipedia.org/wiki/False_discovery_rate [Easy overview of FDR]
* cf. http://dx.doi.org/10.1093/bioinformatics/bti448   [False Negative Rate]
* cf. http://dx.doi.org/10.1093/biomet/93.3.491         [m0 adjustment idea]
* cf. C implementation in mri_fdrize.c                  [trust in the Source]
* cf. https://afni.nimh.nih.gov/pub/dist/doc/misc/FDR/FDR_Jan2008.pdf

++ Compile date = Aug 10 2020 {AFNI_20.2.11:linux_ubuntu_16_64}