text:
function net = insertLossLayers(net, varargin)

import dagnn.*
opts.numClasses = 101;
opts = vl_argparse(opts, varargin);

for p = 1:numel(net.params)
  name = net.params(p).name ;
  sz = size(net.params(p).value);
  if numel(sz) == 4
    if sz(4) == opts.numClasses
      lossAlreadyExists = false;
        for l=1:numel(net.layers)
          if ~isempty(net.layers(l).params)
            if any(strcmp(name, net.layers(l).params(:))), break; end
          end
        end
        pred_layer = l;
        for l = 1:numel(net.layers)
          if any(strcmp(net.layers(pred_layer).outputs, net.layers(l).inputs))  && ...
              isa(net.layers(l).block, 'dagnn.Loss') && isempty(strfind(net.layers(l).name, 'err'))
              lossAlreadyExists = true;
              break;
          end
        end
        if lossAlreadyExists, continue; end;
        name = sprintf('loss_%s',net.layers(pred_layer).name);
        net.addLayer(name, ...
          dagnn.Loss( 'loss', 'softmaxlog'), ...
             [net.layers(pred_layer).outputs{:} {'label'}], sprintf('objective_%s',net.layers(pred_layer).name)) ;
    end
  end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function deploy_each_dagnet(epoch, datapath)

for i = 1:epoch
    fileName = strcat('net-epoch-',num2str(i),'.mat');
    fileName = fullfile(datapath,fileName);
    [net, info] = loadState_deploy(fileName);
    net = net_deploy(net,1) ;
    deployname = strcat('net-deployed-',num2str(i),'.mat');
    save(fullfile(datapath, 'fine-tuned-model', deployname), 'net', 'info', '-v7.3');
end

repo: DBCNN


text:
%%% This program provides a numerical verification that equations (15.2)
%%% and (15.6) coincide in a 2x2 case. The program generates a random 2x2
%%% GOE matrix H and uses it to compute the numerical value of Z(x)
%%% according to the two formulas. You will be asked to choose the real and
%%% imaginary parts of the number x (imaginary part has to be positive).

clear all
close all

%%% Reads the real part from the Command Window
prompt = '\n Choose real part of argument x: ';
re = input(prompt);

%%% Reads the imaginary part from the Command Window
prompt = '\n Choose imaginary part of argument x: ';
im = input(prompt);

if im <= 0
    sprintf('ERROR: Imaginary part has to be positive')
    return;
end

x = re - i*im;

%%% Generating 2x2 GOE matrix
H = randn(2)/sqrt(2);
H = (H + H')/2;

%%% Eigenvalues of matrix H
E = eig(H);

%%% Definition of integrand function for Z
f = @(y1,y2) exp(-i*((x-H(1,1)).*y1.^2 + 2.*H(1,2).*y1.*y2 + (x-H(2,2)).*y2.^2)/2);

%%% Z function
Z = integral2(f,-Inf,Inf,-Inf,Inf);

%%% Explicit form of Z function
Z_exp = 2*pi*exp(-(log(E(1)-x) + log(E(2)-x))/2 + i*pi/2);

sprintf('Value of Z function computed as integral: %6.4f%+6.4fi',real(Z),imag(Z))
sprintf('Value of Z function computed explicitly via eigenvalues of H: %6.4f%+6.4fi',real(Z_exp),imag(Z_exp))

author: RMT-TheoryAndPractice

repo: RMT


text:
function []=anim8_DIC3D_reconstructedPairs_faceMeasure(DIC3DAllPairsResults,faceMeasureString,varargin)
%% function for plotting 3D-DIC results of face measures in STEP3.
% plotting 3D surfaces from camera pairs, animation changing
% with time, and the faces colored according to faceMeasureString
% this function is called in plotMultiDICPairResults
%
% Options:
% anim8_DIC_3D_pairs_faceMeasure(DIC3DAllPairsResults,faceMeasureString)
% anim8_DIC_3D_pairs_faceMeasure(DIC3DAllPairsResults,faceMeasureString,optStruct)
%
% Inputs:
% * DIC3DAllPairsResults
% * faceMeasureString: can be any of the following:
%   'dispMgn','dispX','dispY','dispZ','FaceColors','FaceIsoInd','pairInd'
% * optStruct: optional structure for plotting options which may include any of the following fields:
%   - smoothLogic: logical variable for smoothing (true)/not smoothing (false) the face measure
%   - FaceAlpha: transparacy of the faces (scalar between 0 and 1, where zero is transparent and 1 is opaque)
%   - colorBarLimits: a 2x1 scalar vector for the colobar limits. if not set, it's automatic
%   - dataLimits: a 2x1 scalar vector for the data limits of the face measure. if a face measure is outside these limits, it is set to NaN. if not set no face is set to NaN
%   - colorMap
%   - zDirection: 1 for z up and -1 for z down
%   - lineColor: line color for the mesh. can be for example 'b','k','none',etc...
%   - TitleString=faceMeasureString;

%% Assign plot options
Narg=numel(varargin);
switch Narg
    case 1
        optStruct=varargin{1};
    case 0
        optStruct=struct;
    otherwise
        ('wrong number of input arguments');
end

% complete the struct fields
if ~isfield(optStruct,'smoothLogic')
    optStruct.smoothLogic=0;
end
if ~isfield(optStruct,'dataLimits')
    optStruct.dataLimits=[-inf inf];
end
if ~isfield(optStruct,'zDirection') % 1 or -1
    optStruct.zDirection=1;
end
if ~isfield(optStruct,'lineColor') % 'none' or 'k'
    optStruct.lineColor='none';
end
if ~isfield(optStruct,'TitleString')
    optStruct.TitleString=faceMeasureString;
end
if ~isfield(optStruct,'maxCorrCoeff')
    optStruct.maxCorrCoeff=[];
end

%%
nPairs=numel(DIC3DAllPairsResults);
nFrames=numel(DIC3DAllPairsResults{1}.Points3D);

%% Assign the right face measure into FC
FC=cell(nPairs,nFrames);
switch faceMeasureString

    case {'DispMgn'}
        for ip=1:nPairs
            for it=1:nFrames
                Fnow=DIC3DAllPairsResults{ip}.Faces;
                dispNow=DIC3DAllPairsResults{ip}.Disp.(faceMeasureString){it}; % point measure
                FC{ip,it}=mean(dispNow(Fnow),2); % turn into face measure
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        FCmat = cell2mat(FC);
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[0 prctile(FCmat(:),100)];
        end
        colorBarLogic=1;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='parula';
        end
                if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=0.8;
        end
        legendLogic=0;

    case {'DispX'}
        for ip=1:nPairs
            for it=1:nFrames
                Fnow=DIC3DAllPairsResults{ip}.Faces;
                dispNow=DIC3DAllPairsResults{ip}.Disp.DispVec{it}(:,1); % point measure
                FC{ip,it}=mean(dispNow(Fnow),2); % turn into face measure
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        FCmat = cell2mat(FC);
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[min(FCmat(:)) max(FCmat(:))];
        end
        colorBarLogic=1;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='parula';
        end
                if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=0.8;
        end
        legendLogic=0;

    case {'DispY'}
        for ip=1:nPairs
            for it=1:nFrames
                Fnow=DIC3DAllPairsResults{ip}.Faces;
                dispNow=DIC3DAllPairsResults{ip}.Disp.DispVec{it}(:,2); % point measure
                FC{ip,it}=mean(dispNow(Fnow),2); % turn into face measure
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        FCmat = cell2mat(FC);
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[min(FCmat(:)) max(FCmat(:))];
        end
        colorBarLogic=1;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='parula';
        end
                if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=0.8;
        end
        legendLogic=0;

    case {'DispZ'}
        for ip=1:nPairs
            for it=1:nFrames
                Fnow=DIC3DAllPairsResults{ip}.Faces;
                dispNow=DIC3DAllPairsResults{ip}.Disp.DispVec{it}(:,3); % point measure
                FC{ip,it}=mean(dispNow(Fnow),2); % turn into face measure
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        FCmat = cell2mat(FC);
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[min(FCmat(:)) max(FCmat(:))];
        end
        colorBarLogic=1;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='parula';
        end
                if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=0.8;
        end
        legendLogic=0;

    case {'FaceCorrComb'}
        for ip=1:nPairs
            for it=1:nFrames
                FC{ip,it}=DIC3DAllPairsResults{ip}.(faceMeasureString){it};
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        FCmat = cell2mat(FC);
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[0 prctile(FCmat(:),100)];
        end
        colorBarLogic=1;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='parula';
        end
                if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=0.8;
        end
        legendLogic=0;

    case {'FaceColors'}
        for ip=1:nPairs
            for it=1:nFrames
                FC{ip,it}=DIC3DAllPairsResults{ip}.(faceMeasureString);
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[0 255];
        end
        colorBarLogic=1;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='gray';
        end
        legendLogic=0;
        if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=1;
        end
    case 'pairInd'
        for ip=1:nPairs
            for it=1:nFrames
                FC{ip,it}=ip;
                if ~isempty(optStruct.maxCorrCoeff)
                    corrNow=DIC3DAllPairsResults{ip}.FaceCorrComb{it};
                    FC{ip,it}(corrNow>optStruct.maxCorrCoeff,:)=NaN;
                end
            end
        end
        if ~isfield(optStruct,'colorBarLimits')
            optStruct.colorBarLimits=[1 nPairs];
        end
        colorBarLogic=0;
        if ~isfield(optStruct,'colorMap')
            optStruct.colorMap='gjet';
        end
                if ~isfield(optStruct,'FaceAlpha')
            optStruct.FaceAlpha=0.8;
        end
        legendLogic=1;
    otherwise
        error('unexpected face measure string. plots not created');

end

%% Plot
DIC_3Djoined_results=join3DreconstructedPairs(DIC3DAllPairsResults);
[xl,yl,zl]=axesLimits(DIC_3Djoined_results.Points3D);

animStruct=struct;

hf=cFigure;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

axisGeom;
ax=gca;
ax.CameraUpVector=[0 0 optStruct.zDirection];
Colors=gjet(nPairs);

colormap(optStruct.colorMap);
if colorBarLogic
    cbh=colorbar;
    caxis(optStruct.colorBarLimits);
end

gtitle(optStruct.TitleString,20);
% axis off
% camlight headlight

it=1;
for ip=1:nPairs
    Fnow=DIC3DAllPairsResults{ip}.Faces;
    Pnow=DIC3DAllPairsResults{ip}.Points3D{it};
    CFnow=FC{ip,it};
    if optStruct.smoothLogic
        [CFnow]=triSmoothFaceMeasure(CFnow,Fnow,Pnow,[],[]);
    end
    CFnow(CFnow<optStruct.dataLimits(1))=NaN;
    CFnow(CFnow>optStruct.dataLimits(2))=NaN;
    if strcmp(faceMeasureString,'pairInd')
        hp(ip)=gpatch(Fnow,Pnow,Colors(ip,:),Colors(ip,:),optStruct.FaceAlpha); hold on
    else
        hp(ip)=gpatch(Fnow,Pnow,CFnow,optStruct.lineColor,optStruct.FaceAlpha); hold on
    end
end

if legendLogic
    pairStrings=arrayfun(@num2str,1:nPairs,'unif',0);
    legendStrings=strcat('pair',pairStrings);
    legend(legendStrings);
end

h_ax=gca;
h_ax.XLim = xl; h_ax.YLim = yl; h_ax.ZLim = zl;

animStruct.Time=1:nFrames;
animStruct.Handles=cell(1,nFrames);
animStruct.Props=cell(1,nFrames);
animStruct.Set=cell(1,nFrames);

ic=1;
for it=1:nFrames
    animStruct.Handles{ic}=[];
    animStruct.Props{ic}=cell(1,2*nPairs);
    animStruct.Set{ic}=cell(1,2*nPairs);

    for ip=1:nPairs
        Fnow=DIC3DAllPairsResults{ip}.Faces;
        Pnow=DIC3DAllPairsResults{ip}.Points3D{it};
        CFnow=FC{ip,it};
        if optStruct.smoothLogic
            [CFnow]=triSmoothFaceMeasure(CFnow,Fnow,Pnow,[],[]);
        end
        CFnow(CFnow<optStruct.dataLimits(1))=NaN;
        CFnow(CFnow>optStruct.dataLimits(2))=NaN;
        animStruct.Handles{ic}=[animStruct.Handles{ic} hp(ip) hp(ip)]; %Handles of objects to animate
        animStruct.Props{ic}{2*ip-1}='CData';
        animStruct.Props{ic}{2*ip}='Vertices'; %Properties of objects to animate
        animStruct.Set{ic}{2*ip-1}=CFnow;
        animStruct.Set{ic}{2*ip}=Pnow; %Property values for to set in order to animate
    end

    h_ax.XLim = xl; h_ax.YLim = yl; h_ax.ZLim = zl;

    ic=ic+1;
end
anim8(hf,animStruct);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [] = PSDSInSAR_estimator(Coh,  slcstack, slclist, interfstack, interflist, SHP, reference_ind,  InSAR_path, BroNumthre, Cohthre, Cohthre_slc_filt, InSAR_processor)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This file is part of TomoSAR.
%
%   TomoSAR is distributed in the hope that it will be useful,
%   but without warranty of any kind; without even the implied
%   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%   See the Apache License for more details.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author : Dinh Ho Tong Minh (INRAE) and Yen Nhi Ngo, Jan. 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This function provides the combination of PS and DS targets
% see more detail in section 2 of [1]:
% [1] Dinh Ho Tong Minh and Yen Nhi Ngo.
% "Compressed SAR Interferometry in the Big Data Era". Remote Sensing.
% 2022, 14, 390. https://doi.org/10.3390/rs14020390
%
% This file can be used only for research purposes, you should cite
% the aforementioned papers in any resulting publication.
%

% Phase Linking
[phi_PL,Coh_cal] =  Intf_PL(Coh, 10,reference_ind);

% Phase filtering
[infstack_filt] = Intf_filt(interfstack,SHP,phi_PL,Coh_cal,reference_ind,BroNumthre,Cohthre);

% DeSpeckle
mli_despeckle = Image_DeSpeckle(abs(slcstack),SHP);
[SLCstack] = SLC_filt(mli_despeckle,slcstack,SHP,Coh_cal,BroNumthre,Cohthre_slc_filt);

% Export
switch InSAR_processor
    case 'snap' %
        Intf_export(infstack_filt,interflist,[InSAR_path,'/diff0'],'.psds');
        SLC_export(SLCstack,slclist,[InSAR_path,'/rslc'],'.psar');
    case 'isce'
        Intf_export(infstack_filt,interflist,InSAR_path,'.psds',InSAR_processor);
        SLC_export(SLCstack,slclist,InSAR_path,'.psar',InSAR_processor,reference_ind);
    otherwise
        disp('not yet support')
end

return

author: DinhHoTongMinh

repo: TomoSAR


text:
function R=estExpDecayReliability(t,T,method)
%%ESTEXPDECAYRELIABILITY Determine the probability that a particular
%       example of something will be functional after time T given
%       independent the failure times of n other examples under the
%       assumption that the time of failure given something is functional
%       at time zero is exponentially distributed, so the probability of
%       success (it lasting longer than T) given that it works at time zero
%       is Pr{T>t}=exp(-t/theta) for some parameter theta.
%
%INPUTS: t A 1Xn or nX1 vector of the times when something failed. All
%          elements should be >=0.
%        T The desired time after which a failure is acceptable. T>=0.
%   method The estimation approach. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) Use the
%            expected value estimate.
%          1 Use the maximum likelihood estimate. This estimate is biased.
%
%OUTPUTS: R The probability that  the system will still be operational at
%           some point after time T.
%
%This implements the methods of [1]. The basic reliability model is
%R=p*P
%where p is the probability of a failure at time 0 and P is the probability
%of a failure between time 0 and T. THe exponential decay model applies to
%P.
%
%EXAMPLE:
%We demonstrate here that the value of R at time T is more reliable when
%computing using method 0 than method 1 given n=30 samples of the data.
% p=0.9;
% T=2;
% theta=3;
% lambda=1/theta;
% RTrue=p*(1-ExponentialD.CDF(T,lambda))
%
% numRuns=1e4;
% R0=0;
% R1=0;
% for curRun=1:numRuns
%     n=30;
%     t=zeros(n,1);
%     for k=1:n
%         if(rand()<p)%If working at time 0.
%            t(k)=ExponentialD.rand(1,lambda);
%         else
%            t(k)=0;
%         end
%     end
%     R0=R0+estExpDecayReliability(t,T,0);
%     R1=R1+estExpDecayReliability(t,T,1);
% end
% R0=R0/numRuns
% R1=R1/numRuns
%Typically, R0 will be closet to RTrue than R1.
%
%REFERENCES:
%[1] E. L. Pugh, "The best estimate of reliability in the exponential
%    case," Operations Research, vol. 11, no. 1, pp. 57-61, Jan.-Feb. 1963.
%
%June 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(method))
    method=0;
end

n=length(t);
sel=(t~=0);
p=sum(sel)/n;

%Keep only failure times that are >0.
t=t(sel);
n=length(t);

if(n==0)
   R=0;
   return;
end

theta=sum(t)/n;%Equation 7.
switch(method)
    case 0%The expected value solution.
        R=p*(1-T/(n*theta))^(n-1);%Equation 17 (with p inserted).
    case 1%The ML solution.
        R=p*exp(-T/theta);%Equation 6 (with p inserted).
    otherwise
        error('Unknown method specified.')
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [U,G] = shell(V,F,th,varargin)
  % SHELL Compute a thin shell around a triangle mesh (V,F) with thickness th.
  %
  % Inputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of triangle indices into V
  %   th  desired thickness of shell
  %   Optional:
  %     'Normals' followed by #V by 3 list of vertex normals to use to
  %     construct shell (will be multiplied by th): {area-weighted}
  % Outputs:
  %   U  2*#V by 3 list of output mesh vertex positions: V always comes first
  %   G  #G by 3 list of triangle indices into U
  %

  N = [];
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Normals'}, ...
    {'N'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if isempty(N)
    N = normalizerow(per_vertex_normals(V,F));
  end

  % scale by thickness
  N = N*th;

  O = outline(F);
  n = size(V,1);
  U = [V;V+N];
  G = [ ...
    fliplr(F); ...
    n+F; ...
    0+O(:,1) n+O(:,[2 1]); ...
    n+O(:,2) 0+O(:,1:2)];

end

author: alecjacobson

repo: gptoolbox


text:
cwpos = 1*ones(length(wordsbi), 1);
cwneg = 1*ones(length(wordsbi), 1);

nbind = train_ind;
nbsent = allSNumBi(nbind);
nblbl = labels(nbind);
ppos = sum(nblbl)/length(nblbl); pneg = 1 - ppos;

for i = 1:length(nbsent)
    nbsentraw = nbsent{i};
    nbsentraw=nbsentraw(nbsentraw>0);
    updates = unique(nbsentraw);
    counts = power(histc(nbsentraw, updates)', trainp);
     % counts = ones(size(histc(nbsentraw, updates)));
    if nblbl(i) == 1
        cwpos(updates) = cwpos(updates) + counts * (0.5/ppos);
    else
        cwneg(updates) = cwneg(updates) + counts * (0.5/pneg);
    end
end
% cwpos = cwpos * sum(cwneg) ./ sum(cwpos);
pwpos = log(cwpos / sum(cwpos));
pwneg = log(cwneg / sum(cwneg));

% cwpos = cwpos ./ norm(cwpos);
% cwneg = cwneg ./ norm(cwneg);

nbindtst = test_ind;
nbsenttst = allSNumBi(nbindtst);
nblbltst = labels(nbindtst);

pred = zeros(1, length(nbsenttst));
% predvote = zeros(1, length(nbsenttst));
% load('../data/sentiLex.mat', 'sentlex', 'sentlexval');
%
% sentimap = containers.Map(sentlex, sentlexval);
% wordSent = zeros(length(words),1);
% for i = 1:length(words)
%         if sentimap.isKey(words{i}) && sentimap(words{i}) <= 1
%            wordSent(i) = 2*sentimap(words{i}) - 1;
%         end
% end

softpred = zeros(size(pred));
bipred = zeros(1, length(nbsenttst));
unipred = zeros(1, length(nbsenttst));
for i = 1:length(nbsenttst)
    nbsentraw = nbsenttst{i};
    nbsentraw=nbsentraw(nbsentraw>0);

    updates = unique(nbsentraw);
    counts = power(histc(nbsentraw, updates), testp)';

    sentend = floor(length(nbsentraw)/2);
    sentword = nbsentraw(1:sentend);
    sentbi = nbsentraw(sentend+1:end);

    predsent = cwpos([sentword]) ./ cwneg([sentword]);
    totsent = cwpos([sentword]) + cwneg([sentword]);
    predbi = cwpos([sentbi]) ./ cwneg([sentbi]);

    weightsent = 1 - flukeprob(cwpos([sentword]),  cwneg([sentword]));
    weightbi = 1 - flukeprob(cwpos([sentbi]),  cwneg([sentbi]));

    predval = sum(log(predsent) ) + ...
        sum(log(predbi)) + log(ppos) - log(pneg) ;

    bipred(i) =  prod(predbi) > 1;
    unipred(i) =  prod(predsent) > 1;
    pred(i) = predval > 0;

%     softpredp(i) = poslog;
%     softpredn(i) = neglog;
    vis = 0;

    if vis && pred(i) ~= nblbltst(i)
        if nblbltst(i) ==0
            i;
        end
    for s = 1:length(sentword)
    fprintf('%s\t\t%d\t%d\t%f\t%f\n', wordsbi{sentword(s)},...
    cwpos(sentword(s)) , cwneg(sentword(s)), predsent(s), weightsent(s))
    end

    for s = 1:length(sentbi)
    fprintf('%s\t\t%d\t%d\t%f\t%f\n', wordsbi{sentbi(s)},...
    cwpos(sentbi(s)) , cwneg(sentbi(s)), predbi(s), weightbi(s))
    end

    subplot(2,1,1)
    plot(log(cwpos(sentword) ./ cwneg(sentword)) .* weightsent);
    set(gca,'XTickLabel',wordsbi(sentword))
    set(gca,'XTick',1:sentend)

    subplot(2,1,2)
    plot(log(cwpos(sentbi) ./ cwneg(sentbi)) .* weightbi);
    set(gca,'XTickLabel',wordsbi(sentbi))
    set(gca,'XTick',1:sentend+1)
    end
%     predvote(i) = sum(wordSent(nbsentraw));
end

sum(nblbltst == pred)
% sum(nblbltst == (predvote > 0))

% indvoted = (predvote ~= 0);
% sum(nblbltst(indvoted) == (predvote(indvoted)>0))

%
test_sents = allSStr(test_ind);
fid = fopen('../output/nbtestexamplesc.txt','w');
for i = 1:length(test_sents)
    nbsentraw = nbsenttst{i};
    nbsentraw=nbsentraw(nbsentraw>0);
    poslog = sum(pwpos(nbsentraw));
    neglog = sum(pwneg(nbsentraw));

    fprintf(fid,'%s\n', cell2str(test_sents{i}) );
end
fclose(fid);

c = sum(nblbltst == pred)
l = length(pred)
rate = c/l

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
function [tc,relres,iter,xrec] = gabgrouplasso(x,g,a,M,lambda,varargin)
%GABGROUPLASSO  Group LASSO regression in Gabor domain
%   Usage: [tc,xrec] = gabgrouplasso(x,g,a,M,group,lambda,C,maxit,tol)
%
%   `gabgrouplasso` has been deprecated. Please use |franagrouplasso| instead.
%
%   A call to `gabgrouplasso(x,g,a,M,lambda)` can be replaced by ::
%
%     F=frame('dgt',g,a,M);
%     tc=franagrouplasso(F,lambda);
%
%   Any additional parameters passed to `gabgrouplasso` can be passed to
%   |franagrouplasso| in the same manner.
%
%   See also: frame, franagrouplasso

warning(['LTFAT: GABGROUPLASSO has been deprecated, please use FRANAGROUPLASSO ' ...
          'instead. See the help on FRANAGROUPLASSO for more details.']);

F=newframe('dgt',[],g,a,M);
[tc,relres,iter,xrec] = framegrouplasso(F,lambda,varargin{:});

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [mu, sigma, sigma_points] = prediction_step(mu, sigma, u)
% Updates the belief concerning the robot pose according to the motion model.
% mu: state vector containing robot pose and poses of landmarks obeserved so far
% Current robot pose = mu(1:3)
% Note that the landmark poses in mu are stacked in the order by which they were observed
% sigma: the covariance matrix of the system.
% u: odometry reading (r1, t, r2)
% Use u.r1, u.t, and u.r2 to access the rotation and translation values

% For computing lambda.
global scale;

% Compute sigma points
sigma_points = compute_sigma_points(mu, sigma);

% Dimensionality
n = length(mu);
% lambda
lambda = scale - n;

% TODO: Transform all sigma points according to the odometry command
% Remember to vectorize your operations and normalize angles
% Tip: the function normalize_angle also works on a vector (row) of angles
for i=1:2*n+1
    sigma_points(1:3,i) = sigma_points(1:3,i) + [u.t*cos(sigma_points(3,i) + u.r1); u.t*sin(sigma_points(3,i) + u.r1); u.r1 + u.r2];
    sigma_points(3,i) = normalize_angle(sigma_points(3,i));
endfor
% Computing the weights for recovering the mean
wm = [lambda/scale, repmat(1/(2*scale),1,2*n)];
wc = wm;

% -------- My initialization ----------- %
xbar = 0;
ybar = 0;
mu = zeros(n,1);
sigma = zeros(n);
% TODO: recover mu.
% Be careful when computing the robot's orientation (sum up the sines and
% cosines and recover the 'average' angle via atan2)
for i=1:2*n+1
    mu = mu + wm(i)*sigma_points(:,i);
    xbar = xbar + wm(i)*cos(sigma_points(3,i));
    ybar = ybar + wm(i)*sin(sigma_points(3,i));
endfor
mu(3) = atan2(ybar,xbar);
mu(3) = normalize_angle(mu(3));

% TODO: Recover sigma. Again, normalize the angular difference
for i=1:2*n+1
    tmp = sigma_points(:,i) - mu;
    tmp(3) = normalize_angle(tmp(3));
    sigma = sigma + wc(i)*tmp*tmp';
endfor

% Motion noise
motionNoise = 0.1;
R3 = [motionNoise, 0, 0;
     0, motionNoise, 0;
     0, 0, motionNoise/10];
R = zeros(size(sigma,1));
R(1:3,1:3) = R3;

% TODO: Add motion noise to sigma
sigma = sigma + R;

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
% Point Distribution

colors = 'brgkcm';

if ~exist('n_ima')|~exist('fc'),
   fprintf(1,'No calibration data available.\n');
   return;
end;

check_active_images;

if ~exist(['ex_' num2str(ind_active(1)) ]),
   fprintf(1,'Need to calibrate before analysing reprojection error. Maybe need to load Calib_Results.mat file.\n');
   return;
end;

figure(6);

for kk=1:n_ima

	if exist(['x_' num2str(kk)]),

	if active_images(kk) & eval(['~isnan(x_' num2str(kk) '(1,1))']),

		eval(['plot(x_' num2str(kk) '(1,:),x_' num2str(kk) '(2,:),''' colors(rem(kk-1,6)+1) '+'');']);

		hold on;
	end;

	end;

end;

axis('equal');

axis([0 nx 0 ny]);

title1=pwd;
title1=strrep(title1,'_','\_');

title({'Point Distribution in Images',title1});

xlabel('x');

ylabel('y');

repo: ekfmonoslam


text:
function D = spm_eeg_combineplanar(S)
% Combine data from MEGPLANAR sensors
% FORMAT D = spm_eeg_combineplanar(S)
%
% S        - optional input struct
%  fields of S:
%   D        - MEEG object or filename
%   mode     -
%              'append'     - add combined channels to the origal channels
%              'replace'    - replace MEGPLANAR with combined [default]
%              'replacemeg' - replace all MEG channels with combined but
%                             keep non-MEG
%              'keep'       - only write out the combined channels
%
%   prefix   - prefix for the output file [default: 'P']
%
% Output:
% D        - MEEG object (also written on disk)
%
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_combineplanar.m 7132 2017-07-10 16:22:58Z guillaume $

SVNrev = '$Rev: 7132 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','MEG Combine planar'); spm('Pointer','Watch');

if ~isfield(S, 'prefix'),       S.prefix   = 'P';           end
if ~isfield(S, 'mode'),         S.mode     = 'replace';     end

D = spm_eeg_load(S.D);

isTF = strncmpi(D.transformtype,'TF',2);

chanset = spm_eeg_planarchannelset(D.chanlabels);

chanind = [];
labelnew = {};
for i = 1:size(chanset, 1)
    cind = D.indchannel(chanset(i, 1:2));
    if length(cind) == 2
        chanind  = [chanind cind];
        labelnew = [labelnew; chanset(i, end)];
    end
end

megind    = D.indchantype('MEG');
planarind = D.indchantype('MEGPLANAR');

% add a row of zeros;
chanind(3, end) = 0;
copyind = [];
switch S.mode
    case 'append'
        copyind = 1:D.nchannels;
        copyind = [copyind; copyind];
        chanind(3, :) = (D.nchannels+1):(D.nchannels+size(chanind, 2));
    case {'replace', 'replacemeg'}
        ind = 1;
        for i = 1:D.nchannels
            if any(i == planarind)
                k = find((chanind(1, :) == i) | (chanind(2, :) == i));
                if ~isempty(k) && (chanind(3, k) == 0)
                    chanind(3, k) = ind;
                    ind = ind + 1;
                end
            elseif isequal(S.mode, 'replace') || ~any(i == megind)
                copyind = [copyind [i ind]'];
                ind = ind + 1;
            end
        end
    case 'keep'
        copyind       = [];
        chanind(3, :) = 1:size(chanind, 2);
end

if isempty(copyind)
    Nchannels = max(chanind(3, :));
else
    Nchannels = max([chanind(3, :), copyind(2, :)]);
end

%-Generate new MEEG object with new files
%--------------------------------------------------------------------------
if isTF
    Dnew = clone(D, [S.prefix fname(D)], [Nchannels D.nfrequencies D.nsamples D.ntrials], 1);
else
    Dnew = clone(D, [S.prefix fname(D)], [Nchannels D.nsamples D.ntrials], 1);
end

if strcmp(D.type, 'continuous')
    %-Continuous data
    %----------------------------------------------------------------------
    blksz  = D.fsample;
    blknum = floor(D.nsamples/blksz);

    spm_progress_bar('Init', blknum, 'Data blocks processed');
    if blknum > 100, Ibar = floor(linspace(1, blknum,100));
    else Ibar = 1:blknum; end

    if isTF
        i = 1;
        while 1
            if i<=blknum
                Iblock  = ((i-1)*blksz + 1):(i*blksz);
            elseif (i == (blknum + 1)) && blknum*blksz < D.nsamples
                Iblock  = (blknum*blksz + 1):D.nsamples;
            else
                break;
            end

            planar1 = D(chanind(1, :), :, Iblock);
            planar2 = D(chanind(2, :), :, Iblock);

            Dnew(chanind(3, :), :, Iblock) = planar1 + planar2;

            if ~isempty(copyind)
                Dnew(copyind(2, :), :, Iblock) =  D(copyind(1, :), :, Iblock);
            end

            if any(Ibar == i), spm_progress_bar('Set', i); end

            i = i+1;
        end
    else
        i = 1;
        while 1
            if i<=blknum
                Iblock  = ((i-1)*blksz + 1):(i*blksz);
            elseif (i == (blknum + 1)) && blknum*blksz < D.nsamples
                Iblock  = (blknum*blksz + 1):D.nsamples;
            else
                break;
            end

            planar1 = D(chanind(1, :),  Iblock);
            planar2 = D(chanind(2, :),  Iblock);

            Dnew(chanind(3, :), Iblock) = sqrt(planar1.^2 + planar2.^2);

            if ~isempty(copyind)
                Dnew(copyind(2, :), Iblock) =  D(copyind(1, :), Iblock);
            end

            if any(Ibar == i), spm_progress_bar('Set', i); end

            i = i+1;
        end
    end
else
    %-Epoched data
    %----------------------------------------------------------------------
    spm_progress_bar('Init', D.ntrials, 'Trials processed');
    if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials, 100));
    else Ibar = 1:D.ntrials; end

    if isTF
        for i=1:D.ntrials
            planar1 = D(chanind(1, :), :, :, i);
            planar2 = D(chanind(2, :), :, :, i);

            Dnew(chanind(3, :),:, :, i) = planar1 + planar2;

            if ~isempty(copyind)
                Dnew(copyind(2, :), :, :, i) =  D(copyind(1, :), :, :, i);
            end

            if any(Ibar == i), spm_progress_bar('Set', i); end
        end
    else
        for i=1:D.ntrials
            planar1 = D(chanind(1, :), :, i);
            planar2 = D(chanind(2, :), :, i);

            Dnew(chanind(3, :), :, i) = sqrt(planar1.^2 + planar2.^2);

            if ~isempty(copyind)
                Dnew(copyind(2, :), :, i) =  D(copyind(1, :), :, i);
            end

            if any(Ibar == i), spm_progress_bar('Set', i); end
        end
    end
end

spm_progress_bar('Clear');

Dnew = chanlabels(Dnew, chanind(3, :), labelnew);
Dnew = badchannels(Dnew, chanind(3, :), badchannels(D, chanind(1, :)) | badchannels(D, chanind(2, :)));
Dnew = chantype(Dnew, chanind(3, :), 'MEGCOMB');
Dnew = units(Dnew, chanind(3, :), units(D, chanind(1, :)));
Dnew = coor2D(Dnew, chanind(3, :), 0.5*(coor2D(D, chanind(1,:))+coor2D(D, chanind(1,:))));

if ~isempty(copyind)
    Dnew = chanlabels(Dnew, copyind(2, :), D.chanlabels(copyind(1, :)));
    Dnew = badchannels(Dnew, copyind(2, :), badchannels(D, copyind(1, :)));
    Dnew = chantype(Dnew, copyind(2, :), D.chantype(copyind(1, :)));
    Dnew = units(Dnew, copyind(2, :), D.units(copyind(1, :)));
    Dnew = coor2D(Dnew, copyind(2, :),coor2D(D, copyind(1, :)));
end

%-Save the new M/EEG dataset
%--------------------------------------------------------------------------
Dnew = Dnew.history(mfilename, S);
save(Dnew);

D = Dnew;

%-Cleanup
%--------------------------------------------------------------------------
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','MEG Combine planar: done'); spm('Pointer','Arrow');

author: spm


text:
% Multistage omega scheme.

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory in Section 5.10 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program generates Figs. 5.23, 5.24 in the book

% Solution of (5.89)
% dy/dt + u dy/dx - D d^2u/dx^2 = (hom-1)*beta^2*D*cos(beta*(x-ut))
%	x = 0: Dirichlet,  x = 1: Neumann or absorbing boundary conditions
% Central or upwind differences
% Defect correction optional
% Uniform cell-centered grid

% Functions called: exact_solution, fL, fR

		% ...............Input.........................................
u = 1.1;	% Velocity
D = 0.0005;	% Diffusion coefficient
dt = 1/10;	% Time step
n = 100;	% Number of time steps
J = 30;		% Number of grid cells
central = 1;	% Enter 1 for central scheme or something else for upwind scheme
bc = 2;		% Enter 1 for inhomogeneous Neumann or 2 for homogeneous Neumann
hom = 1;	% Enter 1 for   homogeneous right-hand side
		%    or 2 for inhomogeneous right-hand side
defcor = 1;	% Enter 1 for defect correction or else something else
nalpha = 3;	% Parameter in exact solution: alpha = nalpha*pi
nbeta = 2;	% Parameter in exact solution: beta  = nbeta*pi
		%...............End of input...............................

omega = 1 - sqrt(0.5);			% Parameters in multistage omega scheme
palpha = (1 - 2*omega)/(1 - omega);
palphap = 1 - palpha;
omegap = 1 - 2*omega;

T = n*dt;	% Final time

		% Definition of indexing
%     x=0
% grid |---o---|---o---|---o---|---o---|
%      1   1   2   2   3           J

h = 1.0/J;			% Cell size
d = D*dt/(h*h);			% Diffusion number
sigma = u*dt/h;			% Courant number
meshpeclet = sigma/d;		% Mesh Peclet number

	% Preallocation of J*J tridiagonal matrix A
Ac = spdiags([ ones(J,1)  ones(J,1)  ones(J,1)], [-1 0 1]',J,J);  % Central
Au = spdiags([ ones(J,1)  ones(J,1)  ones(J,1)], [-1 0 1]',J,J);  % Upwind

% Central scheme
Ac(1,1) = sigma/2 + 3*d;  Ac(1,2) = sigma/2 - d;
for j = 2:J-1,
  Ac(j,j-1) = -sigma/2 - d;    Ac(j,j+1) = sigma/2 - d;
  Ac(j,j)   = - Ac(j,j-1) - Ac(j,j+1);
end
j = J; Ac(j,j-1) = -sigma/2 - d; Ac(j,j) = sigma/2 + d;

% Upwind scheme
Au(1,1) = sigma + 3*d;  Au(1,2) =  - d;
for j = 2:J-1,
  Au(j,j-1) = -sigma - d;  Au(j,j+1) = - d;
  Au(j,j)   = - Au(j,j-1) - Au(j,j+1);
end
j = J; Au(j,j-1) = -sigma - d;  Au(j,j) = sigma + d;

	% System to be solved is E*ynew = C*yold + rhs
z = zeros(J,1); e = ones(J,1);   ID = spdiags([z e z],[-1 0 1],J,J);
Ec1 = ID + palpha*omega*Ac;	Cc1 = ID - palphap*omega*Ac;
Ec2 = ID + palphap*omegap*Ac;	Cc2 = ID - palpha*omegap*Ac;
Eu1 = ID + palpha*omega*Au;	Cu1 = ID - palphap*omega*Au;
Eu2 = ID + palphap*omegap*Au;	Cu2 = ID - palpha*omegap*Au;

if defcor ~= 1		% No defect correction
  if central == 1		% Central scheme
    E1 = Ec1;    C1 = Cc1;    E2 = Ec2;    C2 = Cc2;
  else				% Upwind scheme
    E1 = Eu1;    C1 = Cu1;    E2 = Eu2;    C2 = Cu2;
  end
end

	% Preallocations
yold = zeros(J,1);		% Numerical solution at previous time
ynew = zeros(J,1);		% Numerical solution at new time

x = [1:J]'*h - h/2;		% Coordinates of cell centers

%for j = 1:J,  x(j)=(j-0.5)*h; end

t = 0; yold = exact_solution(t,x,D,nalpha,nbeta,u,hom-1);
rhs = zeros(J,1);		% Right-hand side
beta = nbeta*pi;

if defcor ~= 1		% No defect correction
for i = 1:n

% Stage 1
  rhs = (hom-1)*beta*beta*D*dt*omega*cos(beta*(x-u*t));
  rhs(1) = rhs(1) +  palphap*omega*(sigma + 2*d)*fL(t,D,nalpha,nbeta,u,hom-1)...
   + palpha*omega*(sigma + 2*d)*fL(t+omega*dt,D,nalpha,nbeta,u,hom-1);
  if central == 1		% Central scheme
    rhs(J) = rhs(J) -...
     palphap*omega*(sigma/2 - d)*h*fR(t,D,nalpha,nbeta,u,bc,hom-1)...
     - palpha*omega*(sigma/2 - d)*fR(t+omega*dt,D,nalpha,nbeta,u,bc,hom-1);
  else				% Upwind scheme
    rhs(J) = rhs(J) - palphap*omega*( - d)*h*fR(t,D,nalpha,nbeta,u,bc,hom-1)...
     - palpha*omega*( - d)*fR(t+ omega*dt,D,nalpha,nbeta,u,bc,hom-1);
  end
  rhs = rhs + C1*yold;  ynew = E1\rhs;

% Stage 2
  yold = ynew;
  rhs = (hom-1)*beta*beta*D*dt*omegap*cos(beta*(x-u*(t + dt - omega*dt)));
  rhs(1) = rhs(1) + ...
     palpha*omegap*(sigma + 2*d)*fL(t + omega*dt,D,nalpha,nbeta,u,hom-1)...
   + palphap*omegap*(sigma + 2*d)*fL(t + dt - omega*dt,D,nalpha,nbeta,u,hom-1);
  if central == 1		% Central scheme
    rhs(J) = rhs(J) -...
     palpha*omegap*(sigma/2 - d)*h*fR(t + omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
     - palphap*omegap*(sigma/2 - d)*fR(t+dt-omega*dt,D,nalpha,nbeta,u,bc,hom-1);
  else				% Upwind scheme
    rhs(J) = rhs(J) -...
      palpha*omegap*( - d)*h*fR(t + omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
     - palphap*omegap*( - d)*fR(t + dt- omega*dt,D,nalpha,nbeta,u,bc,hom-1);
  end
  rhs = rhs + C2*yold;  ynew = E2\rhs;

% Stage 3
  yold = ynew;
  rhs = (hom-1)*beta*beta*D*dt*omega*cos(beta*(x - u*(t + dt)));
  rhs(1) = rhs(1) +...
    palphap*omega*(sigma + 2*d)*fL(t+ dt- omega*dt,D,nalpha,nbeta,u,hom-1)...
   + palpha*omega*(sigma + 2*d)*fL(t+dt,D,nalpha,nbeta,u,hom-1);
  if central == 1		% Central scheme
    rhs(J) = rhs(J) -...
    palphap*omega*(sigma/2-d)*h*fR(t+ dt- omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
     - palpha*omega*(sigma/2 - d)*fR(t+dt,D,nalpha,nbeta,u,bc,hom-1);
  else				% Upwind scheme
    rhs(J) = rhs(J) -...
   palphap*omega*( - d)*h*fR(t+ dt- omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
     - palpha*omega*( - d)*fR(t+ dt,D,nalpha,nbeta,u,bc,hom-1);
  end
  rhs = rhs + C1*yold;  ynew = E1\rhs;
  t = t + dt;  yold = ynew;
end
else		% One step of defect correction
  dy = zeros(J,1);		% Correction
  drhs = zeros(J,1);		% Correction of right-hand side

  for i = 1:n
		% Stage 1: upwind step
    rhs = (hom-1)*beta*beta*D*dt*omega*cos(beta*(x-u*t));
    rhs(1) = rhs(1) +  palphap*omega*(sigma + 2*d)*fL(t,D,nalpha,nbeta,u,hom-1)...
     + palpha*omega*(sigma + 2*d)*fL(t+omega*dt,D,nalpha,nbeta,u,hom-1);

    drhs = rhs;		% Preparation of rhs for correction step
    drhs(J) = drhs(J) -...
     palphap*omega*(sigma/2 - d)*h*fR(t,D,nalpha,nbeta,u,bc,hom-1)...
     - palpha*omega*(sigma/2 - d)*fR(t+omega*dt,D,nalpha,nbeta,u,bc,hom-1);
    rhs(J) = rhs(J) - palphap*omega*( - d)*h*fR(t,D,nalpha,nbeta,u,bc,hom-1)...
       - palpha*omega*( - d)*fR(t+ omega*dt,D,nalpha,nbeta,u,bc,hom-1);
    rhs = rhs + Cu1*yold;    ynew = Eu1\rhs;

		% Stage 1: correction step
    drhs = drhs - Ec1*ynew + Cc1*yold;    dy = Eu1\drhs;    ynew = ynew + dy;

		% Stage 2: upwind step
    yold = ynew;
    rhs = (hom-1)*beta*beta*D*dt*omegap*cos(beta*(x-u*(t + dt - omega*dt)));
    rhs(1) = rhs(1) + ...
      palpha*omegap*(sigma + 2*d)*fL(t + omega*dt,D,nalpha,nbeta,u,hom-1)...
    + palphap*omegap*(sigma + 2*d)*fL(t + dt - omega*dt,D,nalpha,nbeta,u,hom-1);

    drhs = rhs;		% Preparation of rhs for correction step
    drhs(J) = drhs(J) -...
       palpha*omegap*(sigma/2 - d)*h*fR(t + omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
       - palphap*omegap*(sigma/2 - d)*fR(t+dt-omega*dt,D,nalpha,nbeta,u,bc,hom-1);
    rhs(J) = rhs(J) -...
       palpha*omegap*( - d)*h*fR(t + omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
       - palphap*omegap*( - d)*fR(t + dt- omega*dt,D,nalpha,nbeta,u,bc,hom-1);
    rhs = rhs + Cu2*yold;    ynew = Eu2\rhs;

		% Stage 2: correction step
    drhs = drhs - Ec2*ynew + Cc2*yold;    dy = Eu2\drhs;    ynew = ynew + dy;

		% Stage 3: upwind step
    yold = ynew;
    rhs = (hom-1)*beta*beta*D*dt*omega*cos(beta*(x - u*(t + dt)));
    rhs(1) = rhs(1) +...
    palphap*omega*(sigma + 2*d)*fL(t+ dt- omega*dt,D,nalpha,nbeta,u,hom-1)...
    + palpha*omega*(sigma + 2*d)*fL(t+dt,D,nalpha,nbeta,u,hom-1);

    drhs = rhs;		% Preparation of rhs for correction step
    drhs(J) =drhs(J) -...
      palphap*omega*(sigma/2-d)*h*fR(t+ dt- omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
       - palpha*omega*(sigma/2 - d)*fR(t+dt,D,nalpha,nbeta,u,bc,hom-1);
      rhs(J) = rhs(J) -...
      palphap*omega*( - d)*h*fR(t+ dt- omega*dt,D,nalpha,nbeta,u,bc,hom-1)...
       - palpha*omega*( - d)*fR(t+ dt,D,nalpha,nbeta,u,bc,hom-1);
    rhs = rhs + Cu1*yold;
    ynew = Eu1\rhs;

		% Stage 3: correction step
    drhs = drhs - Ec1*ynew + Cc1*yold;    dy = Eu1\drhs;    ynew = ynew + dy;
    t = t + dt;    yold = ynew;
  end
end

solex =  exact_solution(t,x,D,nalpha,nbeta,u,hom-1);

error = ynew - solex;		% Compute error norms
norm1 = norm(error,1)/J
norm2 = norm(error,2)/sqrt(J)
norminf = norm(error,inf)

figure(1), clf, hold on
xx = 0:0.02:1; plot (xx,exact_solution(t,xx,D,nalpha,nbeta,u,hom-1),'-');
axis([0 1 -1 1]); plot (x,ynew,'o');
if central == 1,  s1 = ['Central fractional step omega-scheme'];
else,             s1 = ['Upwind fractional step omega-scheme'];
end
if defcor == 1
  s1 = ['Fractional step omega-scheme with defect correction'];
end
if bc == 1,  s1 = [s1, ', inhomogeneous Neumann'];
else         s1 = [s1, ', homogeneous Neumann'];
end
title(s1,'fontsize',12);
if hom ==1,  s2 = ['Homogeneous differential equation']
else,        s2 = ['Inhomogeneous differential equation']
end
s3 = ['omega=',num2str(omega),'  D=',num2str(D),' u=',num2str(u),...
 ' h=',num2str(h),'  dt=',num2str(dt),'  t=',num2str(T),'  n=',num2str(n)]
s4 = ['meshpeclet=',num2str(meshpeclet),' Courant=',num2str(sigma),...
 '  d=',num2str(2*d) ]

repo: Aero-matlab


text:
%Image Compression

%compdct.m:
%main code to compress an image we only run compdct.m
%decompdct.m:
% to decompress an image we only run compdct.m
%resize.m:
%here we give an esample: if we have value big than  65792 like 19071001 in
%8 bits we can not save it so we use resize:
%X=19071001;
%Y=resize(X)
%>>Y
%1 35 0 25
%here X=1*256^3+35*256^2+0*256^1+25*256^0
%proba.m:
%like hist.m
%zigzag.m:
%zigzag scan for bloc [8 8]
%zigzaginv.m:
%inverse zigzag scan for bloc [8 8]
%zigzag16.m:
%zigzag scan for bloc [16 16]
%zigzinv16.m:
%inverse zigzag scan for bloc [16 16]
%zigzag32.m:
%zigzag scan for bloc [32 32]
%zigzinv32.m:
%inverse zigzag scan for bloc [32 32]
%rle.m:
%Run length encoding
%irle.m:
% Inverse Run length encoding
%abais.m:
%here we give an esample: if we have value big than 255 like 275 in
%8 bits we can not save it so we use resize:
%X=275;
%Y=resize(X)
%>>Y
%1 19
%here X=1*256^1+19*256^0
%Iabais.m:
%X=Iabais(y);

author: Sable

repo: mcbench-benchmarks


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

g = 1.0 ./ (1.0 + exp(-z));
end


text:
function M = stiefelgeneralizedfactory(n, p, B)
% Returns a manifold structure of "scaled" orthonormal matrices.
%
% function M = stiefelgeneralizedfactory(n, p)
% function M = stiefelgeneralizedfactory(n, p, B)
%
% The generalized Stiefel manifold is the set of "scaled" orthonormal
% nxp matrices X such that X'*B*X is identity. B must be positive definite.
% If B is identity, then this is the standard Stiefel manifold.
%
% The generalized Stiefel manifold is endowed with a scaled metric
% by making it a Riemannian submanifold of the Euclidean space,
% again endowed with the scaled inner product.
%
% Some notions (not all) are from Section 4.5 of the paper
% "The geometry of algorithms with orthogonality constraints",
% A. Edelman, T. A. Arias, S. T. Smith, SIMAX, 1998.
%
% Paper link: http://arxiv.org/abs/physics/9806030.
%
% Note: egrad2rgrad and ehess2rhess involve solving linear systems in B. If
% this is a bottleneck for a specific application, then a way forward is to
% create a modified version of this file which preprocesses B to speed this
% up (typically, by computing a Cholesky factorization of it, then calling
% an appropriate solver; or by exploiting sparsity.)
%
% See also: stiefelfactory  grassmannfactory  grassmanngeneralizedfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, June 30, 2015.
% Contributors: Hiroyuki Sato and Kensuke Aihara, Dec. 27, 2018.
%
% Change log:
%
%    Sep.  6, 2018 (NB):
%        Removed M.exp() as it was not implemented.
%
%    Dec. 27, 2018 (NB):
%        Added retraction M.retr_qr (based on QR factorization) as an
%        alternative to the previous retraction, which is now accessible as
%        M.retr_polar. The new retraction should be more efficient: it is
%        now the default, as M.retr = M.retr_qr. This new retraction is a
%        contribution of H. Sato and K. Aihara, as per their paper:
%        'Cholesky QR-based retraction on the generalized Stiefel manifold'
%        https://link.springer.com/article/10.1007%2Fs10589-018-0046-7

    if ~exist('B', 'var') || isempty(B)
        B = speye(n); % Standard Stiefel manifold.
    end

    M.name = @() sprintf('Generalized Stiefel manifold St(%d, %d)', n, p);

    M.dim = @() (n*p - .5*p*(p+1));

    M.inner = @(X, eta, zeta) trace(eta'*(B*zeta)); % Scaled metric.

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(X, Y) error('stiefelgeneralizedfactory.dist not implemented yet.');

    M.typicaldist = @() sqrt(p);

    % Orthogonal projection of an ambient vector U to the tangent space
    % at X.
    M.proj = @projection;
    function Up = projection(X, U)
        BX = B*X;

        % Projection onto the tangent space
        Up = U - X*symm(BX'*U);
    end

    M.tangent = M.proj;

    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)

        % First, scale egrad according the to the scaled metric in the
        % Euclidean space. Ideally, B should be preprocessed to ease
        % solving linear systems, e.g., via Cholesky factorization.
        egrad_scaled = B\egrad;

        % Second, project onto the tangent space.
        % rgrad = egrad_scaled - X*symm((B*X)'*egrad_scaled);
        %
        % Verify that symm(BX'*egrad_scaled) = symm(X'*egrad).

        rgrad = egrad_scaled - X*symm(X'*egrad);
    end

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)
        egraddot = ehess;
        Xdot = H;

        % Directional derivative of the Riemannian gradient.
        egrad_scaleddot = B\egraddot;
        rgraddot = egrad_scaleddot - Xdot*symm(X'*egrad) ...
                                   - X*symm(Xdot'*egrad) ...
                                   - X*symm(X'*egraddot);

        % Project onto the tangent space.
        rhess = M.proj(X, rgraddot);
    end

    M.retr_polar = @retraction_polar;
    function Y = retraction_polar(X, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = guf(X + t*U); % Ensures Y'*B*Y is identity.
    end

    M.retr_qr = @retraction_qr;
    function Y = retraction_qr(X, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = gqf(X + t*U); % Ensures Y'*B*Y is identity.
    end

    % By default, we use the QR retraction
    M.retr = M.retr_qr;

    M.hash = @(X) ['z' hashmd5(X(:))];

    M.rand = @random;
    function X = random()
        X = guf(randn(n, p)); % Ensures X'*B*X is identity;
        % TODO: test if this can be replaced by gqf.
    end

    M.randvec = @randomvec;
    function U = randomvec(X)
        U = projection(X, randn(n, p));
        U = U / norm(U(:));
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(X) zeros(n, p);

    % This transport is compatible with the generalized polar retraction.
    M.transp = @(X1, X2, d) projection(X2, d);

    M.vec = @(X, u_mat) u_mat(:);
    M.mat = @(X, u_vec) reshape(u_vec, [n, p]);
    M.vecmatareisometries = @() false;

    % Some auxiliary functions
    symm = @(D) (D + D')/2;

    function X = gqf(Y)
        % Generalized QR decomposition of an n-by-p matrix Y.
        % X'*B*X is identity. See Sato and Aihara 2018,
        % Cholesky QR-based retraction on the generalized Stiefel manifold,
        % https://link.springer.com/article/10.1007%2Fs10589-018-0046-7
        %
        % Comment by Nicolas Boumal, Dec. 27, 2018, following discussions
        % with Hiroyuki Sato, Kensuke Aihara and Bamdev Mishra:
        %
        % In principle, to orthonormalize the columns of Y against the
        % inner product defined by B, it would be better to implement (for
        % example) a modified Gram-Schmidt algorithm, possibly run twice.
        % Indeed, the latter would be affected by the condition number of
        % sqrtm(B)*Y. In contrast, the method below first squares the data,
        % hence is affected by the condition number of Y'*B*Y, which is the
        % square of that of sqrtm(B)*Y. Fortunately, it is easily shown
        % that, when Y = X + U for a tangent vector U at X, the condition
        % number of sqrtm(B)*Y is upper bounded by the square root of
        % 1 + ||U||_X^2. In Riemannian optimization, it is seldom necessary
        % to retract very large vectors, hence it is expected that the
        % condition numbers encountered in practice will be reasonable. As
        % a result, the implementation below which calls directly upon
        % low-level routines (Cholesky and a small triangular system solve)
        % is expected to run faster in practice than a home-made
        % implementation of modified Gram-Schmidt, which would involve a
        % loop, and still be sufficiently accurate.
        R = chol(symm(Y'*(B*Y)));
        X = Y / R;
    end

    function X = guf(Y)
        % Generalized polar decomposition of an n-by-p matrix Y.
        % X'*B*X is identity.

        % Method 1
        [u, ~, v] = svd(Y, 0);

        % Instead of the following three steps, an equivalent, but an
        % expensive way is to do X = u*(sqrtm(u'*(B*u))\(v')).
        [q, ssquare] = eig(u'*(B*u));
        qsinv = q/sparse(diag(sqrt(diag(ssquare))));
        X = u*((qsinv*q')*v'); % X'*B*X is identity.

        % Another computation using restricted_svd
        % [u, ~, v] = restricted_svd(Y);
        % X = u*v'; % X'*B*X is identity.

    end

    function [u, s, v] = restricted_svd(Y) %#ok<DEFNU>
        % We compute a thin svd-like decomposition of an n-by-p matrix Y
        % into matrices u, s, and v such that u is an n-by-p matrix
        % with u'*B*u being identity, s is a p-by-p diagonal matrix
        % with positive entries, and v is a p-by-p orthogonal matrix.
        % Y = u*s*v'.
        [v, ssquare] = eig(symm(Y'*(B*Y))); % Y*B*Y is positive definite
        ssquarevec = diag(ssquare);

        s = sparse(diag(abs(sqrt(ssquarevec))));
        u = Y*(v/s); % u'*B*u is identity.
    end

end

author: NicolasBoumal

repo: manopt


text:
classdef InstanceSegMMLoss_randSample < dagnn.Loss
    properties
        marginAlpha_
        adaptiveMM = false;
    end
    properties (Transient)
        curSimMat
        weightMat
        marginMat
        mass_
        size_
	lastLayerName
    end

    methods
        function outputs = forward(obj, inputs, params)
            [h, w, ch, bs] = size(inputs{1});
            sz = [h, w, ch, bs];
            mass = sz(1) * sz(2) + 1;

            obj.size_ = sz;
            obj.mass_ = mass;

            layerName = obj.lastLayerName; %'res7_cosSim';
            if isnan(obj.net.getLayerIndex(layerName))
                layerName = 'obj_instSeg_reg';
            end
            obj.curSimMat = obj.net.layers(obj.net.getLayerIndex(layerName)).block.curSimMat;
            obj.weightMat = obj.net.layers(obj.net.getLayerIndex(layerName)).block.weightMat;

            if obj.adaptiveMM
                obj.marginMat = obj.net.layers(obj.net.getLayerIndex(layerName)).block.marginMat;
                obj.marginMat = max(obj.marginMat*obj.marginAlpha_, obj.marginAlpha_*0.2);
                outputs{1} = vl_nnloss_modified(inputs{1}, obj.curSimMat, 'marginAlpha_', obj.marginAlpha_, ...
                'loss', obj.loss, ...
                'instanceWeights', obj.weightMat, 'marginMat', obj.marginMat ) ; % 1./mass
            else
                outputs{1} = vl_nnloss_modified(inputs{1}, obj.curSimMat, 'marginAlpha_', obj.marginAlpha_, ...
                    'loss', obj.loss, ...
                    'instanceWeights', obj.weightMat ) ; % 1./mass
            end
            n = obj.numAveraged ;
            m = n + size(inputs{1},4) ;
            obj.average = (n * obj.average + double(gather(outputs{1}))) / m ;
            obj.numAveraged = m ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            derInputs{1} = vl_nnloss_modified(inputs{1}, obj.curSimMat, derOutputs{1}, 'marginAlpha_', obj.marginAlpha_, ...
                'loss', obj.loss, ...
                'instanceWeights', obj.weightMat ) ; % 1./mass
            derInputs{2} = [] ;
            derParams = {} ;
        end

        function obj = InstanceSegMMLoss_randSample(varargin)
            obj.load(varargin) ;
        end
    end
end

author: aimerykong


text:
function [SPM] = cl_ext_spm_spm(SPM)
% [Re]ML Estimation of a General Linear Model
%
% :Usage:
% ::
%
%     FORMAT [SPM] = spm_spm(SPM)
%
% :Required fields of SPM:
%
% xY.VY - nScan x 1 struct array of image handles (see spm_vol)
%         Images must have the same orientation, voxel size and data type
%       - Any scaling should have already been applied via the image handle
%         scalefactors.
%
% xX    - Structure containing design matrix information
%       - Required fields are:
%         xX.X      - Design matrix (raw, not temporally smoothed)
%         xX.name   - cellstr of parameter names corresponding to columns
%                     of design matrix
%       - Optional fields are:
%         xX.K      - cell of session-specific structures (see spm_filter)
%                   - Design & data are pre-multiplied by K
%                     (K*Y = K*X*beta + K*e)
%                   - Note that K should not smooth across block boundaries
%                   - defaults to speye(size(xX.X,1))
%         xX.W      - Optional whitening/weighting matrix used to give
%                     weighted least squares estimates (WLS). If not specified
%                     spm_spm will set this to whiten the data and render
%                     the OLS estimates maximum likelihood
%                     i.e. W*W' = inv(xVi.V).
%
% xVi   - Structure describing intrinsic temporal non-sphericity
%       - Required fields are:
%         xVi.Vi    - array of non-sphericity components
%                   - defaults to {speye(size(xX.X,1))} - i.i.d.
%                   - specifying a cell array of constraints (Qi)
%                     These constraints invoke spm_reml to estimate
%                     hyperparameters assuming V is constant over voxels.
%                     that provide a high precise estimate of xX.V
%       - Optional fields are:
%         xX.V      - Optional non-sphericity matrix.  Cov(e) = sigma^2*V
%                     If not specified spm_spm will compute this using
%                     a 1st pass to identify significant voxels over which
%                     to estimate V.  A 2nd pass is then used to re-estimate
%                     the parameters with WLS and save the ML estimates
%                     (unless xX.W is already specified).
%
% xM    - Structure containing masking information, or a simple column vector
%         of thresholds corresponding to the images in VY [default: -Inf]
%       - If a structure, the required fields are:
%         xM.TH - nVar x nScan matrix of analysis thresholds, one per image
%         xM.I  - Implicit masking (0=>none, 1 => implicit zero/NaN mask)
%         xM.VM - struct array of explicit mask image handles
%       - (empty if no explicit masks)
%               - Explicit mask images are >0 for valid voxels to assess.
%               - Mask images can have any orientation, voxel size or data
%                 type. They are interpolated using nearest neighbour
%                 interpolation to the voxel locations of the data Y.
%       - Note that voxels with constant data (i.e. the same value across
%         scans) are also automatically masked out.
%
% swd   - Directory where the output files will be saved [default: pwd]
%         If exists, it becomes the current working directory.
%
% In addition, global SPM "defaults" variable is used (see spm_defaults):
%
% stats.<modality>.UFp - critical F-threshold for selecting voxels over
%                        which the non-sphericity is estimated (if
%                        required) [default: 0.001]
%
% stats.maxres         - maximum number of residual images for smoothness
%                        estimation
%
% stats.maxmem         - maximum amount of data processed at a time (in bytes)
%
% modality             - SPM modality {'PET','FMRI','EEG'}
%
%__________________________________________________________________________
%
% spm_spm is the heart of the SPM package. Given image files and a
% General Linear Model, it estimates the model parameters, variance
% hyperparameters, and smoothness of standardised residual fields, writing
% these out to disk in the current working directory for later
% interrogation in the results section. (NB: Existing analyses in the
% current working directory are overwritten).  This directory
% now becomes the working directory for this analysis and all saved
% images are relative to this directory.
%
% The model is expressed via the design matrix (xX.X). The basic model
% at each voxel is of the form is Y = X*B + e, for data Y, design
% matrix X, (unknown) parameters B and residual errors e. The errors
% are assumed to have a normal distribution.
%
% Sometimes confounds (e.g. drift terms in fMRI) are necessary. These
% can be specified directly in the design matrix or implicitly, in terms
% of a residual forming matrix K to give a generalised linear model
% K*Y = K*X*B + K*e.  In fact K can be any matrix (e.g. a convolution
% matrix).
%
% In some instances i.i.d. assumptions about errors do not hold. For
% example, with serially correlated (fMRI) data or correlations among the
% levels of a factor in repeated measures designs. This non-sphericity
% can be specified in terms of components (SPM.xVi.Vi{i}). If specified
% these covariance components will then be estimated with ReML (restricted
% maximum likelihood) hyperparameters. This estimation assumes the same
% non-sphericity for voxels that exceed the global F-threshold. The ReML
% estimates can then be used to whiten the data giving maximum likelihood
% (ML) or Gauss-Markov estimators. This entails a second pass of the data
% with an augmented model K*W*Y = K*W*X*B + K*W*e where W*W' = inv(xVi.V).
% xVi.V is the non-sphericity based on the hyperparameter estimates.
% W is stored in xX.W and cov(K*W*e) in xX.V. The covariance of the
% parameter estimates is then xX.Bcov = pinv(K*W*X)*xX.V*pinv(K*W*X)'.
%
% If you do not want ML estimates but want to use ordinary least squares
% (OLS) then simply set SPM.xX.W to the identity matrix. Any non-sphericity
% V will still be estimated but will be used to adjust the degrees of freedom
% of the ensuing statistics using the Satterthwaite approximation (c.f.
% the Greenhouse-Geisser corrections).
%
% If [non-spherical] variance components Vi are not specified xVi.Vi and
% xVi.V default to the identity matrix (i.e. i.i.d). The parameters are
% then estimated by OLS.  In this instance the OLS and ML estimates are
% the same.
%
% Note that only a single voxel-specific hyperparameter (i.e. variance
% component) is estimated, even if V is not i.i.d.  This means spm_spm
% always implements a fixed-effects model.
% Random effects models can be emulated using a multi-stage procedure:
% This entails summarising the data with contrasts such that the fixed
% effects in a second model on the summary data are those effects of
% interest (i.e. the population effects). This means contrasts are
% re-entered into spm_spm to make an inference (SPM) at the next
% level. At this higher hierarchical level the residual variance for the
% model contains the appropriate variance components from lower levels.
% See spm_RandFX.man for further details and below.
%
% Under the additional assumption that the standardised error fields
% are non-stationary standard Gaussian random fields, results from
% Random field theory can be applied to estimate the significance
% statistic images (SPM's) adjusting p values for the multiple tests
% at all voxels in the search volume. The parameters required for
% this random field correction are the volume, and Lambda, the covariance
% matrix of partial derivatives of the standardised error fields, estimated
% by spm_est_smoothness.
%
%                           ----------------
%
% The volume analysed is the intersection of the threshold masks,
% explicit masks and implicit masks. See spm_spm_ui for further details
% on masking options.
%
%--------------------------------------------------------------------------
%
% The output of spm_spm takes the form of an SPM.mat file of the analysis
% parameters, and 'float' flat-file images of the parameter and variance
% [hyperparameter] estimates. An 8bit zero-one mask image indicating the
% voxels assessed is also written out, with zero indicating voxels outside
% tha analysed volume.
%
%                           ----------------
%
% The following SPM.fields are set by spm_spm (unless specified)
%
%     xVi.V      - estimated non-sphericity trace(V) = rank(V)
%     xVi.h      - hyperparameters  xVi.V = xVi.h(1)*xVi.Vi{1} + ...
%     xVi.Cy     - spatially whitened <Y*Y'> (used by ReML to estimate h)
%     xVi.CY     - <(Y - <Y>)*(Y - <Y>)'>    (used by spm_spm_Bayes)
%
%                           ----------------
%
%     Vbeta     - struct array of beta image handles (relative)
%     VResMS    - file struct of ResMS image handle  (relative)
%     VM        - file struct of Mask  image handle  (relative)
%
%                           ----------------
%
%     xX.W      - if not specified W*W' = inv(x.Vi.V)
%     xX.V      - V matrix (K*W*Vi*W'*K') = correlations after K*W is applied
%     xX.xKXs   - space structure for K*W*X, the 'filtered and whitened'
%                 design matrix
%               - given as spm_sp('Set',xX.K*xX.W*xX.X) - see spm_sp
%     xX.pKX    - pseudoinverse of K*W*X, computed by spm_sp
%     xX.Bcov   - xX.pKX*xX.V*xX.pKX - variance-covariance matrix of
%                 parameter estimates
%                 (when multiplied by the voxel-specific hyperparameter ResMS
%                 of the parameter estimates (ResSS/xX.trRV = ResMS) )
%     xX.trRV   - trace of R*V
%     xX.trRVRV - trace of RVRV
%     xX.erdf   - effective residual degrees of freedom (trRV^2/trRVRV)
%     xX.nKX    - design matrix (xX.xKXs.X) scaled for display
%                 (see spm_DesMtx('sca',... for details)
%
%                           ----------------
%
%     xVol.M    - 4x4 voxel->mm transformation matrix
%     xVol.iM   - 4x4 mm->voxel transformation matrix
%     xVol.DIM  - image dimensions - column vector (in voxels)
%     xVol.XYZ  - 3 x S vector of in-mask voxel coordinates
%     xVol.S    - Lebesgue measure or volume       (in voxels)
%     xVol.R    - vector of resel counts           (in resels)
%     xVol.FWHM - Smoothness of components - FWHM, (in voxels)
%
%                           ----------------
%
%     xCon      - Contrast structure (created by spm_FcUtil.m)
%     xCon.name - Name of contrast
%     xCon.STAT - 'F', 'T' or 'P' - for F/T-contrast ('P' for PPMs)
%     xCon.c    - (F) Contrast weights
%     xCon.X0   - Reduced design matrix (spans design space under Ho)
%                 It is in the form of a matrix (spm99b) or the
%                 coordinates of this matrix in the orthogonal basis
%                 of xX.X defined in spm_sp.
%     xCon.iX0  - Indicates how contrast was specified:
%                 If by columns for reduced design matrix then iX0 contains
%                 the column indices. Otherwise, it's a string containing
%                 the spm_FcUtil 'Set' action: Usually one of {'c','c+','X0'}
%                 (Usually this is the input argument F_iX0.)
%     xCon.X1o  - Remaining design space (orthogonal to X0).
%                 It is in the form of a matrix (spm99b) or the
%                 coordinates of this matrix in the orthogonal basis
%                 of xX.X defined in spm_sp.
%     xCon.eidf - Effective interest degrees of freedom (numerator df)
%     xCon.Vcon - ...for handle of contrast/ESS image (empty at this stage)
%     xCon.Vspm - ...for handle of SPM image (empty at this stage)
%
%                           ----------------
%
%
% The following images are written to file
%
% mask.{img,hdr}                                   - analysis mask image
% 8-bit (uint8) image of zero-s & one's indicating which voxels were
% included in the analysis. This mask image is the intersection of the
% explicit, implicit and threshold masks specified in the xM argument.
% The XYZ matrix contains the voxel coordinates of all voxels in the
% analysis mask. The mask image is included for reference, but is not
% explicitly used by the results section.
%
%                           ----------------
%
% beta_????.{img,hdr}                                 - parameter images
% These are 32-bit (float32) images of the parameter estimates. The image
% files are numbered according to the corresponding column of the
% design matrix. Voxels outside the analysis mask (mask.img) are given
% value NaN.
%
%                           ----------------
%
% ResMS.{img,hdr}                    - estimated residual variance image
% This is a 64-bit (float64) image of the residual variance estimate.
% Voxels outside the analysis mask are given value NaN.
%
%                           ----------------
%
% RPV.{img,hdr}                      - estimated resels per voxel image
% This is a 64-bit (float64) image of the RESELs per voxel estimate.
% Voxels outside the analysis mask are given value 0.  These images
% reflect the nonstationary aspects the spatial autocorrelations.
%
%                           ----------------
%
% ResI_????.{img,hdr}        - standardised residual (temporary) images
% These are 64-bit (float64) images of standardised residuals. At most
% maxres images will be saved and used by spm_est_smoothness, after which
% they will be deleted.
%
%--------------------------------------------------------------------------
%
% References:
%
% Christensen R (1996) Plane Answers to Complex Questions
%       Springer Verlag
%
% Friston KJ, Holmes AP, Worsley KJ, Poline JB, Frith CD, Frackowiak RSJ (1995)
% ``Statistical Parametric Maps in Functional Imaging:
%   A General Linear Approach''
%       Human Brain Mapping 2:189-210
%
% Worsley KJ, Friston KJ (1995)
% ``Analysis of fMRI Time-Series Revisited - Again''
%       NeuroImage 2:173-181
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
%
% ..
%    % Andrew Holmes, Jean-Baptiste Poline & Karl Friston
%    $Id: spm_spm.m 3960 2010-06-30 17:41:24Z ged $
% ..

SVNid     = '$Rev: 3960 $';

%-Say hello
%--------------------------------------------------------------------------
SPMid     = spm('FnBanner',mfilename,SVNid);
Finter    = spm('FigName','Stats: estimation...'); spm('Pointer','Watch');

%-Get SPM.mat[s] if necessary
%--------------------------------------------------------------------------
if nargin == 0
    P     = cellstr(spm_select(Inf,'^SPM\.mat$','Select SPM.mat[s]'));
    for i = 1:length(P)
        swd     = fileparts(P{i});
        load(fullfile(swd,'SPM.mat'));
        SPM.swd = swd;
        spm_spm(SPM);
    end
    return
end

%-Change to SPM.swd if specified
%--------------------------------------------------------------------------
try
    cd(SPM.swd);
catch
    SPM.swd = pwd;
end

%-Ensure data are assigned
%--------------------------------------------------------------------------
try
    SPM.xY.VY;
catch
    spm('alert!','Please assign data to this design', mfilename);
    spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
    return
end

%-Delete files from previous analyses
%--------------------------------------------------------------------------
if exist(fullfile(SPM.swd,'mask.img'),'file') == 2

    str = {'Current directory contains SPM estimation files:',...
        'pwd = ',SPM.swd,...
        'Existing results will be overwritten!'};
    if spm_input(str,1,'bd','stop|continue',[1,0],1)
        spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
        return
    else
        warning('Overwriting old results\n\t (pwd = %s) ',SPM.swd);
        try, SPM = rmfield(SPM,'xVol'); end
    end
end

files = {'^mask\..{3}$','^ResMS\..{3}$','^RPV\..{3}$',...
         '^beta_.{4}\..{3}$','^con_.{4}\..{3}$','^ResI_.{4}\..{3}$',...
         '^ess_.{4}\..{3}$', '^spm\w{1}_.{4}\..{3}$'};

for i = 1:length(files)
    j = spm_select('List',SPM.swd,files{i});
    for k = 1:size(j,1)
        spm_unlink(deblank(j(k,:)));
    end
end

%==========================================================================
% - A N A L Y S I S   P R E L I M I N A R I E S
%==========================================================================

%-Initialise
%==========================================================================
fprintf('%-40s: %30s','Initialising parameters','...computing');        %-#
xX            = SPM.xX;
[nScan nBeta] = size(xX.X);

%-If xM is not a structure then assume it's a vector of thresholds
%--------------------------------------------------------------------------
try
    xM = SPM.xM;
catch
    xM = -Inf(nScan,1);
end
if ~isstruct(xM)
    xM = struct('T',    [],...
                'TH',   xM,...
                'I',    0,...
                'VM',   {[]},...
                'xs',   struct('Masking','analysis threshold'));
end

%-Check confounds (xX.K) and non-sphericity (xVi)
%--------------------------------------------------------------------------
if ~isfield(xX,'K')
    xX.K  = 1;
end
try
    %-If covariance components are specified use them
    %----------------------------------------------------------------------
    xVi   = SPM.xVi;
catch

    %-otherwise assume i.i.d.
    %----------------------------------------------------------------------
    xVi   = struct( 'form',  'i.i.d.',...
                    'V',     speye(nScan,nScan));
end

%-Get non-sphericity V
%==========================================================================
try
    %-If xVi.V is specified proceed directly to parameter estimation
    %----------------------------------------------------------------------
    V     = xVi.V;
    str   = 'parameter estimation';

catch

    % otherwise invoke ReML selecting voxels under i.i.d assumptions
    %----------------------------------------------------------------------
    V     = speye(nScan,nScan);
    str   = '[hyper]parameter estimation';
end

%-Get whitening/Weighting matrix: If xX.W exists we will save WLS estimates
%--------------------------------------------------------------------------
try
    %-If W is specified, use it
    %----------------------------------------------------------------------
    W     = xX.W;
catch

    if isfield(xVi,'V')

        % otherwise make W a whitening filter W*W' = inv(V)
        %------------------------------------------------------------------
        W     = spm_sqrtm(spm_inv(xVi.V));
        W     = W.*(abs(W) > 1e-6);
        xX.W  = sparse(W);

    else
        % unless xVi.V has not been estimated - requiring 2 passes
        %------------------------------------------------------------------
        W     = speye(nScan,nScan);
        str   = 'hyperparameter estimation (1st pass)';
    end
end

%-Design space and projector matrix [pseudoinverse] for WLS
%==========================================================================
xX.xKXs   = spm_sp('Set',spm_filter(xX.K,W*xX.X));       % KWX
xX.xKXs.X = full(xX.xKXs.X);
xX.pKX    = spm_sp('x-',xX.xKXs);                        % projector
erdf      = spm_SpUtil('trRV',xX.xKXs);                  % Working error df

%-If xVi.V is not defined compute Hsqr and F-threshold under i.i.d.
%--------------------------------------------------------------------------
if ~isfield(xVi,'V')

    Fcname = 'effects of interest';
    iX0    = [SPM.xX.iB SPM.xX.iG];
    xCon   = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
    X1o    = spm_FcUtil('X1o', xCon(1),xX.xKXs);
    Hsqr   = spm_FcUtil('Hsqr',xCon(1),xX.xKXs);
    trRV   = spm_SpUtil('trRV',xX.xKXs);
    trMV   = spm_SpUtil('trMV',X1o);

    % Threshold for voxels entering non-sphericity estimates
    %----------------------------------------------------------------------
    try
        modality = lower(spm_get_defaults('modality'));
        UFp      = spm_get_defaults(['stats.' modality '.ufp']);
    catch
        UFp      = 0.001;
    end
    UF           = spm_invFcdf(1 - UFp,[trMV,trRV]);
end

%-Image dimensions and data
%==========================================================================
VY       = SPM.xY.VY;
spm_check_orientations(VY);

% check files exists and try pwd
%--------------------------------------------------------------------------
for i = 1:numel(VY)
    if ~spm_existfile(VY(i).fname)
        [p,n,e]     = fileparts(VY(i).fname);
        VY(i).fname = [n,e];
    end
end

M        = VY(1).mat;
DIM      = VY(1).dim(1:3)';
xdim     = DIM(1); ydim = DIM(2); zdim = DIM(3);
YNaNrep  = spm_type(VY(1).dt(1),'nanrep');

%-Maximum number of residual images for smoothness estimation
%--------------------------------------------------------------------------
MAXRES   = spm_get_defaults('stats.maxres');
nSres    = min(nScan,MAXRES);

fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done');               %-#

%-Initialise output images (unless this is a 1st pass for ReML)
%==========================================================================
if isfield(xX,'W')
    fprintf('%-40s: %30s','Output images','...initialising');           %-#

    %-Initialise new mask name: current mask & conditions on voxels
    %----------------------------------------------------------------------
    VM    = struct('fname',  'mask.img',...
                   'dim',    DIM',...
                   'dt',     [spm_type('uint8') spm_platform('bigend')],...
                   'mat',    M,...
                   'pinfo',  [1 0 0]',...
                   'descrip','spm_spm:resultant analysis mask');
    VM    = spm_create_vol(VM);

    %-Initialise beta image files
    %----------------------------------------------------------------------
    Vbeta(1:nBeta) = deal(struct(...
        'fname',    [],...
        'dim',      DIM',...
        'dt',       [spm_type('float32') spm_platform('bigend')],...
        'mat',      M,...
        'pinfo',    [1 0 0]',...
        'descrip',  ''));

    for i = 1:nBeta
        Vbeta(i).fname   = sprintf('beta_%04d.img',i);
        Vbeta(i).descrip = sprintf('spm_spm:beta (%04d) - %s',i,xX.name{i});
    end
    Vbeta = spm_create_vol(Vbeta);

    %-Initialise residual sum of squares image file
    %----------------------------------------------------------------------
    VResMS = struct('fname',    'ResMS.img',...
        'dim',      DIM',...
        'dt',       [spm_type('float64') spm_platform('bigend')],...
        'mat',      M,...
        'pinfo',    [1 0 0]',...
        'descrip',  'spm_spm:Residual sum-of-squares');
    VResMS = spm_create_vol(VResMS);

    %-Initialise standardised residual images
    %----------------------------------------------------------------------
    VResI(1:nSres) = deal(struct(...
        'fname',    [],...
        'dim',      DIM',...
        'dt',       [spm_type('float64') spm_platform('bigend')],...
        'mat',      M,...
        'pinfo',    [1 0 0]',...
        'descrip',  'spm_spm:StandardisedResiduals'));

    for i = 1:nSres
        VResI(i).fname   = sprintf('ResI_%04d.img', i);
        VResI(i).descrip = sprintf('spm_spm:ResI (%04d)', i);
    end
    VResI = spm_create_vol(VResI);
    fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised');    %-#
end % (xX,'W')

%==========================================================================
% - F I T   M O D E L   &   W R I T E   P A R A M E T E R    I M A G E S
%==========================================================================

%-MAXMEM is the maximum amount of data processed at a time (bytes)
%--------------------------------------------------------------------------
MAXMEM = spm_get_defaults('stats.maxmem');
mmv    = MAXMEM/8/nScan;
blksz  = min(xdim*ydim,ceil(mmv));                             %-block size
nbch   = ceil(xdim*ydim/blksz);                                %-# blocks
nbz    = max(1,min(zdim,floor(mmv/(xdim*ydim))));   nbz = 1;   %-# planes
blksz  = blksz * nbz;

%-Initialise variables used in the loop
%==========================================================================
[xords, yords] = ndgrid(1:xdim, 1:ydim);
xords = xords(:)'; yords = yords(:)';           % plane X,Y coordinates
S     = 0;                                      % Volume (voxels)
s     = 0;                                      % Volume (voxels > UF)
Cy    = 0;                                      % <Y*Y'> spatially whitened
CY    = 0;                                      % <(Y - <Y>) * (Y - <Y>)'>
EY    = 0;                                      % <Y>    for ReML
i_res = round(linspace(1,nScan,nSres))';        % Indices for residual

%-Initialise XYZ matrix of in-mask voxel co-ordinates (real space)
%--------------------------------------------------------------------------
XYZ   = zeros(3,xdim*ydim*zdim);

%-Cycle over bunches blocks within planes to avoid memory problems
%==========================================================================
spm_progress_bar('Init',100,str,'');

for z = 1:nbz:zdim                       %-loop over planes (2D or 3D data)

    % current plane-specific parameters
    %----------------------------------------------------------------------
    CrPl    = z:min(z+nbz-1,zdim);       %-plane list
    zords   = CrPl(:)*ones(1,xdim*ydim); %-plane Z coordinates
    CrBl    = [];                        %-parameter estimates
    CrResI  = [];                        %-residuals
    CrResSS = [];                        %-residual sum of squares
    Q       = [];                        %-in mask indices for this plane

    for bch = 1:nbch                     %-loop over blocks

        %-Print progress information in command window
        %------------------------------------------------------------------
        if numel(CrPl) == 1
            str = sprintf('Plane %3d/%-3d, block %3d/%-3d',...
                z,zdim,bch,nbch);
        else
            str = sprintf('Planes %3d-%-3d/%-3d',z,CrPl(end),zdim);
        end
        if z == 1 && bch == 1
            str2 = '';
        else
            str2 = repmat(sprintf('\b'),1,72);
        end
        fprintf('%s%-40s: %30s',str2,str,' ');

        %-construct list of voxels in this block
        %------------------------------------------------------------------
        I     = (1:blksz) + (bch - 1)*blksz;       %-voxel indices
        I     = I(I <= numel(CrPl)*xdim*ydim);     %-truncate
        xyz   = [repmat(xords,1,numel(CrPl)); ...
                 repmat(yords,1,numel(CrPl)); ...
                 reshape(zords',1,[])];
        xyz   = xyz(:,I);                          %-voxel coordinates
        nVox  = size(xyz,2);                       %-number of voxels

        %-Get data & construct analysis mask
        %=================================================================
        fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...read & mask data')
        Cm    = true(1,nVox);                      %-current mask

        %-Compute explicit mask
        % (note that these may not have same orientations)
        %------------------------------------------------------------------
        for i = 1:length(xM.VM)

            %-Coordinates in mask image
            %--------------------------------------------------------------
            j = xM.VM(i).mat\M*[xyz;ones(1,nVox)];

            %-Load mask image within current mask & update mask
            %--------------------------------------------------------------
            Cm(Cm) = spm_get_data(xM.VM(i),j(:,Cm),false) > 0;
        end

        %-Get the data in mask, compute threshold & implicit masks
        %------------------------------------------------------------------
        Y     = zeros(nScan,nVox);
        for i = 1:nScan

            %-Load data in mask
            %--------------------------------------------------------------
            if ~any(Cm), break, end                %-Break if empty mask
            Y(i,Cm)  = spm_get_data(VY(i),xyz(:,Cm),false);

            Cm(Cm)   = Y(i,Cm) > xM.TH(i);         %-Threshold (& NaN) mask
            if xM.I && ~YNaNrep && xM.TH(i) < 0    %-Use implicit mask
                Cm(Cm) = abs(Y(i,Cm)) > eps;
            end
        end

        %-Mask out voxels where data is constant
        %------------------------------------------------------------------
        Cm(Cm) = any(diff(Y(:,Cm),1));
        Y      = Y(:,Cm);                          %-Data within mask
        CrS    = sum(Cm);                          %-# current voxels

        %==================================================================
        %-Proceed with General Linear Model (if there are voxels)
        %==================================================================
        if CrS

            %-Whiten/Weight data and remove filter confounds
            %--------------------------------------------------------------
            fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...filtering');%-#

            KWY   = spm_filter(xX.K,W*Y);

            %-General linear model: Weighted least squares estimation
            %--------------------------------------------------------------
            fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...estimation');%-#

            beta  = xX.pKX*KWY;                    %-Parameter estimates
            res   = spm_sp('r',xX.xKXs,KWY);       %-Residuals
            ResSS = sum(res.^2);                   %-Residual SSQ
            clear KWY                              %-Clear to save memory

            %-If ReML hyperparameters are needed for xVi.V
            %--------------------------------------------------------------
            if ~isfield(xVi,'V')

                %-F-threshold & accumulate spatially whitened Y*Y'
                %----------------------------------------------------------
                j   = sum((Hsqr*beta).^2,1)/trMV > UF*ResSS/trRV;
                j   = find(j);
                if ~isempty(j)
                    q  = size(j,2);
                    s  = s + q;
                    q  = spdiags(sqrt(trRV./ResSS(j)'),0,q,q);
                    Y  = Y(:,j)*q;
                    Cy = Cy + Y*Y';
                end

            end % (xVi,'V')

            %-if we are saving the WLS (ML) parameters
            %--------------------------------------------------------------
            if isfield(xX,'W')

                %-sample covariance and mean of Y (all voxels)
                %----------------------------------------------------------
                CY         = CY + Y*Y';
                EY         = EY + sum(Y,2);

                %-Save betas etc. for current plane as we go along
                %----------------------------------------------------------
                CrBl       = [CrBl,    beta];
                CrResI     = [CrResI,  res(i_res,:)];
                CrResSS    = [CrResSS, ResSS];

            end % (xX,'W')
            clear Y                         %-Clear to save memory

        end % (CrS)

        %-Append new inmask voxel locations and volumes
        %------------------------------------------------------------------
        XYZ(:,S + (1:CrS)) = xyz(:,Cm);     %-InMask XYZ voxel coords
        Q                  = [Q I(Cm)];     %-InMask XYZ voxel indices
        S                  = S + CrS;       %-Volume analysed (voxels)

    end % (bch)

    %-Plane complete, write plane to image files (unless 1st pass)
    %======================================================================
    if isfield(xX,'W')

        fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...saving plane'); %-#

        jj = NaN(xdim,ydim,numel(CrPl));

        %-Write Mask image
        %------------------------------------------------------------------
        if ~isempty(Q), jj(Q) = 1; end
        VM    = spm_write_plane(VM, ~isnan(jj), CrPl);

        %-Write beta images
        %------------------------------------------------------------------
        for i = 1:nBeta
            if ~isempty(Q), jj(Q) = CrBl(i,:); end
            Vbeta(i) = spm_write_plane(Vbeta(i), jj, CrPl);
        end

        %-Write standardised residual images
        %------------------------------------------------------------------
        for i = 1:nSres
            if ~isempty(Q), jj(Q) = CrResI(i,:)./sqrt(CrResSS/erdf); end
            VResI(i) = spm_write_plane(VResI(i), jj, CrPl);
        end

        %-Write ResSS into ResMS (variance) image scaled by tr(RV) above
        %------------------------------------------------------------------
        if ~isempty(Q), jj(Q) = CrResSS; end
        VResMS  = spm_write_plane(VResMS, jj, CrPl);

    end % (xX,'W')

    %-Report progress
    %----------------------------------------------------------------------
    fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...done');             %-#
    spm_progress_bar('Set',100*(bch + nbch*(z - 1))/(nbch*zdim));

end % (for z = 1:zdim)
fprintf('\n');                                                          %-#
spm_progress_bar('Clear')

%==========================================================================
% - P O S T   E S T I M A T I O N   C L E A N U P
%==========================================================================
if S == 0, spm('alert!','No inmask voxels - empty analysis!'); return; end

%-average sample covariance and mean of Y (over voxels)
%--------------------------------------------------------------------------
CY   = CY/S;
EY   = EY/S;
CY   = CY - EY*EY';

%-If not defined, compute non-sphericity V using ReML Hyperparameters
%==========================================================================
if ~isfield(xVi,'V')

    %-check there are signficant voxels
    %----------------------------------------------------------------------
    if s == 0
        spm('FigName','Stats: no significant voxels',Finter);
        spm('Pointer','Arrow');
        if isfield(SPM.xGX,'rg')&&~isempty(SPM.xGX.rg)
            figure(Finter);
            plot(SPM.xGX.rg);
            spm('alert*',{'Please check your data'; ...
                'There are no significant voxels';...
                'The globals are plotted for diagnosis'});
        else
            spm('alert*',{'Please check your data'; ...
                'There are no significant voxels'});
        end
        warning('Please check your data: There are no significant voxels.');
        return
    end

    %-ReML estimate of residual correlations through hyperparameters (h)
    %----------------------------------------------------------------------
    str    = 'Temporal non-sphericity (over voxels)';
    fprintf('%-40s: %30s\n',str,'...ReML estimation');                  %-#
    Cy     = Cy/s;

    % ReML for separable designs and covariance components
    %----------------------------------------------------------------------
    if isstruct(xX.K)
        m     = length(xVi.Vi);
        h     = zeros(m,1);
        V     = sparse(nScan,nScan);
        for i = 1:length(xX.K)

            % extract blocks from bases
            %--------------------------------------------------------------
            q     = xX.K(i).row;
            p     = [];
            Qp    = {};
            for j = 1:m
                if nnz(xVi.Vi{j}(q,q))
                    Qp{end + 1} = xVi.Vi{j}(q,q);
                    p           = [p j];
                end
            end

            % design space for ReML (with confounds in filter)
            %--------------------------------------------------------------
            Xp     = xX.X(q,:);
            try
                Xp = [Xp xX.K(i).X0];
            end

            % ReML
            %--------------------------------------------------------------
            fprintf('%-30s\n',sprintf('  ReML Block %i',i));
            [Vp,hp] = spm_reml(Cy(q,q),Xp,Qp);
            V(q,q)  = V(q,q) + Vp;
            h(p)    = hp;
        end
    else
        [V,h] = spm_reml(Cy,xX.X,xVi.Vi);
    end

    % normalize non-sphericity and save hyperparameters
    %----------------------------------------------------------------------
    V         = V*nScan/trace(V);
    xVi.h     = h;
    xVi.V     = V;                  % Save non-sphericity xVi.V
    xVi.Cy    = Cy;                 % spatially whitened <Y*Y'>
    SPM.xVi   = xVi;                % non-sphericity structure

    % If xX.W is not specified use W*W' = inv(V) to give ML estimators
    %----------------------------------------------------------------------
    if ~isfield(xX,'W')
        if spm_matlab_version_chk('7') >=0
            save('SPM','SPM','-V6');
        else
            save('SPM','SPM');
        end
        clear
        load SPM
        SPM = spm_spm(SPM);
        return
    end
end

%-Use non-sphericity xVi.V to compute [effective] degrees of freedom
%==========================================================================
xX.V            = spm_filter(xX.K,spm_filter(xX.K,W*V*W')');% KWVW'K'
[trRV trRVRV]   = spm_SpUtil('trRV',xX.xKXs,xX.V);          % trRV (for X)
xX.trRV         = trRV;                                     % <R'*y'*y*R>
xX.trRVRV       = trRVRV;                                   %-Satterthwaite
xX.erdf         = trRV^2/trRVRV;                            % approximation
xX.Bcov         = xX.pKX*xX.V*xX.pKX';                      % Cov(beta)

%-Set VResMS scalefactor as 1/trRV (raw voxel data is ResSS)
%--------------------------------------------------------------------------
VResMS.pinfo(1) = 1/xX.trRV;
VResMS          = spm_create_vol(VResMS);

%-Smoothness estimates of component fields and RESEL counts for volume
%==========================================================================
try
    FWHM = SPM.xVol.FWHM;
    VRpv = SPM.xVol.VRpv;
    R    = SPM.xVol.R;
catch
    [FWHM,VRpv,R] = spm_est_smoothness(VResI,VM,[nScan erdf]);
end

%-Delete the residuals images
%==========================================================================
% j = spm_select('List',SPM.swd,'^ResI_.{4}\..{3}$');
% for  k = 1:size(j,1)
%     spm_unlink(deblank(j(k,:)));
% end

%-Compute scaled design matrix for display purposes
%--------------------------------------------------------------------------
xX.nKX        = spm_DesMtx('sca',xX.xKXs.X,xX.name);

%-Save remaining results files and analysis parameters
%==========================================================================
fprintf('%-40s: %30s','Saving results','...writing');                   %-#

%-place fields in SPM
%--------------------------------------------------------------------------
SPM.xVol.XYZ   = XYZ(:,1:S);        %-InMask XYZ coords (voxels)
SPM.xVol.M     = M;                 %-voxels -> mm
SPM.xVol.iM    = inv(M);            %-mm -> voxels
SPM.xVol.DIM   = DIM;               %-image dimensions
SPM.xVol.FWHM  = FWHM;              %-Smoothness data
SPM.xVol.R     = R;                 %-Resel counts
SPM.xVol.S     = S;                 %-Volume (voxels)
SPM.xVol.VRpv  = VRpv;              %-Filehandle - Resels per voxel

SPM.Vbeta      = Vbeta;             %-Filehandle - Beta
SPM.VResMS     = VResMS;            %-Filehandle - Hyperparameter
SPM.VM         = VM;                %-Filehandle - Mask

SPM.xVi        = xVi;               % non-sphericity structure
SPM.xVi.CY     = CY;                %-<(Y - <Y>)*(Y - <Y>)'>

SPM.xX         = xX;                %-design structure

SPM.xM         = xM;                %-mask structure

SPM.xCon       = struct([]);        %-contrast structure

SPM.SPMid      = SPMid;
SPM.swd        = pwd;

%-Save analysis parameters in SPM.mat file
%--------------------------------------------------------------------------
if spm_matlab_version_chk('7') >=0
    save('SPM','SPM','-V6');
else
    save('SPM','SPM');
end

%==========================================================================
%- E N D: Cleanup GUI
%==========================================================================
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')                %-#
spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
fprintf('%-40s: %30s\n','Completed',spm('time'))                        %-#
fprintf('...use the results section for assessment\n\n')                %-#

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function experiments = stack_vot2019_rgbd()

set_global_variable('bundle', 'http://data.votchallenge.net/vot2019/rgbd/description.json');

unsupervised.name = 'rgbd-unsupervised';
unsupervised.converter = [];
unsupervised.type = 'unsupervised';
unsupervised.tags = {};

unsupervised.parameters.repetitions = 1;
unsupervised.parameters.burnin = 0;

unsupervised.analysis = {'overlap'};

experiments = {unsupervised};

end

author: votchallenge

repo: toolkit-legacy


text:
function f=peven(f,dim)
%PEVEN   Even part of periodic function
%   Usage:  fe=peven(f);
%           fe=peven(f,dim);
%
%   `peven(f)` returns the even part of the periodic sequence *f*.
%
%   `peven(f,dim)` does the same along dimension *dim*.
%
%   See also:  podd, dft, involute, pconv

if nargin==1
  f=(f+involute(f))/2;
else
  f=(f+involute(f,dim))/2;
end;


text:
classdef nnpdist < nntest
  properties (TestParameter)
    oneToOne = {false, true}
    noRoot = {false, true}
    p = {.5 1 2 3}
    aggregate = {false, true}
  end
  methods (Test)
    function basic(test,oneToOne, noRoot, p, aggregate)
      if aggregate
        % make it smaller to avoid numerical derivative issues with
        % float
        h = 3 ;
        w = 2 ;
      else
        h = 13 ;
        w = 17 ;
      end
      d = 4 ;
      n = 5 ;
      x = test.randn(h,w,d,n) ;
      if oneToOne
        x0 = test.randn(h,w,d,n) ;
      else
        x0 = test.randn(1,1,d,n) ;
      end
      opts = {'noRoot', noRoot, 'aggregate', aggregate} ;

      y = vl_nnpdist(x, x0, p, opts{:}) ;

      % make sure they are not too close in any dimension as this may be a
      % problem for the finite difference dereivatives as one could
      % approach 0 which is not differentiable for some p-norms

      s = abs(bsxfun(@minus, x, x0)) < test.range*1e-1 ;
      x(s) = x(s) + 5*test.range ;

      dzdy = test.rand(size(y)) ;
      dzdx = vl_nnpdist(x,x0,p,dzdy,opts{:}) ;
      test.der(@(x) vl_nnpdist(x,x0,p,opts{:}), x, dzdy, dzdx, test.range * 1e-3) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval