text:
classdef FilterP0 < handle

   properties (Access = private)
       levelSet
       levelSet0
       dens0
   end

   methods (Access = public)

       function obj = FilterP0(ls,d)
           obj.levelSet = ls;
           obj.computeElementalLevelSet(d);
           obj.computeDensityP0();
       end

       function d = getDensity(obj)
           d = obj.dens0;
       end

   end

   methods (Access = private)

       function computeElementalLevelSet(obj,d)
            shape = d.shape;
            conec = d.conec;
            quadr = d.quadr;
            ngaus = quadr.ngaus;
            nelem = size(conec,1);
            nnode = size(shape,1);

            phiP0 = zeros(nelem,ngaus);
            phi   = obj.levelSet;
            for igaus = 1:ngaus
                for inode = 1:nnode
                    nodes = conec(:,inode);
                    phiN = phi(nodes);
                    phiP0(:,igaus) = phiP0(:,igaus) + shape(inode,igaus)*phiN;
                end
            end
            obj.levelSet0 = phiP0;
       end

        function computeDensityP0(obj)
            ls = obj.levelSet0;
            obj.dens0 = obj.computeDensity(ls);
        end

   end

   methods (Access = private, Static)

        function dens = computeDensity(phi)
            dens = 1 - heaviside(phi);
        end

   end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function pop = InitializeFun(pop_size, chromosome_size)
pop = round(rand(pop_size, chromosome_size));
end

author: LiYangSir

repo: Smart-Algorithm


text:
function [net1, avgImg] = initVGG16Net( vgg_16_path )
%VGG16
opts.gpus=1;
cold=true;
%prepareGPUs(opts,cold);

net2=load(vgg_16_path);
net2.layers(31:end) = [];

avgImg=net2.meta.normalization.averageImage;

net1=dagnn.DagNN();

net1.addLayer('conv1_1', dagnn.Conv('size', [3,3,3,64],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'input', 'conv_11', {'conv11_f', 'conv11_b'});
net1.addLayer('relu1_1', dagnn.ReLU(),'conv_11','relu_11');

f = net1.getParamIndex('conv11_f') ;
net1.params(f).value=net2.layers{1}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv11_b') ;
net1.params(f).value=net2.layers{1}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

net1.addLayer('conv1_2', dagnn.Conv('size', [3,3,64,64],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_11', 'conv_12', {'conv12_f', 'conv12_b'});
net1.addLayer('relu1_2', dagnn.ReLU(),'conv_12','relu_12');
net1.addLayer('pool1', dagnn.Pooling('poolSize', [2,2], 'pad',...
    [0,1,0,1], 'stride', [2,2]),'relu_12','pool_1');

f = net1.getParamIndex('conv12_f') ;
net1.params(f).value=net2.layers{3}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv12_b') ;
net1.params(f).value=net2.layers{3}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

%-------------------------covn 2--------------------
net1.addLayer('conv2_1', dagnn.Conv('size', [3,3,64,128],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'pool_1', 'conv_21', {'conv21_f', 'conv21_b'});
net1.addLayer('relu2_1', dagnn.ReLU(),'conv_21','relu_21');

f = net1.getParamIndex('conv21_f') ;
net1.params(f).value=net2.layers{6}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv21_b') ;
net1.params(f).value=net2.layers{6}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

net1.addLayer('conv2_2', dagnn.Conv('size', [3,3,128,128],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_21', 'conv_22', {'conv22_f', 'conv22_b'});
net1.addLayer('relu2_2', dagnn.ReLU(),'conv_22','relu_22');
net1.addLayer('pool2', dagnn.Pooling('poolSize', [2,2], 'pad',...
    [0,1,0,1], 'stride', [2,2]),'relu_22','pool_2');

f = net1.getParamIndex('conv22_f') ;
net1.params(f).value=net2.layers{8}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv22_b') ;
net1.params(f).value=net2.layers{8}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

%-------------------------covn 3--------------------
net1.addLayer('conv3_1', dagnn.Conv('size', [3,3,128,256],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'pool_2', 'conv_31', {'conv31_f', 'conv31_b'});
net1.addLayer('relu3_1', dagnn.ReLU(),'conv_31','relu_31');

f = net1.getParamIndex('conv31_f') ;
net1.params(f).value=net2.layers{11}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv31_b') ;
net1.params(f).value=net2.layers{11}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

net1.addLayer('conv3_2', dagnn.Conv('size', [3,3,256,256],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_31', 'conv_32', {'conv32_f', 'conv32_b'});
net1.addLayer('relu3_2', dagnn.ReLU(),'conv_32','relu_32');

f = net1.getParamIndex('conv32_f') ;
net1.params(f).value=net2.layers{13}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv32_b') ;
net1.params(f).value=net2.layers{13}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

net1.addLayer('conv3_3', dagnn.Conv('size', [3,3,256,256],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_32', 'conv_33', {'conv33_f', 'conv33_b'});
net1.addLayer('relu3_3', dagnn.ReLU(),'conv_33','relu_33');

f = net1.getParamIndex('conv33_f') ;
net1.params(f).value=net2.layers{15}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv33_b') ;
net1.params(f).value=net2.layers{15}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

%-------------------------covn 4--------------------
net1.addLayer('conv4_1', dagnn.Conv('size', [3,3,256,512],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_33', 'conv_41', {'conv41_f', 'conv41_b'});
net1.addLayer('relu4_1', dagnn.ReLU(),'conv_41','relu_41');

f = net1.getParamIndex('conv41_f') ;
net1.params(f).value=net2.layers{18}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv41_b') ;
net1.params(f).value=net2.layers{18}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

net1.addLayer('conv4_2', dagnn.Conv('size', [3,3,512,512],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_41', 'conv_42', {'conv42_f', 'conv42_b'});
net1.addLayer('relu4_2', dagnn.ReLU(),'conv_42','relu_42');

f = net1.getParamIndex('conv42_f') ;
net1.params(f).value=net2.layers{20}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv42_b') ;
net1.params(f).value=net2.layers{20}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

net1.addLayer('conv4_3', dagnn.Conv('size', [3,3,512,512],...
    'hasBias', true, 'pad',...
    [1,1,1,1], 'stride', [1,1]), 'relu_42', 'conv_43', {'conv43_f', 'conv43_b'});
net1.addLayer('relu4_3', dagnn.ReLU(),'conv_43','relu_43');

f = net1.getParamIndex('conv43_f') ;
net1.params(f).value=net2.layers{22}.weights{1};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1e3;

f = net1.getParamIndex('conv43_b') ;
net1.params(f).value=net2.layers{22}.weights{2};
net1.params(f).learningRate=0;
net1.params(f).weightDecay=1;

% net1.addLayer('pool2_show', dagnn.Crop('crop',[1,1]),...
%     {'pool_2','pool_2'}, 'pool2show');
% net1.addLayer('relu33_show', dagnn.Crop('crop',[1,1]),...
%     {'relu_33','relu_33'}, 'relu33_show');

net1.move('gpu');
clear net2;
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tflow vl_imreadjpeg ;
end

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end
end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function export_mri( BstMriFile, OutputMriFile )
% EXPORT_MRI: Export a MRI to one of the supported file formats.
%
% USAGE:  export_mri( BstMriFile, OutputMriFile=[ask] )
%         export_mri( sMri,       OutputMriFile=[ask] )
% INPUT:
%     - BstMriFile    : Full path to input Brainstorm MRI file to be exported
%     - sMri          : Brainstorm MRI structure
%     - OutputMriFile : Full path to target file (extension will determine the format)

% @=============================================================================
% 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, 2008-2016

% ===== PARSE INPUTS =====
if (nargin < 1) || isempty(BstMriFile)
    error('Brainstorm:InvalidCall', 'Invalid use of export_mri()');
end
if (nargin < 2)
    OutputMriFile = [];
end

% ===== LOAD MRI FILE =====
% Show progress bar
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'Export MRI', 'Loading input file');
end
% Load MRI
if ischar(BstMriFile)
    sMri = in_mri_bst(BstMriFile);
else
    sMri = BstMriFile;
    BstMriFile = [];
end
if ~isProgress
    bst_progress('stop');
end

% ===== SELECT OUTPUT FILE =====
if isempty(OutputMriFile)
    % Get default directories and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    % Export extension
    if isempty(DefaultFormats.MriOut)
        DefaultFormats.MriOut = 'Analyze';
    end
    switch DefaultFormats.MriOut
        case 'GIS',     ExportExt = '.ima';
        case 'Analyze', ExportExt = '.img';
        case 'CTF',     ExportExt = '.mri';
        case 'Nifti1',  ExportExt = '.nii';
        case 'FT-MRI',  ExportExt = '.mat';
        otherwise,      ExportExt = '.nii';
    end
    % Build default output filename
    if ~isempty(BstMriFile)
        [BstPath, BstBase, BstExt] = bst_fileparts(BstMriFile);
    else
        BstBase = 'export';
    end
    DefaultOutputFile = bst_fullfile(LastUsedDirs.ExportAnat, [BstBase, ExportExt]);
    DefaultOutputFile = strrep(DefaultOutputFile, '_subjectimage', '');
    DefaultOutputFile = strrep(DefaultOutputFile, 'subjectimage_', '');
    % Put MRI file
    [OutputMriFile, FileFormat, FileFilter] = java_getfile( 'save', ...
        'Export MRI...', ...             % Window title
        DefaultOutputFile, ...       % Default directory
        'single', 'files', ...           % Selection mode
        bst_get('FileFilters', 'mriout'), ...
        DefaultFormats.MriOut);
    % If no file was selected: exit
    if isempty(OutputMriFile)
        return
    end
    % Save new default export path
    LastUsedDirs.ExportAnat = bst_fileparts(OutputMriFile);
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default export format
    DefaultFormats.MriOut = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
end

% ===== SAVE MRI =====
[OutputPath, OutputBase, OutputExt] = bst_fileparts(OutputMriFile);
% Show progress bar
if ~isProgress
    bst_progress('start', 'Export MRI', ['Export MRI to file "' [OutputBase, OutputExt] '"...']);
end
% Switch between file formats
switch lower(OutputExt)
    case '.ima'
        out_mri_gis(sMri, OutputMriFile);
    case {'.img', '.nii'}
        out_mri_nii(sMri, OutputMriFile);
    case '.mri'
        out_mri_ctf(sMri, OutputMriFile);
    case '.mat'
        ftMri = out_fieldtrip_mri(sMri);
        bst_save(OutputMriFile, ftMri, 'v7');
    otherwise
        error(['Unsupported file extension : "' OutputExt '"']);
end
% Hide progress bar
if ~isProgress
    bst_progress('stop');
end

author: brainstorm-tools


text:
function [inchi, metList] = sdf2inchi(sdfFileName, options)
% Converts metabolite structures in an SDF to a cell array of InChI
% strings with OpenBabel.
%
% USAGE:
%
%    [inchi, metList] = sdf2inchi(sdfFileName, options)
%
% INPUT:
%    sdfFileName:    Path to SDF file.
%
% OPTIONAL INPUTS:
%    options:        Write options for InChI strings. See InChI documentation
%                    for details. If no options are specified the function will
%                    output standard InChI.
%
% OUTPUTS:
%    inchi:          Cell array of InChI strings for metabolites in the SDF file.
%    metList:        Cell array of metabolite identifiers (first line of each
%                    molfile in SDF). Will be empty unless write option `t` is
%                    used (i.e., options >= '-xt').
%
% .. Author: - Hulda SH, Nov. 2012

if ~strcmp(sdfFileName(end-3:end),'.sdf') % Check inputs
    sdfFileName = [sdfFileName '.sdf'];
end

if ~exist('options','var')
    options = [];
end
if ~isempty(options)
   options = [' ' strtrim(options)];
end

babelcmd = 'babel';

% Convert to InChI with OpenBabel
[success,result] = system([babelcmd ' ' sdfFileName ' -oinchi' options]);

% Parse output from OpenBabel
if success == 0
    result = regexp(result,'InChI=[^\n]*\n','match');
    result = result';
    result = strtrim(result);

    [inchi,metList] = strtok(result);
    inchi = strtrim(inchi);
    metList = strtrim(metList);
else
    [success,result] = system([babelcmd ' ' sdfFileName ' -oinchi' options])
    fprintf('%s\n','If you get a ''not found'' message from the call to Babel, make sure that Matlab''s LD_LIBRARY_PATH is edited to include correct system libraries. See initVonBertylanffy')
    error('Conversion to InChI not successful. Make sure OpenBabel is installed correctly.\n')
end

if size(inchi,2) > size(inchi,1)
    inchi = inchi';
end
if size(metList,2) > size(metList,1)
    metList = metList';
end

author: opencobra

repo: cobratoolbox


text:
function ExecFlag=DoSpinWatchExec(handles)

global VObj;
global VMag;
global VVar;
global VSeq;
global VCtl;
global VSig;
global VCoi;

%preserve VObj VMag & VCoi
VTmpObj=VObj;
VTmpMag=VMag;
VTmpCoi=VCoi;

handles.Simuh=guidata(handles.Simuh.SimuPanel_figure);
DoDisableButton([],[],handles.Simuh);
%% Do spin execution
try
    % Read tab parameters
    fieldname=fieldnames(handles.Attrh1);
    for i=1:length(fieldname)/2
        try
            eval(['SW.' fieldname{i*2} '=[' get(handles.Attrh1.(fieldname{i*2}),'String') '];']);
        catch me
            TAttr=get(handles.Attrh1.(fieldname{i*2}),'String');
            eval(['SW.' fieldname{i*2} '=''' TAttr{get(handles.Attrh1.(fieldname{i*2}),'Value')}  ''';']);
        end
    end

    % Prescan config
    DoPreScan(handles.Simuh);

    % Create SpinWatcher VOsw, VMsw & VCsw
    VOsw=VObj;
    VOsw.Gyro=SW.Gyro;
    VOsw.ChemShift=SW.ChemShift;
    VOsw.SpinNum=SW.SpinPerVoxel;
    VOsw.TypeNum=SW.TypeNum;
    VOsw.XDim=1;
    VOsw.YDim=1;
    VOsw.ZDim=1;
    VOsw=rmfield(VOsw, 'Rho');
    VOsw=rmfield(VOsw, 'T1');
    VOsw=rmfield(VOsw, 'T2');
    VOsw=rmfield(VOsw, 'T2Star');
    VOsw=rmfield(VOsw, 'Type');
    for i=1:VOsw.TypeNum
        VOsw.Rho(:,:,:,i)=SW.Rho(i);
        VOsw.T1(:,:,:,i)=SW.T1(i);
        VOsw.T2(:,:,:,i)=SW.T2(i);
        VOsw.T2Star(:,:,:,i)=SW.T2Star(i);
    end

    if isfield(VCtl, 'MT_Flag')
        if strcmp(VCtl.MT_Flag, 'on')
            K = VOsw.K;
            VOsw=rmfield(VOsw, 'K');
            VOsw.K=K(SW.LocY,SW.LocX,SW.LocZ,:);
        end
    end

    if isfield(VCtl, 'ME_Flag')
        if strcmp(VCtl.ME_Flag, 'on')
            K = VOsw.K;
            VOsw=rmfield(VOsw, 'K');
            VOsw.K=K(SW.LocY,SW.LocX,SW.LocZ,:);
        end
    end

    if isfield(VCtl, 'CEST_Flag')
        if strcmp(VCtl.CEST_Flag, 'on')
            K = VOsw.K;
            VOsw=rmfield(VOsw, 'K');
            VOsw.K=K(SW.LocY,SW.LocX,SW.LocZ,:);
        end
    end

    if isfield(VCtl, 'GM_Flag')
        if strcmp(VCtl.GM_Flag, 'on')
            K = VOsw.K;
            VOsw=rmfield(VOsw, 'K');
            VOsw.K=K(SW.LocY,SW.LocX,SW.LocZ,:);
        end
    end

    VOsw=rmfield(VOsw, 'Mx');
    VOsw=rmfield(VOsw, 'My');
    VOsw=rmfield(VOsw, 'Mz');
    VOsw.Mx= zeros(1,1,1,VOsw.SpinNum,VOsw.TypeNum);
    VOsw.My= zeros(1,1,1,VOsw.SpinNum,VOsw.TypeNum);
    for i=1:VOsw.TypeNum
        VOsw.Mz(:,:,:,:,i)= ones(1,1,1,VOsw.SpinNum,1).*(VOsw.Rho(i)./VOsw.SpinNum);
    end

    % Gradient Grid
    VMsw=VMag;
    VMsw=rmfield(VMsw, 'Gxgrid');
    VMsw=rmfield(VMsw, 'Gygrid');
    VMsw=rmfield(VMsw, 'Gzgrid');
    VMsw=rmfield(VMsw, 'dB0');
    VMsw=rmfield(VMsw, 'dWRnd');
    VMsw=rmfield(VMsw, 'FRange');
    VMsw.Gxgrid=VMag.Gxgrid(SW.LocY,SW.LocX,SW.LocZ);
    VMsw.Gygrid=VMag.Gygrid(SW.LocY,SW.LocX,SW.LocZ);
    VMsw.Gzgrid=VMag.Gzgrid(SW.LocY,SW.LocX,SW.LocZ);
    if VMag.dB0==0
        VMsw.dB0=0;
    else
        VMsw.dB0=VMag.dB0(SW.LocY,SW.LocX,SW.LocZ);
    end
    if VOsw.SpinNum > 1
        InddWRnd=linspace(0.01,0.99,VOsw.SpinNum);
        for j=1:VOsw.TypeNum
            for i=1:VOsw.SpinNum
                VMsw.dWRnd(:,:,:,i,j)=(1./VOsw.T2Star(:,:,:,j)-1./VOsw.T2(:,:,:,j)).*tan(pi.*(InddWRnd(i)-1/2));
                % need large number of spins for stimulating T2* effect,
                % insufficient number of spins may cause in-accurate simulation
            end
        end
    else
        for j=1:VOsw.TypeNum
            for i=1:VOsw.SpinNum
                VMsw.dWRnd(:,:,:,i,j)= 0;
            end
        end
    end
    VMsw.FRange=1;

    % Coil
    VCsw=VCoi;
    if isfield(handles.Simuh,'CoilTxXMLFile')
        VCsw.TxCoilmg=VCoi.TxCoilmg(SW.LocY - max(1, str2num(get(handles.Simuh.AP1_text,'String'))) + 1, ...   % Align postion with respect to FOV
                                    SW.LocX - max(1, str2num(get(handles.Simuh.LR1_text,'String'))) + 1, ...
                                    SW.LocZ - max(1, str2num(get(handles.Simuh.SI1_text,'String'))) + 1, :);
        VCsw.TxCoilpe=VCoi.TxCoilpe(SW.LocY - max(1, str2num(get(handles.Simuh.AP1_text,'String'))) + 1, ...   % Align postion with respect to FOV
                                    SW.LocX - max(1, str2num(get(handles.Simuh.LR1_text,'String'))) + 1, ...
                                    SW.LocZ - max(1, str2num(get(handles.Simuh.SI1_text,'String'))) + 1, :);
    else
        VCsw.TxCoilmg=VCoi.TxCoilmg(SW.LocY,SW.LocX,SW.LocZ,:);
        VCsw.TxCoilpe=VCoi.TxCoilpe(SW.LocY,SW.LocX,SW.LocZ,:);
    end
    %% Spin execution
    VObj=VOsw;
    VMag=VMsw;
    VCoi=VCsw;

    % Generate Pulse line
    DoPulseGen(handles.Simuh);

    % Initialize Mx, My, Mz, Muts in VSig
    VSig.Mx= zeros(1,1,1,VOsw.SpinNum,VOsw.TypeNum, max(size(VSeq.utsLine)));
    VSig.My= zeros(1,1,1,VOsw.SpinNum,VOsw.TypeNum, max(size(VSeq.utsLine)));
    VSig.Mz= zeros(1,1,1,VOsw.SpinNum,VOsw.TypeNum, max(size(VSeq.utsLine)));
    VSig.Muts=zeros(max(size(VSeq.utsLine)),1);

    VSig.Mx(:,:,:,:,:,1)=VOsw.Mx;
    VSig.My(:,:,:,:,:,1)=VOsw.My;
    VSig.Mz(:,:,:,:,:,1)=VOsw.Mz;
    VSig.Muts(1)=0;

    % Simulation Process
    try
        VCtl.CS=double(VObj.ChemShift*VCtl.B0);
        VCtl.RunMode=int32(1); % Spin scan
        VCtl.MaxThreadNum=int32(handles.Simuh.CPUInfo.NumThreads);
        VCtl.ActiveThreadNum=int32(0);
        DoDataTypeConv(handles.Simuh);

        if isfield(VCtl, 'MT_Flag')
            if strcmp(VCtl.MT_Flag, 'on')
                VCtl.GPUIndex=int32(0);
                DoMTScanAtGPU;  % beta MT kernel, only support two-pool MT
            else
                DoScanAtCPU;  % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed
            end
        elseif isfield(VCtl, 'ME_Flag')
            if strcmp(VCtl.ME_Flag, 'on')
                VCtl.GPUIndex=int32(0);
                DoMEScanAtGPU;  % beta ME kernel
            else
                DoScanAtCPU;  % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed
            end
        elseif isfield(VCtl, 'CEST_Flag')
            if strcmp(VCtl.CEST_Flag, 'on')
                VCtl.GPUIndex=int32(0);
                DoCESTScanAtGPU;  % beta CEST kernel
            else
                DoScanAtCPU;  % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed
            end
        elseif isfield(VCtl, 'GM_Flag')
            if strcmp(VCtl.GM_Flag, 'on')
                VCtl.GPUIndex=int32(0);
                DoGMScanAtGPU;  % beta GM kernel
            else
                DoScanAtCPU;  % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed
            end
        else
            DoScanAtCPU;  % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed
        end

    catch me
        error_msg{1,1}='ERROR!!! Spin execution process aborted.';
        error_msg{2,1}=me.message;
        errordlg(error_msg);
        ExecFlag=0;
        %recover VObj
        VObj=VTmpObj;
        VMag=VTmpMag;
        VCoi=VTmpCoi;
        return;
    end

    % Convert VSig.Mx, VSig.My, VSig.Mz
    VSig.Mx = double(VSig.Mx);
    VSig.My = double(VSig.My);
    VSig.Mz = double(VSig.Mz);

    handles.Mx=VSig.Mx.*double(VObj.SpinNum);
    handles.My=VSig.My.*double(VObj.SpinNum);
    handles.Mz=VSig.Mz.*double(VObj.SpinNum);
    handles.Muts=VSig.Muts;
    handles.MutsStep=length(VSig.Muts);

    Mx=sum(reshape(VSig.Mx,[VObj.SpinNum,VObj.TypeNum, max(size(VSeq.utsLine))]),1);
    My=sum(reshape(VSig.My,[VObj.SpinNum,VObj.TypeNum, max(size(VSeq.utsLine))]),1);
    Mz=sum(reshape(VSig.Mz,[VObj.SpinNum,VObj.TypeNum, max(size(VSeq.utsLine))]),1);
    Mx=reshape(Mx,[VObj.TypeNum max(size(VSeq.utsLine))]);
    My=reshape(My,[VObj.TypeNum max(size(VSeq.utsLine))]);
    Mz=reshape(Mz,[VObj.TypeNum max(size(VSeq.utsLine))]);
    Mx=Mx';
    My=My';
    Mz=Mz';
    handles.MxySum=sqrt(Mx.^2+My.^2);
    handles.MzSum=Mz;
    handles.MaxMxySum=max(max(handles.MxySum));
    handles.MaxMzSum=max(max(handles.MzSum));
    handles.MinMxySum=min(min(handles.MxySum));
    handles.MinMzSum=min(min(handles.MzSum));
    handles.Gzgrid=zeros(1,1,1,VObj.SpinNum);
    handles.Gygrid=zeros(1,1,1,VObj.SpinNum);
    handles.Gxgrid=zeros(1,1,1,VObj.SpinNum);
catch me
    error_msg{1,1}='ERROR!!! Spin execution process aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    ExecFlag=0;
    %recover VObj
    VObj=VTmpObj;
    VMag=VTmpMag;
    VCoi=VTmpCoi;
    return;
end

%recover VObj
VObj=VTmpObj;
VMag=VTmpMag;
VCoi=VTmpCoi;

guidata(handles.SpinWatcherPanel_figure, handles);
ExecFlag=1;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function D = sqdistance(A, B, M)
% Square Euclidean or Mahalanobis distances between all sample pairs
% A: d x n1 data matrix
% B: d x n2 data matrix
% M: d x d  Mahalanobis matrix
% D: n1 x n2 pairwise square distance matrix
% Written by Michael Chen (sth4nth@gmail.com). July 2009.

if nargin == 1
    A = bsxfun(@minus,A,mean(A,2));
    S = full(sum(A.^2,1));
    D = full((-2)*(A'*A));
    D = bsxfun(@plus,D,S);
    D = bsxfun(@plus,D,S');
elseif nargin == 2
    assert(size(A,1)==size(B,1));

    m = (sum(A,2)+sum(B,2))/(size(A,2)+size(B,2));
    A = bsxfun(@minus,A,m);
    B = bsxfun(@minus,B,m);
    D = full((-2)*(A'*B));
    D = bsxfun(@plus,D,full(sum(B.^2,1)));
    D = bsxfun(@plus,D,full(sum(A.^2,1)'));
elseif nargin == 3
    assert(size(A,1)==size(B,1));

    m = (sum(A,2)+sum(B,2))/(size(A,2)+size(B,2));
    A = bsxfun(@minus,A,m);
    B = bsxfun(@minus,B,m);
    D = full((-2)*(A'*M*B));
    D = bsxfun(@plus,D,full(sum(B.*(M*B),1)));
    D = bsxfun(@plus,D,full(sum(A.*(M*A),1)'));
end

author: ZJULearning

repo: MatlabFunc


text:
function dlap_io_test01 ( )

%*****************************************************************************80
%
%% DLAP_IO_TEST01 tests DLAP_FILE_WRITE.
%
%  Discussion:
%
%    The matrix is:
%
%      11  12   0   0  15
%      21  22   0   0   0
%       0   0  33   0  35
%       0   0   0  44   0
%      51   0  53   0  55
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2014
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  nelt = 11;

  a = [ 51.0, 12.0, 11.0, 33.0, 15.0, ...
        53.0, 55.0, 22.0, 35.0, 44.0, ...
        21.0 ];
  ia = [ 5,  1,  1,  3,  1,  5,  5,  2,  3,  4,  2 ];
  irhs = 1;
  isoln = 1;
  isym = 0;
  output_filename = 'a5by5.dlap';
  ja = [ 1,  2,  1,  3,  5,  3,  5,  2,  5,  4,  1 ];
  rhs = [ 110.0, 65.0, 274.0, 176.0, 485.0 ];
  soln = [ 1.0, 2.0, 3.0, 4.0, 5.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DLAP_IO_TEST01\n' );
  fprintf ( 1, '  DLAP_FILE_WRITE writes a matrix in SLAP Triad format\n' );
  fprintf ( 1, '  to a DLAP sparse matrix file.\n' );

  dlap_file_print ( n, nelt, isym, irhs, isoln, ia, ja, a, rhs, ...
    soln, '  The DLAP data to be written to the file.' );

  output = fopen ( output_filename, 'wt' );

  dlap_file_write ( n, nelt, isym, irhs, isoln, ia, ja, a, rhs, soln, output );

  fclose ( output );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the matrix data to "%s"\n', output_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function k = findrows(A,b)
%FINDROWS finds a row vector within a matrix.
%
%   Usage: idx = findrows(A,b)
%
%   Input parameters:
%       A       - matrix
%       b       - vector to find in A
%
%   Output parameters:
%       idx     - indices of found columns in A
%
%   FINDROWS(A,b) returns a column vector with the indices of the rows
%   in the matrix A that are identical to the row vector b.  If no rows
%   in A are identical to b, an empty vector is returned.
%
%   See also: find, findcols

%*****************************************************************************
% 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 *
%*****************************************************************************

% AUTHOR: Peter John Acklam

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

%% ===== Computation =====================================================
k = find( A(:,1)==b(1) );
for jj = 2:size(A,2)
    k = k( A(k,jj)==b(jj) );
    if isempty(k)
        return
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function [gx,dgdx,dgdp] = g_DCMwHRFext(Xt,P,ut,in)

%- initializations
idX1 = 1:in.n5(end);
nreg=length(in.n5);
nresp=length(in.r);

gxr=zeros(nresp,1);
dgdxr=zeros(length(Xt),nresp);
dgdpr=zeros(length(P),nresp);

%== compute gradients
[gxn,dgdxn,dgdpn] = g_HRF3(Xt(idX1),P,ut,in) ;
for iR = 1:numel(in.sourceRespIdx)
    respIdx = in.sourceRespIdx{iR};
    idRinX = in.r(respIdx);
    idRinPhi = in.indr(iR);
    [gxr(respIdx),dgdxr(idRinX,respIdx),dgdpr(idRinPhi,respIdx)] = g_softmax4decoding(Xt(idRinX),P(idRinPhi),ut,in);
end
%== concatenate BOLD and behavioral responses
%- state
gx = [gxn ; gxr] ;

%- Jacobian
dgdx = zeros(length(Xt),length(gx));
dgdx(idX1,1:nreg) = dgdxn ;
dgdx(:,nreg+(1:nresp)) = dgdxr ;

%- gradient wrt parameters
dgdp = zeros(length(P),nreg+nresp);
dgdp(1:size(dgdpn,1),1:nreg) = dgdpn ;
dgdp(:,nreg+(1:nresp)) = dgdpr ;

end

author: MBB-team

repo: VBA-toolbox


text:
%% figuremax
% Below is a demonstration of the features of the |figuremax| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=figuremax(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
clear
% close all
clc;

addpath(genpath('../libs'))
addpath('./myFunc')
% path_to_matconvnet = '../matconvnet';
path_to_matconvnet = '../matconvnet_std';
path_to_matconvnet = '/home/skong2/scratch/matconvnet-1.0-beta23_modifiedDagnn';
path_to_model = '../models/';

run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));
%% read matconvnet model
load('imdb_toydata_v3_from_mnist.mat');
% imdb.path_to_dataset = '/run/shm/toydata';

% set GPU
gpuId = 2; %[1, 2];
gpuDevice(gpuId);
flagSaveFig = true; % {true false} whether to store the result

saveFolder = 'main007_instSeg_v1_absEucMM';
modelName = 'softmax_net-epoch-83.mat';

saveFolder = 'main007_instSeg_v2_absEucMM_invAlpha';
modelName = 'softmax_net-epoch-83.mat';

%saveFolder = 'main007_instSeg_v3_absEucMM_adaptMM';
%modelName = 'softmax_net-epoch-92.mat';

%% modify network for testing
netMat = load( fullfile('./exp', saveFolder, modelName) );
netMat = netMat.net;
netMat = dagnn.DagNN.loadobj(netMat);

% RFinfo = netMat.getVarReceptiveFields('input');

rmLayerName = 'obj_instSeg_MM';
if ~isnan(netMat.getLayerIndex(rmLayerName))
    netMat.removeLayer(rmLayerName); % remove layer
end
rmLayerName = 'obj_instSeg';
if ~isnan(netMat.getLayerIndex(rmLayerName))
    netMat.removeLayer(rmLayerName); % remove layer
end
rmLayerName = 'obj_instSeg_reg';
if ~isnan(netMat.getLayerIndex(rmLayerName))
    netMat.removeLayer(rmLayerName); % remove layer
end
netMat.vars(netMat.layers(netMat.getLayerIndex('res7_l2norm')).outputIndexes).precious = 1;
%% configure training environment
saveFolder = [strrep(saveFolder, '/', '') '_visualization'];
if ~isdir(saveFolder) && flagSaveFig
    mkdir(saveFolder);
end
netMat.move('gpu');
% netMat.mode = 'test';
netMat.mode = 'normal';
testList = find(imdb.set==2);
rawCountsAll = {};

imgMat = [];
semanticMaskMat = [];
instanceMaskMat = [];
predInstanceMaskMat = [];
st = 0;
for i = 1:5%length(testList)
    curBatch = fullfile(imdb.path_to_dataset, imdb.imgList{testList(i)});
    datamat = load(curBatch);

%     datamat.imgMat = reshape(datamat.imgMat, [size(datamat.imgMat,1), size(datamat.imgMat,2), 1, size(datamat.imgMat,3)]);
    datamat.semanticMaskMat = reshape(datamat.semanticMaskMat, [size(datamat.semanticMaskMat,1), size(datamat.semanticMaskMat,2), 1, size(datamat.semanticMaskMat,3)]);
    datamat.instanceMaskMat = reshape(datamat.instanceMaskMat, [size(datamat.instanceMaskMat,1), size(datamat.instanceMaskMat,2), 1, size(datamat.instanceMaskMat,3)]);

    imFeed = bsxfun(@minus, datamat.imgMat, imdb.meta.meanvalue);
    inputs = {'input', gpuArray(imFeed)};
    netMat.eval(inputs) ;
    res7_l2norm = gather(netMat.vars(netMat.layers(netMat.getLayerIndex('res7_l2norm')).outputIndexes).value);

    imgMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.imgMat;
    semanticMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.semanticMaskMat;
    instanceMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = datamat.instanceMaskMat;
    predInstanceMaskMat(:,:,:,st+1:st+size(datamat.imgMat,4)) = res7_l2norm;
    fprintf('\t%d/%d\n', i, length(testList));

    st = st + size(datamat.imgMat,4);
end
%% visualize
num = 25;
M = 5;
N = 5;
h = 28;
w = 28;
panelSZ = round(64/2)*2;
stepSize = (panelSZ-h)/2;

showcaseSamples_R = imgMat(:,:,1,1:num);
showcaseSamples_R = reshape(showcaseSamples_R, [numel(showcaseSamples_R(:,:,1,1)), num]);
showcaseSamples_R = showStitch(showcaseSamples_R, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);
showcaseSamples_G = imgMat(:,:,2,1:num);
showcaseSamples_G = reshape(showcaseSamples_G, [numel(showcaseSamples_G(:,:,1,1)), num]);
showcaseSamples_G = showStitch(showcaseSamples_G, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);
showcaseSamples_B = imgMat(:,:,3,1:num);
showcaseSamples_B = reshape(showcaseSamples_B, [numel(showcaseSamples_B(:,:,1,1)), num]);
showcaseSamples_B = showStitch(showcaseSamples_B, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);
showcaseSamples = cat(3, showcaseSamples_R, showcaseSamples_G, showcaseSamples_B);

showcaseSemanticMask= semanticMaskMat(:,:,:,1:num);
showcaseSemanticMask = reshape(showcaseSemanticMask, [numel(showcaseSemanticMask(:,:,1,1)), num]);
showcaseSemanticMask = showStitch(showcaseSemanticMask, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);

showcaseInstanceMask= instanceMaskMat(:,:,:,1:num);
showcaseInstanceMask = reshape(showcaseInstanceMask, [numel(showcaseInstanceMask(:,:,1,1)), num]);
showcaseInstanceMask = showStitch(showcaseInstanceMask, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);

predInstanceMask_R = predInstanceMaskMat(:,:,1,1:num);
predInstanceMask_R = reshape(predInstanceMask_R, [numel(predInstanceMask_R(:,:,1,1)), num]);
predInstanceMask_R = showStitch(predInstanceMask_R, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);
predInstanceMask_G = predInstanceMaskMat(:,:,2,1:num);
predInstanceMask_G = reshape(predInstanceMask_G, [numel(predInstanceMask_G(:,:,1,1)), num]);
predInstanceMask_G = showStitch(predInstanceMask_G, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);
predInstanceMask_B = predInstanceMaskMat(:,:,3,1:num);
predInstanceMask_B = reshape(predInstanceMask_B, [numel(predInstanceMask_B(:,:,1,1)), num]);
predInstanceMask_B = showStitch(predInstanceMask_B, [panelSZ panelSZ], M, N, 'whitelines', 'linewidth', 5);
predInstanceMask = cat(3, predInstanceMask_R, predInstanceMask_G, predInstanceMask_B);

predInstanceMask = predInstanceMask - min(predInstanceMask(:));
predInstanceMask = predInstanceMask ./ max(predInstanceMask(:));

imgFig2 = figure(2);
subWindowH = 2;
subWindowW = 2;
set(imgFig2, 'Position', [100 100 1400 900]) % [1 1 width height]
windowID = 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseSamples); axis off image; title('random samples'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseSemanticMask); axis off image; title('showcase SemanticMask'); caxis([0, 11]); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(showcaseInstanceMask); axis off image; title('showcase instanceMask'); colorbar; windowID = windowID + 1;
subplot(subWindowH, subWindowW, windowID); imagesc(predInstanceMask); axis off image; title('predInstanceMask'); colorbar; windowID = windowID + 1;

if flagSaveFig
    showcaseSemanticMask = uint8(showcaseSemanticMask * 255/11);
    showcaseInstanceMask = showcaseInstanceMask*255/5;

    imwrite(showcaseSamples, fullfile(saveFolder, sprintf('%04d_showcaseSamples.bmp',i)));
    imwrite(showcaseSemanticMask, fullfile(saveFolder, sprintf('%04d_showcaseSemanticMask.bmp',i)));
    imwrite(showcaseInstanceMask, fullfile(saveFolder, sprintf('%04d_showcaseInstanceMask.bmp',i)));
    imwrite(predInstanceMask, fullfile(saveFolder, sprintf('%04d_predInstanceMask.bmp',i)));

    export_fig( sprintf('%s/%04d_visualization.jpg', saveFolder, i) );
    save(sprintf('%s/results.mat', saveFolder), 'imgMat', 'instanceMaskMat', 'semanticMaskMat', 'predInstanceMaskMat');
end

%% leaving blank

author: aimerykong


text:
function [state_red,inv_state_red,n_g_red] = krylov_reduction(g0,g1,n_v,n_g,n_arnoldi,varargin)
% Reduce the dimensionality of state space using Krylov subspace method
%
% by SeHyoun Ahn, Sept 2016
%
% REFERENCE: Ahn, SeHyoun, Greg Kaplan, Benjamin Moll, Thomas Winberry, and
%    Christian Wolf. "When Inequality Matters for Macro and Macro Matters
%    for Inequality."
%
% PARAMETERS:
%    g0 = LHS matrix, only used to check it satisfies require form
%    g1 = Dynamics matrix
%    n_v = number of jump variables
%    n_g = number of state variables to be reduced
%    n_Z = number of state variables not being reduced
%    n_arnoldi = dimension of Krylov subspace
%    varargin(1) = other observables
%    varargin(2) = updating parameter for B_gv
%    varargin(3) = number of variables to be NOT reduced; it is assumed
%                  that the last varargin(3) variables are not being
%                  reduced, so reorder as necessary
%
% OUTPUTS:
%    state_red = tranformation to get from full grid to reduced states
%    inv_state_red = inverse transform
%    n_g_red = number of state variables after reduction
%
% SYNTAX:
% [state_red,inv_state_red,n_g_red] = krylov_reduction(g0,g1,n_v,n_g,n_arnoldi,varargin)

if nargin == 8
    observable = varargin{1};
    F = varargin{2};
    n_Z = varargin{3};
elseif nargin == 7
    observable = varargin{1};
    F = varargin{2};
    n_Z = 0;
elseif nargin == 6
    observable = varargin{1};
    F = [];
    n_Z = 0;
else
    observable = [];
    F = [];
    n_Z = 0;
end
n_p = size(g1,1) - n_v - n_g;
n_total = n_v + n_g;
n_g = n_g - n_Z;

%% Check to make sure that g0 is an identity matrix
%     This reducetion step assumes that g0 has been cleaned, so if
%     this part fails, make sure to run <clean_g.m>
if (max(max(abs(g0(1:n_total,1:n_total)-speye(n_total))))~=0)
    error('Make sure that g0 is normalized.');
end

% Slice Dynamics Equation into Different Parts
B_pv = -g1(n_total+1:n_total+n_p,n_total+1:n_total+n_p)\g1(n_total+1:n_total+n_p,1:n_v);
B_pg = -g1(n_total+1:n_total+n_p,n_total+1:n_total+n_p)\g1(n_total+1:n_total+n_p,n_v+1:n_v+n_g);
B_pZ = -g1(n_total+1:n_total+n_p,n_total+1:n_total+n_p)\g1(n_total+1:n_total+n_p,n_v+n_g+1:n_v+n_g+n_Z);
B_gg = g1(n_v+1:n_v+n_g,n_v+1:n_v+n_g);
B_gv = g1(n_v+1:n_v+n_g,1:n_v);
B_gp = g1(n_v+1:n_v+n_g,n_total+1:n_total+n_p);
% B_vp = g1(1:n_v,n_v+n_g+1:n_v+n_g+n_p);

%% Drop redundant Directions in B_pg
% In theory, since we do deflated block arnoldi, this step
%    might be redundant. Will be tested later.
obs = [B_pg;observable];
[~,d0,V_g] = svd(full(obs),'econ');
aux = diag(d0);
aux = aux/aux(1);
n_Bpg = sum(aux>10*eps);
V_g = bsxfun(@times,V_g(:,1:n_Bpg),aux(1:n_Bpg)');

% Compute Krylov Subspace
if isa(F,'function_handle')
    A = @(x) F(B_gv'*x)+B_gg'*x + B_pg'*(B_gp'*x);
elseif F
    A = @(x) F'*(B_gv'*x)+B_gg'*x + B_pg'*(B_gp'*x);
else
    A = @(x) B_gg'*x + B_pg'*(B_gp'*x);
end
[V_g,~] = deflated_block_arnoldi(A,V_g,n_arnoldi);
n_g_red = size(V_g,2);

% Build State-Space Reduction transform
state_red = sparse(n_v+n_g_red,n_total+n_p);
state_red(1:n_v,1:n_v) = speye(n_v);
state_red(n_v+1:n_v+n_g_red,n_v+1:n_v+n_g) = V_g';
state_red(n_v+n_g_red+1:n_v+n_g_red+n_Z,n_v+n_g+1:n_v+n_g+n_Z) = speye(n_Z);
%state_red(:,n_total+1:n_total+n_p) = 0;

% Build inverse transform
inv_state_red = sparse(n_total+n_p,n_v+n_g_red);
inv_state_red(1:n_v,1:n_v) = speye(n_v);
inv_state_red(n_v+1:n_v+n_g,n_v+1:n_g_red+n_v) = V_g;
inv_state_red(n_total+1:n_total+n_p,1:n_v) = B_pv;
inv_state_red(n_total+1:n_total+n_p,n_v+1:n_v+n_g_red) = B_pg*V_g;
inv_state_red(n_total+1:n_total+n_p,n_v+n_g_red+1:n_v+n_g_red+n_Z) = B_pZ;
inv_state_red(n_v+n_g+1:n_total,n_v+n_g_red+1:n_v+n_g_red+n_Z) = speye(n_Z);

n_g_red = n_g_red + n_Z;

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
if ~isempty(strfind(fsub,'.mat'))
    if exist(fsub,'file')
        save(fsub,'X','-append')
    else
        save(fsub,'X')
    end
elseif ~isempty(strfind(fsub,'.txt'))
    dlmwrite(fsub,X);
else % no extension - assumed to be an SPM file
    error('For now, we can only deal with .mat and .txt files')
end

author: OHBA-analysis

repo: HMM-MAR


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes the angular displacements of a simple gravity pendulum
%
% PURPOSE: to illustrate no change in oscillatory amplitude under
% variations in the radius of the pendulum bob.
%
% Author: N.A. Battista
% Date: 12/13/2019
% Institution: The College of New Jersey
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function simple_pendulum()

%
% Physical Parameters
%
m = 100;    % Mass of Pendulum bob
g=9.80665;  % Gravitational Acceleration;
L=0.2;      % Length of Pendulum

%
% Vector of Radii
%
rVec = 100*[0.005 0.010 0.015 0.02 0.025];

%
% Temporal Parameters
%
t=0;        % initial time
tFinal = 8; % final-time
dt = 1e-3;  % time-step
ct=1;       % counter
tVec(ct)=t; % time vector storage

%
% Moment of Inertia
%
IVec = m*(0.5*rVec.^2+L^2);

%
% Initial Angular Displacement
%
ang = -(pi/2-pi/5)*ones(1,5);
phi = zeros( size( ang ) );

%
% Perform the time-stepping, e.g., solve the ODEs
%
while t<=tFinal

    % increment the storage counter
    ct = ct+1;

    % increment time
    t = t+dt;

    % save current time
    tVec(ct) = t;

    % solve 1st order ODE for angular displacement
    ang(ct,:) = ang(ct-1,:) + dt * ( phi(ct-1,:) );

    % solve 1st order ODE for angular velocity
    phi(ct,:) = phi(ct-1,:) + dt * ( -m*g*L ./IVec .* sin( ang(ct-1,:) ) );

end

%
% Defining Colors for Plotting Data
%
color2 = [0.635 0.078 0.184];  % dark red
color5 = [0.9 0.525 0.098];    % orange
color6 = [0.929 0.840 0.1250]; % dark yellow
color7 = [0.7 0.7 0];          % 'ugly' green

%
% Note discrepancy is a bit off due to small angle approximation:
% sin(theta) ~ theta gives below period of oscillation
%
T1 = 2*pi*sqrt( IVec(1) / (m*g*L) );
%
% Actually computed value of period is:
T1 = 1.93;

%
% Make Figure of Angular Displacement vs. Non-dimensional Time (# of oscillations of 'r' case)
%
lw=6;
fs=18;
plot(tVec/T1,ang(:,1),'k-','LineWidth',lw); hold on;
plot(tVec/T1,ang(:,2),'b-','LineWidth',lw,'Color',color2); hold on;
plot(tVec/T1,ang(:,3),'g-','LineWidth',lw,'Color',color5); hold on;
plot(tVec/T1,ang(:,4),'r-','LineWidth',lw,'Color',color6); hold on;
plot(tVec/T1,ang(:,5),'c-','LineWidth',lw,'Color',color7); hold on;
xlabel('Non-Dimensional Time (# periods of r-case)');
ylabel('Angular Displacement (Radians)');
axis([0 4*T1 -1 1.25]);
leg=legend('r','2r','3r','4r','5r','Orientation','horizontal');
set(gca,'FontSize',fs);
set(leg,'FontSize',fs);

author: nickabattista


text:
function [mu, inv_sigma] = fit_gaussian(x,class)

N = max(class);
mu = zeros(N,size(x,2));
inv_sigma = zeros(size(x,2),size(x,2),N);

for i=1:N,
    mu(i,:) = mean(x(class==i,:));
    inv_sigma(:,:,i) = inv(cov(x(class==i,:)));
end

author: csn-le


text:
function varargout = uplus(varargin)
%UPLUS   Unary plus for a DISKFUN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = uplus@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function kern = rbfperiodic2KernExpandParam(kern, params)

% RBFPERIODIC2KERNEXPANDPARAM Create kernel structure from RBFPERIODIC2 kernel's parameters.
% FORMAT
% DESC returns a RBF periodic covariance with variying period kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : rbfperiodic2KernParamInit, rbfperiodic2KernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2009
%
% MODIFICATIONS : Andreas C. Damianou, 2011
%
% MODIFICATIONS : Michalis K. Titsias, 2011

% KERN

kern.inverseWidth = params(1);
kern.variance = params(2);
kern.factor = params(3);
kern.period = 2*pi/kern.factor;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function xdot = Equation(t,x)
% Set the input values in order to pass onto ode45
%
n = length(x) ;
thita = x(1) ;          % Angular Position of the bead
dthita = x(2) ;         % Angular Velocity of the bead
g = x(3) ;              % Acceleration due to gravity
M = x(4) ;              % Mass of the bead
R = x(5) ;              % Radius of the hoop
V = x(6) ;              % Frictional coefficient of the bead on the hoop
w0 = x(7) ;             % Frequency of rotation of the bead
%
xdot=zeros(n,1);
xdot(1) = dthita;
xdot(2) = -sin(thita)*(g/R-w0^2*cos(thita))-V/M*dthita ;

author: Sable

repo: mcbench-benchmarks


text:
function fem2d_pack_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests GRID_TEST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  Test the grid routines.\n' );

  grid_test ( 'Q4' );

  grid_test ( 'Q8' );

  grid_test ( 'Q9' );

  grid_test ( 'Q12' );

  grid_test ( 'Q16' );

  grid_test ( 'QL' );

  grid_test ( 'T3' );

  grid_test ( 'T4' );

  grid_test ( 'T6' );

  grid_test ( 'T10' );

  return
end


text:
%
function [map, pose] = AddAKeyScan(map,...
                                   gridMap,...
                                   scan,...
                                   pose,...
                                   hits,...
                                   pixelSize,...
                                   bruteResolution,...
                                   tmax,...
                                   rmax)

% First, evaluate the pose & hits, make sure that there is no large error
lastKeyPose = map.keyscans(end).pose;
dp = DiffPose(lastKeyPose, pose);
if abs(dp(1)) > 0.5 || abs(dp(2)) > 0.5 || abs(dp(3)) > pi
    disp('Oh no no no nobody but you : So Large Error!');
    pose = lastKeyPose;
end

% ToDo: refine the relative pose, estimate the pose's covariance.
% And push them into the map.connections, which will be neede when we close
% a loop (pose graph optimization).

scan_w = Transform(scan, pose);
newPoints = scan_w(hits>1.1, :);
%
if isempty(newPoints)
    return;
end

% keyscans
k = length(map.keyscans);
map.keyscans(k+1).pose = pose;
map.keyscans(k+1).iBegin = size(map.points, 1) + 1;
map.keyscans(k+1).iEnd = size(map.points, 1) + size(newPoints, 1);
map.keyscans(k+1).loopTried = false;
map.keyscans(k+1).loopClosed = false;

% points
map.points = [map.points; newPoints];

% connections
% ToDo: Estimate the relative pose and covariance, they will be useful
% when we close a loop (pose graph optimization).
c = length(map.connections);
map.connections(c+1).keyIdPair = [k, k+1];
map.connections(c+1).relativePose = zeros(3,1);
map.connections(c+1).covariance = zeros(3);

author: meyiao

repo: LaserSLAM

save_path: github-repos/MATLAB/meyiao-LaserSLAM


text:
classdef AD9361TRx
    %% properties applicable to simulation-based study
    properties (Constant)
        % WLAN Toolbox functions operate over 1x-rate data.
        % Oversample signal prior to passing it through AD9361 Rx Simulink model.
        ovx = [2 1] % ovx factor = ovx(1)/ovx(2)
    end

    properties
        % AGC mode (applicable only when SIM_MODE is set to 1 or 3)
        % 0 - manual
        % 1 - fast attack
        % 2 - slow attack
        % 3 - fast during preamble, manual during data payload
        AGC_MODE = 3
        % Log output of ADC block? 1/0 - Y/N
        % The size of log-file will be significantly larger if set to 1
        % since ADC runs at a higher rate than the sample clock.
        % Applicable when SIM_STUDY is set to true.
        LOG_ADC_OUTPUT = 0
        % Save log-data to .mat file? 1/0 - Y/N
        % Applicable when SIM_STUDY is set to true.
        SAVE_LOG_DATA = 1
    end

    properties (Access = private)
        num_sims
        sig_filename
        testwaveform_fileloc
        folder_name
        model
        block
        rx_block
        logged_blocks
        error_message
        log_data
        log_data_indices
    end

    properties
        rxNonHT_2x
        input_waveform
    end

    %% properties applicable to physical radio-based study
    properties
        RXBufferSize = 2^13;
        TransmitRepeat = false
        TimeOut = 2
    end

    methods
        function obj = AD9361TRx(sim_settings, varargin)%ad9361_settings, agc_settings)
            if ~isempty(varargin{1})
                ad9361_settings = varargin{1}{1};
                agc_settings = varargin{1}{2};
            end
            % local variables
            fs = obj.fs;
            M = obj.ovx(1);
            N = obj.ovx(2);
            SampleTime_rx = obj.SampleTime_rx;
            obj.folder_name = 'ad9361_rx_input_files';

            % if (sim_settings.SIM_STUDY)
                %% Simulation-based study
                if ((sim_settings.SIM_MODE == 0) || (sim_settings.SIM_MODE == 2))
                    % Simulation mode #0 or #2 - no simulink-model-in-loop
                    obj.input_waveform = obj.rxWaveform;
                else
                    % Simulation mode #1 or #3 - simulink-model-in-loop
                    % generate all combinations of AGC settings
                    [local_agc_settings, obj.num_sims] = obj.agc_settings_all_combos(agc_settings);

                    % Enable logging for desired signals
                    obj.model = 'ad9361_rx_wlan_testbench';
                    obj.block = 'ad9361_rx';
                    load_system(obj.model); % load model into memory
                    st_param = get_param(obj.model, 'ModelWorkspace');
                    assignin(st_param,'SampleTime_rx', SampleTime_rx); % assign value to sample-time
                    obj = obj.gen_testWaveform_file(); % set Baseband file name
                    obj.rx_block = strcat(obj.model, '/', obj.block);

                    obj.logged_blocks = {['ADC ' newline 'Overload ' newline 'Detector'],...
                        ['LMT ' newline 'Peak ' newline 'Detector'],...
                        'DDC_Filters_RX', 'In', 'Average Power Meter', 'AGC'};
                    if (ad9361_settings.LOG_ADC_OUTPUT)
                        obj.logged_blocks{end+1} = 'ADC_RX';
                    end
                    num_logged_signals = obj.enable_signal_logging();
                    save_system(obj.model);

                    % Tune AGC settings and run the AD9361 simulink model, ...
                    sim_time = ceil(1e4*length(obj.rxNonHT_2x)/(M*fs/N))/1e4;
                    WLAN_ad9361_simin = Simulink.SimulationInput.empty(obj.num_sims, 0);
                    for idx = 1:obj.num_sims
                        WLAN_ad9361_simin(idx) = Simulink.SimulationInput(obj.model);
                        WLAN_ad9361_simin(idx) = WLAN_ad9361_simin(idx).setModelParameter('StopTime',sprintf('%f',sim_time));
                        WLAN_ad9361_simin(idx) = obj.tune_agc_settings(ad9361_settings.AGC_MODE, idx, WLAN_ad9361_simin(idx), local_agc_settings);
                    end
                    WLAN_ad9361_simout = parsim(WLAN_ad9361_simin, 'ShowProgress', 'on');
                    % then, disable signal logging and save model
                    obj.disable_signal_logging();
                    save_system(obj.model);

                    % Extract log-data
                    for idx = 1:obj.num_sims
                        obj.error_message{idx} = WLAN_ad9361_simout(idx).ErrorMessage;
                        [obj.log_data{idx}, obj.log_data_indices{idx}] = obj.extract_log_data(num_logged_signals, WLAN_ad9361_simout(idx));
                        % Resample 9361 output to 1x rate for consumption by WLAN receiver
                        obj.input_waveform{idx} = resample(obj.log_data{idx}{obj.log_data_indices{idx}.AD9361_Output}.Data, N, M);
                    end
                end
%             else
%                 %% Physical radio-based study
%                 obj.TransmitRepeat = false;
%                 obj.TXGain = TXGains;
%
%                 results = RunDeployedDesign(testCase);
%
%
%
%             end

            if ~isempty(varargin{1})
                obj.AGC_MODE = ad9361_settings.AGC_MODE;
                obj.LOG_ADC_OUTPUT = ad9361_settings.LOG_ADC_OUTPUT;
                obj.SAVE_LOG_DATA = ad9361_settings.SAVE_LOG_DATA;
            end
        end

        %% functions applicable to simulation-based study
        % determine all combinations of AGC settings to apply
        function [local_agc_settings, num_combos] = agc_settings_all_combos(obj, agc_settings)
            AGC_fields = fields(agc_settings);
            tuned_settings = [];
            for ii = 1:length(AGC_fields)
                if (length(agc_settings.(AGC_fields{ii})) > 1)
                    tuned_settings = [tuned_settings; ii length(agc_settings.(AGC_fields{ii}))];
                end
            end
            if isempty(tuned_settings)
                local_agc_settings = agc_settings;
                num_combos = 1;
            else
                num_combos = prod(tuned_settings(:, 2));
                num_params = length(tuned_settings(:, 2));

                count = 1;
                sets = cell(num_params, 1);
                combos = zeros(num_combos, num_params);
                for ii = tuned_settings(:, 1).'
                    sets{count} = agc_settings.(AGC_fields{ii});
                    count = count+1;
                end
                combos_temp = cell(1,num_params);
                [combos_temp{:}] = ndgrid(sets{:});
                for ii = 1:num_params
                    combos(:, ii) = combos_temp{ii}(:);
                end

                for ii = 1:length(AGC_fields)
                    if ~ismember(ii, tuned_settings(:, 1))
                        local_agc_settings.(AGC_fields{ii}) = agc_settings.(AGC_fields{ii})*ones(1, num_combos);
                    end
                end
                count = 1;
                for ii = tuned_settings(:, 1).'
                    local_agc_settings.(AGC_fields{ii}) = combos(:, count).';
                    count = count+1;
                end
            end
        end

        % save the test waveform in a baseband file format
        function obj = gen_testWaveform_file(obj)
            fc = obj.fc;
            fs = obj.fs;
            M = obj.ovx(1);
            N = obj.ovx(2);

            warning('off', 'MATLAB:MKDIR:DirectoryExists');
            mkdir(obj.folder_name);
            obj.sig_filename = char('testWaveform_'+regexprep(string(datetime),'[:-\s]','_')+'.bb');
            obj.testwaveform_fileloc = [pwd '\' obj.folder_name '\' obj.sig_filename];
            obj.rxNonHT_2x = resample(obj.rxWaveform, M, N);
            bbw = comm.BasebandFileWriter(obj.testwaveform_fileloc,(M/N)*fs,fc);
            bbw.Metadata = struct('Date',date);
            bbw(obj.rxNonHT_2x);
            release(bbw);

            top_level_blocks = find_system(obj.model);
            for ii = 1:length(top_level_blocks)
                blocks_path = strsplit(top_level_blocks{ii}, '/');
                if (length(blocks_path) < 2)
                    continue;
                end
                str = regexprep(blocks_path{2},'[\n\r-\s_]','');
                if strcmpi(str, 'BasebandFileReader')
                    set_param(top_level_blocks{ii},'Filename',obj.testwaveform_fileloc);
                    break;
                end
            end

            % configure filter parameters
            set_param([obj.model '/' obj.block],'RF',num2str(fc));
            set_param([obj.model '/' obj.block],'LO',num2str(fc));
            set_param([obj.model '/' obj.block],'Trf', [num2str(obj.SampleTime_rx) '/' num2str(M/N)]);
        end

        % programmatically enable logging selected blocks in the Simulink
        % model
        function num_logged_signals = enable_signal_logging(obj)
            num_logged_signals = 0;
            for ii = 1:length(obj.logged_blocks)
                rx_under_mask = [obj.rx_block '/' obj.logged_blocks{ii}];
                ph = get_param(rx_under_mask,'PortHandles');
                str = regexprep(obj.logged_blocks{ii},'[\n\r-\s_]','');
                switch (str)
                    case {'LMTPeakDetector', 'ADCOverloadDetector', 'DDCFiltersRX'}
                        set_param(ph.Outport(1),'DataLogging','on');
                        set_param(ph.Outport(2),'DataLogging','on');
                        num_logged_signals = num_logged_signals+2;
                    otherwise
                        set_param(ph.Outport(1),'DataLogging','on');
                        num_logged_signals = num_logged_signals+1;
                end
            end
        end

        % programmatically disable logging selected blocks in the Simulink
        % model
        function disable_signal_logging(obj)
            for ii = 1:length(obj.logged_blocks)
                rx_under_mask = [obj.rx_block '/' obj.logged_blocks{ii}];
                ph = get_param(rx_under_mask,'PortHandles');
                str = regexprep(obj.logged_blocks{ii},'[\n\r-\s_]','');
                switch (str)
                    case {'LMTPeakDetector', 'ADCOverloadDetector', 'DDCFiltersRX'}
                        set_param(ph.Outport(1),'DataLogging','off');
                        set_param(ph.Outport(2),'DataLogging','off');
                    otherwise
                        set_param(ph.Outport(1),'DataLogging','off');
                end
            end
        end

        % programmatically set AGC parameters in the Simulink model
        function in = tune_agc_settings(obj, AGC_MODE, idx, in, AGC_settings)
            % Toggle manual switch to select AGC mode control input
            rx_top_level = find_system(obj.model);
            for ii = 1:length(rx_top_level)
                rx_top_level_path = strsplit(rx_top_level{ii}, '/');
                if (length(rx_top_level_path) < 2)
                    continue;
                end

                str = regexprep(rx_top_level_path{2},'[\n\r-\s_]','');
                if strcmpi(str, 'ManualSwitch')
                    if (AGC_MODE == 3)
                        in = in.setBlockParameter(rx_top_level{ii}, 'sw', '0');
                    else
                        in = in.setBlockParameter(rx_top_level{ii}, 'sw', '1');
                    end
                elseif strcmpi(str, 'Attackmode')
                    if (AGC_MODE ~= 3)
                        if (AGC_MODE == 0)
                            in = in.setBlockParameter(rx_top_level{ii}, 'Attackmode', 'Manual');
                        elseif (AGC_MODE == 1)
                            in = in.setBlockParameter(rx_top_level{ii}, 'Attackmode', 'Fast');
                        elseif (AGC_MODE == 2)
                            in = in.setBlockParameter(rx_top_level{ii}, 'Attackmode', 'Slow');
                        end
                    end
                end
            end

            rx_under_mask = find_system(obj.rx_block, 'LookUnderMasks', 'on', 'SearchDepth', 1);

            for ii = 1:length(rx_under_mask)
                logged_block_path = strsplit(rx_under_mask{ii}, '/');
                if (length(logged_block_path) < 3)
                    continue;
                end

                str = regexprep(logged_block_path{3},'[\n\r-\s_]','');
                if strcmpi(str, 'AGC')
                    % set AGC mask parameters
                    in = in.setBlockParameter(rx_under_mask{ii}, 'AGC_attack_delay', '1');
                    in = in.setBlockParameter(rx_under_mask{ii}, 'Max_Increase', num2str(AGC_settings.MaxIncrease(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'EnergyLost', num2str(AGC_settings.EnergyLostLevel(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'LowPower', num2str(AGC_settings.LowPwrTh(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'AGC_G_inc_fast', num2str(AGC_settings.AvgPwrLInc(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'AGC_G_inc_slow', num2str(AGC_settings.AvgPwrSInc(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'AGC_G_dec_fast', num2str(AGC_settings.AvgPwrLDec(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'AGC_G_dec_slow', num2str(AGC_settings.AvgPwrSDec(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'AGC_lock_level', num2str(AGC_settings.AGCLockLevel(idx)));
                elseif strcmpi(str, 'LMTPeakDetector')
                    in = in.setBlockParameter(rx_under_mask{ii}, 'HThreshold', num2str(AGC_settings.LMT_Hth(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'LThreshold', num2str(AGC_settings.LMT_Lth(idx)));
                elseif strcmpi(str, 'ADCOverloadDetector')
                    in = in.setBlockParameter(rx_under_mask{ii}, 'Hthreshold', num2str(AGC_settings.ADC_Hth(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'Lthreshold', num2str(AGC_settings.ADC_Lth(idx)));
                    in = in.setBlockParameter(rx_under_mask{ii}, 'Ncycles', num2str(AGC_settings.ADC_Ncycles(idx)));
                elseif strcmpi(str, 'AveragePowerMeter')
                    in = in.setBlockParameter(rx_under_mask{ii}, 'Ncycles', num2str(AGC_settings.AvgPwrMtr_Ncycles(idx)));
                end
            end
        end

        % extract log data from each of the logged blocks in the Simulink
        % model
        function [log_data, indices] = extract_log_data(obj, num_logged_signals, WLAN_ad9361_sim)
            log_data = cell(num_logged_signals, 1);
            for ii = 1:num_logged_signals
                if (isprop(WLAN_ad9361_sim,'logsout'))
                    simData_sigObj = WLAN_ad9361_sim.logsout{ii};
                else
                    continue;
                end
                logged_block_path = strsplit(simData_sigObj.BlockPath.getBlock(1), '/');
                block_name = regexprep(logged_block_path{3},'[\n\r-\s_]','');

                % Enable both output ports for some blocks and for the remaining,
                % enable the only output port
                switch block_name
                    case 'ADCOverloadDetector'
                        if (simData_sigObj.PortIndex == 1)
                           log_data{ii} = struct('Name', strcat(block_name,'Low'), 'Time', simData_sigObj.Values.Time);
                           indices.ADC_Peak_Detector_Low = ii;
                        elseif (simData_sigObj.PortIndex == 2)
                           log_data{ii} = struct('Name', strcat(block_name,'High'), 'Time', simData_sigObj.Values.Time);
                           indices.ADC_Peak_Detector_High = ii;
                        end
                    case 'LMTPeakDetector'
                        if (simData_sigObj.PortIndex == 1)
                           log_data{ii} = struct('Name', strcat(block_name,'High'), 'Time', simData_sigObj.Values.Time);
                           indices.LMT_Peak_Detector_High = ii;
                        elseif (simData_sigObj.PortIndex == 2)
                           log_data{ii} = struct('Name', strcat(block_name,'Low'), 'Time', simData_sigObj.Values.Time);
                           indices.LMT_Peak_Detector_Low = ii;
                        end
                    case 'DDCFiltersRX'
                        if (simData_sigObj.PortIndex == 1)
                           log_data{ii} = struct('Name', strcat(block_name,'OutputPwr'), 'Time', simData_sigObj.Values.Time);
                           indices.AD9361_Output_Power = ii;
                        elseif (simData_sigObj.PortIndex == 2)
                           log_data{ii} = struct('Name', strcat(block_name,'OutputSig'), 'Time', simData_sigObj.Values.Time);
                           indices.AD9361_Output = ii;
                        end
                    otherwise
                        log_data{ii} = struct('Name', block_name, 'Time', simData_sigObj.Values.Time);
                        if (strcmp(block_name, 'ADCRX'))
                            indices.ADC_Output = ii;
                        elseif strcmp(block_name, 'LPFRX')
                            indices.ADC_Input = ii;
                        elseif strcmp(block_name, 'In')
                            indices.AD9361_Input = ii;
                        elseif strcmp(block_name, 'RealImagtoComplex1')
                            indices.LMT_Output = ii;
                        elseif strcmp(block_name, 'AGC')
                            indices.AGC_Gain = ii;
                        elseif strcmp(block_name, 'AveragePowerMeter')
                            indices.Digital_Power = ii;
                        end
                end

                % if log-data is in (1x1xN) matrix format, extract data along the third dimension
                if (numel(size(simData_sigObj.Values.Data)) == 3)
                    log_data{ii}.Data = permute(simData_sigObj.Values.Data(1,1,:),[3 1 2]);
                else
                    log_data{ii}.Data = simData_sigObj.Values.Data;
                end
            end
        end

        %% functions applicable to physical radio-based study
        function results = RunDeployedDesign(obj)
            % waveform transmitted by the radio
            txWaveform = int16(2^15.*obj.txFrms./max(abs(obj.txFrms)));

            % HW setup
            [rx, tx] = obj.SetupHardware(obj.fs,obj.fc);
            obj.reTuneRadio(rx,tx,txWaveform);
            [rx, tx] = testCase.SetupHardware(obj.fs,obj.fc);

            % Transmit and Receive
            fprintf('Initializing RX and TX devices\n');
            if (obj.TransmitRepeat)
                tx.EnableCyclicBuffers = true;
                tx(txWaveform);
                pause(0.1);
            end
            % Clean buffer
            fprintf('Cleaning buffer \n');
            for k = 10:-1:0
                if ~obj.TransmitRepeat
                    % Send packets
                    pause(1);
                    tx(txWaveform);
                end
                rx();
                fprintf('%d ',k);
                pause(0.1);
            end

        end

        function reTuneRadio(testCase,rx,tx,data)
            rx.SamplingRate = fix(rx.SamplingRate*0.9);
            tx.SamplingRate = fix(tx.SamplingRate*0.9);
            rx.CenterFrequency = fix(rx.CenterFrequency*0.9);
            tx.CenterFrequency = fix(tx.CenterFrequency*0.9);
            tx(data);
            for l=1:10
                rx();
            end
            clear rx tx;
            system(['ssh -t root@',testCase.radioIP(4:end),' /root/reg/reg']);
        end

        function [rx,tx] = SetupHardware(obj,Rs,fc)
            % Receiver and transmitter objects
            rx = adi.AD9361.Rx;
            tx = adi.AD9361.Tx;
            rx.uri = testCase.radioIP;
            tx.uri = testCase.radioIP;
            rx.channelCount = 4;
            rx.SamplingRate = Rs;
            tx.SamplingRate = Rs;
            tx.AttenuationChannel0 = obj.TXGain;
            rx.RFBandwidth = fix(Rs*1.2);
            tx.RFBandwidth = fix(Rs*1.2);
            if (AGC_MODE == 0)
                rx.GainControlModeChannel0 = 'manual';
            elseif (AGC_MODE == 1)
                rx.GainControlModeChannel0 = 'slow_attack';
            elseif (AGC_MODE == 2)
                rx.GainControlModeChannel0 = 'fast_attack';
            end

            % Receiver Setup
            rx.DataTimeout = obj.TimeOut;
            rx.CenterFrequency = obj.fc;
            rx.SamplesPerFrame = obj.RXBufferSize;

            % Transmitter Setup
            tx.CenterFrequency = obj.fc;
        end

    end
end

author: analogdevicesinc


text:
function integralImage = computeIntegralImage(table)

integralImage = cumsum(table,1); integralImage = cumsum(integralImage,2);
[height width] = size(table);
%set the first row and the first column 0 in the integral image
integralImage =[zeros(height,1) integralImage];
integralImage=[zeros(1,width+1); integralImage];

author: Cloud-CV

repo: object-proposals


text:
 function [smap x y] = mri_sensemap_sim(varargin)
%function [smap x y] = mri_sensemap_sim(varargin)
%|
%| Simulate 2D sensitivity maps for sensitivity-encoded MRI
%| based on grivich:00:tmf doi:10.1119/1.19461
%|
%| option
%|	nx, ny, dx, dy, ncoil, rcoil, coil_distance, orbit (see below)
%|	chat	0|1	show images?
%|
%| out
%|	smap	[nx ny ncoil]	simulated sensitivity maps (complex!)
%|
%| See also ir_mri_sensemap_sim.m for new 3D version!
%|
%| Copyright 2005-6-20, Jeff Fessler and Amanda Funai, University of Michigan
%| 2014-08-19 JF more testing, verifying phase is correct

if nargin == 1 % tests
	if streq(varargin{1}, 'test')
		mri_sensemap_sim_test
	elseif streq(varargin{1}, 'test1')
		mri_sensemap_sim_test1
	end
return
end

arg.nx = 64;
arg.ny = [];
arg.dx = 3; % pixel size in mm
arg.dy = [];
arg.ncoil = 4; % # of coils
arg.rcoil = 100; % coil radius
arg.orbit = 360;
arg.coil_distance = 1.2; % multiplies fov/2
arg.flag_old = false; % old way based on wang:00:dop
arg.chat = nargout == 0;

arg = vararg_pair(arg, varargin);

if isempty(arg.dy), arg.dy = arg.dx; end
if isempty(arg.ny), arg.ny = arg.nx; end
if isempty(arg.rcoil), arg.rcoil = arg.dx * arg.nx / 2 * 0.50; end

[smap x y] = mri_sensemap_sim_do(arg.nx, arg.ny, arg.dx, arg.dy, ...
	arg.ncoil, arg.rcoil, arg.orbit, arg.coil_distance, ...
	arg.flag_old, arg.chat);

if ~nargout, clear, end

% mri_sensemap_sim_do()
function [smap x y] = mri_sensemap_sim_do(nx, ny, dx, dy, ncoil, rcoil, ...
		orbit, coil_distance, flag_old, chat)

rlist = rcoil * ones(ncoil, 1); % coil radii

plist = zeros(ncoil,3); % position of coil center [x y 0]
nlist = zeros(ncoil,3); % normal vector (inward) from coil center

% circular coil configuration, like head coils
alist = deg2rad(orbit)/ncoil * [0:(ncoil-1)]; % list of coil angles in radians
for ii=1:ncoil
	phi = alist(ii);
	Rad = max(nx/2 * dx, ny/2 * dy) * coil_distance;
	plist(ii,:) = Rad * [cos(phi) sin(phi) 0];
	nlist(ii,:) = -[cos(phi) sin(phi) 0];
	olist(ii,:) = [-sin(phi) cos(phi) 0]; % unit vector orthogonal to nlist
end

% object coordinates for slice z=0
x = ([1:nx] - (nx+1)/2) * dx;
y = ([1:ny] - (ny+1)/2) * dy;
z = 0;
[xx yy zz] = ndgrid(x,y,z);

smap = zeros(nx, ny, ncoil, 'single');
for ii=1:ncoil
	% rotate coordinates to correspond to coil orientation
	zr =	(xx - plist(ii,1)) .* nlist(ii,1) + ...
		(yy - plist(ii,2)) .* nlist(ii,2) + ...
		(zz - plist(ii,3)) .* nlist(ii,3);
	xr =	xx .* nlist(ii,2) - yy .* nlist(ii,1);

	if 0 % see coordinates
		im plc 1 2
		im(1, x, y, xr), xlabel x, ylabel y
		im(2, x, y, zr)
		keyboard
	end

	[sx sy sz] = mri_smap1(xr, 0, zr, rlist(ii)); % in coil coordinates

	if 0 % see field components
		im plc 2 2
		im(1, x, y, sx), cbar
		im(2, x, y, sy), cbar
		im(3, x, y, sz), cbar
		im subplot 2
		tmp = sqrt(sx.^2 + sz.^2);
		quiver(x, y, (sx./tmp)', (sz./tmp)', 0), axis square
	end

	if flag_old
		smap(:,:,ii) = sz; % old way (wrong!) based only on smap_z

	else
		if nlist(ii,3) || olist(ii,3)
			fail 'unsupported'
		end
		% assume z component of plist and nlist are 0
		bx = sz * nlist(ii,1) + sx * olist(ii,1);
		by = sz * nlist(ii,2) + sx * olist(ii,2);
		smap(:,:,ii) = bx + 1i * by;

		if 0 % see final field components vs phase
			im subplot 4
			bb = sqrt(bx.^2 + by.^2);
			quiver(x, y, (bx./bb)', (by./bb)', 0), axis square
			im(2, x, y, angle(smap(:,:,ii))), cbar
			keyboard
		end
	end
end
smap = smap * rlist(1) / (2*pi); % trick: scale so maximum is near unity

if chat && im % show smap and array geometry in z=0 plane
	mri_sensemap_sim_show(smap, x, y, dx, dy, nlist, plist, rlist)
end

% mri_sensemap_sim_show()
function mri_sensemap_sim_show(smap, x, y, dx, dy, nlist, plist, rlist)
[nx ny ncoil] = size(smap);
nshow = min(max(ncoil,2),4);
im('plc', 3, nshow)
for ii=1:min(ncoil,nshow)
	clim = [0 max(abs(smap(:)))];
	tmp = smap(:,:,ii);
	im(ii, x, y, abs(tmp), clim, 'Magnitude'), cbar
%	xmax = max(max(abs(x)), max(plist(:,1)));
%	ymax = max(max(abs(y)), max(plist(:,2)));
	xmax = max([max(abs(x)) max(abs(y)) max(col(plist(:,[1 2])))]);
%	axis([-xmax xmax -ymax ymax]*1.05)
	axis(xmax * [-1 1 -1 1] * 1.1)
	xtick([-1 0 1] * nx/2 * dx)
	ytick([-1 0 1] * ny/2 * dy)

	hold on
	plot(0,0,'.', plist(:,1), plist(:,2), 'o')
	xdir = nlist(ii,2);
	ydir = nlist(ii,1);
	r = rlist(ii);
	plot(plist(ii,1)+r*xdir*[-1 1], plist(ii,2)+r*ydir*[1 -1], '-')
	hold off

%	ph = ir_unwrap(angle(tmp)); % trick: unwrap phase for pretty display
	ph = angle(tmp); % show raw phase (understandable with hsv colormap)
	im(ii+nshow, 'hsv', x, y, ph, [-pi pi], 'Phase'), cbar
	axis(xmax*[-1 1 -1 1]*1.1)
	xtick([-1 0 1] * nx/2 * dx)
	ytick([-1 0 1] * ny/2 * dy)
end

ssos = sqrt(sum(abs(smap), 3)); % bug! missing .^2 !
minmax(ssos)
ssos = ssos / ssos(end/2,end/2);
im(2*nshow+1, x, y, ssos, 'SSoS (normalized)'), cbar
xtick([-1 0 1] * nx/2 * dx)
ytick([-1 0 1] * ny/2 * dy)

if ncoil == 1
	subplot(122)
	bx = real(smap);
	by = imag(smap);
	quiver(x, y, bx', by'), title 'Field pattern in x-y plane'
	axis equal, axis tight
end

% clf, im(angle(smap(:,:,2))), cbar

% mri_smap_r(r, z)
% function for testing near 0
function out = mri_smap_r(r, z)
M = 4 * r ./ ((1 + r).^2 + z.^2); % = k^2, see ellipke
[K E] = ellipke(M);
out = 2 * z ./ r .* ((1 + r).^2 + z.^2).^(-0.5) .* ...
	((1 + r.^2 + z.^2) ./ ((1 - r).^2 + z.^2) .* E - K);

% mri_smap1()
% based on grivich:00:tmf
% for a circular coil in "x-y plane" of radius a
% note that coil x-y plane is not same as object x-y plane!
function [smap_x smap_y smap_z] = mri_smap1(x, y, z, a)
x = x ./ a; % normalized units
y = y ./ a;
z = z ./ a;
r = sqrt(x.^2 + y.^2);
M = 4 * r ./ ((1 + r).^2 + z.^2); % = k^2, see ellipke
[K E] = ellipke(M);
if ir_is_octave
	K = reshape(K, size(M)); % ellipke shape bug in octave
	E = reshape(E, size(M));
end

% the following is B_z in eqn (18) in grivich:00:tmf
% and same as eqn [10] in wang:00:dop to within constant scale factor
smap_z = 2 * ((1 + r).^2 + z.^2).^(-0.5) .* ...
	(K + (1 - r.^2 - z.^2) ./ ((1 - r).^2 + z.^2) .* E);
smap_z = smap_z / a;

if 0 && any(r(:) == 0) % test code to explore when r is near 0
	r0 = linspace(0,5e-7,101);
	z0 = 0.4;
	t0 = mri_smap_r(r0, z0);
	slope = 3*pi * z0 / ((1+z0^2)^2.5);
	clf, plot(r0, t0, '-', r0, slope * r0, '--'); grid, prompt
end

% the following is B_r in eqn (17) in grivich:00:tmf
smap_r = 2 * z ./ r .* ((1+r).^2 + z.^2).^(-0.5) .* ...
	((1 + r.^2 + z.^2) ./ ((1-r).^2 + z.^2) .* E - K);
bad = abs(r) < 1e-6;
smap_r(bad) = 3 * pi * z(bad) ./ ((1 + z(bad).^2).^2.5) .* r(bad);
smap_r = smap_r / a;

if any(isnan(smap_r(:))) || any(isnan(smap_z(:)))
	keyboard
end

smap_x = smap_r .* div0(x, r);
smap_y = smap_r .* div0(y, r);

%phi = atan2(y, x);
%smap_x = smap_r .* cos(phi);
%smap_y = smap_r .* sin(phi);

% mri_sensemap_sim_test1
% test mri_smap1 routine, cf Fig. 4 of grivich:00:tmf
function mri_sensemap_sim_test1
a = 1;
x = linspace(-2,2,99);
y = linspace(-2,2,97);
zlist = [0.1 0.2 0.5 1.0];
[xx yy zz] = ndgrid(x, y, zlist);
[smap_x smap_y smap_z] = mri_smap1(xx, yy, zz, a);
im('plc', 4, numel(zlist))
mri_sensemap_sim_test1_show(smap_x, x, y, 0, zlist, 'x')
mri_sensemap_sim_test1_show(smap_y, x, y, 4, zlist, 'y')
mri_sensemap_sim_test1_show(smap_z, x, y, 8, zlist, 'z')
smap_b = sqrt(smap_x.^2 + smap_y.^2);
mri_sensemap_sim_test1_show(smap_b, x, y, 12, zlist, 'b')

% mri_sensemap_sim_test1_show()
function mri_sensemap_sim_test1_show(map, x, y, offset, zlist, leg)
clim = [-20 20];
for iz = 1:numel(zlist)
	p = (iz-1) * 4;
	im(offset+iz, x, y, map(:,:,iz), leg, clim), cbar
	axis equal
	if iz == 1
		xtick([-2 2]), ytick([-2 2])
	else
		xtick off, ytick off
	end
	if zlist(iz) < 0.5
		blim = [7 12 19];
	else
		blim = [1 3 5];
	end
	hold on
	contour(x, y, abs(map(:,:,iz))', blim, 'b-')
	contour(x, y, abs(map(:,:,iz))', [0 0]+0.001, 'g-')
	hold off
end
drawnow

% mri_sensemap_sim_test
function mri_sensemap_sim_test

[smap x y] = mri_sensemap_sim('chat', 1, 'nx', 32, ...
	'rcoil', [], 'ncoil', 4, 'coil_distance', 1.2);

if 1 % check rotational symmetry in 4-coil case
	for ic=2:4
		tmp = rot90(smap(:,:,1), (ic-1));
		equivs(abs(tmp), abs(smap(:,:,ic)))
	%	im(8+ic, x, y, abs(tmp)), cbar
		p1 = angle(tmp) + (ic-1) * pi/2; % add pi/2 to rotated
		p2 = angle(smap(:,:,ic));
		tmp = exp(1i * (p2 - p1));
		equivs(tmp, ones(size(tmp))) % trick: equivs mod 2*pi
	end
end

if 0 && im % see ellipke
	m = linspace(0,1,101);
	[k e] = ellipke(m);
	clf, plot(m, k, '-', m, e, '--'), legend('k', 'e')
	yaxis_pi('0 p/2 p 3*p/2')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function submit()
  addpath('./lib');

  conf.assignmentSlug = 'anomaly-detection-and-recommender-systems';
  conf.itemName = 'Anomaly Detection and Recommender Systems';
  conf.partArrays = { ...
    { ...
      '1', ...
      { 'estimateGaussian.m' }, ...
      'Estimate Gaussian Parameters', ...
    }, ...
    { ...
      '2', ...
      { 'selectThreshold.m' }, ...
      'Select Threshold', ...
    }, ...
    { ...
      '3', ...
      { 'cofiCostFunc.m' }, ...
      'Collaborative Filtering Cost', ...
    }, ...
    { ...
      '4', ...
      { 'cofiCostFunc.m' }, ...
      'Collaborative Filtering Gradient', ...
    }, ...
    { ...
      '5', ...
      { 'cofiCostFunc.m' }, ...
      'Regularized Cost', ...
    }, ...
    { ...
      '6', ...
      { 'cofiCostFunc.m' }, ...
      'Regularized Gradient', ...
    }, ...
  };
  conf.output = @output;

  submitWithConfiguration(conf);
end

function out = output(partId, auxstring)
  % Random Test Cases
  n_u = 3; n_m = 4; n = 5;
  X = reshape(sin(1:n_m*n), n_m, n);
  Theta = reshape(cos(1:n_u*n), n_u, n);
  Y = reshape(sin(1:2:2*n_m*n_u), n_m, n_u);
  R = Y > 0.5;
  pval = [abs(Y(:)) ; 0.001; 1];
  Y = (Y .* double(R));  % set 'Y' values to 0 for movies not reviewed
  yval = [R(:) ; 1; 0];
  params = [X(:); Theta(:)];
  if partId == '1'
    [mu sigma2] = estimateGaussian(X);
    out = sprintf('%0.5f ', [mu(:); sigma2(:)]);
  elseif partId == '2'
    [bestEpsilon bestF1] = selectThreshold(yval, pval);
    out = sprintf('%0.5f ', [bestEpsilon(:); bestF1(:)]);
  elseif partId == '3'
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 0);
    out = sprintf('%0.5f ', J(:));
  elseif partId == '4'
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 0);
    out = sprintf('%0.5f ', grad(:));
  elseif partId == '5'
    [J] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                       n, 1.5);
    out = sprintf('%0.5f ', J(:));
  elseif partId == '6'
    [J, grad] = cofiCostFunc(params, Y, R, n_u, n_m, ...
                             n, 1.5);
    out = sprintf('%0.5f ', grad(:));
  end
end

author: Ayatans

repo: Machine-Learning-homework


text:
classdef DesignVariableAcademic < handle

    properties (Access = public)
        value
        valueOld
    end

    properties (Access = private)

    end

    properties (Access = private)

    end

    methods (Access = public)

        function obj = DesignVariableAcademic()

        end

        function init(obj,x0)
            obj.value = x0;
        end

    end

    methods (Access = public)

        function update(obj,x)
            obj.value = x;
        end

        function updateOld(obj)
            obj.valueOld = obj.value;
        end

    end

end

author: SwanLab

repo: Swan


text:
function [ patch_expert, corr, rms_error] = CreatePatchExpert( samples, labels, unnormed_samples, normalisation_options)
%CREATEPATCHEXPERT Summary of this function goes here
%   Detailed explanation goes here

    num_examples = size(samples, 1);

    region_length = normalisation_options.normalisationRegion - normalisation_options.patchSize + 1;
    region_length = region_length(1) * region_length(2);

    % this part sets the split boundaries for SVR training, logit training and test subsets
    train_SVR_start = 1;
    train_SVR_end = int32(normalisation_options.svmRatio * num_examples - 1);
    train_SVR_end = train_SVR_end - mod(train_SVR_end, region_length);

    train_logit_start = train_SVR_end + 1;
    train_logit_end = train_logit_start + int32(normalisation_options.logitRatio * num_examples) - 1;

    % make sure it ends within same area of interest (region)
    train_logit_end = train_logit_end - mod(train_logit_end, region_length);

    test_start = train_logit_end + 1;
    test_end = size(samples,1);

    % picking training data for SVR
    examples_train_SVR = samples(train_SVR_start:train_SVR_end, :);
    labels_train_SVR = labels(train_SVR_start:train_SVR_end);

    % Train the SVR using liblinear (older version used libSVM,
    % but liblinear is much faster)
    [w, b] = Train_SVR(examples_train_SVR, labels_train_SVR);

    % Now create the test dataset
    examples_test = samples(test_start:test_end, :);
    labels_test = labels(test_start:test_end);

    % Training the logistic regressor now
    examples_train_logit = samples(train_logit_start:train_logit_end, :);
    labels_train_logit = labels(train_logit_start:train_logit_end, :);

    patch_expert = CreateLogisticRegressor(examples_train_logit, labels_train_logit, w, b, normalisation_options);

    % Evaluate the patch expert now
    [rms_error, corr, ~] = EvaluatePatchExpert(examples_test, labels_test, patch_expert, false);

    fprintf('Rms error %.3f, correlation %.3f\n', rms_error, corr);

    % Assert that our implementation and the convolution based one are equivalent
    [~, ~, responses_svm] = EvaluatePatchExpert(samples(1:size(unnormed_samples,1)*region_length,:), labels(1:size(unnormed_samples,1)*region_length), patch_expert, false);
    [responses_ncc] = SVR_expert_ncc_response(unnormed_samples, patch_expert, normalisation_options, normalisation_options.normalisationRegion, region_length);
    assert(mean(abs(responses_svm-responses_ncc))< 1e-2);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function i4_modp_test ( )

%*****************************************************************************80
%
%% I4_MODP_TEST tests I4_MODP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 4;

  ndivid = [ 50, -50, 50, -50 ];
  number = [ 107, 107, -107, -107 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_MODP_TEST\n' );
  fprintf ( 1, '  I4_MODP factors a number\n' );
  fprintf ( 1, '  into a multiple and a remainder.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Number   Divisor  Multiple Remainder\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    nrem = i4_modp ( number(test), ndivid(test) );
    nmult = floor ( ( number(test) - nrem ) / ndivid(test) );
    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
      number(test), ndivid(test), nmult, nrem );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat using MOD:\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    nrem = mod ( number(test), ndivid(test) );
    nmult = floor ( number(test) / ndivid(test) );
    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
      number(test), ndivid(test), nmult, nrem );
  end

  return
end

author: johannesgerer


text:
function [stateLogCol] = ma_resNumToStateLogCol(resNum)
%ma_resNumToStateLogCol Summary of this function goes here
%   Detailed explanation goes here

    switch resNum
        case 1
            stateLogCol = 10;
        case 2
            stateLogCol = 11;
        case 3
            stateLogCol = 12;
        otherwise
            error('Resource number %u has no mapping to a state log column number!', resNum);
    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function zscoredData = BF_zscore(inputData)
% BF_zscore     z-score the input data vector (without a Statistics Toolbox licence)
%
%---INPUT:
% inputData, the input time series (or any vector).
%
%---OUTPUT:
% zscoredData, the z-scored transformation of the input.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% Check for NaNs:
if any(isnan(inputData))
    error('inputData contains NaNs');
end

% By default, z-score twice to reduce the numerical error:
zscoredData = (inputData - mean(inputData)) / std(inputData);
zscoredData = (zscoredData - mean(zscoredData)) / std(zscoredData);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [Psi, Astats, Kstats] = sampleHMMhypers( Psi, algParams )
% Sample hyperparameters (alpha,kappa) that control HMM transition distr.
% Uses a Metropolis-Hastings random walk proposal centered on current state
% INPUT:
%   Psi : input Markov Chain state
% OUTPUT
%   Psi : new Markov chain state
% DETAILS:
%  Given hyperparams alpha (non-sticky) and kappa (sticky) params,
%   each sequence draws transition probabilities as follows:
%     eta(jj,:) ~ Gamma( alpha, alpha, ... alpha+kappa, alpha, alpha )
%     pi(jj,:) = eta(jj,:)/sum( eta(jj,:)

%------------------------ Unpack
F = Psi.F;
[N,K] = size(F);
alpha0 = Psi.TransM.prior.alpha;
kappa0 = Psi.TransM.prior.kappa;
TransM   = Psi.TransM;

% Hyperparameters for prior on kappa:
a_kappa = TransM.prior.a_kappa;
b_kappa = TransM.prior.b_kappa;
% Variance of gamma proposal:
var_kappa = algParams.HMM.var_kappa;

% Hyperparameters for prior on alpha:
a_alpha = TransM.prior.a_alpha;
b_alpha = TransM.prior.b_alpha;
% Variance of gamma proposal:
var_alpha = algParams.HMM.var_alpha;

% Build sufficient stats
%    Ki( ii ) : # total features on for object ii
%    Skk( ii ) : sum( log(   diag entries of Pz(ii)  )  )
%    Sall( ii )    : sum( log(   all entries of Pz(ii)  ) )
Ki = zeros(1,N);
Skk = zeros(1,N);
Sall = zeros(1,N);
for ii=1:N
    Ki(ii) = sum(F(ii,:));
    pi_ii = TransM.pi( ii );

    Skk( ii ) = sum( log( diag( pi_ii ) ) );
    Sall( ii ) = sum( log(  pi_ii(:) ) );
end

Astats.nAccept = 0;
Astats.nTotal   = algParams.HMM.Niter;

Kstats.nAccept = 0;
Kstats.nTotal   = algParams.HMM.Niter;

for nn=1:algParams.HMM.Niter

    %%%%%%% Sample kappa given alpha %%%%%%%

    % (a,b) hyperparameters of gamma prior based on fixed variance and setting
    % mean equal to previous kappa value:
    aa_kappa0 = (kappa0^2)/var_kappa;
    bb_kappa0 = kappa0/var_kappa;

    % Sample a proposed kappa:
    kappaP = randgamma(aa_kappa0) / bb_kappa0;

    % Determine log-likelihood of transition distributions given previous kappa
    % value and proposed kappa value:
    log_diff_Z = sum( ...
              Ki .* ( gammaln(alpha0*Ki+kappaP) - gammaln(alpha0*Ki+kappa0) )...
            - Ki .* ( gammaln(alpha0+kappaP) - gammaln(alpha0+kappa0)) ...
                            + (kappaP-kappa0)*Skk ...
                    );
    % Add in prior probability of previous and proposed kappa values:
    log_diff_Prior = (a_kappa-1)*(log(kappaP)-log(kappa0))-(kappaP-kappa0)*b_kappa;

    % (a,b) hyperparameters of gamma prior based on fixed variance and setting
    % mean equal to proposed kappa value:
    aa_kappaP = (kappaP^2)/var_kappa;

    log_diff_Q  = (gammaln(aa_kappa0) - gammaln(aa_kappaP))...
                  + (aa_kappaP-aa_kappa0-1)*log(kappa0) - (aa_kappa0-aa_kappaP-1)*log(kappaP)...
                  + (aa_kappa0-aa_kappaP)*log(var_kappa);

    % Log accept-reject ratio:
    log_rho = log_diff_Z + log_diff_Prior + log_diff_Q;

    if isinf(log_rho)
        log_rho = -Inf;
    end
    rho = exp(log_rho);

    if rand < rho
        kappa0 = kappaP;
        Kstats.nAccept = Kstats.nAccept + 1;
    end
    % otherwise, just keep kappa0 to current value

    %%%%%%% Sample alpha given kappa %%%%%%%

    % (a,b) hyperparameters of gamma prior based on fixed variance and setting
    % mean equal to previous alpha value:
    aa_alpha0 = (alpha0^2)/var_alpha;
    bb_alpha0 = alpha0/var_alpha;

    % Sample a proposed alpha:
    alphaP = randgamma(aa_alpha0) / bb_alpha0;

    % Determine log-likelihood of transition distributions given previous alpha
    % value and proposed alpha value:
    log_diff_Z = sum( ...
                          Ki .* ( gammaln( alphaP*Ki + kappa0 ) - gammaln( alpha0*Ki + kappa0 )) ...
                        - Ki .* ( gammaln( alphaP+kappa0 )      - gammaln( alpha0+kappa0 ) ) ...
                        - Ki .* (Ki-1) .* ( gammaln( alphaP ) - gammaln( alpha0 )  ) ...
                        + ( alphaP - alpha0 ) .* Sall ...
                      );

    % Add in prior probability of previous and proposed alpha values:
    log_diff_Prior = (a_alpha-1)*(log(alphaP)-log(alpha0))-(alphaP-alpha0)*b_alpha;

    % (a,b) hyperparameters of gamma prior based on fixed variance and setting
    % mean equal to proposed kappa value:
    aa_alphaP = (alphaP^2)/var_alpha;
    %bb_alpha = alphaP/var_alpha; % seems to be unused

    log_diff_Q = (gammaln(aa_alpha0) - gammaln(aa_alphaP))...
        + (aa_alphaP-aa_alpha0-1)*log(alpha0) - (aa_alpha0-aa_alphaP-1)*log(alphaP)...
        + (aa_alpha0-aa_alphaP)*log(var_alpha);

    % Log accept-reject ratio:
    log_rho = log_diff_Z + log_diff_Prior + log_diff_Q;

    if isinf(log_rho)
        log_rho = -Inf;
    end
    rho = exp(log_rho);

    if rand < rho
        alpha0 = alphaP;
        Astats.nAccept = Astats.nAccept + 1;
    end
    % otherwise, just keep previous value of alpha0

end % end iterative loop over params

% -------------------------------  Repack
Psi.TransM.prior.alpha = alpha0;
Psi.TransM.prior.kappa = kappa0;

author: michaelchughes

repo: NPBayesHMM


text:
%% Copyright (C) 2016-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun logint (@var{x})
%% Numerical logint function.
%%
%% Example:
%% @example
%% @group
%% logint (1.1)
%%   @result{} ans = -1.6758
%% @end group
%% @end example
%%
%% @strong{Note} this function may be slow for large numbers of inputs.
%% This is because it is not a native double-precision implementation
%% but rather the numerical evaluation of the Python @code{mpmath} function
%% @code{li}.
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @seealso{@@sym/logint}
%% @end defun

function y = logint (x)
  if (nargin ~= 1)
    print_usage ();
  end
  cmd = { 'L = _ins[0]'
          'A = [complex(mpmath.li(x)) for x in L]'
          'return A,' };
  c = pycall_sympy__ (cmd, num2cell (x(:)));
  y = reshape (cell2mat (c), size (x));
end

%!error logint (1, 2)

%!test
%! x = 1.1;
%! y = sym(11)/10;
%! A = logint (x);
%! B = double (logint (y));
%! assert (A, B, -4*eps);

%!test
%! y = [2 3 sym(pi); exp(sym(1)) 5 6];
%! x = double (y);
%! A = logint (x);
%! B = double (logint (y));
%! assert (A, B, -4*eps);

%!test
%! % maple:
%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I];
%! % > for a in A do evalf(Li(a)) end do;
%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
%! A = [  1.3876787420229375511 + 2.5087546988592328752*1i
%!        1.6987684473874802274 + 4.5936366057115204667*1i
%!        30.126141584079629926
%!        3.4936715673748995398 + 5.5260023797127391973*1i
%!        0.90264689772681592152e-5 + 3.1415953634267361942*1i
%!       -2.3996350854560916779 - 7.6971739096353664559*1i ];
%! B = logint (x);
%! assert (A, B, -eps)

%!assert (logint (inf), inf)
%!assert (isnan (logint (-inf)))

repo: octsympy