text:
function MapFig = createMapFig(Rob,Sen,Lmk, Trj, Frm, Fac,SimRob,SimSen,SimLmk,FigOpt)

% CREATEMAPFIG  Create 3D map figure and handles.
%   MAPFIG = CREATEMAPFIG(Rob,Sen,Lmk,SimRob,SimSen,SimLmk,MapFigure)
%   creates the map figure in figure 1, containing the following graphics
%   objects:
%       - a grid representing the ground
%       - simulated robots and sensors - in green
%       - estimated robots and sensors - in blue
%       - simulated landmarks (virtual world) - in red
%       - estimated landmarks, containing mean and uncertainty ellipsoid - in
%       colors depending on the landmark type.
%
%   The output MAPFIG is a structure of handles to all these graphics
%   objects. See the Matlab documentation for information about graphics
%   handles and the way to efficiently manipulate graphics. MAPFIG has the
%   following fields:
%       .fig    handle to the figure
%       .axes   handle to the axes
%       .simRob array of handles to the simulated robots 'patch' objects
%       .simSen array of handles to the simulated sensors 'patch' objects
%       .simLmk handle to the simulated landmarks 'line' object
%       .Rob array of handles to the estimated robots 'patch' objects
%       .Sen array of handles to the estimated sensors 'patch' objects
%       .Lmk array of structures with handles to the estimated landmarks, with
%           .mean       handle to the ellipsoid's center 'line' object
%           .ellipse    handle to the ellipsoid's contour 'line' object
%
%   The figure is updated using drawMapFig.
%
%   See also DRAWMAPFIG, CREATESENFIG, MAPOBSERVER, LINE, PATCH, SURFACE,
%   SET, GET.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

global Map

% Figure
if FigOpt.createVideo
    MapFig.fig = figure(99);
    set(MapFig.fig, 'WindowStyle', 'normal');
    figPos     = get(MapFig.fig,'position');
    figSize    = FigOpt.map.size;
    newFigPos  = [0 figPos(2)  figSize];
    set(MapFig.fig,'position',newFigPos);
else
    if ishandle(99)
        MapFig.fig = figure(99);
    else
        MapFig.fig = figure(99);
        if ~strcmp( get(MapFig.fig, 'WindowStyle'), 'docked')
            figPos     = get(MapFig.fig,'position');
            figSize    = FigOpt.map.size;
            newFigPos  = [0 figPos(2)  figSize];
            set(MapFig.fig,'position',newFigPos);
        end
    end
end
clf
moreindatatip

set(MapFig.fig,...
    'numbertitle',   'off',...
    'name',          '3D Map',...
    'doublebuffer',  'off',...
    'renderer',      FigOpt.renderer,...
    'toolbar',       'none',...
    'color',         FigOpt.map.colors.bckgnd);
cameratoolbar('show');
cameratoolbar('setmode','orbit');

% World dimensions
world.lims         = FigOpt.map.lims;
world.dims.l       = FigOpt.map.lims.xMax - FigOpt.map.lims.xMin;
world.dims.w       = FigOpt.map.lims.yMax - FigOpt.map.lims.yMin;
world.dims.h       = FigOpt.map.lims.zMax - FigOpt.map.lims.zMin;
world.center.xMean = (FigOpt.map.lims.xMax + FigOpt.map.lims.xMin)/2;
world.center.yMean = (FigOpt.map.lims.yMax + FigOpt.map.lims.yMin)/2;
world.center.zMean = (FigOpt.map.lims.zMax + FigOpt.map.lims.zMin)/2;

% Map viewpoint
viewPnt = mapObserver(world,FigOpt.map.view);

% Axes
axis equal
MapFig.axes = gca;
set(MapFig.axes,...
    'parent',              MapFig.fig,...
    'nextplot',            'replacechildren',...
    'position',            [ 0 0 1 1],...
    'SortMethod',          'childorder',... % {'depth'  'childorder'}
    'projection',          FigOpt.map.proj,...
    'CameraPosition',      viewPnt.X,...
    'CameraPositionMode',  'manual',...
    'CameraUpVector',      viewPnt.upvec,...
    'CameraUpVectorMode',  'manual',...
    'CameraViewAngle',     viewPnt.fov,...
    'CameraViewAngleMode', 'manual',...
    'CameraTarget',        viewPnt.tgt,...
    'CameraTargetMode',    'manual',...
    'xlim',                [FigOpt.map.lims.xMin FigOpt.map.lims.xMax],...
    'ylim',                [FigOpt.map.lims.yMin FigOpt.map.lims.yMax],...
    'zlim',                [FigOpt.map.lims.zMin FigOpt.map.lims.zMax],...
    'alimmode',            'manual',...
    'climmode',            'manual',...
    'vis',                 'off');%,...

% OBJECTS COMMON TO SIMULATION AND ESTIMATION
% Ground
MapFig.ground = createGround(FigOpt.map,MapFig.axes,FigOpt.map.colors.ground);

% ESTIMATED OBJECTS
% robots
for rob = 1:numel(Rob)

    % create and draw robot - with ellipse
    MapFig.Rob(rob).patch = createObjPatch(Rob(rob),FigOpt.map.colors.est,MapFig.axes);
    MapFig.Rob(rob).ellipse = line(...
        'parent', MapFig.axes,...
        'xdata',  [],    ...
        'ydata',  [],    ...
        'zdata',  [],    ...
        'color',  'r',   ...
        'marker', 'none');

    % sensors
    for sen = Rob(rob).sensors

        % create and draw sensor
        MapFig.Sen(sen) = createObjPatch(Sen(sen),FigOpt.map.colors.est,MapFig.axes);

        % redraw sensor in robot frame
        F = composeFrames(Rob(rob).frame,Sen(sen).frame);
        drawObject(MapFig.Sen(sen),Sen(sen),F);

    end

    % trajectory
    n = size(Frm,2);
    Z = zeros(2,n);
    MapFig.Rob(rob).trj = line(...
        Z,    ...
        Z,    ...
        Z,    ...
        'parent', MapFig.axes,...
        'visible','off', ...
        'color',  FigOpt.map.colors.graph.motion,   ...
        'marker', 'o');

    % measurement factors
    n = numel(Fac);
    Z = zeros(2,n);
    MapFig.Rob(rob).factors = line(...
        Z,    ...
        Z,    ...
        Z,    ...
        'parent', MapFig.axes,...
        'visible','off', ...
        'color',  FigOpt.map.colors.graph.meas,   ...
        'marker', 'none');

end

% landmarks
for lmk = 1:numel(Lmk)

    MapFig.Lmk(lmk) = createLmkGraphics(Lmk(lmk),FigOpt.map.colors.label,MapFig.axes);

end

% SIMULATED OBJECTS
if ~isempty(SimRob) && ~isempty(SimSen) && ~isempty(SimLmk)

    % Landmarks - do not loop, draw all at once
    MapFig.simLmk = createSimLmkGraphics(SimLmk,FigOpt.map.colors.simLmk,MapFig.axes,FigOpt.map.showSimLmk);

    % Robots
    for rob = 1:numel(SimRob)

        % create and draw robot
        MapFig.simRob(rob) = createObjPatch(SimRob(rob),FigOpt.map.colors.simu,MapFig.axes);

        % Sensors
        for sen = SimRob(rob).sensors

            % create and draw sensor
            MapFig.simSen(sen) = createObjPatch(SimSen(sen),FigOpt.map.colors.simu,MapFig.axes);

            % redraw sensor in robot frame
            F = composeFrames(SimRob(rob).frame,SimSen(sen).frame);
            drawObject(MapFig.simSen(sen),SimSen(sen),F);
        end
    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function planC = sumDose(doseNums,wtfactor,lqParamS,assocScan,newDoseName,planC)
%function planC = sumDose(doseNumsV,weightsV,assocScan)
%
%This function creates a new dose distribution by adding up doseNums
%according to factors wtfactor. assocScan is the scan number to associate
%new dose with. newDoseName is the name of new dose.
%
%APA, 12/23/2009
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

wtfactorTmp = wtfactor;
wtfactor = ones(1,max(doseNums));
wtfactor(doseNums) = wtfactorTmp;

indexS = planC{end};
if assocScan > 0
    assocScanUID = planC{indexS.scan}(assocScan).scanUID;
    %Get associated transM
    assocTransM = planC{indexS.scan}(assocScan).transM;
    if isempty(assocTransM)
        assocTransM = eye(4);
    end
else %No Association
    assocScanUID = '';
    assocTransM = eye(4);
end

% doseNums = checkedDoses;

%Get the x,y,z grid for new dose
for i = 1:length(doseNums)
    doseNum = doseNums(i);
    %Get x,y,z values for doseNum
    [xV, yV, zV] = getDoseXYZVals(planC{indexS.dose}(doseNum));

    %Get the corners of the original dataset.
    [xCorn, yCorn, zCorn] = meshgrid([min(xV) max(xV)], [min(yV) max(yV)], [min(zV) max(zV)]);

    %Add ones to the corners so we can apply a transformation matrix.
    corners = [xCorn(:) yCorn(:) zCorn(:) ones(prod(size(xCorn)), 1)];

    %Apply transform to corners, so we know boundary of the slice.
    transM = getTransM(planC{indexS.dose}(doseNum),planC);
    if isempty(transM) || isequal(transM,eye(4))
        [xV,yV,zV] = getDoseXYZVals(planC{indexS.dose}(doseNum));
        xGrid{doseNum} = xV(:)';
        yGrid{doseNum} = yV(:)';
        zGrid{doseNum} = zV(:)';
    else
        newCorners = inv(assocTransM) * transM * corners';
        xGrid{doseNum} = linspace(min(newCorners(1,:)), max(newCorners(1,:)), length(xV));
        yGrid{doseNum} = linspace(max(newCorners(2,:)), min(newCorners(2,:)), length(yV));
        zGrid{doseNum} = linspace(min(newCorners(3,:)), max(newCorners(3,:)), length(zV));
    end
    xRes{doseNum} = length(xV);
    yRes{doseNum} = length(yV);
    zRes{doseNum} = length(zV);

    %Get associated scan
    % assocScanV{doseNum} = getAssociatedScan(planC{indexS.dose}(doseNum).assocScanUID);
    assocScanV{doseNum} = getDoseAssociatedScan(doseNum,planC);
end
newXgrid = linspace(min(cell2mat(xGrid)),max(cell2mat(xGrid)),max(cell2mat(xRes)));
newYgrid = linspace(max(cell2mat(yGrid)),min(cell2mat(yGrid)),max(cell2mat(yRes)));
newZgrid = linspace(min(cell2mat(zGrid)),max(cell2mat(zGrid)),max(cell2mat(zRes)));

%Obtain doses with same grid
if isempty(assocTransM)
    assocTransM = eye(4);
end
doseIndC = {};
doseNumsTmp = doseNums;
for iSortAll = 1:length(doseNums)
    iSort = doseNums(iSortAll);
    indRemaining = doseNums;
    indRemaining(iSortAll) = [];
    doseSortM = [iSort];
    if ~ismember(iSort,[doseIndC{:}])
        for jSortAll = 1:length(indRemaining)
            jSort = indRemaining(jSortAll);
            if ~isempty(getTransM('dose',iSort,planC))
                doseITM = inv(assocTransM) * getTransM('dose',iSort,planC);
            else
                doseITM = inv(assocTransM);
            end
            if ~isempty(getTransM('dose',jSort,planC))
                doseJTM = inv(assocTransM) * getTransM('dose',jSort,planC);
            else
                doseJTM = inv(assocTransM);
            end
            if isequal([xRes{iSort},yRes{iSort},zRes{iSort}],[xRes{jSort},yRes{jSort},zRes{jSort}]) && isequal(doseITM,doseJTM)
                doseSortM(end+1) = jSort;
                indJsort = find(doseNumsTmp == jSort);
                doseNumsTmp(indJsort) = [];
            end
        end
    end
    doseIndC{iSort} = doseSortM;
end

doseNums = doseNumsTmp;

%Loop over doses and add over new grid
hWait = waitbar(0,'Summing Dose distributions');
doseSumM = zeros([length(newYgrid),length(newXgrid),length(newZgrid)],'single');
doseEmptyM = zeros([length(newYgrid),length(newXgrid)],'single');

%Assume dose units are same as that of 1st dose
doseUnits = getDoseUnitsStr(doseNums(1),planC);
for i = 1:length(doseNums)

    doseNum = doseNums(i);

    %Check for transM
    if ~isempty(assocTransM) && ~isequal(assocTransM,eye(4))
        doseTtransM = getTransM('dose',doseNum,planC);
        if isempty(doseTtransM)
            doseTtransM = eye(4);
        end
        inputTM = inv(assocTransM) * doseTtransM;
    else
        inputTM = getTransM('dose',doseNum,planC);
        if isempty(inputTM)
            inputTM = eye(4);
        end
    end

    SOPInstanceUIDv = {planC{indexS.beams}.SOPInstanceUID};
    if ~isempty(lqParamS)
        paramS.Tk.val = inf;         %Kick-off time of repopulation (days)
        paramS.Tp.val = NaN;        %Potential tumor doubling time (days)
        paramS.alpha.val = NaN;
        paramS.abRatio.val = lqParamS.abRatio; %10;  %alpha/beta
        paramS.stdFractionSize = lqParamS.stdFractionSize; % 2;
    end

    %Get the summation for this grid
    doseCombinedM = [];
    for iDoseAll = 1:length(doseIndC{doseNum})
        iDose = doseIndC{doseNum}(iDoseAll);
        doseUnits2 = getDoseUnitsStr(iDose,planC);
        if strcmpi(doseUnits, 'Gy') && strcmpi(doseUnits2, 'cGy')
            multFact = 0.01;
        elseif strcmpi(doseUnits, 'cGy') && strcmpi(doseUnits2, 'Gy')
            multFact = 100;
        else
            multFact = 1;
        end
        doseOffset = planC{indexS.dose}(iDose).doseOffset;
        if isempty(doseOffset)
            doseOffset = 0;
        end
        doseArray = single(getDoseArray(planC{indexS.dose}(iDose)) - doseOffset);

        % Apply BED/EQD2 correction
        if ~isempty(lqParamS)
            ReferencedSOPInstanceUID = planC{indexS.dose}(iDose)...
                .DICOMHeaders.ReferencedRTPlanSequence.Item_1.ReferencedSOPInstanceUID;
            planNum = find(strcmpi(ReferencedSOPInstanceUID,SOPInstanceUIDv));
            paramS.numFractions.val = planC{indexS.beams}(planNum(1)).FractionGroupSequence...
                .Item_1.NumberOfFractionsPlanned;
            paramS.numFractions.val = double(paramS.numFractions.val);
            paramS.frxSize.val = doseArray / paramS.numFractions.val;
            doseArray = calc_BED(paramS) / (1+paramS.stdFractionSize/paramS.abRatio.val);
        end

        if ~isempty(doseCombinedM)
            doseCombinedM = doseCombinedM + multFact * wtfactor(iDose) * doseArray;
        else
            doseCombinedM = multFact * wtfactor(iDose) * doseArray;
        end
    end

    %Check if this dose is on same grid as the new one
    %if isequal(xGrid{doseNum},newXgrid) && isequal(yGrid{doseNum},newYgrid) && isequal(zGrid{doseNum},newZgrid)
    chkLength = length(xGrid{doseNum})==length(newXgrid) && length(yGrid{doseNum})==length(newYgrid) && length(zGrid{doseNum})==length(newZgrid);
    if sum(sum((inputTM-eye(4)).^2)) < 1e-3 && chkLength && max(abs(xGrid{doseNum}-newXgrid)) < 1e-3 && max(abs(yGrid{doseNum}-newYgrid)) < 1e-3 && max(abs(zGrid{doseNum}-newZgrid)) < 1e-3

        doseSumM = doseSumM + doseCombinedM;
        waitbar((i-1)/length(doseNums),hWait,['Calculating contribution from Dose ', num2str(doseNum)])

    else %interpolation required

        %Transform this dose
        doseTmpM = [];

        for slcNum=1:length(newZgrid)
            [xV, yV, zV] = getDoseXYZVals(planC{indexS.dose}(doseNum));
            doseTmp = slice3DVol(doseCombinedM, xV, yV, zV, newZgrid(slcNum), 3, 'linear', inputTM, [], newXgrid, newYgrid);
            if isempty(doseTmp)
                doseTmpM(:,:,slcNum) = doseEmptyM;
            else
                doseTmpM(:,:,slcNum) = doseTmp;
            end
            waitbar((i-1)/length(doseNums) + (slcNum-1)/length(newZgrid)/length(doseNums) ,hWait,['Calculating contribution from Dose ', num2str(doseNum)])
        end

        doseSumM = doseSumM + doseTmpM;

    end

end

delete(hWait)

%Create new dose distribution
newDoseNum = length(planC{indexS.dose}) + 1;
planC{indexS.dose}(newDoseNum).doseArray = doseSumM;
clear doseSumM
planC{indexS.dose}(newDoseNum).doseUID = createUID('dose');
planC{indexS.dose}(newDoseNum).assocScanUID = assocScanUID;

%Find minimum value in 3d array, use its negative as the offset
maxDose = max(max(max(planC{indexS.dose}(newDoseNum).doseArray)));
offset = -min(min(min(planC{indexS.dose}(newDoseNum).doseArray)));
if offset > 0
    planC{indexS.dose}(newDoseNum).doseOffset = offset;
    planC{indexS.dose}(newDoseNum).doseArray = planC{indexS.dose}(newDoseNum).doseArray + offset;
end

%set labels on new dose, overwriting some of the copied labels **Check for more labels that need to be replaced
planC{indexS.dose}(newDoseNum).doseNumber = newDoseNum;
planC{indexS.dose}(newDoseNum).fractionGroupID = newDoseName;

%Remove old caching info.
planC{indexS.dose}(newDoseNum).cachedMask = [];
planC{indexS.dose}(newDoseNum).cachedColor = [];
planC{indexS.dose}(newDoseNum).cachedTime = [];

%Set coordinates.
planC{indexS.dose}(newDoseNum).sizeOfDimension1 = length(newXgrid);
planC{indexS.dose}(newDoseNum).sizeOfDimension2 = length(newYgrid);
planC{indexS.dose}(newDoseNum).sizeOfDimension3 = length(newZgrid);
planC{indexS.dose}(newDoseNum).horizontalGridInterval = newXgrid(2)-newXgrid(1);
planC{indexS.dose}(newDoseNum).verticalGridInterval = newYgrid(2)-newYgrid(1);
planC{indexS.dose}(newDoseNum).depthGridInterval = newZgrid(2)-newZgrid(1);
planC{indexS.dose}(newDoseNum).coord1OFFirstPoint = newXgrid(1);
planC{indexS.dose}(newDoseNum).coord2OFFirstPoint = newYgrid(1);
planC{indexS.dose}(newDoseNum).coord3OfFirstPoint = newZgrid(1);
planC{indexS.dose}(newDoseNum).zValues = newZgrid;
planC{indexS.dose}(newDoseNum).doseUnits = doseUnits;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%	nufft_tune_gauss
%	tune the FWHM for gaussian interpolation in 1D NUFFT
%	to minimize the worst-case error

if ~isvar('err_gauss')
	N = 2^8;
	K_N = 2;				% oversampling factor
	Jlist = [2:16];
	Plist = linspace(0.2, 0.4, 101);	% sig = parameter * sqrt(J)
	rlist = [0:20]'/40;			% fractions of gamma

	err_gauss.zn = zeros([length(rlist) length(Jlist) length(Plist)]);
	err_gauss.ft = zeros([length(rlist) length(Jlist) length(Plist)]);

	for jj=1:length(Jlist)
		J = Jlist(jj);	printf('J=%d', J)

		K = K_N * N;
		gam = 2*pi/K;
		om = gam * rlist;

		%
		%	gaussian with various FWHM
		%
		for ii=1:length(Plist)
			sig = Plist(ii) * sqrt(J);
			[Gfunc, Gfunc_ft] = nufft_gauss('inline', J, sig);

			%
			%	interpolator worst-case error
			%
		        err_gauss.ft(:,jj,ii) = nufft1_error(om, N, J, K, ...
				Gfunc, Gfunc_ft);
		        err_gauss.zn(:,jj,ii) = nufft1_error(om, N, J, K, ...
				Gfunc);
		end
	end
end

emax.zn = squeeze( max(err_gauss.zn, [], 1) );		% [M,J,F] -> [J,F]
[emax.zn, ibest.zn] = min(emax.zn, [], 2);		% [J,F] -> J
pbest.zn = Plist(ibest.zn);

emax.ft = squeeze( max(err_gauss.ft, [], 1) );		% [M,J,F] -> [J,F]
[emax.ft, ibest.ft] = min(emax.ft, [], 2);		% [J,F] -> J
pbest.ft = Plist(ibest.ft);

%	plot best sig vs J (on log scale!)
clf
if 1
	log_J = log(Jlist);
	log_sig.zn = log(pbest.zn .* sqrt(Jlist));
	log_sig.ft = log(pbest.ft .* sqrt(Jlist));

	subplot(131)
	plot(log_J, log_sig.zn, 'co', log_J, log_sig.ft, 'yx'), axis tight
	lsline, legend('zn', 'ft')
	xlabel log(J), ylabel 'log(sig best)'
	pol = polyfit(log_J, log_sig.zn, 1)
	printf('\\sigma \\approx %g * J^{%g}', exp(pol(2)), pol(1))
	pol = polyfit(log_J, log_sig.ft, 1)
	printf('\\sigma \\approx %g * J^{%g}', exp(pol(2)), pol(1))
end

%
%	plot emax vs J
%
if 1
	subplot(132)
	semilogy(Jlist, emax.zn, 'c-o', Jlist, emax.ft, 'y-x')
	xlabel J, ylabel E_{max}
        title(sprintf('Maximum error for K/N=%g', K/N))
	legend('zn', 'ft')

	subplot(133)
	plot(emax.zn ./ emax.ft), axis tight, title 'emax: zn / ft'
%	savefig c 'fig_?'
end

%
%	save both sets of tuned sigma's to "private" file
%
if 1
	Jgauss2 = Jlist;
	Sgauss2.zn = pbest.zn .* sqrt(Jlist);		% tuned sig zn
	Sgauss2.ft = pbest.ft .* sqrt(Jlist);		% tuned sig ft
	save nufft_gauss2 Jgauss2 Sgauss2
	disp('MOVE TO private!')
end

author: thomaskuestner


text:
function results = vl_test_kdtree(varargin)
% VL_TEST_KDTREE
vl_test_init ;

function s = setup()
randn('state',0) ;
s.X = single(randn(10, 1000)) ;
s.Q = single(randn(10, 10)) ;

function test_nearest(s)
for tmethod = {'median', 'mean'}
  for type = {@single, @double}
    conv = type{1} ;
    tmethod = char(tmethod) ;

    X = conv(s.X) ;
    Q = conv(s.Q) ;
    tree = vl_kdtreebuild(X,'ThresholdMethod', tmethod) ;
    [nn, d2] = vl_kdtreequery(tree, X,Q) ;

    D2 = vl_alldist2(X, Q, 'l2') ;
    [d2_, nn_] = min(D2) ;

    vl_assert_equal(...
      nn,uint32(nn_),...
      'incorrect nns: type=%s th. method=%s', func2str(conv), tmethod) ;
    vl_assert_almost_equal(...
      d2,d2_,...
      'incorrect distances: type=%s th. method=%s', func2str(conv), tmethod) ;
  end
end

function test_nearests(s)
numNeighbors = 7 ;
tree = vl_kdtreebuild(s.X) ;
[nn, d2] = vl_kdtreequery(tree, s.X, s.Q, ...
                          'numNeighbors', numNeighbors) ;

D2 = vl_alldist2(s.X, s.Q, 'l2') ;
[d2_, nn_] = sort(D2) ;
d2_ = d2_(1:numNeighbors, :) ;
nn_ = nn_(1:numNeighbors, :) ;

vl_assert_equal(nn,uint32(nn_)) ;
vl_assert_almost_equal(d2,d2_) ;

function test_ann(s)
vl_twister('state', 1) ;
numNeighbors = 7 ;
maxComparisons = numNeighbors * 50 ;
tree = vl_kdtreebuild(s.X) ;
[nn, d2] = vl_kdtreequery(tree, s.X, s.Q, ...
                          'numNeighbors', numNeighbors, ...
                          'maxComparisons', maxComparisons) ;

D2 = vl_alldist2(s.X, s.Q, 'l2') ;
[d2_, nn_] = sort(D2) ;
d2_ = d2_(1:numNeighbors, :) ;
nn_ = nn_(1:numNeighbors, :) ;

for i=1:size(s.Q,2)
  overlap = numel(intersect(nn(:,i), nn_(:,i))) / ...
            numel(union(nn(:,i), nn_(:,i))) ;
  assert(overlap > 0.6, 'ANN did not return enough correct nearest neighbors') ;
end

function test_ann_forest(s)
vl_twister('state', 1) ;
numNeighbors = 7 ;
maxComparisons = numNeighbors * 25 ;
numTrees = 5 ;
tree = vl_kdtreebuild(s.X, 'numTrees', 5) ;
[nn, d2] = vl_kdtreequery(tree, s.X, s.Q, ...
                          'numNeighbors', numNeighbors, ...
                          'maxComparisons', maxComparisons) ;

D2 = vl_alldist2(s.X, s.Q, 'l2') ;
[d2_, nn_] = sort(D2) ;
d2_ = d2_(1:numNeighbors, :) ;
nn_ = nn_(1:numNeighbors, :) ;

for i=1:size(s.Q,2)
  overlap = numel(intersect(nn(:,i), nn_(:,i))) / ...
            numel(union(nn(:,i), nn_(:,i))) ;
  assert(overlap > 0.6, 'ANN did not return enough correct nearest neighbors') ;
end

author: Cloud-CV

repo: object-proposals


text:
function linplus_test154 ( )

%*****************************************************************************80
%
%% TEST154 tests R8BTO_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  l = 3;
  m = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST154\n' );
  fprintf ( 1, '  For a real block Toeplitz matrix,\n' );
  fprintf ( 1, '  R8BTO_INDICATOR sets up an indicator matrix\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Block order M =  %d\n', m );
  fprintf ( 1, '  Block number L = %d\n', l );
  fprintf ( 1, '  Matrix order N = %d\n', m * l );

  a = r8bto_indicator ( m, l );

  r8bto_print ( m, l, a, '  The block Toeplitz matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ImgOut = fftsf(ImgIn,degree)
% FFTSF Frequency domain auto filtering (enhances powerfull frequencies)
%
% Input :	ImgIn  - grayscale image
%             degree - emphasis given to regular patterns
%                    values: 1, 2, 3, {4} (default)
%
%                     Filtering    |   soft    strong
%                     ------------------------------
%                     emphasis     |   1       2
%                     attenuation  |   3       4
%
% Output: class double data
% Display:
%
% im = imread('kb818c.jpg');
% figure
% subplot(2,3,2)
% imshow(im,[])
% title('Input')
% subplot(2,3,3)
% imshow(fftsf(im,1),[])
% title('1. Soft emphazis')
% subplot(2,3,6)
% imshow(fftsf(im,2),[])
% title('2. Strong emphazis')
% subplot(2,3,5)
% title([{'Self-filtering'},{'emphasizes or attenuates'},{'regular patterns'}])
% axis off
% subplot(2,3,1)
% imshow(fftsf(im,3),[])
% title('3. Soft attenuation')
% subplot(2,3,4)
% imshow(fftsf(im,4),[])
% title('4. Strong attenuation')
%
% Proceedure: The most powerfull frequencies in an image are enhanced /
%   attenuated by multiplying / dividing the complex FFT signal with the
%   magnitude of the same FFT.
%
% Reference: This code is based on the technique described in:
%   D.G. Bailey - Detecting regular patterns using frequency domain
%   self-filtering, 1997 Intl. Conf. on Image Processing, 1:440-3.
%
% Written by Vlad Atanasiu 2003.09.25
% 2008.11.07 - modified sample code
% 2007.07.07 - added new magnitude enhancement options
%            - removed bug in 'switch'
%            - changed name from 'fftselffilter' to 'fftsf' (is shorter)

F = fft2(double(ImgIn));	% FFT
R = abs(F);                     % magnitude
Z = angle(F);                   % phase
switch degree                  % filtering
case(2)
    F2 = F.*R.*sqrt(R.^2+Z.^2);
case(3)
    R = R + eps;
    F2 = sqrt(R).*exp(i*Z);
case(4)
    R = R + eps;
    F2 = F./R;
otherwise
    F2 = F.*R;
end
ImgOut = real(ifft2(F2));		% filtered image

author: Sable

repo: mcbench-benchmarks


text:
function [mix] = spm_mix (y,m,verbose)
% Fit a multivariate Gaussian Mixture model using VB
% FORMAT [mix] = spm_mix (y,m,verbose)
%
% y          [N x d] data matrix containing N samples of d-dim data
% m          Number of mixture components
% verbose    Set to 1 to see evolution of free energy, 0 otherwise
%            (default=1)
%
% mix        Returned model
%
%--------------------------------------------------------------------------
% The fields in mix are:
%
% m                The number of components
% fm               The negative free energy. This decomposes into
%                  fm=acc-kl_proportions-kl_covs-kl_centres
%
% acc              model accuracy
% kl_proportions   complexity penalty for cluster proportions
% kl_covs          complexity penalty for cluster covariances
% kl_centres       complexity penalty for cluster centres
%
% Fields:
%
% lambda           Post mixers, q(pi|D) = D(lambda)
% gamma            [m x N] matrix of belonging probabilities
% state(s).a       Post precisions, q(Gamma|D)=W(a,B)
% state(s).B
% state(s).C       Post covariance
% state(s).m       Post mean, q(mu|D)=N(m_s,beta_s Gamma_s)
% state(s).beta
% state(s).prior   Estimated mixing proportions
%
%                  In the field prior:
%
% lambda_0         Prior mixers, p(pi) = D(lambda_0)
% a_0,B_0          Prior precisions, p(Gamma)=W(a_0,B_0)
% m_0,beta_0       Prior means, p(mu)=N(m_0,beta_0 Gamma_s)
%
%__________________________________________________________________________
% Copyright (C) 2007-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mix.m 5962 2014-04-17 12:47:43Z spm $

% This code implements the algorithm in:
%
% See Attias, H. (2000) A Variational
% Bayesian framework for Graphical Models, NIPS.
%
% See also:
%
% W.D. Penny (2001) Variational Bayes for d-dimensional Gaussian mixture models
% Wellcome Department of Imaging Neuroscience, University College London.
%
% for Negative Free Energy expression and validation.

if nargin < 3
    verbose=1;
end
[N,d]=size(y);

% Put model info into data structure format
mix.nin=d;
mix.m=m;

% PRIORS
% Set mixing priors
lambda_0=1;
% Set mean priors
m_0=mean(y);
beta_0=1;
% Set precision priors
% These settings set the prior covariance matrix to be
% equal to the 0.01*d*identity matrix - if original data is
% zmuv, this is reasonable
a_0=d;
B_0=0.01*d*eye(d);
%B_0=d*eye(d);

% Also, store in data structure ready for output
mix.prior.lambda_0=lambda_0;
mix.prior.m_0=m_0;
mix.prior.beta_0=beta_0;
mix.prior.a_0=a_0;
mix.prior.B_0=B_0;

sd=0;
for i=1:d,
  sd=sd+psi((a_0+1-i)/2);
end
log_tilde_gamma_0=sd-log(det(B_0))+d*log(2);

if (m==1)
    % For a single component we just have a Gaussian model

    % Set posterios
    state(1).a=N+a_0;
    state(1).beta=N+beta_0;
    state(1).m=m_0;
    Cy=cov(y);
    state(1).B=N*Cy+B_0;

    fm=spm_lg_gamma (d,0.5*state(1).a);
    fm=fm-spm_lg_gamma (d,0.5*a_0);
    fm=fm-0.5*d*N*log(pi);
    fm=fm+0.5*d*(log(beta_0)-log(state(1).beta));
    fm=fm+0.5*a_0*log(det(B_0));
    fm=fm-0.5*state(1).a*log(det(state(1).B));

    mix.m=m;
    mix.fm=fm;
    state(1).C=state(1).B/state(1).a;
    mix.state=state;

    mix.priors=1;
    return;
end

% Run kmeans on the data
% to initialise posterior
[priors,means,covs] = spm_kmeans(y,m);

% Add pseudo-counts to ML priors
priors=priors+(1/N);
lambda=priors;
for s=1:m,
  state(s).m=means(s,:);
  state(s).beta=priors(s)*N+beta_0;
  state(s).a=priors(s)*N+a_0;
  Cov=covs(:,:,s);
  % If component cov is rank deficient set to prior
  if rank(Cov) < size(Cov,1)
    state(s).B=B_0;
  else
    state(s).B=priors(s)*N*Cov;
  end
  N_bar(s)=priors(s)*N;
end

% Start algorithm
lik=[];
tol=0.0001;
max_loops=32;
for loops=1:max_loops,

    % E-step
    lambda_tot=sum(lambda);
    for s=1:m,
        state(s).bar_gamma=state(s).a*inv(state(s).B);
        log_tilde_pi(s)=psi(lambda(s))-psi(lambda_tot);
        sd=0;
        for i=1:d,
            sd=sd+psi((state(s).a+1-i)/2);
        end
        log_tilde_gamma(s)=sd-log(det(state(s).B))+d*log(2);

        tilde_pi(s)=exp(log_tilde_pi(s));
        tilde_gamma(s)=exp(log_tilde_gamma(s));
        for n=1:N,
            gamma(s,n)=tilde_pi(s)*tilde_gamma(s)^0.5;
            dy=(y(n,:)-state(s).m);
            gamma(s,n)=gamma(s,n)*(exp(-0.5*dy*state(s).bar_gamma*dy')+eps)*exp(-d/(2*state(s).beta));
        end
    end

    gamma_n=sum(gamma);
    for s=1:m,
        if mean(gamma_n) > eps
            % If component still exists
            gamma(s,:)=gamma(s,:)./gamma_n;
        end
    end

    % M-step

    % Part-I
    for s=1:m,
        pi_bar(s)=mean(gamma(s,:))+eps;
        N_bar(s)=N*pi_bar(s)+eps;
        bar_mu(s,:)=(1/N_bar(s))*sum(gamma(s,:)'*ones(1,d).*y);

    end

    % get weighted means and covariances
    for s=1:m,
        state(s).bar_sigma=zeros(d,d);
        for n=1:N,
            dy=y(n,:)-bar_mu(s,:);
            state(s).bar_sigma=state(s).bar_sigma+gamma(s,n).*(dy'*dy);
        end
        state(s).bar_sigma=(state(s).bar_sigma)/N_bar(s);

        if isnan(state(s).bar_sigma)
            state(s).bar_sigma
        end
    end

    % Now compute the free energy
    f1=-spm_kl_dirichlet(lambda,lambda_0*ones(1,m),log_tilde_pi);
    for s=1:m,
        f2(s)=-spm_kl_wishart(state(s).a,state(s).B,a_0,B_0);

        % KL-method for computing f3(s)
        Cs=state(s).B/(state(s).beta*state(s).a);
        C0=state(s).B/(beta_0*state(s).a);
        f3(s)=-spm_kl_normal(state(s).m,Cs,m_0,C0);

        % Check f3(s)
        check_f3=-0.5*d*(log(state(s).beta)-log(beta_0)+(beta_0/state(s).beta)-1);
        dm=state(s).m-m_0;
        check_f3=check_f3-0.5*dm*beta_0*state(s).a*inv(state(s).B)*dm';
        f3(s)=check_f3;

        f4(s)=N_bar(s)*log_tilde_pi(s)-sum(gamma(s,:).*log(gamma(s,:)+eps));
        LaB=0;
        for i=1:d,
            LaB=LaB+psi((state(s).a+1-i)/2);
        end
        LaB=LaB+d*log(2)-log(det(state(s).B));
        f5(s)=0.5*N_bar(s)*(-d*log(2*pi)+LaB-trace(state(s).bar_gamma*state(s).bar_sigma)-(d/state(s).beta));
        fkl(s)=f2(s)+f3(s)+f4(s)+f5(s);
        fkl_adj(s)=f2(s)+f3(s)+(f4(s)+f5(s))/N_bar(s);
    end
    fm=f1+sum(f2)+sum(f3)+sum(f4)+sum(f5);

    acc=sum(f4)+sum(f5);
    kl_proportions=-f1;
    kl_covs=-sum(f2);
    kl_centres=-sum(f3);

    if verbose
        disp(sprintf('Iter=%d, F1=%1.2f, F2=%1.2f, F3=%1.2f, F4=%1.2f, F5=%1.2f, Fm=%1.2f',loops,f1,sum(f2),sum(f3),sum(f4),sum(f5),fm));
    end

    % Convergence criterion
    oldlik=lik;
    lik=fm;
    if (loops>1)
        if abs((lik-oldlik)/lik) < tol
            break;
        end
    end;

    % Part-II

    for s=1:m,
        % Posterior mixing coefficients
        lambda(s)=N_bar(s)+lambda_0;
        % Posterior means
        state(s).m=(N_bar(s)*bar_mu(s,:)+beta_0*m_0)/(N_bar(s)+beta_0);
        state(s).beta=N_bar(s)+beta_0;
        % Posterior precisions
        state(s).a=N_bar(s)+a_0;
        dy=bar_mu(s,:)-m_0;
        state(s).B=N_bar(s)*state(s).bar_sigma + (N_bar(s)*beta_0*dy'*dy)/(N_bar(s)+beta_0)+B_0;
    end
end

% Put variables into data structure
mix.m=m;
mix.fm=fm;
mix.acc=acc;
mix.kl_proportions=kl_proportions;
mix.kl_covs=kl_covs;
mix.kl_centres=kl_centres;
mix.state=state;
mix.lambda=lambda;
mix.gamma=gamma;

% Put info into data structure
for j=1:m,
    mix.state(j).prior=pi_bar(j);
    mix.state(j).C=mix.state(j).B/mix.state(j).a;
end

author: spm


text:
function response_structure_description

%RESPONSE_STRUCTURE_DESCRIPTION Describe response structure fields
%   The fields of a response structure are as follows:
%       scnl: scnlobject associated with the response
%       time: time at which the response is valid
%frequencies: Nx1 vector of frequencies (Hz)
%     values: Nx1 vector of complex responses at each frequency.
%               Values should be normalized such that the maximum amplitude
%               of the response (max(abs(values))) is near 1.
%      calib: Scaler alibration value used to scale the normalized responses
%               to their real value. NOTE that many waveform objects already
%               have this value applied.
%      units: Character string describing the units
% sampleRate: Sample rate in Hz
%   respFile: Character string of name (and path) of response file.
%     status: Open-ended character string describing any issues with file
%
%  NOTE ABOUT CALIBRATION VALUES AND ANTELOPE DATABASES: In many cases the
%  calib value and units will already have been applied to a waveform if it
%  was loaded with WAVEFORM from an Antelope database. WAVEFORM attempts to
%  read the CALIB from the wfdisc table. DB_GET_RESPONSE retrieves the
%  calib value from the calibration table. Its inclusion here allows the
%  user to verify that the two values are the same. If they are not, it may
%  be necessary to reconcile the two. The units variable (translated from
%  the segtype field) should match as well.

author: geoscience-community-codes

repo: GISMO


text:
function [R,S]=rga(G)
% RGA   Relative Gain Array
%
% R=rga(G) returns the Relative Gain Array, R based on the gain matrix, G.
%
% For square matrix, G, R is the normal RGA, defined as follows:
%   R = inv(G') .* G
% If G is nonsquare, R is the so called General RGA:
%   R = pinv(G') .* G
%
% The RGA is a very useful tool for control structure design. For a square
% system (the number of inputs equals to the number of output), the RGA can
% be used to find input-output pairs (one input to control one output).
% There are many publications on how to pair input and output using the
% RGA. The basic rule is:
% 1. Avoid input and output pairs which have negative relative gains.
% 2. Avoid input and output pairs which have large relative gains.
% 3. Select input and output pairs which have the relative gain close to 1.
%
% Interestingly, if the RGA is repeatedly calculated based on R, i.e.
%
% R = rga(R);
%
% Eventualy, the RGA will converge to a selection (0,1)-matrix. With the
% second output,
%
% [R,S] = rga(G) produces a selection matrix, S based on the repeated RGA.
%
% where S is (0,1) matrix with the same dimension as G. Each row and each
% column of S have only one element of 1, but all others are 0. The element
% of 1 indicates which input (column) should be used to pair which output
% (row).
% Note: The selection matrix, S, is only a recommendation. It may not
% satisfy the rules stated above and may not be appropriate for some
% systems.
%
% For nonsquare systems, the general RGA can also be used to select inputs
% if the number of inputs is larger than the number of outputs, or to
% choose outputs, if the system has more outputs than inputs. The second
% output in this case gives the input (output) effectiveness, E.
%
% [R,E] = rga(G)
%
% Where E is a column vector corresponding to the input (output). Inputs
% (outputs) with largest input (output) effectiveness are recommented to be
% selected for the control system.
%
% By Yi Cao at Cranfield University on 7th March 2008
%
% References:
%
% 1. Bristol, E.H., On a new measure of interactions for multivariable process
% control. IEEE Trans. Automatic Control, AC-11:133-134, 1966.
% 2. Cao, Y and Rossiter, D, An input pre-screening technique for control
% structure selection, Computers and Chemical Engineering, 21(6), pp.
% 563-569, 1997.

% Examples:
%
% Example 1: A random 5 x 5 matrix
%{
G = randn(5);
[R,S] = rga(G)
%}
% Example 2: A random 5 x 10 matrix for input selection
%{
G = randn(5,10);
[R,IE] = rga(G)
%}
% Example 3: A random 10 x 5 matrix for output selection
%{
G = randn(10,5);
[R,OE] = rga(G)
%}
%

% Check input and output
error(nargchk(1,1,nargin));
error(nargoutchk(0,2,nargout));

% The RGA
R=pinv(G.').*G;
if nargout<2
    return
end

[m,n]=size(G);
% The square case
if n==m
    S=R;
    while norm(round(S)-S)>1e-9
        S(S<0)=0;
        S = inv(S.').*S;
    end
    return
end

% The nonsquare case
S = sqrt(sum(R))';
if m>n
    S = sqrt(sum(R,2));
end


text:
function newstr = untabify(str, n)
%UNTABIFY Convert tabs to spaces.
%
%   NEWSTR = UNTABIFY(STR) converts each tab character in STR to the
%   appropriate number of space characters.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 11:17:52 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 2, nargin));
   if ~ischar(str)
      error('First argument must be a string.');
   end
   if nargin < 2
      n = 8;
   else
      if any(size(n) ~= 1) | ~isnumeric(n) | n ~= round(n) | n < 0;
         error('Second argument must be a non-negative integer.');
      end
   end

   tab = sprintf('\t');                 % tab character

   strlen = length(str);                % length of input string

   newstr = '';
   if strlen == 0
      return
   end

   tabpos = find(str == tab);           % position of tabs in string
   tabpos = [ tabpos strlen+1 ];        % add "end of line mark"

   ntabs = length(tabpos);              % number of tabs

   newstr = str(1:tabpos(1)-1);
   newstrlen = length(newstr);

   for i = 1:ntabs-1
      nspc = n*ceil((newstrlen+1)/n) - newstrlen;
      newstr = [ newstr blanks(nspc) str(tabpos(i)+1:tabpos(i+1)-1) ];
      newstrlen = length(newstr);
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function X = lyap(A, B, C)
%LYAP  Solve continuous-time Lyapunov equations.
%
%   X = LYAP(A,C) solves the special form of the Lyapunov matrix
%   equation:
%
%           A*X + X*A' = -C
%
%   X = LYAP(A,B,C) solves the general form of the Lyapunov matrix
%   equation (also called Sylvester equation):
%
%           A*X + X*B = -C
%
%   See also  DLYAP.

%	S.N. Bangert 1-10-86
%	Copyright (c) 1986-98 by The MathWorks, Inc.
%	$Revision: 1.4 $  $Date: 1997/12/01 22:11:32 $
%	Last revised JNL 3-24-88, AFP 9-3-95

ni = nargin;

if ni==2,
   C = B;
   B = A';
end

[ma,na] = size(A);
[mb,nb] = size(B);
[mc,nc] = size(C);

% A and B must be square and C must have the rows of A and columns of B
if (ma ~= na) | (mb ~= nb) | (mc ~= ma) | (nc ~= mb)
   error('Dimensions do not agree.')
elseif ma==0,
   X = [];
   return
end

% Perform schur decomposition on A (and convert to complex form)
[ua,ta] = schur(A);
[ua,ta] = rsf2csf(ua,ta);
if ni==2,
   % Schur decomposition of A' can be calculated from that of A.
   j = ma:-1:1;
   ub = ua(:,j);
   tb = ta(j,j)';
else
   % Perform schur decomposition on B (and convert to complex form)
   [ub,tb] = schur(B);
   [ub,tb] = rsf2csf(ub,tb);
end

% Check all combinations of ta(i,i)+tb(j,j) for zero
p1 = diag(ta).'; % Use .' instead of ' in case A and B are not real
p1 = p1(ones(mb,1),:);
p2 = diag(tb);
p2 = p2(:,ones(ma,1));
sum = abs(p1) + abs(p2);
if any(any(sum == 0)) | any(any(abs(p1 + p2) < 1000*eps*sum))
   error('Solution does not exist or is not unique.');
end

% Transform C
ucu = -ua'*C*ub;

% Solve for first column of transformed solution
y = zeros(ma,mb);
ema = eye(ma);
y(:,1) = (ta+ema*tb(1,1))\ucu(:,1);

% Solve for remaining columns of transformed solution
for k=2:mb,
   km1 = 1:(k-1);
   y(:,k) = (ta+ema*tb(k,k))\(ucu(:,k)-y(:,km1)*tb(km1,k));
end

% Find untransformed solution
X = ua*y*ub';

% Ignore complex part if real inputs (better be small)
if isreal(A) & isreal(B) & isreal(C),
   X = real(X);
end

% Force X to be symmetric if ni==2 and C is symmetric
if (ni==2) & isequal(C,C'),
   X = (X+X')/2;
end

% end lyap


text:
function [ sindex ] = tournament(pop, index1, index2)
%   Apply tournament selection to row indexed by index1 and index2
%   Vectorization on this function may not give better speed, so beter
%   to keep as it is.

global nreal ;
global nobj ;
global ncon ;

obj_col = nreal+1:nreal+nobj;
cv_col = nreal+nobj+ncon+1;

% check the dominance
flag = check_dominance([pop(index1,obj_col), pop(index1,cv_col)], ...
                       [pop(index2,obj_col), pop(index2,cv_col)]);

cdist1 = pop(index1, end);
cdist2 = pop(index2, end);
if (flag == 1)
    sindex = index1 ;
    return ;
end
if(flag == -1)
    sindex = index2 ;
    return ;
end
if(cdist1 > cdist2)
    sindex = index1 ;
    return ;
end
if(cdist2 > cdist1)
    sindex = index2 ;
    return ;
end
if(rand(1) < 0.5)
% if(randomperc() < 0.5) % SLOW !!!
    sindex = index1 ;
    return ;
else
    sindex = index2 ;
    return ;
end

author: chudur-budur


text:
% function hewma_timeseries_plot
%
% graphic display of significant voxels in whole-brain hewma analysis
% and plotting of timeseries points you click on in the image.
%
% ..
%    tor wager
% ..
if ~exist('overlayimg') == 1, overlayimg = 'hewma_sig.img';,end
disp(['Overlay image is: ' overlayimg])

% used in button-up fcn callback
E2 = EXPT;
clear EXPT

global VOL
global f
global f2
global EXPT
EXPT = E2;

    % -------------------------------------------------------------------
    % * essential stuff for the viewing
    % -------------------------------------------------------------------

    cl = mask2clusters(overlayimg);

% view clusters
cluster_orthviews(cl,'unique');
set(gcf,'WindowButtonUpFcn','[dat,files,stats,mycov] = hewma_plot_coord_btnupfcn;')

% get coordinate mapping matrix
VOL = struct('M',cl(1).M);

% prepare figure
f1 = figure('Color','w','Name','Hewma plots');
f = f1;     % button callback uses figure f

cd ..

%return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function M1 = glassbrain_avi(fps,len,clusters)
% :Usage:
% ::
%
%    M = glassbrain_avi(fps,len,clusters)
%
% works with cluster_surf.m
% assumes a gray brain where no activation, or creates one
%
% SEE ALSO cluster_surf_movie.m, which works well too.
%
% ..
%    tor wager
% ..

dofadein = 0;
dozionimg = 0;
if isempty(clusters), doclustersurf = 0;, else, doclustersurf = 1;,end

if dozionimg
% ---------------------------------------
% Special Script
% ---------------------------------------
hh = open('Brook Zion 2.jpg');
hh = hh.Brook_Zion_2;
for i = 1:3, hh2(:,:,i) = flipud(hh(:,:,i));,end
hh = hh2;
f1 = figure('Color',[0 0 .7]); image(-size(hh,2)./2,-size(hh,1)./2,hh); hold on; axis image; axis auto
set(gca,'Color',[0 0 1],'YDir','Normal'); axis off
Ps = which('surf_single_subj_T1_gray.mat'); %'c:\tor_scripts\3DheadUtility\surf_single_subj_T1_gray.mat';
load(Ps)
 h = patch('Faces',faces,'Vertices',vertices,'FaceColor',[.5 .5 .5], ...
  'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',0.01,'SpecularExponent',200);

%camzoom(3)
camzoom(1.2)
drawnow

 O = struct('name','glassZion.avi','fps',5,'length',10,'timeres',1, 'H',h, ...
    'azOffset',0,'elOffset',0,'zoom',2,'add2movie',[],'closemovie',1);

    finalTrans = .5;
    O.timecourse{1} = [0:(finalTrans/(O.fps*O.length)):finalTrans] + .01;

    M2 = make3Davi(O);  % close this one
    O.name = 'glassZion_act.avi'; O.closemovie = 0;
    camzoom(1/3); set(h,'FaceAlpha',.01);
    M1 = make3Davi(O);  % leave this one open

    axis vis3d

% ---------------------------------------
% First part: transparent brain fade-in
% ---------------------------------------
elseif dofadein
    f1 = figure('Color','w'); h = addbrain;
    view(0,90); axis off;

    set(h,'Tag','cluster1');

    O = struct('name','glassbrainX.avi','fps',fps,'length',len(1),'timeres',1, 'H',h, ...
    'azOffset',0,'elOffset',0,'zoom',1,'add2movie',[],'closemovie',0);

    finalTrans = .5;
    O.timecourse{1} = [0:(finalTrans/(O.fps*O.length)):finalTrans];
    M1 = make3Davi(O);

else

    f1 = gcf;
    h3 = get(gca,'Children');for i = 1:length(h3),h4=get(h3(i));if isfield(h4,'FaceVertexCData'),h5 = h3(i);,end,end
    h = h5; %gco;
    tmp = get(h,'FaceVertexCData'); if size(tmp,2) ~= 3 | ~(size(tmp,1) > 5000), error('Make surface object current object'),end

    O = struct('name','glassbrainX.avi','fps',fps,'length',len(1),'timeres',1, 'H',h, ...
    'azOffset',0,'elOffset',0,'zoom',1,'add2movie',[],'closemovie',0);
    M1 = avifile('glassbrainX.avi','Quality',75,'Compression','none','Fps',O.fps);

end

% ---------------------------------------
% Second part: change cluster surface
% ---------------------------------------
if doclustersurf
    cl = cat(2,clusters.Z); refZ = [0 max(cl) min(cl(cl < 0)) 0];,

    %for j = 1:length(clusters)
    %      clusters(j).origZ = clusters(j).Z;
    %  end
    %cluster_surf(clusters,h,refZ,10,{[1 0 0]},'colorscale','heatmap');

    cluster_surf(clusters,h,10,{[1 0 0]},'colorscale','heatmap');
    figure(f1)
end

% ---------------------------------------
% Second part: cluster activations
% ---------------------------------------
axis vis3d
mov = avifile('glassbrainX_act1.avi','Quality',75,'Compression','none','Fps',O.fps);

% add frames

tmp = get(h,'FaceVertexCData');
t1 = find(~(all(tmp == .5,2)));
ind = 1;

for i = [10:-1:1 1:10 10:-1:1 1:10]

    if ind > 10, ff = get(h,'FaceAlpha'); set(h,'FaceAlpha',min(1,ff+.1));, end
    fprintf(1,'%3.0f . ',i)
    M1 = get_frame(M1,i,tmp,h,t1);
    mov = addframe(mov,gca);
    ind = ind + 1;

end

mov = close(mov);

% ---------------------------------------
% Third part: tip down
% ---------------------------------------
[az,el] = view; ind = 1;
ii = repmat([10:-1:1 1:10 10:-1:1 1:10 10:-1:1 1:10],1,10);

mov = avifile('glassbrainX_act2.avi','Quality',75,'Compression','none','Fps',O.fps);

for i = 5:5:90

    view(az-2*i,el - i)

    fprintf(1,'%3.0f . ',i)
    M1 = get_frame(M1,ii(ind),tmp,h,t1);
    mov = addframe(mov,gca);

    ind = ind + 1;
    lightFollowView
end

mov = close(mov);

% ---------------------------------------
% Fourth part: spin
% ---------------------------------------
[az,el] = view; ind = 1;
ii = repmat([10:-1:1 1:10 10:-1:1 1:10 10:-1:1 1:10],1,10);

mov = avifile('glassbrainX_act3.avi','Quality',75,'Compression','none','Fps',O.fps);

for i = 0:10:360

    if i < 70
        view(az-i,el+i./2)
    elseif i == 70
        [dummy,el] = view;
        view(az-i,el)
    else
        view(az-i,el);
    end

    fprintf(1,'%3.0f . ',i)
    M1 = get_frame(M1,ii(ind),tmp,h,t1);
    mov = addframe(mov,gca);

    ind = ind + 1;
    lightFollowView
end

mov = close(mov);

M1 = close(M1);

return

% ---------------------------------------
% Alternative: Random blobs
% ---------------------------------------
n = noisevector(size(tmp,1),[1:-.1:0],1);
gray = .5 * ones(size(tmp));
set(h,'FaceVertexCData',gray);
t2 = gray;

ind = 1;
ni = find(n > 1);
 bins = round(1:length(ni)./30:length(ni));
 cols = rand(length(bins),3);
cols(:,3) = 0;

for i = 1:length(bins)-1
    whh{ind} = ni(ni > bins(i) & ni < bins(i+1)) ;
    t2(whh{ind},:) = repmat(cols(ind,:),length(whh{ind}),1);
    ind = ind+1;
end
set(h,'FaceVertexCData',t2);

function M1 = get_frame(M1,i,tmp,h,t1)

    fprintf(1,'%3.0f . ',i)
    t2 = tmp; t2(t1,:) = 2.5 .* t2(t1,:) ./ (i .^ .7); % scale down

    t2(t1,:) = t2(t1,:) + .05 * randn(length(t2(t1,1)),3); % add a little noise
        t2(t2 > 1) = 1;
        t2(t1,3) = 0;        % no blue

    t2(sum(t2,2) < .2,:) = .5;                  % replace low values with gray
    set(h,'FaceVertexCData',t2);
    drawnow

    H = gca;
    %mov = addframe(mov,H);
    M1 = addframe(M1,H);

 return


text:
function y = cbrt(x)
%CBRT   The real cube root
%
%   CBRT(X) is the real cube root of X (assuming X is real).  X
%   can be any shape.

  y = abs(x).^(1/3);
  y(x < 0) = -y(x < 0);
end


text:
function grid_params = bp_parse_grid_params( per_pulse_metadata, num_samples_used, varargin )
%BP_PARSE_GRID_PARAMS Compute grid for backprojection
% BP_PARSE_GRID_PARAMS(PER_PULSE_METADATA, NUM_SAMPLES_USED, ...
%    'PropertyName', PropertyValue, ...)
%
%       Property name     Description
%       center            Center of image to form in ECF coordinates.
%                         Default is scene reference point.
%       grid_type         'slant', 'ground', or 'DEM'.  Default is ground.
%       image_size_meters Scene size for image formation ([range_size
%                         azimuth_size]) in meters.  Default is total size
%                         supported by collect.
%       image_size_pixels Size of the image in pixels (only used for images
%                         formed to a DEM)
%       sample_rate       Samples per IPR.  Default is 1.5.
%       grid              The grid of image points (ECEF) onto which the
%                         image should be formed.  The format of the grid
%                         parameters is identical to the (MxNx3) grid
%                         returned from this routine if no grid is
%                         supplied.
%
% Output is a structure with the following fields:
%       center            Same as input property
%       grid_type         Same as input property
%       grid              Same as input property
%       sample_rate       Same as input property
%       range_vect_hat    Unit vector point from center of aperture to
%                         image center
%       image_size_meters Extent of the image in meters
%       image_size_pixels Size of the image in pixels
%       row_coords        Vector describing position of each row.  For
%                         planar grids, this is the distance from the
%                         center of the image in meters.  For DEMs, this is
%                         the latitude position in degrees of each row.
%       col_coords        Vector describing position of each column.
%       row_unit_vector   For planar grids, the unit vector in direction of
%                         increasing row
%       col_unit_vector   For planar grids, the unit vector in direction of
%                         increasing column
%
% Authors: Wade Schwartzkopf and Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Calculate default image extent and supported resolution (in slant plane)
pulse_bandwidth = per_pulse_metadata.SCSS * (num_samples_used-1);
[max_resolution, max_extent] = pulse_info_to_resolution_extent(...
    per_pulse_metadata.TxPos - per_pulse_metadata.SRPPos, ... % Line-of-sight vector between ARP and ORP
    per_pulse_metadata.SC0 + (pulse_bandwidth/2), per_pulse_metadata.SCSS, ...
    pulse_bandwidth);

% Parse input parameters
p1 = inputParser;
p1.KeepUnmatched=true;
p1.addParamValue('grid_type',  'ground',                     @(x) any(strcmp(x,{'slant','ground','DEM'})));
p1.addParamValue('grid',       []); % Single grid that fits in memory.  Perhaps in the future, this could also allow for a file input.
p1.addParamValue('center',     mean(per_pulse_metadata.SRPPos), @(x) isempty(x) | (isvector(x)&&(length(x)==3))); % Centerpoint of formed image (ECEF, meters)
p1.addParamValue('image_size_meters', max_extent,            @(x) (isvector(x)&&(length(x)==2))); % Scene extent [x y] (m)
p1.addParamValue('image_size_pixels', [],                    @(x) (isvector(x)&&(length(x)==2))); % Scene extent [x y] (pixels)
p1.addParamValue('sample_rate',1.5,                          @(x) (isvector(x)&&(length(x)<=2)));
p1.FunctionName = mfilename;
p1.parse(varargin{:});
grid_params = p1.Results;

% Some geometry info
range_vect = grid_params.center - ... % Range vector points from center of aperture to image center
    ((per_pulse_metadata.TxPos(ceil(length(per_pulse_metadata.TxPos)/2),:) + ... % Use middle pulse
    per_pulse_metadata.RcvPos(ceil(length(per_pulse_metadata.RcvPos)/2),:))/2); % Bistatic definition, average of tx and rcv
% range_vect = grid_params.center - ... % Range vector points from center of aperture to image center
%     per_pulse_metadata.TxPos(ceil(length(per_pulse_metadata.TxPos)/2),:); % Use middle pulse
grid_params.range_vect_hat = range_vect/norm(range_vect); % Unit vector
vel_vector = per_pulse_metadata.TxPos(end,:)-per_pulse_metadata.TxPos(1,:); % Velocity vector
spn=cross(grid_params.range_vect_hat,vel_vector); % Slant plane normal
spn=spn/norm(spn); % Unit vector
% Slant plane normal to point "up" (away from origin).  This is affected by left/right looking.
spn = spn * sign(grid_params.center*cross(grid_params.range_vect_hat,vel_vector).'); % Left/right fix

% Determine final image size.  IMAGE_SIZE_PIXELS input parameter is ignored
% in all cases except for DEM.  We only allow it for the DEM case, since we
% don't currently have a robust automated way to determine it there.
if ~isempty(grid_params.grid) % Arbitrary grid was passed in.  Nothing to compute.
    grid_params.grid_type = 'ignore';
    grid_params.center = squeeze(mean(mean(grid_params.grid,1),2)); % Centroid of points in grid
    % grid_params.image_size_meters = ???
    grid_params.image_size_pixels = [size(grid_params.grid,1), size(grid_params.grid,2)];
    grid_params = rmfield(grid_params,'sample_rate'); % Not valid
    return; % Nothing left to do for explicitly passed grid
elseif strcmpi(grid_params.grid_type,'DEM')
    if isempty(grid_params.image_size_pixels) % Default
        % Assures that at least sample rate is achieved in all directions.
        grid_params.image_size_pixels = ceil(grid_params.image_size_meters.*...
            max(grid_params.sample_rate)./min(max_resolution));
    end
else % Planar image (slant or ground)
    % Adjust distances from slant to ground if necessary
    if strcmpi(grid_params.grid_type,'ground')
        gpn = wgs_84_norm( grid_params.center ); % Ground plane normal, based on (inflated) wgs_84 ellipsoid
        graze = asin(gpn.'*(-grid_params.range_vect_hat.')); % radians
        twist = asin(cross(gpn.',spn)*(-grid_params.range_vect_hat.')/...
            norm(cross(-grid_params.range_vect_hat.',gpn.'))); % radians
        if any(strcmp(p1.UsingDefaults,'image_size_meters'))
            grid_params.image_size_meters(1) = grid_params.image_size_meters(1)/cos(graze);
            grid_params.image_size_meters(2) = grid_params.image_size_meters(2)/cos(twist);
        end
        max_resolution(1) = max_resolution(1)/cos(graze);
        max_resolution(2) = max_resolution(2)/cos(twist);
    end
    grid_params.image_size_pixels = ceil(grid_params.image_size_meters.*...
        grid_params.sample_rate./max_resolution);
end

% Compute grid onto which the image will be formed
switch grid_params.grid_type
    % For planar images (ground and slant) compute the unit vectors in the
    % direction of increasing row/column direction
    case 'ground'
        row_vector = grid_params.range_vect_hat-(grid_params.range_vect_hat*gpn)*gpn.'; % Project range vector to ground
        grid_params.row_unit_vector = row_vector/norm(row_vector); % Unit vector
        col_vector = -cross(grid_params.row_unit_vector,gpn); % Vector in direction of increasing column (ECEF)
        grid_params.col_unit_vector = col_vector/norm(col_vector); % Unit vector
    case 'slant'
        % Slant plane really has no meaning in backprojection.  Any
        % arbitrary grid can be defined.  We define it here merely for
        % comparison against PFA.
        grid_params.row_unit_vector = grid_params.range_vect_hat; % Unit vector in direction of increasing row (ECEF)
        grid_params.col_unit_vector = -cross(grid_params.row_unit_vector,spn); % Unit vector in direction of increasing column (ECEF)
    case 'DEM'
        % We'll also get the DEM data for the image.  We assume that the
        % entire image fits within .2 degrees of lat and long and that the
        % entire DEM can fit into memory (reasonable for SRTM, but probably
        % not LIDAR).
        center_lla  = ecf_to_geodetic(grid_params.center);
        [lats,lons,elevations] = get_DEM_heights_region([center_lla(1)-0.1,center_lla(2)-0.1], ...
            [center_lla(1)+0.1,center_lla(2)+0.1]);
        [lats, lons] = meshgrid(lats, lons);
        % Return DEM interpolation function, rather than the grid itself,
        % since for image formation we will need values not lying exactly
        % on the DEM points.
        grid_params.F_1 = TriScatteredInterp(lats(:), lons(:), elevations(:));

        % Convert scene extent from meters to arc-degrees
        e2 = 6.6943799901377997e-3;  % eccentricity squared of Earth (WGS 84 value)
        a = 6378137.0;               % semimajor radius of the Earth (WGS 84 value)
        b = a .* sqrt( 1.0 - e2);    % semiminor radius of the Earth
        lat = pi/180*center_lla(1);
        M = (a*b)^2/((a*cos(lat))^2+(b*sin(lat))^2)^(3/2);
        N = a^2/sqrt((a*cos(lat))^2+(b*sin(lat))^2);
        arcradius_lat = pi/180*M;          % length of arc-degree in latitude direction (m)
        arcradius_lon = pi/180*cos(lat)*N; % length of arc-degree in longitude direction (m)
        scene_extent_lat = grid_params.image_size_meters(1)/arcradius_lat;
        scene_extent_lon = grid_params.image_size_meters(2)/arcradius_lon;
end
% Compute row/column coordinates
switch grid_params.grid_type
    % For a plane, coordinates are distances (in meters) from the center
    % of image for each row/column
    case {'ground','slant'}
        grid_params.row_coords = grid_params.image_size_meters(1)*...
            linspace(-1,1,grid_params.image_size_pixels(1))/2;
        grid_params.col_coords = grid_params.image_size_meters(2)*...
            linspace(-1,1,grid_params.image_size_pixels(2))/2;
    % For a DEM, coordinates of each row/column are in lat/lon degrees;
    % that is, angularly, rather than along a plane.
    case 'DEM'
        grid_params.row_coords = center_lla(1) + (scene_extent_lat*...
            linspace(-1,1,grid_params.image_size_pixels(1))/2);
        grid_params.col_coords = center_lla(2) + (scene_extent_lon*...
            linspace(-1,1,grid_params.image_size_pixels(2))/2);
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
rng default;
% Signal space
N = 256;
% Number of measurements
M = 64;
% Sparsity level
K = 10;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K);
% Generate bi-uniform signals
x = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
y = Phi.apply(x);
% OMP MMV solver instance
result = spx.pursuit.single.omp_chol(double(Phi), y, K, 1e-6);
% Solution vector
z = result.z;

stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
spx.commons.sparse.print_recovery_performance(stats);

mf = spx.graphics.Figures();
mf.new_figure('OMP solution');
subplot(411);
stem(x, '.');
title('Sparse vector');
subplot(412);
stem(z, '.');
title('Recovered sparse vector');
subplot(413);
stem(abs(x - z), '.');
title('Recovery error');
subplot(414);
stem(y, '.');
title('Measurement vector');

author: indigits

repo: sparse-plex


text:
function nMM=grMaxMatch(E)
% Function nMM=grMaxMath(E) solve the maximal matching problem.
% Input parameter:
%   E(m,2) or (m,3) - the edges of graph and their weight;
%     1st and 2nd elements of each row is numbers of vertexes;
%     3rd elements of each row is weight of edge;
%     m - number of edges.
%     If we set the array E(m,2), then all weights is 1.
% Output parameter:
%   nMM - the list of the numbers of edges included
%     in the maximal (weighted) matching.
% Uses the reduction to integer LP-problem.
% Required the Optimization Toolbox v.3.0.1 or over.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

% ============= Input data validation ==================
if nargin<1,
  error('There are no input data!')
end
[m,n,E] = grValidation(E); % E data validation

% ============= Parameters of integer LP problem ==========
A=zeros(n,m); % for incidence matrix
A(E(:,1:2)+repmat(([1:m]'-1)*n,1,2))=1; % we fill the incidence matrix
options=optimset('bintprog'); % the default options
options.Display='off'; % we change the output

% ============= We solve the integer LP problem ==========
xmin=bintprog(-E(:,3),A,ones(n,1),[],[],[],options);
nMM=find(round(xmin)); % the answer - numbers of edges
return

author: HuangCongQing


text:
function pass = test_nonlinSys1_C1(pref)
% Test 2x2 system (sin/cos). This is nonlinearification of the test
%       test_linearSystem1
%
% Asgeir Birkisson, April 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

% Smooth domain:
d = [-pi pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% Colloc1
pref.discretization = @chebcolloc1;

A = chebop(@(x,u,v) [u - diff(v,2) + u.^2; diff(u) + sin(v)],d);
A.lbc = @(u,v) u-1;
A.rbc = @(u,v) [v-1/2; diff(v)];

u12 = mldivide(A, f, pref);
u1 = u12{1}; u2 = u12{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = chebfun(A.rbc(u1,u2));

pass(1) = norm( chebfun(A(x, u1, u2))) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
filename='Cantilever_triangle_coarse';
  ptype = '';
method = '';
materialType = '';
initial_case = '';
cost = {''};
weights = '';
constraint = {''};
optimizer = '';
incrementFactor = 1;
designVariable = '';
filterType = '';

nsteps = 1;
Vfrac_final = '';
Perimeter_target = 1;
optimality_final = '';
constr_final = '';

Vfrac_initial = '';
optimality_initial = '';
constr_initial = '';

TOL.rho_plus = '';
TOL.rho_minus = '';
TOL.E_plus = '';
TOL.E_minus = '';
TOL.nu_plus = '';
TOL.nu_minus = '';

plotting = '';
showBC = '';
monitoring = '';
monitoring_interval = 5;

isOld = false;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function weight = randInitializeWeights(inputs, outputs)
    % Initialize the weights randomly.
    % This will break the symmetry while training the neural network.
    epsilon = 0.12;
    weight = rand(outputs, inputs + 1) * (2 * epsilon) - epsilon;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function FemMat = fem_remove_elem(FemMat, iRemoveElem)
% FEM_REMOVE_ELEM: Remove some elements from a 3D mesh

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2020

% Remove elements
FemMat.Elements(iRemoveElem, :) = [];
FemMat.Tissue(iRemoveElem) = [];
if isfield(FemMat, 'Tensors') && ~isempty(FemMat.Tensors)
    FemMat.Tensors(iRemoveElem, :) = [];
end

% Find vertices to remove
nVert = size(FemMat.Vertices, 1);
iVertCut = setdiff(1:nVert, unique(FemMat.Elements(:)));
% Re-numbering matrix
iVertKept = setdiff(1:nVert, iVertCut);
iVertMap = zeros(1, nVert);
iVertMap(iVertKept) = 1:length(iVertKept);
% Remove vertices
FemMat.Vertices(iVertCut,:) = [];
% Renumber vertices in elements list
FemMat.Elements = iVertMap(FemMat.Elements);

author: brainstorm-tools


text:
function status = test_wfs_iir_prefilter(modus)
%TEST_WFS_IIR_PREFILTER tests the IIR WFS pre-equalization filter
%
%   Usage: status = test_wfs_iir_prefilter(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual (not available)
%
%   Output parameters:
%       status - true or false
%
%   TEST_WFS_IIR_PREFILTER(modus) test the WFS pre-euqalization IIR filter
%   design. This works only in Matlab as the Signal Processing Toolbox is used.
%   See wfs_iir_prefilter.m for details

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

status = false;

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Configuration ===================================================
conf = SFS_config;

%% ===== Calculation =====================================================
% call with default values
hpre1 = wfs_iir_prefilter(conf)
conf.fs = 44100;
conf.hpreflow = 200;
conf.hprefhigh = 1500;
conf.hpreBandwidth_in_Oct = 2;
conf.hpreIIRorder = 4;
hpre2 = wfs_iir_prefilter(conf)

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
function [ Y ] = TMED( X, h )
% TMED  - The fuzzy filter with Triangular function with MEDian center
% x     - input (noisy image)
% y     - output (de-noised image)
% F     - window function, weighting function

if nargin < 2
    h=3;    % local search window size: h=3; h=5; h=7; h=9; h=11;
    f=1;    % padding value for window: f=1; f=2; f=3; f=4; f=5;
else
    f=(h-1)/2;    % padding value for window: f=1; f=2; f=3; f=4; f=5;
end
F = zeros(h*h,1);
Y  = zeros(size(X));

X = double(padarray(X,[f f],'symmetric'));
[m n ~] = size(X);

for i=1+f:1:m-f
   for j=1+f:1:n-f
       x=reshape(X(i-f:i+f, j-f:j+f),[],1);
       xmax = max(x);
       xmin = min(x);
       xmed = median(x);
       xmm = max(xmax-xmed, xmed-xmin);

       F(:,:) = 0;
       if xmm==0
           F(:,:) = 1;
       else
           F = (abs(x-xmed)/xmm);
       end

       Y(i-f,j-f) = sum(sum(F.*x))/sum(sum(F));
       clear xmax xmin xmed xmm;

   end
end
end


text:
function obj = addConstraint(obj, label, nodes, cstr_array)
    % Add NLP constraints constraints
    %
    %
    % Parameters:
    % label: the label name (column) of the constraints @type char
    % nodes: the node list of the variable @type rowvec
    % cstr_array: an array of structures that could be used to construct
    % the NlpFunction object. @type struct
    %
    % @see NlpFunction, removeConstraint
    %
    % @note

    cstr_names = obj.ConstrTable.Properties.VariableNames;

%     if ismember(label,cstr_names)
%         warning('The NLP constraint (%s) already exist.\n Overwriting the existing constraint.', label);
%     end

    if ischar(nodes)
        switch nodes
            case 'first'
                node_list = 1;
            case 'last'
                node_list = obj.NumNode;
            case 'terminal'
                node_list = [1 obj.NumNode];
            case 'except-first'
                node_list = 2:obj.NumNode;
            case 'except-last'
                node_list = 1:obj.NumNode-1;
            case 'except-terminal'
                node_list = 2:obj.NumNode-1;
            case 'all'
                node_list = 1:obj.NumNode;
            case 'cardinal'
                node_list = 1:2:obj.NumNode;
            case 'interior'
                node_list = 2:2:obj.NumNode-1;
            otherwise
                error('Unknown node type.');
        end
    else
        if ~isnumeric(nodes)
            error(['The node must be specified as a list or following supported characters:\n',...
                '%s'],implode({'first','last','all','except-first','except-last','except-terminal', 'cardinal', 'interior'},','));
        else
            node_list = nodes;
        end
    end

    assert(length(node_list)==length(cstr_array),...
        'The length of the array (%d) of constraint structures must be the same as the number of nodes (%d) to be imposed.',...
        length(cstr_array),length(node_list));

    % create empty NlpVariable array
    if isa(cstr_array,'NlpFunction')
        constr = repmat(NlpFunction(),obj.NumNode,1);
        for j=1:numel(node_list)
            constr(node_list(j)) = cstr_array(j);
        end
    elseif isstruct(cstr_array)

        constr = repmat(NlpFunction(),obj.NumNode,1);
        for j=1:numel(node_list)
            constr(node_list(j)) = NlpFunction(cstr_array(j));
        end
    end

    % add to the constraints table
    if ismember(label,cstr_names)
        obj.ConstrTable.(label)(node_list) = constr(node_list);
    else
        obj.ConstrTable.(label) = constr;
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function MDP = DEMO_MDP_questions
% Demo of active inference for visual salience
%__________________________________________________________________________
%
% This routine provide simulations of reading to demonstrate deep temporal
% generative models. It builds upon the scene construction simulations to
% equip the generative model with a second hierarchical level. In effect,
% this creates an agent that can accumulate evidence at the second level
% based upon epistemic foraging at the first. In brief, the agent has to
% categorise a sentence or narrative into one of two categories (happy or
% sad), where it entertains six possible sentences. Each sentence comprises
% four words, which are themselves constituted by two pictures or graphemes
% These are the same visual outcomes used in previous illustrations of
% scene construction and saccadic searches.
%
% Here, the agent has policies at two levels. The second level policy (with
% just one step into the future) allows it to either look at the next word
% or stay on the current page and make a decision. Concurrently, a first
% level policy entails one of four saccadic eye movements to each quadrant
% of the current page, where it will sample a particular grapheme.
%
% This provides a rough simulation of reading - that can be made more
% realistic by terminating first level active inference, when there can be
% no further increase in expected free energy (i.e., all uncertainty about
% the current word has been resolved). The subsequent inferred hidden
% states then become the outcome for the level above.
%
% To illustrate the schemes biological plausibility, one can change the
% agent's prior beliefs and repeat the reading sequence under violations of
% either local (whether the graphemes are flipped vertically) or globally
% (whether the sentence is surprising) expectations. This produces a
% mismatch negativity (MMN) under local violations) and a MMN with a
% P300 with global violations.
%
% see also: DEM_demo_MDP_habits.m and spm_MPD_VB_X.m
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEMO_MDP_questions.m 7766 2020-01-05 21:37:39Z karl $

% set up and preliminaries: first level
%==========================================================================
%
%--------------------------------------------------------------------------

rng('default')

% first level (lexical)
%==========================================================================
% question 1: {'noun (7)'}
% question 2: {'noun (7)','adverb (9)'}
% question 3: {'adjective (8)','noun (7)','adverb (9)'}
%--------------------------------------------------------------------------
% probabilistic mapping from hidden states to outcomes: A
%--------------------------------------------------------------------------
sentence{1} = {'Is there a ','&1_1','?'};
sentence{2} = {'Is any ','&2_1','&2_3','?'};
sentence{3} = {'Is a ','&3_2','&3_1','&3_3','?'};
sentence{4} = {'Yes','!'};
sentence{5} = {'No','!'};
sentence{6} = {'Not sure','!'};
sentence{7} = {'Ready','?'};

syn{1}      = {'Ready','OK','please begin'};
syn{2}      = {'Not sure','I dont know','Cannot say'};
syn{3}      = {'No','Sorry, no','There is not'};
syn{4}      = {'Yes','Yes there is','Well done'};

% assemble hidden states and transition probabilities
%--------------------------------------------------------------------------
name  = {};
for i = 1:numel(sentence)
    for j = 1:numel(sentence{i})
        phrase = sentence{i}(j);
        if ~any(ismember(name,phrase))
        name(end + 1) = phrase;
        end
    end
end

label.factor{1} = 'noun';   label.name{1}  = {'square ','triangle '};
label.factor{2} = 'adj.';   label.name{2}  = {'green ','red '};
label.factor{3} = 'adverb'; label.name{3}  = {'above','below'};
label.factor{4} = 'syntax'; label.name{4}  = name;

% prior beliefs about initial states D
%--------------------------------------------------------------------------
for i = 1:numel(label.factor)
    n    = numel(label.name{i});
    D{i} = ones(n,1)/n;
end

% restricted initial states to the beginning of a sentence
%--------------------------------------------------------------------------
state = label.name{4};
D{4}  = spm_zeros(D{4});
for i = 1:numel(sentence)
    j = ismember(state,sentence{i}(1));
    D{4}(j) = 1;
end

% probabilistic mapping from hidden states to outcomes: A
%--------------------------------------------------------------------------
Nf    = numel(D);
for f = 1:Nf
    Ns(f) = numel(D{f});
end

outcome = {};
for i = 1:numel(label.name)
    outcome = [outcome,label.name{i}];
end
for i = 1:numel(syn)
    outcome = [outcome,syn{i}];
end
outcome = unique(outcome);
j     = [];
for i = 1:numel(outcome)
    if outcome{i}(1) ~= '&'
        j = [j,i];
    end
end
outcome = outcome(j);

% single outcome modality with multiple phrases
%--------------------------------------------------------------------------
label.modality{1} = 'phrase';   label.outcome{1} = outcome;

for f1 = 1:Ns(1)
    for f2 = 1:Ns(2)
        for f3 = 1:Ns(3)
            for f4 = 1:Ns(4)

                % indices of the state
                %----------------------------------------------------------
                j = {f1,f2,f3,f4};

                % assemble phrases under this state
                %----------------------------------------------------------
                name = label.name{4}{f4};
                if name(1) == '&'
                    p   = eval(name(end));
                    q   = eval(['f' name(end)]);
                    out = label.name{p}(q);
                else
                    out = name;
                    for i = 1:numel(syn)
                        if any(ismember(syn{i},name))
                            out = syn{i};
                        end
                    end
                end

                % placing likelihood matrix
                %----------------------------------------------------------
                i = ismember(outcome,out);
                A{1}(i,j{:}) = 1/sum(i);

            end
        end
    end
end

% transitions: B{f} for each factor
%--------------------------------------------------------------------------
for f = 1:Nf
    B{f} = eye(Ns(f));
end

% specifies syntax
%--------------------------------------------------------------------------
B{4}  = spm_zeros(B{4});
for s = 1:numel(sentence)
    for t = 2:numel(sentence{s})
        i = find(ismember(state,sentence{s}(t - 1)));
        j = find(ismember(state,sentence{s}(t)));
        B{4}(j,i) = 1;
    end
    B{4}(j,j) = 1;
end

% MDP Structure
%--------------------------------------------------------------------------
mdp.T = 5;                      % number of updates
mdp.A = A;                      % observation model
mdp.B = B;                      % transition probabilities
mdp.D = D;                      % prior over initial states

mdp.label = label;
mdp.chi   = 1/32;

clear A B D

MDP   = spm_MDP_check(mdp);

% % check belief updates (and behaviour)
% %------------------------------------------------------------------------
% MDP.s = [1 2 1 8]';
% MDP   = spm_MDP_VB_X(MDP);
%
% spm_figure('GetWin','Figure 1'); clf
% spm_MDP_VB_trial(MDP,[3 4],1);
%
% % illustrate phase-precession and responses
% %------------------------------------------------------------------------
% spm_figure('GetWin','Figure 2'); clf
% spm_MDP_VB_LFP(MDP,[],4);
%
% cell2mat(MDP.label.outcome{1}(MDP.o))

% set up and preliminaries: first level
%==========================================================================
%
%--------------------------------------------------------------------------

%% second level (narrative)
%==========================================================================
% question 1: {'noun (7)'}
% question 2: {'noun (7)','adverb (9)'}
% question 3: {'adjective (8)','noun (7)','adverb (9)'}
%--------------------------------------------------------------------------
label.factor{1}  = 'narrative';    label.name{1}  = {'ready','question','answer'};
label.factor{2}  = 'question';     label.name{2}  = {'is?','where?','what?'};
label.factor{3}  = 'upper colour'; label.name{3}  = {'green','red'};
label.factor{4}  = 'lower colour'; label.name{4}  = {'green','red'};
label.factor{5}  = 'upper shape';  label.name{5}  = {'square','triangle'};
label.factor{6}  = 'lower shape';  label.name{6}  = {'square','triangle'};
label.factor{7}  = 'noun';         label.name{7}  = {'square','triangle'};
label.factor{8}  = 'adjective';    label.name{8}  = {'green','red'};
label.factor{9}  = 'adverb';       label.name{9}  = {'above','below'};

% prior beliefs about initial states D
%--------------------------------------------------------------------------
for i = 1:numel(label.factor)
    n    = numel(label.name{i});
    D{i} = ones(n,1)/n;
end

% known initial states
%--------------------------------------------------------------------------
D{1}(1)  = 128;

% probabilistic mapping from hidden states to outcomes: A
%--------------------------------------------------------------------------
label.modality{1} = 'noun';   label.outcome{1}  = MDP.label.name{1};
label.modality{2} = 'adj.';   label.outcome{2}  = MDP.label.name{2};
label.modality{3} = 'adverb'; label.outcome{3}  = MDP.label.name{3};
label.modality{4} = 'syntax'; label.outcome{4}  = MDP.label.name{4}(find(MDP.D{4}));

Nf    = numel(D);
for f = 1:Nf
    Ns(f) = numel(D{f});
end
for f1 = 1:Ns(1)
    for f2 = 1:Ns(2)
        for f3 = 1:Ns(3)
            for f4 = 1:Ns(4)
                for f5 = 1:Ns(5)
                    for f6 = 1:Ns(6)
                        for f7 = 1:Ns(7)
                            for f8 = 1:Ns(8)
                                for f9 = 1:Ns(9)

                                    % indices
                                    %--------------------------------------
                                    j = {f1,f2,f3,f4,f5,f6,f7,f8,f9};

                                    % answer: depending on question and beliefs
                                    %--------------------------------------
                                    Y  = 0;
                                    if f1 == 3
                                        if f2 == 1
                                            Y = (f7 == f5) | (f7 == f6);
                                        elseif f2 == 2
                                            if f9 == 1
                                                Y = (f7 == f5);
                                            elseif f9 == 2
                                                Y = (f7 == f6);
                                            end
                                        elseif f2 == 3
                                            if f9 == 1
                                                Y = (f7 == f5) & (f8 == f3);
                                            elseif f9 == 2
                                                Y = (f7 == f6) & (f8 == f4);
                                            end
                                        end
                                    end

                                    % A{1} noun:
                                    %======================================
                                    A{1}(f7,j{:}) = 1;

                                    % A{2} adjective:
                                    %======================================
                                    A{2}(f8,j{:}) = 1;

                                    % A{3} adverb:
                                    %======================================
                                    A{3}(f9,j{:}) = 1;

                                    % A{4} syntax: {'1','2','3','Y','N','?'}
                                    %======================================
                                    if f1 == 1
                                        A{4}(7,j{:})     = 1;
                                    elseif f1 == 2
                                        A{4}(f2,j{:})    = 1;
                                    elseif f1 == 3
                                        if Y
                                            A{4}(4,j{:}) = .98;
                                        else
                                            A{4}(5,j{:}) = .98;
                                        end
                                         A{4}(6,j{:})    = .02;
                                    end

                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

% likelihood mappings
%--------------------------------------------------------------------------
Ng    = numel(A);
for g = 1:Ng
    No(g) = size(A{g},1);
    if any(spm_vec(sum(A{g})) ~= 1), disp('check A'), return, end
end

% controlled transitions: B{f} for each factor
%--------------------------------------------------------------------------
for f = 1:Nf
    B{f} = eye(Ns(f));
    label.action{f} = {'stay'};
end

% transitions B(1): {'ready','question','answer'}
%--------------------------------------------------------------------------
B{1}(:,:,1) = spm_speye(Ns(1),Ns(1),-1); B{1}(1,Ns(1),1) = 1;

% control states B(2): question {'1','2' or '3'} & control states B(7):B{9}
%--------------------------------------------------------------------------
% D{2} ;      % question;     {'1','2','3'};
% D{7} ;      % noun:         {'square','triangle'}
% D{8} ;      % adjective:    {'green','red'}
% D{9} ;      % adverb:       {'above','below'}
%--------------------------------------------------------------------------
for f = [2 7 8 9]
    for k = 1:Ns(f)
        B{f}(:,:,k + 1) = 0;
        B{f}(k,:,k + 1) = 1;
    end
    label.action{f} = ['stay', label.name{f}];
end

% allowable policies (time x polcy x factor)
%--------------------------------------------------------------------------
% question 1: {'noun (7)'}
% question 2: {'noun (7)','adverb (9)'}
% question 3: {'adjective (8)','noun (7)','adverb (9)'}
%--------------------------------------------------------------------------
V         = ones(2,14,Nf);
V(1,:,2)  = 1 + [1 1 2 2 2 2 3 3 3 3 3 3 3 3];
V(1,:,7)  = 1 + [1 2 1 1 2 2 1 1 1 1 2 2 2 2];
V(1,:,8)  = 1 + [1 1 1 1 1 1 1 1 2 2 1 1 2 2];
V(1,:,9)  = 1 + [1 1 1 2 1 2 1 2 1 2 1 2 1 2];
V(2,:,:)  = 1;

% priors: (utility) C: A{4} syntax: {'1','2','3','Y','N','S'}
%--------------------------------------------------------------------------
for g = 1:Ng
    C{g}  = zeros(No(g),1);
end

% the agent expects affirmative answers
%--------------------------------------------------------------------------
C{4}(4,:) =  1/4;
C{4}(5,:) = -1/4;

% actual state of the world
%--------------------------------------------------------------------------
s    = ones(Nf,1);
s(4) = 2;

% MDP Structure
%--------------------------------------------------------------------------
mdp.MDP    = MDP;
mdp.label  = label;             % names of factors and outcomes
mdp.tau    = 4;                 % time constant of belief updating
mdp.erp    = 4;                 % initialization
mdp.chi    = 0;                 % initialization

mdp.V = V;                      % allowable policies
mdp.A = A;                      % observation model
mdp.B = B;                      % transition probabilities
mdp.C = C;                      % preferred outcomes
mdp.D = D;                      % prior over initial states (context)
mdp.s = s;                      % initial state
mdp.o = [];                     % outcomes

mdp.link = spm_MDP_link(mdp);
MDP      = spm_MDP_check(mdp);
%   spm_MDP_factor_graph(mdp);

%% illustrate questioning
%==========================================================================
clear MDP
OPTIONS.D    = 1;
[MDP(1,1:6)] = deal(mdp);

MDP   = spm_MDP_VB_X(MDP,OPTIONS);

% show belief updates (and behaviour)
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1'); clf
spm_MDP_VB_trial(MDP(1),[1 2 4],[1 3 4]);

% illustrate phase-precession and responses
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf
spm_MDP_VB_LFP(MDP,[],3);

spm_figure('GetWin','Figure 2A'); clf
spm_MDP_VB_LFP(MDP,[],3,1);

spm_figure('GetWin','Figure 3'); clf
spm_MDP_VB_ERP(MDP(4:6),[3,2]);

spm_figure('GetWin','Figure 4'); clf
for i = 1:size(MDP,2)
    subplot(4,3,i)
    spm_questions_plot(MDP(1,i))
end

% illustrate violations
%==========================================================================
j      = 5;                               % which answer
NDP    = MDP(j);                          % get states and outcomes
if NDP.o(4,3) == 4                        % switch the answer
    NDP.o(4,3) = 5;
else
    NDP.o(4,3) = 4;
end
NDP = rmfield(NDP,'link');                % remove link (outomes are given)
NDP = spm_MDP_VB_X(NDP);

% find greatest effect on belief updating (about the scene)
%--------------------------------------------------------------------------
for f = 3:6
    v(f) =  norm(spm_vec(MDP(j).X{f}) - spm_vec(NDP.X{f}),'inf');
end
[v,f] = max(v);

% responses to appropriate and inappropriate answers
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 6 expected' ); clf, spm_MDP_VB_LFP(MDP(5),[],f);
spm_figure('GetWin','Figure 7 violation'); clf, spm_MDP_VB_LFP(NDP   ,[],f);

% illustrate 'communication'
%==========================================================================

% increase efficiency (i.e., suppress neurophysiological correlates)
%--------------------------------------------------------------------------
mdp.tau    = 3;                         % time constant of belief updating
mdp.erp    = 1;                         % initialization
mdp.chi    = 1/64;                      % initialization

% create two agents
%--------------------------------------------------------------------------
clear MDP
[MDP(1:2,1:6)] = deal(mdp);

% give a second subject veridical beliefs about the scene
%--------------------------------------------------------------------------
for i = [3 4 5 6]
    MDP(2,1).D{i} = sparse(s(i),1,1,Ns(i),1);
end

% confident player answers then asks
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 8'); clf
for i = 1:length(MDP)

    if i < 5
        % first model asks and the second answers
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[2 1 2];
        MDP(2,i).n = ones(Ng,1)*[2 1 2];
    else
        % switch roles
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[1 2 1];
        MDP(2,i).n = ones(Ng,1)*[1 2 1];
    end
end

TDP   = spm_MDP_VB_X(MDP,OPTIONS);
for i = 1:size(TDP,2)
    subplot(4,3,i), spm_questions_plot(TDP(:,i))
end

% confident player asks then answers
%--------------------------------------------------------------------------
for i = 1:length(MDP)

    if i < 5
        % first model asks and the second answers
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[1 2 1];
        MDP(2,i).n = ones(Ng,1)*[1 2 1];
    else
        % switch roles
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[2 1 2];
        MDP(2,i).n = ones(Ng,1)*[2 1 2];
    end
end

TDP   = spm_MDP_VB_X(MDP,OPTIONS);
for i = 1:size(TDP,2)
    subplot(4,3,i + 6)
    spm_questions_plot(TDP(:,i))
end

% confident player tells then asks
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 9'); clf
for i = 1:length(MDP)

    if i < 7
        % first model asks and the second answers
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[1 2 2];
        MDP(2,i).n = ones(Ng,1)*[1 2 2];
    else
        % switch roles
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[1 2 1];
        MDP(2,i).n = ones(Ng,1)*[1 2 1];
    end
end

TDP   = spm_MDP_VB_X(MDP,OPTIONS);
for i = 1:size(TDP,2)
    subplot(4,3,i)
    spm_questions_plot(TDP(:,i))
end

% building fantasies
%==========================================================================
spm_figure('GetWin','Figure 10'); clf

% make both players uncertain about the scene at hand
%--------------------------------------------------------------------------
MDP(2,1).D = MDP(1,1).D;

for i = 1:length(MDP)

    if i < 3
        % first model asks and the second answers
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[1 2 1];
        MDP(2,i).n = ones(Ng,1)*[1 2 1];
    else
        % switch roles
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[2 1 2];
        MDP(2,i).n = ones(Ng,1)*[2 1 2];
    end
end

TDP   = spm_MDP_VB_X(MDP,OPTIONS);
for i = 1:size(TDP,2)
    subplot(4,3,i)
    spm_questions_plot(TDP(:,i))
end

% now repeat when precluding uncertain responses
%--------------------------------------------------------------------------
for i = 1:length(MDP)

    MDP(1,i).A{4} = MDP(1,i).A{4} > 1/2;
    MDP(2,i).A{4} = MDP(2,i).A{4} > 1/2;

    if i < 3
        % first model asks and the second answers
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[1 2 1];
        MDP(2,i).n = ones(Ng,1)*[1 2 1];
    else
        % switch roles
        %------------------------------------------------------------------
        MDP(1,i).n = ones(Ng,1)*[2 1 2];
        MDP(2,i).n = ones(Ng,1)*[2 1 2];
    end
end

TDP   = spm_MDP_VB_X(MDP,OPTIONS);
for i = 1:size(TDP,2)
    subplot(4,3,6 + i)
    spm_questions_plot(TDP(:,i))
end

return

function spm_questions_plot(MDP)
%% illustrate beliefs
%--------------------------------------------------------------------------

% probabilistic mapping from hidden states to outcomes: A
%--------------------------------------------------------------------------
% label.modality{1} = 'noun';   label.outcome{1}  = {'square','triangle'};
% label.modality{2} = 'adj.';   label.outcome{2}  = {'green','red'};
% label.modality{3} = 'adverb'; label.outcome{3}  = {'above','below'};
% label.modality{4} = 'syntax'; label.outcome{4}  = {'1','2','3','Y','N','Sil.'};

% plot question, answer and posterior beliefs
%==========================================================================
cla;

for m = 1:numel(MDP)

    % Assemble question-and-answer
    %----------------------------------------------------------------------
    question = MDP(m).o(4,2);
    answer   = MDP(m).o(4,3);

    try
        qstr = cell2mat(MDP(m).MDP.label.outcome{1}(MDP(m).mdp(2).o));
        astr = cell2mat(MDP(m).MDP.label.outcome{1}(MDP(m).mdp(3).o));
    catch
        noun     = {'square','triangle'};
        adj      = {'green','red'};
        adverb   = {'above','below'};
        if question == 1
            qstr = ['is there a ' noun{MDP(m).o(1,2)} ' ?'];
        elseif question == 2
            qstr = ['Is any ' noun{MDP(m).o(1,2)} ' ' adverb{MDP(m).o(3,2)} ' ?'];
        elseif question == 3
            qstr = ['Is a ' adj{MDP(m).o(2,2)} ' ' noun{MDP(m).o(1,2)} ' ' adverb{MDP(m).o(3,2)} ' ?'];
        else
            qstr = '!';
        end
        if answer == 4
            astr = 'Yes !';
        elseif answer == 5
            astr = 'No !';
        else
            astr = 'I''m not sure';
        end
    end

    % is the answer right (for a single player)?
    %----------------------------------------------------------------------
    ind      = num2cell(MDP(m).s(:,3));
    if answer == find(MDP(m).A{4}(:,ind{:}),1) || m > 1
        cor = spm_softmax(2*[0;1;0]);
    else
        cor = spm_softmax(2*[1;0;0]);
    end

    %  plot posterior beliefs
    %----------------------------------------------------------------------
    % label.factor{3}  = 'upper colour'; label.name{3}  = {'green','red'};
    % label.factor{4}  = 'lower colour'; label.name{4}  = {'green','red'};
    % label.factor{5}  = 'upper shape';  label.name{5}  = {'square','triangle'};
    % label.factor{6}  = 'lower shape';  label.name{6}  = {'square','triangle'};

    % upper and lower object
    %----------------------------------------------------------------------
    T       = 2;
    col{1}  = [MDP(m).X{3}(2,T) MDP(m).X{3}(1,T) 0];
    col{1}  = col{1}*MDP(m).X{5}(1,T) + (1 - MDP(m).X{5}(1,T));
    col{2}  = [MDP(m).X{3}(2,T) MDP(m).X{3}(1,T) 0];
    col{2}  = col{2}*MDP(m).X{5}(2,T) + (1 - MDP(m).X{5}(2,T));
    col{3}  = [MDP(m).X{4}(2,T) MDP(m).X{4}(1,T) 0];
    col{3}  = col{3}*MDP(m).X{6}(1,T) + (1 - MDP(m).X{6}(1,T));
    col{4}  = [MDP(m).X{4}(2,T) MDP(m).X{4}(1,T) 0];
    col{4}  = col{4}*MDP(m).X{6}(2,T) + (1 - MDP(m).X{6}(2,T));

    plot(m,0,'^','MarkerSize',24,'LineWidth',4,'Color',col{4}), hold on
    plot(m,0,'s','MarkerSize',24,'LineWidth',4,'Color',col{3})
    plot(m,1,'^','MarkerSize',24,'LineWidth',4,'Color',col{2})
    plot(m,1,'s','MarkerSize',24,'LineWidth',4,'Color',col{1})

end

try, nq = MDP(1).n(4,2); catch, nq = 1; end
try, na = MDP(1).n(4,3); catch, na = 1; end
if ~nq, nq = 1; end
if ~na, na = 1; end
text(nq, 2,qstr,'HorizontalAlignment','Center')
text(na,-1,astr,'HorizontalAlignment','Center','FontWeight','bold','Color',cor)
axis([0 (m + 1) -1.5 2.5]), axis off, axis square

% return if multiple agents
%--------------------------------------------------------------------------
if numel(MDP) > 1, return, end

% upper and lower object
%--------------------------------------------------------------------------
rgb     = {[0 1 0],[1 0 0]};
shape   = {'s','^'};
plot(1 + 1/2,1,shape{MDP.s(5)},'MarkerSize',8,'LineWidth',1,'Color',rgb{MDP.s(3)})
plot(1 + 1/2,0,shape{MDP.s(6)},'MarkerSize',8,'LineWidth',1,'Color',rgb{MDP.s(4)})