diff --git a/Makefile b/Makefile
index 940ec2d3af7611a8a6685e100c8914e1c33ca7e4..744fb651452b52d246487b1c5311cb894aaac5d9 100644
--- a/Makefile
+++ b/Makefile
@@ -16,8 +16,6 @@ TEST_OBJS = test.o
GGMIX_OBJS = ggmix.o
-FSL_DR_OBJS = melhlprfns.o fsl_dualreg.o
-
FSL_GLM_OBJS = melhlprfns.o fsl_glm.o
FSL_SBCA_OBJS = melhlprfns.o fsl_sbca.o
@@ -34,7 +32,7 @@ MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o
TESTXFILES = test
-XFILES = fsl_dualreg fsl_glm fsl_sbca fsl_mvlm fsl_regfilt fsl_schurprod fsl_bwtf melodic
+XFILES = fsl_glm fsl_sbca fsl_mvlm fsl_regfilt fsl_schurprod fsl_bwtf melodic
OTHERFILES =
@@ -65,7 +63,7 @@ RUNTCLS = Melodic
SCRIPTS = melodicreport dual_regression
-all: ggmix fsl_regfilt fsl_glm fsl_dr melodic ${OTHERFILES}
+all: ggmix fsl_regfilt fsl_glm melodic ${OTHERFILES}
ggmix: ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS}
@@ -76,9 +74,6 @@ melodic: ${MELODIC_OBJS}
test: ${TEST_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${TEST_OBJS} ${LIBS}
-fsl_dr: ${FSL_DR_OBJS}
- $(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_DR_OBJS} ${LIBS}
-
fsl_glm: ${FSL_GLM_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_GLM_OBJS} ${LIBS}
diff --git a/doc/anim_gui.gif b/doc/anim_gui.gif
deleted file mode 100644
index e8be7a74f85b186ee9e0d9d8a4d46bde700c4b17..0000000000000000000000000000000000000000
Binary files a/doc/anim_gui.gif and /dev/null differ
diff --git a/doc/concat_diag.png b/doc/concat_diag.png
deleted file mode 100644
index 252b34034c426ffa797762ae538f090312dfb208..0000000000000000000000000000000000000000
Binary files a/doc/concat_diag.png and /dev/null differ
diff --git a/doc/gui.png b/doc/gui.png
deleted file mode 100644
index dfcab87a5fcd4683f8b5d33ac692e4bb53f39036..0000000000000000000000000000000000000000
Binary files a/doc/gui.png and /dev/null differ
diff --git a/doc/index.html b/doc/index.html
deleted file mode 100644
index 66af2ead58b0a20c57d489e78a2e6d8291784a21..0000000000000000000000000000000000000000
--- a/doc/index.html
+++ /dev/null
@@ -1,260 +0,0 @@
-<HTML><HEAD><link REL="stylesheet" TYPE="text/css" href="../fsl.css"><TITLE>FSL</TITLE></HEAD><BODY><TABLE BORDER=0 WIDTH="100%"><TR><TD ALIGN=CENTER><H1>
-MELODIC v3.0
-</H1>
-Multivariate Exploratory Linear Optimized Decomposition into Independent Components
-<TD ALIGN=RIGHT><a href="../index.html"><IMG BORDER=0 SRC="../images/fsl-logo.jpg"></a></TR></TABLE><HR>
-
-<IMG ALIGN=RIGHT hspace=20 vspace=20 width=90% SRC="tica_diag.png"
-ALT="TICA diagram">
-
-<H2>INTRODUCTION</H2>
-
-<P>MELODIC 3.0 uses Independent Component Analysis to decompose a
-single or multiple 4D data sets into different spatial and temporal
-components. For ICA group analysis, MELODIC uses either Tensorial
-Independent Component Analysis (TICA, where data is decomposed into
-spatial maps, time courses and subject/session modes) or a simpler
-temporal concatenation approach. MELODIC can pick out different
-activation and artefactual components without any explicit time series
-model being specified.
-
-<P>For more detail on MELODIC and an updated journal reference, see
-the <A
-HREF="http://www.fmrib.ox.ac.uk/analysis/research/melodic/">MELODIC
-research web page</A>. If you use MELODIC in your research, please
-quote the journal reference listed there.
-
-<p>The different MELODIC programs are:
-<UL>
-<LI><a href="#MelodicGUI">Melodic</a> - MELODIC GUI<br>
-<LI><a href="#melodic">melodic</a> - command-line MELODIC program<br>
-<LI><a href="#regfilt">fsl_regfilt</a> - command-line tool for removing regressors from data (melodic denoising)<br>
-</UL>
-
-
-<A NAME="MelodicGUI"></A><hr><H2>Melodic GUI</H2>
-
-<IMG ALIGN=RIGHT hspace=20 vspace=20 SRC="gui.png"
-ALT="Example GUI view">
-
-<p>To call the MELODIC GUI, either type <b>Melodic</b> in a terminal
- (type <b>Melodic_gui</b> on Mac), or run <b>fsl</b> and press
- the <b>MELODIC</b> button.
-
-<p>Before calling the GUI, you need to prepare each session's
-data as a 4D NIFTI or Analyze format image; there are utilities in
-fsl/bin called <a href="../avwutils/index.html" target="_top">fslmerge
-and fslsplit</a> to convert between multiple 3D images and a single 4D
-(3D+time) image.
-
-<p>Structural images for use as "highres" images in registration
-should normally be brain-extracted using <a href="../bet2/index.html"
-target="_top">BET</a>.
-
-
-<H3>GUI details:</H3>
-<UL>
-<LI><a href="#misc">Misc</a><br>
-<LI><a href="#data">Data</a><br>
-<LI><a href="#prestats">Pre-Stats</a><br>
-<LI><a href="#reg">Registration</a><br>
-<LI><a href="#stats">Stats</a><br>
-<LI><a href="#poststats">Post-Stats</a><br>
-<LI><a href="#buttons">Bottom Row of Buttons</a><br>
-</UL>
-<UL>
-<LI><a href="#output">MELODIC report output</a><br>
-</UL>
-
-<a name="misc"></a>
-<hr><H3>Misc</H3>
-
-<p><b>Balloon help</b> (the popup help messages in the MELODIC GUI) can
-be turned off once you are familiar with the GUI.
-
-<p>The <b>Progress watcher</b> button allows you to tell Melodic not to
-start a web browser to watch the analysis progress. If you are
-running lots of analyses you probably want to turn this off; you can
-view the same logging information by looking at the report_log.html or log.txt
-files in any MELODIC directories instead.
-
-<a name="data"></a>
-<hr><H3>Data</H3>
-
-<p>First, set the filename of the 4D input image
-(e.g. <b>/users/sibelius/origfunc.nii.gz</b>) by pressing <b>Select 4D
-data</b>. You can select multiple files if you want MELODIC to perform
-a group analysis or if you want to run separate ICAs with the same
-setup. Results for each input file will be saved in separate .ica
-directories, the name of which is based on the input data's filename
-(unless you enter an <b>Output directory</b> name).
-
-<p><b>Delete volumes</b> controls the number of initial FMRI volumes
-to delete before any further processing.
-
-<p><b>TR</b> controls the time (in seconds) between scanning
-successive FMRI volumes. Changes here will not affect the analysis and only change the x-axis units of the final time series plots.
-
-<p>The <b>High pass filter cutoff</b> controls the longest temporal
-period that you will allow.
-
-<a name="prestats"></a>
-<hr><H3>Pre-Stats</H3>
-
-<p>Low-frequency drifts and motion in the data can adversely affect
-the decomposition. In most cases, you would want to motion-correct the
-data, remove these drifts first or perform other types of typical data pre-processing before running the analysis. This can be done from within the
-Melodic GUI <code>Pre-stats</code> section.
-
-<a name="reg"></a>
-<hr><H3>Registration</H3>
-
-<p>Before any multi-session or multi-subject analyses can be carried
-out, the different sessions need to be registered to each other. This
-is made easy within MELODIC which performs registration on input data as part of an analysis using FEAT functionality. Unlike registration step in FEAT this here needs to be performed before the statistical analysis so that the filtered functional data is transformed into the standard space. For information on using multi-stage registration please consult the <a href="../feat5/detail.html#reg" target="_top">FEAT</a> manual.
-
-<p><b>Standard space</b> refers to the standard (reference) image; it
-should be an image already in standard space, ideally with the
-non-brain structures already removed.
-
-<p><b>Resampling resolution (mm)</b> refers to the desired isotropic voxel dimension of the resampled data. In order to save on disk space and on required memory during the analysis it is advisable to resample the filtered data into standard space but keeping the resampled resolution at the FMRI resolution (typically 4mm or 5mm).
-
-<a name="stats"></a>
-<hr><H3>Stats</H3>
-
-<p>The Stats section lets you control some of the options for the decomposition. The default setting will most probably already be set to what you would want most of the time.
-
-<p> By default, MELODIC will variance-normalise timecourses.
-
-<p> By default, Melodic will automatically estimate the number of
- components from the data - you can switch this option off and then can specify the number of components explicitly.<p>
-
-<p> You can now select the type of analysis. MELODIC currently offers three options:
-<p>
-<table border=0>
- <TR><TD width=50%>
- <UL>
- <LI><b>Single-session ICA:</b> This will perform standard 2D ICA on each of the input files. The input data will each be represented as a 2D time x space matrix. MELODIC then de-composes each matrix separately into pairs of time courses and spatial maps. The original data is assumed to be the sum of outer products of time courses and spatial maps. All the different time courses (one per component) will be saved in the <i>mixing matrix</i> <code>melodic_mix</code> and all the spatial maps (one per component) will be saved in the 4D file <code>melodic_IC</code>.
- <p>When using separate analyses, MELODIC will attempt to find components which are relevant and non-Gaussian relative to the residual fixed-effects within session/subject variation. It is recommended to use this option in order to check for session-specific effects (such as MR artefacts). You will need to use this option if you want to perform MELODIC denoising using <a href="#regfilt">fsl_regfilt</a>. When using single-session ICA the component are ordered in order of decreasing amounts of uniquely explained variance.
-
- </UL></TD><TD valign=top>
- <IMG ALIGN=RIGHT hspace=20 vspace=20 width=80% SRC="pica_diag.png" ALT="PICA diag">
- </TD></TR><TR><TD width=50%>
- <UL>
- <LI><b>Multi-session temporal concatenation:</b> This will perform a single 2D ICA run on the concatenated data matrix (obtained by stacking all 2D data matrices of every single data set on top of each other).
-<p>
- It is recommended to use this approach in cases where one is looking for common spatial patterns but can not assume that the associated temporal response is consistent between sessions/subjects. Examples include activation studies where the design was randomised between sessions or the analysis of data acquired without stimulation (<i>resting-state FMRI</i>).
- <p>This approach does not assume that the temporal response pattern is the same across the population, though the final web report will contain the first Eigenvector of all different temporal responses as a summary time course. Access to all time courses is available: the time series plot is linked to a text file (<code>tXX.txt</code>) which contains the first Eigenvector, the best model fit in case a time series design was specified and all different subject/session-specific time courses as columns.
-
-For each component the final mixing matrix <code>melodic_mix</code> contains the temporal response of all different data sets concatenated into a single column vector. The final reported time course will be the best rank-1 approximation to these different responses. <BR>
- </UL>
- </TD> <TD valign=top>
- <IMG ALIGN=RIGHT hspace=20 vspace=20 width =80% SRC="concat_diag.png" ALT="CONCAT diag">
- </TD></TR><TR><TD width=50%>
- <UL>
- <LI><b>Multi-session Tensor-ICA:</b> This will perform a 3D Tensor-ICA decomposition of the data. All individual data sets will be represented as a single time x space x sessions/subjects block of data. Tensor-ICA will decompose this block of data into triplets of time courses, spatial maps and session/subject modes, which - for each component - characterise the signal variation across the temporal, spatial and subject/session domain.
-<p>It is recommended to use this approach for data where the stimulus paradigm is consistent between session/subjects. Tensor-ICA assumes that the temporal response pattern is the same across the population and provides a single decomposition for all original data sets. MELODIC will attempt to find components which are highly non-Gaussian relative to the full mixed-effects variance of the residuals.
-<p>Estimated components typically fall into 2 classes: components which describe effects common to all or most subjects/sessions, and components which describe effects only contained in a small number of subjects/sessions. The former will have a non-zero estimated effect size while the latter will have an effect size around 0 for most subjects/sessions and only few high non-zero values. These different types of components can be identified easily by looking at the boxplots provided. When using Tensor-ICA the components are ordered in order of decreasing amount of median response amplitude. For details on the decomposition see the technical report <a href="http://www.fmrib.ox.ac.uk/analysis/techrep/"> TR04CB1 </a>.
-</UL> </TD><TD valign=top>
- <IMG ALIGN=RIGHT hspace=20 vspace=20 width =80% SRC="tica_diag.png" ALT="TICA diag">
-</TD></TR>
-</table>
-
-
-<a name="poststats"></a>
-<hr><H3>Post-Stats</H3>
-
-<p> Melodic will also by default carry out inference on the estimated maps
- using a mixture model and an alternative hypothesis testing approach. A threshold level of 0.5
- in the case of alternative hypothesis testing means that a
- voxel 'survives' thresholding as soon as the probability
- of being in the 'active' class (as modelled by the Gamma
- densities) exceeds the probability of being in the
- 'background' noise class. This threshold level assumes that
- you are placing an equal loss on false-positives and
- false-negatives. If, however, you consider e.g. false-positives as
- being twice as bad as false-negatives you should change this value
- to 0.66...
-
-<p> You can select the background image used for the generation of the
- spatial map overlay images.
-
-<p> If you select the <b>Output full stats folder</b> option, MELODIC will save thresholded maps and probability maps in a <code>/stats</code> subdirectory within its output folder.
-
- <p>You can specify a temporal design matrix (and in the case
- of a group analysis also, a session/subject design matrix) as well as corresponding contrast matrices. If these matrices are set in the GUI, MELODIC will perform a post-hoc regression analysis on estimated time courses and session/subject modes. This can be a helpful tool in order to identify whether or not a given component is task related. The matrices themselves can be created easily using the <a href=../feat5/programs.html><b>Glm</b></a> GUI.
-
- <a name="buttons"></a>
- <hr><H3>Bottom Row of Buttons</H3>
-
- <p>When you have finished setting up MELODIC, press <b>Go</b> to run the
- analysis. Once MELODIC is running, you can either <b>Exit</b> the GUI, or
- setup further analyses.
-
- <p>The <b>Save</b> and <b>Load</b> buttons enable you to save and load
- the complete MELODIC setup to and from file.
-
-<a name="output"></a>
-<hr><H3>MELODIC report output</H3>
-
-Melodic will then generate the results and
-your terminal window will tell you where to find the web report.
-
-Each IC_XX.html webpage shows one spatial map thresholded and rendered on top of a background image
-followed by the relevant time-course of the ICA decomposition and the power-spectrum of the time-course.
-
-<p> If you click on the thresholded map, you can inspect the raw IC output together with probability maps and the Mixture Model fit.
-
-
-<p>In the case of TICA or simple time series concatenation the time course plotted is the rank-1 approximation to all the different time courses that correspond to the given spatial map within the population.
-
-<p>If a temporal design was specified in the <a href="#poststats" target="_top">Post-Stats</a> section then the time series plot will also contain a plot of the total model fit. In addition, a simple GLM table will describe the fit in detail, providing information of the regression parameter estimates (PEs). Furthermore, MELODIC will perform a simple F-test on the estimated time course and the total model fit. For task-related components the model fit will explain large amounts of the variation contained in the estimated time couse. In addition, if a contrast matrix was specified, the table will also contain Z-statistics and p-values for all the contrasts.
-
-If a group analysis was carried out then the report page will also include information on the distribution of the effect size across the population. A simple plot and a boxplot show the relative effect size across the different sessions/subjects. If a design matrix was specified in the GUI setup then MELODIC will also include a GLM regression fit table.
-
-<A NAME="melodic"></A><HR><H2>melodic command-line program</H2>
-
-<p>Type <b>melodic --help</b> to get usage.
-
-<A name="regfilt"></A><HR><H2>fsl_regfilt command-line program</H2>
-
-<p>Running MELODIC can be a useful tool for gaining insight into unexpected artefacts or activation in your data.
-
-<p>As well as being a good way to find structured noise (or unexpected
- activation) in your data, ICA can also be used to remove chosen
- components (normally obvious scanner-related or physiological
- artefacts) from your data in order, for example, in order to improve
- the FEAT results. In order to do this:
-
- <UL>
-
- <LI> Run MELODIC single-session ICA on a 4D image file
-
- <LI> Open the MELODIC report
- (<code>melodic_output_directory.ica/filtered_func_data.ica/report/00index.html</code>)
- in a web browser and look through the components to identify those
- that you wish to remove; record the list of component numbers to
- remove.
-
- <LI> In a terminal, run the MELODIC denoising, using the
- commands:<pre>
-cd melodic_output_directory.ica
-fsl_regfilt -i filtered_func_data -o denoised_data -d filtered_func_data.ica/melodic_mix -f "2,5,9"</pre>
- where you should replace the comma-separated list of component numbers with the list that you previously recorded when viewing the MELODIC report.<br>
- </UL>
- The output file <code> denoised_data.nii.gz</code> then contains the filtered and denoised data set which can be used e.g. within FEAT. When running FEAT on this data make sure that the analysis is set to <code>Stats + Post-stats </code> as you do not want to run the other filtering steps (smoothing etc.) again on this data.
-
-
-<p><HR><FONT SIZE=1>Copyright © 2001-2007, University of
-Oxford. Written by <A
- HREF="http://www.fmrib.ox.ac.uk/~beckmann/">C.F. Beckmann </A> and <A
-HREF="http://www.fmrib.ox.ac.uk/~steve/index.html">S. Smith</A>.</FONT>
-
-<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
-<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
-<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
-<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
-
-</BODY></HTML>
-
-
diff --git a/doc/mel1.png b/doc/mel1.png
deleted file mode 100644
index 6de5835db8a62e54cdd441c61893dc087cc32b03..0000000000000000000000000000000000000000
Binary files a/doc/mel1.png and /dev/null differ
diff --git a/doc/mel2.png b/doc/mel2.png
deleted file mode 100644
index e61c591c062687ff32674e385b8a035e481b28d5..0000000000000000000000000000000000000000
Binary files a/doc/mel2.png and /dev/null differ
diff --git a/doc/mel3.png b/doc/mel3.png
deleted file mode 100644
index c3850b8f214535bcc7f2980494a48d2ec7965ae0..0000000000000000000000000000000000000000
Binary files a/doc/mel3.png and /dev/null differ
diff --git a/doc/pica_diag.png b/doc/pica_diag.png
deleted file mode 100644
index 120863cc9c8b60c640ad1a7521d954eeace1cc6a..0000000000000000000000000000000000000000
Binary files a/doc/pica_diag.png and /dev/null differ
diff --git a/doc/tica_diag.png b/doc/tica_diag.png
deleted file mode 100644
index 0e9199cc5d594731aa618cdd1735764ff9894b7d..0000000000000000000000000000000000000000
Binary files a/doc/tica_diag.png and /dev/null differ
diff --git a/dual_regression b/dual_regression
deleted file mode 100755
index 4806200f6b165d79075835bb4b16b5e7a56ca5ea..0000000000000000000000000000000000000000
--- a/dual_regression
+++ /dev/null
@@ -1,128 +0,0 @@
-#!/bin/sh
-
-# dual_regression - take group-ICA maps (etc) and get subject-specific versions of them (and associated timecourses)
-#
-# Stephen Smith and Christian Beckmann, FMRIB Image Analysis Group
-#
-# Copyright (C) 2011-2012 University of Oxford
-#
-# SHCOPYRIGHT
-
-
-Usage() {
- cat <<EOF
-
-dual_regression v0.5 (beta)
-
-***NOTE*** ORDER OF COMMAND-LINE ARGUMENTS IS DIFFERENT FROM PREVIOUS VERSION
-
-Usage: dual_regression <group_IC_maps> <des_norm> <design.mat> <design.con> <n_perm> <output_directory> <input1> <input2> <input3> .........
-e.g. dual_regression groupICA.gica/groupmelodic.ica/melodic_IC 1 design.mat design.con 500 grot \`cat groupICA.gica/.filelist\`
-
-<group_IC_maps_4D> 4D image containing spatial IC maps (melodic_IC) from the whole-group ICA analysis
-<des_norm> 0 or 1 (1 is recommended). Whether to variance-normalise the timecourses used as the stage-2 regressors
-<design.mat> Design matrix for final cross-subject modelling with randomise
-<design.con> Design contrasts for final cross-subject modelling with randomise
-<n_perm> Number of permutations for randomise; set to 1 for just raw tstat output, set to 0 to not run randomise at all.
-<output_directory> This directory will be created to hold all output and logfiles
-<input1> <input2> ... List all subjects' preprocessed, standard-space 4D datasets
-
-<design.mat> <design.con> can be replaced with just
--1 for group-mean (one-group t-test) modelling.
-If you need to add other randomise option then just edit the line after "EDIT HERE" below
-
-EOF
- exit 1
-}
-
-############################################################################
-
-[ "$6" = "" ] && Usage
-
-ORIG_COMMAND=$*
-
-ICA_MAPS=`${FSLDIR}/bin/remove_ext $1` ; shift
-
-DES_NORM=--des_norm
-if [ $1 = 0 ] ; then
- DES_NORM=""
-fi ; shift
-
-if [ $1 = "-1" ] ; then
- DESIGN="-1"
- shift
-else
- dm=$1
- dc=$2
- DESIGN="-d $1 -t $2"
- shift 2
-fi
-
-NPERM=$1 ; shift
-
-OUTPUT=`${FSLDIR}/bin/remove_ext $1` ; shift
-
-while [ _$1 != _ ] ; do
- INPUTS="$INPUTS `${FSLDIR}/bin/remove_ext $1`"
- shift
-done
-
-############################################################################
-
-mkdir $OUTPUT
-LOGDIR=${OUTPUT}/scripts+logs
-mkdir $LOGDIR
-echo $ORIG_COMMAND > $LOGDIR/command
-if [ "$DESIGN" != -1 ] ; then
- /bin/cp $dm $OUTPUT/design.mat
- /bin/cp $dc $OUTPUT/design.con
-fi
-
-echo "creating common mask"
-j=0
-for i in $INPUTS ; do
- echo "$FSLDIR/bin/fslmaths $i -Tstd -bin ${OUTPUT}/mask_`${FSLDIR}/bin/zeropad $j 5` -odt char" >> ${LOGDIR}/drA
- j=`echo "$j 1 + p" | dc -`
-done
-ID_drA=`$FSLDIR/bin/fsl_sub -T 10 -N drA -l $LOGDIR -t ${LOGDIR}/drA`
-cat <<EOF > ${LOGDIR}/drB
-#!/bin/sh
-\$FSLDIR/bin/fslmerge -t ${OUTPUT}/maskALL \`\$FSLDIR/bin/imglob ${OUTPUT}/mask_*\`
-\$FSLDIR/bin/fslmaths $OUTPUT/maskALL -Tmin $OUTPUT/mask
-\$FSLDIR/bin/imrm $OUTPUT/mask_*
-EOF
-chmod a+x ${LOGDIR}/drB
-ID_drB=`$FSLDIR/bin/fsl_sub -j $ID_drA -T 5 -N drB -l $LOGDIR ${LOGDIR}/drB`
-
-echo "doing the dual regressions"
-j=0
-for i in $INPUTS ; do
- s=subject`${FSLDIR}/bin/zeropad $j 5`
- echo "$FSLDIR/bin/fsl_glm -i $i -d $ICA_MAPS -o $OUTPUT/dr_stage1_${s}.txt --demean -m $OUTPUT/mask ; \
- $FSLDIR/bin/fsl_glm -i $i -d $OUTPUT/dr_stage1_${s}.txt -o $OUTPUT/dr_stage2_$s --out_z=$OUTPUT/dr_stage2_${s}_Z --demean -m $OUTPUT/mask $DES_NORM ; \
- $FSLDIR/bin/fslsplit $OUTPUT/dr_stage2_$s $OUTPUT/dr_stage2_${s}_ic" >> ${LOGDIR}/drC
- j=`echo "$j 1 + p" | dc -`
-done
-ID_drC=`$FSLDIR/bin/fsl_sub -j $ID_drB -T 30 -N drC -l $LOGDIR -t ${LOGDIR}/drC`
-
-echo "sorting maps and running randomise"
-j=0
-Nics=`$FSLDIR/bin/fslnvols $ICA_MAPS`
-while [ $j -lt $Nics ] ; do
- jj=`$FSLDIR/bin/zeropad $j 4`
-
- RAND=""
- if [ $NPERM -eq 1 ] ; then
- RAND="$FSLDIR/bin/randomise -i $OUTPUT/dr_stage2_ic$jj -o $OUTPUT/dr_stage3_ic$jj -m $OUTPUT/mask $DESIGN -n 1 -V -R"
- fi
- if [ $NPERM -gt 1 ] ; then
- # EDIT HERE
- RAND="$FSLDIR/bin/randomise -i $OUTPUT/dr_stage2_ic$jj -o $OUTPUT/dr_stage3_ic$jj -m $OUTPUT/mask $DESIGN -n $NPERM -T -V"
- fi
-
- echo "$FSLDIR/bin/fslmerge -t $OUTPUT/dr_stage2_ic$jj \`\$FSLDIR/bin/imglob $OUTPUT/dr_stage2_subject*_ic${jj}.*\` ; \
- $FSLDIR/bin/imrm \`\$FSLDIR/bin/imglob $OUTPUT/dr_stage2_subject*_ic${jj}.*\` ; $RAND" >> ${LOGDIR}/drD
- j=`echo "$j 1 + p" | dc -`
-done
-ID_drD=`$FSLDIR/bin/fsl_sub -j $ID_drC -T 60 -N drD -l $LOGDIR -t ${LOGDIR}/drD`
-
diff --git a/fsl_dualreg.cc b/fsl_dualreg.cc
deleted file mode 100644
index bc0fb918f2217522e8e5f8d25f6772b86ec58428..0000000000000000000000000000000000000000
--- a/fsl_dualreg.cc
+++ /dev/null
@@ -1,254 +0,0 @@
-/* test.cc
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-20013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "libvis/miscplot.h"
-#include "miscmaths/miscmaths.h"
-#include "miscmaths/miscprob.h"
-#include "utils/options.h"
-#include <vector>
-#include <ctime>
-#include "newimage/newimageall.h"
-#include "melhlprfns.h"
-#include <iostream>
-
-#ifdef __APPLE__
-#include <mach/mach.h>
-#define memmsg(msg) { \
- struct task_basic_info t_info; \
- mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
- if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
- { \
- cout << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
- cout.flush(); \
- } \
-}
-#else
-#define memmsg(msg) { \
- cout << msg; \
-}
-#endif
-
-// a simple message macro that takes care of cout and log
-#define message(msg) { \
- cout << msg; \
- cout.flush(); \
-}
-
-#define outMsize(msg,Mat) { \
- cerr << " " << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
-}
-
-
-using namespace MISCPLOT;
-using namespace MISCMATHS;
-using namespace Utilities;
-using namespace std;
-using namespace Melodic;
-
-// GLOBALS
-
-clock_t tictime;
-
-
-// The two strings below specify the title and example usage that is
-// printed out as the help or usage message
-
- string title=string("fsl_BLAH")+
- string("\nAuthor: Christian F. Beckmann \nCopyright(c) 2008-2013 University of Oxford\n")+
- string(" \n \n")+
- string(" \n");
- string examples="fsl_BLAH [options]";
-
-//Command line Options {
- Option<string> fnin(string("-i,--in"), string(""),
- string("input file name (matrix 3D or 4D image)"),
- false, requires_argument);
- Option<string> fnmask(string("-m"), string(""),
- string("mask file name "),
- false, requires_argument);
- Option<int> help(string("-h,--help"), 0,
- string("display this help text"),
- false,no_argument);
- Option<int> xdim(string("-x,--xdim"), 0,
- string("xdim"),
- false,requires_argument);
- Option<int> ydim(string("-y,--ydim"), 0,
- string("ydim"),
- false,requires_argument);
- Option<int> econ(string("-e,--econ"), 0,
- string("econ: how to liump stuff"),
- false,requires_argument);
- /*
-}
-*/
-////////////////////////////////////////////////////////////////////////////
-
-// Local functions
-
-void tic(){
- tictime = clock();
-}
-
-void toc(){
- cerr << endl << "TOC: " << float(clock()-tictime)/CLOCKS_PER_SEC << " seconds" << endl<<endl;
-}
-
-Matrix calccorr(const Matrix& in, int econ)
- {
- Matrix Res;
- int nrows=in.Nrows();
- int ncols=in.Ncols();
- Res = zeros(nrows,nrows);
-
- if(econ>0){
- RowVector colmeans(ncols);
- for (int n=1; n<=ncols; n++) {
- colmeans(n)=0;
- for (int m=1; m<=nrows; m++) {
- colmeans(n)+=in(m,n);
- }
- colmeans(n)/=nrows;
- }
- int dcol = econ;
- Matrix suba;
-
- for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
- suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
- int scolmax = suba.Ncols();
-
- for (int n=1; n<=scolmax; n++) {
- double cmean=colmeans(ctr + n - 1);
- for (int m=1; m<=nrows; m++) {
- suba(m,n)-=cmean;
- }
- }
-
- Res += suba*suba.t() / ncols;
- }
- }
- else
- Res = cov(in.t());
- return Res;
- } //Matrix calccorr
-
-void convert_to_mat(volume4D<float>& in, volume<float>& mask, Matrix& out)
-{
- out = in[0].vec(mask).t();
- in.deletevolume(0);
- while(in.tsize()>0){
- out &= in[0].vec(mask).t();
- in.deletevolume(0);
- }
-}
-
-int do_work(int argc, char* argv[]) {
-
-
- tic();
- Matrix MatrixData;
- volume<float> Mean;
-
- if(xdim.value()==0 && ydim.value()==0)
- {
- volume4D<float> RawData;
- volume<float> theMask;
- toc();
- //read data
- message("Reading data file " << (string)fnin.value() << " ... ");
- read_volume4D(RawData,fnin.value());
- message(" done" << endl);
-
- Mean = meanvol(RawData);
- toc();
- message("Reading mask file " << (string)fnmask.value() << " ... ");
- read_volume(theMask,fnmask.value());
-
- MatrixData = RawData.matrix();
-
- memmsg("start");
- Matrix res;
- convert_to_mat(RawData,theMask, res);
-
- memmsg(" Before reshape ");
-
- write_ascii_matrix("mat1", res);
- write_ascii_matrix("mat2", MatrixData);
-
-
- }
- else{
- Matrix data = unifrnd(xdim.value(),ydim.value());
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
- data = unifrnd(10,100);
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
- data = unifrnd(100,1000);
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
- data = unifrnd(100,10000);
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
- data = unifrnd(300,200000);
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
- data = unifrnd(500,20000);
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
- data = unifrnd(500,200000);
- outMsize("data", data);
- tic(); calccorr(data,econ.value()); toc();
-
-
- }
-
-
- return 0;
-}
-
-////////////////////////////////////////////////////////////////////////////
-
- int main(int argc,char *argv[]){
- Tracer tr("main");
- OptionParser options(title, examples);
- try{
- // must include all wanted options here (the order determines how
- // the help message is printed)
-
- options.add(fnin);
- options.add(fnmask);
- options.add(help);
- options.add(xdim);
- options.add(ydim);
- options.add(econ);
- options.parse_command_line(argc, argv);
-
- // line below stops the program if the help was requested or
- // a compulsory option was not set
- if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
- options.usage();
- exit(EXIT_FAILURE);
- }else{
- // Call the local functions
- return do_work(argc,argv);
- }
- }catch(X_OptionError& e) {
- options.usage();
- cerr << endl << e.what() << endl;
- exit(EXIT_FAILURE);
- }catch(std::exception &e) {
- cerr << e.what() << endl;
- }
- }
-
diff --git a/meldata.cc b/meldata.cc
index e1fd9b4a8d3e1175f7eef3558e92b09c9d1a02c7..fc1ee3e8d79f5196ce679fff2abd98310aef0d0d 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -56,6 +56,16 @@ namespace Melodic{
tmpData = RawData.matrix(Mask);
memmsg(" after reshape ");
}
+
+ // If a time series model design was specified, check
+ // that the data dimensions match the model dimensions
+ if (Tdes.Storage() && (tmpData.Nrows() != Tdes.Nrows())) {
+
+ cerr << "ERROR: " << fname << " " <<
+ "- data dimensions (" << tmpData.Nrows() << ") " <<
+ "do not match model dimensions (" << Tdes.Nrows() << ")" << endl;
+ exit(2);
+ }
//convert to percent BOLD signal change
if(opts.pbsc.value()){
@@ -79,7 +89,7 @@ namespace Melodic{
if(opts.remove_meantc.value()){
remmean(tmpData,meanC,2);
}
-
+
//convert to power spectra
if(opts.pspec.value()){
message(" Converting data to powerspectra ...");
@@ -115,48 +125,18 @@ namespace Melodic{
message(" done" << endl);
}
- //INSTACORRS
- Matrix tmpTC;
- tmpTC = tmpData * insta_maps.t();
-
- if(opts.insta_fn.value().length()>0){
- dbgmsg(string("BEGIN: INSTACORR") << endl);
- if(opts.insta_idx.value()<1 || opts.insta_idx.value()>tmpTC.Ncols()){
- cerr << "ERROR:: INSTACORR index is wrong \n\n";
- exit(2);
- }
-
- Matrix tmpRef = tmpTC.Column(opts.insta_idx.value());
- if(opts.insta_idx.value()>1){
- // swap columns
- dbgmsg(string("INSTACORR: swap columns") << endl);
- tmpTC.Column(opts.insta_idx.value()) << tmpTC.Column(1);
- tmpTC.Column(1) << tmpRef;
- }
-
- if(opts.insta_partial.value() && tmpTC.Ncols()>1){
- // partal correlations
- dbgmsg(string("INSTACORR: partial analysis") << endl);
- Matrix tmpConf = tmpTC.Columns(2,tmpTC.Ncols());
- tmpData -= tmpConf * (pinv(tmpConf) * tmpData);
- tmpRef -= tmpConf * (pinv(tmpConf) * tmpRef);
- }
-
- if(opts.insta_varnorm.value()){
- dbgmsg(string("INSTACORR: varnorm") << endl);
- Matrix vscales = pow(stdev(tmpData,1),-1);
- varnorm(tmpData,vscales);
- varnorm(tmpRef,pow(stdev(tmpRef,1),-1));
- }
+ //convert to instacorrs
+ if(opts.insta_fn.value()>""){
+ Matrix vscales = pow(stdev(tmpData,1),-1);
+ varnorm(tmpData,vscales);
- // Shur product
- dbgmsg(string("INSTACORR: SP") << endl);
- SP4(tmpData,tmpRef);
+ Matrix tmpTC = tmpData * insta_mask.t();
+ varnorm(tmpTC,pow(stdev(tmpTC),-1));
- dbgmsg(string("END: INSTACORR") << endl);
- }
- //END INSTACORRS
+ for(int ctr=1; ctr <=tmpData.Ncols();ctr++)
+ tmpData.Column(ctr) = SP(tmpData.Column(ctr),tmpTC);
+ }
tmpData.Release();
dbgmsg(string("END: process_file") << endl);
@@ -650,6 +630,19 @@ namespace Melodic{
create_RXweight();
}
+ //set up instacorr mask image
+ if(opts.insta_fn.value()>""){
+ dbgmsg(string(" Setting up instacorr mask") << endl);
+ volume4D<float> tmp_im;
+ read_volume4D(tmp_im,opts.insta_fn.value());
+
+ if(!samesize(Mean,tmp_im[0])){
+ cerr << "ERROR:: instacorr mask and data have different voxel dimensions \n\n";
+ exit(2);
+ }
+ insta_mask = tmp_im.matrix(Mask);
+ }
+
//seed the random number generator
double tmptime = time(NULL);
if ( opts.seed.value() != -1 ) {
@@ -687,8 +680,20 @@ namespace Melodic{
TconF = read_ascii_matrix(opts.fn_TconF.value());
if(opts.fn_SconF.value().length()>0)
SconF = read_ascii_matrix(opts.fn_SconF.value());
-
- if(numfiles>1 && Sdes.Storage() == 0){
+
+ // Check that the number of input
+ // files matches the session design
+ if (Sdes.Storage()) {
+ if (Sdes.Nrows() != numfiles) {
+ cerr << "ERROR: Number of input files (" << numfiles << ") " <<
+ "does not match subject/session design (" << Sdes.Nrows() << ")" << endl;
+ exit(2);
+ }
+ }
+
+ // Or create a default session design
+ // if one was not specified
+ else if(numfiles>1){
Sdes = ones(numfiles,1);
if(Scon.Storage() == 0){
Scon = ones(1,1);
@@ -697,20 +702,6 @@ namespace Melodic{
}
remmean(Tdes);
- //INSTACORRS
- if(opts.insta_fn.value().length()>0){
- message(" Reading in " << opts.insta_fn.value()
- << " for instantaneous correlation analysis" <<endl);
- volume4D<float> tmp_im;
- read_volume4D(tmp_im,opts.insta_fn.value());
-
- if(!samesize(Mean,tmp_im[0])){
- cerr << "ERROR:: instacorr mask and data have different voxel dimensions \n\n";
- exit(2);
- }
- insta_maps = tmp_im.matrix(Mask);
- }
-
dbgmsg(string("END: setup_misc") << endl);
}
diff --git a/melgmix.cc b/melgmix.cc
deleted file mode 100644
index c8a8707570c6ee5b793fb820f59fd33383a61252..0000000000000000000000000000000000000000
--- a/melgmix.cc
+++ /dev/null
@@ -1,687 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melgmix.cc - Gaussian Mixture Model
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "newimage/newimageall.h"
-//#include "melmmopts.h"
-#include "melgmix.h"
-//#include "melmm.h"
-#include "utils/log.h"
-#include "miscmaths/miscprob.h"
-#include <time.h>
-#include "libvis/miscplot.h"
-#include "libvis/miscpic.h"
-
-using namespace Utilities;
-using namespace NEWIMAGE;
-using namespace MISCPLOT;
-using namespace MISCPIC;
-
-string float2str(float f,int width, int prec, bool scientif){
- ostringstream os;
- int redw = int(std::abs(std::log10(std::abs(f))))+1;
- if(width>0)
- os.width(width);
- if(scientif)
- os.setf(ios::scientific);
- os.precision(redw+std::abs(prec));
- os.setf(ios::internal, ios::adjustfield);
- os << f;
- return os.str();
-}
-
-namespace Melodic{
-
- void MelGMix::setup(const RowVector& dat,
- const string dirname,
- int cnum, volume<float> themask,
- volume<float> themean,
- int num_mix, float eps, bool fixit){
- cnumber = cnum;
- Mask = themask;
- Mean = themean;
- prefix = string("IC_")+num2str(cnum);
-
- fitted = false;
- nummix = num_mix;
- numdata = dat.Ncols();
-
- //normalise data
- datamean = mean(dat,2).AsScalar();
- datastdev= stdev(dat,2).AsScalar();
- data=(dat - datamean)/datastdev;
-
- dbgmsg(" mapdat; mean: " << datamean << " std: " <<datastdev << endl);
-
- props=zeros(1,nummix);
- vars=zeros(1,nummix);
- means=zeros(1,nummix);
- Params=zeros(1,nummix);
- logprobY = 1.0;
-
- props = std::pow(float(nummix),float(-1.0));
-
- Matrix tmp1;
- tmp1 = data * data.t() / numdata;
- vars = tmp1.AsScalar();
-
- float Dmin, Dmax, IntSize;
- Dmin = min(data).AsScalar(); Dmax = max(data).AsScalar();
- IntSize = Dmax / nummix;
-
- means(1)=mean(data,2).AsScalar();
- for (int ctr=2; ctr <= means.Ncols(); ctr++){
- means(ctr) = Dmax - (ctr - 1.5) * IntSize;
- }
- means(2)=means(1)+2*sqrt(vars(1));
- if(nummix>2)
- means(3)=means(1)-2*sqrt(vars(1));
-
- epsilon = eps;
- if(epsilon >=0 && epsilon < 0.0000001)
- {epsilon = std::log(float(dat.Ncols()))/1000 ;}
- fixdim = fixit;
- dbgmsg(" parameters; " << means << " : " << vars << " : " << props << endl);
- }
-
- Matrix MelGMix::threshold(const RowVector& dat, string levels){
- Matrix Res;
- Res = 1.0;
- string tmpstr;
- tmpstr = levels;
- //cerr << " Levels : " << levels << endl << endl;
- Matrix levls(1,4);
- levls = 0;
- Matrix fpr;
- Matrix fdr;
- Matrix nht;
-
- char *p;
- char t[1024];
- const char *discard = ", [];{(})abceghijklmoqstuvwxyzABCEGHIJKLMNOQSTUVWXYZ~!@#$%^&*_-=+|\':></?";
-
- strcpy(t, tmpstr.c_str());
- p=strtok(t,discard);
- while(p){
- Matrix New(1,1);
- New(1,1) = atof(p);
- if(strchr(p,'d')){
- levls(1,3)++;
- if(fdr.Storage()>0){
- fdr = fdr | New;
- }else{
- fdr = New;
- }
- }else{
- if(strchr(p,'p')){
- levls(1,2)++;
- if(fpr.Storage()>0){
- fpr = fpr | New;
- }else{
- fpr = New;
- }
- }else{
- if(strchr(p,'n')){
- levls(1,4)++;
- if(nht.Storage()>0){
- nht = nht | New;
- }else{
- nht = New;
- }
- }else{
- levls(1,1)++;
- levls = levls | New;
- }
- }
- }
- p=strtok(NULL,discard);
- }
-
- if(fpr.Storage()>0){levls = levls | fpr;}
- if(fdr.Storage()>0){levls = levls | fdr;}
- if(nht.Storage()>0){levls = levls | nht;}
-
- // cerr << " levles : " << levls << endl << endl;
- Res = threshold(data, levls);
- set_threshmaps(Res);
-
- return Res;
- }
-
- Matrix MelGMix::threshold(const RowVector& dat, Matrix& levels){
- Matrix tests;
- tests=levels;
- Matrix Nprobs;
-
- //if only single Gaussian: resort to nht
- if(nummix==1||props(1)>0.999||probmap.Sum()<0.05){
- if(levels(1,4)==0){
- Matrix New(1,6);
- New(1,5)=0.05;
- New(1,6)=0.01;
- New(1,4)=2;New(1,1)=0;New(1,2)=0;New(1,3)=0;tests=New;
- }else{
- Matrix New;
- New = levels.Columns(int(1+levels(1,1)+levels(1,2)
- +levels(1,3)),levels.Ncols());
- New(1,4) = levels(1,4);
- New(1,1)=0;New(1,1)=0;New(1,3)=0;
- tests=New;
- }
- }
-
- int numtests = int(tests(1,1)+tests(1,2)+tests(1,3)+tests(1,4));
- Matrix Res(numtests,numdata);
- Res = 0.0;
- int next = 1;
-
- for(int ctr1=1;ctr1<=tests(1,1);ctr1++){
- if(4+next <= tests.Ncols()){
- message(" alternative hypothesis test at p > " << tests(1,4+next) << endl);
- add_infstr(" alternative hypothesis test at p > "+float2str(tests(1,4+next),0,2,false));
- Matrix tmpNull;
- tmpNull = dat;
-/*
- float cutoffpos, cutoffneg;
- cutoffpos = means(1)+100*std::sqrt(vars(1)+0.0000001);
- cutoffneg = means(1)-100*std::sqrt(vars(1)+0.0000001);
-
- for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
- if( probmap(ctr) > tests(1,4+next) ){
- if( dat(ctr) > means(1) )
- cutoffpos = std::min(cutoffpos, float(dat(ctr)));
- else
- cutoffneg = std::max(cutoffneg, float(dat(ctr)));
- }
-
- for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
- if( (dat(ctr) > cutoffneg) && (dat(ctr) < cutoffpos) )
- tmpNull(1,ctr)=0.0;
-*/
- for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
- if( probmap(ctr) < tests(1,4+next) ){
- tmpNull(1,ctr)=0.0;
- }
-
- Res.Row(next) << tmpNull;
- }
- next++;
- }
-
- for(int ctr1=1;ctr1<=tests(1,2);ctr1++){
- if(4+next <=tests.Ncols()){
- cerr << " false positives control " << tests(1,4+next)<<endl;
- Matrix tmp1;
- tmp1 = normcdf(dat,means(1),vars(1));
- Matrix tmpNull;
- tmpNull = dat;
- for(int ctr=1; ctr<=tmp1.Ncols(); ctr++)
- if(tmp1(1,ctr) < tests(1,4+next))
- tmpNull(1,ctr)=0.0;
- Res.Row(next) << tmpNull;
- }
- next++;
- }
-
- for(int ctr1=1;ctr1<=tests(1,3);ctr1++){
- if(4+next <=tests.Ncols()){
- message(" Local False Discovery Rate control at p < " << tests(1,4+next) << endl);
- add_infstr(" Local False Discovery Rate control at p < "+float2str(tests(1,4+next),0,2,false));
- RowVector tmp=dat;
- SortAscending(tmp);
- RowVector newcdf(tmp);
- newcdf << normcdf(tmp,means(1),vars(1));
-
- float thrp = tmp(tmp.Ncols())+0.01;
- float thrn = tmp(1)-0.01;
- int ctr=tmp.Ncols();
- do{
- thrp = tmp(ctr);
- ctr-=1;
- }while(ctr>0 && ( (1.0-newcdf(ctr))*tmp.Ncols() < (tests(1,4+next)*(tmp.Ncols()-ctr+1)) ));
-
- ctr=1;
- do{
- thrn = tmp(ctr);
- ctr+=1;
- }while(ctr<=tmp.Ncols() && ( (newcdf(ctr))*tmp.Ncols() < (tests(1,4+next)*ctr)));
-
- tmp = dat;
- for(ctr=1; ctr<=tmp.Ncols();ctr++)
- if((tmp(ctr) < thrp)&&(tmp(ctr) > thrn))
- tmp(ctr) = 0.0;
- Res.Row(next) << tmp;
- }
- next++;
- }
-
- for(int ctr1=1;ctr1<=tests(1,4);ctr1++){
- if(4+next <=tests.Ncols()){
- message(" 2-sided null hypothesis test at " << tests(1,4+next)<<endl);
- add_infstr(" 2-sided null hypothesis test at "+float2str(tests(1,4+next),0,2,false));
- double mu, sig;
- mu = dat.Sum()/numdata;
- sig = var(dat,2).AsScalar();
- Matrix tmp1;
- tmp1 = normcdf(dat,mu,std::abs(sig));
- Matrix tmpNull;
- tmpNull = dat;
- for(int ctr=1; ctr<=tmp1.Ncols(); ctr++)
- if((tmp1(1,ctr) < 1-0.5*(tests(1,4+next))&&
- (tmp1(1,ctr) > 0.5*(tests(1,4+next)))))
- tmpNull(1,ctr)=0.0;
- Res.Row(next) << tmpNull;
- }
- next++;
- }
-
- return Res;
- }
-
- /* GMM fitting */
-
- void MelGMix::gmmupdate(){
- int it_ctr = 1;
- bool exitloop = false;
- float oldll;
-
- Matrix tmp0;Matrix tmp1;Matrix prob_K__y_theta;
- Matrix kdata;
- RowVector prob_Y__theta;RowVector Nbar;
- RowVector mubar;RowVector sigmabar;RowVector pibar;
-
- do{
- oldll = logprobY;
-
- //make sure all variances are acceptable
- for(int ctr1=1; ctr1 <=vars.Ncols(); ctr1++)
- if(vars(ctr1)<0.0001){
- vars(ctr1) = 0.0001;
- }
-
- tmp0 = normpdf(data,means,vars);
- tmp1 = SP(props.t()*ones(1,numdata),tmp0);
- prob_Y__theta << sum(tmp1,1);
- logprobY = log(prob_Y__theta).Sum();
- prob_K__y_theta = SP(tmp1,pow(ones(nummix,1)*prob_Y__theta,-1));
- Nbar << sum(prob_K__y_theta,2).t();
- pibar = Nbar / numdata;
- kdata = ones(nummix,1)*data;
- mubar <<SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));
- kdata -= mubar.t()*ones(1,numdata);
- kdata = pow(kdata,2);
- sigmabar << SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));
-
- means = mubar;
- vars = sigmabar;
- props = pibar;
-
- if(epsilon<0){exitloop = it_ctr >= -epsilon;}
- else{exitloop = (((logprobY-oldll < epsilon)&&(it_ctr>20))
- ||(it_ctr>400));}
- it_ctr++;
- }while(!exitloop);
- }
-
- void MelGMix::gmmfit(){
- int i,j;
-
- if(fixdim){
- if(nummix>1){
- gmmupdate();
- add_params(means,vars,props,logprobY,MDL,Evi,true);
- }else{
- means.ReSize(1);
- means = data.Sum()/numdata;
- vars.ReSize(1);
- vars = var(data,2);
- props.ReSize(1);
- props = 1.0;
- gmmevidence();
- }
- }else{
- RowVector Score(Params.Ncols());
- do{
- gmmupdate();
- Score(nummix) = gmmevidence();
- int idx1,idx2;
- RowVector pitmp = props;
-
- pitmp.MaximumAbsoluteValue1(idx1);
- pitmp(idx1)=0.0;
- pitmp.MaximumAbsoluteValue1(idx2);
-
- if(props.MaximumAbsoluteValue1(i)<0.9){
- if((vars(idx2)>0.15)&&
- (std::abs(means(idx2)-means(idx1))<0.5*vars(idx1))){
- Score(nummix) = Score(nummix)/(2*(means(idx1)));
- }
- }
-
- add_params(means,vars,props,logprobY,MDL,Evi,true);
-
- gmmreducemm();
- means = means.Columns(1,nummix);
- vars = vars.Columns(1,nummix);
- props = props.Columns(1,nummix);
-
- }while(nummix>1);
- means.ReSize(1);
- means = data.Sum()/numdata;
- vars.ReSize(1);
- vars = var(data,2);
- props.ReSize(1);
- props = 1.0;
- Score(nummix) = gmmevidence();
- add_params(means,vars,props,logprobY,MDL,Evi,true);
- //identify best MM
- Score.MinimumAbsoluteValue2(i,j);
- means.ReSize(j);
- vars.ReSize(j);
- props.ReSize(j);
- nummix = j;
- int index; index = 3 + (j-1)*5;
- means = Params.SubMatrix(index,index,1,j);
- vars = Params.SubMatrix(index+1,index+1,1,j);
- props = Params.SubMatrix(index+2,index+2,1,j);
- }
-
- props.MaximumAbsoluteValue2(i,j);
- if(j>1){
- float tmp;
- tmp = means(1);means(1) = means(j);means(j)=tmp;
- tmp = vars(1);vars(1) = vars(j);vars(j)=tmp;
- tmp = props(1);props(1) = props(j);props(j)=tmp;
- }
-
- add_params(means,vars,props,logprobY,MDL,Evi,true);
-
- if(nummix==1)
- probmap << normcdf(data,means(1),vars(1));
- else{
- Matrix Nprobs;
- Matrix tmp0;
- tmp0 = normpdf(data,means,vars);
- Nprobs = SP(props.t()*ones(1,numdata),tmp0);
- tmp0 = ones(Nprobs.Nrows(),1)*pow(sum(Nprobs,1),-1);
- Nprobs = SP(tmp0,Nprobs);
- probmap << SP(sum(Nprobs.Rows(2,Nprobs.Nrows()),1),
- pow(sum(Nprobs,1),-1));
- }
- }
-
- float MelGMix::gmmevidence(){
- Matrix tmp0;
- if(means.Ncols()>1){
- tmp0 = normpdf(data,means,vars);
- }else{
- tmp0 = normpdf(data,means.AsScalar(),vars.AsScalar());
- }
- Matrix tmp1;
- tmp1 = SP(props.t()*ones(1,numdata),tmp0);
- tmp0 = SP(tmp0,pow(ones(nummix,1)*sum(tmp1,1),-1));
- tmp0 = pow(tmp0-ones(nummix,1)*tmp0.Row(nummix),2);
- float logH = 0;
- if(means.Ncols()>1){
- logH = sum(log(sum(tmp0.Rows(1,nummix-1),2)),1).AsScalar();
- }
- logH = logH + 2*sum(log(std::sqrt(2.0)*numdata*props),2).AsScalar();
- logH = logH - 2*sum(props,2).AsScalar();
-
- RowVector prob_Y__theta;
- prob_Y__theta << sum(tmp1,1);
- logprobY = log(prob_Y__theta).Sum();
- MDL = -logprobY + (1.5*nummix-0.5)* std::log(float(numdata));
- Evi = -logprobY +nummix*std::log(2.0)-std::log(MISCMATHS::gamma(nummix))
- -3*nummix/2*std::log(M_PI)+0.5*logH;
-
- return Evi;
- }
-
- void MelGMix::gmmreducemm(){
- Matrix dlm(nummix,nummix);
- Matrix mus(nummix,nummix);
- Matrix rs(nummix,nummix);
-
- for(int ctri=1;ctri<=nummix; ctri++){
- for(int ctrj=1;ctrj<=nummix; ctrj++){
- mus(ctri,ctrj) = (props(ctri)*means(ctri)+props(ctrj)*means(ctrj))
- /(props(ctri)+props(ctrj));
- rs(ctri,ctrj) = (props(ctri)*(vars(ctri)+
- std::pow(means(ctri)-mus(ctri,ctrj),2) )
- + props(ctrj)*(vars(ctrj)
- + std::pow(means(ctrj)-mus(ctri,ctrj),2)))
- / (props(ctri)+props(ctrj));
- dlm(ctri,ctrj) = 0.5*numdata*(
- props(ctri)*std::log(
- std::abs(rs(ctri,ctrj))/std::abs(vars(ctri)))
- + props(ctrj)*std::log(std::abs(rs(ctri,ctrj))
- / std::abs(vars(ctrj))));
- }
- }
-
- dlm += IdentityMatrix(nummix)*dlm.Maximum();
-
- int i,j;
- float val;
- val=dlm.MinimumAbsoluteValue2(i,j);
-
- nummix--;
-
- RowVector newmean(nummix);
- RowVector newvars(nummix);
- RowVector newprop(nummix);
-
- int ctr0 = 1;
- for(int ctr=1; ctr<=nummix+1; ctr++){
- if(ctr!=i&&ctr!=j){
- newmean(ctr0) = means(ctr);
- newvars(ctr0) = vars(ctr);
- newprop(ctr0) = props(ctr);
- ctr0++;
- }
- }
- //cerr << "ctr0 " << ctr0 << endl;
- if(ctr0<=nummix){
- newmean(ctr0) = mus(i,j);
- newvars(ctr0) = rs(i,j);
- newprop(ctr0) = props(i)+props(j);
-
- means = newmean;
- vars=newvars;
- props=newprop;
- }
- }
-
- void MelGMix::ggmfit(){
- // fit a mixture of a Gaussian and multiple Gamma functions to the histogram
-
- float log_p_y_theta = 1.0;
- float old_ll = 2.0;
- float g_eps = 0.000001;
- int it_ctr = 0;
- double Dmax, Dmin;
-
- Dmax = 2 * data.Maximum();
- Dmin = -2 * data.Minimum();
-
- //resize means, vars and props
- if(nummix > 3)
- nummix = 3;
- means = means.Columns(1,nummix);
- vars = vars.Columns(1,nummix);
- props = props.Columns(1,nummix);
-
- means(1) = -2*mean(data,2).AsScalar();
-
- Matrix tmp1;Matrix prob_K__y_theta;
- Matrix kdata;
- RowVector prob_Y__theta;RowVector Nbar;
- RowVector mubar;RowVector sigmabar;RowVector pibar;
- Matrix p_ygx(nummix,numdata);
- offset = 0.0;
- float const2;
- Matrix negdata(data);
- negdata = -data;
-
- while((it_ctr<30) ||
- ((std::abs(old_ll - log_p_y_theta)>g_eps) && (it_ctr<500))){ // fit GGM
- it_ctr++;
- //offset = (std::min(0.2,1-props(1)))*std::sqrt(vars(1));
-
- //make sure all variances are acceptable
- for(int ctr1=1; ctr1 <=nummix; ctr1++)
- if(vars(ctr1)<0.0001){
- vars(ctr1) = 0.0001;
- }
-
- p_ygx = 0.0;
- p_ygx.Row(1) << normpdf(data,means(1),vars(1));
-
- const2 = (2.6-props(1))*sqrt(vars(1))+means(1); //min. nht level
-
- means(2) = (std::max(means(2), std::max(0.001,
- 0.5 * ( const2 + std::sqrt( const2*const2 + 4*vars(2))))));
- vars(2) = std::max(std::min(vars(2), 0.5*std::pow(means(2),2)),0.0001);
- p_ygx.Row(2) << gammapdf(data,means(2),vars(2));
-
- if(nummix>2){
- const2 = (2.6-props(1))*sqrt(vars(1))-means(1);
-
- means(3) = -(std::max(-means(3), std::max(0.001,
- 0.5 * ( const2 + std::sqrt( const2*const2 + 4*vars(3))))));
- vars(3) = std::max(std::min(vars(3), 0.5*std::pow(means(3),2)),0.0001);
- p_ygx.Row(3) << gammapdf(negdata,-means(3),vars(3));
- }
-
- tmp1 = SP(props.t()*ones(1,numdata),p_ygx);
- prob_Y__theta << sum(tmp1,1);
-
- //deal with non-modelled voxels
- for(int ctr=1; ctr<=tmp1.Ncols(); ctr++)
- if(prob_Y__theta(ctr) < 0.0001)
- prob_Y__theta(ctr) = 0.0001;
-
- old_ll = log_p_y_theta;
- log_p_y_theta = log(prob_Y__theta).Sum();
- if((it_ctr<30) ||
- ((std::abs(old_ll - log_p_y_theta)>g_eps) && (it_ctr<300))){//update
-
- prob_K__y_theta = SP(tmp1,pow(ones(nummix,1)*prob_Y__theta,-1));
- Nbar << sum(prob_K__y_theta,2).t();
- for(int ctr=1; ctr<=Nbar.Ncols(); ctr++)
- if(Nbar(ctr) < 0.0001 * numdata)
- Nbar = Nbar + 0.0001;
- pibar= Nbar / numdata;
- // cerr << "pibar :" << pibar << endl;
- kdata = ones(nummix,1)*data;
- mubar <<SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));
- // cerr << "mubar :" << mubar << endl;
-
- kdata -= mubar.t()*ones(1,numdata);
- kdata = pow(kdata,2);
- sigmabar << SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));
-
- means = mubar;
- vars = sigmabar;
- props = pibar;
- }//update
- } //while loop
-
- props = props / sum(props,2).AsScalar();
- add_params(means,vars,props,logprobY,MDL,Evi,true);
-
- probmap << SP(sum(tmp1.Rows(2,tmp1.Nrows()),1),
- pow(sum(tmp1,1),-1));
-
- dbgmsg(" mu: " << means << " sig: " << vars << " prop: " << props << endl);
-
- if(props(1)<0.4 ){
- //set up GMM
- message(" try Gaussian Mixture Model " << endl);
- props=zeros(1,nummix);
- vars=zeros(1,nummix);
- means=zeros(1,nummix);
- Params=zeros(1,nummix);
- logprobY = 1.0;
- props = std::pow(float(nummix),float(-1.0));
-
- tmp1 = data * data.t() / numdata;
- vars = tmp1.AsScalar();
- float Dmin, Dmax, IntSize;
- Dmin = min(data).AsScalar(); Dmax = max(data).AsScalar();
- IntSize = Dmax / nummix;
- means(1)=mean(data,2).AsScalar();
- for (int ctr=2; ctr <= means.Ncols(); ctr++){
- means(ctr) = Dmax - (ctr - 1.5) * IntSize;
- }
- means(2)=means(1)+sqrt(vars(1));
- if(nummix>2)
- means(3)=means(1)-sqrt(vars(1));
-
- fit(string("GMM"));
- }
-
- }
-
- /* INPUT / OUTPUT */
-
- void MelGMix::add_params(Matrix& mu, Matrix& sig, Matrix& pi,
- float logLH, float MDL, float Evi, bool advance){
- int size = Params.Ncols();
- if(size<2){size=2;}
- Matrix New(5,size);
- New = 0;
-
- New.SubMatrix(3,3,1,mu.Ncols())=mu;
- New.SubMatrix(4,4,1,mu.Ncols())=sig;
- New.SubMatrix(5,5,1,mu.Ncols())=pi;
- New(1,1)=nummix;
-
- New(1,2)=-logLH;
- New(2,1)=Evi;
- New(2,2)=MDL;
- if(Params.Storage()>nummix){
- Params = New & Params;
- }else{
- Params = New;
- }
- }
-
- void MelGMix::get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi,
- float logLH, float MDL, float Evi){
-
- }
-
- void MelGMix::save(){
-
- }
-
- void MelGMix::status(const string &txt){
- cerr << txt << "epsilon " << epsilon << endl;
- cerr << txt << "nummix " << nummix << endl;
- cerr << txt << "numdata " << numdata << endl;
- cerr << txt << "means " << means << endl;
- cerr << txt << "vars " << vars << endl;
- cerr << txt << "props " << props << endl;
- //write_ascii_matrix(logger.appendDir(string(txt + "means")),means);
- //write_ascii_matrix(logger.appendDir(string(txt + "vars")),vars);
- //write_ascii_matrix(logger.appendDir(string(txt + "props")),props);
- }
-
- void MelGMix::create_rep(){
-
- }
-
-
-}
-
-
-
diff --git a/melgmix.h b/melgmix.h
deleted file mode 100644
index 357eaae72e1e6fcff3734272572f23fe1b21ed2a..0000000000000000000000000000000000000000
--- a/melgmix.h
+++ /dev/null
@@ -1,207 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melgmix.h - class for Gaussian/Gamma Mixture Model
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#ifndef __MELGMIX_h
-#define __MELGMIX_h
-
-#include "newimage/newimageall.h"
-#include "utils/log.h"
-#include "melodic.h"
-#include "utils/options.h"
-#include "meloptions.h"
-//#include "melreport.h"
-
-using namespace Utilities;
-using namespace NEWIMAGE;
-
-namespace Melodic{
-
- class MelGMix{
- public:
-
- MelGMix(MelodicOptions &popts, Log &plogger):
- opts(popts),
- logger(plogger){}
-
- ~MelGMix() {
- //mainhtml << endl << "<hr></CENTER></BODY></HTML>" << endl;
- }
-
- void save();
-
- void setup(const RowVector& dat, const string dirname,
- int here, volume<float> themask,
- volume<float> themean, int num_mix = 3,
- float eps = 0.0, bool fixdim = false);
-
- void gmmfit();
- void ggmfit();
-
- inline void fit(string mtype = string("GGM")){
- mmtype = mtype;
- if(mmtype==string("GGM"))
- this->ggmfit();
- else
- this->gmmfit();
-
- //re-insert mean and stdev
- data = data*datastdev + datamean;
- //threshmaps = threshmaps*datastdev + datamean;
- means = means*datastdev + datamean;
- vars = vars*datastdev*datastdev;
- }
-
- inline Matrix threshold(string levels){
- return this->threshold(data, levels);
- }
- inline Matrix threshold(RowVector& levels){
- return this->threshold(data, levels);
- }
- Matrix threshold(const RowVector& dat, Matrix& levels);
- Matrix threshold(const RowVector& dat, string levels);
-
- void status(const string &txt);
-
- inline RowVector& get_means() {return means;}
- inline void set_means(RowVector& Arg) {means = Arg;}
-
- inline RowVector& get_vars() {return vars;}
- inline void set_vars(RowVector& Arg) {vars = Arg;}
-
- inline RowVector& get_pi() {return props;}
- inline void set_pi(RowVector& Arg) {props = Arg;}
-
- inline RowVector& get_data() {return data;}
- inline void set_data(RowVector& Arg) {data = Arg;}
-
- inline RowVector& get_prob() {return probmap;}
-
- inline float get_eps() {return epsilon;}
- inline void set_eps(float Arg) {epsilon = Arg;}
-
- inline Matrix& get_threshmaps() {return threshmaps;}
- inline void set_threshmaps(Matrix& Arg) {threshmaps = Arg;}
-
- inline bool isfitted(){return fitted;}
-
- inline int mixtures(){return nummix;}
-
- inline string get_type() { return mmtype;}
- inline void set_type(string Arg) { mmtype = Arg;}
-
- inline string get_prefix() { return prefix;}
- inline void set_prefix(string Arg) { prefix = Arg;}
-
- inline RowVector get_probmap() {return probmap;}
-
- inline float get_offset() {return offset;}
- inline void set_offset(float Arg) {offset = Arg;}
-
- inline void flipres(int num){
- means = -means;
- data = -data;
- threshmaps = -threshmaps;
- if(mmtype=="GGM"){
- float tmp;
- tmp= means(2);means(2)=means(3);means(3)=tmp;
- tmp=vars(2);vars(2)=vars(3);vars(3)=tmp;
- tmp=props(2);props(2)=props(3);props(3)=tmp;
- }
- }
-
- void create_rep();
-
- inline void add_infstr(string what){
- threshinfo.push_back(what);
- }
-
- inline string get_infstr(int num){
- if((threshinfo.size()<(unsigned int)(num-1))||(num<1))
- return string("");
- else
- return threshinfo[num-1];
- }
-
- inline int size_infstr(){
- return threshinfo.size();
- }
-
- inline void clear_infstr(){
- threshinfo.clear();
- }
-
- inline void smooth_probs(float howmuch){
- volume4D<float> tempVol;
- tempVol.setmatrix(probmap,Mask);
- tempVol[0]= smooth(tempVol[0],howmuch);
- probmap = tempVol.matrix(Mask);
- }
-
- inline Matrix get_params(){
- Matrix tmp = zeros(3,means.Ncols());
- tmp.Row(1) = means;
- tmp.Row(2) = vars;
- tmp.Row(3) = props;
- return tmp;
- }
-
- double datamean;
- double datastdev;
-
- private:
- __attribute__((unused)) MelodicOptions &opts;
- __attribute__((unused)) Log &logger; //global log file
-
- //Log mainhtml;
-
- void gmmupdate();
- float gmmevidence();
- void gmmreducemm();
- void add_params(Matrix& mu, Matrix& sig, Matrix& pi,
- float logLH, float MDL, float Evi, bool advance = false);
- void get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi,
- float logLH, float MDL, float Evi);
-
- Matrix Params;
- Matrix threshmaps;
-
- RowVector means;
- RowVector vars;
- RowVector props;
- RowVector data;
- RowVector probmap;
-
- volume<float> Mean;
- volume<float> Mask;
-
- float epsilon;
- float logprobY;
- float MDL;
- float Evi;
- float offset;
-
- int nummix;
- int numdata;
- int cnumber;
-
- bool fitted;
- bool fixdim;
-
- string prefix;
- string mmtype;
- string dirname;
-
- vector<string> threshinfo;
-
- };
-}
-
-#endif
diff --git a/melhlprfns.cc b/melhlprfns.cc
deleted file mode 100644
index 76fa04c4fbd9bc59c5546251819b7032afa1b3cc..0000000000000000000000000000000000000000
--- a/melhlprfns.cc
+++ /dev/null
@@ -1,924 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melhlprfns.cc - misc functions
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "melhlprfns.h"
-#include "libprob.h"
-#include "miscmaths/miscprob.h"
-#include "miscmaths/t2z.h"
-#include "miscmaths/f2z.h"
-
-namespace Melodic{
-
- void update_mask(volume<float>& mask, Matrix& Data)
- {
- Matrix DStDev=stdev(Data);
- volume4D<float> tmpMask, RawData;
- tmpMask.setmatrix(DStDev,mask);
-
- float tMmax;
- volume<float> tmpMask2;
- tmpMask2 = tmpMask[0];
-
- tMmax = tmpMask2.max();
- double st_mean = DStDev.Sum()/DStDev.Ncols();
- double st_std = stdev(DStDev.t()).AsScalar();
-
- volume4D<float> newmask;
- newmask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,(float) 0.01*st_mean),tMmax);
-
- Matrix newmaskM,newData;
- newmaskM = newmask.matrix(mask);
- int N = Data.Nrows();
-
- if(int(newmaskM.Row(1).SumAbsoluteValue() + 0.3) < newmaskM.Ncols()){
- RawData.setmatrix(Data.Row(1),mask);
- newData = RawData.matrix(newmask[0]);
- for(int r=2; r <= N; r++){
- RawData.setmatrix(Data.Row(r),mask);
- newData &= RawData.matrix(newmask[0]);
- }
- Data = newData;
- mask = newmask[0];
- }
- }
-
- void del_vols(volume4D<float>& in, int howmany)
- {
- for(int ctr=1; ctr<=howmany; ctr++){
- in.deletevolume(ctr);
- }
- }
-
- Matrix calc_FFT(const Matrix& Mat, const bool logpwr)
- {
- Matrix res;
- for(int ctr=1; ctr <= Mat.Ncols(); ctr++){
- ColumnVector tmpCol;
- tmpCol=Mat.Column(ctr);
- ColumnVector FtmpCol_real;
- ColumnVector FtmpCol_imag;
- ColumnVector tmpPow;
- if(tmpCol.Nrows()%2 != 0){
- Matrix empty(1,1); empty=0;
- tmpCol &= empty;
- }
- RealFFT(tmpCol,FtmpCol_real,FtmpCol_imag);
- tmpPow = pow(FtmpCol_real,2)+pow(FtmpCol_imag,2);
- tmpPow = tmpPow.Rows(2,tmpPow.Nrows());
- if(logpwr) tmpPow = log(tmpPow);
- if(res.Storage()==0){res= tmpPow;}else{res|=tmpPow;}
- }
- return res;
- } //Matrix calc_FFT()
-
- Matrix smoothColumns(const Matrix& inp)
- {
- Matrix temp(inp);
- int ctr1 = temp.Nrows();
- Matrix temp2(temp);
- temp2=0;
-
- temp = temp.Row(4) & temp.Row(3) & temp.Row(2) & temp & temp.Row(ctr1-1)
- & temp.Row(ctr1-2) &temp.Row(ctr1-3);
-
- double kern[] ={0.0045 , 0.055, 0.25, 0.4, 0.25, 0.055, 0.0045};
- double fac = 0.9090909;
-
-
- for(int cc=1;cc<=temp2.Ncols();cc++){
- for(int cr=1;cr<=temp2.Nrows();cr++){
- temp2(cr,cc) = fac*( kern[0] * temp(cr,cc) + kern[1] * temp(cr+1,cc) +
- kern[2] * temp(cr+2,cc) + kern[3] * temp(cr+3,cc) +
- kern[4] * temp(cr+4,cc) + kern[5] * temp(cr+5,cc) +
- kern[6] * temp(cr+6,cc));
- }
- }
- return temp2;
- } //Matrix smoothColumns()
-
- Matrix convert_to_pbsc(Matrix& inp)
- {
- Matrix meanimg;
- meanimg = mean(inp);
- float eps = 0.00001;
-
- for(int ctr=1; ctr <= inp.Ncols(); ctr++){
- if(meanimg(1,ctr) < eps)
- meanimg(1,ctr) = eps;
- }
-
- for(int ctr=1; ctr <= inp.Nrows(); ctr++){
- Matrix tmp;
- tmp << inp.Row(ctr);
- inp.Row(ctr) << 100 * SP((tmp - meanimg),pow(meanimg,-1));
- }
-
- inp = remmean(inp);
- return meanimg;
- } //void convert_to_pbsc
-
- RowVector varnorm(Matrix& in, int dim, float level, int econ)
- {
- SymmetricMatrix Corr(cov_r(in,false,econ));
- RowVector out;
- out = varnorm(in,Corr,dim,level, econ);
- return out;
- } //RowVector varnorm
-
- void varnorm(Matrix& in, const RowVector& vars)
- {
- for(int ctr=1; ctr <=in.Nrows();ctr++)
- in.Row(ctr) = SD(in.Row(ctr),vars);
- }
-
- RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim, float level, int econ)
- {
-
- Matrix tmpE, white, dewhite;
- RowVector tmpD, tmpD2;
-
- std_pca(remmean(in,2), Corr, tmpE, tmpD, econ);
- calc_white(tmpE,tmpD, dim, white, dewhite);
-
- Matrix ws = white * in;
- for(int ctr1 = 1; ctr1<=ws.Ncols(); ctr1++)
- for(int ctr2 = 1; ctr2<=ws.Nrows(); ctr2++)
- if(std::abs(ws(ctr2,ctr1)) < level)
- ws(ctr2,ctr1)=0;
- tmpD = stdev(in - dewhite*ws);
- for(int ctr1 = 1; ctr1<=tmpD.Ncols(); ctr1++)
- if(tmpD(ctr1) < 0.01){
- tmpD(ctr1) = 1.0;
- in.Column(ctr1) = 0.0*in.Column(ctr1);
- }
- varnorm(in,tmpD);
-
- return tmpD;
- } //RowVector varnorm
-
-
-
- Matrix SP2(const Matrix& in, const Matrix& weights, int econ)
- {
- Matrix Res;
- Res = in;
- if(econ>0){
- if(weights.Ncols() == in.Ncols()){
- ColumnVector tmp;
- for(int ctr=1; ctr <= in.Ncols(); ctr++){
- tmp = in.Column(ctr);
- tmp = tmp * weights(1,ctr);
- Res.Column(ctr) = tmp;
- }
- }
- }
- else{
- Res = ones(in.Nrows(),1)*weights.Row(1);
- Res = SP(in,Res);
- }
- return Res;
- } //Matrix SP2
-
- void SP3(Matrix& in, const Matrix& weights)
- {
- for(int ctr=1; ctr <= in.Nrows(); ctr++){
- in.Row(ctr) << SP(in.Row(ctr),weights.AsRow());
- }
- }
-
- void SP4(Matrix& in, const Matrix& weights){
- RowVector tmp;
- for(int ctr=1; ctr <= in.Nrows(); ctr++){
- tmp = in.Row(ctr);
- tmp = tmp * weights(ctr,1);
- in.Row(ctr) << tmp;
- }
- }
-
- Matrix corrcoef(const Matrix& in1, const Matrix& in2){
- Matrix tmp = in1;
- tmp |= in2;
- Matrix out;
- out=MISCMATHS::corrcoef(tmp,0);
- return out.SubMatrix(1,in1.Ncols(),in1.Ncols()+1,out.Ncols());
- }
-
- Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part){
- Matrix tmp1 = in1, tmp2 = in2, out;
- if(part.Storage()){
- tmp1 = tmp1 - part * pinv(part) * tmp1;
- tmp2 = tmp2 - part * pinv(part) * tmp2;
- }
-
- out = corrcoef(tmp1,tmp2);
- return out;
- }
-
- float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
- {
-
-// tmpD2= tmpD | tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols());
-// cerr << tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols()) << endl;
-
-//
-
-// Matrix tmpPE;
-// tmpPE = SP(param,ones(param.Nrows(),1)*pow(stdev(param,1)*std::sqrt((float)param.Nrows()),-1));
-
-// RE |= tmpPE;
-// RowVector tmpD2;
-// tmpD2 = tmpD | stdev(param,1).AsRow()*std::sqrt((float)param.Nrows());
-// RD << abs(diag(tmpD2.t()));
-// RD << RD.SymSubMatrix(N-dim+1,N);
-
- Matrix RE;
- DiagonalMatrix RD;
- int N = tmpE.Ncols();
- dim = std::min(dim,N);
-
-// cerr << stdev(param) << endl;
- RE = tmpE.Columns(std::min(N-dim+1+param.Ncols(),N-2),N);
- RE |= param;
-
-// cerr << paramS.Nrows() << " x " << paramS.Ncols() << endl;
-// cerr << paramS << endl;
- RowVector tmpD2;
- tmpD2 = tmpD | pow(paramS,2).AsRow();
- RD << abs(diag(tmpD2.t()));
-
-// cerr << " " <<tmpD2.Ncols() << " " << N << " " << dim << endl;
- RD << RD.SymSubMatrix(N-dim+1+param.Ncols(),N+param.Ncols());
-
- float res = 1.0;
- white = sqrt(abs(RD.i()))*RE.t();
- dewhite = RE*sqrt(RD);
-
- if(dim < N)
- res = PercEV(dim);
- return res;
- } //Matrix calc_white
-
- float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite)
- {
- Matrix RE;
- DiagonalMatrix RD;
- int N = tmpE.Ncols();
- dim = std::min(dim,N);
- RE = tmpE.Columns(N-dim+1,N);
- RD << abs(diag(tmpD.t()));
- RD << RD.SymSubMatrix(N-dim+1,N);
-
- float res = 1.0;
- white = sqrt(abs(RD.i()))*RE.t();
- dewhite = RE*sqrt(RD);
-
- if(dim < N)
- res = PercEV(dim);
- return res;
- } //Matrix calc_white
-
- void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite)
- {
- RowVector tmp(tmpE.Ncols());
- float tmp2;
- tmp2 = calc_white(tmpE,tmpD, tmp, dim, white, dewhite);
- } //Matrix calc_white
-
- void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
- {
- RowVector tmp(tmpE.Ncols());
- float tmp2;
- tmp2 = calc_white(tmpE,tmpD, tmp, dim, param, paramS, white, dewhite);
- } //Matrix calc_white
-
- void calc_white(const SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite)
- {
- Matrix RE;
- DiagonalMatrix RD;
- RowVector tmp2;
- EigenValues(Corr,RD,RE);
- tmp2 = diag(RD).t();
- calc_white(RE,tmp2, dim, white, dewhite);
- } //Matrix calc_white
-
-
- void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
- {
- if(weights.Storage()>0)
- Corr = cov_r(Mat, weights, econ);
- else
- Corr = cov_r(Mat,false,econ);
-
- DiagonalMatrix tmpD;
- EigenValues(Corr,tmpD,evecs);
- evals = tmpD.AsRow();
- } //void std_pca
-
- void std_pca(const Matrix& Mat, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
- {
- Matrix weights;
- std_pca(Mat,weights,Corr,evecs,evals, econ);
- } //void std_pca
-
- void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc, int iter)
- {
- Matrix guess;
- guess = normrnd(Mat.Nrows(),num_pc);
- em_pca(Mat,guess,evecs,evals,num_pc,iter);
- } //void em_pca
-
- void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc, int iter)
- {
- Matrix C;
- if(guess.Ncols() < num_pc){
- C=normrnd(Mat.Nrows(),num_pc);
- C.Columns(1,guess.Ncols()) = guess;
- }
- else
- C = guess;
-
- Matrix tmp, tmp2;
- for(int ctr=1; ctr <= iter; ctr++){
- // E-Step
- tmp = C.t()*C;
- tmp = tmp.i();
- tmp = tmp * C.t();
- tmp = tmp * Mat;
- // M-Step
- tmp2 = tmp * tmp.t();
- tmp2 = tmp2.i();
- tmp2 = Mat*tmp.t()*tmp2;
- C = tmp2;
- }
-
- symm_orth(C);
- Matrix Evc;
- SymmetricMatrix tmpC;
- RowVector Evl;
- tmp = C.t() * Mat;
- std_pca(tmp,tmpC,Evc,Evl);
- evals = Evl;
- evecs = C * Evc;
- } //void em_pca
-
- float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
- {
- SymmetricMatrix Corr;
- Matrix Evecs, tmpWM, tmpDWM, tmp;
- RowVector Evals;
- std_pca(Mat.t(), Corr, Evecs, Evals);
- calc_white(Corr, dim, tmpWM, tmpDWM);
- tmp = tmpWM * Mat.t();
- cols = tmp.t();
- rows << tmpDWM;
- float res;
- Evals=fliplr(Evals);
- res = sum(Evals.Columns(1,dim),2).AsScalar()/sum(Evals,2).AsScalar()*100;
- return res;
- } // rankapprox
-
- RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows)
- {
- Matrix out; RowVector res(Mat.Ncols());
- for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++){
- Matrix tmpVals(cols.Nrows(),rows.Nrows());
- for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
- tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() *
- (ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
-
- Matrix tmpcols, tmprows;
- res(ctr1) =rankapprox(tmpVals, tmpcols, tmprows);
- cols.Column(ctr1) = tmpcols;
- rows.Column(ctr1) = tmprows;
- }
- return res;
- } // krfact
-
- RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows)
- {
- RowVector res;
- cols = zeros(colnum,Mat.Ncols());
- rows = zeros(int(Mat.Nrows() / colnum),Mat.Ncols());
- res = krfact(Mat,cols,rows);
- return res;
- } // krfact
-
- Matrix krprod(const Matrix& cols, const Matrix& rows)
- {
- Matrix out;
- out = zeros(cols.Nrows()*rows.Nrows(),cols.Ncols());
- for(int ctr1 = 1; ctr1 <= cols.Ncols(); ctr1++)
- for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
- {
- out.SubMatrix(cols.Nrows()*(ctr2-1)+1,cols.Nrows()*ctr2,ctr1,ctr1) << cols.Column(ctr1) * rows(ctr2,ctr1);
- }
- return out;
- } // krprod
-
- Matrix krapprox(const Matrix& Mat, int size_cols, int dim)
- {
- Matrix out, cols, rows;
- out = zeros(Mat.Nrows(), Mat.Ncols());
- cols = zeros(size_cols,Mat.Ncols());
- rows = zeros(int(Mat.Nrows() / size_cols), Mat.Ncols());
- krfact(Mat,cols,rows);
- out = krprod(cols, rows);
- return out;
- } // krapprox
-
- void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2, RowVector& out3, int& out4, int num_vox, float resels)
- {
- RowVector AdjEV;
- AdjEV << in.AsRow();
- AdjEV = AdjEV.Columns(3,AdjEV.Ncols());
- AdjEV = AdjEV.Reverse();
-
- RowVector CircleLaw;
- int NumVox = (int) floor(num_vox/(2.5*resels));
-
- CircleLaw = Feta(int(AdjEV.Ncols()), NumVox);
-
- for(int ctr=1;ctr<=CircleLaw.Ncols(); ctr++){
- if(CircleLaw(ctr)<5*10e-10){CircleLaw(ctr) = 5*10e-10;}
- }
-
- /* float slope;
- slope = CircleLaw.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
- int(AdjEV.Ncols()/4)).Sum() /
- AdjEV.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
- int(AdjEV.Ncols()/4)).Sum();*/
-
- RowVector PercEV(AdjEV);
- PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
-
- AdjEV << SP(AdjEV,pow(CircleLaw.Columns(1,AdjEV.Ncols()),-1));
-
- SortDescending(AdjEV);
- int maxEV = 1;
- float threshold = 0.98;
- for(int ctr_i = 1; ctr_i < PercEV.Ncols(); ctr_i++){
- if((PercEV(ctr_i)<threshold)&&(PercEV(ctr_i+1)>=threshold)){maxEV=ctr_i;}
- }
-
- if(maxEV<3){maxEV=PercEV.Ncols()/2;}
- RowVector NewEV;
- Matrix temp1;
- temp1 = abs(AdjEV);
- NewEV << temp1.SubMatrix(1,1,1,maxEV);
-
- AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
-
- out1 = AdjEV;
- out2 = PercEV;
- out3 = NewEV;
- out4 = maxEV;
- } //adj_eigspec
-
- void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2)
- {
- RowVector AdjEV, PercEV;
- AdjEV = in.Reverse();
- SortDescending(AdjEV);
-
- PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
- AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
- out1 = AdjEV;
- out2 = PercEV;
- } //adj_eigspec
-
- RowVector Feta(int n1, int n2)
- {
- float nu = (float) n1/n2;
- float bm = pow((1-sqrt(nu)),2.0);
- float bp = pow((1+sqrt(nu)),2.0);
-
- float lrange = 0.9*bm;
- float urange = 1.1*bp;
-
- RowVector eta(30*n1);
- float rangestepsize = (urange - lrange) / eta.Ncols();
- for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
- eta(ctr_i) = lrange + rangestepsize * (ctr_i);
- }
-
- RowVector teta(10*n1);
- teta = 0;
- float stepsize = (bp - bm) / teta.Ncols();
- for(int ctr_i = 1; ctr_i <= teta.Ncols(); ctr_i++){
- teta(ctr_i) = stepsize*(ctr_i);
- }
-
- RowVector feta(teta);
- feta = SP(pow(2*M_PI*nu*(teta + bm),-1), pow(SP(teta, bp-bm-teta),0.5));
-
- teta = teta + bm;
-
- RowVector claw(eta);
- claw = 0;
- for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
- double tmpval = 0.0;
- for(int ctr_j = 1; ctr_j <= teta.Ncols(); ctr_j++){
- if(( double(teta(ctr_j))/double(eta(ctr_i)) )<1)
- tmpval += feta(ctr_j);
- }
- claw(ctr_i) = n1*(1-stepsize*tmpval);
- }
-
- RowVector Res(n1); //invert the CDF
- Res = 0;
- for(int ctr_i = 1; ctr_i < eta.Ncols(); ctr_i++){ //Should this be <= instead of <?
- if(floor(claw(ctr_i))>floor(claw(ctr_i+1))){
- Res(int(floor(claw(ctr_i)))) = eta(ctr_i);
- }
- }
-
- return Res;
- } //RowVector Feta
-
- RowVector cumsum(const RowVector& Inp)
- {
- UpperTriangularMatrix UT(Inp.Ncols());
- UT=1.0;
- RowVector Res;
- Res = Inp * UT;
- return Res;
- } //RowVector cumsum
-
- int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, SymmetricMatrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which)
- {
- std_pca(in,weights,Corr,tmpE,tmpD);
-
- int maxEV = 1;
- RowVector NewEV;
- adj_eigspec(tmpD.AsRow(),AdjEV,PercEV,NewEV,maxEV,in.Ncols(),resels);
-
- int res;
- PPCA = ppca_est(NewEV, in.Ncols(),resels);
- ColumnVector tmp = ppca_select(PPCA, res, maxEV, which);
-
- PPCA = tmp | PPCA;
- return res;
- } //int ppca_dim
-
- int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which)
- {
- RowVector tmpD;
- Matrix tmpE;
- SymmetricMatrix Corr;
-
- int res = ppca_dim(in, weights, PPCA, AdjEV, PercEV, Corr, tmpE, tmpD, resels, which);
- return res;
- } //int ppca_dim
-
- int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which)
- {
- ColumnVector PPCA;
- RowVector AdjEV, PercEV;
- int res = ppca_dim(in,weights,PPCA,AdjEV,PercEV,resels,which);
- return res;
- } //int ppca_dim
-
- ColumnVector ppca_select(Matrix& PPCAest, int& dim, int maxEV, string which)
- {
- RowVector estimators(5);
- estimators = 1.0;
-
- for(int ctr=1; ctr<=PPCAest.Ncols(); ctr++){
- PPCAest.Column(ctr) = (PPCAest.Column(ctr) -
- min(PPCAest.Column(ctr)).AsScalar()) /
- ( max(PPCAest.Column(ctr)).AsScalar() -
- min(PPCAest.Column(ctr)).AsScalar());
- }
-
- int ctr_i = 1;
- while((ctr_i< PPCAest.Nrows()-1)&&
- (PPCAest(ctr_i,2) < PPCAest(ctr_i+1,2))&&(ctr_i<maxEV))
- {estimators(1)=ctr_i+1;ctr_i++;}
- ctr_i = 1;
- while((ctr_i< PPCAest.Nrows()-1)&&
- (PPCAest(ctr_i,3) < PPCAest(ctr_i+1,3))&&(ctr_i<maxEV))
- {estimators(2)=ctr_i+1;ctr_i++;}
- ctr_i = 1;
- while((ctr_i< PPCAest.Nrows()-1)&&
- (PPCAest(ctr_i,4) < PPCAest(ctr_i+1,4))&&(ctr_i<maxEV))
- {estimators(3)=ctr_i+1;ctr_i++;}
- ctr_i = 1;
- while((ctr_i< PPCAest.Nrows()-1)&&
- (PPCAest(ctr_i,5) < PPCAest(ctr_i+1,5))&&(ctr_i<maxEV))
- {estimators(4)=ctr_i+1;ctr_i++;}
- ctr_i = 1;
- while((ctr_i< PPCAest.Nrows()-1)&&
- (PPCAest(ctr_i,6) < PPCAest(ctr_i+1,6))&&(ctr_i<maxEV))
- {estimators(5)=ctr_i+1;ctr_i++;}
-
- int res = 0;
- ColumnVector PPCA;
- RowVector PercEV(PPCAest.Column(1).t());
- PercEV = cumsum(PercEV / sum(PercEV,2).AsScalar());
-
- if(which == string("aut")) {
- if(int(estimators(2)) < int(estimators(1)) &&
- float(PercEV(int(estimators(2))))>0.8){
- res=int(estimators(2));
- PPCA << PPCAest.Column(3);
- }else{
- res = int(estimators(1));
- PPCA << PPCAest.Column(2);
- }
- }
- if(which == string("lap")){
- res = int(estimators(1));
- PPCA << PPCAest.Column(2);
- }
- if(which == string("bic")){
- res = int(estimators(2));
- PPCA << PPCAest.Column(3);
- }
- if(which == string("mdl")){
- res = int(estimators(3));
- PPCA << PPCAest.Column(4);
- }
- if(which == string("rrn")){
- res = int(estimators(4));
- PPCA << PPCAest.Column(5);
- }
- if(which == string("aic")){
- res = int(estimators(5));
- PPCA << PPCAest.Column(6);
- }
- if(which == string("median")){
- RowVector unsorted(estimators);
- SortAscending(unsorted);
- ctr_i=1;
- res=int(unsorted(3));
- while(res != int(estimators(ctr_i)))
- ctr_i++;
- PPCA << PPCAest.Column(ctr_i);
- }
- if(res==0 || which == string("mean")){//median estimator
- PPCA = mean(PPCAest.Columns(2,6),2);
- res=int(mean(estimators,2).AsScalar());
- }
-
- dim = res;
- return PPCA;
- } //RowVector ppca_select
-
- Matrix ppca_est(const RowVector& eigenvalues, const int N1, const float N2)
- {
- Matrix Res;
- Res = ppca_est(eigenvalues, (int) floor(N1/(2.5*N2)));
- return Res;
- } //Matrix ppca_est
-
- Matrix ppca_est(const RowVector& eigenvalues, const int N)
- {
- RowVector logLambda(eigenvalues);
- logLambda = log(eigenvalues);
-
- int d = logLambda.Ncols();
-
- RowVector k(d);
- for(int ctr_i = 1; ctr_i <=d; ctr_i++){
- k(ctr_i)=ctr_i;
- }
-
- RowVector m(d);
- m=d*k-0.5*SP(k,k+1);
-
- RowVector loggam(d);
- loggam = 0.5*k.Reverse();
- for(int ctr_i = 1; ctr_i <=d; ctr_i++){
- loggam(ctr_i)=lgam(loggam(ctr_i));
- }
- loggam = cumsum(loggam);
-
- RowVector l_probU(d);
- l_probU = -log(2)*k + loggam - cumsum(0.5*log(M_PI)*k.Reverse());
-
- RowVector tmp1;
- tmp1 = -cumsum(eigenvalues).Reverse()+sum(eigenvalues,2).AsScalar();
- tmp1(1) = 0.95*tmp1(2);
- tmp1=tmp1.Reverse();
-
- RowVector tmp2;
- tmp2 = -cumsum(logLambda).Reverse()+sum(logLambda,2).AsScalar();
- tmp2(1)=tmp2(2);
- tmp2=tmp2.Reverse();
-
- RowVector tmp3;
- tmp3 = d-k;
- tmp3(d) = 1.0;
-
- RowVector tmp4;
- tmp4 = SP(tmp1,pow(tmp3,-1));
- for(int ctr_i = 1; ctr_i <=d; ctr_i++){
- if(tmp4(ctr_i)<0.01){tmp4(ctr_i)=0.01;}
- if(tmp3(ctr_i)<0.01){tmp3(ctr_i)=0.01;}
- if(tmp1(ctr_i)<0.01){tmp1(ctr_i)=0.01;}
- }
-
- RowVector l_nu;
- l_nu = SP(-N/2*(d-k),log(tmp4));
- l_nu(d) = 0;
-
- RowVector l_lam;
- l_lam = -(N/2)*cumsum(logLambda);
-
- RowVector l_lhood;
- l_lhood = SP(pow(tmp3,-1),tmp2)-log(SP(pow(tmp3,-1),tmp1));
-
- Matrix t1,t2, t3;
- UpperTriangularMatrix triu(d);
- triu = 1.0;
- for(int ctr_i = 1; ctr_i <= triu.Ncols(); ctr_i++){
- triu(ctr_i,ctr_i)=0.0;
- }
- t1 = (ones(d,1) * eigenvalues);
- t1 = SP(triu,t1.t() - t1);
- t2 = pow(tmp4,-1).t()*ones(1,d);
- t3 = ones(d,1)*pow(eigenvalues,-1);
- t2 = SP(triu, t2.t()-t3.t());
- for(int ctr_i = 1; ctr_i <= t1.Ncols(); ctr_i++){
- for(int ctr_j = 1; ctr_j <= t1.Nrows(); ctr_j++){
- if(t1(ctr_j,ctr_i)<=0){t1(ctr_j,ctr_i)=1;}
- }
- }
- for(int ctr_i = 1; ctr_i <= t2.Ncols(); ctr_i++){
- for(int ctr_j = 1; ctr_j <= t2.Nrows(); ctr_j++){
- if(t2(ctr_j,ctr_i)<=0){t2(ctr_j,ctr_i)=1;}
- }
- }
- t1 = cumsum(sum(log(t1),2).AsRow());
- t2 = cumsum(sum(log(t2),2).AsRow());
-
- RowVector l_Az(d);
- l_Az << (t1+t2);
-
- RowVector l_lap;
- l_lap = l_probU + l_nu +l_Az + l_lam + 0.5*log(2*M_PI)*(m+k)-0.5*log(N)*k;
-
- RowVector l_BIC;
- l_BIC = l_lam + l_nu - 0.5*log(N)*(m+k);
-
- RowVector l_RRN;
- l_RRN = -0.5*N*SP(k,log(SP(cumsum(eigenvalues),pow(k,-1))))+l_nu;
-
- RowVector l_AIC;
- l_AIC = -2*N*SP(tmp3,l_lhood)+ 2*(1+d*k+0.5*(k-1));
- l_AIC = -l_AIC;
-
- RowVector l_MDL;
- l_MDL = -N*SP(tmp3,l_lhood)+ 0.5*(1+d*k+0.5*(k-1))*log(N);
- l_MDL = -l_MDL;
-
- Matrix Res;
-
- Res = eigenvalues.t();
- Res |= l_lap.t();
- Res |= l_BIC.t();
- Res |= l_MDL.t();
- Res |= l_RRN.t();
- Res |= l_AIC.t();
-
-
- return Res;
- } //Matrix ppca_est
-
- ColumnVector acf(const ColumnVector& in, int order)
- {
- double meanval;
- meanval = mean(in).AsScalar();
- int tpoints = in.Nrows();
-
- ColumnVector y, res;
- Matrix X, tmp;
-
- y = in.Rows(order+1,tpoints) - meanval;
- X = zeros(order + 1, order);
- for(int ctr1 = 1; ctr1 <= order; ctr1++)
- X.Column(ctr1) = in.Rows(order + 1 - ctr1, tpoints - ctr1) - meanval;
- tmp = X.t()*X;
- tmp = tmp.i();
- tmp = tmp * X.t();
- res << tmp * y;
- return res;
- } //ColumnVector acf
-
- ColumnVector pacf(const ColumnVector& in, int maxorder)
- {
- int tpoint = in.Nrows();
- ColumnVector res;
- res = acf(in, maxorder);
- for(int ctr1 = 1; ctr1 <= maxorder; ctr1++)
- if ( res.Column(ctr1).AsScalar() < (1/tpoint) + 2/(float)std::pow(tpoint,0.5))
- res.Column(ctr1) = 0;
- return res;
- } //ColumnVector pacf
-
- Matrix est_ar(const Matrix& Mat, int maxorder)
- {
- Matrix res;
- res = zeros(maxorder, Mat.Ncols());
- ColumnVector tmp;
- for (int ctr = 1; ctr <= Mat.Ncols(); ctr++){
- tmp = pacf(Mat.Column(ctr));
- res.Column(ctr) = tmp;
- }
- return res;
- } //Matrix est_ar
-
- ColumnVector gen_ar(const ColumnVector& in, int maxorder)
- {
- float sdev;
- sdev = stdev(in).AsScalar();
- ColumnVector x, arcoeff, scaled;
- scaled = in / sdev;
- arcoeff = pacf( scaled, maxorder);
- x = normrnd(in.Nrows(),1).AsColumn() * sdev;
- for(int ctr1=2; ctr1 <= in.Nrows(); ctr1++)
- for(int ctr2 = 1; ctr2 <= maxorder; ctr2++)
- x(ctr1) = arcoeff(ctr2) * x(std::max(1,int(ctr1-ctr2))) + x(ctr1);
- return x;
- } //ColumnVector gen_ar
-
- Matrix gen_ar(const Matrix& in, int maxorder)
- {
- Matrix res;
- res = in;
- ColumnVector tmp;
- for(int ctr=1; ctr <= in.Ncols(); ctr++){
- tmp = in.Column(ctr);
- res.Column(ctr) = gen_ar(tmp, maxorder);
- }
- return res;
- } //Matrix gen_ar
-
- Matrix gen_arCorr(const Matrix& in, int maxorder)
- {
- Matrix res;
- res = zeros(in.Nrows(), in.Nrows());
- ColumnVector tmp;
- for(int ctr=1; ctr<= in.Ncols(); ctr++){
- tmp = in.Column(ctr);
- tmp = gen_ar(tmp, maxorder);
- res += tmp * tmp.t();
- }
- return res;
- } //Matrix gen_arCorr
-
- void basicGLM::olsfit(const Matrix& data, const Matrix& design,
- const Matrix& contrasts, int requestedDOF)
- {
- beta = zeros(design.Ncols(),1);
- residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
- t = zeros(1); z = zeros(1); p=-1.0*ones(1);
- dof = (int)-1; cbeta = -1.0*ones(1);
-
- if(data.Nrows()==design.Nrows()){
- Matrix dat = data;
- Matrix tmp = design.t()*design;
- Matrix pinvdes = tmp.i()*design.t();
-
- beta = pinvdes * dat;
- residu = dat - design*beta;
- dof = ols_dof(design);
- if ( requestedDOF>0)
- dof = requestedDOF;
- sigsq = sum(SP(residu,residu))/dof;
-
- float fact = float(dof) / design.Ncols();
- f_fmf = SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact;
-
- pf_fmf = f_fmf.Row(1);
- for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
- pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),dof,f_fmf.Column(ctr1).AsScalar());
-
- if(contrasts.Storage()>0 && contrasts.Ncols()==beta.Nrows()){
- cbeta = contrasts*beta;
- Matrix tmp = contrasts*pinvdes*pinvdes.t()*contrasts.t();
- varcb = diag(tmp)*sigsq;
- t = SP(cbeta,pow(varcb,-0.5));
- z = t; p=t;
- for(int ctr1=1;ctr1<=t.Ncols();ctr1++){
- ColumnVector tmp = t.Column(ctr1);
- T2z::ComputeZStats(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
- z.Column(ctr1) << tmp;
- T2z::ComputePs(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
- p.Column(ctr1) << exp(tmp);
- }
- }
- }
-
- }
-
-
-}
diff --git a/melhlprfns.h b/melhlprfns.h
deleted file mode 100644
index 1ff66df9b639d2c80db736b777c8309f087d8e06..0000000000000000000000000000000000000000
--- a/melhlprfns.h
+++ /dev/null
@@ -1,133 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melhlprfns.cc - misc functions
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#ifndef __MELODICHLPR_h
-#define __MELODICHLPR_h
-
-#include "newimage/newimageall.h"
-#include "newmatap.h"
-#include "newmatio.h"
-
- #ifdef __APPLE__
- #include <mach/mach.h>
- #define mmsg(msg) { \
- struct task_basic_info t_info; \
- mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
- if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
- { \
- cout << " MEM: " << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
- } \
- }
- #else
- #define mmsg(msg) { \
- cout << msg; \
- }
- #endif
-
-using namespace NEWIMAGE;
-
-namespace Melodic{
-
- void update_mask(volume<float>& mask, Matrix& Data);
- void del_vols(volume4D<float>& in, int howmany);
-
- Matrix smoothColumns(const Matrix& inp);
- Matrix calc_FFT(const Matrix& Mat, const bool logpwr = 0);
-
- Matrix convert_to_pbsc(Matrix& Mat);
-
- RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6, int econ = 20000);
- void varnorm(Matrix& in, const RowVector& vars);
- RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim = 30, float level = 1.6, int econ = 20000);
-
- Matrix SP2(const Matrix& in, const Matrix& weights, int econ = 20000);
- void SP3(Matrix& in, const Matrix& weights);
- void SP4(Matrix& in, const Matrix& weights);
-
- RowVector Feta(int n1,int n2);
- RowVector cumsum(const RowVector& Inp);
-
- Matrix corrcoef(const Matrix& in1, const Matrix& in2);
- Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part);
- float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
- float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite);
- void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
- void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite);
- void calc_white(const SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite);
-
- void std_pca(const Matrix& Mat, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
- void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
- void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
- void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
-
- float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim = 1);
- RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows);
- RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows);
- Matrix krprod(const Matrix& cols, const Matrix& rows);
- Matrix krapprox(const Matrix& Mat, int size_col, int dim = 1);
-
- void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2, RowVector& out3, int& out4, int num_vox, float resels);
- void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2);
-
- int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, SymmetricMatrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which);
- int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which);
- int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which);
- ColumnVector ppca_select(Matrix& PPCAest, int& dim, int maxEV, string which);
- Matrix ppca_est(const RowVector& eigenvalues, const int N1, const float N2);
- Matrix ppca_est(const RowVector& eigenvalues, const int N);
-
- ColumnVector acf(const ColumnVector& in, int order);
- ColumnVector pacf(const ColumnVector& in, int maxorder = 1);
- Matrix est_ar(const Matrix& Mat, int maxorder);
- ColumnVector gen_ar(const ColumnVector& in, int maxorder = 1);
- Matrix gen_ar(const Matrix& in, int maxorder);
- Matrix gen_arCorr(const Matrix& in, int maxorder);
-
- class basicGLM{
- public:
-
- //constructor
- basicGLM(){}
-
- //destructor
- ~basicGLM(){}
-
- void olsfit(const Matrix& data, const Matrix& design,
- const Matrix& contrasts, int DOFadjust = -1);
-
- inline Matrix& get_t(){return t;}
- inline Matrix& get_z(){return z;}
- inline Matrix& get_p(){return p;}
- inline Matrix& get_f_fmf(){return f_fmf;}
- inline Matrix& get_pf_fmf(){return pf_fmf;}
- inline Matrix& get_cbeta(){return cbeta;}
- inline Matrix& get_beta(){return beta;}
- inline Matrix& get_varcb(){return varcb;}
- inline Matrix& get_sigsq(){return sigsq;}
- inline Matrix& get_residu(){return residu;}
- inline int get_dof(){return dof;}
-
- private:
- Matrix beta;
- Matrix residu;
- Matrix sigsq;
- Matrix varcb;
- Matrix cbeta;
- Matrix f_fmf, pf_fmf;
- int dof;
- Matrix t;
- Matrix z;
- Matrix p;
- };
-// Matrix glm_ols(const Matrix& dat, const Matrix& design);
-}
-
-#endif
diff --git a/melica.cc b/melica.cc
deleted file mode 100644
index 3d6ed123f22a93bcd9a05d09eab0131f3923f79d..0000000000000000000000000000000000000000
--- a/melica.cc
+++ /dev/null
@@ -1,484 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melica.cc - ICA estimation
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include <stdlib.h>
-#include "newimage/newimageall.h"
-#include "utils/log.h"
-#include "meloptions.h"
-#include "meldata.h"
-#include "melodic.h"
-#include "newmatap.h"
-#include "newmatio.h"
-#include "melica.h"
-#include "melpca.h"
-#include "melhlprfns.h"
-#include "miscmaths/miscprob.h"
-
-using namespace Utilities;
-using namespace NEWIMAGE;
-
-namespace Melodic {
-
- void MelodicICA::ica_fastica_symm(const Matrix &Data){
- // based on Aapo Hyvärinen's fastica method
- // see www.cis.hut.fi/projects/ica/fastica/
-
- //initialize matrices
- Matrix redUMM_old, rank1_old;
- Matrix tmpU;
- //srand((unsigned int)timer(NULL));
- redUMM = melodat.get_white()*
- unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
-
- if(opts.debug.value())
- cerr << "redUMM init submatrix : " << endl << redUMM.SubMatrix(1,2,1,2) << endl;
-
- if(opts.guessfname.value().size()>1){
- message(" Use columns in " << opts.guessfname.value()
- << " as initial values for the mixing matrix " <<endl);
- Matrix guess ;
- guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
- redUMM.Columns(1,guess.Ncols()) = guess;
- }
-
- symm_orth(redUMM);
-
- int itt_ctr,itt_ctr2=1,cum_itt=0,newmaxitts = opts.maxNumItt.value();
- double minAbsSin = 1.0, minAbsSin2 = 1.0;
- if(opts.approach.value() == string("tica"))
- opts.maxNumItt.set_T(opts.rank1interval.value());
-
- rank1_old = melodat.get_dewhite()*redUMM;
- rank1_old = melodat.expand_dimred(rank1_old);
- rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
- do{// TICA loop
- itt_ctr = 1;
- do{ // da loop!!!
- redUMM_old = redUMM;
- //calculate IC estimates
- tmpU = Data.t() * redUMM;
-
- //update redUMM depending on nonlinearity
- if(opts.nonlinearity.value()=="pow4"){
- redUMM = (Data * pow(tmpU,3.0)) / samples - 3 * redUMM;
- }
- if(opts.nonlinearity.value()=="pow3"){
- tmpU /= opts.nlconst1.value();
- redUMM = 3 * (Data * pow(tmpU,2.0)) / samples -
- (SP(ones(dim,1)*sum(tmpU,1),redUMM))/ samples;
- }
- if(opts.nonlinearity.value()=="tanh"){
- Matrix hyptanh;
- hyptanh = tanh(opts.nlconst1.value()*tmpU);
- redUMM = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
- sum(1-pow(hyptanh,2),1),redUMM))/samples;
- }
- if(opts.nonlinearity.value()=="gauss"){
- Matrix tmpUsq;
- Matrix tmpU2;
- Matrix gauss;
- Matrix dgauss;
- tmpUsq = pow(tmpU,2);
- tmpU2 = exp(-(opts.nlconst2.value()/2) * tmpUsq);
- gauss = SP(tmpU,tmpU2);
- dgauss = SP(1-opts.nlconst2.value()*tmpUsq,tmpU2);
- redUMM = (Data * gauss - SP(ones(dim,1)*
- sum(dgauss,1),redUMM))/samples;
- }
-
- // orthogonalize the unmixing-matrix
- symm_orth(redUMM);
- if(opts.approach.value() == string("tica")){
- message("");
- }
-
- //termination condition : angle between old and new < epsilon
- minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum();
-
- message(" Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);
- // if((abs(minAbsSin) < opts.epsilon.value())&&
- // (opts.approach.value()!=string("tica"))){ break;}
- if((abs(minAbsSin) < opts.epsilon.value())){ break;}
-
- itt_ctr++;
- } while(itt_ctr < opts.maxNumItt.value());
- cum_itt += itt_ctr;
- itt_ctr2++;
- if(opts.approach.value() == string("tica")){
- message(" Rank-1 approximation of the time courses; ");
- Matrix temp(melodat.get_dewhite() * redUMM);
- temp = melodat.expand_dimred(temp);
- temp = krapprox(temp,int(temp.Nrows()/melodat.get_numfiles()));
- minAbsSin2 = 1 - diag(abs(corrcoef(temp,rank1_old))).Minimum();
- rank1_old = temp;
- temp = melodat.reduce_dimred(temp);
- redUMM = melodat.get_white() * temp;
-
- message(" change : " << minAbsSin2 << endl);
- if(abs(minAbsSin2) < opts.epsilonS.value() && abs(minAbsSin) < opts.epsilon.value()){ break;}
- }
- } while(
- (itt_ctr2 < newmaxitts/opts.maxNumItt.value()) &&
- (opts.approach.value() == string("tica")) &&
- cum_itt < newmaxitts);
-
- if((itt_ctr>=opts.maxNumItt.value() && (opts.approach.value()!=string("tica")))
- || (cum_itt >= newmaxitts && opts.approach.value()==string("tica"))){
- cerr << " No convergence after " << cum_itt <<" steps "<<endl;
- } else {
- message(" Convergence after " << cum_itt <<" steps " << endl << endl);
- no_convergence = false;
-
- {Matrix temp(melodat.get_dewhite() * redUMM);
- melodat.set_mix(temp);}
- {Matrix temp(redUMM.t()*melodat.get_white());
- melodat.set_unmix(temp);}
- }
- }
-
- void MelodicICA::ica_fastica_defl(const Matrix &Data){
- if(!opts.explicitnums || opts.numICs.value()>dim){
- opts.numICs.set_T(dim);
- message(" Using numICs:" << opts.numICs.value() << endl);
- }
-
- //redUMM = zeros(dim); // got to start somewhere
- redUMM = melodat.get_white()*
- unifrnd(melodat.get_white().Ncols(),opts.numICs.value());
- redUMM = zeros(melodat.get_white().Nrows(),opts.numICs.value());
- Matrix guess;
- int guesses=0;
- if(opts.guessfname.value().size()>1){
- message(" Use columns in " << opts.guessfname.value() << " as initial values for the mixing matrix " <<endl);
- guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
- guesses = guess.Ncols();
- }
-
- int ctrIC = 1;
- int numRestart = 0;
- while(ctrIC<=opts.numICs.value()){
- message(" Extracting IC " << ctrIC << " ... ");
- ColumnVector w;
- ColumnVector w_old;
- ColumnVector tmpU;
- if(ctrIC <= guesses){
- w = w - redUMM * redUMM.t() * w;
- w = w / norm2(w);
- w_old = zeros(w.Nrows(),1);
- int itt_ctr = 1;
- do{
- w_old = w;
- tmpU = Data.t() * w;
- if(opts.nonlinearity.value()=="pow4"){
- w = (Data * pow(tmpU,3.0)) / samples - 3 * w;
- }
- if(opts.nonlinearity.value()=="tanh"){
- ColumnVector hyptanh;
- hyptanh = tanh(opts.nlconst1.value()*tmpU);
- w = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
- sum(1-pow(hyptanh,2),1),w))/samples;
- }
- if(opts.nonlinearity.value()=="pow3"){
- tmpU /= opts.nlconst1.value();
- w = 3*(Data * pow(tmpU,2.0)) / samples - 2*(SP(ones(dim,1)*
- sum(tmpU,1),w))/samples;
- }
- if(opts.nonlinearity.value()=="gauss"){
- ColumnVector tmpUsq;
- ColumnVector tmpU2;
- ColumnVector gauss;
- ColumnVector dgauss;
- tmpUsq = pow(tmpU,2);
- tmpU2 = exp(-(opts.nlconst2.value()/2) * tmpUsq);
- gauss = SP(tmpU,tmpU2);
- dgauss = SP(1-opts.nlconst2.value()*tmpUsq,tmpU2);
- w = (Data * gauss - SP(ones(dim,1)*
- sum(dgauss,1),w))/samples;
- }
-
- // orthogonalize w
- w = w - redUMM * redUMM.t() * w;
- w = w / norm2(w);
-
- //termination condition : angle between old and new < epsilon
- if(((norm2(w-w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10)) ||
- ((norm2(w+w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10))){
- break;
- }
- //cout << norm2(w-w_old) << " " << norm2(w+w_old) << endl;
- itt_ctr++;
- } while(itt_ctr <= opts.maxNumItt.value());
-
- if(itt_ctr<opts.maxNumItt.value()){
- redUMM.Column(ctrIC) = w;
- message(" estimated using " << itt_ctr << " iterations " << endl);
- ctrIC++;
- numRestart = 0;
- } else{
- if(numRestart > opts.maxRestart.value()){
- message(endl << " Estimation failed after "
- << numRestart << " attempts "
- << " giving up " << endl);
- break;
- }else{
- numRestart++;
- message(endl <<" Estimation failed - restart "
- << numRestart << endl);
- }
- }
- }
- if(numRestart < opts.maxRestart.value()){
- no_convergence = false;
- {Matrix temp(melodat.get_dewhite() * redUMM);
- melodat.set_mix(temp);}
- {Matrix temp(redUMM.t()*melodat.get_white());
- melodat.set_unmix(temp);}
- }
- }
- }
-
- void MelodicICA::ica_maxent(const Matrix &Data){
- // based on Aapo Hyvärinen's fastica method
- // see www.cis.hut.fi/projects/ica/fastica/
- message(" MAXENT " << endl);
- //initialize matrices
- Matrix redUMM_old;
- Matrix tmpU;
- Matrix gtmpU;
- double lambda = 0.015/std::log((double)melodat.get_white().Ncols());
-
- //srand((unsigned int)timer(NULL));
- redUMM = melodat.get_white()*
- unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
-
- if(opts.guessfname.value().size()>1){
- message(" Use columns in " << opts.guessfname.value()
- << " as initial values for the mixing matrix " <<endl);
- Matrix guess ;
- guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
- redUMM.Columns(1,guess.Ncols()) = guess;
- }
-
- // symm_orth(redUMM);
- int itt_ctr=1;
- double minAbsSin = 1.0;
- Matrix Id;
- Id = IdentityMatrix(redUMM.Ncols());
- //cerr << " nonlinearity : " << opts.nonlinearity.value() << endl;
-
- do{ // da loop!!!
- redUMM_old = redUMM;
- //calculate IC estimates
- tmpU = Data.t() * redUMM;
- if(opts.nonlinearity.value()=="tanh"){
- //Matrix hyptanh;
- //hyptanh = tanh(opts.nlconst1.value()*tmpU);
- //redUMM = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
- //sum(1-pow(hyptanh,2),1),redUMM))/samples;
- gtmpU = tanh(opts.nlconst1.value()*tmpU);
- redUMM = redUMM + lambda*(Id+(1-2*gtmpU.t()*tmpU))*redUMM;
- }
- if(opts.nonlinearity.value()=="gauss"){
- gtmpU = pow(1+exp(-(opts.nlconst2.value()/2) * tmpU),-1);
- redUMM = redUMM + lambda*(Id - (gtmpU.t()-tmpU.t())*tmpU)*redUMM;
- }
-
- //termination condition : angle between old and new < epsilon
- minAbsSin = abs(1 - diag(abs(redUMM.t()*redUMM_old)).Minimum());
- message(" Step no. " << itt_ctr << " change : " << minAbsSin << endl);
- if(abs(minAbsSin) < opts.epsilon.value()){ break;}
-
- itt_ctr++;
- } while(itt_ctr < opts.maxNumItt.value());
-
- if(itt_ctr>=opts.maxNumItt.value()){
- cerr << " No convergence after " << itt_ctr <<" steps "<<endl;
- } else {
- message(" Convergence after " << itt_ctr <<" steps " << endl << endl);
- no_convergence = false;
- {Matrix temp(melodat.get_dewhite() * redUMM);
- melodat.set_mix(temp);}
- {Matrix temp(redUMM.t()*melodat.get_white());
- melodat.set_unmix(temp);}
- }
- }
-
- void MelodicICA::ica_jade(const Matrix &Data){
- int dim_sym = (int) (dim*(dim+1))/2;
- int num_CM = dim_sym;
- Matrix CM;
- Matrix R; R = IdentityMatrix(dim);
- Matrix Qij; Qij = zeros(dim);
- Matrix Xim;
- Matrix Xjm;
- Matrix scale; scale = ones(dim,1)/samples;
-
- for (int im =1; im <= dim; im++){
- Xim = Data.Row(im);
- write_ascii_matrix("Xim",Data.Row(1));
- //Qij = SP((scale * pow(Xim,2)),Data) * Data.t();//- R - 2*R.Column(im)*R.Column(im).t();
- Qij = (pow(Xim,2)) * Data.t();//- R - 2*R.Column(im)*R.Column(im).t();
- if(im==1){CM = Qij; write_ascii_matrix("CM",CM);exit(2);}else{CM |= Qij;}
- for (int jm = 1; jm < im; jm++){
- Xjm = Data.Row(jm);
- Qij = (SP((scale * SP(Xim,Xjm)),Data) * Data.t() -
- R.Column(im)*R.Column(jm).t() - R.Column(jm)*R.Column(im).t());
- Qij *= sqrt(2);
- CM |= Qij;
- }
- }
-
- write_ascii_matrix("CM",CM);
- Matrix redUMM; redUMM = IdentityMatrix(dim);
-
- bool exitloop = false;
- int ctr_itt = 0;
- int ctr_updates = 0;
- Matrix Givens; Givens = zeros(2,num_CM);
- Matrix Givens_ip; Givens_ip = zeros(2);
- Matrix Givens_ro; Givens_ro = zeros(2);
- double det_on, det_off;
- double cos_theta, sin_theta, theta;
-
- while(!exitloop && ctr_itt <= opts.maxNumItt.value()){
- ctr_itt++;
- cout << "loop" <<endl;
- for(int ctr_p = 1; ctr_p < dim; ctr_p++){
- for(int ctr_q = ctr_p+1; ctr_q <= dim; ctr_q++){
-
- for(int ctr_i = 0; ctr_i < num_CM; ctr_i++){
- int Ip = ctr_p + ctr_i * dim;
- int Iq = ctr_q + ctr_i * dim;
- Givens(1,ctr_i + 1) = CM(ctr_p,Ip) - CM(ctr_q,Iq);
- Givens(2,ctr_i + 1) = CM(ctr_p,Iq) - CM(ctr_q,Ip);
- }
-
- Givens_ip = Givens * Givens.t();
- det_on = Givens_ip(1,1) - Givens_ip(2,2);
- det_off = Givens_ip(1,2) + Givens_ip(2,1);
- theta = 0.5 * atan2(det_off, det_on + sqrt(det_on*det_on + det_off*det_off));
-
- cout << theta << endl;
-
- if(abs(theta) > opts.epsilon.value()){
- ctr_updates++;
- message(" Step No. "<< ctr_itt << " change : " << theta << endl);
-
- //create Givens rotation matrix
- cos_theta = cos(theta);
- sin_theta = sin(theta);
- Givens_ro(1,1) = cos_theta;
- Givens_ro(1,2) = -sin_theta;
- Givens_ro(2,1) = sin_theta;
- Givens_ro(2,2) = cos_theta;
-
- //update 2 entries of redUMM
- {Matrix tmp_redUMM;
- tmp_redUMM = redUMM.Column(ctr_p);
- tmp_redUMM |= redUMM.Column(ctr_q);
- tmp_redUMM *= Givens_ro;
- redUMM.Column(ctr_p) = tmp_redUMM.Column(1);
- redUMM.Column(ctr_q) = tmp_redUMM.Column(2);}
-
- //update Cumulant matrix
- {Matrix tmp_CM;
- tmp_CM = CM.Row(ctr_p);
- tmp_CM &= CM.Row(ctr_q);
- tmp_CM = Givens_ro.t() * tmp_CM;
- CM.Row(ctr_p) = tmp_CM.Row(1);
- CM.Row(ctr_q) = tmp_CM.Row(2);}
-
- //update Cumulant matrices
- for(int ctr_i = 0; ctr_i < num_CM; ctr_i++){
- int Ip = ctr_p + ctr_i * dim;
- int Iq = ctr_q + ctr_i * dim;
- CM.Column(Ip) = cos_theta*CM.Column(Ip)+sin_theta*CM.Column(Iq);
- CM.Column(Iq) = cos_theta*CM.Column(Iq)-sin_theta*CM.Column(Ip);
- }
- }else{
- exitloop = true;
- }
- }
- }
- }//while loop
- if(ctr_itt > opts.maxNumItt.value()){
- cerr << " No convergence after " << ctr_itt <<" steps "<<endl;
- } else {
- message(" Convergence after " << ctr_itt <<" steps " << endl << endl);
- no_convergence = false;
- {Matrix temp(melodat.get_dewhite() * redUMM);
- melodat.set_mix(temp);}
- {Matrix temp(redUMM.t()*melodat.get_white());
- melodat.set_unmix(temp);}
- }
- }
-
- Matrix MelodicICA::sign(const Matrix &Inp){
- Matrix Res = Inp;
- Res = 1;
- for(int ctr_i = 1; ctr_i <= Inp.Ncols(); ctr_i++){
- for(int ctr_j = 1; ctr_j <= Inp.Nrows(); ctr_j++){
- if(Inp(ctr_j,ctr_i)<0){Res(ctr_j,ctr_i)=-1;}
- }
- }
- return Res;
- }
-
- void MelodicICA::perf_ica(const Matrix &Data){
- message("Starting ICA estimation using " << opts.approach.value()
- << endl << endl);
- dim = Data.Nrows();
- samples = Data.Ncols();
- no_convergence = true;
- //switch to the chosen method
- if(opts.approach.value()==string("symm") ||
- opts.approach.value()==string("tica") ||
- opts.approach.value()==string("parafac") ||
- opts.approach.value()==string("concat"))
- ica_fastica_symm(Data);
- if(opts.approach.value()==string("defl"))
- ica_fastica_defl(Data);
- if(opts.approach.value()==string("jade"))
- ica_jade(Data);
- if(opts.approach.value()==string("maxent"))
- ica_maxent(Data);
-
- if(!no_convergence){//calculate the IC
- Matrix temp(melodat.get_unmix()*melodat.get_Data());
- // Add the mean time course again
- // temp += melodat.get_unmix()*melodat.get_meanC()*ones(1,temp.Ncols());
-
- //re-normalise the decomposition to std(mix)=1
- Matrix scales;
- scales = stdev(melodat.get_mix());
-
- //cerr << " SCALES 1 " << scales << endl;
- Matrix tmp, tmp2;
- tmp = SP(melodat.get_mix(),ones(melodat.get_mix().Nrows(),1)*pow(scales,-1));
- temp = SP(temp,scales.t()*ones(1,temp.Ncols()));
-
- scales = scales.t();
-
- melodat.set_mix(tmp);
-
- melodat.set_IC(temp);
-
- melodat.set_ICstats(scales);
- melodat.sort();
-
- //message("Calculating T- and S-modes " << endl);
- melodat.set_TSmode();
-
- }
- }
-}
-
-
diff --git a/melica.h b/melica.h
deleted file mode 100644
index dc9594033322e511a137af4d90e28eaedd890f09..0000000000000000000000000000000000000000
--- a/melica.h
+++ /dev/null
@@ -1,66 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melica.h - ICA estimation
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#ifndef __MELODICICA_h
-#define __MELODICICA_h
-#include "newimage/newimageall.h"
-#include "utils/log.h"
-#include "melpca.h"
-#include "meloptions.h"
-#include "meldata.h"
-#include "melodic.h"
-#include "newmatap.h"
-#include "newmatio.h"
-#include "melreport.h"
-#include "ggmix.h"
-
-using namespace Utilities;
-using namespace NEWIMAGE;
-
-namespace Melodic{
-
- class MelodicICA{
- public:
- MelodicICA(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger, MelodicReport &preport):
- melodat(pmelodat),
- opts(popts),
- logger(plogger),
- report(preport){}
-
- void perf_ica(const Matrix &Data);
-
- bool no_convergence;
-
- private:
- MelodicData &melodat;
- MelodicOptions &opts;
- Log &logger;
- MelodicReport &report;
-
- int dim;
- int samples;
- Matrix redUMM;
-
- void ica_fastica_symm(const Matrix &Data);
- void ica_fastica_defl(const Matrix &Data);
- void ica_maxent(const Matrix &Data);
- void ica_jade(const Matrix &Data);
- //void tica();
- //void tica_concat();
- //void parafac();
- Matrix randm(const int dim1, const int dim2);
- void sort();
- Matrix sign(const Matrix &Inp);
-
- };
-}
-
-#endif
diff --git a/melodic.h b/melodic.h
deleted file mode 100644
index 3e3e32d9f8f2a9fab19e4a7816f33757260750e6..0000000000000000000000000000000000000000
--- a/melodic.h
+++ /dev/null
@@ -1,88 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melodic.h - main program header
-
- Christian F. Beckmann and Matthew Webster, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#ifndef __MELODIC_h
-#define __MELODIC_h
-
-#include <strstream>
-
-#ifdef __APPLE__
-#include <mach/mach.h>
-#define memmsg(msg) { \
- MelodicOptions&opt = MelodicOptions::getInstance(); \
- struct task_basic_info t_info; \
- mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
- if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
- { \
- if(opt.debug.value()) {\
- cout << " MEM: " << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
- cout.flush(); \
- } \
- } \
-}
-#else
-#define memmsg(msg) { \
- MelodicOptions&opt = MelodicOptions::getInstance(); \
- if(opt.debug.value()) {\
- cout << msg; \
- cout.flush(); \
- } \
-}
-#endif
-
-
-
-// a simple message macro that takes care of cout and log
-#define message(msg) { \
- MelodicOptions& opt = MelodicOptions::getInstance(); \
- if(opt.verbose.value()) \
- { \
- cout << msg; \
- } \
- Log& logger = LogSingleton::getInstance(); \
- logger.str() << msg; \
- cout.flush(); \
-}
-
-#define dbgmsg(msg) { \
- MelodicOptions&opt = MelodicOptions::getInstance(); \
- if(opt.debug.value()) {\
- cout << msg; } \
-}
-
-#define outMsize(msg,Mat) { \
- MelodicOptions& opt = MelodicOptions::getInstance(); \
- if(opt.debug.value()) \
- cerr << " " << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
-}
-
-namespace Melodic{
-
-const string version = "3.141";
-
-// The two strings below specify the title and example usage that is
-// printed out as the help or usage message
-const string title=string("MELODIC (Version ")+version+")"+
- string("\n Multivariate Exploratory Linear Optimised Decomposition into Independent Components\n")+
- string("\nAuthor: Christian F. Beckmann \n Copyright(c) 2001-2013 University of Oxford");
-
-const string usageexmpl=string(" melodic -i <filename> <options>")+
- string("\n \t \t to run melodic")+
- // string("\n melodic -i <filename> --mix=melodic_mix")+
- // string(" --filter=\"string of component numbers\"")+
- // string("\n \t \t to remove estimated ICs from input")+
- string("\n melodic -i <filename> --ICs=melodic_IC")+
- string(" --mix=melodic_mix <options>")+
- string("\n \t \t to run Mixture Model based inference on estimated ICs")+
- string("\n melodic --help ");
-}
-
-#endif
diff --git a/meloptions.cc b/meloptions.cc
index 0d96dbed5c1125a8208e543f3bfe1703c78b285e..72835c345eaa1ac34a7dde20e642b76de164f2cc 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -157,10 +157,6 @@ MelodicOptions* MelodicOptions::gopt = NULL;
}
if (insta_fn.value() > ""){
varnorm.set_T(false);
- cout << " inc_fn = " << insta_fn.value() << endl;
- cout << " inc_idx = " << insta_idx.value() << endl;
- cout << " inc_partial = " << insta_partial.value() << endl;
- cout << " inc_varnorm = " << insta_varnorm.value() << endl;
}
//in the case of indirect inputs, create the vector of input names here
diff --git a/meloptions.h b/meloptions.h
deleted file mode 100644
index 020884ff7e5e23fcfb264d37b188d321c6644ab9..0000000000000000000000000000000000000000
--- a/meloptions.h
+++ /dev/null
@@ -1,543 +0,0 @@
- /* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- meloptions.h - class for command line options
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#ifndef __MELODICOPTIONS_h
-#define __MELODICOPTIONS_h
-
-#include <string>
-#include <strstream>
-#include <iostream>
-#include <iomanip>
-#include <fstream>
-#include <stdlib.h>
-#include <stdio.h>
-#include "utils/options.h"
-#include "utils/log.h"
-#include "melodic.h"
-
-using namespace Utilities;
-
-namespace Melodic {
-
-class MelodicOptions {
- public:
- static MelodicOptions& getInstance();
- ~MelodicOptions() { delete gopt; }
-
- string version;
- string binpath;
- string logfname;
- bool filtermode;
- bool explicitnums;
-
- Option<string> logdir;
- Option< std::vector<string> > inputfname;
-
- Option<string> outputfname;
-
- Option<string> maskfname;
- Option<bool> use_mask;
- Option<bool> update_mask;
- Option<bool> perf_bet;
- Option<float> threshold;
-
- Option<int> pca_dim;
- Option<string> pca_est;
- Option<bool> joined_whiten;
- Option<bool> joined_vn;
- Option<bool> dr_pca;
- Option<bool> migp;
- Option<int> migpN;
- Option<bool> migp_shuffle;
- Option<int> migp_factor;
- Option<bool> dr;
- Option<bool> dr_out;
- Option<float> vn_level;
- Option<int> numICs;
- Option<string> approach;
- Option<string> nonlinearity;
-
- Option<bool> varnorm;
- Option<bool> varnorm2;
- Option<bool> pbsc;
- Option<bool> pspec;
- Option<string> segment;
- Option<bool> tsmooth;
- Option<float> epsilon;
- Option<float> epsilonS;
- Option<int> maxNumItt;
- Option<int> maxRestart;
- Option<int> rank1interval;
-
- Option<string> mmthresh;
- Option<bool> perf_mm;
- Option<string> ICsfname;
- Option<string> filtermix;
- Option<string> smodename;
- Option<string> filter;
-
- Option<bool> genreport;
- Option<string> guireport;
- Option<string> bgimage;
- Option<float> tr;
- Option<bool> logPower;
- Option<bool> addsigchng;
- Option<bool> allPPCA;
- Option<bool> varplots;
- Option<bool> varvals;
-
- Option<string> fn_Tdesign;
- Option<string> fn_Tcon;
- Option<string> fn_TconF;
- Option<string> fn_Sdesign;
- Option<string> fn_Scon;
- Option<string> fn_SconF;
-
- Option<bool> output_all;
- Option<bool> output_unmix;
- Option<bool> output_MMstats;
- Option<bool> output_pca;
- Option<bool> output_white;
- Option<bool> output_origIC;
- Option<bool> output_mean;
-
- Option<bool> verbose;
- Option<bool> vers;
- Option<bool> copyright;
- Option<bool> help;
- Option<bool> debug;
-
- Option<string> guessfname;
- Option<string> paradigmfname;
- Option<string> axials_str;
-
- Option<int> dummy;
- Option<int> repeats;
- Option<int> seed;
- Option<float> nlconst1;
- Option<float> nlconst2;
- Option<float> smooth_probmap;
-
- Option<string> insta_fn;
- Option<int> insta_idx;
- Option<bool> insta_partial;
- Option<bool> insta_varnorm;
-
- Option<bool> remove_meanvol;
- Option<bool> remove_meantc;
- Option<bool> remove_endslices;
- Option<bool> rescale_nht;
-
- Option<bool> guess_remderiv;
- Option<bool> temporal;
-
- Option<float> retryfactor;
- Option<int> econ;
-
- void parse_command_line(int argc, char** argv, Log& logger, const string &p_version);
-
- private:
- MelodicOptions();
- const MelodicOptions& operator=(MelodicOptions&);
- MelodicOptions(MelodicOptions&);
-
- OptionParser options;
-
- static MelodicOptions* gopt;
-
- void print_usage(int argc, char *argv[]);
- void print_version();
- void print_copyright();
- void status();
-
-};
-
- inline MelodicOptions& MelodicOptions::getInstance(){
- if(gopt == NULL)
- gopt = new MelodicOptions();
-
- return *gopt;
- }
-
- inline MelodicOptions::MelodicOptions() :
- logdir(string("-o,--outdir"), string("log.txt"),
- string("output directory name\n"),
- false, requires_argument),
- inputfname(string("-i,--in"), std::vector<string>(),
- string("input file names (either single file name or comma-separated list or text file)"),
- true, requires_argument),
- outputfname(string("-O,--out"), string("melodic"),
- string("output file name"),
- false, requires_argument,false),
- maskfname(string("-m,--mask"), string(""),
- string("file name of mask for thresholding"),
- false, requires_argument),
- use_mask(string("--nomask"), true,
- string("switch off masking"),
- false, no_argument),
- update_mask(string("--update_mask"), true,
- string("switch off mask updating"),
- false, no_argument),
- perf_bet(string("--nobet"), true,
- string("\tswitch off BET"),
- false, no_argument),
- threshold(string("--bgthreshold"), 0.01,
- string("brain / non-brain threshold (only if --nobet selected)\n"),
- false, requires_argument),
- pca_dim(string("-d,--dim"), 0,
- string("dimensionality reduction into #num dimensions (default: automatic estimation)"),
- false, requires_argument),
- pca_est(string("--dimest"), string("lap"),
- string("use specific dim. estimation technique: lap, bic, mdl, aic, mean (default: lap)"),
- false, requires_argument),
- joined_whiten(string("--sep_whiten"), false,
- string("switch on separate whitening"),
- false, no_argument, false),
- joined_vn(string("--sep_vn"), true,
- string("switch on separate variance nomalisation (as opposed to separate VN)"),
- false, no_argument),
- dr_pca(string("--mod_pca"), true,
- string("switch off modified PCA for concat ICA"),
- false, no_argument, false),
- migp(string("--migp"), false,
- string("switch on MIGP data reduction"),
- false, no_argument, false),
- migpN(string("--migpN"), 0,
- string("Number of internal Eigenmaps"),
- false, requires_argument, false),
- migp_shuffle(string("--migp_shuffle"), true,
- string("Randomise MIGP file order (default: TRUE)"),
- false, no_argument, false),
- migp_factor(string("--migp_factor"), 2,
- string("Internal Factor of mem-threshold relative to number of Eigenmaps (default: 2)"),
- false, requires_argument, false),
- dr(string("--dr"), false,
- string("Dual Regression (default: false)"),
- false, no_argument, false),
- dr_out(string("--dr_out"), false,
- string("Dual Regression output for MIGP/concat ICA"),
- false, no_argument, false),
- vn_level(string("--vn_level"), float(2.3),
- string("variance nomalisation threshold level (Z> value is ignored)"),
- false, requires_argument, false),
- numICs(string("-n,--numICs"), -1,
- string("numer of IC's to extract (for deflation approach)"),
- false, requires_argument),
- approach(string("-a,--approach"), string("symm"),
- string("approach for decomposition, 2D: defl, symm (default), 3D: tica, concat (default)"),
- false, requires_argument),
- nonlinearity(string("--nl"), string("pow3"),
- string("\tnonlinearity: gauss, tanh, pow3, pow4"),
- false, requires_argument),
- varnorm(string("--vn,--varnorm"), true,
- string("switch off variance normalisation"),
- false, no_argument),
- varnorm2(string("--vn2"), true,
- string("switch off 2nd level variance normalisation"),
- false, no_argument, false),
- pbsc(string("--pbsc"), false,
- string(" switch on conversion to percent BOLD signal change"),
- false, no_argument, false),
- pspec(string("--pspec"), false,
- string(" switch on conversion to powerspectra"),
- false, no_argument, false),
- segment(string("--covarweight"), string(""),
- string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
- false, requires_argument),
- tsmooth(string("--spca"), false,
- string("smooth the eigenvectors prior to IC decomposition"),
- false, no_argument, false),
- epsilon(string("--eps"), 0.00005,
- string("minimum error change"),
- false, requires_argument),
- epsilonS(string("--epsS"), 0.03,
- string("minimum error change for rank-1 approximation in TICA"),
- false, requires_argument),
- maxNumItt(string("--maxit"), 500,
- string("\tmaximum number of iterations before restart"),
- false, requires_argument),
- maxRestart(string("--maxrestart"), -1,
- string("maximum number of restarts\n"),
- false, requires_argument),
- rank1interval(string("--rank1interval"), 10,
- string("number of iterations between rank-1 approximation (TICA)\n"),
- false, requires_argument,false),
- mmthresh(string("--mmthresh"), string("0.5"),
- string("threshold for Mixture Model based inference"),
- false, requires_argument),
- perf_mm(string("--no_mm"), true,
- string("\tswitch off mixture modelling on IC maps\n "),
- false, no_argument),
- ICsfname(string("--ICs"), string(""),
- string("\tinput filename of the IC components file for mixture modelling"),
- false, requires_argument),
- filtermix(string("--mix"), string(""),
- string("\tinput filename of mixing matrix for mixture modelling / filtering"),
- false, requires_argument),
- smodename(string("--smode"), string(""),
- string("\tinput filename of matrix of session modes for report generation"),
- false, requires_argument),
- filter(string("-f,--filter"), string(""),
- string("list of component numbers to remove\n "),
- false, requires_argument),
- genreport(string("--report"), false,
- string("generate Melodic web report"),
- false, no_argument),
- guireport(string("--guireport"), string(""),
- string("modify report for GUI use"),
- false, requires_argument, false),
- bgimage(string("--bgimage"), string(""),
- string("specify background image for report (default: mean image)\n "),
- false, requires_argument),
- tr(string("--tr"), 0.0,
- string("\tTR in seconds"),
- false, requires_argument),
- logPower(string("--logPower"), false,
- string("calculate log of power for frequency spectrum\n"),
- false, no_argument),
- addsigchng(string("--sigchng"), false,
- string("add signal change estimates to report pages\n"),
- false, no_argument, false),
- allPPCA(string("--allPPCA"), false,
- string("add all PPCA plots\n"),
- false, no_argument, false),
- varplots(string("--varplots"), false,
- string("add std error envelopes to time course plots\n"),
- false, no_argument, false),
- varvals(string("--varvals"), false,
- string("add rank1 values after plots\n"),
- false, no_argument, false),
- fn_Tdesign(string("--Tdes"), string(""),
- string(" design matrix across time-domain"),
- false, requires_argument),
- fn_Tcon(string("--Tcon"), string(""),
- string(" t-contrast matrix across time-domain"),
- false, requires_argument),
- fn_TconF(string("--Tconf"), string(""),
- string(" F-contrast matrix across time-domain"),
- false, requires_argument, false),
- fn_Sdesign(string("--Sdes"), string(""),
- string(" design matrix across subject-domain"),
- false, requires_argument),
- fn_Scon(string("--Scon"), string(""),
- string(" t-contrast matrix across subject-domain"),
- false, requires_argument),
- fn_SconF(string("--Sconf"), string(""),
- string(" F-contrast matrix across subject-domain"),
- false, requires_argument,false),
- output_all(string("--Oall"), false,
- string(" output everything"),
- false, no_argument),
- output_unmix(string("--Ounmix"), false,
- string("output unmixing matrix"),
- false, no_argument),
- output_MMstats(string("--Ostats"), false,
- string("output thresholded maps and probability maps"),
- false, no_argument),
- output_pca(string("--Opca"), false,
- string("\toutput PCA results"),
- false, no_argument),
- output_white(string("--Owhite"), false,
- string("output whitening/dewhitening matrices"),
- false, no_argument),
- output_origIC(string("--Oorig"), false,
- string("\toutput the original ICs"),
- false, no_argument),
- output_mean(string("--Omean"), false,
- string("\toutput mean volume\n"),
- false, no_argument),
- verbose(string("-v,--verbose"), false,
- string("switch on diagnostic messages"),
- false, no_argument),
- vers(string("-V,--version"), false,
- string("prints version information"),
- false, no_argument),
- copyright(string("--copyright"), false,
- string("prints copyright information"),
- false, no_argument),
- help(string("-h,--help"), false,
- string("prints this help message"),
- false, no_argument),
- debug(string("--debug"), false,
- string(" switch on debug messages"),
- false, no_argument),
- guessfname(string("--init_ica"), string(""),
- string("file name of FEAT paradigm file (design.mat) for ICA initialisation"),
- false, requires_argument, false),
- paradigmfname(string("--init_pca"), string(""),
- string("file name of FEAT paradigm file (design.mat) for PCA initialisation"),
- false, requires_argument, false),
- axials_str(string("--report_maps"), string(" -s 2 -A 950 "),
- string("control string for spatial map images (see slicer)"),
- false, requires_argument),
- dummy(string("--dummy"), 0,
- string("number of dummy volumes"),
- false, requires_argument,false),
- repeats(string("--repeats"), 1,
- string("number of repeats (multistart)"),
- false, requires_argument, false),
- seed(string("--seed"), -1,
- string("integer seed for random number generator within melodic"),
- false, requires_argument, false),
- nlconst1(string("--nl1,--nlconst1"), 1.0,
- string("nonlinear constant 1"),
- false, requires_argument, false),
- nlconst2(string("--nl2,--nlconst2"), 1.0,
- string("nonlinear constant 2"),
- false, requires_argument, false),
- smooth_probmap(string("--smooth_pm"), 0.0,
- string("width of smoothing kernel for probability maps"),
- false, requires_argument, false),
- insta_fn(string("--inc_fn"), string(""),
- string(" file name for instacorr volumes"),
- false, requires_argument, false),
- insta_idx(string("--inc_idx"), 1,
- string(" index for instacorr calculation"),
- false, requires_argument, false),
- insta_partial(string("--inc_partial"), true,
- string(" switch off partial correlation analysis for instacorrs"),
- false, requires_argument, false),
- insta_varnorm(string("--inc_varnorm"), true,
- string(" switch off varnorm for instantaneous correlations"),
- false, no_argument, false),
- remove_meanvol(string("--keep_meanvol"), true,
- string("do not subtract mean volume"),
- false, no_argument, false),
- remove_meantc(string("--remove_meantc"), false,
- string("remove mean time course"),
- false, no_argument, false),
- remove_endslices(string("--remEndslices"), false,
- string("delete end slices (motion correction artefacts)"),
- false, no_argument,false),
- rescale_nht(string("--rescale_nht"), true,
- string("switch off map rescaling after mixture-modelling"),
- false, no_argument,false),
- guess_remderiv(string("--remove_deriv"), false,
- string("removes every second entry in paradigm file (EV derivatives)"),
- false, no_argument, false),
- temporal(string("--temporal"), false,
- string("perform temporal ICA"),
- false, no_argument, false),
- retryfactor(string("--retryfactor"), float(0.95),
- string("multiplicative factor for determining new dim if estimated dim fails to converge"),
- false, requires_argument, false),
- econ(string("--econ"), 20000,
- string("set ctrl parameter for helperfns econ mode"),
- false, requires_argument, false),
- options(title, usageexmpl)
- {
- try {
- options.add(logdir);
- options.add(inputfname);
- options.add(outputfname);
- options.add(guessfname);
- options.add(maskfname);
- options.add(use_mask);
- options.add(update_mask);
- options.add(perf_bet);
- options.add(threshold);
- options.add(pca_dim);
- options.add(pca_est);
- options.add(joined_whiten);
- options.add(joined_vn);
- options.add(dr_pca);
- options.add(migp);
- options.add(migpN);
- options.add(migp_shuffle);
- options.add(migp_factor);
- options.add(dr);
- options.add(dr_out);
- options.add(vn_level);
- options.add(numICs);
- options.add(approach);
- options.add(nonlinearity);
- options.add(varnorm);
- options.add(varnorm2);
- options.add(pbsc);
- options.add(pspec);
- options.add(segment);
- options.add(tsmooth);
- options.add(epsilon);
- options.add(epsilonS);
- options.add(maxNumItt);
- options.add(maxRestart);
- options.add(rank1interval);
- options.add(mmthresh);
- options.add(perf_mm);
- options.add(ICsfname);
- options.add(filtermix);
- options.add(smodename);
- options.add(filter);
- options.add(genreport);
- options.add(guireport);
- options.add(bgimage);
- options.add(tr);
- options.add(logPower);
- options.add(addsigchng);
- options.add(allPPCA);
- options.add(varplots);
- options.add(varvals);
- options.add(fn_Tdesign);
- options.add(fn_Tcon);
- options.add(fn_TconF);
- options.add(fn_Sdesign);
- options.add(fn_Scon);
- options.add(fn_SconF);
- options.add(output_all);
- options.add(output_unmix);
- options.add(output_MMstats);
- options.add(output_pca);
- options.add(output_white);
- options.add(output_origIC);
- options.add(output_mean);
- options.add(verbose);
- options.add(vers);
- options.add(copyright);
- options.add(help);
- options.add(debug);
-
- options.add(guessfname);
- options.add(paradigmfname);
- options.add(axials_str);
- options.add(dummy);
- options.add(repeats);
- options.add(seed);
- options.add(nlconst1);
- options.add(nlconst2);
- options.add(smooth_probmap);
- options.add(insta_fn);
- options.add(insta_idx);
- options.add(insta_partial);
- options.add(insta_varnorm);
- options.add(remove_meanvol);
- options.add(remove_meantc);
- options.add(remove_endslices);
- options.add(rescale_nht);
- options.add(guess_remderiv);
- options.add(temporal);
- options.add(retryfactor);
- options.add(econ);
- }
- catch(X_OptionError& e) {
- options.usage();
- cerr << endl << e.what() << endl;
- }
- catch(std::exception &e) {
- cerr << e.what() << endl;
- }
-
- }
-}
-
-#endif
-
-
-
diff --git a/melpca.cc b/melpca.cc
deleted file mode 100644
index 1b77439832a1b3f5ad10c25ec0a1b48d7bfb65cb..0000000000000000000000000000000000000000
--- a/melpca.cc
+++ /dev/null
@@ -1,127 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melpca.cc - PCA and whitening
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "newimage/newimageall.h"
-#include "utils/log.h"
-#include "meloptions.h"
-#include "meldata.h"
-#include "melodic.h"
-#include "newmatap.h"
-#include "newmatio.h"
-#include "melpca.h"
-#include "melhlprfns.h"
-#include "libprob.h"
-
-using namespace Utilities;
-using namespace NEWIMAGE;
-
-namespace Melodic{
-
- void MelodicPCA::perf_pca(Matrix& in, Matrix& weights){
- message("Starting PCA ... ");
-
- SymmetricMatrix Corr;
- Matrix tmpE;
- RowVector tmpD, AdjEV, PercEV;
-
- if(opts.paradigmfname.value().length()>0)
- {
- basicGLM tmpglm;
- tmpglm.olsfit(in,melodat.get_param(),IdentityMatrix(melodat.get_param().Ncols()));
- std_pca(tmpglm.get_residu(),weights,Corr,tmpE,tmpD);
- }
- else{
- std_pca(in,weights,Corr,tmpE,tmpD);
- }
-
- if(opts.tsmooth.value()){
- message(endl << " temporal smoothing of Eigenvectors " << endl);
- tmpE=smoothColumns(tmpE);
- }
-
- adj_eigspec(tmpD,AdjEV,PercEV);
- melodat.set_pcaE(tmpE);
- melodat.set_pcaD(tmpD);
- melodat.set_EVP(PercEV);
- melodat.set_EV((AdjEV));
- write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV);
-
- message("done" << endl);
- }
-
- void MelodicPCA::perf_white(){
- int N = melodat.get_pcaE().Ncols();
- if(opts.pca_dim.value() > N){
- message("dimensionality too large - using -dim " << N <<
- " instead " << endl);
- opts.pca_dim.set_T(N);
- }
- if(opts.pca_dim.value() < 0){
- if(opts.remove_meanvol.value()){
- opts.pca_dim.set_T(N-2);
- }else{
- opts.pca_dim.set_T(N-1);
- }
- }
- if(opts.pca_dim.value() ==0){
- opts.pca_dim.set_T(pcadim());
- if(melodat.get_Data().Nrows()<20){
- opts.pca_dim.set_T(N-2);
- message("too few data points for full estimation, using "
- << opts.pca_dim.value() << " instead"<< endl);
- }
- }
- if(opts.approach.value()==string("jade") && opts.pca_dim.value() > 30){
- message("dimensionality too large for jade estimation - using --dim 30 instead" << endl);
- opts.pca_dim.set_T(30);
- }
-
- message("Start whitening using "<< opts.pca_dim.value()<<" dimensions ... " << endl);
- Matrix tmpWhite;
- Matrix tmpDeWhite;
-
- float varp = 1.0;
-
- if (opts.paradigmfname.value().length()>0)
- varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
- melodat.get_EVP(),opts.pca_dim.value(),melodat.get_param(),melodat.get_paramS(),tmpWhite,tmpDeWhite);
- else
- varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
- melodat.get_EVP(),opts.pca_dim.value(),tmpWhite,tmpDeWhite);
-
- melodat.set_white(tmpWhite);
- melodat.set_dewhite(tmpDeWhite);
- message(" retaining "<< varp*100 <<" percent of the variability " << endl);
- message(" ... done"<< endl << endl);
- }
-
- int MelodicPCA::pcadim()
- {
- message("Estimating the number of sources from the data (PPCA) ..." << endl);
-
- ColumnVector PPCA; RowVector AdjEV, PercEV;
- int res = ppca_dim(melodat.get_Data(),melodat.get_RXweight(), PPCA,
- AdjEV, PercEV, melodat.get_resels(), opts.pca_est.value());
-
- write_ascii_matrix(logger.appendDir("PPCA"),PPCA);
- write_ascii_matrix(logger.appendDir("eigenvalues_adjusted"),AdjEV.t());
- write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV.t());
-
- melodat.set_EVP(PercEV);
- melodat.set_EV(AdjEV);
- melodat.set_PPCA(PPCA);
-
- return res;
- }
-
-}
-
-
diff --git a/melpca.h b/melpca.h
deleted file mode 100644
index 346e9a5ba4d24aac6652b895b2d12a12b22f0b10..0000000000000000000000000000000000000000
--- a/melpca.h
+++ /dev/null
@@ -1,55 +0,0 @@
- /* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melpca.h - PCA and whitening
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#ifndef __MELODICPCA_h
-#define __MELODICPCA_h
-#include "newimage/newimageall.h"
-#include "utils/log.h"
-#include "meloptions.h"
-#include "meldata.h"
-#include "melodic.h"
-//#include "melreport.h"
-#include "newmatap.h"
-#include "newmatio.h"
-
-using namespace Utilities;
-using namespace NEWIMAGE;
-
-namespace Melodic{
-
- class MelodicReport;
-
- class MelodicPCA{
- public:
- MelodicPCA(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger,
- MelodicReport &preport):
- melodat(pmelodat),
- opts(popts),
- logger(plogger),
- report(preport){}
-
- void perf_pca(Matrix& in, Matrix& weights);
- inline void perf_pca(){
- perf_pca(melodat.get_Data(),melodat.get_RXweight());
- }
- void perf_white();
-
- private:
- MelodicData &melodat;
- MelodicOptions &opts;
- Log &logger;
- __attribute__((unused)) MelodicReport &report;
-
- int pcadim();
- };
-}
-
-#endif
diff --git a/melreport.cc b/melreport.cc
deleted file mode 100644
index 417f3ef78497e2d93826c2f13e58526a5f8e276e..0000000000000000000000000000000000000000
--- a/melreport.cc
+++ /dev/null
@@ -1,807 +0,0 @@
-/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
-
- melreport.cc - report generation
-
- Christian F. Beckmann, FMRIB Analysis Group
-
- Copyright (C) 1999-2013 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "newimage/newimageall.h"
-#include "utils/log.h"
-#include "melreport.h"
-#include "melhlprfns.h"
-#include "miscmaths/miscprob.h"
-
-namespace Melodic{
-
- void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim, Matrix ICstats){
-
- if( bool(opts.genreport.value()) ){
- addlink(mmodel.get_prefix()+".html",num2str(cnum));
- IChtml.setDir(report.getDir(),mmodel.get_prefix()+".html");
-
- {//start IC page
- IChtml << "<HTML><HEAD><link REL=stylesheet TYPE=text/css href=file:" +
- (string) getenv("FSLDIR") +"/doc/fsl.css>" << endl
- << "<style type=\"text/css\">OBJECT { width: 100% }</style>"
- << "<TITLE>FSL</TITLE></HEAD>" << endl
- << "<IFRAME height=" << int(melodat.get_numfiles()/30 + 1)*50
- <<"px width=100% src=nav.html frameborder=0></IFRAME><BR>"<< endl
- << "<Center>" << endl;
-
- if(cnum>1)
- IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
- <<".html\"><</a> - ";
- else
- IChtml << "<a href=\"" << string("IC_")+num2str(melodat.get_mix().Ncols())
- <<".html\"><</a> - ";
-
- if(cnum<dim)
- IChtml << "<a href=\"" << string("IC_")+num2str(cnum+1)
- <<".html\">></a>";
- else
- IChtml << "<a href=\"" << string("IC_")+num2str(1)
- <<".html\">></a>";
- IChtml << "<p><H3>MELODIC Component " << num2str(cnum)
- << "<br></b></H3><hr><p>" << endl;
- }
- {//output IC stats
- if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
- IChtml << fixed << setprecision(2) << std::abs(ICstats(cnum,1)) << " % of explained variance";
- if(ICstats.Ncols()>1)
- IChtml << "; " << std::abs(ICstats(cnum,2)) << " % of total variance";
- if(ICstats.Ncols()>2&&opts.addsigchng.value()){
- IChtml << "<p>" <<endl;
- IChtml << " " << ICstats(cnum,3) << " % signal change (pos peak voxel); " << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
- }
- IChtml << "<hr><p>" << endl;
- }
- }
-
- volume4D<float> tempVol;
- if(mmodel.get_threshmaps().Storage()>0&&
- (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols()))
- {//Output thresholded IC map
- tempVol.setmatrix(mmodel.get_threshmaps().Row(1),melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
-
- volume<float> newvol;
-
- miscpic newpic;
-
- float map1min = std::max((map1 + binarise(tempVol[0],tempVol[0].min(),
- float(0.0)) * map1.max()).min(),float(0.001));
- float map1max = std::max(map1.max(),float(0.001));
- float map2min = std::min(map2.min(),float(-0.001));
- float map2max = std::min((map2 + binarise(tempVol[0],float(0.0),
- tempVol[0].max()) * map2.min()).max(),float(-0.001));
-
- newpic.overlay(newvol, melodat.get_bg(), map1, map2,
- melodat.get_bg().percentile(0.02),
- melodat.get_bg().percentile(0.98),
- map1min, map1max, map2max, map2min,
- 0, 0);
- char instr[10000];
-
- sprintf(instr," ");
- strcat(instr,axials_instr.c_str());
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+
- "_thresh.png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - thresholded IC map ") + mmodel.get_infstr(1));
- newpic.set_cbar(string("ysb"));
-
- if((std::abs(map1.max()-map1.min())>0.01) || (std::abs(map2.max()-map2.min())>0.01))
- newpic.slicer(newvol, instr);
- else
- newpic.slicer(newvol, instr);
- IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
- IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
- +"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
- }
-
- {//plot time course
- IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
- miscplot newplot;
-
- Matrix tmptc = melodat.get_Tmodes(cnum-1).Column(1).t();
-
- newplot.col_replace(0,0xFF0000);
-
- newplot.add_label(string("IC ")+num2str(cnum)+" time course");
- //add GLM OLS fit
- if(melodat.Tdes.Storage()>0 &&
- melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){
- tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t();
- newplot.add_label("full model fit");
- }
-
- //add deviation around time course
- if(melodat.get_Tmodes(cnum-1).Ncols()>1 && opts.varplots.value()){
- Matrix tmp = stdev(melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t(),1);
- tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()+tmp;
- tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()-tmp;
- newplot.add_label("std error across subjects");
- newplot.col_replace(tmptc.Nrows()-1,0x808080);
- newplot.col_replace(tmptc.Nrows()-2,0x808080);
- }
-
- if(opts.tr.value()>0.0)
- newplot.add_xlabel(string("Time (seconds); TR = ")+
- float2str(opts.tr.value(),0,2,0)+" s");
- else
- newplot.add_xlabel(string("Time (TRs)"));
-
- newplot.add_ylabel("Normalised Response");
- newplot.set_yrange(tmptc.Row(1).Minimum()-0.05*(tmptc.Row(1).Maximum() -
- tmptc.Row(1).Minimum()),tmptc.Row(1).Maximum()+
- 0.05*(tmptc.Row(1).Maximum()-tmptc.Row(1).Minimum()));
- newplot.grid_swapdefault();
-
- newplot.set_xysize(750,150);
-
- newplot.timeseries(tmptc,
- report.appendDir(string("t")+num2str(cnum)+".png"),
- string("Timecourse No. ")+num2str(cnum),
- opts.tr.value(),150,12,1,false);
-
- if(melodat.get_Tmodes(cnum-1).Ncols()>1)
- tmptc &= melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t();
- write_ascii_matrix(report.appendDir(string("t")
- +num2str(cnum)+".txt"),tmptc.t());
- IChtml << "<A HREF=\"" << string("t")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
-
- if(melodat.get_numfiles()>1 && melodat.explained_var.Storage()>0
- && melodat.explained_var.Ncols()>=cnum && opts.varvals.value())
- IChtml << "Rank-1 approximation of individual time courses explains "
- << std::abs(melodat.explained_var(cnum)) << "% of variance.<p>" << endl;
- }//time series plot
-
- if(!opts.pspec.value())
- {//plot frequency
- miscplot newplot;
- RowVector empty(1);
- empty = 0.0;
- int fact = int(std::pow(10.0,int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
- if(opts.logPower.value())
- newplot.add_ylabel(string("log-Power"));
- else
- newplot.add_ylabel(string("Power"));
-
- Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1).Column(1), opts.logPower.value());
-
- newplot.set_Ylabel_fmt("%.0f");
- newplot.set_yrange(0.0,1.02*fmixtc.Maximum());
- newplot.grid_swapdefault();
- newplot.set_xysize(750,150);
-
- if(opts.tr.value()>0.0){
- newplot.add_xlabel(string("Frequency (in Hz / ")+num2str(fact)+ " )");
- newplot.timeseries(empty | fmixtc.t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string("Powerspectrum of timecourse"),
- fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()), 150,0,2);
- }else{
- newplot.add_xlabel(string("Frequency (in cycles); ")
- +"frequency(Hz)=cycles/("
- +num2str(melodat.get_Tmodes(0).Nrows())
- +"* TR); period(s)=("
- +num2str(melodat.get_Tmodes(0).Nrows())
- +"* TR)/cycles"
- );
- newplot.timeseries(fmixtc.t(),
- report.appendDir(string("f")+num2str(cnum)+".png"),
- string("Powerspectrum of timecourse"));
- }
- write_ascii_matrix(report.appendDir(string("f")
- +num2str(cnum)+".txt"), fmixtc);
- IChtml << "<A HREF=\"" << string("f")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
- }//frequency plot
- {//add T-mode GLM F-stats for full model fit & contrasts
- if(melodat.Tdes.Storage() > 0 &&
- melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){
- IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
- "<CAPTION><EM> <b>GLM (OLS) on time series </b></EM></CAPTION>" << endl
- << "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
- if(melodat.Tcon.Storage() > 0)
- IChtml << "<TH ><EM>Contrasts</EM>"<<endl;
- IChtml << "<TR><TD><TABLE border=0><TR><TD align=right>" << endl;
- for(int ctr=1;ctr <= melodat.Tdes.Ncols();ctr++)
- IChtml << " PE(" <<num2str(ctr)+"): <br>" << endl;
- IChtml << "<TD align=right>" << endl;
- for(int ctr=1;ctr <= melodat.Tdes.Ncols();ctr++)
- IChtml << melodat.glmT.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
- IChtml << "</TABLE>" <<
- " <TD align=center> F = "<< melodat.glmT.get_f_fmf().Column(cnum) <<
- " <BR> dof1 = " << melodat.Tdes.Ncols() << "; dof2 = "
- << melodat.glmT.get_dof() << "<BR>" <<endl;
- if(melodat.glmT.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
- IChtml << fixed << setprecision(5) <<"<b> p < " << melodat.glmT.get_pf_fmf().Column(cnum) <<
- "<BR> (uncorrected for #comp.)<b></TD>" << endl;
- else
- IChtml << fixed << setprecision(5) << " p < " <<
- melodat.glmT.get_pf_fmf().Column(cnum) <<
- "<BR> (uncorrected for #comp.)</TD>" << endl;
- if(melodat.Tcon.Storage() > 0 &&
- melodat.Tdes.Ncols() == melodat.Tcon.Ncols()){
- IChtml << fixed << setprecision(2) << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
- for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
- IChtml << "COPE(" << num2str(ctr) << "): <br>" << endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
- IChtml <<" z = <BR>" <<endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
- IChtml << melodat.glmT.get_z().Column(cnum).Row(ctr) <<";<BR>" <<endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
- if(melodat.glmT.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
- IChtml << fixed << setprecision(5) << "<b> p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
- "</b><BR>" << endl;
- else
- IChtml << fixed << setprecision(5) <<" p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
- "<BR>" << endl;
- IChtml << "</TABLE></td></tr>" << endl;
- }
- }
- IChtml << "</TABLE><p>" << endl;
- }
-
- if(cnum <= (int)melodat.get_Smodes().size())
- {//plot subject mode
-
- Matrix smode;
- smode = melodat.get_Smodes(cnum-1);
-
- if(smode.Nrows() > 1){
- IChtml << "<hr><H3> Sessions/Subjects mode </H3><p>" << endl <<endl;
- miscplot newplot;
-
- //add GLM OLS fit
- //newplot.setscatter(smode,2);
-
- if(melodat.Sdes.Storage() > 0&&
- melodat.glmS.get_beta().Nrows() == melodat.Sdes.Ncols()){
- smode |= melodat.Sdes * melodat.glmS.get_beta().Column(cnum);
- newplot.add_label(string("IC ")+num2str(cnum)+" subject/session-mode");
- newplot.add_label("full model fit");
- }
- newplot.grid_swapdefault();
- newplot.set_Ylabel_fmt("%.2f");
- newplot.add_xlabel(" Subject Number");
-
-
- newplot.set_xysize(750,150);
- newplot.timeseries(smode.t(),
- report.appendDir(string("s")+num2str(cnum)+".png"),
- string("Subject/Session mode No. ") + num2str(cnum));
- newplot.clear_xlabel();
- newplot.clear_labels();
- newplot.set_xysize(120,200);
- newplot.set_minmaxscale(1.1);
-
- newplot.boxplot((Matrix)smode.Column(1),
- report.appendDir(string("b")+num2str(cnum)+".png"),
- string("Subject/Session mode"));
- write_ascii_matrix(report.appendDir(string("s")
- +num2str(cnum)+".txt"), smode);
- IChtml << "<A HREF=\"" << string("s")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("s")+num2str(cnum)+".png\"></A>" << endl;
- IChtml << "<A HREF=\"" << string("s")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("b")+num2str(cnum)+".png\"></A><p>" << endl;
- }
- {//add S-mode GLM F-stats for full model fit & contrasts
- if(melodat.Sdes.Storage() > 0 &&
- melodat.glmS.get_beta().Nrows() == melodat.Sdes.Ncols()){
- IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
- "<CAPTION><EM> <b>GLM (OLS) on subject/session-mode </b></EM></CAPTION>" << endl
- << "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
- if(melodat.Scon.Storage() > 0)
- IChtml << "<TH ><EM>Contrasts</EM>"<<endl;
- IChtml << "<TR><TD><TABLE border=0><TR><TD align=right>" << endl;
- for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
- IChtml << " PE(" <<num2str(ctr)+"): <br>" << endl;
- IChtml << "<TD align=right>" << endl;
- for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
- IChtml << melodat.glmS.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
- IChtml << "</TABLE>" <<
- " <TD align=center> F = "<< melodat.glmS.get_f_fmf().Column(cnum) <<
- " <BR> dof1 = " << melodat.Sdes.Ncols() << "; dof2 = "
- << melodat.glmS.get_dof() << "<BR>" <<endl;
- if(melodat.glmS.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
- IChtml << fixed << setprecision(5) <<"<b> p < " << melodat.glmS.get_pf_fmf().Column(cnum) <<
- "<BR> (uncorrected for #comp.)<b></TD>" << endl;
- else
- IChtml << fixed << setprecision(5) << " p < " <<
- melodat.glmS.get_pf_fmf().Column(cnum) <<
- "<BR> (uncorrected for #comp.)</TD>" << endl;
- if(melodat.Scon.Storage() > 0 && melodat.Sdes.Storage() > 0 &&
- melodat.Sdes.Ncols() == melodat.Scon.Ncols()){
- IChtml << fixed << setprecision(2) << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- IChtml << "COPE(" << num2str(ctr) << "): <br>" << endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- IChtml <<" z = <BR>" <<endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- IChtml << melodat.glmS.get_z().Column(cnum).Row(ctr) <<";<BR>" <<endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- if(melodat.glmS.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
- IChtml << fixed << setprecision(5) << "<b> p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
- "</b><BR>" << endl;
- else
- IChtml << fixed << setprecision(5) <<" p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
- "<BR>" << endl;
- IChtml << "</TABLE></td></tr>" << endl;
- }
- }
- IChtml << "</TABLE><p>" << endl;
- }
- }//subject mode plot
-
- if(mmodel.get_threshmaps().Storage()>0&&
- (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
- (mmodel.get_threshmaps().Nrows()>1))
- {//Output other thresholded IC map
-
- for(int tctr=2; tctr<=mmodel.get_threshmaps().Nrows(); tctr++){
- tempVol.setmatrix(mmodel.get_threshmaps().Row(tctr),melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
-
- volume<float> newvol;
- miscpic newpic;
-
- float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
- float(0.0)) * map1.max()).min();
- float map1max = map1.max();
- float map2min = map2.min();
- float map2max = (map2 + binarise(tempVol[0],float(0.0),
- tempVol[0].max()) * map2.min()).max();
-
- //cerr << endl << map1min << " " << map1max << endl
- // << map2min << " " << map2max << endl;
-
- // if(map1.max()-map1.min()>0.01)
- newpic.overlay(newvol, melodat.get_bg(), map1, map2,
- melodat.get_bg().percentile(0.02),
- melodat.get_bg().percentile(0.98),
- map1min, map1max, map2max, map2min,
- 0, 0);
-
- char instr[10000];
-
- sprintf(instr," ");
- strcat(instr,axials_instr.c_str());
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+"_thresh"+
- num2str(tctr)+".png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - thresholded IC map ("+
- num2str(tctr)+") ")+
- mmodel.get_infstr(tctr));
- newpic.set_cbar(string("ysb"));
- //cerr << instr << endl;
- newpic.slicer(newvol, instr);
-
- IC_rep_det(mmodel, cnum, dim);
- IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
- IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
- +"_thresh"+num2str(tctr)+".png\" ALT=\"MMfit\"></A><p>" << endl;
- }
- }
-
- { //finish IC page
- IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
- << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
- << "FMRIB Software Library</A>.</FONT></Center>" << endl
- << "</BODY></HTML>" << endl;
- } //finish IC page
- IC_rep_det(mmodel, cnum, dim);
- }
- }
-
- void MelodicReport::IC_rep_det(MelGMix &mmodel, int cnum, int dim){
- if( bool(opts.genreport.value()) ){
-
- {//start IC2 page
- IChtml2.setDir(report.getDir(),mmodel.get_prefix()+"_MM.html");
- IChtml2 << "<HTML><HEAD><link REL=stylesheet TYPE=text/css href=file:" +
- (string) getenv("FSLDIR") +"/doc/fsl.css>" << endl
- << "<style type=\"text/css\">OBJECT { width: 100% }</style>"
- << "<TITLE>FSL</TITLE></HEAD>" << endl
- << "<IFRAME height="<< int(melodat.get_numfiles()/30 + 1)*50
- <<"px width=100% src=nav.html frameborder=0></IFRAME><p>"<< endl
- << "<CENTER>";
- if(cnum>1)
- IChtml2 << "<b><a href=\"" << string("IC_")+num2str(cnum-1)
- <<"_MM.html\"><</a> - ";
- else
- IChtml2 << "<b><a href=\"" << string("IC_")+num2str(melodat.get_mix().Ncols())
- <<"_MM.html\"><</a> - ";
- // IChtml << "<a href=\"00index.html\"> index </a>" ;
-
- if(cnum<dim)
- IChtml2 << "<a href=\"" << string("IC_")+num2str(cnum+1)
- <<"_MM.html\">></a>";
- else
- IChtml2 << "<a href=\"" << string("IC_")+num2str(1)
- <<"_MM.html\">></a>";
- IChtml2 << "<p><H3>Component " << num2str(cnum)
- << " Mixture Model fit <br></b></H3><hr><p>" << endl;
- }
-
- volume4D<float> tempVol;
-
- if(melodat.get_IC().Storage()>0)
- {//Output raw IC map
-
- // tempVol.setmatrix(melodat.get_IC().Row(cnum),
- // melodat.get_mask());
- tempVol.setmatrix(mmodel.get_data(),
- melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0),
- tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(),
- float(-0.0));
-
- volume<float> newvol;
- miscpic newpic;
-
- // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
- // float(0.0)) * map1.max()).robustmin();
- float map1max = map1.percentile(0.99);
- float map2min = map2.percentile(0.01);
- //float map2max = (map2 + binarise(tempVol[0],float(0.0),
- // tempVol[0].max()) * map2.min()).robustmax();
-
- newpic.overlay(newvol, melodat.get_bg(), map1, map2,
- float(0.0),
- float(0.0),
- float(0.01), map1max, float(-0.01), map2min,
- 0, 0);
-
- char instr[10000];
-
- sprintf(instr," ");
- strcat(instr,axials_instr.c_str());
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+
- ".png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - raw Z transformed IC map (1 - 99 percentile)"));
- newpic.set_cbar(string("ysb"));
-
- newpic.slicer(newvol, instr);
- }
- IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
- IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
- ".png\"><A><p>" << endl;
-
- if(mmodel.get_probmap().Storage()>0&&
- (mmodel.get_probmap().Ncols() == mmodel.get_data().Ncols())&&
- (mmodel.get_probmap().Nrows() == mmodel.get_data().Nrows()))
- {//Output probmap
- tempVol.setmatrix(mmodel.get_probmap(),melodat.get_mask());
-
- volume<float> map;
- map = tempVol[0];
-
- volume<float> newvol;
- miscpic newpic;
-
- newpic.overlay(newvol, melodat.get_bg(), map, map,
- melodat.get_bg().percentile(0.02),
- melodat.get_bg().percentile(0.98),
- float(0.1), float(1.0), float(0.0), float(0.0),
- 0, 0);
-
- char instr[10000];
-
- sprintf(instr," ");
- strcat(instr,"-l render1 ");
- strcat(instr,axials_instr.c_str());
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+
- "_prob.png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - Mixture Model probability map"));
-
- newpic.set_cbar(string("y"));
- newpic.slicer(newvol, instr);
-
- IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
- IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
- "_prob.png\">" << endl;
-
- IChtml2 << "</A><p>" << endl;
- }
-
- RowVector dat = mmodel.get_data().Row(1);
- if(dat.Maximum()>dat.Minimum())
- {//Output GGM/GMM fit
- miscplot newplot;
-
- if(mmodel.get_type()=="GGM"){
- newplot.add_label("IC map histogram");
- newplot.add_label("full GGM fit");
- newplot.add_label("background Gaussian");
- newplot.add_label("Gamma distributions");
- newplot.gmmfit(mmodel.get_data().Row(1),
- mmodel.get_means(),
- mmodel.get_vars(),
- mmodel.get_pi(),
- report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
- string(mmodel.get_prefix() +
- " Gaussian/Gamma Mixture Model("+num2str(mmodel.mixtures())+") fit"),
- true, float(0.0), float(0.0));
-
- }
- else{
- newplot.add_label("IC map histogram");
- newplot.add_label("full GMM fit");
- newplot.add_label("individual Gaussians");
- newplot.gmmfit(mmodel.get_data().Row(1),
- mmodel.get_means(),
- mmodel.get_vars(),
- mmodel.get_pi(),
- report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
- string(mmodel.get_prefix() +
- " Gaussian Mixture Model("+num2str(mmodel.mixtures())+") fit"),
- false, float(0.0), float(2.0));
-
- }
-
- // IChtml2 << "<A HREF=\"" +mmodel.get_prefix()+"_MMfit_detail.png\"> ";
- IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
- "_MMfit.png\"><p>" << endl;
- } //GGM/GMM plot
-
- {//MM parameters
- IChtml2 << "<br> " << mmodel.get_prefix()
- << " Mixture Model fit <br>" << endl
- << "<br> Means : " << mmodel.get_means() << endl
- << "<br> Vars : " << mmodel.get_vars() << endl
- << "<br> Prop. : " << mmodel.get_pi() << endl;
- }
-
- { //finish IC2 page
- IChtml2<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
- << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
- << "FMRIB Software Library</A>.</FONT></CENTER>" << endl
- << "</BODY></HTML>" << endl;
- } //finish IC2 page
- }
- }
-
- void MelodicReport::IC_simplerep(string prefix, int cnum, int dim){
- if( bool(opts.genreport.value()) ){
- addlink(prefix+".html",num2str(cnum));
- IChtml.setDir(report.getDir(),prefix+".html");
-
- {//start IC page
-
- IChtml << "<HTML> " << endl
- << "<TITLE>MELODIC Component " << num2str(cnum)
- << "</TITLE>" << endl
- << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
- << "/doc/images/fsl-bg.jpg\">" << endl
- << "<hr><CENTER><H1>MELODIC Component " << num2str(cnum)
- << "</H1>"<< endl;
-
- if(cnum>1)
- IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
- <<".html\">previous</a> - ";
-
- IChtml << "<a href=\"00index.html\"> index </a>" ;
-
- if(cnum<dim)
- IChtml << " - <a href=\"" << string("IC_")+num2str(cnum+1)
- <<".html\">next</a><p>";
-
- IChtml << "<hr><p>" << endl;
- }
-
- volume4D<float> tempVol;
-
- if(melodat.get_IC().Storage()>0)
- {//Output raw IC map
-
- tempVol.setmatrix(melodat.get_IC().Row(cnum),
- melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0),
- tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(),
- float(-0.0));
-
- volume<float> newvol;
- miscpic newpic;
-
- // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
- // float(0.0)) * map1.max()).robustmin();
- float map1max = map1.percentile(0.99);
- float map2min = map2.percentile(0.01);
- //float map2max = (map2 + binarise(tempVol[0],float(0.0),
- // tempVol[0].max()) * map2.min()).robustmax();
-
- newpic.overlay(newvol, melodat.get_bg(), map1, map2,
- float(0.0),
- float(0.0),
- float(0.01), map1max, float(-0.01), map2min,
- 0, 0);
-
- char instr[10000];
-
- sprintf(instr," ");
- strcat(instr,axials_instr.c_str());
- strcat(instr,string(report.appendDir(prefix+
- ".png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - raw Z transformed IC map (1 - 99 percentile)"));
- newpic.set_cbar(string("ysb"));
-
- newpic.slicer(newvol, instr);
- }
-
- IChtml << "<img BORDER=0 SRC=\""+ prefix+
- ".png\"><p>" << endl;
-
- {//plot time course
- miscplot newplot;
-
- if(opts.tr.value()>0.0)
- newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
- report.appendDir(string("t")+
- num2str(cnum)+".png"),
- string("Timecourse (in seconds); TR = ")+
- float2str(opts.tr.value(),0,2,0)+" s",
- opts.tr.value(),150,4,1);
- else
- newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
- report.appendDir(string("t")+
- num2str(cnum)+".png"),
- string("Timecourse (in TRs)"));
- write_ascii_matrix(report.appendDir(string("t")
- +num2str(cnum)+".txt"),
- melodat.get_Tmodes(cnum-1));
- IChtml << "<A HREF=\"" << string("t")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
- }//time series plot
-
- {//plot frequency
- miscplot newplot;
- int fact = int(std::pow(10.0,
- int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
-
- if(opts.tr.value()>0.0)
- newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string("FFT of timecourse (in Hz / ") +
- num2str(fact)+")",
- fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()),
- 150,0,2);
- else
- newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string(string("FFT of timecourse (in cycles); ")
- +"frequency(Hz)=cycles/("
- +num2str(melodat.get_Tmodes(0).Nrows())
- +"* TR); period(s)=("
- +num2str(melodat.get_Tmodes(0).Nrows())
- +"* TR)/cycles"));
- write_ascii_matrix(report.appendDir(string("f")
- +num2str(cnum)+".txt"),
- melodat.get_Tmodes(cnum-1));
- IChtml << "<A HREF=\"" << string("f")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
- }//frequency plot
-
- { //finish IC page
- IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
- << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
- << "FMRIB Software Library</A>.</FONT>" << endl
- << "</BODY></HTML>" << endl;
- } //finish IC page
- }
- }
-
- void MelodicReport::PPCA_rep(){
-
- {//plot time course
- report << "<p><hr><b>PCA estimates </b> <p>" << endl;
-
- Matrix what;
- miscplot newplot;
-
- what = melodat.get_EV();
- what &= melodat.get_EVP();
-
- newplot.add_label("ordered Eigenvalues");
- newplot.add_label("% of expl. variance");
-
- if(melodat.get_PPCA().Storage()>0){
- what = what.Columns(1,melodat.get_PPCA().Nrows());
- if(opts.allPPCA.value()&&melodat.get_PPCA().Ncols()==7){
- what &= melodat.get_PPCA().Columns(3,7).t();
- newplot.add_label("Laplace");
- newplot.add_label("BIC");
- newplot.add_label("MDL");
- newplot.add_label("RRN");
- newplot.add_label("AIC");
- }else{
- what &= melodat.get_PPCA().Column(1).t();
- newplot.add_label("dim. estimate");
- }
- }
-
- newplot.set_Ylabel_fmt("%.2f");
- newplot.add_xlabel("Number of included components");
- newplot.set_yrange(0.0,1.02);
- newplot.grid_swapdefault();
- newplot.timeseries(what,
- report.appendDir("EVplot.png"),
- string("Eigenspectrum Analysis"),
- 0,450,4,0);
-
- report << "<img BORDER=0 SRC=\"EVplot.png\"><p>" << endl;
- }//time series plot
- }
-
- void MelodicReport::Smode_rep(){
- if(melodat.get_Smodes().size()>0){
- report << "<p><hr><b>TICA Subject/Session modes </b> <br>" << endl;
- miscplot newplot;
- report << "Boxplots show the relative response amplitudes across the "
- << "session/subject domain (" << melodat.get_numfiles()
- << " input files). Components have been sorted in decreasing order of "
- << " the median response per component. <br><br>";
-
- outMsize("Smode.at(0)", melodat.get_Smodes().at(0));
- Matrix allmodes = melodat.get_Smodes().at(0);
- for(int ctr = 1; ctr < (int)melodat.get_Smodes().size();++ctr)
- allmodes |= melodat.get_Smodes().at(ctr);
-
- outMsize("allmodes", allmodes);
- newplot.add_xlabel("Component No.");
- newplot.add_ylabel("");
- if(allmodes.Ncols()<100)
- newplot.set_xysize(100+30*allmodes.Ncols(),300);
- else
- newplot.set_xysize(1200,300);
- newplot.boxplot(allmodes,report.appendDir(string("bp_Smodes.png")),
- string("Subject/Session modes"));
- report << "<img BORDER=0 SRC=\"bp_Smodes.png\"><p>" << endl;
- }
- }
-}