diff --git a/+nla/+edge/+test/SandwichEstimator.m b/+nla/+edge/+test/SandwichEstimator.m index a16d90b5..b1bbf1b4 100755 --- a/+nla/+edge/+test/SandwichEstimator.m +++ b/+nla/+edge/+test/SandwichEstimator.m @@ -150,13 +150,14 @@ stdErrInput.scanMetadata = scanMetadata; stdErrInput.residual = residual; stdErrInput.pinvDesignMtx = pinv(designMtx); - - %sweRes.stdError = stdErrCalcObj.calculate(stdErrInput); - stdError = stdErrCalcObj.calculate(stdErrInput); + stdErrInput.contrasts = input.contrasts; + + %change stdError to compute contrast SE + contrastSE = stdErrCalcObj.calculate(stdErrInput); contrastCalc = input.contrasts * regressCoeffs; - contrastSE = sqrt((input.contrasts.^2) * (stdError.^2)); + dof = obj.calcDegreesOfFreedom(designMtx); diff --git a/+nla/+gfx/PlotValue.m b/+nla/+gfx/PlotValue.m new file mode 100644 index 00000000..86b2f8f9 --- /dev/null +++ b/+nla/+gfx/PlotValue.m @@ -0,0 +1,5 @@ +classdef PlotValue + enumeration + PVALUE, STATISTIC + end +end \ No newline at end of file diff --git a/+nla/+gfx/ProbPlotMethod.m b/+nla/+gfx/ProbPlotMethod.m index d57c18f1..e646aece 100755 --- a/+nla/+gfx/ProbPlotMethod.m +++ b/+nla/+gfx/ProbPlotMethod.m @@ -1,5 +1,5 @@ classdef ProbPlotMethod enumeration - DEFAULT, LOG, NEGATIVE_LOG_10, STATISTIC + DEFAULT, LOG, NEGATIVE_LOG_10 end end \ No newline at end of file diff --git a/+nla/+helpers/+stdError/Guillaume.m b/+nla/+helpers/+stdError/Guillaume.m index 616f4501..a35e4101 100755 --- a/+nla/+helpers/+stdError/Guillaume.m +++ b/+nla/+helpers/+stdError/Guillaume.m @@ -15,7 +15,7 @@ methods - function stdError = calculate(obj, sweStdErrInput) + function contrastStdError = calculate(obj, sweStdErrInput) @@ -61,7 +61,9 @@ end - stdError = sqrt(betaCovar.v(betaCovar.getDiagElemIdxs,:)); + stdErr = sqrt(betaCovar.v(betaCovar.getDiagElemIdxs,:)); + + contrastStdError = sqrt((sweStdErrInput.contrasts.^2) * (stdErr.^2)); end diff --git a/+nla/+helpers/+stdError/Heteroskedastic.m b/+nla/+helpers/+stdError/Heteroskedastic.m index 2dc7174b..47e06899 100755 --- a/+nla/+helpers/+stdError/Heteroskedastic.m +++ b/+nla/+helpers/+stdError/Heteroskedastic.m @@ -6,7 +6,16 @@ methods - function stdErr = calculate(obj, sweStdErrInput) + function contrastStdErr = calculate(obj, sweStdErrInput) + + FORCE_USE_FAST_ALGO = true; + if FORCE_USE_FAST_ALGO + %There is a faster, but possibly larger memory + %implementation of this algorithm. Don't + fastAlgoObj = nla.helpers.stdError.Heteroskedastic_FAST(); + contrastStdErr = fastAlgoObj.calculate(sweStdErrInput); + return; + end %Calculation of standard error assuming heteroskedascticity @@ -31,6 +40,8 @@ stdErr(:,fcEdgeIdx) = sqrt(correctionFactor * diag(betaCovariance)); end + + contrastStdErr = sqrt((sweStdErrInput.contrasts.^2) * (stdErr.^2)); end diff --git a/+nla/+helpers/+stdError/Heteroskedastic_FAST.m b/+nla/+helpers/+stdError/Heteroskedastic_FAST.m index ef672a86..5490ba69 100755 --- a/+nla/+helpers/+stdError/Heteroskedastic_FAST.m +++ b/+nla/+helpers/+stdError/Heteroskedastic_FAST.m @@ -6,7 +6,7 @@ methods - function stdErr = calculate(obj, sweStdErrInput) + function contrastStdErr = calculate(obj, sweStdErrInput) %Computes Standard Error, but accelerated using assumption of %heteroskeadisticity for quicker computation % @@ -51,6 +51,8 @@ diagElemIdxsInFlatArr = 1:(numCovariates+1):numCovariates^2; stdErr = sqrt(betaCovarianceFlat(diagElemIdxsInFlatArr,:)); + contrastStdErr = sqrt((sweStdErrInput.contrasts.^2) * (stdErr.^2)); + end end diff --git a/+nla/+helpers/+stdError/Homoskedastic.m b/+nla/+helpers/+stdError/Homoskedastic.m index 0c7b0c1b..9c3851e0 100755 --- a/+nla/+helpers/+stdError/Homoskedastic.m +++ b/+nla/+helpers/+stdError/Homoskedastic.m @@ -6,7 +6,7 @@ methods - function stdErr = calculate(obj, sweStdErrInput) + function contrastStdErr = calculate(obj, sweStdErrInput) %Calculation of standard error assuming homoskedasticity %(errors are independent and identically distributed iid) @@ -20,7 +20,9 @@ meanSqErr = sum(residual.^2) ./ degOfFree; %In regression, divide by dof instead of number of data points (per wikipedia) - stdErr = sqrt(diag(pinvDesignMtx * pinvDesignMtx')*meanSqErr); + stdErr = sqrt(diag(pinvDesignMtx * pinvDesignMtx')*meanSqErr); + + contrastStdErr = sqrt((sweStdErrInput.contrasts.^2) * (stdErr.^2)); end diff --git a/+nla/+helpers/+stdError/UnconstrainedBlocks.m b/+nla/+helpers/+stdError/UnconstrainedBlocks.m index 682848e9..1c22ce06 100755 --- a/+nla/+helpers/+stdError/UnconstrainedBlocks.m +++ b/+nla/+helpers/+stdError/UnconstrainedBlocks.m @@ -1,69 +1,88 @@ classdef UnconstrainedBlocks < nla.helpers.stdError.AbstractSwEStdErrStrategy - properties - - SPARSITY_THRESHOLD = 0.2; - - end properties (SetAccess = protected) REQUIRES_GROUP = true; end methods - function stdErr = calculate(obj, sweStdErrInput) - %Computes Standard Error assuming unconstrained blocks - %Uses standard matrix multiplication to compute standard error, - %since it does not make assumption that V is sparse. - + function contrastStdErr = calculate(obj, sweStdErrInput) + + + %rename variables for readability + pinvDesignMtx = sweStdErrInput.pinvDesignMtx; + residual = sweStdErrInput.residual; groupIds = sweStdErrInput.scanMetadata.groupId; unqGrps = unique(groupIds); obj.throwErrorIfVEntirelyFull(unqGrps); - vSparsity = obj.computeVSparsity(groupIds); + [numCovariates, ~] = size(pinvDesignMtx); + [numObs, numFcEdges] = size(residual); - %Ben Kay 'Half Sandwich' algorithm seems to be at least as good - %or better than any of the other clever approaches so far. - %Might be able to beat it by using the clever approach and only - %computing the diagonal of covBat - FORCE_HALF_SW_ALGO = true; + stdErr = zeros(numCovariates, numFcEdges); - if FORCE_HALF_SW_ALGO - stdErrStrategy = nla.helpers.stdError.UnconstrainedBlocks_BenKay(); - elseif vSparsity <= obj.SPARSITY_THRESHOLD - stdErrStrategy = nla.helpers.stdError.UnconstrainedBlocks_Sparse(); + numNonzeroValuesInContrast = sum(sweStdErrInput.contrasts~=0); + if numNonzeroValuesInContrast > 1 + WALD_TEST = true; else - stdErrStrategy = nla.helpers.stdError.UnconstrainedBlocks_Dense(); + WALD_TEST = false; end - stdErr = stdErrStrategy.calculate(sweStdErrInput); - + if ~WALD_TEST + covB = zeros(numCovariates,numFcEdges); + else + covB = zeros(numCovariates,numCovariates,numFcEdges); + end + %NOTE, optimized from swe_block.m from Benjamin Kay + %if NOT WALD_TEST, can just compute diagonals of cov(B), which + %makes the original halfsandwich of + %pinvDesignMtx(:,subjThisGrp) * residual(subjThisGrp, fcIdx) + %into a [numCovars x 1] matrix, and then squaring for the + %elements and multiplying by the one non-zero value of the + %contrast + + if ~WALD_TEST + for grpIdx = 1:length(unqGrps) + thisGrpId = unqGrps(grpIdx); + subjThisGrp = groupIds == thisGrpId; + halfSandwich = pinvDesignMtx(:, subjThisGrp) * residual(subjThisGrp,:); + + covB = covB + halfSandwich .* halfSandwich; + + end + + stdErr(:) = sqrt(covB); + contrastStdErr = sqrt((sweStdErrInput.contrasts.^2) * (stdErr.^2)); + else + for grpIdx = 1:length(unqGrps) + + thisGrpId = unqGrps(grpIdx); + subjThisGrp = groupIds == thisGrpId; + halfSandwich = pinvDesignMtx(:, subjThisGrp) * residual(subjThisGrp,:); + + for fcEdgeIdx = 1:numFcEdges + covB(:,:,fcEdgeIdx) = covB(:,:,fcEdgeIdx) + ... + (halfSandwich(:,fcEdgeIdx) * halfSandwich(:,fcEdgeIdx)'); + end + end + + %Computation of contrast StdErr here + contrasts = sweStdErrInput.contrasts; + contrastStdErr = zeros(1,numFcEdges); + for fcEdgeIdx = 1:numFcEdges + contrastStdErr(fcEdgeIdx) = contrasts * covB(:,:,fcEdgeIdx) * contrasts'; + end + + end + end - - end - - methods (Access = protected) - - function fractionNonZero = computeVSparsity(obj, groupIds) - %Do quick check of how many elements of V we expect to be - %nonzero given the group Ids of our observations. - %This calculation will only be fast if V is sparse, so we - %should determine how full V will be and warn user if this - %method will be slow. - unqGrps = unique(groupIds); - countInGrps = histcounts(groupIds,[unqGrps;Inf]); - - numNonzeroElems = sum(countInGrps.^2); - totalElems = length(groupIds)*length(groupIds); - - fractionNonZero = numNonzeroElems / totalElems; - - end - + + methods (Access = protected) + function throwErrorIfVEntirelyFull(obj, uniqueGroupIds) %If V is entirely full, throw an error if length(uniqueGroupIds) == 1 @@ -76,30 +95,9 @@ function throwErrorIfVEntirelyFull(obj, uniqueGroupIds) end end - - function [groupedPinv, groupedResidual, groupIds] = reorderDataByGroup(obj, origPinvDesignMtx, origResidual, origGrps) - - [numCovariates, ~] = size(origPinvDesignMtx); - - %Group all data in one matrix and sort by group - %NOTE: need to use transpose of pinvDesignMatrix!!! - allData = [origGrps, origPinvDesignMtx', origResidual]; - sortedData= sortrows(allData,1); - - pinvColRange = 2:(numCovariates+1); - residualColRange = (numCovariates+2):(size(allData,2)); - - groupedPinvTpose = sortedData(:,pinvColRange); - groupedPinv = groupedPinvTpose'; - groupedResidual = sortedData(:,residualColRange); - groupIds = sortedData(:,1); - - end - - - + end - - + + end \ No newline at end of file diff --git a/+nla/+helpers/+stdError/UnconstrainedBlocks_BenKay.m b/+nla/+helpers/+stdError/UnconstrainedBlocks_BenKay.m deleted file mode 100755 index d743e76b..00000000 --- a/+nla/+helpers/+stdError/UnconstrainedBlocks_BenKay.m +++ /dev/null @@ -1,71 +0,0 @@ -classdef UnconstrainedBlocks_BenKay < nla.helpers.stdError.UnconstrainedBlocks - - methods - - function stdErr = calculate(obj, sweStdErrInput) - - - %rename variables for readability - pinvDesignMtx = sweStdErrInput.pinvDesignMtx; - residual = sweStdErrInput.residual; - groupIds = sweStdErrInput.scanMetadata.groupId; - unqGrps = unique(groupIds); - - obj.throwErrorIfVEntirelyFull(unqGrps); - - [numCovariates, ~] = size(pinvDesignMtx); - [numObs, numFcEdges] = size(residual); - - stdErr = zeros(numCovariates, numFcEdges); - - WALD_TEST = false; - - if ~WALD_TEST - covB = zeros(numCovariates,numFcEdges); - else - covB = zeros(numCovariates,numCovariates,numFcEdges); - end - - %NOTE, optimized from swe_block.m from Benjamin Kay - %if NOT WALD_TEST, can just compute diagonals of cov(B), which - %makes the original halfsandwich of - %pinvDesignMtx(:,subjThisGrp) * residual(subjThisGrp, fcIdx) - %into a [numCovars x 1] matrix, and then squaring for the - - - for grpIdx = 1:length(unqGrps) - - thisGrpId = unqGrps(grpIdx); - subjThisGrp = groupIds == thisGrpId; - halfSandwich = pinvDesignMtx(:, subjThisGrp) * residual(subjThisGrp,:); - - if WALD_TEST - - for fcEdgeIdx = 1:numFcEdges - covB(:,:,fcEdgeIdx) = covB(:,:,fcEdgeIdx) + ... - (halfSandwich(:,fcEdgeIdx) * halfSandwich(:,fcEdgeIdx)'); - end - - else - - covB = covB + halfSandwich .* halfSandwich; - - end - - end - - - if WALD_TEST - waldRunTime = toc - stdErr = rand(numCovariates, numFcEdges); - error('Not Implemented!'); - else - stdErr(:) = sqrt(covB); - end - - - end - - end - -end \ No newline at end of file diff --git a/+nla/+helpers/+stdError/UnconstrainedBlocks_Dense.m b/+nla/+helpers/+stdError/UnconstrainedBlocks_Dense.m deleted file mode 100755 index bf0474a5..00000000 --- a/+nla/+helpers/+stdError/UnconstrainedBlocks_Dense.m +++ /dev/null @@ -1,62 +0,0 @@ -classdef UnconstrainedBlocks_Dense < nla.helpers.stdError.UnconstrainedBlocks - - methods - - function stdErr = calculate(obj, sweStdErrInput) - %Computes Standard Error assuming unconstrained blocks - %Uses standard matrix multiplication to compute standard error, - %since it does not make assumption that V is sparse. - - - %rename variables for readability - pinvDesignMtx = sweStdErrInput.pinvDesignMtx; - residual = sweStdErrInput.residual; - groupIds = sweStdErrInput.scanMetadata.groupId; - unqGrps = unique(groupIds); - numGrps = length(unqGrps); - - obj.throwErrorIfVEntirelyFull(unqGrps); - - [numCovariates, numObs] = size(pinvDesignMtx); - [~,numFcEdges] = size(residual); - - %Preallocate size of stdErr output - stdErr = zeros(numCovariates, numFcEdges); - - %Reorder data so that grouped observations are together - [pinvDesignMtx_grp, residual_grp, groupIds_grp] = obj.reorderDataByGroup(pinvDesignMtx, residual, groupIds); - - %Determine which observations fit in which group once, outside - %of loop for fc edges - inGroupFlags = zeros(numObs, numGrps); - for grpIdx = 1:numGrps - thisGrpId = unqGrps(grpIdx); - inGroupFlags(:,grpIdx) = (groupIds_grp == thisGrpId); - end - - %For each fc edge, build the V and multiply it by - %pinvDesignMatrix on both sides - thisV = zeros(numObs, numObs); - - for fcIdx = 1:numFcEdges - - thisV(:,:) = 0; - for grpIdx = 1:numGrps - thisGrpFlags = logical(inGroupFlags(:,grpIdx)); - residThisGrp = residual_grp(thisGrpFlags,fcIdx); - thisGrpVBlock = residThisGrp * residThisGrp'; - thisV(thisGrpFlags,thisGrpFlags) = thisGrpVBlock; - end - - thisBetaCovar = pinvDesignMtx_grp * thisV * pinvDesignMtx_grp'; - stdErr(:,fcIdx) = sqrt(diag(thisBetaCovar)); - - end - - - end - - end - - -end \ No newline at end of file diff --git a/+nla/+helpers/+stdError/UnconstrainedBlocks_MatlabSparse.m b/+nla/+helpers/+stdError/UnconstrainedBlocks_MatlabSparse.m deleted file mode 100755 index 483d2936..00000000 --- a/+nla/+helpers/+stdError/UnconstrainedBlocks_MatlabSparse.m +++ /dev/null @@ -1,83 +0,0 @@ -classdef UnconstrainedBlocks_MatlabSparse < nla.helpers.stdError.UnconstrainedBlocks - - methods - - function stdErr = calculate(obj, sweStdErrInput) - %Computes Standard Error assuming unconstrained blocks - %Uses standard matrix multiplication to compute standard error, - %since it does not make assumption that V is sparse. - - - %rename variables for readability - pinvDesignMtx = sweStdErrInput.pinvDesignMtx; - residual = sweStdErrInput.residual; - groupIds = sweStdErrInput.scanMetadata.groupId; - unqGrps = unique(groupIds); - numGrps = length(unqGrps); - - obj.throwErrorIfVEntirelyFull(unqGrps); - - [numCovariates, numObs] = size(pinvDesignMtx); - [~,numFcEdges] = size(residual); - - %Preallocate size of stdErr output - stdErr = zeros(numCovariates, numFcEdges); - - %Reorder data so that grouped observations are together - [pinvDesignMtx_grp, residual_grp, groupIds_grp] = obj.reorderDataByGroup(pinvDesignMtx, residual, groupIds); - - %Determine which observations fit in which group once, outside - %of loop for fc edges - inGroupFlags = zeros(numObs, numGrps); - for grpIdx = 1:numGrps - thisGrpId = unqGrps(grpIdx); - inGroupFlags(:,grpIdx) = (groupIds_grp == thisGrpId); - end - - %To prepare sparse matrix, find number of entries for each group - [subjPerGrp,~] = groupcounts(groupIds_grp); - totalNonzerosInV = sum(subjPerGrp.*subjPerGrp); - thisRowInds = zeros(totalNonzerosInV,1); - thisColInds = zeros(totalNonzerosInV,1); - - %precompute row and col inds of nonzero idxs in V - idxOfNonzero = 1; - for grpIdx = 1:numGrps - thisGrpInds = find(inGroupFlags(:,grpIdx)); - - for sparseRow = 1:length(thisGrpInds) - for sparseCol = 1:length(thisGrpInds) - thisRowInds(idxOfNonzero) = thisGrpInds(sparseRow); - thisColInds(idxOfNonzero) = thisGrpInds(sparseCol); - idxOfNonzero = idxOfNonzero+1; - end - end - - end - - - %For each fc edge, build the V and multiply it by - %pinvDesignMatrix on both sides - pinvDesignMtx_grp_tpose = pinvDesignMtx_grp'; - for fcIdx = 1:numFcEdges - - - thisSparseVal = residual_grp(thisRowInds, fcIdx).*residual_grp(thisColInds, fcIdx); - - sparseV = sparse(thisRowInds, thisColInds, thisSparseVal); - - thisBetaCovar = pinvDesignMtx_grp * ... - sparseV * ... - pinvDesignMtx_grp_tpose; - - stdErr(:,fcIdx) = sqrt(diag(thisBetaCovar)); - - end - - - end - - end - - -end \ No newline at end of file diff --git a/+nla/+helpers/+stdError/UnconstrainedBlocks_Sparse.m b/+nla/+helpers/+stdError/UnconstrainedBlocks_Sparse.m deleted file mode 100755 index 477cb756..00000000 --- a/+nla/+helpers/+stdError/UnconstrainedBlocks_Sparse.m +++ /dev/null @@ -1,212 +0,0 @@ -classdef UnconstrainedBlocks_Sparse < nla.helpers.stdError.UnconstrainedBlocks - - properties - APPROX_MAX_MEMORY_GB = 20; %Limit how much data this algorithm should have in memory at a time - end - - methods - - function stdErr = calculate(obj, sweStdErrInput) - %Computes Standard Error, but uses assumption that V matrix is - %sparse to speed up calculation. If V is not sparse, will - %probably take longer than normal matrix multiplication - %V matrix will be calculated in blocks, where each block along - %the diagonal corresponds to a group. Each block will be the - %residual error of all observations in the group multiplied by - %the transpose of the error, ie residual * residual'; - % - %Efficient algo adapted from - %https://www.mathworks.com/matlabcentral/answers/87629-efficiently-multiplying-diagonal-and-general-matrices - %(And in case that page goes away, copying text in file - %/data/wheelock/data1/people/ecka/fastDiagMatrixMultAlgo.txt) - - %Wrinkle here is that when off diagonal elements of V can be - %non-zero, memory considerations must be applied. - - - %rename variables for readability - pinvDesignMtx = sweStdErrInput.pinvDesignMtx; - residual = sweStdErrInput.residual; - groupIds = sweStdErrInput.scanMetadata.groupId; - unqGrps = unique(groupIds); - - obj.throwErrorIfVEntirelyFull(unqGrps); - - - [numCovariates, ~] = size(pinvDesignMtx); - [numObs, numFcEdges] = size(residual); - - %determine which elements of V will be non-zero based on - %grouping. - %First reorganize data so that observations within groups are - %next to eachother - [pinvDesignMtx_grp, residual_grp, groupIds_grp] = obj.reorderDataByGroup(pinvDesignMtx, residual, groupIds); - - - %Find row and column locations where non-zero values in V will - %be - [rowIdxs, colIdxs] = obj.getCoordsOfNonzerosInV(groupIds); - numNonzeroElems = length(rowIdxs); - - %Split data into blocks so that memory used by this algorithm - %will attempt to be kept below APPROX_MAX_MEMORY_GB - edgesPerMemBlock = obj.calcNumEdgesPerMemoryBlock(numNonzeroElems, numCovariates, obj.APPROX_MAX_MEMORY_GB); - - if edgesPerMemBlock < 1 - error(['%s: Number of nonzero V elements will exceed memory ceiling allocated ',... - 'to this class (APPROX_MAX_MEMORY_GB property currently set to %i GB'],class(obj), obj.APPROX_MAX_MEMORY_GB) - end - - %Make precomputed matrix that will multiply all - invDesMtxSelfMultPreCompute = zeros(numCovariates,numCovariates,numNonzeroElems); - - for i = 1:numNonzeroElems - invDesMtxSelfMultPreCompute(:,:,i) = pinvDesignMtx_grp(:,colIdxs(i)) * pinvDesignMtx_grp(:,rowIdxs(i))'; - end - invDesMtxSelfMultPreCompute = reshape(invDesMtxSelfMultPreCompute,numCovariates^2, numNonzeroElems); - - %Use pregenerated matrix to compute covariance of our estimates - %of the regressors. - - - %Precompute which elems of flattened beta covariance matrix - %correspond to diagonal elements if matrix were square - diagElemIdxsInFlatArr = 1:(numCovariates+1):numCovariates^2; - - %Compute stdError one memory block at a time - flatVTheseFcEdges = zeros(numNonzeroElems, edgesPerMemBlock); - stdErr = zeros(numCovariates, numFcEdges); - fcEdgeBlockStart = 1; - - while fcEdgeBlockStart <= numFcEdges - fcEdgeBlockEnd = min(fcEdgeBlockStart + edgesPerMemBlock - 1, numFcEdges); - - %If this is the last block, it may not have as many fcEdges - %as the other blocks and therefore won't fill up the - %preallocated flatVTheseFcEdges matrix. - %Clear the old one and create a new, smaller one. - %This is because I have not found a way to efficiently - %subindex a large matrix in MATLAB. See this MATLAB forum - %post for the weirdness: - % https://www.mathworks.com/matlabcentral/answers/54522-why-is-indexing-vectors-matrices-in-matlab-very-inefficient - %As part of their test, they showed for a large vector a, - %that "a(1:end) = 1" is 5 times slower than "a(:) = 1" - %Best actual speed solution is probably to incorporate - %these algos into a C MEX file. ADE 20220526 - fcEdgesThisBlock = fcEdgeBlockEnd - fcEdgeBlockStart + 1; - if fcEdgesThisBlock < edgesPerMemBlock - clear flatVTheseFcEdges - flatVTheseFcEdges = zeros(numNonzeroElems, fcEdgesThisBlock); - end - - - flatVTheseFcEdges = obj.makeNonZeroFlatV(... - residual_grp(:,fcEdgeBlockStart:fcEdgeBlockEnd),... - groupIds_grp, flatVTheseFcEdges ); - - betaCovarianceFlat = (invDesMtxSelfMultPreCompute * flatVTheseFcEdges); - - - stdErr(:,fcEdgeBlockStart:fcEdgeBlockEnd) = sqrt(betaCovarianceFlat(diagElemIdxsInFlatArr,:)); - - fcEdgeBlockStart = fcEdgeBlockEnd + 1; - end - - end - - end - - methods (Access = private) - - - function [rowIdxs, colIdxs] = getCoordsOfNonzerosInV(obj, groupIdVec) - %Find the row and column indices that will have non-zero values - %if we were to build the whole V matrix for correlated blocks - %based on the groupIds of the elements - - numObs = length(groupIdVec); - vTemplate = zeros(numObs, numObs); - rowIndexMtx = (1:numObs)' * ones(1,numObs); - colIndexMtx = ones(numObs,1) * (1:numObs); - - unqGrps = unique(groupIdVec); - countInGrps = histcounts(groupIdVec,[unqGrps;Inf]); - - tmpRowIdx = 1; - for grpIdx = 1:length(unqGrps) - startIdx = tmpRowIdx; - endIdx = tmpRowIdx + countInGrps(grpIdx)-1; - - vTemplate([startIdx:endIdx],[startIdx:endIdx]) = 1; - tmpRowIdx = endIdx + 1; - end - - vTemplateFlat = reshape(vTemplate,numel(vTemplate),1); - rowMtxFlat = reshape(rowIndexMtx,numel(vTemplate),1); - colMtxFlat = reshape(colIndexMtx,numel(vTemplate),1); - - rowIdxs = rowMtxFlat(vTemplateFlat==1); - colIdxs = colMtxFlat(vTemplateFlat==1); - - end - - function flatV = makeNonZeroFlatV(obj, residual, groupIds, flatV) - %assumes residuals have already been sorted by groupId - %computes the nonzero elements of V given correlated errors - %within groups - % - %Accepts a previously allocated block of flatV as input and - %overwrites values during algo for speed / memory concerns - % - %If flatV is larger than needed for residual, write only to the - %first columns needed - - - unqGrps = unique(groupIds); - countInGrps = histcounts(groupIds,[unqGrps;Inf]); - - [~, numFcEdgesThisBlock] = size(residual); - [~, numFcEdgesInFlatV] = size(flatV); - - if numFcEdgesInFlatV > numFcEdgesThisBlock - flatV(:,(numFcEdgesThisBlock+1):end) = []; - end - - rowIdx = 1; - - for grpIdx = 1:length(unqGrps) - - thisGrpId = unqGrps(grpIdx); - numObsThisGrp = countInGrps(grpIdx); - thisGrpFlag = groupIds == thisGrpId; - residThisGrp = residual(thisGrpFlag,:); - - - for i = 1:numObsThisGrp - for j = 1:numObsThisGrp - flatVAllEdgesThisRow = residThisGrp(i,:) .* residThisGrp(j,:); - flatV(rowIdx,:) = flatVAllEdgesThisRow; - rowIdx = rowIdx + 1; - end - end - - end - - end - - function edgesPerMemBlock = calcNumEdgesPerMemoryBlock(obj, numNonzeroElemsInV, numCovariates, memBlockSize_GB) - - GBPerMtxElem = 8 / (1024*1024*1024); - - memPrecomputedMtx_GB = numNonzeroElemsInV * (numCovariates^2) * GBPerMtxElem; - memForFcEdgeBlocks_GB = memBlockSize_GB - memPrecomputedMtx_GB; - - memPerEdge_GB = numNonzeroElemsInV * GBPerMtxElem; - - edgesPerMemBlock = floor(memForFcEdgeBlocks_GB / memPerEdge_GB); - - end - - end - -end \ No newline at end of file diff --git a/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp b/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp index 9fdb4e49..420fd283 100644 Binary files a/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp and b/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp differ diff --git a/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m b/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m index dc5363df..bc2c8c0c 100644 --- a/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m +++ b/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m @@ -34,6 +34,8 @@ ConvergencePlotColorDropDownLabel matlab.ui.control.Label ConvergencePlotColorDropDown matlab.ui.control.DropDown ApplyButton matlab.ui.control.Button + PlotValueDropDownLabel matlab.ui.control.Label + PlotValueDropDown matlab.ui.control.DropDown Panel_2 matlab.ui.container.Panel end @@ -102,18 +104,7 @@ function getPlotTitle(app) import nla.net.result.NetworkTestResult if ~isequal(app.matrix_plot, false) - app.settings = struct(); - app.settings.upperLimit = app.UpperLimitEditField.Value; - app.settings.lowerLimit = app.LowerLimitEditField.Value; - app.settings.prevPlotScale = app.matrix_plot.plot_scale; - app.settings.plotScale = app.PlotScaleDropDown.Value; - app.settings.pValueThreshold = app.pvalueThresholdEditField.Value; - app.settings.cohensD = app.CohensDThresholdCheckBox.Value; - app.settings.cohensDValue = app.CohensDThresholdEditField.Value; - app.matrix_plot.display_legend.Visible = app.LegendVisibleDropDown.Value; - app.matrix_plot.plot_title.String = {}; - app.network_test_options.prob_max = app.pvalueThresholdEditField.Value; - app.settings.ranking = app.RankingDropDown.Value; + app.save_plot_settings(); end if isfield(app.settings, "ranking") @@ -146,30 +137,24 @@ function getPlotTitle(app) end probability = NetworkTestResult().getPValueNames(app.test_method, app.network_test_result.test_name); - - app.network_test_options.show_ROI_centroids = app.ROIcentroidsonbrainplotsCheckBox.Value; - - app.parameters = nla.net.result.NetworkResultPlotParameter(app.network_test_result, app.edge_test_options.net_atlas,... - app.network_test_options); - + app.parameters = nla.net.result.NetworkResultPlotParameter(app.network_test_result,... + app.edge_test_options.net_atlas, app.network_test_options); probability_parameters = app.parameters.plotProbabilityParameters(app.edge_test_options, app.edge_test_result,... - app.test_method, probability, sprintf(app.title), mcc, app.createSignificanceFilter(),... - app.RankingDropDown.Value); + app.test_method, probability, sprintf(app.title), mcc, app.createEffectSizeFilter(),... + app.RankingDropDown.Value, app.PlotValueDropDown.Value); % if ~isequal(app.UpperLimitEditField.Value, 0.3) && ~isequal(app.LowerLimitEditField.Value, 0.3) - probability_parameters.p_value_plot_max = app.pvalueThresholdEditField.Value; + probability_parameters.plot_max = app.pvalueThresholdEditField.Value; % end + app.network_test_options.show_ROI_centroids = app.ROIcentroidsonbrainplotsCheckBox.Value; + plotter = nla.net.result.plot.PermutationTestPlotter(app.edge_test_options.net_atlas); - [width, height, app.matrix_plot] = plotter.plotProbability(app.Panel_2, probability_parameters, nla.inputField.LABEL_GAP, -50); + [width, height, app.matrix_plot] = plotter.plotProbability(app.Panel_2, probability_parameters,... + nla.inputField.LABEL_GAP, -50); + if ~isequal(app.settings, false) - app.PlotScaleDropDown.Value = app.settings.plotScale; - app.UpperLimitEditField.Value = app.settings.upperLimit; - app.LowerLimitEditField.Value = app.settings.lowerLimit; - app.pvalueThresholdEditField.Value = app.settings.pValueThreshold; - app.CohensDThresholdCheckBox.Value = app.settings.cohensD; - app.CohensDThresholdEditField.Value = app.settings.cohensDValue; - app.matrix_plot.display_legend.Visible = app.LegendVisibleDropDown.Value; + app.load_plot_settings(); app.applyScaleChange(); end end @@ -177,7 +162,53 @@ function getPlotTitle(app) methods (Access = private) - function cohens_d_filter = createSignificanceFilter(app) + function save_plot_settings(app) + if isequal(app.settings, false) + app.settings = struct(); + app.settings.statValueLimits = struct(); + app.settings.pValueLimits = struct(); + end + if ~isfield(app.settings, "plotValue") + app.settings.pValueLimits.upperLimit = app.UpperLimitEditField.Value; + app.settings.pValueLimits.lowerLimit = app.LowerLimitEditField.Value; + app.settings.statValueLimits.upperLimit = max(app.network_test_result.no_permutations.(app.network_test_result.ranking_statistic).v); + app.settings.statValueLimits.lowerLimit = min(app.network_test_result.no_permutations.(app.network_test_result.ranking_statistic).v); + elseif isequal(app.settings.plotValue, "nla.gfx.PlotValue.PVALUE") + app.settings.pValueLimits.upperLimit = app.UpperLimitEditField.Value; + app.settings.pValueLimits.lowerLimit = app.LowerLimitEditField.Value; + else + app.settings.statValueLimits.upperLimit = app.UpperLimitEditField.Value; + app.settings.statValueLimits.lowerLimit = app.LowerLimitEditField.Value; + end + app.settings.plotValue = app.PlotValueDropDown.Value; + app.settings.prevPlotScale = app.matrix_plot.plot_scale; + app.settings.plotScale = app.PlotScaleDropDown.Value; + app.settings.pValueThreshold = app.pvalueThresholdEditField.Value; + app.settings.cohensD = app.CohensDThresholdCheckBox.Value; + app.settings.cohensDValue = app.CohensDThresholdEditField.Value; + app.matrix_plot.display_legend.Visible = app.LegendVisibleDropDown.Value; + app.matrix_plot.plot_title.String = {}; + app.network_test_options.prob_max = app.pvalueThresholdEditField.Value; + app.settings.ranking = app.RankingDropDown.Value; + end + + function load_plot_settings(app) + app.PlotScaleDropDown.Value = app.settings.plotScale; + if isequal(app.PlotValueDropDown.Value, "nla.gfx.PlotValue.PVALUE") + app.UpperLimitEditField.Value = app.settings.pValueLimits.upperLimit; + app.LowerLimitEditField.Value = app.settings.pValueLimits.lowerLimit; + else + app.UpperLimitEditField.Value = app.settings.statValueLimits.upperLimit; + app.LowerLimitEditField.Value = app.settings.statValueLimits.lowerLimit; + end + app.PlotValueDropDown.Value = app.settings.plotValue; + app.pvalueThresholdEditField.Value = app.settings.pValueThreshold; + app.CohensDThresholdCheckBox.Value = app.settings.cohensD; + app.CohensDThresholdEditField.Value = app.settings.cohensDValue; + app.matrix_plot.display_legend.Visible = app.LegendVisibleDropDown.Value; + end + + function cohens_d_filter = createEffectSizeFilter(app) %REMOVE COHENS D FILTERING UNTIL WE DETERMINE CORRECT CALCULATION FOR IT - ADE 2025MAR24 num_nets = app.edge_test_options.net_atlas.numNets; cohens_d_filter = nla.TriMatrix(num_nets, "logical", nla.TriMatrixDiag.KEEP_DIAGONAL); @@ -266,7 +297,7 @@ function drawChords(app, event) probability = NetworkTestResult().getPValueNames(app.test_method, app.network_test_result.test_name); p_value = strcat("uncorrected_", probability); probability_parameters = app.parameters.plotProbabilityParameters(app.edge_test_options, app.edge_test_result,... - app.test_method, p_value, sprintf(app.title), app.MultipleComparisonCorrectionDropDown.Value, app.createSignificanceFilter(),... + app.test_method, p_value, sprintf(app.title), app.MultipleComparisonCorrectionDropDown.Value, app.createEffectSizeFilter(),... app.RankingDropDown.Value); chord_plotter = nla.net.result.chord.ChordPlotter(app.edge_test_options.net_atlas, app.edge_test_result); @@ -282,9 +313,17 @@ function drawChords(app, event) function PlotScaleValueChanged(app, event) if isequal(app.settings, false) app.settings = struct(); + app.settings.pValueLimits = struct(); + app.settings.statValueLimits = struct(); end - app.settings.upperLimit = app.UpperLimitEditField.Value; - app.settings.lowerLimit = app.LowerLimitEditField.Value; + if isequal(app.PlotValueDropDown.Value, "nla.gfx.PlotValue.PVALUE") + app.settings.pValueLimits.upperLimit = app.UpperLimitEditField.Value; + app.settings.pValueLimits.lowerLimit = app.LowerLimitEditField.Value; + else + app.settings.statValueLimits.upperLimit = app.UpperLimitEditField.Value; + app.settings.statValueLimits.lowerLimit = app.LowerLimitEditField.Value; + end + app.settings.plotValue = app.PlotValueDropDown.Value; app.settings.plotScale = app.PlotScaleDropDown.Value; app.settings.pValueThreshold = app.pvalueThresholdEditField.Value; app.settings.cohensD = app.CohensDThresholdCheckBox.Value; @@ -448,36 +487,36 @@ function createComponents(app) % Create UpperLimitEditFieldLabel app.UpperLimitEditFieldLabel = uilabel(app.Panel); app.UpperLimitEditFieldLabel.HorizontalAlignment = 'right'; - app.UpperLimitEditFieldLabel.Position = [46 267 70 22]; + app.UpperLimitEditFieldLabel.Position = [46 236 70 22]; app.UpperLimitEditFieldLabel.Text = 'Upper Limit'; % Create UpperLimitEditField app.UpperLimitEditField = uieditfield(app.Panel, 'numeric'); app.UpperLimitEditField.ValueChangedFcn = createCallbackFcn(app, @PlotScaleValueChanged, true); - app.UpperLimitEditField.Position = [126 267 52 22]; + app.UpperLimitEditField.Position = [126 236 52 22]; app.UpperLimitEditField.Value = 0.05; % Create LowerLimitEditFieldLabel app.LowerLimitEditFieldLabel = uilabel(app.Panel); app.LowerLimitEditFieldLabel.HorizontalAlignment = 'right'; - app.LowerLimitEditFieldLabel.Position = [259 267 70 22]; + app.LowerLimitEditFieldLabel.Position = [46 206 70 22]; app.LowerLimitEditFieldLabel.Text = 'Lower Limit'; % Create LowerLimitEditField app.LowerLimitEditField = uieditfield(app.Panel, 'numeric'); app.LowerLimitEditField.ValueChangedFcn = createCallbackFcn(app, @PlotScaleValueChanged, true); - app.LowerLimitEditField.Position = [339 267 52 22]; + app.LowerLimitEditField.Position = [126 206 52 22]; % Create pvalueThresholdEditFieldLabel app.pvalueThresholdEditFieldLabel = uilabel(app.Panel); app.pvalueThresholdEditFieldLabel.HorizontalAlignment = 'right'; - app.pvalueThresholdEditFieldLabel.Position = [14 237 102 22]; + app.pvalueThresholdEditFieldLabel.Position = [227 267 102 22]; app.pvalueThresholdEditFieldLabel.Text = 'p-value Threshold'; % Create pvalueThresholdEditField app.pvalueThresholdEditField = uieditfield(app.Panel, 'numeric'); app.pvalueThresholdEditField.ValueChangedFcn = createCallbackFcn(app, @pvalueThresholdEditFieldValueChanged, true); - app.pvalueThresholdEditField.Position = [126 237 52 22]; + app.pvalueThresholdEditField.Position = [339 267 52 22]; app.pvalueThresholdEditField.Value = 0.05; % Create CohensDThresholdEditFieldLabel @@ -542,7 +581,7 @@ function createComponents(app) app.ROIcentroidsonbrainplotsCheckBox = uicheckbox(app.Panel); app.ROIcentroidsonbrainplotsCheckBox.ValueChangedFcn = createCallbackFcn(app, @ROIcentroidsonbrainplotsCheckBoxValueChanged, true); app.ROIcentroidsonbrainplotsCheckBox.Text = 'ROI centroids on brain plots'; - app.ROIcentroidsonbrainplotsCheckBox.Position = [7 207 171 22]; + app.ROIcentroidsonbrainplotsCheckBox.Position = [221 27 171 22]; % Create ViewChordPlotsButton app.ViewChordPlotsButton = uibutton(app.Panel, 'push'); @@ -593,6 +632,19 @@ function createComponents(app) app.ApplyButton.Position = [21 27 100 22]; app.ApplyButton.Text = 'Apply'; + % Create PlotValueDropDownLabel + app.PlotValueDropDownLabel = uilabel(app.Panel); + app.PlotValueDropDownLabel.HorizontalAlignment = 'right'; + app.PlotValueDropDownLabel.Position = [-2 270 65 22]; + app.PlotValueDropDownLabel.Text = 'Plot Value'; + + % Create PlotValueDropDown + app.PlotValueDropDown = uidropdown(app.Panel); + app.PlotValueDropDown.Items = {'p-value', 'Test Statistic'}; + app.PlotValueDropDown.ItemsData = {'nla.gfx.PlotValue.PVALUE', 'nla.gfx.PlotValue.STATISTIC', ''}; + app.PlotValueDropDown.Position = [80 270 98 20]; + app.PlotValueDropDown.Value = 'nla.gfx.PlotValue.PVALUE'; + % Create Panel_2 app.Panel_2 = uipanel(app.UIFigure); app.Panel_2.AutoResizeChildren = 'off'; diff --git a/+nla/+net/+result/NetworkResultPlotParameter.m b/+nla/+net/+result/NetworkResultPlotParameter.m index 401a067f..9a7bc49b 100644 --- a/+nla/+net/+result/NetworkResultPlotParameter.m +++ b/+nla/+net/+result/NetworkResultPlotParameter.m @@ -27,12 +27,13 @@ end function result = plotProbabilityParameters(obj, edge_test_options, edge_test_result, test_method, plot_statistic,... - plot_title, fdr_correction, effect_size_filter, ranking_method) + plot_title, fdr_correction, effect_size_filter, ranking_method, plot_value) % plot_title - this will be a string % plot_statistic - this is the stat that will be plotted - % significance filter - this will be a boolean or some sort of object (like Cohen's D > D-value) + % effect_size_filter - this will be a boolean or some sort of object (like Cohen's D > D-value) % fdr_correction - a struct of fdr_correction (found in nla.net.mcc) % test_method - 'no permutations', 'within network pair', 'full connectome' + % plot_value - nla.gfx.PlotValue.PVALUE or nla.gfx.PlotValue.STATISTIC import nla.TriMatrix nla.TriMatrixDiag % We're going to use a default filter here @@ -76,33 +77,49 @@ % scale values very slightly for display so numbers just below % the threshold don't show up white but marked significant statistic_input_scaled = TriMatrix(obj.number_of_networks, "double", TriMatrixDiag.KEEP_DIAGONAL); - statistic_input_scaled.v = statistic_input.v .* (obj.default_discrete_colors / (obj.default_discrete_colors + 1)); - - % default values for plotting - statistic_plot_matrix = statistic_input_scaled; p_value_plot_max = p_value_max; - significance_type = "nla.gfx.SigType.DECREASING"; - % determine colormap and operate on values if it's -log10 - switch obj.updated_test_options.prob_plot_method - case "LOG" % FUCK Matlab and their enums - color_map = nla.net.result.NetworkResultPlotParameter.getLogColormap(obj.default_discrete_colors,... - statistic_input, p_value_max); - % Here we take a -log10 and change the maximum value to show on the plot - case "NEGATIVE_LOG_10" - color_map = parula(obj.default_discrete_colors); - - statistic_matrix = nla.TriMatrix(obj.number_of_networks, "double", nla.TriMatrixDiag.KEEP_DIAGONAL); - statistic_matrix.v = -log10(statistic_input.v); - statistic_plot_matrix = statistic_matrix; - if strcmp(test_method, "full_connectome") || strcmp(test_method, "within_network_pair") - p_value_plot_max = 2; - else - p_value_plot_max = 40; - end - significance_type = "nla.gfx.SigType.INCREASING"; - otherwise - color_map = nla.net.result.NetworkResultPlotParameter.getColormap(obj.default_discrete_colors,... - p_value_max); + + if isequal(plot_value, "nla.gfx.PlotValue.PVALUE") + statistic_input_scaled.v = statistic_input.v .* (obj.default_discrete_colors / (obj.default_discrete_colors + 1)); + + % default values for plotting + statistic_plot_matrix = statistic_input_scaled; + significance_type = "nla.gfx.SigType.DECREASING"; + % determine colormap and operate on values if it's -log10 + switch obj.updated_test_options.prob_plot_method + case "LOG" % FUCK Matlab and their enums + color_map = nla.net.result.NetworkResultPlotParameter.getLogColormap(obj.default_discrete_colors,... + statistic_input, p_value_max); + % Here we take a -log10 and change the maximum value to show on the plot + case "NEGATIVE_LOG_10" + color_map = parula(obj.default_discrete_colors); + + statistic_matrix = nla.TriMatrix(obj.number_of_networks, "double", nla.TriMatrixDiag.KEEP_DIAGONAL); + statistic_matrix.v = -log10(statistic_input.v); + statistic_plot_matrix = statistic_matrix; + if strcmp(test_method, "full_connectome") || strcmp(test_method, "within_network_pair") + p_value_plot_max = 2; + else + p_value_plot_max = 40; + end + significance_type = "nla.gfx.SigType.INCREASING"; + otherwise + color_map = nla.net.result.NetworkResultPlotParameter.getColormap(obj.default_discrete_colors,... + p_value_max); + end + else + ranking_statistic = obj.network_test_results.ranking_statistic; + statistic_input_scaled.v = obj.network_test_results.no_permutations.(ranking_statistic).v .* (... + obj.default_discrete_colors / (obj.default_discrete_colors + 1)... + ); + + statistic_plot_matrix = statistic_input_scaled; + % switch obj.updated_test_options.prob_plot_method + % case "LOG" + % color_map + significance_type = "nla.gfx.SigType.DECREASING"; + color_map = nla.net.result.NetworkResultPlotParameter.getColormap(obj.default_discrete_colors,... + p_value_plot_max); end % callback function for brain image. diff --git a/+nla/+net/unittests/NetworkResultPlotParameterTestCase.m b/+nla/+net/unittests/NetworkResultPlotParameterTestCase.m index 2842a68f..09f8221a 100644 --- a/+nla/+net/unittests/NetworkResultPlotParameterTestCase.m +++ b/+nla/+net/unittests/NetworkResultPlotParameterTestCase.m @@ -92,7 +92,7 @@ function plotProbabilityParametersDefaultPlottingTest(testCase) probability_parameters = plot_parameters.plotProbabilityParameters(testCase.edge_test_options,... testCase.edge_test_result, "full_connectome", probability, 'Title', nla.net.mcc.Bonferroni(),... - false, nla.RankingMethod.UNCORRECTED); + false, nla.RankingMethod.UNCORRECTED, "nla.gfx.PlotValue.PVALUE"); expected_p_value_max = testCase.network_test_options.prob_max / testCase.network_atlas.numNetPairs(); expected_plot = nla.TriMatrix(plot_parameters.number_of_networks, "double", nla.TriMatrixDiag.KEEP_DIAGONAL); @@ -120,7 +120,7 @@ function plotProbabilityParametersLogPlottingTest(testCase) probability_parameters = plot_parameters.plotProbabilityParameters(testCase.edge_test_options,... testCase.edge_test_result, "full_connectome", probability, 'Title', nla.net.mcc.Bonferroni(),... - false, nla.RankingMethod.UNCORRECTED); + false, nla.RankingMethod.UNCORRECTED, "nla.gfx.PlotValue.PVALUE"); expected_p_value_max = testCase.network_test_options.prob_max / testCase.network_atlas.numNetPairs(); expected_plot = nla.TriMatrix(plot_parameters.number_of_networks, "double", nla.TriMatrixDiag.KEEP_DIAGONAL); @@ -158,7 +158,7 @@ function plotProbabilityParametersNegLogTest(testCase) probability_parameters = plot_parameters.plotProbabilityParameters(testCase.edge_test_options,... testCase.edge_test_result, "full_connectome", probability, 'Title', nla.net.mcc.Bonferroni(),... - false, nla.RankingMethod.UNCORRECTED); + false, nla.RankingMethod.UNCORRECTED, "nla.gfx.PlotValue.PVALUE"); expected_p_value_max = 2; expected_plot = nla.TriMatrix(plot_parameters.number_of_networks, "double", nla.TriMatrixDiag.KEEP_DIAGONAL); diff --git a/+nla/unittests/TestPoolTest.m b/+nla/unittests/TestPoolTest.m index 1f2a9814..e228a217 100644 --- a/+nla/unittests/TestPoolTest.m +++ b/+nla/unittests/TestPoolTest.m @@ -133,18 +133,6 @@ function hyperGeometricTest(testCase) end end - function kolmogorovSmirnovTest(testCase) - edge_result = load(strcat(testCase.root_path, fullfile("+nla", "unittests", "pearson_result.mat"))); - testCase.tests.net_tests = {nla.net.test.KolmogorovSmirnovTest()}; - network_result = testCase.tests.runNetTestsPerm(testCase.network_test_options, testCase.edge_test_options.net_atlas, edge_result.edge_result); - - expected_result = load(strcat(testCase.root_path, fullfile("+nla", "unittests", "kolmogorov_smirnov_result.mat"))); - property_names = properties(network_result{1}); - for prop_name = property_names - testCase.verifyEqual(expected_result.network_result.(prop_name{1}), network_result{1}.(prop_name{1})); - end - end - function studentTTest(testCase) edge_result = load(strcat(testCase.root_path, fullfile("+nla", "unittests", "pearson_result.mat"))); testCase.tests.net_tests = {nla.net.test.StudentTTest()};