text:
function hf=anim8(varargin)

% function hf=anim8(hf,animStruct)
% ------------------------------------------------------------------------
%
% Change log:
% 2019/05/20 Added loading of saved anim8 figure by using the path to the
% figure as sole input (or no input which triggers uigetfile)
% 2019/08/09 Changed to use uicontrol slider rather than java slider due
% to future removal of javacomponent
% 2019/10/13 Fixed bug in relation to handling a single time step (no
% annimation), fixed by copying state twice and animating clones states.
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case {0,1} %Load previous
        if nargin==1 %If path to figure is specified
            figureFullPath=varargin{1};
        else %If no path is specified
            [figureFileName,figurePath]=uigetfile;
            figureFullPath=fullfile(figurePath,figureFileName);
        end
        hf=open(figureFullPath); %Open figure and keep handle

        %Add view control widget
        hp=vcw(hf);
        hf.UserData.cFigure.Handles.vcw=hp;

        %Add export figure widget
        efw(hf);

        if ~isfield(hf.UserData,'anim8')
            error('Loaded figure lacks an anim8 structure as UserData');
        else
            animStruct=hf.UserData.anim8.animStruct; %Get anim8 structure
        end

        try %Try to remove a slider if present
            hSlider=hf.UserData.anim8.sliderHandles{1};
            if isa(hSlider,'matlab.ui.control.UIControl')
                delete(hSlider)
            else %Old jSlider type
                try %javacomponent is scheduled to be removed but try it here for old slider bar type
                    warning off %Disable all warnings
                    [~,hc]=javacomponent(hSlider);
                    delete(hc);
                    warning on %Enable warnings again
                catch
                    %No fix if the above throws an error
                end
            end
        catch %Nothing to delete
        end

        try %Try to remove text label above slider if present
            delete(hf.UserData.anim8.ButtonHandles.hTextTime)
        catch %Nothing to delete
        end

    case 2 %Create new
        hf=varargin{1}; %Figure handle
        animStruct=varargin{2}; %The anim8 structure
end

if ~isfield(animStruct,'Time')
    animStruct.Time=linspace(0,numel(animStruct.Handles),numel(animStruct.Handles));
end

%% Visualisation settings

% fontColor='w';
fontSize=15;
% cMap=gjet(250);
scrollBarWidth=30;

%% Defining slider
figure(hf); drawnow;

if numel(animStruct.Time)==1
    warning('Only 1 step for animation');
    animStruct.Time(end+1)=animStruct.Time(end);
    animStruct.Handles{end+1}=animStruct.Handles{end}; %Handles of objects to animate
    animStruct.Props{end+1}=animStruct.Props{end}; %Properties of objects to animate
    animStruct.Set{end+1}=animStruct.Set{end};
end

animTime=animStruct.Time(:);
sliceIndexI=numel(animTime); %Initial index at end
sliderStep=[1/(numel(animTime)-1) 1/(numel(animTime)-1)]; %Slider step sizes

%Initialize slider
hSlider= uicontrol(hf,'Style','slider','Position',[0,0,round(hf.Position(3)),scrollBarWidth]);
set(hSlider,'Value',sliceIndexI,'Min',1,'Max',numel(animTime),'SliderStep',sliderStep);
hSlider.Callback={@updateViewFunc,hf};
% hSlider.KeyPressFcn={@updateViewFunc,hf};
addlistener(hSlider,'ContinuousValueChange',@(hObject, event) updateViewFunc(hObject, event,hf));
addlistener(hSlider,'Value','PostSet',@(hObject, event) updateViewFunc(hObject, event,hf));

%% Set resize function

% set(hf,'ResizeFcn',{@figResize,{hf,scrollBarWidth,hSlider}});
% set(hf,'ResizeFcn',@(h,e)figResize(h,e,{hf,scrollBarWidth,hSlider}));

hFunc=get(hf,'ResizeFcn');

if iscell(hFunc)
    warning('anim8 changed the ResizeFcn function for this figure. Specify your ResizeFcn in the form @(h,e)figResize(h,e,c) to avoid this behavior.');
    set(hf,'ResizeFcn',@(a,b)figResize(a,b,{hf,scrollBarWidth,hSlider}));
else
    if isempty(hFunc)
        set(hf,'ResizeFcn',@(a,b)figResize(a,b,{hf,scrollBarWidth,hSlider}));
    else
        set(hf,'ResizeFcn',@(a,b)(cellfun(@(x)feval(x,a,b),{hFunc,@(a,b)figResize(a,b,{hf,scrollBarWidth,hSlider})})));
    end
end

%% Initialize figure callbacks

set(hf,'KeyPressFcn', {@figKeyPressFunc,{hf}});%'WindowButtonDownFcn', {@figMouseDown,{hf}},'WindowButtonUpFcn', {@mouseup,hf});

%% Initialise tool bar and buttons

filePath=mfilename('fullpath');
toolboxPath=fileparts(fileparts(filePath));
iconPath=fullfile(toolboxPath,'icons');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Toolbar
toolbarTag='anim8_toolbar';

%Remove existing anim8 toolbar
hb = findall(hf,'Tag',toolbarTag); % hb = findall(hf,'Type','uitoolbar');
if ~isempty(hb)
    delete(hb)
else
    %Find all toolbars
    hb = findall(hf,'Type','uitoolbar');

    %Remove them if they are not the default bar
    if ~isempty(hb)
        for q=1:1:numel(hb)
            if ~strcmp(hb(q).Tag,'FigureToolBar')
                delete(hb(q));
            end
        end
    end
end

%Create new anim8 toolbar
hb = uitoolbar(hf);
hb.Tag=toolbarTag;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Help button

%get icon
D=importdata(fullfile(iconPath,'help.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
S(S==1)=NaN;
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end
% Create a uipushtool in the toolbar
uipushtool(hb,'TooltipString','Help','CData',S,'Tag','help_button','ClickedCallback',@helpFunc,'Separator','on');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Play button

%get icon 1
D=importdata(fullfile(iconPath,'play.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
S(S==1)=NaN;
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end
iconPlay=S;

%get icon 2
D=importdata(fullfile(iconPath,'stop.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
S(S==1)=NaN;
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end
iconStop=S;

% Create a uitoggletool in the toolbar
hPlay=uitoggletool(hb,'TooltipString','Play','CData',iconPlay,'Tag','play_button','ClickedCallback',{@playFunc,{hf}});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Time button

%get icon
D=importdata(fullfile(iconPath,'time.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
S(S==1)=NaN;
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end
% Create a uipushtool in the toolbar
uipushtool(hb,'TooltipString','Time','CData',S,'Tag','time_button','ClickedCallback',{@timeFunc,{hf}}); %,'Separator','on'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Cycle button

%get icon
D=importdata(fullfile(iconPath,'cycle.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
S(S==1)=NaN;
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end
% Create a uitoggletool in the toolbar
hCycle=uitoggletool(hb,'TooltipString','Cycle forward-backward','CData',S,'Tag','cycle_button');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Save animation button

%get icon
D=importdata(fullfile(iconPath,'saveAnimation.jpg'));
S=double(D);
S=S-min(S(:));
S=S./max(S(:));
S(S==1)=NaN;
if size(S,3)==1
    S=repmat(S,[1 1 3]);
end
% Create a uipushtool in the toolbar
uipushtool(hb,'TooltipString','Save .gif animation','CData',S,'Tag','saveAnimation_button','ClickedCallback',{@saveAnimationFunc,{hf}}); %,'Separator','on'

%% Text fields

% title(sprintf('%6.16e ',T),'FontSize',hf.UserData.anim8.fontSize);
% title(sprintf('%f',T),'FontSize',hf.UserData.anim8.fontSize);

hTextTime = uicontrol(hf,'Style','text',...
    'String',[' Time: ',sprintf('%f',animTime(1))],...
    'Position',[0 scrollBarWidth round(hf.Position(3)) scrollBarWidth],...
    'BackgroundColor',[1 1 1],'HorizontalAlignment','Left',...
    'fontSize',fontSize);

%% Set figure UserData

hf.UserData.anim8.sliderHandles={hSlider};
% hf.UserData.anim8.ButtonHandles.Sample=hSample;
hf.UserData.anim8.animStruct=animStruct;
hf.UserData.anim8.fontSize=fontSize;
% animHandles=animStruct.Handles; %Handles of objects to animate
% animProps=animStruct.Props; %Properties of objects to animate
% animSets=animStruct.Set; %Property values for to set in order to animate
% animTime=animStruct.Time;
hf.UserData.anim8.playDir=1;
hf.UserData.anim8.ButtonHandles.Play=hPlay;
hf.UserData.anim8.ButtonHandles.hCycle=hCycle;
hf.UserData.anim8.ButtonHandles.hTextTime=hTextTime;
hf.UserData.anim8.pauseTime=1/numel(animStruct.Time);%mean(diff(animStruct.Time));
hf.UserData.anim8.shiftMag=1;%ceil(numel(animStruct.Time)/20);
hf.UserData.anim8.icons.play=iconPlay;
hf.UserData.anim8.icons.stop=iconStop;

% Store current settings
hf.UserData.anim8.WindowButtonDownFcn=hf.WindowButtonDownFcn;
hf.UserData.anim8.WindowButtonUpFcn=hf.WindowButtonUpFcn;
hf.UserData.anim8.KeyPressFcn=hf.KeyPressFcn;
hf.UserData.anim8.WindowScrollWheelFcn=hf.WindowScrollWheelFcn;
hf.UserData.anim8.BusyAction=hf.BusyAction;

%Export figure widget settings
hf.UserData.efw.defaultPath=fullfile(cd,'efw');
hf.UserData.efw.imName=['figure',num2str(get(hf,'Number'))];
hf.UserData.efw.imExt='png';
hf.UserData.efw.imRes='50';
hf.UserData.efw.exportFigOpt='-nocrop';
hf.UserData.efw.exportGifOpt='1';

%% Initialize slider locations
set(hSlider,'Value',round(sliceIndexI/2));
set(hSlider,'Value',sliceIndexI);
drawnow;

end

%% Scroll bar resizing

function figResize(~,~,inputCell)
hf=inputCell{1};
scrollBarWidth=inputCell{2};

hSlider=inputCell{3};
if isa(hSlider,'matlab.ui.control.UIControl') %Check if its a MATLAB style slider
    set(hSlider,'Position',[0,0,round(hf.Position(3)), scrollBarWidth]);
    set(hf.UserData.anim8.ButtonHandles.hTextTime,'Position',[0 scrollBarWidth round(hf.Position(3)) scrollBarWidth]);
else %Old JAVA jSlider
%     try
% %         delete(hSlider)
%     catch
%     end
%     try %javacomponent is scheduled to be removed but try it here for old slider bar type
%         jc=javacomponent(hSlider,[0,0,round(hf.Position(3)), scrollBarWidth]);
%         jc.delete
%     catch
%         %No fix if the above throws an error
%     end
end

end

%% Help

function helpFunc(~,~)

msgText={'* Play button or press the space bar -> Start animation',...
    '* Stop button or press the space bar -> Stop animation',...
    '* Time button -> Change time stepping settings',...
    '* Cycle button -> Toggle forward / forward-backward mode',...
    '* Save button -> Export snap shorts (e.g. jpg) and animated .gif files',...
    '* Press the v key to activate the view control widget (vcw), note that vcw changes key press functions until vcw is deactivated',...
    };
helpdlg(msgText,'Help information');

end

%% play

function playFunc(~,~,inputCell)
hf=inputCell{1};
shiftMag=hf.UserData.anim8.shiftMag;

set(hf.UserData.anim8.ButtonHandles.Play,'CData',hf.UserData.anim8.icons.stop,'TooltipString','Stop');

while strcmp(get(hf.UserData.anim8.ButtonHandles.Play,'State'),'on')
    tic
    if strcmp(get(hf.UserData.anim8.ButtonHandles.hCycle,'State'),'on')
        hSlider=hf.UserData.anim8.sliderHandles{1};

        sliderValue=get(hSlider,'Value');
        sliderValueNew=sliderValue+(shiftMag*hf.UserData.anim8.playDir);

        sliderMax=get(hSlider,'Max');
        sliderMin=get(hSlider,'Min');

        if sliderValueNew<sliderMin
            hf.UserData.anim8.playDir=hf.UserData.anim8.playDir*-1;
            sliderValueNew=sliderValue+(shiftMag*hf.UserData.anim8.playDir);
        elseif sliderValueNew>sliderMax
            hf.UserData.anim8.playDir=hf.UserData.anim8.playDir*-1;
            sliderValueNew=sliderValue+(shiftMag*hf.UserData.anim8.playDir);
        end
        set(hSlider,'Value',sliderValueNew);
    else
        shiftSlider(hf.UserData.anim8.sliderHandles{1},shiftMag);
    end
    updateViewFunc([],[],hf);
    drawnow;

    t=toc;
    tPause=hf.UserData.anim8.pauseTime-t;
    if tPause>0
        pause(tPause);
    end
end

set(hf.UserData.anim8.ButtonHandles.Play,'CData',hf.UserData.anim8.icons.play,'TooltipString','Play');

end

%% Time
function timeFunc(~,~,inputCell)
hf=inputCell{1};

prompt = {'Enter pause time:','Enter step size:'};
dlg_title = 'Time stepping settings';
defaultOptions = {num2str(hf.UserData.anim8.pauseTime),num2str(hf.UserData.anim8.shiftMag)};
s=40+max([cellfun(@numel,prompt) cellfun(@numel,defaultOptions)]);
Q = inputdlg(prompt,dlg_title,[1 s],defaultOptions);
if ~isempty(Q)
    if ~isempty(Q{1})
        hf.UserData.anim8.pauseTime=str2double(Q{1});
    end
    if ~isempty(Q{2})
        hf.UserData.anim8.shiftMag=str2double(Q{2});
    end
end

end

%% Save Animation

function saveAnimationFunc(~,~,inputCell)

hf=inputCell{1};
defStruct=hf.UserData.efw; %Default option structure
exportGifAnim8(hf,defStruct,1); %Export gif data

end

%% Figure key press

function figKeyPressFunc(~,eventData,inputCell)

hf=inputCell{1}; %Figure handle

% Key input options
switch eventData.Key
    case {'leftarrow','downarrow'}
        shiftSlider(hf.UserData.anim8.sliderHandles{1},-1*hf.UserData.anim8.shiftMag);
        updateViewFunc([],[],hf);
    case {'rightarrow','uparrow'}
        shiftSlider(hf.UserData.anim8.sliderHandles{1},1*hf.UserData.anim8.shiftMag);
        updateViewFunc([],[],hf);
    case 'home' % Go to the start
        hSlider=hf.UserData.anim8.sliderHandles{1};
        sliderVal=get(hSlider,'Min');
        set(hSlider,'Value',sliderVal);
        updateViewFunc([],[],hf);
    case 'end' % Go to the end
        hSlider=hf.UserData.anim8.sliderHandles{1};
        sliderVal=get(hSlider,'Max');
        set(hSlider,'Value',sliderVal);
        updateViewFunc([],[],hf);
    case 'v' % Activate vcw
        set(hf.UserData.cFigure.Handles.vcw,'State','On');
    case'space'
        if strcmp(get(hf.UserData.anim8.ButtonHandles.Play,'State'),'off')
            set(hf.UserData.anim8.ButtonHandles.Play,'State','on');
            playFunc([],[],{hf});
        else
            set(hf.UserData.anim8.ButtonHandles.Play,'State','off');
        end
end

end

%% updateViewFunc

function updateViewFunc(~,~,hf)

animStruct=hf.UserData.anim8.animStruct;
hSlider=hf.UserData.anim8.sliderHandles{1};
sliderValue=round(get(hSlider,'Value'));
set(hSlider,'Value',sliderValue);

T=animStruct.Time(sliderValue);
H=animStruct.Handles{sliderValue};%e.g. [hp,hp]; %Handles of objects to animate
P=animStruct.Props{sliderValue};% e.g. {'FaceColor','Vertices'}; %Properties of objects to animate
S=animStruct.Set{sliderValue};%e.g.{c(q,:),V}; %Property values for to set in order to animate

set(hf.UserData.anim8.ButtonHandles.hTextTime,'String',[' Time: ',sprintf('%f',T)],'FontSize',hf.UserData.anim8.fontSize);

for q=1:1:numel(H)
    h=H(q); % Current graphics handle
    p=P{q}; % Current graphics property
    s=S{q}; % Current property setting
    h.(p)=s;% Setting the property
end

end

%% Shift slider

function shiftSlider(hSlider,shiftMag)

sliderValue=get(hSlider,'Value');
sliderValueNew=sliderValue+shiftMag;

sliderMax=get(hSlider,'Max');
sliderMin=get(hSlider,'Min');

if sliderValueNew<sliderMin
    set(hSlider,'Value',sliderMax);
elseif sliderValueNew>sliderMax
    set(hSlider,'Value',1);
else
    set(hSlider,'Value',sliderValueNew);
end

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function setup_vocpart()
%% download images
if(~exist('../data/VOC_part','dir'))
    system('wget -O ../VOCPart.tar.gz --no-check-certificate http://www.stat.ucla.edu/~xianjie.chen/pascal_part_dataset/trainval.tar.gz');
    mkdir('../data/VOC_part/')
    system('tar -xvzf ../VOCPart.tar.gz -C ../data/VOC_part/');
    system('rm ../VOCPart.tar.gz');
    system('wget -O ../img.tar http://host.robots.ox.ac.uk/pascal/VOC/voc2010/VOCtrainval_03-May-2010.tar');
    system('tar -xvf ../img.tar -C ../data/VOC_part/');
    system('rm ../img.tar');
end

%% settings
conf.data.catedir='../data/VOC_part/';
conf.data.imgdir='../data/VOC_part/VOCdevkit/VOC2010/JPEGImages/';
conf.data.readCode='./data_input/data_input_VOC/';

conf.data.minArea=50^2;

conf.output.dir='./mat/';

Name_batch={'bird','cat','cow','dog','horse','sheep'};
for i=1:numel(Name_batch)
    mkdir([conf.output.dir,Name_batch{i}]);
    save([conf.output.dir,Name_batch{i},'/conf.mat'],'conf');
end

addpath(genpath('./tool'));
addpath(conf.data.readCode);

repo: interpretableCNN


text:
function K = sdlfmXsdlfmKernCompute(sdlfmKern1, sdlfmKern2, t1, t2, covIC, type)

% SDLFMXSDLFMKERNCOMPUTE Compute a cross kernel between two SDLFM kernels.
% FORMAT
% DESC computes cross kernel terms between two switching dynamical
% LFM kernels for the multiple output kernel.
% ARG sdlfmKern1 : the kernel structure associated with the first SDLFM
% kernel.
% ARG sdlfmKern2 : the kernel structure associated with the second SDLFM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two SDLFM kernels for
% the multiple output kernel.
% ARG sdlfmKern1 : the kernel structure associated with the first SDLFM
% kernel.
% ARG sdlfmKern2 : the kernel structure associated with the second SDLFM
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two SDLFM kernels for
% the multiple output kernel. The SDLFM kernels can correspond to Position
% X Position (default), Velocity X Position, Velocity X Velocity,
% Acceleration X Position, Acceleration X Velocity, Acceleration X
% Acceleration. The type of kernel to be computed is specified in 'type'.
% ARG sdlfmKern1 : the kernel structure associated with the first SDLFM
% kernel.
% ARG sdlfmKern2 : the kernel structure associated with the second SDLFM
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% ARG type : specifies the type of kerne to be computed
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : sdlfmKernParamInit, sdlfmKernCompute, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 6
    type = 'PosPos';
    if nargin < 5
        covIC = t2;
        t2 = t1;
    end
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end

compareInverseWidth = sdlfmKern1.inverseWidth == sdlfmKern2.inverseWidth;
if sum(sum(compareInverseWidth))~=(sdlfmKern1.nIntervals*sdlfmKern1.nlfPerInt)
    error('Kernels cannot be cross combined if they have different inverse widths.')
end
compareSwitchingTimes = sdlfmKern1.switchingTimes == sdlfmKern2.switchingTimes;
if sum(sum(compareSwitchingTimes))~=sdlfmKern1.nIntervals
    error('Kernels cannot be cross combined if they have different switching points.')
end

switch type
    case 'PosPos'
        fhandle = 'sdlfmXsdlfmKernComputeBlock';
    case 'VelPos'
        fhandle = 'sdlfmvXsdlfmKernComputeBlock';
    case 'VelVel'
        fhandle = 'sdlfmvXsdlfmvKernComputeBlock';
    case 'AccelPos'
        fhandle = 'sdlfmaXsdlfmKernComputeBlock';
    case 'AccelVel'
        fhandle = 'sdlfmaXsdlfmvKernComputeBlock';
    case 'AccelAccel'
        fhandle = 'sdlfmaXsdlfmaKernComputeBlock';
end

fhandle = str2func(fhandle);

%Form the basic kernels
lfmKern1 = struct();
lfmKern2 = struct();

% Create structures that will make easy the computation of the kernels

spVector = [cumsum(sdlfmKern1.switchingTimes) t1(end)+50];

dim1 = zeros(1, sdlfmKern1.nIntervals); dim2 = zeros(1, sdlfmKern1.nIntervals);

for i=1:sdlfmKern1.nIntervals
    for j =1:sdlfmKern1.nlfPerInt
        % Create the appropriate set of kernel structures
        lfmKern1(i,j).mass = sdlfmKern1.mass;
        lfmKern1(i,j).spring = sdlfmKern1.spring;
        lfmKern1(i,j).damper = sdlfmKern1.damper;
        lfmKern1(i,j).inverseWidth = sdlfmKern1.inverseWidth(j,i);
        lfmKern1(i,j).sensitivity = sdlfmKern1.sensitivity(j,i);
        lfmKern2(i,j).mass = sdlfmKern2.mass;
        lfmKern2(i,j).spring = sdlfmKern2.spring;
        lfmKern2(i,j).damper = sdlfmKern2.damper;
        lfmKern2(i,j).inverseWidth = sdlfmKern2.inverseWidth(j,i);
        lfmKern2(i,j).sensitivity = sdlfmKern2.sensitivity(j,i);
        lfmKern1(i,j).limit = spVector(i+1) - spVector(i);
        lfmKern2(i,j).limit = spVector(i+1) - spVector(i);
        lfmKern1(i,j).isNormalised = sdlfmKern1.isNormalised;
        lfmKern2(i,j).isNormalised = sdlfmKern2.isNormalised;
    end
    newt1 = t1(t1> spVector(i) & t1<spVector(i+1));
    newt2 = t2(t2> spVector(i) & t2<spVector(i+1));
    dim1(i) = length(newt1);
    dim2(i) = length(newt2);
end

if sum(dim1)~=length(t1) || sum(dim2)~=length(t2)
    error('A problem with the dimensions of the switching intervals occured')
end

% Compute some necessary constants

generalConst = sdlfmKernComputeConstant(sdlfmKern1.nIntervals, ...
    lfmKern1(1), lfmKern2(1), spVector);

kyy = zeros(sdlfmKern1.nIntervals);kyv = zeros(sdlfmKern1.nIntervals);
kvy = zeros(sdlfmKern1.nIntervals);kvv = zeros(sdlfmKern1.nIntervals);
kyy(1,1) = covIC(1,1);
kvy(1,1) = covIC(2,1);
kyv(1,1) = covIC(1,2);
kvv(1,1) = covIC(2,2);

%%% Compute initial conditions for intervals (1-2), (2-2) and (2-1)

tInit1 = [spVector(1) - spVector(1);spVector(2) - spVector(1)];
tInit2 = [spVector(1) - spVector(1);spVector(2) - spVector(1)];

% Pos -- Pos
temp = sdlfmXsdlfmKernComputeBlock(lfmKern1(1), lfmKern2(1), ...
    tInit1, tInit2, kyy(1,1), kyv(1,1), kvy(1,1), kvv(1,1), 1, 1, generalConst);
kyy(2,2) = temp(2,2); kyy(1,2) = temp(1,2); kyy(2,1)  = temp(2,1);
% Vel -- Pos
temp = sdlfmvXsdlfmKernComputeBlock(lfmKern1(1), lfmKern2(1), ...
    tInit1, tInit2, kyy(1,1), kyv(1,1), kvy(1,1), kvv(1,1), 1, 1, generalConst);
kvy(2,2) = temp(2,2); kvy(1,2) = temp(1,2); kvy(2,1)  = temp(2,1);
% Pos -- Vel
temp = sdlfmXsdlfmvKernComputeBlock(lfmKern1(1), lfmKern2(1), ...
    tInit1, tInit2, kyy(1,1), kyv(1,1), kvy(1,1), kvv(1,1), 1, 1, generalConst);
kyv(2,2) = temp(2,2); kyv(1,2) = temp(1,2); kyv(2,1)  = temp(2,1);
% Vel -- Vel
temp = sdlfmvXsdlfmvKernComputeBlock(lfmKern1(1), lfmKern2(1), ...
    tInit1, tInit2, kyy(1,1), kyv(1,1), kvy(1,1), kvv(1,1), 1, 1, generalConst);
kvv(2,2) = temp(2,2); kvv(1,2) = temp(1,2); kvv(2,1)  = temp(2,1);

tempPosPos = cell(sdlfmKern1.nIntervals);
tempVelPos = cell(sdlfmKern1.nIntervals);
tempPosVel = cell(sdlfmKern1.nIntervals);
tempVelVel = cell(sdlfmKern1.nIntervals);

startValOne = 1;
endValOne   = 0;

for i=1:sdlfmKern1.nIntervals
    endValOne = endValOne + dim1(i);
    startValThree = 1;
    endValThree = 0;
    for j=1:sdlfmKern1.nIntervals
        if i>j
            lfmKern1Local = lfmKern1(j,:);
            lfmKern2Local = lfmKern2(j,:);
        else
            lfmKern1Local = lfmKern1(i,:);
            lfmKern2Local = lfmKern2(i,:);
        end
        endValThree = endValThree + dim2(j);
        % POS -- POS (Kernel and initial positions)
        K(startValOne:endValOne, startValThree:endValThree) = fhandle(lfmKern1Local, ...
            lfmKern2Local, t1(startValOne:endValOne) - spVector(i), t2(startValThree:endValThree) - spVector(j), ...
            kyy(i,j), kyv(i,j), kvy(i,j), kvv(i,j), i, j, generalConst);
        % Time vector initial conditions
        tInit1 = [spVector(i) - spVector(i);spVector(i+1) - spVector(i)];
        tInit2 = [spVector(j) - spVector(j);spVector(j+1) - spVector(j)];
        tempPosPos{i,j} = sdlfmXsdlfmKernComputeBlock(lfmKern1Local, ...
            lfmKern2Local, tInit1, tInit2, ...
            kyy(i,j), kyv(i,j), kvy(i,j), kvv(i,j), i, j, generalConst);
        kyy = organizeIC(kyy, tempPosPos, i, j);
        % VEL -- POS
        tempVelPos{i,j} = sdlfmvXsdlfmKernComputeBlock(lfmKern1Local, ...
            lfmKern2Local, tInit1, tInit2 , ...
            kyy(i,j), kyv(i,j), kvy(i,j), kvv(i,j), i, j, generalConst);
        kvy = organizeIC(kvy, tempVelPos, i, j);
        % POS -- VEL
        tempPosVel{i,j} = sdlfmXsdlfmvKernComputeBlock(lfmKern1Local, ...
            lfmKern2Local, tInit1, tInit2, ...
            kyy(i,j), kyv(i,j), kvy(i,j), kvv(i,j), i, j, generalConst);
        kyv = organizeIC(kyv, tempPosVel, i, j);
        % POS -- VEL
        tempVelVel{i,j} = sdlfmvXsdlfmvKernComputeBlock(lfmKern1Local, ...
            lfmKern2Local, tInit1, tInit2, ...
            kyy(i,j), kyv(i,j), kvy(i,j), kvv(i,j), i, j, generalConst);
        kvv = organizeIC(kvv, tempVelVel, i, j);
        startValThree = endValThree + 1;
    end
    startValOne = endValOne + 1;
end

K = real(K);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% clc;
clear;
warning off
addpath('support_methods/');
addpath(genpath('tools'));
vl_setupnn;

%% general parameters
% 4 supported datasets: ChallengeDB_release LIVE TID2013 CSIQ
testDatabase = 'ChallengeDB_release';
% 3 supported CNN architectures: ResNet AlexNet S_CNN
ModelType = 'S_CNN';
repetitions = 10;
trainingPropertion = 0.8;
%% parameters for PQR
quantizationMethod = 'uniform'; % LloydMax
beta = 64;
bins = 5;  % set bins=1 for scalar quality score regression

switch ModelType
    case 'AlexNet'
        patchSize = 227;
        patchStep = 16;
        epoch = 20;
    case 'ResNet'
        patchSize = 224;
        patchStep = 16;
        epoch = 10;
    case 'S_CNN'
        patchSize = 64;
        patchStep = 8;
        epoch = 40;
end

for patchNum = [50]
    if strcmp(ModelType,'S_CNN')
        patchNum = patchNum * 10;
    end

%% training
for seed = 1:repetitions
    TrainModel('Model',ModelType,'Database',testDatabase,'PatchNum', patchNum,...
               'seed',seed, 'trainingpropertion', trainingPropertion,...
               'quantizationMethod',quantizationMethod,...
               'bins',bins, 'beta', beta,'patchSize',patchSize,...
               'patchStep',patchStep,'epoch',epoch);
end

%% testing
for seed = 1:repetitions
    fprintf('seed = %d....\n', seed);
    for i = 1:epoch
        netStruct = load(fullfile('data',[ModelType '_' testDatabase '_TrainPropertion' num2str(trainingPropertion)...
                    '_PatchSize' num2str(patchSize) '_PatchNum' num2str(patchNum) '_Seed' num2str(seed) '_bins' num2str(bins)...
                    '_beta' num2str(beta)],['net-epoch-' num2str(i) '.mat']));
        net = dagnn.DagNN.loadobj(netStruct.net) ;
        move(net, 'gpu')
        net.mode = 'test';
        [SRCC(seed,i),PLCC(seed,i)] = testModel(ModelType, net, testDatabase,seed,trainingPropertion);
    end
end

%% save results
median_SRCC =[]; median_PLCC=[]; std_SRCC=[]; std_PLCC=[];
median_SRCC = median(SRCC);
median_PLCC = median(PLCC);
std_SRCC = std(SRCC);
std_PLCC = std(PLCC);
[max_SRCC,idx] = max(median_SRCC);
max_PLCC = median_PLCC(idx);
file = fopen(fullfile('result','results.txt'),'a');
fprintf(file,'Parameters: %s; %s; bins = %d; patches = %d;\n',testDatabase,ModelType,bins,patchNum);
fprintf(file,'SRCC:     ');fprintf(file,'%.4f; ',median_SRCC);fprintf(file,'\n');
fprintf(file,'std_SRCC: ');fprintf(file,'%.4f; ',std_SRCC);fprintf(file,'\n');
fprintf(file,'PLCC:     ');fprintf(file,'%.4f; ',median_PLCC);fprintf(file,'\n');
fprintf(file,'std_PLCC: ');fprintf(file,'%.4f; ',std_PLCC);fprintf(file,'\n');
fprintf(file,'best_SRCC: %.4f; best_PLCC: %.4f \n',max_SRCC,max_PLCC);
fclose(file);

end

author: HuiZeng


text:
% PRTBRVVB - PRT BRV Variational Bayes Object
%   Contains properties and methods for prtBrv objects that impliment
%   iterative variational Bayesian algorithms.
%
% Properties:
%   vbConvergenceThreshold = 1e-6; % This is a percentage of the total change in NFE
%   vbMaxIterations = 1e3; % Maximum number of iterations
%   vbVerboseText = false; % display text
%   vbVerbosePlot = false; % plot after each iteration
%   vbVerboseMovie = false; % make a movie with each frame of plotting
%   vbVerboseMovieFrames = []; % Where we store the frames
% Methods:
%   vbCheckConvergence - Used to check convergence at the end of each VB
%       iteration. Checks the percentage change in negative free energy.

classdef prtBrvVb

    methods (Abstract)
        [self, training] = vbBatch(self, x);
    end

    properties
        vbCheckConvergence = true; % Actually bother to calculate NFE
        vbConvergenceThreshold = 1e-6; % This is total change in NFE
        vbConvergenceDecreaseThreshold = 1e-3; % This is total change in NFE that will cause an exit
        vbMaxIterations = 1e3; % Maximum number of iterations
        vbVerboseText = false; % display text
        vbVerbosePlot = false; % plot after each iteration
        vbVerboseMovie = false; % make a movie with each frame of plotting
        vbVerboseMovieFrames = []; % Where we store the frames
    end

    methods

        function [self, training] = vb(self, x)
            [self, training] = vbBatch(self, x);
        end

        function [converged, err] = vbIsConverged(self, priorSelf, x, training) %#ok<INUSL>
            F = training.negativeFreeEnergy;
            pF = training.previousNegativeFreeEnergy;

            signedPercentage = (F-pF);
            converged = signedPercentage < self.vbConvergenceThreshold;

            if self.vbVerboseText
                fprintf('\t\t\t Negative Free Energy: %0.2f, %+.2e Change',F,signedPercentage)
            end

            if signedPercentage < 0
                if self.vbVerboseText
                    fprintf('\tDecrease in negative free energy occured!!!')
                end
                converged = false;
                if abs(signedPercentage)>self.vbConvergenceDecreaseThreshold
                    if self.vbVerboseText
                        fprintf('\n\t\tDecrease in negative free energy is beyond numerical error expectations. Exiting.')
                    end
                    err = true;
                else
                    err = false;
                end
            else
                err = false;
            end

            if self.vbVerboseText
                fprintf('\n');
            end

        end
        function [converged, err] = vbIsConvergedAbs(self, priorSelf, x, training) %#ok<INUSL>
            F = training.negativeFreeEnergy;
            pF = training.previousNegativeFreeEnergy;

            signedPercentage = (F-pF);
            converged = abs(signedPercentage) < self.vbConvergenceThreshold;

            if self.vbVerboseText
                fprintf('\t\t\t Negative Free Energy: %0.2f, %+.2e Change\n',F,signedPercentage)
            end

            if signedPercentage < 0
                %if self.vbVerboseText
                %    fprintf('\t\tDecrease in negative free energy occured!!! Exiting.\n')
                %end
                %converged = false;
                err = true;
            else
                err = false;
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function saveresults(res,peval,p, name_res_file)
% saveresults(res,peval,p, name_res_file)

if ~exist(peval.res_dir,'dir')
    fprintf('Creating directory:\n%s\n', peval.res_dir)
    mkdir (peval.res_dir);
end

if exist('peval.dir_res','var');
    savehere = [peval.dir_res '/' name_res_file];
else
    savehere = name_res_file;
end

if ~isfield (peval,'fid') peval.fid=1;end
mfprintf(peval.fid, ['Saving data to:\n%s\n'], savehere)
save (savehere,'res', 'peval', 'p');

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Script to illustrate the computation of the hydration layer fraction.
%  s_HLF
%
% For reference see:
%

%% Load and align the data
%
dataDir = 'C:\u\brian\Matlab\mrDataExample\hydrationLayer';
refIm = 't1.nii.gz';
rawDir = 'raw';

% 1 for trilinear, 7 for b-spline
interp = 7;
if(interp==1), outDirName = 'trilin';
else           outDirName = 'bsplin'; end

% Fit method:
%   'f' for fminsearch
%   'l' for lsqnolin
%   'pf','pl' for parallel
% versions with disttoolbox
% What is disttoolbox?
fitMethod = 'ls';
outMontage = ['./wst1.png'];

outBaseName = outDirName;
outFile = fullfile(dataDir,[outBaseName '_aligned.mat']);
if(exist(outFile,'file'))
    disp(['Loading aligned data from ' outFile '...']);
    load(outFile);
else
    outDir = dataDir;
    % Load all the series in the struct 's'
    s = dicomLoadAllSeries(rawDir);
    ref = niftiRead(refIm);
    %mmPerVox = ref.pixdim(1:3);
    mmPerVox = [0.9 0.9 0.9];
    % Align all the series to this subject's reference volume
    [s,xform,alignInds] = relaxAlignAll(s,ref,mmPerVox,true,interp);
    s = s(alignInds);
    save(outFile,'s', 'xform', 'dataDir', 'refIm', 'outDir', 'mmPerVox');
end

% Drop the first measure- too much motion
s = s(2:end);

outDir = fullfile(outDir, [outDirName fitMethod(end)]);
if(~exist(outDir,'dir')), mkdir(outDir); end

%
% Sort out image types
%
sequenceNames = {s(:).sequenceName};
spgrInds = false(size(sequenceNames));

% find series that are *either* SPGR or FSPGR
spgrInds(strmatch('3DGRASS',sequenceNames)) = 1;
spgrInds(strmatch('EFGRE3D',sequenceNames)) = 1;

tiInds = spgrInds & [s(:).inversionTime]>0;
t1Inds = spgrInds & ~tiInds;

% Mean image, use only t1Inds to calculate it
mn = mean(cat(4,s([2,4,5]).imData),4);
[brainMask,checkSlices] = mrAnatExtractBrain(mn, mmPerVox, 0.40);
%brainMask = dtiCleanImageMask(brainMask,3,1,.2);
%brainMask = dtiCleanImageMask(mrAnatHistogramClip(mn,0.3,0.9)>0.4,10,1,.3);
brainMask(any(cat(4,s(t1Inds).imData)<=0,4)) = 0;

tr = [s(t1Inds).TR];
if(~all(tr==tr(1))), error('TR''s do not match!'); end
tr = tr(1);

flipAngles = [s(t1Inds).flipAngle];

% compute t1 with no b1 bias correction
[t1,gPD_T2s] = relaxFitT1(cat(4,s(t1Inds).imData),flipAngles,tr);
t1(~brainMask) = 0; gPD_T2s(~brainMask) = 0;
dtiWriteNiftiWrapper(single(t1), xform, fullfile(outDir,'T1_nob1.nii.gz'));
dtiWriteNiftiWrapper(single(pd), xform, fullfile(outDir,'gPD_T2s_nob1.nii.gz'));

%% Compute a b1 map using the DESPOT1-HIFI method
%b1Map = 1.0;
if(~exist(fullfile(outDir,'b1.nii.gz')))
    b1Map = relaxFitDespoT1(s, brainMask);
    dtiWriteNiftiWrapper(single(b1Map), xform, fullfile(outDir,'b1.nii.gz'));
else
    ni = niftiRead(fullfile(outDir,'b1.nii.gz'));
    b1Map = double(ni.data);
end

[t1,gPD_T2s] = relaxFitT1(cat(4,s(t1Inds).imData),flipAngles,tr,b1Map);
t1(~brainMask) = 0; gPD_T2s(~brainMask) = 0;
t1(t1>5) = 5;

dtiWriteNiftiWrapper(single(t1), xform, fullfile(outDir,'T1.nii.gz'));
dtiWriteNiftiWrapper(single(gPD_T2s), xform, fullfile(outDir,'gPD_T2s.nii.gz'));

% Compute R1
t1(t1>3) = 3;
t1(brainMask & t1<0.5) = 0.5;
r1 = zeros(size(t1));
r1(brainMask) = 1./t1(brainMask);
dtiWriteNiftiWrapper(single(r1), xform, fullfile(outDir,'R1.nii.gz'));

% FSL commands
% fast R1.nii.gz
% run_first_all -b -i R1.nii.gz -o first

%% GridFit to the gPDt2* image the T2* in neglected hopefully with good reason
[PDc CoilG]=BaisGridFit(t1,PD,xform,outDir);

%% Final HLF computation
mField = 3;
[HLF,Wf] = hlfCompute(t1,pd,xform,outDir,mField)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% NEWTIMEFTRIALBASELN - Remove baseline power values for single trials in
%        newtimef. This only remove single trial baseline (not the average
%        baseline).
%
% Usage:
%   >>  tf = newtimefbaseln(tf, tvals, 'key', val);
%
% Inputs:
%   tf       - [3-D or 4-D array] single-trial spectral estimates
%           [freqs x times x trials] or [channels x freqs x times x trials)
%           The function may also process cell arrays of such inputs.
%   tvals    - [array] time values
%
% Optional inputs: 'baseline', 'basenorm' and 'trialbase'. Same definition
%  as for newtimef. If trialbase is 'off' this function does nothing.
%
% Outputs:
%   tf        - Baseline correct power (same size as input)
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, August 2016

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, 2016, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function PP = newtimeftrialbaseln(PPori, timesout, varargin)

if nargin < 3
    help newtimefbaseln;
    return;
end

[g timefreqopts ] = finputcheck(varargin, ...
    {'basenorm'      'string'    {'on','off'} 'off';
    'baseline'      'real'      []          0;
    'trialbase'     'string'    {'on','off','full'} 'off';
    'verbose'       'string'    {'on','off'} 'on';
    }, 'newtimeftrialbaseln', 'ignore');
if ischar(g) error(g); return; end
PP = PPori; if ~iscell(PP), PP = { PP }; end

% ---------------
% baseline length
% ---------------
if size(g.baseline,2) == 2
    baseln = [];
    for index = 1:size(g.baseline,1)
        tmptime   = find(timesout >= g.baseline(index,1) & timesout <= g.baseline(index,2));
        baseln = union_bc(baseln, tmptime);
    end
    if length(baseln)==0
        error( [ 'There are no sample points found in the default baseline.' 10 ...
            'This may happen even though data time limits overlap with' 10 ...
            'the baseline period (because of the time-freq. window width).' 10 ...
            'Either disable the baseline, change the baseline limits.' ] );
    end
else
    if ~isempty(find(timesout < g.baseline))
         baseln = find(timesout < g.baseline); % subtract means of pre-0 (centered) windows
    else baseln = 1:length(timesout); % use all times as baseline
    end
end

for ind = 1:length(PP(:))

    P = PP{ind};

    % -----------------------------------------
    % remove baseline on a trial by trial basis
    % -----------------------------------------
    if strcmpi(g.trialbase, 'on'), tmpbase = baseln;
    else                           tmpbase = 1:size(P,2); % full baseline
    end
    if ~strcmpi(g.trialbase, 'off')
        if ndims(P) == 4
            mbase = mean(P(:,:,tmpbase,:),3);
            if strcmpi(g.basenorm, 'on')
                mstd = std(P(:,:,tmpbase,:),[],3);
                P = bsxfun(@rdivide, bsxfun(@minus, P, mbase), mstd);
            else P = bsxfun(@rdivide, P, mbase);
            end
        else
            mbase = mean(P(:,tmpbase,:),2);
            if strcmpi(g.basenorm, 'on')
                mstd = std(P(:,tmpbase,:),[],2);
                P = (P-repmat(mbase,[1 size(P,2) 1]))./repmat(mstd,[1 size(P,2) 1]); % convert to log then back to normal
            else
                P = P./repmat(mbase,[1 size(P,2) 1]);
                %P = 10 .^ (log10(P) - repmat(log10(mbase),[1 size(P,2) 1])); % same as above
            end
        end
    end

    PP{ind} = P;
end
if ~iscell(PPori) PP = PP{1}; end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [EC,ec,degij] = edge_nei_overlap_bu(CIJ)
% EDGE_NEI_OVERLAP_BU        Overlap amongst neighbors of two adjacent nodes
%
%   [EC,ec,degij] = edge_nei_bu(CIJ);
%
%   This function determines the neighbors of two nodes that are linked by
%   an edge, and then computes their overlap.  Connection matrix must be
%   binary and directed.  Entries of 'EC' that are 'inf' indicate that no
%   edge is present.  Entries of 'EC' that are 0 denote "local bridges", i.e.
%   edges that link completely non-overlapping neighborhoods.  Low values
%   of EC indicate edges that are "weak ties".
%
%   If CIJ is weighted, the weights are ignored.
%
%   Inputs:     CIJ,    undirected (binary/weighted) connection matrix
%
%   Outputs:    EC,     edge neighborhood overlap matrix
%               ec,     edge neighborhood overlap per edge, in vector format
%               degij,  degrees of node pairs connected by each edge
%
%   Reference: Easley and Kleinberg (2010) Networks, Crowds, and Markets.
%              Cambridge University Press, Chapter 3.
%
%   Olaf Sporns, Indiana University, 2012

[ik,jk,ck] = find(CIJ);
lel = length(ck);
N = size(CIJ,1);

[deg] = degrees_und(CIJ);

ec = zeros(1,lel);
degij = zeros(2,lel);
for e=1:lel
    neiik = setdiff(union(find(CIJ(ik(e),:)),find(CIJ(:,ik(e))')),[ik(e) jk(e)]);
    neijk = setdiff(union(find(CIJ(jk(e),:)),find(CIJ(:,jk(e))')),[ik(e) jk(e)]);
    ec(e) = length(intersect(neiik,neijk))/length(union(neiik,neijk));
    degij(:,e) = [deg(ik(e)) deg(jk(e))];
end;

ff = find(CIJ);
EC = 1./zeros(N);
EC(ff) = ec;                        %#ok<FNDSB>

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [mi2] = yd22mi2(yd2)
% Convert area from square yards to square miles.
% Chad A. Greene 2012
mi2 = yd2*3.228305785124E-7;

author: Sable

repo: mcbench-benchmarks


text:
function [d, pt1, pt2] = distanceLines3d(line1, line2)
%DISTANCELINES3D Minimal distance between two 3D lines.
%
%   D = distanceLines3d(LINE1, LINE2);
%   Returns the distance between line LINE1 and the line LINE2, given as:
%   LINE1 : [x0 y0 z0 dx dy dz] (or M-by-6 array)
%   LINE2 : [x0 y0 z0 dx dy dz] (or N-by-6 array)
%   D     : (positive) array M-by-N
%
%   [D, PT1, PT2] = distanceLines3d(LINE1, LINE2);
%   Also returns the points located on LINE1 and LINE2 corresponding to the
%   shortest distance.
%   One should get the following:
%   distancePoints3d(PT1, PT2) - D == 0
%
%
%   Example
%     line1 = [2 3 4 0 1 0];
%     line2 = [8 8 8 0 0 1];
%     distanceLines3d(line1, line2)
%     ans =
%         6.0000
%
%   See also:
%   lines3d, distancePoints3d
%
%   ---------
%   authors: Brandon Baker, oqilipo, David Legland
%   created January 19, 2011
%

% number of points of each array
n1 = size(line1, 1);
n2 = size(line2, 1);

if nargout <= 1
    % express line coordinate as n1-by-n2 arrays
    v1x = repmat(line1(:,4), [1 n2]);
    v1y = repmat(line1(:,5), [1 n2]);
    v1z = repmat(line1(:,6), [1 n2]);
    p1x = repmat(line1(:,1), [1 n2]);
    p1y = repmat(line1(:,2), [1 n2]);
    p1z = repmat(line1(:,3), [1 n2]);

    v2x = repmat(line2(:,4)', [n1 1]);
    v2y = repmat(line2(:,5)', [n1 1]);
    v2z = repmat(line2(:,6)', [n1 1]);
    p2x = repmat(line2(:,1)', [n1 1]);
    p2y = repmat(line2(:,2)', [n1 1]);
    p2z = repmat(line2(:,3)', [n1 1]);

    % calculates distance for each set of lines
    vcross = cross([v1x(:) v1y(:) v1z(:)], [v2x(:) v2y(:) v2z(:)]);
    num = ([p1x(:) p1y(:) p1z(:)] - [p2x(:) p2y(:) p2z(:)]) .* vcross;
    t1 = sum(num,2);
    d = abs(t1) ./ (vectorNorm3d(vcross) + eps);

    % returns result as n1-by-n2 array
    d = reshape(d, n1, n2);

else
    % check input dimension, as we need to be able to match each pair of
    % lines
    if n1 ~= n2
        error('geom3d:distanceLines3d:IllegalInputArgument', ...
            'when output points are requested, number of lines should be the same');
    end

    p1 = line1(:, 1:3);
    p2 = line2(:, 1:3);
    dp = p2 - p1;
    v1 = line1(:, 4:6);
    v2 = line2(:, 4:6);

    % compute distance
    vcross = cross(v1, v2, 2);
    num = dp .* vcross;
    t1 = sum(num, 2);
    d = abs(t1) ./ (vectorNorm3d(vcross) + eps);

    % precomputations
    a = dot(v1, v1, 2);
    b = dot(v1, v2, 2);
    e = dot(v2, v2, 2);
    den = a.*e - b.*b; % 0, if lines are parallel

    % vector between origin of both lines
    r = line1(:,1:3) - line2(:,1:3);

    % solve linear system
    c = dot(v1, r, 2);
    f = dot(v2, r, 2);
    s = (b .* f - c .* e) ./ den;
    t = (a .* f - c .* b) ./ den;

    % convert to coordinates of points on lines
    pt1 = line1(:,1:3) + v1 .* s;
    pt2 = line2(:,1:3) + v2 .* t;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% Background Subtraction demo
% An example using drawContours to clean up a background segmentation result.
%
% This program demonstrates a simple method of connected components clean up
% of background subtraction.
% When the program starts, it begins learning the background.
% You can toggle background learning on and off using the checkbox.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/segment_objects.cpp>
%

%% Set up video source: video file or camera
fname = fullfile(mexopencv.root(),'test','768x576.avi');
if exist(fname, 'file') ~= 2
    % download video from Github
    disp('Downloading video...')
    url = 'https://cdn.rawgit.com/opencv/opencv/3.1.0/samples/data/768x576.avi';
    urlwrite(url, fname);
end
cap = cv.VideoCapture(fname);
%cap = cv.VideoCapture(0); pause(2);
assert(cap.isOpened(), 'Cannot open video source');

frame = cap.read();
assert(~isempty(frame), 'Cannot read data from the video source');

%% Create a background subtractor
bs = cv.BackgroundSubtractorMOG2('VarThreshold',10);

%%
ZR = zeros(size(frame,1), size(frame,2), 'uint8');
niters = 3;

%% Prepare window
hFig = figure('Position',[100 100 840 630], ...
    'KeyPressFcn',@(o,e)setappdata(o,'flag',true));
setappdata(hFig, 'flag',false);
subplot(221), hImg(1) = imshow(frame); title('video')
subplot(222), hImg(2) = imshow(frame); title('BG model')
subplot(223), hImg(3) = imshow(ZR); title('FG segmented')
subplot(224), hImg(4) = imshow(ZR); title('FG largest component')
hCB = uicontrol('Style','checkbox', 'Position',[20 20 200 20], ...
    'String','Update Background Model', 'Value',true);

%% Main loop
while ishghandle(hFig)
    % get next frame
    frame = cap.read();
    if isempty(frame), break; end

    % determine whether to update BG model or not, then compute FG mask
    update_bg_model = get(hCB,'Value');
    if update_bg_model
        learnRate = -1; % automatically chosen learning rate
    else
        learnRate = 0;  % dont update BG model
    end
    fgmask = bs.apply(frame, 'LearningRate',learnRate);

    % process mask and extract largest connected component
    bw = cv.dilate(fgmask, 'Iterations',niters);
    bw = cv.erode(bw, 'Iterations',niters*2);
    bw = cv.dilate(bw, 'Iterations',niters);
    [contours, hierarchy] = cv.findContours(bw, ...
        'Mode','CComp', 'Method','Simple');
    if ~isempty(contours)
        maxArea = 0;
        largestComp = 0;
        idx = 0;
        while idx >= 0
            % compute contour area
            a = abs(cv.contourArea(contours{idx+1}));
            if a > maxArea
                maxArea = a;
                largestComp = idx;
            end
            % iterate through all the top-level contours
            idx = hierarchy{idx+1}(1);
        end
        out_frame = cv.drawContours(ZR, contours, ...
            'ContourIdx',largestComp, 'Hierarchy',hierarchy, ...
            'Color',255, 'Thickness','Filled');
    else
        out_frame = ZR;
    end

    % update images
    set(hImg(1), 'CData',frame)     % video frame
    if update_bg_model
        % get the current BG model
        set(hImg(2), 'CData',bs.getBackgroundImage())
    end
    set(hImg(3), 'CData',fgmask)    % FG mask
    set(hImg(4), 'CData',out_frame) % largest contour in FG mask

    % Terminate if any user input
    flag = getappdata(hFig, 'flag');
    if isempty(flag)||flag, break; end
    pause(0.1);
end

%%
% release camera
cap.release()

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function fx1 = p04_fx1 ( x )

%*****************************************************************************80
%
%% P04_FX1 evaluates the derivative of the function for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX1, the first derivative of the function at X.
%
  fx1 = exp ( x ) + 2.0 / ( 100.0 * x * x * x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef (Abstract) FilterX < BaseX % Extends trackingX.BaseX
% FilterX Abstract class
%
% Summary of FilterX:
% This is the base class for all TrackingX filters.
% Any custom defined Filter should be derived from this FilterX base class.
%
% By default, FilterX objects are designed with the intention to be utilised
% as (Markov-Chain) State Estimators that opererate on a given problem as
% deterministic Finite State Machines (FSMs), with self-contained memory.
% In other words, it is expected that for most problems, FilterX objects will
% go through any heavy parameterisation once, after which they will be used
% to recursively execute a Prediction-Update State Estimation loop.
% FilterX objects aim to preserve all the inter-state information they require
% such that they can execute the next FSM step, without being provided with
% information they already know.
%
%
% FilterX Properties:
%   + Model      - Object handle to StateSpaceModelX class, containing model
%                  parameterisations. Normally set once during filter initialisation and
%                  accessed later on.
%   + Prior      - Structure containing prior information (e.g. pdf)
%                  Normally set once during filter initialisation and
%                  accessed later on.
%   + Prediction - Structure containing prediction information (e.g. pdf)
%                  Data is over-written on every prediction step and is
%                  subsequently read by the update (and later) step(s)
%   + Posterior  - Structure containing posterior information (e.g. pdf)
%                  Data is over-written on every update step and is subsequently
%                  utilised by
%   + MeasurementList   - A (matrix of) column vector(s) representing the
%                         measurement(s) to be used during the filter update step
%
%   (*) Signifies properties necessary to instantiate a class object
%
% FilterX Methods:
%   + FilterX    - Constructor method
%   + predict    - Performs filter prediction step
%   + update     - Performs filter update/correction step
%
% (+) denotes puplic properties/methods
%
% See also TransitionModelX, MeasurementModelX and ControlModelX
%
% November 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
        Model
        MeasurementList
    end

    properties (Access = protected)

    end

    methods (Abstract)
        predict(this);
        update(this);
    end

    methods (Access = protected)
        function initialise_(this, config)
            if (isfield(config,'Model'))
                this.Model = config.Model;
            end
        end
    end

    methods
        function this = FilterX(varargin)
        % FilterX Constructor method
        %
        % Parameters
        % ----------
        % Model: StateSpaceModelX
        %  Object handle to a given state-space model
        %
        % Usage
        % -----
        % * FilterX(__, Name, Value) instantiates an object handle, configured
        %   the parameters specified by one or more Name,Value pair arguments.
        % * FilterX(config) instantiates an object handle configured with
        %   the parameters specified inside the 'config' structure, whose
        %   fieldnames correspond to the parameter names as given in the
        %   Parameters section above.
        %
        %  See also predict, update, iterate, smooth.

            if(nargin==0)
                return;
            end

            % First check to see if a structure was received
            if(nargin==1)
                if(isstruct(varargin{1}))
                    this.initialise_(varargin{1});
                    return;
                end
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            this.initialise_(parser.Unmatched);
        end

        function initialise(this, varargin)
        % initialise Initialisation method. Reset and re-configure the filter.
        %
        % Parameters
        % ----------
        % Model: StateSpaceModelX
        %  Object handle to a given state-space model.
        %
        % Usage
        % -----
        % * filter.initialise(__, Name, Value) instantiates an object handle,
        %   configured the parameters specified by one or more Name,Value
        %   pair arguments.
        % * filter.initialise(config) instantiates an object handle configured
        %   with the parameters specified inside the 'config' structure, whose
        %   fieldnames correspond to the parameter names as given in the
        %   Parameters section above.
        %
        %  See also predict, update, iterate, smooth.

            if(nargin==1)
                return;
            end

            % First check to see if a structure was received
            if(nargin==2)
                if(isstruct(varargin{1}))
                    this.initialise_(varargin{1});
                    return;
                end
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            this.initialise_(parser.Unmatched);
        end

        % ===============================>
        % ACCESS METHODS
        % ===============================>
        function set.Model(this,newModel)
            this.Model = setModel(this,newModel);
        end

        function set.MeasurementList(this,newMeasurementList)
            this.MeasurementList = setMeasurementList(this,newMeasurementList);
        end
    end

    methods (Access = protected)
        % ===============================>
        % ACCESS METHOD HANDLES
        % ===============================>
        function Prior = setPrior(this,newPrior)
            Prior = newPrior;
        end
        function Prediction = setPrediction(this,newPrediction)
            Prediction = newPrediction;
        end
        function Posterior = setPosterior(this,newPosterior)
            Posterior = newPosterior;
        end
        function Model = setModel(this,newModel)
            Model = newModel;
        end
        function measurementList = setMeasurementList(this,newMeasurementList)
            measurementList = newMeasurementList;
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function test_model = merge_model_arr(test_model_arr)
    test_model.class = test_model_arr(1).class;
    test_model.S_bar = test_model_arr(1).S_bar;
    test_model.defBasis = test_model_arr(1).defBasis;
    test_model.part_names = test_model_arr(1).part_names;
    test_model.sigma_sq = test_model_arr(1).sigma_sq;
    test_model.points3 = vertcat(test_model_arr(:).points3);
    test_model.points3 = [test_model.points3(1:3:end,:);test_model.points3(2:3:end,:);...
        test_model.points3(3:3:end,:)];
    test_model.rots = [test_model_arr(:).rots];
    test_model.trs = vertcat(test_model_arr(:).trs);
    test_model.def_weights = vertcat(test_model_arr(:).def_weights);
    test_model.voc_image_id = [test_model_arr(:).voc_image_id]';
    if(isfield(test_model_arr(1),'voc_rec_id'))
        test_model.voc_rec_id = vertcat(test_model_arr(:).voc_rec_id);
    end
    test_model.bbox = vertcat(test_model_arr(:).bbox);
    test_model.points2 = vertcat(test_model_arr(:).points2);
    test_model.points2 = [test_model.points2(1:2:end,:);test_model.points2(2:2:end,:)];
    test_model.params = test_model_arr(1).params;
    test_model.c = vertcat(test_model_arr(:).c);
    if(isfield(test_model_arr(1),'flip'))
        test_model.flip = vertcat(test_model_arr(:).flip);
    end

    try
        for i=1:length(test_model_arr)
            test_model.seg.poly_x{i} = test_model_arr(i).seg.poly_x{1};
            test_model.seg.poly_y{i} = test_model_arr(i).seg.poly_y{1};
        end
        test_model.seg.poly_x = test_model.seg.poly_x';
        test_model.seg.poly_y = test_model.seg.poly_y';
    catch
        test_model.seg = {test_model_arr(:).seg}';
    end
end

repo: CategoryShapes


text:
function [ V VS ] = dicomreadvolume(zipfile)
%DICOMREADVOLUME   Read a zip-file with DICOM slices as volume.
%   [V VS] = DICOMREADVOLUME(FNAME) loads the 3-D volume data V from
%   the volume dicom slices in the zip-file FNAME.
%   The voxel size is returned by VS.
%
%   See also DICOMWRITEVOLUME, DICOMWRITE, DICOMREAD
%
%   Author: medical-image-processing.com

% extract the zip file to the temporary directory
fnames = unzip(zipfile, tempdir);

% number of files
N = numel(fnames);

if (N<1)
    error('Empty zip file');
end

% read first slice for determining slice properties
S = squeeze(dicomread(fnames{1}));
I1 = dicominfo(fnames{1});

% voxel size information
VS = [I1.PixelSpacing(:) ; I1.SliceThickness];

% slice size and datatype
sz = size(S);
tp = class(S);

% pre-allocate data
VT = zeros([sz N], tp);
V = VT;
POS = zeros(N,2);

% load each slice and its properties
for i=1:N
    VT(:,:,i) = squeeze(dicomread(fnames{i}));
    info = dicominfo(fnames{i});
    if isfield(info, 'ImagePositionPatient')
        POS(i,:) = [info.ImagePositionPatient(3) i];
    else
        POS(i,:) = [i i];
    end
    delete(fnames{i});
end

% resort the slices according to the image position
POS = sortrows(POS,1);
for i=1:N
   V(:,:,i) = VT(:,:,POS(i,2));
end


text:
function [rx,sx,tx,ry,sy,ty,rz,sz,tz,J] = GeometricFactors3D(x,y,z,Dr,Ds,Dt)

% function [rx,sx,tx,ry,sy,ty,rz,sz,tz,J] = GeometricFactors3D(x,y,z,Dr,Ds,Dt)
% Purpose  : Compute the metric elements for the local mappings of the elements

% calculate geometric factors
xr = Dr*x; xs = Ds*x; xt = Dt*x;
yr = Dr*y; ys = Ds*y; yt = Dt*y;
zr = Dr*z; zs = Ds*z; zt = Dt*z;

J = xr.*(ys.*zt-zs.*yt) - yr.*(xs.*zt-zs.*xt) + zr.*(xs.*yt-ys.*xt);
rx =  (ys.*zt - zs.*yt)./J; ry = -(xs.*zt - zs.*xt)./J; rz = (xs.*yt - ys.*xt)./J;
sx = -(yr.*zt - zr.*yt)./J; sy =  (xr.*zt - zr.*xt)./J; sz = -(xr.*yt - yr.*xt)./J;
tx =  (yr.*zs - zr.*ys)./J; ty = -(xr.*zs - zr.*xs)./J; tz = (xr.*ys - yr.*xs)./J;
return;

repo: Aero-matlab


text:
function VARopt = VARoption
%========================================================================
% Optional inputs for VAR analysis. This function is run automatically in
% the VARmodel function.
%========================================================================
% VARopt = VARoption
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

VARopt.vnames    = [];      % endogenous variables names
VARopt.vnames_ex = [];      % exogenous variables names
VARopt.snames    = [];      % shocks names
VARopt.nsteps    = 40;      % number of steps for computation of IRFs and FEVDs
VARopt.impact    = 0;       % size of the shock for IRFs: 0=1stdev, 1=unit shock
VARopt.shut      = 0;       % forces the IRF of one variable to zero
VARopt.ident     = 'short'; % identification method for IRFs ('short' zero short-run restr, 'long' zero long-run restr, 'sign' sign restr, 'iv' external instrument)
VARopt.recurs    = 'wold';  % method for computation of recursive stuff ('wold' form MA representation, 'comp' for companion form)
VARopt.ndraws    = 1000;    % number of draws for bootstrap or sign restrictions
VARopt.mult      = 10;      % multiple of draws to be printed at screen.
VARopt.pctg      = 95;      % confidence level for bootstrap
VARopt.method    = 'bs';    % methodology for error bands, 'bs' for standard bootstrap, 'wild' wild bootstrap
VARopt.sr_hor    = 1;       % number of periods that sign restrcitions are imposed on
VARopt.sr_rot    = 500;     % max number of rotations for finding sign restrictions
VARopt.sr_draw   = 100000;  % max number of total draws for finding sign restrictions
VARopt.sr_mod    = 1;       % model uncertainty for sign restrictions (1=yes, 0=no)
VARopt.pick      = 0;       % selects one variable for IRFs and FEVDs plots (0 => plot all)
VARopt.quality   = 1;       % quality of exported figures: 1=high (ghostscript required), 0=low
VARopt.suptitle  = 0;       % title on top of figures
VARopt.datesnum  = [];      % numeric vector of dates in the VAR
VARopt.datestxt  = [];      % cell vector of dates in the VAR
VARopt.datestype = 1;       % 1 smart labels; 2 less smart labels
VARopt.firstdate = [];      % initial date of the sample in format 1999.75 => 1999Q4 (both for annual and quarterly data)
VARopt.frequency = 'q';     % frequency of the data: 'm' monthly, 'q' quarterly, 'y' yearly
VARopt.figname   = [];      % string for name of exported figure
VARopt.FigSize   = [26,24]; % size of window for plots

% In progress
%VARopt.maxvd_N   = 1;       % position of variable to maximize variance decomposition
%VARopt.maxvd_H   = 10;      % horizon over which to maximize variance decomposition
%VARopt.maxvd_rot = 1000;    % max number of rotations for maximization of variance decomposition

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function h=pheaviside(L)
%PHEAVISIDE  Periodic Heaviside function
%   Usage: h=pheaviside(L);
%
%   `pheaviside(L)` returns a periodic Heaviside function. The periodic
%   Heaviside function takes on the value 1 for indices corresponding to
%   positive frequencies, 0 corresponding to negative frequencies and the
%   value .5 for the zero and Nyquist frequencies.
%
%   To get a function that weights the negative frequencies by 1 and the
%   positive by 0, use `involute(pheaviside(L))`
%
%   As an example, the `pheaviside` function can be use to calculate the
%   Hilbert transform for a column vector *f*::
%
%     h=2*ifft(fft(f).*pheaviside(length(f)));
%
%   See also: middlepad, involute, fftindex

%   AUTHOR : Peter L. Søndergaard.
%   REFERENCE: OK
%   TESTING: OK

complainif_argnonotinrange(nargin,1,1,mfilename);

h=zeros(L,1);
if L>0
    % First term is .5
    h(1)=.5;

    % Set positive frequencies to 1.
    h(2:ceil(L/2))=1;

    % Last term (Nyquist frequency) is also .5, if it exists.
    if rem(L,2)==0
        h(L/2+1)=.5;
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% EX_STOKES_DRIVEN_CAVITY_TH: solve the Stokes problem in the driven cavity with generalized Taylor-Hood elements.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined from the aspect ratio using the NURBS toolbox
aspect_ratio=[1.5 0.3];
nrb_section = nrb4surf([0 0], [1 0], [0 aspect_ratio(1)], [1 aspect_ratio(1)]);
problem_data.geo_name = nrbextrude (nrb_section, [0 0 aspect_ratio(2)]);

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:6;
problem_data.nmnn_sides = [];

% Physical parameters
problem_data.viscosity = @(x, y, z) ones (size (x));

% Force and boundary terms
problem_data.f  = @(x, y, z) zeros ([3, size(x)]);
problem_data.h  = @test_stokes_3d_symdrivcav_h_drchlt;

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.element_name = 'sg';        % Element type for discretization
method_data.degree       = [2  2  2];  % Degree of the splines (pressure space)
method_data.regularity   = [1  1  1];  % Regularity of the splines (pressure space)
method_data.nsub         = [2  2  2];  % Number of subdivisions
method_data.nquad        = [4  4  4];  % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = ...
                       solve_stokes (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file = 'Driven_cavity_3d_SG_Deg2_Reg1_Sub2';

fprintf ('The result is saved in the files %s \n and %s \n \n', ...
           [output_file '_vel'], [output_file '_press']);
vtk_pts = {linspace(0, 1, 15), linspace(0, 1, 15), linspace(0, 1, 15)};
sp_to_vtk (press, space_p, geometry, vtk_pts, [output_file '_press'], 'press')
sp_to_vtk (vel,   space_v, geometry, vtk_pts, [output_file '_vel'  ], {'velocity', 'divergence'}, {'value', 'divergence'})

%!test
%! aspect_ratio=[1.5 0.3];
%! nrb_section = nrb4surf([0 0], [1 0], [0 aspect_ratio(1)], [1 aspect_ratio(1)]);
%! problem_data.geo_name = nrbextrude (nrb_section, [0 0 aspect_ratio(2)]);
%! problem_data.drchlt_sides = 1:6;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y, z) ones (size (x));
%! problem_data.f  = @(x, y, z) zeros ([3, size(x)]);
%! problem_data.h  = @test_stokes_3d_symdrivcav_h_drchlt;
%! method_data.element_name = 'sg';        % Element type for discretization
%! method_data.degree       = [2  2  2];  % Degree of the splines (pressure space)
%! method_data.regularity   = [1  1  1];  % Regularity of the splines (pressure space)
%! method_data.nsub         = [2  2  2];  % Number of subdivisions
%! method_data.nquad        = [4  4  4];  % Points for the Gaussian quadrature rule
%! [geometry, msh, space_v, vel, space_p, press] = ...
%!                        solve_stokes (problem_data, method_data);
%! assert (msh.nel, 64)
%! assert (space_p.ndof, 64)
%! assert (space_v.ndof, 1029)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
try
    EEGDATA.marker_buffer(:,1+mod(EEGDATA.mmax:EEGDATA.mmax+size(EEGDATA_marker_chunk,2)-1,EEGDATA.marker_buffer_len)) = EEGDATA_marker_chunk;
catch e
    if any(strcmp(e.identifier,{'MATLAB:heterogeneousStrucAssignment','MATLAB:heterogenousStrucAssignment'}))
        EEGDATA.marker_buffer = EEGDATA_marker_chunk(ones(1,EEGDATA.marker_buffer_len));
        EEGDATA.marker_buffer(:,1+mod(EEGDATA.mmax:EEGDATA.mmax+size(EEGDATA_marker_chunk,2)-1,EEGDATA.marker_buffer_len)) = EEGDATA_marker_chunk;
    else
        rethrow(e);
    end
end

repo: REST


text:
if 0
    rxnForms = {' -> A','A -> B','B -> C', 'B -> D','D -> C','C ->'};
    rxnNames = {'Ain','AB','BC','BD','DC', 'Cout'};
    model = createModel(rxnNames, rxnNames,rxnForms);
    model.lb(3) = 2; %BD
    model.lb(4) = 2;
    model.ub(6) = 2;
else
    rxnForms = {' -> A','A -> B','A -> D','B -> C', 'B -> D','D -> C','C ->'};
    rxnNames = {'Ain','AB','AD','BC','BD','DC', 'Cout'};
    model = createModel(rxnNames, rxnNames,rxnForms);
    model.lb(1) = 5;
    %model.ub(6) = 2;
    model.ub(6) = 2.1;
    model.ub(4) = 1;
end

disp('-')
[nMet,nRxn]=size(model.S);

if exist('param','var')
clear param
end

if 0
    %weighting on zero norm of fluxes
    param.gamma0   = 0*1e-6;
    param.gamma1   = 1e-6;

    % weighting on relaxation of reaction bounds
    param.alpha0   = 10;
    param.alpha1   = 1;

    %weighting on relaxation of steady state constraints S*v = b
    param.lambda0   = 0;
    param.lambda1   = 0;
end

param.printLevel=2;

%stopping criterion
param.epsilon = 1e-6;
%capped l0 parameter
param.theta0   = 0.5;

sol = relaxedFBA(model,param);

model = findSExRxnInd(model);

return
plotRelaxedFBA(sol, model)

author: opencobra

repo: cobratoolbox


text:
%  Figure 10.24      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%   fig10_24.m is a script to generate Fig. 10.24,
%   the transient response of the LQR symmetric rootlocus compensator
%   of the satellite position control, non-colocated case WITH ESTIMATOR

% parameter values
m=[1, 0.1]; k=[0, 0.091] ; d=[0, 0.0036]; k1=[0, 0.4];
% call function
[f,g,h,j] = twomass(m,k,d);
s=[f, g;h, 0];r=[0*g;1]; n=s\r;nx=n(1:4);nu=n(5);
% call function
[f1,g,h,j] = twomass(m,k1,d);

% form G(s)G(-s) model
a=[f, 0*f;
-h'*h, -f'];
b=[g;0*g];
d=[0];
c=[0*h, g'];
hold off; clf
P=eig(a-b*c*0.1621);
pc=P(real(P<0)==1);
K=place(f,g,pc);
nbar=nu+K*nx;
% eig(f-g*K)
P=eig(a-b*c*3.056e7);
pe=P(real(P<0)==1);
L=place(f',h',pe)';
ac=f-g*K-L*h ;bc=L;cc=K;dc=0;
[Aol,Bol,Col,Dol]=series(ac,bc,cc,dc,f,g,h,j);
[acl,bcl,ccl,dcl]=feedback(f,g,h,j,ac,bc,cc,dc);
[acl1,bcl,ccl,dcl]=feedback(f1,g,h,j,ac,bc,cc,dc);
bcl= nbar*[g;g];
t=0:.25:30;
syscl=ss(acl,bcl,ccl,dcl);
step(syscl,t);
hold on;
grid;
gtext('nominal case')
syscl1=ss(acl1,bcl,ccl,dcl);
step(syscl1,t) ;
gtext('stiff-spring case');
title('Fig. 10.24 Closed-loop step response for the SRL design with an estimator')


text:
%
% Usage:   x =mexConjGrad(A,b,x0,tol,itermax)
%
% Name: mexConjGrad
%
% Description: Conjugate gradient algorithm, sometimes faster than the
%    equivalent Matlab function pcg. In order to solve Ax=b;
%
% Inputs: A:  double square n x n matrix. HAS TO BE POSITIVE DEFINITE
%         b:  double vector of length n.
%         x0: double vector of length n. (optional) initial guess.
%         tol: (optional) tolerance.
%         itermax: (optional) maximum number of iterations.
%
% Output: x: double vector of length n.
%
% Author: Julien Mairal, 2009

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%  Copyright (c) 2014, Karen Simonyan
%  All rights reserved.
%  This code is made available under the terms of the BSD license (see COPYING file).

function [res, extra] = eval(config, scores, gt)

    % predicted labels
    class = 2 * (scores >= config.threshold) - 1;

    % class-n accuracy
    res = mean(class == gt) * 100;

    extra = [];
end

author: AlfredXiangWu


text:
% Demonstrates how to use the property pane.
%
% See also: PropertyGrid

% Copyright 2010 Levente Hunyadi
function example_propertygrid

properties = [ ...
    PropertyGridField('double', pi, ...
        'Category', 'Primitive types', ...
        'DisplayName', 'real double', ...
        'Description', 'Standard MatLab type.') ...
    PropertyGridField('single', pi, ...
        'Category', 'Primitive types', ...
        'DisplayName', 'real single', ...
        'Description', 'Single-precision floating point number.') ...
    PropertyGridField('integer', int32(23), ...
        'Category', 'Primitive types', ...
        'DisplayName', 'int32', ...
        'Description', 'A 32-bit integer value.') ...
    PropertyGridField('interval', int32(2), ...
        'Type', PropertyType('int32', 'scalar', [0 6]), ...
        'Category', 'Primitive types', ...
        'DisplayName', 'int32', ...
        'Description', 'A 32-bit integer value with an interval domain.') ...
    PropertyGridField('enumerated', int32(-1), ...
        'Type', PropertyType('int32', 'scalar', {int32(-1), int32(0), int32(1)}), ...
        'Category', 'Primitive types', ...
        'DisplayName', 'int32', ...
        'Description', 'A 32-bit integer value with an enumerated domain.') ...
    PropertyGridField('logical', true, ...
        'Category', 'Primitive types', ...
        'DisplayName', 'logical', ...
        'Description', 'A Boolean value that takes either true or false.') ...
    PropertyGridField('doublematrix', [], ...
        'Type', PropertyType('denserealdouble', 'matrix'), ...
        'Category', 'Compound types', ...
        'DisplayName', 'real double matrix', ...
        'Description', 'Matrix of standard MatLab type with empty initial value.') ...
    PropertyGridField('string', 'a sample string', ...
        'Category', 'Compound types', ...
        'DisplayName', 'string', ...
        'Description', 'A row vector of characters.') ...
    PropertyGridField('rowcellstr', {'a sample string','spanning multiple','lines'}, ...
        'Category', 'Compound types', ...
        'DisplayName', 'cell row of strings', ...
        'Description', 'A row cell array whose every element is a string (char array).') ...
    PropertyGridField('colcellstr', {'a sample string';'spanning multiple';'lines'}, ...
        'Category', 'Compound types', ...
        'DisplayName', 'cell column of strings', ...
        'Description', 'A column cell array whose every element is a string (char array).') ...
    PropertyGridField('season', 'spring', ...
        'Type', PropertyType('char', 'row', {'spring','summer','fall','winter'}), ...
        'Category', 'Compound types', ...
        'DisplayName', 'string', ...
        'Description', 'A row vector of characters that can take any of the predefined set of values.') ...
    PropertyGridField('set', [true false true], ...
        'Type', PropertyType('logical', 'row', {'A','B','C'}), ...
        'Category', 'Compound types', ...
        'DisplayName', 'set', ...
        'Description', 'A logical vector that serves an indicator of which elements from a universe are included in the set.') ...
    PropertyGridField('root', [], ...  % [] (and no type explicitly set) indicates that value is not editable
        'Category', 'Compound types', ...
        'DisplayName', 'root node') ...
    PropertyGridField('root.parent', int32(23), ...
        'Category', 'Compound types', ...
        'DisplayName', 'parent node') ...
    PropertyGridField('root.parent.child', int32(2007), ...
        'Category', 'Compound types', ...
        'DisplayName', 'child node') ...
];

% arrange flat list into a hierarchy based on qualified names
properties = properties.GetHierarchy();

% create figure
f = figure( ...
    'MenuBar', 'none', ...
    'Name', 'Property grid demo - Copyright 2010 Levente Hunyadi', ...
    'NumberTitle', 'off', ...
    'Toolbar', 'none');

% procedural usage
g = PropertyGrid(f, ...            % add property pane to figure
    'Properties', properties, ...  % set properties explicitly
    'Position', [0 0 0.5 1]);
h = PropertyGrid(f, ...
    'Position', [0.5 0 0.5 1]);

% declarative usage, bind object to grid
obj = SampleObject;  % a value object
h.Item = obj;        % bind object, discards any previously set properties

% update the type of a property assigned with type autodiscovery
userproperties = PropertyGridField.GenerateFrom(obj);
userproperties.FindByName('IntegerMatrix').Type = PropertyType('denserealdouble', 'matrix');
disp(userproperties.FindByName('IntegerMatrix').Type);
h.Bind(obj, userproperties);

% wait for figure to close
uiwait(f);

% display all properties and their values on screen
disp('Left-hand property grid');
disp(g.GetPropertyValues());
disp('Right-hand property grid');
disp(h.GetPropertyValues());
disp('SampleObject (modified)');
disp(h.Item);
disp('SampleNestedObject (modified)');
disp(h.Item.NestedObject);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function url = prtGraphDataUrl
% url = prtGraphDataUrl
%   Return the url where the big graph data can be found.

url = 'http://www.duke.edu/~pat7/prtGraphData/';

author: covartech

repo: PRT


text:
function dt6_new = dtiResliceTensorSn3d(dt6, sn3dParams, nonLinearFlag, mmPerVoxOut, xformToTensor, rotTensorMethod)
% dt6_new = dtiResliceTensor(dt6, sn3dParams, [nonLinearFlag], [mmPerVoxOut], [xformToTensor], [rotTensorMethod])
%
% Applies spm-style sn3d parameters to a tensor volume in dt6 format.
% These sn3d parameters include a linear (affine) component and a non-linear component.
%
% INPUT:
%   dt6                 XxYxZx6 tensor array. If last dimension is 1 then data is
%                       treated as scalar volume.
%   sn3dParams          The sn3d parameter structure from SPM.
%   nonLinearFlag       If 1 (default), non-linear component is included.
%                       If 0, then non-linear component is not included.
%   mmPerVoxOut         mm per voxel in resliced array.
%   xformToTensor       Linear xform required with raw data. Not required if
%                       data is in Analyze format (set to identity).
%   rotTensorMethod     Method of tensor rotation. Ignored if data is scalar.
%                           'SPM' (using spm_matrix.m) - default
%                           'FS' (finite strain)
%                           'PPD' (preservation fo principal direction)
%                           '' (no rotation)
%
% OUTPUT:
%   dt6_new             XxYxZx6 tensor array, or XxYxZ if data is scalar.
%
% REQUIRES:
% spm99 or spm2 in the path.
%
% HISTORY:
% 2004.09.13 ASH & RFD Wrote it.
% 2005.01.12 RFD: Cleaning code, adding support for SPM2 spatial norm
% params, and adding code to do a proper tensor correction with nonlinear
% deformations.
%
% Example:
% sn3dParams = load('dti_analyses/may021126_B0_sn3d.mat');
% xformToB0 = computeCannonicalXformFromIfile('dti/B0.001');
% xformToTensor = inv(xformToB0);
% [eigVec,eigVal,mm] = dtiLoadTensor('/snarp/u1/dti/adultData/may021126/dti/Vectors.float');
% dt6 = dtiRebuildTensor(eigVec, eigVal);
% dt6_new = dtiResliceTensorSn3d(dt6, sn3dParams, mmPerVoxOut, xformToTensor);
% [newVec, newVal] = dtiSplitTensor(dt6_new);
% fa = dtiComputeFA(newVal);
% figure; imagesc(fa(:,:,23)); axis equal; colormap gray; axis xy;

if (size(dt6,4)~=6 & size(dt6,4)~=1),
    error('Wrong input format')
end
if(~exist('nonLinearFlag') | isempty(nonLinearFlag))
    nonLinearFlag = 1;
end
if(~exist('mmPerVoxOut') | isempty(mmPerVoxOut))
    mmPerVoxOut = [2 2 2];
end
if(~exist('xformToTensor') | isempty(xformToTensor))
    xformToTensor = eye(4);
end
if(~exist('rotTensorMethod') | isempty(rotTensorMethod))
    rotTensorMethod = 'SPM';
end

% Create a bounding box for the resliced data.
% The following is the spm default bounding box. Note that the box is
% defined in Talairach space (units = mm).
bb = [-78 -112 -50;
       78  76   85];

x = (bb(1,1):mmPerVoxOut(1):bb(2,1));
y = (bb(1,2):mmPerVoxOut(2):bb(2,2));
z = (bb(1,3):mmPerVoxOut(3):bb(2,3));

if(isfield(sn3dParams, 'Dims'))
    mm = sn3dParams.Dims(6,:)';
else
    % support spm2 style sn params
    mm = diag(sn3dParams.VF.mat(1:3,1:3));
end

[X,Y,Z] = meshgrid(x, y, z);
clear x y z;

talCoords = [X(:) Y(:) Z(:)];
newSize = size(X);
clear X Y Z;

imgCoords = mrAnatGetImageCoordsFromSn(sn3dParams, talCoords, nonLinearFlag);
if(isfield(sn3dParams,'VF'))
    [trans,rot,mm,skew] = affineDecompose(sn3dParams.VF.mat);
    imgCoords = mrAnatXformCoords(inv(sn3dParams.VF.mat)*xformToTensor, imgCoords')';
else
    mm = sn3dParams.Dims(6,:)';
    imgCoords = imgCoords./repmat(mm,1,size(imgCoords,2));
    imgOrigin = sn3dParams.Dims(5,:)'/2;
    imgCoords = imgCoords+repmat(imgOrigin,1,size(imgCoords,2));
    % We now have image coords for the space defined by the bounding box (in
    % Talairach space). These image coords are in the space of the original
    % image that was used to compute the sn3d params. But that may not be the
    % same space as the tensor image.
    imgCoords = mrAnatXformCoords(xformToTensor, imgCoords);
    imgCoords(1,:) = imgCoords(1,:).*mm(1);
    imgCoords(2,:) = imgCoords(2,:).*mm(2);
    imgCoords(3,:) = imgCoords(3,:).*mm(3);
end

fprintf('Interpolate coordinate grid...\n')
if (size(dt6,4)==6),
    % convert from matlab 1-indexing to C 0-indexing
	imgCoords = imgCoords - 1;

	% *** HACK so that matlab can find the shared library
	old = pwd;
	cd(fileparts(which('dtiTensorInterp_Pajevic')));
	dt6_new = dtiTensorInterp_Pajevic(dt6, [imgCoords(2,:);imgCoords(1,:);imgCoords(3,:)]', mm, 1, mm./2);
	cd(old);
    dt6_new = reshape(dt6_new, [newSize, 6]);
else
    dt6_new = myCinterp3(dt6, [size(dt6,1) size(dt6,2)], size(dt6,3), imgCoords', 0.0);
    dt6_new = reshape(dt6_new, [newSize, 1]);
end
dt6_new = real(dt6_new);

% ROTATE THE TENSORS
fprintf('Rotate tensors...\n')
xformTalToImg = sn3dParams.MF * sn3dParams.Affine * inv(sn3dParams.MG);
xformImgToTal = inv(xformTalToImg);

% Apply just the rigid rotation component of the xform to the tensors
if (size(dt6,4)==6),
    switch rotTensorMethod,
    case 'SPM',
        % spm_matrix method
		p = spm_imatrix(xformImgToTal);
		p([1:3,10:12]) = 0; p(7:9) = 1;
		rot = spm_matrix(p); rot = rot(1:3,1:3);
     	dt6_new = dtiXformTensors(dt6_new, rot);
    case 'FS',
        % Finite strain method
    	rigidXform = dtiFiniteStrainDecompose(xformImgToTal);
     	dt6_new = dtiXformTensors(dt6_new, rigidXform);
    case 'PPD',
        % Preservation of principal direction method
        dt6_new = dtiXformTensorsPPD(dt6_new, xformImgToTal);
    end
end
% Apply cannonical transform
% xformToTensor is same as inv(xformToTensor)
dt6_new = dtiXformTensors(dt6_new, xformToTensor);

return;

repo: vistasoft


text:
% SCIP  Solve a LP/MILP/QP/MIQP/QCQP/MIQCQP/NLP/MINLP using SCIP
%
% SCIP uses the Solving Constraint Integer Programs library, available
% from: http://scip.zib.de/scip.shtml
%
%   [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, opts)
%
%   Input arguments*:
%       H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU])
%       f - linear objective vector
%       A - linear constraint matrix (sparse)
%       rl - linear constraint lhs
%       ru - linear constraint rhs
%       lb - decision variable lower bounds
%       ub - decision variable upper bounds
%       xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary)
%       sos - SOS structure with fields type, index and weight (see below)
%       qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below)
%       nl - Nonlinear Objective and Constraints structure (see below)
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       stats - statistics structure
%
%   Option Fields (all optional - also see scipset):
%       tolrfun - LP primal convergence tolerance
%       maxiter - maximum LP solver iterations
%       maxnodes - maximum nodes to explore
%       maxtime - maximum execution time [s]
%       display - solver display level [0-5]
%       objbias - constant objective bias term
%
%   Return Status:
%       0 - Unknown
%       1 - User Interrupted
%       2 - Node Limit Reached
%       3 - Total Node Limit Reached
%       4 - Stall Node Limit Reached
%       5 - Time Limit Reached
%       6 - Memory Limit Reached
%       7 - Gap Limit Reached
%       8 - Solution Limit Reached
%       9 - Solution Improvement Limit Reached
%      10 - Restart Limit Reached
%      11 - Problem Solved to Optimality
%      12 - Problem is Infeasible
%      13 - Problem is Unbounded
%      14 - Problem is Either Infeasible or Unbounded
%
%   Special Ordered Sets (SOS):
%       type    - A string containing '1' or '2' for each SOS1 or SOS2 constraint.
%       index   - A double array of the indices of the variables in the SOS,
%                 group multiple SOS index arrays in a cell array.
%       weight  - A double array of the weights of each of the variable
%                 above, indicating variables next to each other. Group
%                 multiple SOS via cell arrays as above.
%
%   Quadratic Constraints (QC) [qrl <= x'Qx + l'x <= qru]:
%       Q       - A sparse double matrix of the quadratic terms for the
%                 constraint, group multiple quadratic constraints via a
%                 cell array of matrices. [NOT TRIL / TRIU]
%       l       - A column vector of the linear terms for the constraint,
%                 group multiple quadratic constraints in a matrix, each
%                 column representing each constraint.
%       qrl     - A scalar representing the quadratic constraint lower
%                 bound. Group multiple quadratic constraints in a
%                 column vector, each row representing each constraint.
%       qru     - A scalar representing the quadratic constraint upper
%                 bound. Group multiple quadratic constraints in a
%                 column vector, each row representing each constraint.
%
%   Nonlinear Objective & Constraints (min f(x) and cl <= c(x) <= cu)
%       See opti_scipnl for usage. Contains multiple fields of instructions
%       lists which the MEX interface parses into SCIP expressions.
%
%   * From OPTI v1.78 scip can also solve AMPL .nl models. Simply pass the FULL path
%	  of the AMPL .nl file to scip as the first argument, and options (if requred) as
%	  the second argument. This interface will also write an AMPL solution file to the
%	  same directory as the input file. Remember SCIP does not solve models with
%	  trigonometric functions.
%
%   Copyright (C) 2012-2013 Jonathan Currie (IPL)


text:
classdef math_model_opf_acps < mp.math_model_opf_acp
%MP.MATH_MODEL_OPF_ACPS  MATPOWER mathematical model for AC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_ACPS ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function tag = form_tag(obj)
            tag = 'acps';
        end

        function name = form_name(obj)
            name = 'AC-polar-power';
        end

        function add_node_balance_constraints(obj, nm, dm, mpopt)
            %% power balance constraints
            nn = nm.node.N;             %% number of nodes
            fcn_mis = @(x)obj.nodal_power_balance_fcn(x, nm);
            hess_mis = @(x, lam)obj.nodal_power_balance_hess(x, lam, nm);
            obj.add_nln_constraint({'Pmis', 'Qmis'}, [nn;nn], 1, fcn_mis, hess_mis);
        end

        function [lam_p, lam_q] = node_power_balance_prices(obj, nm)
            %% shadow prices on node power balance
            nne = obj.get_idx('nle');
            lambda = obj.soln.lambda;
            lam_p = lambda.eqnonlin(nne.i1.Pmis:nne.iN.Pmis);
            lam_q = lambda.eqnonlin(nne.i1.Qmis:nne.iN.Qmis);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [v, template] = struct2vector( s )

  if isnumeric(s)

    v = s(:);
    template = size(s);

  elseif iscell(s)

    [v, template] = cell2vector(s);

  else

    fields = fieldnames(s)';
    vs = {};
    shapes = {};
    for fi = 1:length(fields)
      x = getfield(s, fields{fi});
      if iscell(x)
        [vs{fi}, shapes{fi}] = cell2vector(x);
      elseif isstruct(x)
        [vs{fi}, shapes{fi}] = struct2vector(x);
      else
        vs{fi} = x(:);
        shapes{fi} = size(x);
      end
    end

    v = cat(1,vs{:});

    template.fields = fields;
    template.shapes = shapes;

  end
end


text:
function tn = tagnames(item, dflag)

% function tn = tagnames(item, dflag)
% Return the tags of all children in the job tree of an item. dflag
% indicates whether the filled (false) or defaults (true) part of the
% tree should be searched.
%
% This function is identical for all cfg_intree classes.
% It is not defined for leaf items.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: tagnames.m 4864 2012-08-27 13:57:31Z volkmar $

rev = '$Rev: 4864 $'; %#ok

tp = treepart(item, dflag);
citems = subsref(item, substruct('.', tp));
tn = cell(size(citems));
for k = 1:numel(citems)
    tn{k} = gettag(citems{k});
end;

author: spm


text:
function [params] = pf_calc(params)
% function [params] = pf_calc(params)
% -----------------------------------
% Calculation of the probabilistic forecast test.
%
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bRandom            Perform random simulation (=1) or real calculation (=0)
%   params.nCalculation       Number of random simulations
%   params.bMap               Calculate a map (=1) or a cross-section (=0)
%   params.bNumber            Use constant number (=1) or constant radius (=0)
%   params.nNumberEvents      Number of earthquakes if bNumber == 1
%   params.fRadius            Radius of gridnode if bNumber == 0
%   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
%   params.fForecastPeriod    Forecasting period in years
%   params.bLearningPeriod    Use different learning period than the rest of the catalog
%   params.fLearningPeriod    Learning period in years
%   params.bSignificance      Calculate significance during random simulation
%                             using params.fRealProbability
%   params.fRealProbability   Probability of calculation with real data
%   params.nCalculateMC       Method to calculate magnitude of completeness (see also: help calc_Mc)
%   params.nTestMethod        Method to calculate the Kagan & Jackson test (see also: help kj_poissonian)
%   params.bMinMagMc          Use magnitude of completeness as lower limit of magnitude range for testing (=1)
%                             Use params.fMinMag as lower limit (=0)
%   params.fMinMag            Lower limit of magnitude range for testing
%   params.fMaxMag            Upper limit of magnitude range for testing
%
% Output parameters:
%   Same as input parameters including
%   params.mValueGrid         Matrix of calculated Kagan & Jackson test values
%   params.vRandomMeans       Vector of means of probability differences per simulation run
%   params.vSignificanceLevel Vector of significance levels per simulation run
%   params.fBValueOverall     Calculated overall b-value
%   params.fStdDevOverall     Calculated standard deviation
%   params.fMcOverall         Calculated magnitude of completeness
%
% Danijel Schorlemmer
% April 24, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Init variable
% params.mRandomDist = [];

% Create random probability-ratios
% --------------------------------
if (params.bRandomNode) | (params.bRandomArea)
  if params.bForceRandomCalculation
    % Init the matrix of random probability-ratios (node-wise) for all calculations
    mRandomDist_ = [];
  else
    try
      mRandomDist_ = params.mRandomDist;
    catch
      mRandomDist_ = [];
    end
  end
  if isempty(mRandomDist_)
    % Loop over the calculations
    for nCnt_ = 1:params.nNumberCalculationNode
      % Init the vector of probability-ratios for one simulation
      vRandomDist_ = [];
      % Randomize catalog
      mRandomCatalog_ = params.mCatalog;
      mRandomCatalog_(:,6) = params.mCatalog(randperm(length(params.mCatalog)), 6);
      % Loop over the grid
      for nNode_ = 1:length(params.mPolygon(:,1))
        % Create catalog of given gridnode
        mNodeCatalog_ = mRandomCatalog_(params.caNodeIndices{nNode_}, :);
        % Calculate the probability-ratio
        [rCalcNodeResult_] = pf_CalcNode(mNodeCatalog_, params.fSplitTime, params.bLearningPeriod, params.fLearningPeriod, ...
          params.bForecastPeriod, params.fForecastPeriod, params.nCalculateMC, params.fMcOverall, params.bMinMagMc, ...
          params.fMinMag, params.fMaxMag, params.nTestMethod, params.nMinimumNumber, params.fBValueOverall, params.fStdDevOverall);
        % Store it into the vector
        vRandomDist_ = [vRandomDist_; rCalcNodeResult_.fProbDiff];
      end % of for nNode_
      % Store the vector into the matrix
      mRandomDist_ = [mRandomDist_ vRandomDist_];
    end % of for nCnt_
    params.mRandomDist = mRandomDist_;
  end
end % of if (params.bRandomNode) | (params.bRandomArea)

% Perform the real calculation
% ----------------------------
% Init result matrix
mValueGrid_ = [];
% Loop over all grid nodes
for nNode_ = 1:length(params.mPolygon(:,1))
  % Create node catalog
  mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
  % Calculate the probability-ratio
  [rCalcNodeResult_] = pf_CalcNode(mNodeCatalog_, params.fSplitTime, params.bLearningPeriod, params.fLearningPeriod, ...
    params.bForecastPeriod, params.fForecastPeriod, params.nCalculateMC, params.fMcOverall, params.bMinMagMc, ...
    params.fMinMag, params.fMaxMag, params.nTestMethod, params.nMinimumNumber, params.fBValueOverall, params.fStdDevOverall);
  % Calculate the significnance of the probability-ratio
  if params.bRandomNode
    fSignificanceLevel_ = kj_calcsig(rCalcNodeResult_.fProbDiff, mRandomDist_(nNode_,:), 2);
    fSignificanceLevel_ = fSignificanceLevel_ - 50;
    fSignificanceLevel_ = fSignificanceLevel_ * (-2); % Flip it (we want to be positive)
    fNormSignificanceLevel_ = kj_CalcNormSig(mRandomDist_(nNode_,:), rCalcNodeResult_.fProbDiff);
    fNormSignificanceLevel_ = fNormSignificanceLevel_ - 50;
    fNormSignificanceLevel_ = fNormSignificanceLevel_ * (-2); % Flip it (we want to be positive)
  else % of if params.bRandomNode
    fSignificanceLevel_ = nan;
    fNormSignificanceLevel_ = nan;
  end % of if params.bRandomNode
  % Store the results
  mValueGrid_= [mValueGrid_; fSignificanceLevel_ fNormSignificanceLevel_ rCalcNodeResult_.fProbDiff ...
      rCalcNodeResult_.fProbK rCalcNodeResult_.fProbO rCalcNodeResult_.nEventsLearning rCalcNodeResult_.nEventsObserved ...
      rCalcNodeResult_.fWeightK rCalcNodeResult_.fWeightO params.fBValueOverall rCalcNodeResult_.fBValueO rCalcNodeResult_.fMc];
end % for nNode
params.vcsGridNames = cellstr(char('Significance level', 'Significance level (normal distribution)', 'Probability difference', 'Kagan & Jackson', 'Our model', 'Number events in learning period', ...
  'Number events in forecasting period', 'Weighting of the overall b-value', 'Weighting of the node b-value', ...
  'b-value used for Kagan & Jackson model', 'b-value used for our model', 'Magnitude of completeness'));
params.mValueGrid = mValueGrid_;
% params.mRandomDist = mRandomDist_;

% % Perform the random simulation over the area
% % -------------------------------------------
% if params.bRandomArea
%   % Init matrix for significance values (n-times per gridnode)
%   mSignificanceGrid_ = [];
%   % Loop over the calculations
%   for nCnt_ = 1:params.nNumberCalculationArea
%     % Randomize catalog
%     mRandomCatalog_ = params.mCatalog;
%     mRandomCatalog_(:,6) = params.mCatalog(randperm(length(params.mCatalog)), 6);
%     % Init vector for significance values of one calculation over the entire grid
%     vSignificanceGrid_ = [];
%     % Loop over the grid
%     for nNode_ = 1:length(params.mPolygon(:,1))
%       % Create node catalog
%       mNodeCatalog_ = mRandomCatalog_(params.caNodeIndices{nNode_}, :);
%       % Calculate the probability-ratio
%       [rCalcNodeResult_] = pf_CalcNode(mNodeCatalog_, params.fSplitTime, params.bLearningPeriod, params.fLearningPeriod, ...
%         params.bForecastPeriod, params.fForecastPeriod, params.nCalculateMC, params.fMcOverall, params.bMinMagMc, ...
%         params.fMinMag, params.fMaxMag, params.nTestMethod, params.nMinimumNumber, params.fBValueOverall, params.fStdDevOverall);
%       % Calculate the significnance of the probability-ratio
%       fSignificanceLevel_ = kj_calcsig(rCalcNodeResult_.fProbDiff, mRandomDist_(nNode_,:), 2);
%       fSignificanceLevel_ = fSignificanceLevel_ - 50;
%       fSignificanceLevel_ = fSignificanceLevel_ * (-2); % Flip it (we want to be positive)
%       %       fNormSignificanceLevel_ = kj_CalcNormSig(mRandomDist_(nNode_,:), rCalcNodeResult_.fProbDiff);
%       %       fNormSignificanceLevel_ = fNormSignificanceLevel_ - 50;
%       %       fNormSignificanceLevel_ = fNormSignificanceLevel_ * (-2); % Flip it (we want to be positive)
%       % Store it into the vector
%       vSignificanceGrid_ = [vSignificanceGrid_ fSignificanceLevel_];
%     end % of for nNode_
%     % Store the vector into the matrix
%     mSignificanceGrid_ = [mSignificanceGrid_; vSignificanceGrid_];
%   end % of for nCnt_
% params.mSignificanceGrid = mSignificanceGrid_;
%
% % Evaluate overall significance
% % -----------------------------
%
% % Has to be set in the dialog
params.fSignificanceLevel = 95;
% % Compute the number of grid nodes with significant values
vSel_ = mValueGrid_(:,1) > params.fSignificanceLevel;
params.nRealPositive_ = sum(vSel_);
vSel_ = mValueGrid_(:,1) < -params.fSignificanceLevel;
params.nRealNegative_ = sum(vSel_);
% % Display them
disp(['Number positive: ' num2str(params.nRealPositive_)]);
disp(['Number negative: ' num2str(params.nRealNegative_)]);
%
% mRandomPositive_ = [];
% mRandomNegative_ = [];
%
% for nCnt_ = 1:params.nNumberCalculationArea
%  vSel_ = mSignificanceGrid_(nCnt_,:) > params.fSignificanceLevel;
%   %mTempGrid_ = mSignificanceGrid_(nCnt_,vSel_);
%   mRandomPositive_ = [mRandomPositive_; sum(vSel_)];
%
%  vSel_ = mSignificanceGrid_(nCnt_,:) < -params.fSignificanceLevel;
%   %mTempGrid_ = mSignificanceGrid_(nCnt_,vSel_);
%   mRandomNegative_ = [mRandomNegative_; sum(vSel_)];
% end
%
% fPositiveSignificance = kj_calcsig(nRealPositive_, mRandomPositive_, 3);
% fNegativeSignificance = kj_calcsig(nRealNegative_, mRandomNegative_, 3);
%
% disp(['Positive significance: ' num2str(fPositiveSignificance)]);
% disp(['Negative significance: ' num2str(fNegativeSignificance)]);
% end % of if params.bRandomArea

author: CelsoReyes