text:
function [p,d] = check(X)
% CHECK(F)  Displays/calculates constraint residuals on constraint F
%
% [pres,dres] = CHECK(F)
%
% pres : Primal constraint residuals
% dres : Dual constraint residuals
%
% If no output argument is supplied, tabulated results are displayed
%
% Primal constraint residuals are calculated as:
%
%  Semidefinite constraint F(x)>0 : min(eig(F))
%  Element-wise constraint F(x)>0 : min(min(F))
%  Equality constraint F==0       : -max(max(abs(F)))
%  Second order cone t>||x||      : t-||x||
%  Integrality constraint on x    : max(abs(x-round(x)))
%  Sum-of-square constraint       : Minus value of largest (absolute value) coefficient
%                                   in the polynomial p-v'*v
%
% Dual constraints are evaluated similarily.
%
%  See also  SOLVESDP, SOLVESOS, SOSD, DUAL

switch nargout
    case 0
        check(lmi(X));
    case 1
        p = check(lmi(X));
    case 2
        [p,d] = check(lmi(X));
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%% Transformation Texture
%
%%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function AInv=invert2X2Matrix(A)
%%INVERT2X2MATRIX Inver the 2X2 matrix A. This function demonstrates the
%       algorithm given in [1], which is more efficient than Gaussian
%       elimination. This could be used as a template for implementation in
%       other languages. In Matlab, this function is not faster than the
%       built-in inv function.
%
%INPUTS: A A 2X2 real or complex invertible matrix.
%
%OUTPUTS: AInv The 2X2 matrix inverse of A.
%
%EXAMPLE:
%This just shows that this function produces the same result as Matlab
%within finite precision limitations.
% A=randn(2,2);
% C=invert2X2Matrix(A);
% C1=inv(A);%Matlab's way.
% RelErr=max(max(abs(abs((C-C1)./C1))))
%
%REFERENCES:
%[1] V. Strassen, "Gaussian elimination is not optimal," Numerische
%    Mathematik, vol. 13, pp. 354-356, 1965.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

I=1/A(1,1);
II=A(2,1)*I;
III=I*A(1,2);
IV=A(2,1)*III;
V=IV-A(2,2);
VI=1/V;

AInv=zeros(2,2);
AInv(1,2)=III*VI;
AInv(2,1)=VI*II;
VII=III*AInv(2,1);
AInv(1,1)=I-VII;
AInv(2,2)=-VI;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y = any(x)

% ALL    True if all elements of a vector are nonzero.
%    For vectors, ALL(V) returns logical 1 (TRUE) if none of the elements
%    of the vector are zero.  Otherwise it returns logical 0 (FALSE).  For
%    matrices, ALL(X) operates on the columns of X, returning a row vector
%    of logical 1's and 0's. For N-D arrays, ALL(X) operates on the first
%    non-singleton dimension.

if nargin>1
  error('this implementation is only supported with one input argument');
end

siz = size(x);
if numel(siz)>2
  error('this implementation is only supported with vector or matrix input');
end

if siz(1)==1 || siz(2)==1
  y = true;
  for i=1:prod(siz)
    if ~x(i)
      y = false;
      break
    end
  end

else
  y = true(1,siz(2));
  for j=1:siz(2)
    for i=1:siz(1)
      if ~x(i,j)
        y(1,j) = false;
        break
      end
    end
  end

end

author: PatternRecognition

repo: OpenBMI


text:
classdef WrapperMshFile < FileReader

    properties (Access = private)
        dataBase
        lineNumber0
        tScan
        dimension
        nnode
    end

    properties (Access = private)
        lineNumber
    end

    methods (Access = public)

        function obj = WrapperMshFile(cParams)
            obj.init(cParams)
        end

        function read(obj)
            obj.openFile();
            obj.scanFile();
            obj.readDimensionAndNnode();
            obj.readCoordinates();
            obj.readConnec();
            obj.closeFile();
        end

        function d = getDataBase(obj)
            d = obj.dataBase;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.filePath = cParams.filePath;
            obj.lineNumber0 = 7;
        end

        function scanFile(obj)
            t = textscan(obj.fid,'%s','delimiter','\n', 'headerlines',0);
            obj.tScan = t{1};
        end

        function readDimensionAndNnode(obj)
            lineDimension = 5;
            tLine = obj.readLines(lineDimension);
            tLinesS = split(tLine{1});
            tLinesD = str2double(tLinesS);
            obj.dimension = tLinesD(4);
            obj.nnode = tLinesD(8);
        end

        function readCoordinates(obj)
            iline = 0;
            line0 = obj.lineNumber0;
            isCoord = true;
            xNodes = zeros(1,obj.dimension);
            while isCoord
                iline = iline + 1;
                tLine = obj.readLines(iline + line0);
                if strcmp(tLine{1},'end coordinates ')
                    isCoord = false;
                else
                    isCoord = true;
                    tLinesS = split(tLine{1});
                    tLinesD = str2double(tLinesS);
                    xNodes(iline,:) = tLinesD(2:end-1,1);
                end
            end
            obj.dataBase.coord = xNodes;
        end

        function readConnec(obj)
            line0 = obj.lineNumber0 + size(obj.dataBase.coord,1) + 3;
            iline = 0;
            isCoord = true;
            conn = zeros(1,obj.nnode);
            while isCoord
                iline = iline + 1;
                tLine = obj.readLines(iline + line0);
                if strcmp(tLine{1},'end elements')
                    isCoord = false;
                else
                    isCoord = true;
                    tLinesS = split(tLine{1});
                    tLinesD = str2double(tLinesS);
                    conn(iline,:) = tLinesD(2:obj.nnode+1,1);
                end
            end
            obj.dataBase.connec = conn;
        end

        function tLine = readLines(obj,lineNumber,nLines)
            tLine = obj.tScan(lineNumber,:);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [cfg] = ft_topoplotTFR(cfg, varargin)

% FT_TOPOPLOTTFR plots the topographic distribution over the head
% of a 3-dimensional data representations such as time-frequency
% representation of the power or coherence spectrum.
%
% Use as
%   ft_topoplotTFR(cfg, freq)
%
% The input freq structrure should contain a time-resolved power or
% coherence spectrum from FT_FREQANALYSIS or FT_FREQDESCRIPTIVES.
%
% The configuration can have the following parameters
%   cfg.parameter          = field that contains the data to be plotted as color, for example 'avg', 'powspctrm' or 'cohspctrm' (default is automatic)
%   cfg.maskparameter      = field in the data to be used for masking of data. It should have alues between 0 and 1, where 0 corresponds to transparent.
%   cfg.xlim               = limit for 1st dimension in data (e.g., time), can be 'maxmin' or [xmin xmax] (default = 'maxmin')
%   cfg.ylim               = limit for 2nd dimension in data (e.g., freq), can be 'maxmin' or [ymin ymax] (default = 'maxmin')
%   cfg.zlim               = limits for color dimension, 'maxmin', 'maxabs', 'zeromax', 'minzero', or [zmin zmax] (default = 'maxmin')
%   cfg.channel            = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.refchannel         = name of reference channel for visualising connectivity, can be 'gui'
%   cfg.baseline           = 'yes','no' or [time1 time2] (default = 'no'), see FT_TIMELOCKBASELINE or FT_FREQBASELINE
%   cfg.baselinetype       = 'absolute' or 'relative' (default = 'absolute')
%   cfg.trials             = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.colormap           = string, or Nx3 matrix, see FT_COLORMAP
%   cfg.marker             = 'on', 'labels', 'numbers', 'off'
%   cfg.markersymbol       = channel marker symbol (default = 'o')
%   cfg.markercolor        = channel marker color (default = [0 0 0] (black))
%   cfg.markersize         = channel marker size (default = 2)
%   cfg.markerfontsize     = font size of channel labels (default = 8 pt)
%   cfg.highlight          = 'off', 'on', 'labels', 'numbers'
%   cfg.highlightchannel   =  Nx1 cell-array with selection of channels, or vector containing channel indices see FT_CHANNELSELECTION
%   cfg.highlightsymbol    = highlight marker symbol (default = 'o')
%   cfg.highlightcolor     = highlight marker color (default = [0 0 0] (black))
%   cfg.highlightsize      = highlight marker size (default = 6)
%   cfg.highlightfontsize  = highlight marker size (default = 8)
%   cfg.hotkeys            = enables hotkeys (pageup/pagedown/m) for dynamic zoom and translation (ctrl+) of the color limits
%   cfg.colorbar           = 'yes'
%                            'no' (default)
%                            'North'              inside plot box near top
%                            'South'              inside bottom
%                            'East'               inside right
%                            'West'               inside left
%                            'NorthOutside'       outside plot box near top
%                            'SouthOutside'       outside bottom
%                            'EastOutside'        outside right
%                            'WestOutside'        outside left
%   cfg.colorbartext       =  string indicating the text next to colorbar
%   cfg.interplimits       = limits for interpolation (default = 'head')
%                            'electrodes' to furthest electrode
%                            'head' to edge of head
%   cfg.interpolation      = 'linear','cubic','nearest','v4' (default = 'v4') see GRIDDATA
%   cfg.style              = plot style (default = 'both')
%                            'straight' colormap only
%                            'contour' contour lines only
%                            'both' (default) both colormap and contour lines
%                            'fill' constant color between lines
%                            'blank' only the head shape
%   cfg.gridscale          = scaling grid size (default = 67)
%                            determines resolution of figure
%   cfg.shading            = 'flat' or 'interp' (default = 'flat')
%   cfg.comment            = 'no', 'auto' or 'xlim' (default = 'auto')
%                            'auto': date, xparam, yparam and parameter limits are printed
%                            'xlim': only xparam limits are printed
%                            'ylim': only yparam limits are printed
%   cfg.commentpos         = string or two numbers, position of the comment (default = 'leftbottom')
%                            'lefttop' 'leftbottom' 'middletop' 'middlebottom' 'righttop' 'rightbottom'
%                            'title' to place comment as title
%                            'layout' to place comment as specified for COMNT in layout
%                            [x y] coordinates
%   cfg.interactive        = Interactive plot 'yes' or 'no' (default = 'yes')
%                            In a interactive plot you can select areas and produce a new
%                            interactive plot when a selected area is clicked. Multiple areas
%                            can be selected by holding down the SHIFT key.
%   cfg.directionality     = '', 'inflow' or 'outflow' specifies for
%                            connectivity measures whether the inflow into a
%                            node, or the outflow from a node is plotted. The
%                            (default) behavior of this option depends on the dimor
%                            of the input data (see below).
%   cfg.layout             = specify the channel layout for plotting using one of
%                            the supported ways (see below).
%   cfg.interpolatenan     = string 'yes', 'no' (default = 'yes')
%                            interpolate over channels containing NaNs
%
% For the plotting of directional connectivity data the cfg.directionality
% option determines what is plotted. The default value and the supported
% functionality depend on the dimord of the input data. If the input data
% is of dimord 'chan_chan_XXX', the value of directionality determines
% whether, given the reference channel(s), the columns (inflow), or rows
% (outflow) are selected for plotting. In this situation the default is
% 'inflow'. Note that for undirected measures, inflow and outflow should
% give the same output. If the input data is of dimord 'chancmb_XXX', the
% value of directionality determines whether the rows in data.labelcmb are
% selected. With 'inflow' the rows are selected if the refchannel(s) occur in
% the right column, with 'outflow' the rows are selected if the
% refchannel(s) occur in the left column of the labelcmb-field. Default in
% this case is '', which means that all rows are selected in which the
% refchannel(s) occur. This is to robustly support linearly indexed
% undirected connectivity metrics. In the situation where undirected
% connectivity measures are linearly indexed, specifying 'inflow' or
% 'outflow' can result in unexpected behavior.
%
% The layout defines how the channels are arranged. You can specify the
% layout in a variety of ways:
%  - you can provide a pre-computed layout structure (see prepare_layout)
%  - you can give the name of an ascii layout file with extension *.lay
%  - you can give the name of an electrode file
%  - you can give an electrode definition, i.e. "elec" structure
%  - you can give a gradiometer definition, i.e. "grad" structure
% If you do not specify any of these and the data structure contains an
% electrode or gradiometer structure, that will be used for creating a
% layout. If you want to have more fine-grained control over the layout
% of the subplots, you should create your own layout file.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure. For this particular function, the input should be
% structured as a cell-array.
%
% See also FT_TOPOPLOTER, FT_TOPOPLOTIC, FT_SINGLEPLOTTFR, FT_MULTIPLOTTFR, FT_PREPARE_LAYOUT

% Undocumented options:
%
% It is possible to use multiple highlight-selections (e.g.: multiple
% statistical clusters of channels) To do this, all the content of
% the highlight-options (including cfg.highlight) should be placed
% in a cell-array (even if the normal content was already in a
% cell-array). Specific marker settings (e.g. color, size) are defaulted
% when not present.
%
% Example (3 selections):
% cfg.highlight          = {'labels', 'labels', 'numbers'}
% cfg.highlightchannel   = {{'MZF03','MZC01','MRT54'}, [1:5], 'C*'}
% cfg.highlightsymbol    = {'o',[],'+'}        % the empty option will be defaulted
% cfg.highlightcolor     = {'r',[0 0 1]};      % the missing option will be defaulted
% cfg.highlightsize      = [];                 % will be set to default, as will the missing cfg.highlightfontsize
%
% Other options:
% cfg.labeloffset (offset of labels to their marker, default = 0.005)

% Copyright (C) 2005-2017, F.C. Donders Centre
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPERS NOTE: This code is organized in a similar fashion for multiplot/singleplot/topoplot
% and for ER/TFR and should remain consistent over those 6 functions.
% Section 1: general cfg handling that is independent from the data
% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data
% Section 3: cfg handling that depends on the data
% Section 4: actual plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Section 1: general cfg handling that is independent from the data

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar varargin
ft_preamble provenance varargin

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% make sure figure window titles are labeled appropriately, pass this onto the actual plotting function
cfg.funcname = mfilename;
cfg.dataname = dataname;

% prepare the layout, this should be done only once
tmpcfg = keepfields(cfg, {'layout', 'channel', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
cfg.layout = ft_prepare_layout(tmpcfg, varargin{1});

% call the common function that is shared between ft_topoplotER and ft_topoplotTFR
[cfg] = topoplot_common(cfg, varargin{:});

% remove this field again, it is only used for figure labels
cfg = removefields(cfg, 'funcname');

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous varargin
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function planC = calculateGRE(baseScanNum,movScanNum,planC)
% function planC = calculateGRE(baseScanNum,movScanNum,planC)
%
% APA, 03/21/2017

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

% Absolute Difference between two scans
siz = size(planC{indexS.scan}(baseScanNum).scanArray);
baseMask3M = logical(maskByThresh3D(planC{indexS.scan}(baseScanNum).scanArray));
movMask3M = logical(maskByThresh3D(planC{indexS.scan}(movScanNum).scanArray));
sA1 = zeros(siz,'single');
meanSa1 = mean(single(planC{indexS.scan}(baseScanNum).scanArray(baseMask3M)));
sdSa1 = std(single(planC{indexS.scan}(baseScanNum).scanArray(baseMask3M)));
sA1(baseMask3M) = (single(planC{indexS.scan}(baseScanNum).scanArray(baseMask3M)) - meanSa1)/sdSa1;
sA2 = zeros(siz,'single');
meanSa2 = mean(single(planC{indexS.scan}(movScanNum).scanArray(movMask3M)));
sdSa2 = std(single(planC{indexS.scan}(movScanNum).scanArray(movMask3M)));
sA2(movMask3M) = (single(planC{indexS.scan}(movScanNum).scanArray(movMask3M)) - meanSa2)/sdSa2;
diff3M = abs(sA1 - sA2);

% Window size
slcWindow = 5;
rowWindow = 5;
colWindow = 5;

% Number of levels for histogram
numLevels = 16;
% Create initial imM with dimension of slcWindow
imM = [];
for slc = 1:slcWindow
    imTmpM = im2col(diff3M(:,:,slc),[rowWindow colWindow],'sliding');
    imM = [imM;imTmpM];
end

numNeighbors = rowWindow*colWindow;
numSlcs = size(diff3M,3);
numRows = size(diff3M,1);
numCols = size(diff3M,2);
entropy3M = zeros(siz,'single');
mean3M = zeros(siz,'single');
var3M = zeros(siz,'single');
for slc = 1:numSlcs
    disp(['-------------------', num2str(slc)])
    if slc > floor(slcWindow/2) && slc <= (numSlcs-floor(slcWindow/2))
        imM(1:numNeighbors,:) = [];
        imTmpM = im2col(diff3M(:,:,slc),[rowWindow colWindow],'sliding');
        imM = [imM;imTmpM];
    end
    countsM = hist(imM,numLevels);
    countsM = countsM/numNeighbors/colWindow;
    entrpy2M = col2im(-sum(countsM.*log2(countsM+eps)),[rowWindow colWindow],[numRows numCols],'sliding');
    mean2M = col2im(mean(imM),[rowWindow colWindow],[numRows numCols],'sliding');
    var2M = col2im(var(imM),[rowWindow colWindow],[numRows numCols],'sliding');
    for i = 1:floor(slcWindow/2)
        entrpy2M = [entrpy2M(:,1), entrpy2M, entrpy2M(:,end)];
        entrpy2M = [entrpy2M(1,:); entrpy2M; entrpy2M(end,:)];
        mean2M = [mean2M(:,1), mean2M, mean2M(:,end)];
        mean2M = [mean2M(1,:); mean2M; mean2M(end,:)];
        var2M = [var2M(:,1), var2M, var2M(:,end)];
        var2M = [var2M(1,:); var2M; var2M(end,:)];
    end
    entropy3M(:,:,slc) = entrpy2M;
    mean3M(:,:,slc) = mean2M;
    var3M(:,:,slc) = var2M;
end

gre3M = entropy3M/median(entropy3M(baseMask3M)) + mean3M/median(mean3M(baseMask3M)) + var3M/median(var3M(baseMask3M));
gre3M = gre3M / 3;

%
% showIMDose(diff3M,'Diff',1);
% showIMDose(entropy3M,'Entropy',1);
% showIMDose(mean3M,'Mean',1);
% showIMDose(var3M,'Variance',1);
% showIMDose(gre3M,'GRE',1);

register = 'UniformCT';  %Currently only option supported.  Dose has the same shape as the uniformized CT scan.
doseError = [];
doseEdition = 'Generalized Registration Error';
overWrite = 'no';  %Overwrite the last CERR dose?
if ~exist('assocScanNum','var')
    assocScanNum = 1;
end
fractionGroupID = 'GRE';
assocScanUID = planC{indexS.scan}(assocScanNum).scanUID;
description = '';
planC = dose2CERR(gre3M,doseError,fractionGroupID,doseEdition,description,register,[],overWrite,assocScanUID,planC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%
%   run an example of ordfilt3 on a noisy sphere. The example is the same as
%   that written by Olivier Salvado.

clear
randn('seed',0)

[x,y,z] = meshgrid(1:100,1:100,1:100);
sphere = 20 + 200*( ( sqrt((x-50).^2+(y-50).^2+(z-50).^2) ) < 40 );
clear x y z
sphere = uint8(1*sphere + 50*randn(size(sphere)));

% -- median filter
tic
[Vr] = ordfilt3(sphere,'med',5);
toc
clf
% -- compare to box filter
subplot(221)
p1 = patch(isosurface(sphere,100), ...
   'FaceColor','blue','EdgeColor','none');
p2 = patch(isocaps(sphere,100), ...
    'FaceColor','interp','EdgeColor','none');
isonormals(sphere,p1)
view(3); axis vis3d square
camlight; lighting phong

subplot(222)
p1 = patch(isosurface(Vr,100), ...
   'FaceColor','blue','EdgeColor','none');
p2 = patch(isocaps(Vr,100), ...
    'FaceColor','interp','EdgeColor','none');
isonormals(Vr,p1)
view(3); axis vis3d square
camlight; lighting phong

%%
% show some slice
for k=1:20,
    subplot(223)
    imagesc( sphere(:,:,k) ,[0 255]),axis image
    title('original')

    subplot(224)
    imagesc( Vr(:,:,k) ,[0 255]),axis image
    title('Filtered with a 3D median filter')

    drawnow
end

author: Sable

repo: mcbench-benchmarks


text:
% Test file for chebtech/real.m

function pass = test_real(pref)

if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    % Test a scalar-valued function:
    f = testclass.make(@(x) exp(1i*x) + 1i*sin(x), [], pref);
    g = testclass.make(@(x) cos(x), [], pref);
    h = real(f);
    pass(n, 1) = norm(h.coeffs - g.coeffs, inf) < 10*vscale(h)*eps;

    % Test an array-valued function:
    f = testclass.make(@(x) [exp(1i*x) + 1i*sin(x), -exp(1i*x)], [], pref);
    g = testclass.make(@(x) [cos(x), -real(exp(1i*x))], [], pref);
    h = real(f);
    pass(n, 2) = norm(h.coeffs - g.coeffs, inf) < 10*max(vscale(h)*eps);

    % Test a real function:
    f = testclass.make(@(x) 1i*cos(x), [], pref);
    g = real(f);
    pass(n, 3) = numel(g.coeffs) == 1 && g.coeffs == 0;

    % Test an array-valued real function:
    f = testclass.make(@(x) 1i*[cos(x), sin(x), exp(x)], [], pref);
    g = real(f);
    pass(n, 4) = all(size(g.coeffs) == [1, 3]) && all(g.coeffs == 0);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function linplus_test37 ( )

%*****************************************************************************80
%
%% TEST37 tests R8GE_NP_ML.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST37\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_NP_ML computes A*x or A''*X\n' );
  fprintf ( 1, '    where A has been factored by R8GE_NP_FA.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  for job = 0 : 1
%
%  Set the matrix.
%
    [ a, seed ] = r8ge_random ( n, n, seed );
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r8ge_mxv ( n, n, a, x );
    else
      b = r8ge_vxm ( n, n, a, x );
    end
%
%  Factor the matrix.
%
    [ a_lu, info ] = r8ge_np_fa ( n, a );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TEST37 - Fatal error!\n' );
      fprintf ( 1, '  R8GE_NP_FA declares the matrix is singular!\n' );
      fprintf ( 1, '  The value of INFO is %d\n', info );
      continue;
    end
%
%  Now multiply factored matrix times solution to get right hand side again.
%
    b2 = r8ge_np_ml ( n, a_lu, x, job );

    if ( job == 0 )
      r8vec2_print_some ( n, b, b2, 10, '  A*x and PLU*x' );
    else
      r8vec2_print_some ( n, b, b2, 10, '  A''*x and (PLU)''*x' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [sos,P] = sospolymatrixvar(sos,ZSym,n,matrixstr)
% SOSMATRIXVAR --- Declare a polynomial matrix variable P in the sos program
% SOS of
%
% [SOSP,P] = sosmatrixvar(SOSP,ZSym,n,matrixstr)
%
% SOSP is the sum of squares program.
% P is the new polynomial matrix.
% n is the dimension of the matrix P: n(1) x n(2)
% ZSym is the vector of monomials contained in VAR. Decision
% variables corresponding to those monomials will be assigned
% automatically by SOSPOLYVAR.
% matrixstr is a char string with the option 'symmetric' when required

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.

% JA&GV - 6/13/2013

% Original Code
if nargin == 4
    if matrixstr=='symmetric'
        if n(1)==n(2)
            if isfield(sos,'symvartable')
                P = sym(zeros(n(1),n(2)));
            else
                % Code for multipoly: PJS 9/9/2013
                P = polynomial(zeros(n(1),n(2)));
            end
            for i = 1:n(1)
                for j = i:n(1)
                    [sos,var] = sospolyvar(sos,ZSym);
                    P(i,j) = var;
                    P(j,i) = var;
                    clear var;
                end
            end
        else
            disp(['''symmetric''' ' option used, matrix must be square.']);
            P = [];
            return
        end
    else
        disp(['Matrix structure ' matrixstr ' is not defined.' ]);
        P = [];
        return

    end
else
    if isfield(sos,'symvartable')
        P = sym(zeros(n(1),n(2)));
    else
        % Code for multipoly: PJS 9/9/2013
        P = polynomial(zeros(n(1),n(2)));
    end
    for i = 1:n(1)
        for j = 1:n(2)
            [sos,var] = sospolyvar(sos,ZSym);
            P(i,j) = var;
            clear var;
        end
    end

end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%op_CSIIntegrate.m
%Brenden Kadota, Sunnybrook 2021.
%
% USAGE:
% in=op_CSIintegrate(in);
%
% DESCRIPTION:
% Integrates the spectrum from ppmmin to ppmmax in all csi voxels.
%
% INPUTS:
% in        = CSI FID-A data structure
% pmmmin    = lower integration bounds
% ppmmax    = upper integration bounds
% mode      = mode (optional):
%                        -'re' (integral performed on real part (default)).
%                        -'im' (integral performed on imaginary part).
%                        -'mag' (integral performed on magnitude part).
%
% OUTPUTS:
% out       = map of integrated area
function [map, signal, noised] = op_CSIgetSNR(MRSIStruct, NAAppmmin, NAAppmmax, noiseppmmin, noiseppmmax)
    arguments
        MRSIStruct (1, 1) struct
        NAAppmmin (1, 1) double = nan
        NAAppmmax (1, 1) double = nan
        noiseppmmin (1, 1) double = nan
        noiseppmmax (1, 1) double = nan
    end
    checkArguments(MRSIStruct);
    arguments = removeNanArguments(NAAppmmin, NAAppmmax, noiseppmmin, noiseppmmax);

    MRSIStruct = reshapeDimensions(MRSIStruct, {'y', 'x', 't'});
    map = zeros(getSizeFromDimensions(MRSIStruct, {'y', 'x', 'extras'}));
    signal = zeros(getSizeFromDimensions(MRSIStruct, {'y', 'x', 'extras'}));
    noised = zeros(getSizeFromDimensions(MRSIStruct, {'y', 'x', 'extras'}));
    for e = 1:getSizeFromDimensions(MRSIStruct, {'extras'})
        for x = 1:getSizeFromDimensions(MRSIStruct, {'x'})
            for y = 1:getSizeFromDimensions(MRSIStruct, {'y'})
                voxel = op_CSItoMRS(MRSIStruct, x, y, "Extra", e);
                [voxSNR, voxSignal, voxNoiseD] = op_getSNR(voxel, arguments{:});
                signal(y, x, e) = voxSignal;
                noised(y, x, e) = voxNoiseD;
                map(y, x, e) = voxSNR;
            end
        end
    end
end

function checkArguments(in)
    if(in.flags.spectralFT == 0)
        error('FID-A Error: Input type invalid. Please fourier transform along the spectral dimension');
    end
    if(in.flags.spatialFT == 0)
        error('FID-A Error: Input type invalid. Please fourier transform along the spacial dimension');
    end
end

function args = removeNanArguments(NAAppmmin, NAAppmmax, noiseppmmin, noiseppmmax)
    args = {NAAppmmin, NAAppmmax, noiseppmmin, noiseppmmax};
    emptyIndecies = zeros(1, size(args, 2), 'logical');
    for iArgument = 1:length(args)
        if(isnan(args{iArgument}))
            emptyIndecies(iArgument) = 1;
        end
    end
    args(emptyIndecies) = [];
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function opts = mkopts_bj(T);

opts.algorithm = 'bj';
opts.use_kd_tree = 0;
opts.sis = 20;
opts.initial_K = T;
opts.do_greedy = 0;
opts.do_split = 0;
opts.do_merge = 0;
opts.do_sort = 0;

% Local Variables: ***
% mode: matlab ***
% End: ***

repo: REST


text:
function [dynpcm2] = bar2dynpcm2(bar)
% Convert pressure from bar to dynes/cm^2.
% Chad Greene 2012
dynpcm2 = bar*1.00000e+6;


text:
function out = compute_hash_function(words,hash_key)

if (nargin==1)
  rand('state',0);
  hash_key = rand(1,1000);
end

if ~iscell(words)
  words = {words};
end

for a=1:length(words)

  out(a) = sum( hash_key(1:length(words{a})) .* double(uint8(words{a})) );

end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function spm_dcm_Granger_demo
% Demo routine for induced responses
%==========================================================================
%
% This routine illustrates the relationship between Geweke Granger
% causality (GC) in frequency space and modulation transfer functions
% (MTF).  We first compare and contrast analytic results for GC with
% estimates based on a simulated time series. These synthetic data are
% chosen to show that (analytic) GC can, in principle, detect sparsity
% structure in terms of missing causal connections (however, GC estimates
% are not so efficient). We then demonstrate the behaviour of (analytic)
% GC by varying the strength of forward connections, backward connections
% and intrinsic gain.  There is reasonable behaviour under these
% manipulations. However, when we introduce realistic levels of (power law)
% measurement noise, GC fails. The simulations conclude by showing that DCM
% recovery of the underlying model parameters can furnish  (analytic) GC
% among sources (in the absence of measurement noise). [delete the 'return'
% below to see these simulations].
%
% See also:
%  spm_ccf2csd.m, spm_ccf2mar, spm_csd2ccf.m, spm_csd2mar.m, spm_mar2csd.m,
%  spm_csd2coh.m, spm_ccf2gew, spm_dcm_mtf.m, spm_Q.m, spm_mar.m and
%  spm_mar_spectral.m
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_Granger_demo.m 6059 2014-06-19 11:57:31Z vladimir $

% Model specification
%==========================================================================
rng('default')

% number of regions
%--------------------------------------------------------------------------
Nc    = 2;                                       % number of channels
Ns    = 2;                                       % number of sources
ns    = 2*128;                                   % sampling frequency
dt    = 1/ns;                                    % time bins
Hz    = 1:128;                                   % frequency
p     = 16;                                      % autoregression order
options.spatial  = 'LFP';
options.model    = 'CMC';
options.analysis = 'CSD';
M.dipfit.model = options.model;
M.dipfit.type  = options.spatial;
M.dipfit.Nc    = Nc;
M.dipfit.Ns    = Ns;
M.pF.D         = [1 4];

% extrinsic connections (forward an backward)
%--------------------------------------------------------------------------
A{1} = [0 0; 0 0];
A{2} = [0 0; 0 0];
A{3} = [0 0; 0 0];
B    = {};
C    = sparse(2,0);

% get priors
%--------------------------------------------------------------------------
pE    = spm_dcm_neural_priors(A,B,C,options.model);
pE    = spm_L_priors(M.dipfit,pE);
pE    = spm_ssr_priors(pE);
x     = spm_dcm_x_neural(pE,options.model);

% (log) connectivity parameters
%--------------------------------------------------------------------------
pE.A{1}(2,1) = 2;
pE.S         = 1/8;

% (log) amplitude of fluctations and noise
%--------------------------------------------------------------------------
pE.a(1,:) = -2;
pE.b(1,:) = -8;
pE.c(1,:) = -8;

% orders and model
%==========================================================================
nx    = length(spm_vec(x));

% create forward model
%--------------------------------------------------------------------------
M.f   = 'spm_fx_cmc';
M.g   = 'spm_gx_erp';
M.x   = x;
M.n   = nx;
M.pE  = pE;
M.m   = Ns;
M.l   = Nc;
M.Hz  = Hz;
M.Rft = 4;

% specify M.u - endogenous input (fluctuations) and intial states
%--------------------------------------------------------------------------
M.u   = sparse(Ns,1);

% solve for steady state
%--------------------------------------------------------------------------
M.x   = spm_dcm_neural_x(pE,M);

% Analytic spectral chararacterisation
%==========================================================================
spm_figure('GetWin','Figure 1'); clf

[csd,Hz,mtf] = spm_csd_mtf(pE,M);
csd          = csd{1};
mtf          = mtf{1};
ccf          = spm_csd2ccf(csd,Hz,dt);
mar          = spm_ccf2mar(ccf,p);
mar          = spm_mar_spectra(mar,Hz,ns);

spm_figure('GetWin','Figure 1'); clf
spm_spectral_plot(Hz,csd,  'b', 'frequency','density')
spm_spectral_plot(Hz,mar.P,'r', 'frequency','density')

legend('cross spectral density',...
       'autoregressive model')

%  comparison of expected results
%==========================================================================
spm_figure('GetWin','Figure 2'); clf

dtf  = mar.dtf;
gew  = mar.gew;
new  = spm_csd2gew(csd,Hz);

spm_spectral_plot(Hz,mtf,'b', 'frequency','density')
spm_spectral_plot(Hz,dtf,'g', 'frequency','density')
spm_spectral_plot(Hz,gew,'r', 'frequency','density')
spm_spectral_plot(Hz,new,'r:','frequency','density')

subplot(2,2,3), a = axis; subplot(2,2,2), axis(a);

legend('modulation transfer function',...
       'directed transfer function',...
       'Granger causality (parametric)',...
       'Granger causality (non-parametric)')

% effects of changing various model parameters
%==========================================================================

% (log) scaling, and parameters
%--------------------------------------------------------------------------
logs  = [ ((1:4)/1 - 2);
          ((1:4)/1 - 2);
          ((1:4)/6 + 0);
          ((1:4)/1 - 6)];

param = {'A{1}(2,1)','A{3}(1,2)','S','c(1,2)'};
str   = {     'forward connectivity',
              'backward connectivity',
              'intrinsic gain',
              'amplitude of noise'};

% expected transfer function and Gramger causality
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3'); clf
a     = [0 Hz(end) 0 .35];
ca    = 0;
for i = 1:size(logs,1)
    for j = 1:size(logs,2)

        P     = pE;
        eval(sprintf('P.%s = %i;',param{i},logs(i,j)));

        % create forward model and solve for steady state
        %------------------------------------------------------------------
        M.x   = spm_dcm_neural_x(P,M);

        % Analytic spectral chararacterisation
        %==================================================================
        [csd,Hz,mtf,qew] = spm_csd_mtf(P,M);

        ccf   = spm_csd2ccf(csd{1},Hz,dt);
        mar   = spm_ccf2mar(ccf,p);
        mar   = spm_mar_spectra(mar,Hz,ns);

        % plot forwards and backwards functions
        %------------------------------------------------------------------
        subplot(4,2,ca + 1)
        plot(Hz,abs(mar.gew(:,2,1)),Hz,abs(qew{1}(:,2,1)))
        xlabel('frequency')
        ylabel('absolute value')
        title(sprintf('%s',str{i}),'FontSize',16)
        axis square, hold on, set(gca,'XLim',[0 Hz(end)])
        axis(a);

        subplot(4,2,ca + 2)
        plot(Hz,abs(mar.gew(:,1,2)),Hz,abs(qew{1}(:,1,2)))
        xlabel('frequency')
        ylabel('absolute value')
        title(sprintf('backward'),'FontSize',16)
        axis square, hold on, set(gca,'XLim',[0 Hz(end)])
        axis(a);

    end
    ca  = ca + 2;

end

% Lyapunov functions
%==========================================================================

% (log) scaling, and parameters
%--------------------------------------------------------------------------
param = {'A{1}(2,1)','A{3}(1,2)','S','D(1,1)'};
str   = {     'forward connectivity',
              'backward connectivity',
              'intrinsic gain',
              'intrinsic delays'};

% changing parameters and Lyapunov functions
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4'); clf
k     = linspace(0,1,32);
for i = 1:length(param)
    for j = 1:length(k)

        P     = pE;
        eval(sprintf('P.%s = %i;',param{i},k(j)));

        % create forward model and solve for steady state
        %------------------------------------------------------------------
        M.x   = spm_dcm_neural_x(P,M);

        % Analytic spectral chararacterisation
        %==================================================================
        [f,dfdx,D] = feval(M.f,M.x,M.u,P,M);
        dfdx       = D*dfdx;
        S(i,j)     = max(real(eig(full(dfdx),'nobalance')));

        % recursive check on MAR approximation
        %------------------------------------------------------------------
        csd        = spm_csd_mtf(P,M);
        ccf        = spm_csd2ccf(csd{1},Hz);
        [mar,c]    = spm_ccf2mar(ccf,p);
        pcond(i,j) = c;

    end
end

% plot
%--------------------------------------------------------------------------
subplot(2,2,1)
plot(k,S),hold on
xlabel('parameter (log) scaling')
ylabel('principal eigenvalue')
title('Lyapunov exponents','FontSize',16)
axis square; spm_axis tight
plot(k,k*0,'k:',k,k*0 - 1/(dt*p),'k-.'),hold off

subplot(2,2,2)
semilogy(k,pcond),hold on
xlabel('parameter (log) scaling')
ylabel('condition number')
title('condition of cross covariance','FontSize',16)
axis square;
legend(str)

plot(k,k*0 + 1e4,'k:'),hold off

% MAR approximation
%==========================================================================
param = {'S'};
for i = 1:length(param)
    for j = 1:length(k)

        P     = pE;
        eval(sprintf('P.%s = %i;',param{i},k(j)));

        % create forward model and solve for steady state
        %------------------------------------------------------------------
        M.x   = spm_dcm_neural_x(P,M);

        % check on MAR approximation
        %------------------------------------------------------------------
        csd   = spm_csd_mtf(P,M);
        mar   = spm_csd2mar(csd{1},Hz,p);
        mar   = spm_mar_spectra(mar,Hz);

        CSD(:,j) = csd{1}(:,2,2);
        PSD(:,j) = mar.P(:,2,2);

    end
end

subplot(2,2,3)
imagesc(k,Hz,log(abs(CSD)))
xlabel('parameter (log) scaling')
ylabel('frequency')
title('(log) spectral density','FontSize',16)
axis square

subplot(2,2,4)
imagesc(k,Hz,abs(PSD - CSD).^(1/4))
xlabel('parameter (log) scaling')
ylabel('frequency')
title('MAR approximation error','FontSize',16)
axis square

% a more careful examination of [in]stability on Granger causality
%==========================================================================
spm_figure('GetWin','Figure 5'); clf
k     = linspace(0,2/3,8);
for j = 1:length(k)

    % amplitude of observation noise
    %----------------------------------------------------------------------
    P        = pE;
    P.S      = k(j);

    % create forward model and solve for steady state
    %----------------------------------------------------------------------
    M.x      = spm_dcm_neural_x(P,M);

    % Analytic spectral chararacterisation (parametric)
    %======================================================================
    [csd,Hz,mtf,qew] = spm_csd_mtf(P,M);

    ccf      = spm_csd2ccf(csd{1},Hz,dt);
    mar      = spm_ccf2mar(ccf,p);
    mar      = spm_mar_spectra(mar,Hz,ns);

    % and non-parametric
    %======================================================================
    gew      = spm_csd2gew(csd{1},Hz);

    % save forwards and backwards functions
    %----------------------------------------------------------------------
    GCF(:,j) = abs(mar.gew(:,2,1));
    GCB(:,j) = abs(mar.gew(:,1,2));

    % plot forwards and backwards functions (parametric)
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 5');

    subplot(3,2,1)
    plot(Hz,abs(mar.gew(:,2,1)),Hz,abs(qew{1}(:,2,1)))
    xlabel('frequency')
    ylabel('absolute value')
    title('forward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    a  = axis;

    subplot(3,2,2)
    plot(Hz,abs(mar.gew(:,1,2)),Hz,abs(qew{1}(:,1,2)))
    xlabel('frequency')
    ylabel('absolute value')
    title('backward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    % plot forwards and backwards functions (nonparametric)
    %----------------------------------------------------------------------
    subplot(3,2,3)
    plot(Hz,abs(gew(:,2,1)),Hz,abs(qew{1}(:,2,1)))
    xlabel('frequency')
    ylabel('absolute value')
    title('forward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    subplot(3,2,4)
    plot(Hz,abs(gew(:,1,2)),Hz,abs(qew{1}(:,1,2)))
    xlabel('frequency')
    ylabel('absolute value')
    title('backward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

end

a = .2;
subplot(3,2,5)
image(Hz,k,GCF'*64/a)
xlabel('frequency')
ylabel('log(exponent)')
title('forward','FontSize',16)
axis square

subplot(3,2,6)
image(Hz,k,GCB'*64/a)
xlabel('frequency')
ylabel('log(exponent)')
title('backward','FontSize',16)
axis square

% a more careful examination of measurement noise
%==========================================================================
spm_figure('GetWin','Figure 6' ); clf
spm_figure('GetWin','Figure 6a'); clf
k     = linspace(-8,-0,8);
for j = 1:length(k)

    % amplitude of observation noise
    %----------------------------------------------------------------------
    P        = pE;
    P.b(1)   = k(j);
    P.c(1,:) = k(j) + 1;

    % create forward model and solve for steady state
    %----------------------------------------------------------------------
    M.x      = spm_dcm_neural_x(P,M);

    % Analytic spectral chararacterisation (parametric)
    %======================================================================
    [csd,Hz,mtf,qew] = spm_csd_mtf(P,M);

    ccf      = spm_csd2ccf(csd{1},Hz,dt);
    mar      = spm_ccf2mar(ccf,p);
    mar      = spm_mar_spectra(mar,Hz,ns);

    % and non-parametric
    %======================================================================
    gew      = spm_csd2gew(csd{1},Hz);
    qew      = qew{1};

    % save forwards and backwards functions
    %----------------------------------------------------------------------
    GCF(:,j) = abs(gew(:,2,1));
    GCB(:,j) = abs(gew(:,1,2));

    % plot forwards and backwards functions  (parametric)
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 6');

    subplot(3,2,1)
    plot(Hz,abs(mar.gew(:,2,1)),Hz,abs(qew(:,2,1)))
    xlabel('frequency')
    ylabel('absolute value')
    title('forward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    a  = axis;

    subplot(3,2,2)
    plot(Hz,abs(mar.gew(:,1,2)),Hz,abs(qew(:,1,2)))
    xlabel('frequency')
    ylabel('absolute value')
    title('backward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    % plot forwards and backwards functions  (nonparametric)
    %----------------------------------------------------------------------
    subplot(3,2,3)
    plot(Hz,abs(gew(:,2,1)),Hz,abs(qew(:,2,1)))
    xlabel('frequency')
    ylabel('absolute value')
    title('forward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    subplot(3,2,4)
    plot(Hz,abs(gew(:,1,2)),Hz,abs(qew(:,1,2)))
    xlabel('frequency')
    ylabel('absolute value')
    title('backward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    % plot associated coherence
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 6a' );

    if j > 4, str = ':'; else, str = '-'; end
    coh = spm_csd2coh(csd{1},Hz);
    spm_spectral_plot(Hz,coh,str,'frequency','coherence')

end

spm_figure('GetWin','Figure 6' );

subplot(3,2,5)
imagesc(Hz,k,GCF')
xlabel('frequency')
ylabel('log(exponent)')
title('forward','FontSize',16)
axis square

subplot(3,2,6)
imagesc(Hz,k,GCB')
xlabel('frequency')
ylabel('log(exponent)')
title('backward','FontSize',16)
axis square

% return if in demonstration mode
%--------------------------------------------------------------------------
DEMO = 1;
if DEMO, return, end

% DCM estimates of coupling
%==========================================================================
rng('default')

% get priors and generate data
%--------------------------------------------------------------------------
pE    = spm_dcm_neural_priors(A,B,C,options.model);
pE    = spm_L_priors(M.dipfit,pE);
pE    = spm_ssr_priors(pE);

% (log) connectivity parameters (forward connection only)
%--------------------------------------------------------------------------
pE.A{1}(2,1) = 2;
pE.S         = 1/8;

% (log) amplitude of fluctations and noise
%--------------------------------------------------------------------------
pE.a(1,:) = -2;
pE.b(1,:) = -8;
pE.c(1,:) = -8;

% expected cross spectral density
%--------------------------------------------------------------------------
csd       = spm_csd_mtf(pE,M);

% Get spectral profile of fluctuations and noise
%--------------------------------------------------------------------------
[Gu,Gs,Gn] = spm_csd_mtf_gu(pE,Hz);

% Integrate with power law process (simulate multiple trials)
%--------------------------------------------------------------------------
PSD   = 0;
CSD   = 0;
N     = 1024;
U.dt  = dt;
for t = 1:16

    % neuronal fluctuations
    %----------------------------------------------------------------------
    U.u      = spm_rand_power_law(Gu,Hz,dt,N);
    LFP      = spm_int_L(pE,M,U);

    % and measurement noise
    %----------------------------------------------------------------------
    En       = spm_rand_power_law(Gn,Hz,dt,N);
    Es       = spm_rand_power_law(Gs,Hz,dt,N);
    E        = Es + En*ones(1,Ns);

    % and estimate spectral features under a MAR model
    %----------------------------------------------------------------------
    MAR      = spm_mar(LFP + E,p);
    MAR      = spm_mar_spectra(MAR,Hz,ns);
    CSD      = CSD + MAR.P;

    % and using Welch's method
    %----------------------------------------------------------------------
    PSD      = PSD + spm_csd(LFP + E,Hz,ns);

    CCD(:,t) = abs(CSD(:,1,2)/t);
    PCD(:,t) = abs(CSD(:,1,1)/t);

    % plot
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 7'); clf
    spm_spectral_plot(Hz,csd{1},'r',  'frequency','density')
    spm_spectral_plot(Hz,CSD/t, 'b',  'frequency','density')
    spm_spectral_plot(Hz,PSD/t, 'g',  'frequency','density')
    legend('real','estimated (AR)','estimated (CSD)')
    drawnow

end

%  show convergence of spectral estimators
%--------------------------------------------------------------------------
subplot(2,2,3), hold off
imagesc(Hz,1:t,log(PCD'))
xlabel('frequency')
ylabel('trial number')
title('log auto spectra','FontSize',16)
axis square

subplot(2,2,4), hold off
imagesc(Hz,1:t,log(CCD'))
xlabel('frequency')
ylabel('trial number')
title('log cross spectra','FontSize',16)
axis square

% DCM set up (allow for both forward and backward connections)
%==========================================================================

% (log) connectivity parameters (forward connection only)
%--------------------------------------------------------------------------
pE.A{1}(2,1) = 2;
pE.S         = 1/8;

% (log) amplitude of fluctations and noise
%--------------------------------------------------------------------------
pE.a(1,:) = -2;
pE.b(1,:) = -4;
pE.c(1,:) = -4;

% expected cross spectral density
%--------------------------------------------------------------------------
[csd,Hz,mtf] = spm_csd_mtf(pE,M);

DCM.options.model   = 'CMC';
DCM.options.spatial = 'LFP';
DCM.options.DATA    = 0;

DCM.M.dipfit.Nc  = Nc;
DCM.M.dipfit.Ns  = Ns;
DCM.M.Nmax = 32;
DCM.M.U    = eye(Nc,Nc);

DCM.A      = {[0 0; 1 0],[0 1; 0 0],[0 0; 0 0]};
DCM.B      = {};
DCM.C      = sparse(Ns,0);

% place in data structure
%--------------------------------------------------------------------------
DCM.xY.y  = csd;
DCM.xY.dt = dt;
DCM.xY.Hz = Hz;

% estimate
%--------------------------------------------------------------------------
DCM       = spm_dcm_csd(DCM);

% show results in terms of transfer functions and Granger causality
%==========================================================================
spm_figure('GetWin','Figure 8'); clf

% transfer functions in the absence of measurement noise
%--------------------------------------------------------------------------
Ep    = DCM.Ep;
Ep.b  = Ep.b - 32;              % and suppress non-specific and
Ep.c  = Ep.c - 32;              % specific channel noise

Gu           = spm_csd_mtf_gu(Ep,Hz);
[psd Hz dtf] = spm_csd_mtf(Ep,DCM.M);

mtf   = mtf{1};
csd   = csd{1};
dtf   = dtf{1};

% transfer functions and Granger causality among sources and channels
%--------------------------------------------------------------------------
gwe  = spm_dtf2gew(mtf,Gu);
qew  = spm_dtf2gew(dtf,Gu);
gew  = spm_csd2gew(csd,Hz);

spm_spectral_plot(Hz,mtf,'b:', 'frequency','density')
spm_spectral_plot(Hz,gwe,'b',  'frequency','density')
spm_spectral_plot(Hz,qew,'r',  'frequency','density')
spm_spectral_plot(Hz,gew,'g',  'frequency','density')
legend('modulation transfer function',...
       'Granger causality (true)',...
       'Granger causality (source)',...
       'Granger causality (channel)')

subplot(2,2,3), a = axis; subplot(2,2,2), axis(a);

return

% NOTES: a more careful examination of delays
%==========================================================================
spm_figure('GetWin','Figure 9'); clf

k     = linspace(8,11,8);
for j = 1:length(k)

    % keep total power of fluctuations constant
    %----------------------------------------------------------------------
    M.pF.D(2,1)  = k(j);

    % create forward model and solve for steady state
    %----------------------------------------------------------------------
    M.x          = spm_dcm_neural_x(pE,M);

    % Analytic spectral chararacterisation (parametric)
    %======================================================================
    [csd,Hz,mtf] = spm_csd_mtf(pE,M);
    ccf          = spm_csd2ccf(csd{1},Hz,dt);
    mar          = spm_ccf2mar(ccf,p);
    mar          = spm_mar_spectra(mar,Hz,ns);

    % and non-parametric)
    %======================================================================
    gew          = spm_csd2gew(csd{1},Hz);

    % plot forwards and backwards functions
    %----------------------------------------------------------------------
    subplot(3,2,1)
    plot(Hz,abs(mar.gew(:,2,1)),Hz,abs(mtf{1}(:,2,1)))
    xlabel('frequency')
    ylabel('absolute value')
    title('forward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    a  = axis;

    subplot(3,2,2)
    plot(Hz,abs(mar.gew(:,1,2)),Hz,abs(mtf{1}(:,1,2)))
    xlabel('frequency')
    ylabel('absolute value')
    title('backward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    % plot forwards and backwards functions
    %----------------------------------------------------------------------
    subplot(3,2,3)
    plot(Hz,abs(gew(:,2,1)),Hz,abs(mtf{1}(:,2,1)))
    xlabel('frequency')
    ylabel('absolute value')
    title('forward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

    subplot(3,2,4)
    plot(Hz,abs(gew(:,1,2)),Hz,abs(mtf{1}(:,1,2)))
    xlabel('frequency')
    ylabel('absolute value')
    title('backward','FontSize',16)
    axis square, hold on, set(gca,'XLim',[0 Hz(end)])
    axis(a);

end

author: spm


text:
clear all; close all;
I=imread('cameraman.tif');
J=imnoise(I, 'poisson');
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function CPD = hhmmQ_CPD(bnet, self, Qnodes, d, D, varargin)
% HHMMQ_CPD Make the CPD for a Q node at depth D of a D-level hierarchical HMM
% CPD = hhmmQ_CPD(bnet, self, Qnodes, d, D, ...)
%
%  Fd(t-1) \   Q1:d-1(t)
%           \  |
%            \ v
%  Qd(t-1) -> Qd(t)
%            /
%           /
%  Fd+1(t-1)
%
% We assume parents are ordered (numbered) as follows:
% Qd(t-1), Fd+1(t-1), Fd(t-1), Q1(t), ..., Qd(t)
%
% The parents of Qd(t) can either be just Qd-1(t) or the whole stack Q1:d-1(t) (allQ)
% In either case, we will call them Qps.
% If d=1, Qps does not exist. Also, the F1(t-1) -> Q1(t) arc is optional.
% If the arc is missing, startprob does not need to be specified,
% since the toplevel is assumed to never reset (F1 does not exist).
% If d=D, Fd+1(t-1) does not exist (there is no signal from below).
%
% optional args [defaults]
%
% transprob - transprob(i,k,j) = prob transition from i to j given Qps = k ['leftright']
% selfprob  - prob of a transition from i to i given Qps=k [0.1]
% startprob - startprob(k,j) = prob start in j given Qps = k ['leftstart']
% startargs - other args to be passed to the sub tabular_CPD for learning startprob
% transargs - other args will be passed to the sub tabular_CPD for learning transprob
% allQ      - 1 means use all Q nodes above d as parents, 0 means just level d-1 [0]
% F1toQ1    - 1 means add F1(t-1) -> Q1(t) arc, 0 means level 1 never resets [0]
%
% For d=1, startprob(1,j) is only needed if F1toQ1=1
% Also, transprob(i,j) can be used instead of transprob(i,1,j).
%
% hhmmQ_CPD is a subclass of tabular_CPD so we inherit inference methods like CPD_to_pot, etc.
%
% We create isolated tabular_CPDs with no F parents to learn transprob/startprob
% so we can avail of e.g., entropic or Dirichlet priors.
% In the future, we will be able to represent the transprob using a tree_CPD.
%
% For details, see "Linear-time inference in hierarchical HMMs", Murphy and Paskin, NIPS'01.

ss = bnet.nnodes_per_slice;
%assert(self == Qnodes(d)+ss);
ns = bnet.node_sizes(:);
CPD.Qsizes = ns(Qnodes);
CPD.d = d;
CPD.D = D;
allQ = 0;

% find out which parents to use, to get right size
for i=1:2:length(varargin)
  switch varargin{i},
   case 'allQ', allQ = varargin{i+1};
  end
end

if d==1
  CPD.Qps = [];
else
  if allQ
    CPD.Qps = Qnodes(1:d-1);
  else
    CPD.Qps = Qnodes(d-1);
  end
end

Qsz = ns(self);
Qpsz = prod(ns(CPD.Qps));

% set default arguments
startprob = 'leftstart';
transprob = 'leftright';
startargs = {};
transargs = {};
CPD.F1toQ1 = 0;
selfprob = 0.1;

for i=1:2:length(varargin)
  switch varargin{i},
   case 'transprob', transprob = varargin{i+1};
   case 'selfprob',  selfprob = varargin{i+1};
   case 'startprob', startprob = varargin{i+1};
   case 'startargs', startargs = varargin{i+1};
   case 'transargs', transargs = varargin{i+1};
   case 'F1toQ1',    CPD.F1toQ1 = varargin{i+1};
  end
end

Qps = CPD.Qps + ss;
old_self = self-ss;

if strcmp(transprob, 'leftright')
  LR = mk_leftright_transmat(Qsz, selfprob);
  transprob = repmat(reshape(LR, [1 Qsz Qsz]), [Qpsz 1 1]); % transprob(k,i,j)
  transprob = permute(transprob, [2 1 3]); % now transprob(i,k,j)
end
transargs{end+1} = 'CPT';
transargs{end+1} = transprob;
CPD.sub_CPD_trans = mk_isolated_tabular_CPD([old_self Qps], ns([old_self Qps self]), transargs);
S = struct(CPD.sub_CPD_trans);
CPD.transprob = myreshape(S.CPT, [Qsz Qpsz Qsz]);

if strcmp(startprob, 'leftstart')
  startprob = zeros(Qpsz, Qsz);
  startprob(:,1) = 1;
end

if (d==1) & ~CPD.F1toQ1
  CPD.sub_CPD_start = [];
  CPD.startprob = [];
else
  startargs{end+1} = 'CPT';
  startargs{end+1} = startprob;
  CPD.sub_CPD_start = mk_isolated_tabular_CPD(Qps, ns([Qps self]), startargs);
  S = struct(CPD.sub_CPD_start);
  CPD.startprob = myreshape(S.CPT, [Qpsz Qsz]);
end

CPD = class(CPD, 'hhmmQ_CPD', tabular_CPD(bnet, self));

CPD = update_CPT(CPD);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function test_tutorial_plotting

% WALLTIME 00:15:00
% MEM 4gb
% DEPENDENCY ft_multiplotER ft_singleplotER ft_topoplotER ft_singleplotTFR ft_multiplotTRF ft_megplanar ft_combineplanar ft_volumereslice

% this tutorial reflects the plotting tutorial up to March 2017
% around that time we did a data visualization workshop at the donders for which I created an updated tutorial

% see http://www.fieldtriptoolbox.org/tutorial/plotting
% this testscript corresponds to the version on the wiki at 23 December 2012

% use the tutorial dataset from home/common

load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/avgFC.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/GA_FC.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/TFRhann.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/statERF.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/statTFR.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/sourceDiff.mat'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% Singleplot Functions %%%%%%%%%
% ft_singleplotER
cfg = [];
cfg.xlim = [-0.2 1.0];
cfg.ylim = [-1e-13 3e-13];
cfg.channel = 'MLC24';
clf;
ft_singleplotER(cfg,avgFC);

% MATLAB toolbox alternative of ft_singleplotER
selected_data = avgFC.avg(9,241:601); %MLC24 is the 9th channel, -0.2 to 1.0 is sample 241 to 601
time = avgFC.time(241:601);
figure;
plot(time, selected_data)
xlim([-0.2 1.0])
ylim([-1e-13 3e-13])

% ft_singleplotTFR
cfg = [];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim = [-1.5e-27 1.5e-27];
cfg.channel       = 'MLC*';
cfg.renderer      = 'painters'; % might be needed if colorbar does not print properly
figure;ft_singleplotTFR(cfg, TFRhann); colorbar;

cfg.channel       = 'MRC*';
figure;ft_singleplotTFR(cfg, TFRhann); colorbar;

%%%%%%%%% Topoplot Functions %%%%%%%%%
cfg = [];
cfg.xlim = [0.3 0.5];
cfg.zlim = [0 6e-14];
cfg.layout = 'CTF151.lay';
cfg.parameter = 'individual'; % the default (avg) is not present
figure; ft_topoplotER(cfg,GA_FC); colorbar;

cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [-1e-27 1e-27];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';
figure; ft_topoplotTFR(cfg,TFRhann)

% multiple topoplots in one figure
cfg.xlim = -0.4:0.2:1.4;
cfg.comment = 'xlim';
cfg.commentpos = 'title';
figure; ft_topoplotTFR(cfg,TFRhann)

% data selection options prior to calling topoplot
% not specific to topoplot
cfg = [];
cfg.xlim = [0.9 1.3];
cfg.ylim = [15 20];
cfg.zlim = [-1e-27 1e-27];
cfg.baseline = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';

% specific to topoplot - option 1 (regular colours)
cfg.gridscale = 300;
cfg.style = 'straight';
cfg.marker = 'labels';
figure; ft_topoplotTFR(cfg,TFRhann)

% specific to topoplot - option 2 (gray scale)
cfg.gridscale = 300;
cfg.contournum = 10;
cfg.colormap = gray(10);
figure; ft_topoplotTFR(cfg,TFRhann)

% specific to topoplot - option 3 (spring colormap)
cfg.gridscale = 300;
cfg.contournum = 4;
cfg.colormap = spring(4);
cfg.markersymbol = '.';
cfg.markersize = 12;
cfg.markercolor = [0 0.69 0.94];
figure; ft_topoplotTFR(cfg,TFRhann)

%%%%%%%%% Plotting with interactive mode / multiplots  %%%%%%%%%

cfg = [];
cfg.baseline = [-0.5 -0.1];
cfg.zlim = [-3e-27 3e-27];
cfg.baselinetype = 'absolute';
cfg.layout = 'CTF151.lay';
cfg.interactive = 'yes';
figure; ft_multiplotTFR(cfg,TFRhann)

cfg = [];
cfg.xlim = [-0.2 1.0];
cfg.ylim = [-1e-13 3e-13];
cfg.layout = 'CTF151.lay';
cfg.interactive = 'yes';
figure; ft_multiplotER(cfg,avgFC);

%%%%%%%%% Connectivity plots  %%%%%%%%%

cfg = [];
cfg.covariance = 'yes';
tmp1 = ft_timelockanalysis(cfg, avgFC);

cfg = [];
cfg.method = 'corr';
tmp2 = ft_connectivityanalysis(cfg, tmp1);

cohFC      = tmp2;
cohFC.coh  = tmp2.corr;
cohFC.freq = 1;
cohFC = rmfield(cohFC, 'coh');

cfg = [];
cfg.foi = 1;
cfg.layout = 'CTF151.lay';
figure; ft_topoplotCC(cfg, cohFC);

%%%%%%%%% Cluster plots  %%%%%%%%%

% timelocked data
cfg = [];
cfg.zlim = [-6 6]; %Tvalues
cfg.alpha = 0.05;
cfg.layout = 'CTF151.lay';
figure; ft_clusterplot(cfg,statERF);

% freq data
cfg = [];
cfg.zlim = [-5 5];
cfg.alpha = 0.05;
cfg.layout = 'CTF151.lay';
ft_clusterplot(cfg,statTFR);

%%%%%%%% beamforming plotting  %%%%%%%

% (1) multiple axial slices

%% load MRI and interpolate functional source data to MRI
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/plotting/Subject01.mri'));
mri = ft_volumereslice([], mri);

cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'pow';
sourceDiffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri);

%% plot multiple 2D axial slices
cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
ft_sourceplot(cfg, sourceDiffInt);

%% (2) ortho plot
% volume normalise
cfg = [];
cfg.coordsys      = 'ctf';
cfg.nonlinear     = 'no';
sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);

% plot ortho
cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
figure; ft_sourceplot(cfg, sourceDiffIntNorm);

%% (3) MNI surface
cfg = [];
cfg.method         = 'surface';
cfg.funparameter   = 'pow';
cfg.maskparameter  = cfg.funparameter;
cfg.funcolorlim    = [0.0 1.2];
cfg.funcolormap    = 'jet';
cfg.opacitylim     = [0.0 1.2];
cfg.opacitymap     = 'rampup';
cfg.projmethod     = 'nearest';
cfg.surffile      = 'surface_white_both.mat';
cfg.surfdownsample = 10;
ft_sourceplot(cfg, sourceDiffIntNorm);
view ([90 0])


text:
function prefs = mrPrefs;
% Prompts user to adjust MR object-related preferences. Initializes defaults
% if they don't exist.
%
% prefs = mrPrefs;
%
% The user will usually only be prompted if the preferences don't yet
% exist. However, if you call this function with no return arg, it will
% always prompt.
%
% STILL UNDER CONSTRUCTION
%
% ras, 04/2006.
if(nargout==0)
    promptUser = true;
else
    promptUser = false;
end

% get default values for prefs --
% to edit new prefs, define them in this sub-routine.
dlg = defaultPrefsDialog;
prefNames = {dlg.fieldName}; % names of all prefs to be defined

% check that all methods are defined
if ispref('mr')
    prefs = getpref('mr');
    notAssigned = setdiff(prefNames, fieldnames(prefs));
    if ~isempty(notAssigned)
        promptUser = true;
    end
else
    % need to initialize
    disp('Initializing mrVista2 MR object preferences')
    dlg = defaultPrefsDialog;
    for i = 1:length(dlg)
        setpref('mr', dlg(i).fieldName, dlg(i).value);
    end
    prefs = getpref('mr');
    promptUser = true;
end

% put up a dialog, if requested or necessary
if promptUser
    % plug in any already-defined prefs over the default values
    if exist('prefs') & ~isempty(prefs)
        for f = fieldnames(prefs)
            ii = cellfind(prefNames, f{1});
            if ii>0 % this is a valid pref to set
                dlg(ii).value = prefs.(f{1});
            end
        end
    end

    prefs = generalDialog(dlg, 'mrVista2 MR Preferences');

    % exit quietly if canceled -- otherwise update prefs
    if isempty(prefs)
        return;
    else
        for f = fieldnames(prefs)'
            setpref('mr', f{1}, prefs.(f{1}));
        end
    end
end

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function dlg = defaultPrefsDialog;
% return the default MR prefs names/values as a dialog. This seems
% to be the easiest way to describe new preferences.
dlg(1).fieldName = 'interp';
dlg(1).style = 'popup';
dlg(1).string = 'Interpolation Method for Different Coordinate Spaces?';
dlg(1).list = {'Nearest-Neighbor' 'Linear' 'Cubic' 'Spline'};
dlg(1).value = 'Linear';

dlg(end+1).fieldName = 'alphaMethod';
dlg(end).style = 'popup';
dlg(end).string = 'Alpha Method for Multiple Overlays?';
dlg(end).list = {'Average: [1 0 0] + [.2 0 1] => [.6 0 .5]' ...
                'Add and saturate: [1 0 0] + [.2 0 1] => [1 0 1]' ...
                'Opaque: [1 0 0] + [.2 0 1] => [.2 0 1]'};
dlg(end).value = 1;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% show the eigenface images
% examples are conducted in Yale face database
[vectors,Training_matrix,Y]=train(20);  % 10 is the number of eigenvectors
% the size of vectors is 77760*10
for k=1:20
    col=vectors(:,k);  %pick up every column vector
    im=reshape(col,243,320); % transform this vector into 2D-image
    m_a=max(max(im));
    n_b=min(min(im));
    for i=1:243
       for j=1:320
       im(i,j)=(im(i,j)-n_b)/(m_a-n_b)*255;
       end
    end
    subplot(4,5,k);
    imshow(uint8(im));
end

author: sfvsfv

repo: Mathematical-modeling


text:
function [dataout] = mmse_det(ch_and_rx)
% mmse_det  User-defined function for MMSE detector.
% Input: 'ch_and_rx'
%    - first elements are MIMO parameters (STBC, Nss, etc.)
%    - next 16*NST elements are channel estimates.
%    - next (NSTS+1)*NST elements are reference training symbols.
%    - next (NSTS+1)*NST elements are received training symbols.
%    - remaining x*NST elements are received data symbols.
% Output: 'detout'
%    - x*NST elements are MMSE detector outputs.
%    - x*NST elements are Viterbi metric weights.
global PER_snr snr_idx ch_type;

% Initialize vectors
Clk = ch_and_rx(1);
STBC = ch_and_rx(2);
Nss  = ch_and_rx(3);
Nsts = Nss + STBC;
Nltf = Nsts;
if (Nsts==3) % use 4 TRN symbols for NSTS=3
   Nltf = 4;
end

Nrx  = ch_and_rx(4);
ch_est  = ch_and_rx(4+(1:16*56));
ref_trn = ch_and_rx(4+16*56+(1:4*56*(Nltf)));
rx_trn  = ch_and_rx(4+16*56+4*56*(Nltf)+(1:4*56*(Nltf)));
rx_data = ch_and_rx((5+16*56+2*4*56*(Nltf)):end);

%%%% Form channel matrix
idx = 0;
H = zeros(4,4,56);
for c=1:4
    for r=1:4
        H(r,c,:) = ch_est(idx+(1:56));
        idx=idx+56;
    end
end
H = H(1:Nrx,1:Nsts,:);  %% channel matrix is Nrx x Nsts

%%%% Train data symbols...
idx = 0;
nsymtrn = (length(ref_trn)/4/56);
trnRef = zeros(Nrx,56,nsymtrn);
for n=1:nsymtrn
    for m=1:Nsts
        trnRef(m,:,n) = ref_trn(idx+(m-1)*56+(1:56));
        trnRx(m,:,n)  =  rx_trn(idx+(m-1)*56+(1:56));
    end
    idx=idx+4*56;
end

%%%% Rx data symbols...
idx = 0;
nsymdata = (length(rx_data)/4/56);
rdata = zeros(4,56,nsymdata);
for n=1:nsymdata
    for m=1:Nrx
        rdata(m,:,n) = rx_data(idx+(m-1)*56+(1:56));
    end
    idx=idx+4*56;
end

%%%%
%%%% STBC or MMSE detection...
%%%%
if (STBC)
    %%%% Form space-time decoding matrices
    for k=1:56
        if (STBC==1 && Nss==1)
            % Heff matrix (assumes channel stationary for 2 symbols)
            Heff(:,:,k) = [      H(1,1,k)      -H(1,2,k); ...
                            conj(H(1,2,k)) conj(H(1,1,k)) ];
        elseif (STBC==1 && Nss==2)
            % Heff matrix (assumes channel stationary for 2 symbols)
            Heff(:,:,k) = [      H(1,1,k)      -H(1,2,k)       H(1,3,k)              0;  ...
                            conj(H(1,2,k)) conj(H(1,1,k))             0  conj(H(1,3,k)); ...
                                 H(2,1,k)      -H(2,2,k)       H(2,3,k)              0;  ...
                            conj(H(2,2,k)) conj(H(2,1,k))             0  conj(H(2,3,k)) ];
        elseif (STBC==1 && Nss==3)
            % Heff matrix (assumes channel stationary for 2 symbols)
            Heff(:,:,k) = [      H(1,1,k)      -H(1,2,k)       H(1,3,k)              0        H(1,4,k)              0;  ...
                            conj(H(1,2,k)) conj(H(1,1,k))             0  conj(H(1,3,k))              0  conj(H(1,4,k)); ...
                                 H(2,1,k)      -H(2,2,k)       H(2,3,k)              0        H(2,4,k)              0;  ...
                            conj(H(2,2,k)) conj(H(2,1,k))             0  conj(H(2,3,k))              0  conj(H(2,4,k)); ...
                                 H(3,1,k)      -H(3,2,k)       H(3,3,k)              0        H(3,4,k)              0;  ...
                            conj(H(3,2,k)) conj(H(3,1,k))             0  conj(H(3,3,k))              0  conj(H(3,4,k)) ];
        elseif (STBC==2 && Nss==2)
            % Heff matrix (assumes channel stationary for 2 symbols)
            Heff(:,:,k) = [      H(1,1,k)      -H(1,2,k)       H(1,3,k)      -H(1,4,k);  ...
                            conj(H(1,2,k)) conj(H(1,1,k)) conj(H(1,4,k)) conj(H(1,3,k)); ...
                                 H(2,1,k)      -H(2,2,k)       H(2,3,k)      -H(2,4,k);  ...
                            conj(H(2,2,k)) conj(H(2,1,k)) conj(H(2,4,k)) conj(H(2,3,k)) ];
        end
    end

    %%%% Compute noise power...
    N_o = 0;
    for k=1:56
        He = Heff(:,:,k);
        N_o = N_o + (1/56)*sum(sum(He.*conj(He)))/Nrx;     % Signal power
    end
    snr_val = PER_snr(snr_idx);
    N_o = N_o * 10^(-snr_val/10);     % Adjust by SNR

    %%%% Compute space-time detector coeff's
    C = zeros(Nss*2,Nss*2,56);
    for k=1:56
        %%%% MMSE coeff's
        He = Heff(:,:,k);
        C(:,:,k) = inv(He'*He + N_o*eye(Nss*2))*He';
    end

    %%%% Compute space-time decoding outputs, Viterbi metric weights
    idx = 0;
    detout = zeros(4*56*nsymdata,1);
    metout = zeros(4*56*nsymdata,1);
    for n=1:2:nsymdata
        for k=1:56
            if (STBC==1 && Nss==1)
                % MMSE detector outputs
                detset = C(:,:,k) * [ rdata(1,k,n); ...
                                      conj(rdata(1,k,n+1)) ];
                detout(k+idx+0*56) = detset(1);         % first symbol in time, ss1
                detout(k+idx+4*56) = conj(detset(2));   % second symbol in time, ss1
                % Viterbi metric weights
                metset = eps+abs(1-diag(C(:,:,k)*Heff(:,:,k)));
                metout(k+idx+0*56) = 1/metset(1);
                metout(k+idx+4*56) = 1/metset(2);
            elseif (STBC==1 && Nss==2)
                % MMSE detector outputs
                detset = C(:,:,k) * [ rdata(1,k,n); ...
                                      conj(rdata(1,k,n+1)); ...
                                      rdata(2,k,n); ...
                                      conj(rdata(2,k,n+1)) ];
                detout(k+idx+0*56) = detset(1);         % first symbol in time, ss1
                detout(k+idx+4*56) = conj(detset(2));   % second symbol in time, ss1
                detout(k+idx+1*56) = detset(3);         % first symbol in time, ss2
                detout(k+idx+5*56) = conj(detset(4));   % second symbol in time, ss2
                % Viterbi metric weights
                metset = eps+abs(1-diag(C(:,:,k)*Heff(:,:,k)));
                metout(k+idx+0*56) = 1/metset(1);
                metout(k+idx+4*56) = 1/metset(2);
                metout(k+idx+1*56) = 1/metset(3);
                metout(k+idx+5*56) = 1/metset(4);
            elseif (STBC==1 && Nss==3)
                % MMSE detector outputs
                detset = C(:,:,k) * [ rdata(1,k,n); ...
                                      conj(rdata(1,k,n+1)); ...
                                      rdata(2,k,n); ...
                                      conj(rdata(2,k,n+1)); ...
                                      rdata(3,k,n); ...
                                      conj(rdata(3,k,n+1)) ];
                detout(k+idx+0*56) = detset(1);         % first symbol in time, ss1
                detout(k+idx+4*56) = conj(detset(2));   % second symbol in time, ss1
                detout(k+idx+1*56) = detset(3);         % first symbol in time, ss2
                detout(k+idx+5*56) = conj(detset(4));   % second symbol in time, ss2
                detout(k+idx+2*56) = detset(5);         % first symbol in time, ss3
                detout(k+idx+6*56) = conj(detset(6));   % second symbol in time, ss3
                % Viterbi metric weights
                metset = eps+abs(1-diag(C(:,:,k)*Heff(:,:,k)));
                metout(k+idx+0*56) = 1/metset(1);
                metout(k+idx+4*56) = 1/metset(2);
                metout(k+idx+1*56) = 1/metset(3);
                metout(k+idx+5*56) = 1/metset(4);
                metout(k+idx+2*56) = 1/metset(5);
                metout(k+idx+6*56) = 1/metset(6);
            elseif (STBC==2 && Nss==2)
                % MMSE detector outputs
                detset = C(:,:,k) * [ rdata(1,k,n); ...
                                      conj(rdata(1,k,n+1)); ...
                                      rdata(2,k,n); ...
                                      conj(rdata(2,k,n+1)) ];
                detout(k+idx+0*56) = detset(1);         % first symbol in time, ss1
                detout(k+idx+4*56) = conj(detset(2));   % second symbol in time, ss1
                detout(k+idx+1*56) = detset(3);         % first symbol in time, ss2
                detout(k+idx+5*56) = conj(detset(4));   % second symbol in time, ss2
                % Viterbi metric weights
                metset = eps+abs(1-diag(C(:,:,k)*Heff(:,:,k)));
                metout(k+idx+0*56) = 1/metset(1);
                metout(k+idx+4*56) = 1/metset(2);
                metout(k+idx+1*56) = 1/metset(3);
                metout(k+idx+5*56) = 1/metset(4);
            end
        end
        idx=idx+8*56;
    end

else

    %%%% Compute noise power...
    N_o = 0;
    for k=1:56
        He = H(:,:,k);
        N_o = N_o + (1/56)*sum(sum(He.*conj(He)))/Nrx;     % Signal power
    end
    snr_val = PER_snr(snr_idx);
    N_o = N_o * 10^(-snr_val/10);     % Adjust by SNR

    %%%%% MMSE linear detector (no STBC)
    C = zeros(Nsts,Nrx,56);
    for k=1:56
        %%%% MMSE coeff's
        He = H(:,:,k);
        C(:,:,k) = inv(He'*He + N_o*eye(Nrx)) * He';
    end

    %%%% MMSE detector outputs, Viterbi metric weights
    idx = 0;
    detset = zeros(4,1);
    detout = zeros(4*56*nsymdata,1);
    metset = eps+zeros(4,1);
    metout = zeros(4*56*nsymdata,1);
    for n=1:nsymdata
        for k=1:56
            % MMSE detector outputs
            detset(1:Nsts) = C(:,:,k)*rdata(1:Nrx,k,n);
            detout(k+idx+0*56) = detset(1);
            detout(k+idx+1*56) = detset(2);
            detout(k+idx+2*56) = detset(3);
            detout(k+idx+3*56) = detset(4);
            % Viterbi metric weights
            metset(1:Nsts) = eps+abs(1-diag(C(:,:,k)*H(:,:,k)));
            metout(k+idx+0*56) = 1/metset(1);
            metout(k+idx+1*56) = 1/metset(2);
            metout(k+idx+2*56) = 1/metset(3);
            metout(k+idx+3*56) = 1/metset(4);
        end
        idx=idx+4*56;
    end

end

% Adjust metric scaling range (0 - least confident, 7 - most confident)
% (based on use of 4-bit data, soft decisions, for Viterbi module)
max_metout = max(metout(find(metout<1e14)));
log_metout = log(metout + 1e-6);
log_max_metout = log(max_metout + 1e-6);
metout = floor(7.99 * log_metout/log_max_metout);
metout = (metout >= 7).*7 + (metout < 7).*metout;

% MMSE module output
dataout = [detout; metout];


text:
%
%   FILE NAME:    addManuallyCuratedRxnAssoc.m
%
%   PURPOSE: Incorporate manual curation results of rxn association
%            with MNXref ids to the mergedModel
%

% Move to the target folder of manual curation files
cd('/Users/haowa/Box Sync/HMR3/Curation Files/manualRxnCuration');

% Load the manual curation results
[~, textData1]=xlsread('conflictAssocByIDsAndViaMets.xlsx','ManualCuration_2018May');
[~, textData2]=xlsread('filledRxnAssocViaMets.xlsx','ManualCuration_2018May');
[~, textData3]=xlsread('rxnAssocOnlyByRxnIDs.xlsx','ManualCuration_2018May');

manuallyCuratedRxnAssoc.rxnHMRID=[textData1(2:end,1);textData2(2:end,1);textData3(2:end,1)];
manuallyCuratedRxnAssoc.rxnMNXID=[textData1(2:end,2);textData2(2:end,2);textData3(2:end,2)];

% Move back to rxnAssoc folder, and save this data structure
%save('manuallyCuratedRxnAssoc.mat','manuallyCuratedRxnAssoc');  % 2018-05-25

% Load the HMR2 model after merging compartments
load('mergedModel.mat');    % 2018-05-25

numel(find(~cellfun(@isempty,mergedModel.confirmedMNXID)))  %ans = 2506
% Add manually curated rxn associations to confirmedMNXID field, some of
% which are mostly not assigned (with empty values)
for i=1:numel(manuallyCuratedRxnAssoc.rxnHMRID)
		% No need to check the duplicateRxns field, because all reactions
		% for curation were selected from the rxns field
		[~, index]=ismember(manuallyCuratedRxnAssoc.rxnHMRID{i},mergedModel.rxns);
		if isempty(mergedModel.confirmedMNXID{index})
				mergedModel.confirmedMNXID{index}{1}=manuallyCuratedRxnAssoc.rxnMNXID{i};
		else
				if ~ismember(manuallyCuratedRxnAssoc.rxnMNXID{i},mergedModel.confirmedMNXID{index})
						mergedModel.confirmedMNXID{index}=[mergedModel.confirmedMNXID{index};manuallyCuratedRxnAssoc.rxnMNXID{i}];
				end
		end
end
numel(find(~cellfun(@isempty,mergedModel.confirmedMNXID)))  %ans = 2988
% A total of 492 (10 already have rxnAssoc) reactions in mergedModel were newly associated to MNX

% Re-processing the filed 'confirmedFilteredMNXID'
mergedModel.confirmedFilteredMNXID=filterBalancedRxns(mergedModel.confirmedMNXID);

% Save the mergedModel
%save('mergedModel.mat','mergedModel');  % 2018-05-25

%% sub functions
%Filter MNX rxns according to their balance status
function filteredMNXrxns = filterBalancedRxns(rxnList)

load('MNXRxns.mat');          % load MNX reactions
filteredMNXrxns=cell(numel(rxnList),1);
filteredMNXrxns(:)={''};

for i=1:numel(rxnList)
		if ~isempty(rxnList{i})
				countNum=numel(rxnList{i});
				if countNum==1  % there is only one MNX id, skip balance check
						filteredMNXrxns{i}=rxnList{i};
				elseif countNum>1   % there are several MNX ids
						[a, b]=ismember(rxnList{i},MNXRxns.MNX_ID);
						% The priority order of Balance status: true > NA > ambiguous
						if find(strcmp('true',MNXRxns.Balance(b)));
								trueHits=find(strcmp('true',MNXRxns.Balance(b)));
								filteredMNXrxns{i}=rxnList{i}(trueHits);
						elseif find(strcmp('NA',MNXRxns.Balance(b)))
								NAHits=find(strcmp('NA',MNXRxns.Balance(b)));
								filteredMNXrxns{i}=rxnList{i}(NAHits);
						elseif find(strcmp('ambiguous',MNXRxns.Balance(b)))
								amHits=find(strcmp('ambiguous',MNXRxns.Balance(b)));
								filteredMNXrxns{i}=rxnList{i}(amHits);
						else
								filteredMNXrxns{i}=rxnList{i};
						end
				end
		end
end

end

author: SysBioChalmers

repo: Human-GEM


text:
function showSpaceImages(img, x, y, scaling, limits)
%
% Displays images organized in a 2D space
%
% input:
%   img: array [nrows ncols 3 nimages]
%   x,y  = coordinate for each image

if nargin<5
    x1 = min(x(:));
    y1 = min(y(:));
    x2 = max(x(:));
    y2 = max(y(:));
else
    x1 = limits(1);
    x2 = limits(2);
    y1 = limits(3);
    y2 = limits(4);
end

%xa = x - min(x(:));
%ya = y - min(y(:));

if nargin<4
    scaling = 2;
end

%S = max(xa);
%xa = (0.975-0.05*scaling)*xa/S+0.025;
%S = max(ya);
%ya = (0.975-0.05*scaling)*ya/S+0.025;

xa = (0.975-0.05*scaling)*(x - x1)/(x2-x1) + 0.05*scaling;
ya = (0.975-0.05*scaling)*(y - y1)/(y2-y1) + 0.05*scaling;

figure
for n = length(xa):-1:1
    n
    h=axes('position', [xa(n) ya(n) .05*scaling .05*scaling]);
    image(uint8(img(:,:,:,n)), 'parent', h)
    axis('off'); axis('equal')
end

author: CSAILVision

repo: LabelMeToolbox


text:
function [rtk,stat0]=relpos(rtk,obsr,obsb,nav)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% relative positioning %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%input:  rtk - rtk control struct
%        obsr - rover observations
%        obsb - base observations
%        nav - navigation message
%output: rtk - rtk control struct
%        stat0 - state (0:error 1:ok)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%relative positioning include PPD, PPK and PPS
%PPD:post-processing differenced, based on double-differenced PR
%PPK:post-processing kinematic, based on double-differenced CP and PR
%PPS:post-processing static, based on double-differenced CP and PR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
opt=rtk.opt; nf=rtk.NF; dt=timediff(obsr(1).time,obsb(1).time);
nobsr=size(obsr,1); nobsb=size(obsb,1);

if rtk.opt.mode<=glc.PMODE_DGNSS
    stat=glc.SOLQ_DGNSS;
else
    stat=glc.SOLQ_FLOAT;
end

for i=1:glc.MAXSAT
    [sys,~]=satsys(i);
    rtk.sat(i).sys=sys;
    for j=1:glc.NFREQ
        rtk.sat(i).vsat(j)=0;
    end
    for j=1:glc.NFREQ
        rtk.sat(i).snr(j)=0;
    end
end

svr=satposs(obsr,nav,opt.sateph);
svb=satposs(obsb,nav,opt.sateph);

% zero-differenced residuals for base
[zdb,stat_tmp]=zdres(2,obsb,svb,nav,rtk.basepos,opt,rtk);
if stat_tmp==0,stat0=0;return;end

ind=selcomsat(obsr,obsb,opt,zdb);
if ind.ns==0,stat0=0;return;end

% debug tracing
% fprintf('\n\n\n');
% fprintf('time=%d\n',obsr(1).time.time);
% fprintf('before time update\n');
% printx(rtk.x,rtk);
% printP(rtk.P,rtk);

% time update
rtk=udstate(rtk,obsr,obsb,nav,ind);

% debug tracing
% fprintf('after time update\n');
% printx(rtk.x,rtk);
% printP(rtk.P,rtk);

xp=rtk.x; Pp=zeros(rtk.nx,rtk.nx);
for i=1:opt.niter

    % zero-differenced residuals for rover
    [zdr,stat_tmp]=zdres(1,obsr,svr,nav,xp(1:3),opt,rtk);
    if stat_tmp==0,stat=glc.SOLQ_NONE;break;end

    % double-differenced residuals,measurement matrix and noise matrix
    [rtk,v,H,R,nv,~]=ddres(rtk,nav,dt,xp,Pp,zdr,zdb,ind);
    if nv<=0,stat=glc.SOLQ_NONE;break;end

    % debug tracing
%     printv(v);
%     printH(H);
%     printR(R);

    % measurement update
    Pp=rtk.P;
    [xp,Pp,stat_tmp]=Kfilter_h(v,H,R,xp,Pp);
    if stat_tmp==0
        [week,sow]=time2gpst(obs(1).time);
        fprintf('Warning:GPS week = %d sow = %.3f,filter error!\n',week,sow);
        stat=glc.SOLQ_NONE;
        break;
    end

    % debug tracing
%     fprintf('after measurement update\n');
%     printx(xp,rtk);
%     printP(Pp,rtk);

end

if stat~=glc.SOLQ_NONE
    % zero-differenced residuals for rover
    [zdr,stat_tmp]=zdres(1,obsr,svr,nav,xp(1:3),opt,rtk);
    if stat_tmp~=0
        % post-fit residuals for float solution
        [rtk,v,~,R,nv,vflag]=ddres(rtk,nav,dt,xp,Pp,zdr,zdb,ind);
        % validation of float solution
        if valpos_rel(rtk,v,R,vflag,nv,4)
            rtk.x=xp; rtk.P=Pp;
            rtk.sol.ns=0;
            for i=1:ind.ns
                for f=1:nf
                    sati=ind.sat(i);
                    if ~rtk.sat(sati).vsat(f),continue;end
                    rtk.sat(sati).lock(f)=rtk.sat(sati).lock(f)+1;
                    rtk.sat(sati).outc(f)=0;
                    if f==1,rtk.sol.ns=rtk.sol.ns+1;end
                end
            end
            if rtk.sol.ns<4,stat=glc.SOLQ_NONE;end
        else
            stat=glc.SOLQ_NONE;
        end
    end
end

% ambiguity resolution
if stat~=glc.SOLQ_NONE
    [rtk,~,xa,nb]=resamb(rtk);
    if nb>1
        % zero-differenced residuals for rover
        [zdr,stat_tmp]=zdres(1,obsr,svr,nav,xa(1:3),opt,rtk);
        if stat_tmp~=0
            % post-fit reisiduals for fixed solution
            [rtk,v,~,R,nv,vflag]=ddres(rtk,nav,dt,xa,Pp,zdr,zdb,ind);
            % validation of fixed solution
            if valpos_rel(rtk,v,R,vflag,nv,4)
                rtk.nfix=rtk.nfix+1;
                if rtk.nfix>=opt.minfix&&rtk.opt.modear==glc.ARMODE_FIXHOLD
                    % hold integer ambiguity
                    rtk=holdamb(rtk,xa);
                end
                stat=glc.SOLQ_FIX;
            end
        end
    end
end

% save solution status
if stat==glc.SOLQ_FIX
    rtk.sol.pos=rtk.xa(1:3)';
    rtk.sol.posP(1)=rtk.Pa(1,1); rtk.sol.posP(2)=rtk.Pa(2,2);
    rtk.sol.posP(3)=rtk.Pa(3,3); rtk.sol.posP(4)=rtk.Pa(1,2);
    rtk.sol.posP(5)=rtk.Pa(2,3); rtk.sol.posP(6)=rtk.Pa(1,3);
    if rtk.opt.dynamics==1
        rtk.sol.vel=rtk.xa(4:6)';
        rtk.sol.velP(1)=rtk.Pa(4,4); rtk.sol.velP(2)=rtk.Pa(5,5);
        rtk.sol.velP(3)=rtk.Pa(6,6); rtk.sol.velP(4)=rtk.Pa(4,5);
        rtk.sol.velP(5)=rtk.Pa(5,6); rtk.sol.velP(6)=rtk.Pa(4,6);
    end
else
    rtk.sol.pos=rtk.x(1:3)';
    rtk.sol.posP(1)=rtk.P(1,1); rtk.sol.posP(2)=rtk.P(2,2);
    rtk.sol.posP(3)=rtk.P(3,3); rtk.sol.posP(4)=rtk.P(1,2);
    rtk.sol.posP(5)=rtk.P(2,3); rtk.sol.posP(6)=rtk.P(1,3);
    if rtk.opt.dynamics==1
        rtk.sol.vel=rtk.xa(4:6)';
        rtk.sol.velP(1)=rtk.P(4,4); rtk.sol.velP(2)=rtk.P(5,5);
        rtk.sol.velP(3)=rtk.P(6,6); rtk.sol.velP(4)=rtk.P(4,5);
        rtk.sol.velP(5)=rtk.P(5,6); rtk.sol.velP(6)=rtk.P(4,6);
    end
    rtk.nfix=0;
end

for i=1:nobsr
    for j=1:nf
        if obsr(i).L(j)==0,continue;end
        rtk.sat(obsr(i).sat).pt(1,j)=obsr(i).time;
        rtk.sat(obsr(i).sat).ph(1,j)=obsr(i).L(j);
    end
end
for i=1:nobsb
    for j=1:nf
        if obsb(i).L(j)==0,continue;end
        rtk.sat(obsb(i).sat).pt(2,j)=obsb(i).time;
        rtk.sat(obsb(i).sat).ph(2,j)=obsb(i).L(j);
    end
end

for i=1:ind.ns
    for j=1:nf
        rtk.sat(ind.sat(i)).snr(j)=obsr(ind.ir(i)).S(j);
    end
end

for i=1:glc.MAXSAT
    for j=1:nf
        if rtk.sat(i).fix(j)==2&&stat~=glc.SOLQ_FIX
            rtk.sat(i).fix(j)=1;
        end
        if bitand(rtk.sat(i).slip(j),1)
            rtk.sat(i).slipc(j)=rtk.sat(i).slipc(j)+1;
        end
    end
end

if stat~=glc.SOLQ_NONE,rtk.sol.stat=stat;end

stat0=(stat~=glc.SOLQ_NONE);

return

repo: GINav


text:
function neqn = p32_neqn ( )

%*****************************************************************************80
%
%% P32_NEQN returns the number of equations for problem P32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 3;

  return
end