text:
function data = loadmsgpack(fname,varargin)
%
%   data = loadmsgpack(fname,varargin)
%
%LOADMSGPACK parses a msgpack byte buffer into Matlab data structures
% LOADMSGPACK(BYTES)
%    reads BYTES as msgpack data, and creates Matlab data structures
%    from it.
%    - strings are converted to strings
%    - numbers are converted to appropriate numeric values
%    - true, false are converted to logical 1, 0
%    - nil is converted to []
%    - arrays are converted to cell arrays
%    - maps are converted to containers.Map
%
% (c) 2016 Bastian Bechtold
% modified by Qianqian Fang <q.fang at neu.edu>
%
% license:
%     BSD 3-clause license or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

    if(exist(fname,'file'))
       fid = fopen(fname,'rb');
       bytes = fread(fid,inf,'uint8=>char')';
       fclose(fid);
    else
       bytes=fname;
    end

    opt=varargin2struct(varargin{:});
    opt.simplifycell=jsonopt('SimplifyCell',1,opt);
    opt.simplifycellarray=jsonopt('SimplifyCellArray',0,opt);

    jsoncount=1;
    idx=0;
    while idx <= length(bytes)
        [obj, idx] = parse(uint8(bytes(:)), 1, opt);
        data{jsoncount}=obj;
        jsoncount=jsoncount+1;
    end

    jsoncount=length(data);
    if(jsoncount==1 && iscell(data))
        data=data{1};
    end
    if(iscell(data))
        data=cellfun(@(x) jdatadecode(x),data,'UniformOutput',false);
    elseif(isstruct(data))
        data=jdatadecode(data,'Base64',0, opt);
    end
end

function [obj, idx] = parse(bytes, idx, varargin)
    % masks:
    b10000000 = 128;
    b01111111 = 127;
    b11000000 = 192;
    b00111111 = 63;
    b11100000 = 224;
    b00011111 = 31;
    b11110000 = 240;
    b00001111 = 15;
    % values:
    b00000000 = 0;
    b10010000 = 144;
    b10100000 = 160;

    currentbyte = bytes(idx);

    if bitand(b10000000, currentbyte) == b00000000
        % decode positive fixint
        obj = int8(currentbyte);
        idx = idx + 1;
        return
    elseif bitand(b11100000, currentbyte) == b11100000
        % decode negative fixint
        obj = typecast(currentbyte, 'int8');
        idx = idx + 1;
        return
    elseif bitand(b11110000, currentbyte) == b10000000
        % decode fixmap
        len = double(bitand(b00001111, currentbyte));
        [obj, idx] = parsemap(len, bytes, idx+1, varargin{:});
        return
    elseif bitand(b11110000, currentbyte) == b10010000
        % decode fixarray
        len = double(bitand(b00001111, currentbyte));
        [obj, idx] = parsearray(len, bytes, idx+1, varargin{:});
        return
    elseif bitand(b11100000, currentbyte) == b10100000
        % decode fixstr
        len = double(bitand(b00011111, currentbyte));
        [obj, idx] = parsestring(len, bytes, idx + 1);
        return
    end

    switch currentbyte
        case 192 % nil
            obj = [];
            idx = idx+1;
      % case 193 % unused
        case 194 % false
            obj = false;
            idx = idx+1;
        case 195 % true
            obj = true;
            idx = idx+1;
        case 196 % bin8
            len = double(bytes(idx+1));
            [obj, idx] = parsebytes(len, bytes, idx+2);
        case 197 % bin16
            len = double(bytes2scalar(bytes(idx+1:idx+2), 'uint16'));
            [obj, idx] = parsebytes(len, bytes, idx+3);
        case 198 % bin32
            len = double(bytes2scalar(bytes(idx+1:idx+4), 'uint32'));
            [obj, idx] = parsebytes(len, bytes, idx+5);
        case 199 % ext8
            len = double(bytes(idx+1));
            [obj, idx] = parseext(len, bytes, idx+1);
        case 200 % ext16
            len = double(bytes2scalar(bytes(idx+1:idx+2), 'uint16'));
            [obj, idx] = parseext(len, bytes, idx+3);
        case 201 % ext32
            len = double(bytes2scalar(bytes(idx+1:idx+4), 'uint32'));
            [obj, idx] = parseext(len, bytes, idx+5);
        case 202 % float32
            obj = bytes2scalar(bytes(idx+1:idx+4), 'single');
            idx = idx+5;
        case 203 % float64
            obj = bytes2scalar(bytes(idx+1:idx+8), 'double');
            idx = idx+9;
        case 204 % uint8
            obj = bytes(idx+1);
            idx = idx+2;
        case 205 % uint16
            obj = bytes2scalar(bytes(idx+1:idx+2), 'uint16');
            idx = idx+3;
        case 206 % uint32
            obj = bytes2scalar(bytes(idx+1:idx+4), 'uint32');
            idx = idx+5;
        case 207 % uint64
            obj = bytes2scalar(bytes(idx+1:idx+8), 'uint64');
            idx = idx+9;
        case 208 % int8
            obj = bytes2scalar(bytes(idx+1), 'int8');
            idx = idx+2;
        case 209 % int16
            obj = bytes2scalar(bytes(idx+1:idx+2), 'int16');
            idx = idx+3;
        case 210 % int32
            obj = bytes2scalar(bytes(idx+1:idx+4), 'int32');
            idx = idx+5;
        case 211 % int64
            obj = bytes2scalar(bytes(idx+1:idx+8), 'int64');
            idx = idx+9;
        case 212 % fixext1
            [obj, idx] = parseext(1, bytes, idx+1);
        case 213 % fixext2
            [obj, idx] = parseext(2, bytes, idx+1);
        case 214 % fixext4
            [obj, idx] = parseext(4, bytes, idx+1);
        case 215 % fixext8
            [obj, idx] = parseext(8, bytes, idx+1);
        case 216 % fixext16
            [obj, idx] = parseext(16, bytes, idx+1);
        case 217 % str8
            len = double(bytes(idx+1));
            [obj, idx] = parsestring(len, bytes, idx+2);
        case 218 % str16
            len = double(bytes2scalar(bytes(idx+1:idx+2), 'uint16'));
            [obj, idx] = parsestring(len, bytes, idx+3);
        case 219 % str32
            len = double(bytes2scalar(bytes(idx+1:idx+4), 'uint32'));
            [obj, idx] = parsestring(len, bytes, idx+5);
        case 220 % array16
            len = double(bytes2scalar(bytes(idx+1:idx+2), 'uint16'));
            [obj, idx] = parsearray(len, bytes, idx+3, varargin{:});
        case 221 % array32
            len = double(bytes2scalar(bytes(idx+1:idx+4), 'uint32'));
            [obj, idx] = parsearray(len, bytes, idx+5, varargin{:});
        case 222 % map16
            len = double(bytes2scalar(bytes(idx+1:idx+2), 'uint16'));
            [obj, idx] = parsemap(len, bytes, idx+3, varargin{:});
        case 223 % map32
            len = double(bytes2scalar(bytes(idx+1:idx+4), 'uint32'));
            [obj, idx] = parsemap(len, bytes, idx+5, varargin{:});
        otherwise
            error('JSONLAB:MSGPACK:InvalidFormat', ...
                  ['Unknown type "' dec2bin(currentbyte) '"']);
    end
end

function value = bytes2scalar(bytes, type)
    % reverse byte order to convert from little-endian to big-endian
    value = typecast(bytes(end:-1:1), type);
end

function [str, idx] = parsestring(len, bytes, idx)
    if(~isoctavemesh)
         str = native2unicode(bytes(idx:idx+len-1)', 'utf-8');
    else
         str = char(bytes(idx:idx+len-1)');
    end
    idx = idx + len;
end

function [out, idx] = parsebytes(len, bytes, idx)
    out = bytes(idx:idx+len-1);
    idx = idx + len;
end

function [out, idx] = parseext(len, bytes, idx)
    obj.type = bytes(idx);
    obj.data = bytes(idx+1:idx+len);
    idx = idx + len + 1;
end

function [out, idx] = parsearray(len, bytes, idx, varargin)
    out = cell(1,len);
    for n=1:len
        [out{n}, idx] = parse(bytes, idx, varargin{:});
    end
    if(len==1)
        out=out{1};
    end
    if(varargin{1}.simplifycell)
        if(iscell(out) && ~isempty(out) && isnumeric(out{1}))
          if(all(cellfun(@(e) isequal(size(out{1}), size(e)) , out(2:end))))
              try
                  oldobj=out;
                  if(iscell(out) && length(out)>1 && ndims(out{1})>=2)
                      catdim=size(out{1});
                      catdim=ndims(out{1})-(catdim(end)==1)+1;
                      out=cat(catdim,out{:});
                      out=permute(out,ndims(out):-1:1);
                  else
                      out=cell2mat(out')';
                  end
                  if(iscell(oldobj) && isstruct(out) && numel(out)>1 && varargin{1}.simplifycellarray==0)
                      out=oldobj;
                  end
              catch
              end
          end
        end
        if(~iscell(out) && size(out,2)>1 && ndims(out)==2)
            out=out';
        end
    end
end

function [out, idx] = parsemap(len, bytes, idx, varargin)
    out = struct();
    for n=1:len
        [key, idx] = parse(bytes, idx, varargin{:});
        [out.(encodevarname(char(key))), idx] = parse(bytes, idx, varargin{:});
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [candidates, scores] = get_candidates(candidate_dir,method_config, img_id, num_candidates, allow_filtering)

  if nargin < 4
    allow_filtering = true;
  end

  [candidates, scores, rerun_num_candidates] = read_candidates_mat(candidate_dir, img_id);
  if iscell(candidates) && iscell(scores)
    [~,idx] = min(abs(rerun_num_candidates - num_candidates));
    candidates = candidates{idx};
    scores = scores{idx};
  end

  if allow_filtering
    if strcmp(method_config.opts.order, 'none')
      % nothing to do
    elseif strcmp(method_config.opts.order, 'biggest')
      w = candidates(:,3) - candidates(:,1) + 1;
      h = candidates(:,4) - candidates(:,2) + 1;
      areas = w .* h;
      [~,order] = sort(areas, 'descend');
      candidates = candidates(order,:);
      scores = scores(order,:);
    elseif strcmp(method_config.opts.order, 'smallest')
      w = candidates(:,3) - candidates(:,1) + 1;
      h = candidates(:,4) - candidates(:,2) + 1;
      areas = w .* h;
      [~,order] = sort(areas, 'ascend');
      candidates = candidates(order,:);
      scores = scores(order,:);
    elseif strcmp(method_config.opts.order, 'random')
      s = RandStream('mt19937ar','Seed',0);
      perm = randperm(s, size(candidates,1));
      candidates = candidates(perm,:);
      if numel(scores) > 0
        scores = scores(perm);
      end
    else
      [scores, argsort] = sort(scores, method_config.opts.order);
      candidates = candidates(argsort,:);
    end

    num_candidates = min(num_candidates, size(candidates, 1));
    candidates = candidates(1:num_candidates,:);
    if numel(scores) > 0
      scores = scores(1:num_candidates,:);
    end
  else
    error('this shouldn''t be used');
  end
end

author: Cloud-CV

repo: object-proposals


text:
% Haar transform on Healpix
% Y = HealpixHaarTransform(X)

function Y = HealpixHaarTransform(X)

% nest depth
n = length(X);
m = log2(n / 12) / 2; % = log4(n / 12)

% Haar wavelet basis
A = [ 1  1  1  1;
      1  1 -1 -1;
      1 -1  1 -1;
      1 -1 -1  1 ] / sqrt(4);
%A = [ 1  1  1  1;
%      1  1 -1 -1;
%      sqrt(2) -sqrt(2)  0  0;
%      0  0  sqrt(2) -sqrt(2) ] / sqrt(4);

Y = HealpixNestedLinearTrans(X, A, m);

author: mtex-toolbox

repo: mtex

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


text:
function out_ztimes = genas(cumu,xt,totbin,bin0,bin1)
    % Uses the GenAS algorithm to determine times of maximum Z values (as given by function AS) for a cumulative time curve
    %
    %  Syntax:     ztimes = genas(cumu,xt,totbin,bin0,bin1)
    %
    %  This Matlab function uses the GenAS algorithm of Habermann
    %  to determine times of maximum Z values (as given by function AS)
    %  for a cumulative time curve.
    %
    %  cumu is a histogram of events using a predefined bin length
    %  xt is the total time vector (in decimal years)
    %  bin0 is the cutoff at the beginning and bin1 at end of the analyses
    %  totbin is the total number of bins (including those with 0 z-values)
    %  ztimes is a vector with max-zvalues, its indexes give the bin number
    %  ------------------                                  R. Zuniga, 4/94

    global ztimes
    global sumx

    report_this_filefun();
    as=zeros(1,totbin);
    % if ~exist('sumx', 'var'); sumx = sum(cumu); end
    sumx = [sumx; sum(cumu)];
    sumx = max(sumx);
    par2 = sumx*0.1;
    %
    %
    for i = bin0+3:1:bin1-3          % calculate mean and z value for AS
        mean1 = mean(cumu(bin0:i));
        mean2 = mean(cumu(i+1:bin1));
        var1 = cov(cumu(bin0:i));
        var2 = cov(cumu(i+1:bin1));
        if mean1 && mean2 ~= 0
            as(i) = (mean1 - mean2)/(sqrt(var1/(i-bin0+1)+var2/(bin1-i)));
        else
            as(i) = 0;

        end     %if mean1

    end     % for i

    %   S = sprintf('bin0 %3d bin1 %3d i  %d',bin0, bin1, i);
    %   disp(S)
    %
    % check for threshold  (z = 1.96 -> 95%,  2.57  -> 99%)
    %as
    [xmax,ixs] = max(abs(as));
    if abs(as(ixs)) >= 2.57
        zmax = as(ixs);
        as = as*0 ;
        as(ixs) = zmax;

        ztimes(ixs) = as(ixs);      % form (vector) ztimes
        %  find(ztimes)

        xsum = cumsum(cumu);

        t1(1) = xt(ixs);
        t1(2) = xsum(ixs);
        t1p = [  t1(1)  t1(2); t1(1)   t1(2)+par2 ];
        plot(t1p(:,1),t1p(:,2),'k');
        set(gca,'NextPlot','add');

        S = sprintf('bin0 %d sig-Z at %d bin1 %d ',bin0, ixs, bin1);
        disp(S)

        ztimes = genas(cumu,xt,totbin,bin0,ixs);  %call genas again for both extremes
        ztimes = genas(cumu,xt,totbin,ixs,bin1);

    end     %if abs

    ztimes(1,totbin) = 0;    %   pad the end of ztimes
    as = as*0;
    out_ztimes = ztimes; %return a version that isn't the global

    %
    %  Plot the as(t)
    %
    %figure_w_normalized_uicontrolunits(2)
    %clf
    % orient tall
    % rect = [0.2,  0.20, 0.55, 0.75];
    % axes('position',rect)
    %set(gca,'NextPlot','add')
    %%plotyy(xt,cumu2,xt,as*10,'m')
    %% y2label('z-value')
    %% plot(xt,as*10,'+m')
    %% plot(xt,as*10,'m')
    % %text(0.70,0.5,'+: AS * 10','sc')
end

author: CelsoReyes


text:
function Seqs_SDC = Generate_SDC(Seqs_Long, rate)

Seqs_SDC = Seqs_Long;

parfor n = 1:length(Seqs_Long)
    Start_new = Seqs_Long(n).Start + ...
        rand*(Seqs_Long(n).Stop - Seqs_Long(n).Start)*(1-rate);

    Stop_new = Start_new + rate*(Seqs_Long(n).Stop - Seqs_Long(n).Start);

    Seqs_SDC(n).Start = Start_new;
    Seqs_SDC(n).Stop = Stop_new;

    ind = find(Seqs_SDC(n).Time>=Start_new && Seqs_SDC(n).Time<=Stop_new);
    Seqs_SDC(n).Time = Seqs_SDC(n).Time(ind);
    Seqs_SDC(n).Mark = Seqs_SDC(n).Mark(ind);
    if ~isempty(Seqs_SDC(n).Feature)
        Seqs_SDC(n).Feature = Seqs_SDC(n).Feature(ind,:);
    end
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
classdef DiffReactTests < matlab.unittest.TestCase

    properties (TestParameter)
%         file = {'testDiffReactHexagon', 'testDiffReactTorus', 'testDiffReactCube'}
        file = {'testDiffReactHexagon'}
        file3d = {'testDiffReactTorus', 'testDiffReactCube'}
        LHStype = {'DiffReactNeumann', 'DiffReactRobin'}
    end

    methods (Test, TestTags = {'DiffReact', '2D'})

        function testHexagon(testCase, file, LHStype)
            s   = testCase.createFEMparameters(file, LHStype);
            RHS = testCase.createRHS(s.mesh);
            fem = FEM.create(s);
            fem.computeLHS(0.1857);
            fem.computeVariables(RHS);
%             fem.print(filename)
            err = testCase.computeError(file, LHStype, fem);
            tol = 1e-6;
            testCase.verifyLessThanOrEqual(err, tol)
        end

    end

    methods (Test, TestTags = {'DiffReact', '3D'})

        function test3D(testCase, file3d)
            lhstype = 'DiffReactNeumann';
            s   = testCase.createFEMparameters(file3d, lhstype);
            RHS = testCase.createRHS(s.mesh);
            fem = FEM.create(s);
            fem.computeLHS(0.1857);
            fem.computeVariables(RHS);
%             fem.print(filename)
            err = testCase.computeError(file3d, lhstype, fem);
            tol = 1e-6;
            testCase.verifyLessThanOrEqual(err, tol)
        end

    end

    methods (Access = private)

        function error = computeError(testCase, file, LHStype, fem)
            file2load = append(file, '_', LHStype);
            cV = fem.variables.x;
            sV = load(file2load).x;
            error = norm(sV - cV)/norm(sV);
        end

        function name = loadFile(testCase, file, LHStype)
            name = append(file, '_', LHStype);
        end

        function s = createFEMparameters(testCase, file, LHStype)
            gidParams = testCase.createGiDparameters(file);
            s.dim       = gidParams.pdim;
            s.type      = gidParams.ptype;
            s.scale     = gidParams.scale;
            s.mesh      = gidParams.mesh;
            s.LHStype   = LHStype;
        end

        function gidParams = createGiDparameters(testCase, file)
            gidReader = FemInputReader_GiD();
            gidParams = gidReader.read(file);
        end

        function rhs = createRHS(testCase, mesh)
            M = testCase.computeM(mesh);
            u = testCase.createDisplacement(M);
            rhs = M*u;
        end

        function M = computeM(testCase, mesh)
            a.mesh    = mesh;
            a.fValues = zeros(mesh.nnodes, 1);
            f = P1Function(a);
            s.type         = 'MassMatrix';
            s.quadratureOrder     = 'QUADRATICMASS';
            s.mesh         = mesh;
            s.fun        = f;
            LHS = LHSintegrator.create(s);
            M = LHS.compute();
        end

        function u = createDisplacement(testCase, M)
            sizeM = size(M,1);
            sizeI = floor(sizeM/2);
            sizeJ = sizeM - sizeI;
            u = [ones(sizeI, 1); zeros(sizeJ, 1)];
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [s, cfg] = ft_statfun_indepsamplesF(cfg, dat, design)

% FT_STATFUN_INDEPSAMPLESF calculates the independent samples F-statistic on the
% biological data in dat (the dependent variable), using the information on the
% independent variable (ivar) in design.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option
%   cfg.statistic = 'ft_statfun_indepsamplesF'
%
% Configuration options
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computeprob    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or
% cfg.computeprob='yes'.
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail  = -1, 0, or 1, left, two-sided, or right (default=1)
%               cfg.tail in combination with cfg.computecritval='yes'
%               determines whether the critical value is computed at
%               quantile cfg.alpha (with cfg.tail=-1), at quantiles
%               cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%               quantile (1-cfg.alpha) (with cfg.tail=1).
%
% Design specification
%   cfg.ivar  = row number of the design that contains the labels of the conditions that must be
%               compared (default=1). The labels range from 1 to the number of conditions.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

% Copyright (C) 2006, Eric Maris
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% set the defaults
if ~isfield(cfg, 'computestat'),       cfg.computestat='yes';     end
if ~isfield(cfg, 'computecritval'),    cfg.computecritval='no';   end
if ~isfield(cfg, 'computeprob'),       cfg.computeprob='no';      end
if ~isfield(cfg, 'alpha'),             cfg.alpha=0.05;            end
if ~isfield(cfg, 'tail'),              cfg.tail=1;                end

% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
    ft_error('P-values can only be calculated if the test statistics are calculated.');
end
if isfield(cfg,'uvar') && ~isempty(cfg.uvar)
    ft_error('cfg.uvar should not exist for an independent samples statistic');
end

ncond=length(unique(design(cfg.ivar,:)));
nrepl=0;
for condindx=1:ncond
    nrepl=nrepl+length(find(design(cfg.ivar,:)==condindx));
end
if nrepl<size(design,2)
  ft_error('Invalid specification of the independent variable in the design array.');
end
if nrepl<=ncond
    ft_error('The must be more trials/subjects than levels of the independent variable.');
end
dfnum = ncond - 1;
dfdenom = nrepl - ncond;

nsmpls = size(dat,1);

if strcmp(cfg.computestat, 'yes')
  % compute the statistic
  nobspercell = zeros(1,ncond);
  avgs = zeros(nsmpls,ncond);
  pooledvar = zeros(nsmpls,1);
  for condindx=1:ncond
      sel=find(design(cfg.ivar,:)==condindx);
      nobspercell(condindx)=length(sel);
      avgs(:,condindx)=mean(dat(:,sel),2);
      pooledvar = pooledvar + nobspercell(condindx)*var(dat(:,sel),1,2);
  end
  pooledvar = pooledvar/dfdenom;
  globalavg = mean(dat,2);
  mseffect = ((avgs-repmat(globalavg,1,ncond)).^2)*nobspercell'./dfnum;
  s.stat = mseffect./pooledvar;
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  s.dfnum   = dfnum;
  s.dfdenom = dfdenom;
  if cfg.tail==-1
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a left tail critical value.');
  end
  if cfg.tail==0
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a two-sided critical value.');
  end
  if cfg.tail==1
    s.critval = finv(1-cfg.alpha,s.dfnum,s.dfdenom);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  s.dfnum   = dfnum;
  s.dfdenom = dfdenom;
  if cfg.tail==-1
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a left tail p-value.');
  end
  if cfg.tail==0
      ft_error('For an independent samples F-statistic, it does not make sense to calculate a two-sided p-value.');
  end
  if cfg.tail==1
    s.prob = 1-fcdf(s.stat,s.dfnum,s.dfdenom);
  end
end

author: spm


text:
function niederreiter2_test04 ( )

%*****************************************************************************80
%
%% NIEDERREITER2_TEST04 tests NIEDERREITER2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2008
%
%  Author:
%
%    John Burkardt
%
  dim_max = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'NIEDERREITER2_TEST04\n' );
  fprintf ( 1, '  NIEDERREITER2 computes the next element of\n' );
  fprintf ( 1, '  a Niederreiter quasirandom sequence using base 2.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we call NIEDERREITER2 repeatedly.\n' );

  for dim_num = 2 : dim_max

    seed = 0;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using dimension DIM_NUM =   %d\n', dim_num );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Seed  Seed   Niederreiter2\n' );
    fprintf ( 1, '  In    Out\n' );
    fprintf ( 1, '\n' );

    for i = 0 : 110
      seed_in = seed;
      [ r, seed ] = niederreiter2 ( dim_num, seed );
      seed_out = seed;
      if ( i <= 11 | 95 <= i )
        fprintf ( 1, '%3d  %3d  ', seed_in, seed_out );
        for j = 1 : dim_num
          fprintf ( 1, '%10f', r(j) );
        end
        fprintf ( 1, '\n' );
      elseif ( i == 12 )
        fprintf ( 1, '......................\n' );
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function finetuned_model_path = train_net_bbox_reg(...
    image_db_train, image_db_test, pooler, opts)

data_param                  = opts.data_param;
solver_param.test_iter      = data_param.test_iter;
solver_param.test_interval  = data_param.test_interval;
solver_param.max_iter       = opts.max_iter;

% ------------------------------------------------

work_dir  = opts.finetune_rst_dir;
timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
log_file  = fullfile(work_dir, 'output', ['mrcnn_regression_' timestamp, '.txt']);

mkdir_if_missing(work_dir);
mkdir_if_missing(fileparts(log_file));
diary(log_file);

fprintf('data_param:\n');
disp(data_param);
fprintf('solver_param:\n');
disp(solver_param);
fprintf('opts_param:\n');
disp(opts);
fprintf('pooler:\n');
disp(pooler);

data.feature_paths_train = image_db_train.feature_paths;
data.feature_paths_test  = image_db_test.feature_paths;

[data.fg_windows_train] = setup_data(image_db_train, data_param);
[data.fg_windows_test]  = setup_data(image_db_test,  data_param);

% data.fg_windows_train = setup_data(image_db_train, data_param);
% data.fg_windows_test  = setup_data(image_db_test,  data_param);

data.fg_windows_num_train = cell2mat(cellfun(@(x) size(x, 1), data.fg_windows_train, 'UniformOutput', false));
data.fg_windows_num_test  = cell2mat(cellfun(@(x) size(x, 1), data.fg_windows_test,  'UniformOutput', false));

fprintf('total fg_windows_num_train : %d\n', sum(data.fg_windows_num_train));
fprintf('total fg_windows_num_test  : %d\n', sum(data.fg_windows_num_test));

model.pooler = pooler;

current_dir = pwd;
cd(fileparts(opts.finetune_net_def_file))
% init caffe solver
[solver, iter_] = InitializeSolver(opts);

fprintf('Starting iteration %d\n', iter_);

% fix the random seed for repeatability
data_param_test                = data_param;
data_param_test.iter_per_batch = solver_param.test_iter;
data_param_test.random_scale   = 0; % for test use the best scale
data_param_test.nTimesMoreData = 10;

data.save_data_test_file = fullfile(work_dir, ...
    sprintf('data_test_file_batch%d_test_iter%d.mat', ...
    data_param_test.img_num_per_iter, data_param_test.iter_per_batch));

%% finetune
prev_rng = seed_rand();
mean_test_error = TestNet(solver, model, data_param_test, solver_param, data);

%min_test_error = TestNet(model, data_param_test, solver_param, data);

if isfield(opts, 'test_only') && opts.test_only
    finetuned_model_path = '';
    return;
end
last_finetuned_model_prev_path = '';

rng(prev_rng);
while(iter_ < solver_param.max_iter)
    mean_train_error = TrainNet(solver, model, data_param,      solver_param, data);
    mean_test_error  = TestNet( solver, model, data_param_test, solver_param, data);
    iter_            = solver.iter();
    ShowState(iter_, mean_train_error, mean_test_error);

    last_finetuned_model_path = fullfile(opts.finetune_rst_dir, sprintf('%s_iter_%d.caffemodel', opts.snapshot_prefix, iter_));
    assert(exist(last_finetuned_model_path,'file')>0);
    fprintf('last save as %s\n', last_finetuned_model_path);

    delete_prev_model_if_exist(last_finetuned_model_prev_path);
    last_finetuned_model_prev_path = last_finetuned_model_path;

    diary; diary; % flush diary
end
cd(current_dir);
finetuned_model_path = last_finetuned_model_path;
fprintf('Optimization is finished\n');
end

function [solver, iter_] = InitializeSolver(opts)
solver = caffe.Solver(opts.finetune_net_def_file);
iter_   = 1;

if isfield(opts,'solver_state_file') && ~isempty(opts.solver_state_file)
    solver.restore(opts.solver_state_file);
    [directory,filename] = fileparts(opts.solver_state_file);
    model_file = [directory,filesep, filename, '.caffemodel'];
    solver.net.copy_from(model_file);
    solver.test_nets(1).copy_from(model_file);
    iter_   = solver.iter();
    rng('shuffle')
else
    solver.net.copy_from(opts.net_file);
    solver.test_nets(1).copy_from(opts.net_file);
end
end

function ShowState(iter, train_error, test_error)
fprintf('\n------------------------- Iteration %d -------------------------\n', iter);
fprintf('Error - Training : %.4f - Testing : %.4f\n', train_error, test_error);
end

function delete_prev_model_if_exist(finetuned_model_prev_path)
if exist(finetuned_model_prev_path, 'file')
	delete(finetuned_model_prev_path);
end
[directory,filename] = fileparts(finetuned_model_prev_path);
finetuned_model_prev_path_solver_state = [directory,filesep, filename,'.solverstate'];
if exist(finetuned_model_prev_path_solver_state, 'file')
    delete(finetuned_model_prev_path_solver_state);
end
end

function fg_windows = setup_data(image_db, data_param)

image_paths = image_db.image_paths;
all_regions = image_db.all_regions;
all_bbox_gt = image_db.all_bbox_gt;

num_elems  = length(image_paths);
fg_windows = cell(num_elems, 1);

num_classes  = data_param.num_classes;
fg_threshold = data_param.fg_threshold;

for i = 1:num_elems
    fg_windows_this = setup_img_data(image_paths{i}, all_bbox_gt{i}, all_regions{i}, num_classes, fg_threshold);
    fg_windows{i}   = [ones(size(fg_windows_this,1),1,'single')*i, fg_windows_this];
end

end

function bbox_out = setup_img_data(image_file, bbox_gt, bbox_proposals, num_classes, low_ov_threshold)

bbox_gt_label = bbox_gt(:,5);

bbox_all    = [bbox_gt(:,1:4); bbox_proposals(:,1:4)];
overlap     = zeros(size(bbox_all,1), num_classes,     'single');
bbox_target = zeros(size(bbox_all,1), num_classes * 4, 'single');

for class_idx = 1:num_classes
    if any(bbox_gt_label == class_idx)
        bbox_gt_class = bbox_gt(bbox_gt_label == class_idx,1:4);
        overlap_class = boxoverlap(bbox_all(:,1:4), bbox_gt_class(:,1:4));
        [overlap(:, class_idx), gt_index] = max(overlap_class, [], 2);
        does_overlap = overlap(:, class_idx) > 0.001;
        bbox_target(does_overlap, (class_idx-1)*4 + (1:4)) = bbox_gt_class(gt_index(does_overlap), 1:4);
    end
end

[max_overlap, label] = max(overlap, [], 2);
keep = max_overlap >= low_ov_threshold;

bbox_all         = bbox_all(   keep,:);
overlap          = overlap(    keep,:);
bbox_target      = bbox_target(keep,:);
bbox_target_vals = encode_bbox_targets_to_reg_vals(bbox_all, bbox_target);

max_overlap      = max_overlap(keep);
label            = label(keep);

bbox_out         = single([label, max_overlap, bbox_all, overlap, bbox_target_vals, bbox_target]);
end

function [data_test, label_test, windows_test] = create_test_batch(model, data_param_test, data )
try
    ld = load(data.save_data_test_file);
    data_test    = ld.data_test;
    label_test   = ld.label_test;
    windows_test = ld.windows_test;
catch
    windows_test = sample_window_data(data.fg_windows_test, data.fg_windows_num_test, data_param_test);
    data_test    = prepare_batch_feat_input_data(model.pooler, data.feature_paths_test, windows_test, data_param_test);
    label_test   = prepare_batch_target_data(windows_test, data_param_test);

    save(data.save_data_test_file, 'data_test', 'label_test', 'windows_test', '-v7.3');
end
end

function [windows, select_img_idx] = sample_window_data(fg_windows, fg_windows_num, data_param)

nTimesMoreData = 10;

if isfield(data_param,'nTimesMoreData')
    nTimesMoreData = data_param.nTimesMoreData;
end

fg_num_total = data_param.img_num_per_iter * data_param.iter_per_batch;
num_img      = length(fg_windows);

% random perm image and get top n image with total more than patch_num
% windows
img_idx           = randperm(num_img);
fg_windows_num    = fg_windows_num(img_idx);
fg_windows_cumsum = cumsum(fg_windows_num);
img_idx_end       = find(fg_windows_cumsum > (fg_num_total * nTimesMoreData), 1, 'first');

if isempty( img_idx_end ), img_idx_end = length(img_idx); end

select_img_idx = img_idx(1:img_idx_end);
fg_windows_sel = cell2mat(fg_windows(select_img_idx));

% random perm all windows, and drop redundant data
window_idx     = randperm(size(fg_windows_sel, 1));
fg_indices     = mod((1:fg_num_total) - 1, length(window_idx)) + 1;
fg_windows_sel = fg_windows_sel(window_idx(fg_indices), :);

random_seed = randi(10^6, length(select_img_idx),1);
for i = 1:length(select_img_idx), rng(random_seed(i), 'twister'); end

windows     = fg_windows_sel(randperm(size(fg_windows_sel, 1)),:);
end

function label = prepare_batch_target_data(windows, data_param)
num_classes  = data_param.num_classes;
num_reg_vals = num_classes * 4;
fg_threshold = data_param.fg_threshold;
mask         = windows(:,7+(1:num_classes)) >= fg_threshold;
mask         = repmat( mask, [4, 1]);
mask         = reshape(mask, [size(windows,1), num_classes*4]);
label        = [single(mask), windows(:,(7+num_classes)+(1:num_reg_vals))]';
end

function [mean_ave_recall, mean_test_error] = TestNet(solver, model, data_param, solver_param, data)
fprintf('Testing: ');
th = tic;

[data_test, label_test, windows_test] = create_test_batch(model, data_param, data );

assert(isnumeric(data_test));
assert(isnumeric(label_test));
assert(size(data_test,2)    == solver_param.test_iter * data_param.img_num_per_iter);
assert(size(data_test,1)    == data_param.feat_dim);
assert(size(label_test,2)   == solver_param.test_iter * data_param.img_num_per_iter);
assert(size(windows_test,1) == solver_param.test_iter * data_param.img_num_per_iter);

% windows_test = [fg_img_idx, fg_label, fg_ov, fg_window, fg_overlap_all, fg_bbox_target_vals, fg_bbox_target];
mini_batch_size = data_param.img_num_per_iter;

fprintf('elapsed time %.2f + ', toc(th));
th = tic;
test_error = nan(solver_param.test_iter, 1);
test_pred  = cell(solver_param.test_iter,1);
% fprintf('Testing: ');
for it = 1:solver_param.test_iter
    start_idx = (it-1) * mini_batch_size + 1;
    stop_idx  = it * mini_batch_size;
    assert(stop_idx <= size(data_test,2));

    mini_batch = prepare_mini_batch(data_test( :,start_idx:stop_idx), ...
                               label_test(:,start_idx:stop_idx), data_param);
    caffe_set_input(solver.test_nets(1), mini_batch);
    forward_prefilled(solver.test_nets(1));

    test_error(it) = solver.test_nets(1).blobs('error').get_data();
    test_pred{it}  = solver.test_nets(1).blobs('predictions').get_data();
end

mean_test_error = mean(test_error);

num_classes      = data_param.num_classes;
bbox_init        = windows_test(:,4:7);
all_bbox_pred    = repmat(bbox_init, [1, num_classes]);
all_bbox_targets = windows_test(:, (7+num_classes+num_classes*4) + (1:(num_classes*4)));
all_mask         = label_test(1:(num_classes*4),:)';

fprintf('Initial Candidate Boxes');
mean_ave_recall_init  = compute_mAR_of_bbox_reg(bbox_init, all_bbox_pred, all_bbox_targets, all_mask, num_classes);

all_reg_values  = cell2mat(test_pred(:)')';
all_bbox_pred   = decode_reg_vals_to_bbox_targets(bbox_init, all_reg_values);
fprintf('After Regression Candidate Boxes');
mean_ave_recall = compute_mAR_of_bbox_reg(bbox_init, all_bbox_pred, all_bbox_targets, all_mask, num_classes);
fprintf('%.2f sec --- error %.4f mAR %.5f --> %.5f\n', toc(th), mean_test_error, mean_ave_recall_init, mean_ave_recall);
mean_ave_recall = 1 - mean_ave_recall;

end

function mean_train_error = TrainNet(solver, model, data_param, solver_param, data)

test_interval   = solver_param.test_interval;
num_batches     = data_param.iter_per_batch;
num_epochs      = test_interval / num_batches;
mini_batch_size = data_param.img_num_per_iter;

train_count = 0;
train_error_sum = 0;
for epoch = 1:num_epochs
    fprintf('Training: ');
    th = tic;

    windows_train = sample_window_data(data.fg_windows_train,   data.fg_windows_num_train, data_param);
    data_train    = prepare_batch_feat_input_data(model.pooler, data.feature_paths_train, windows_train, data_param);
    label_train   = prepare_batch_target_data(windows_train, data_param);

    assert(isnumeric(data_train));
    assert(isnumeric(label_train));
    assert(size(data_train,2)  == num_batches * data_param.img_num_per_iter);
    assert(size(data_train,1)  == data_param.feat_dim);
    assert(size(label_train,2) == num_batches * data_param.img_num_per_iter);

    fprintf('elapsed time %.2f + ', toc(th));

    th = tic;
    train_error_sum_this = 0;
    for it = 1:num_batches
        start_idx = (it-1) * mini_batch_size + 1;
        stop_idx  = it * mini_batch_size;
        assert(stop_idx <= size(data_train,2));
        caffe_set_input(solver.net, ...
            prepare_mini_batch(data_train( :,start_idx:stop_idx), ...
                               label_train(:,start_idx:stop_idx), data_param));
        solver.step(1);

        results = caffe_get_output(solver.net);

        train_error_sum_this = train_error_sum_this + results{1}(1);
    end
    train_error_sum = train_error_sum + train_error_sum_this;
    train_count = train_count + num_batches;
    data_train = []; label_train = [];
    fprintf('%.2f sec --- error[total] %.4f - error[this] %.4f\n', toc(th), train_error_sum/train_count, train_error_sum_this/num_batches);
end
mean_train_error = train_error_sum/train_count;
end

function mini_batch = prepare_mini_batch(data, labels, data_param)
data   = prepare_data(    data, data_param);
labels = prepare_labels(labels, data_param);
mini_batch = [data, labels];
end

function labels = prepare_labels(labels, data_param)

num_samples = size(labels,2);

if ~isfield(data_param, 'labels_split_div')
    labels_split_div = size(labels,1);
else
    labels_split_div = data_param.labels_split_div;
    assert(sum(labels_split_div) == size(labels,1));
end

labels = mat2cell(labels, labels_split_div, size(labels,2))';

if isfield(data_param, 'reshape_labels')
    reshape_labels = data_param.reshape_labels;
    num_labels_types = length(labels);
    assert(length(labels) == length(reshape_labels));
    for i = 1:num_labels_types
        assert(length(reshape_labels{i}) == 4 || length(reshape_labels{i}) == 3);
        if length(reshape_labels{i}) == 3
            reshape_labels{i}(4) = num_samples;
        end
        labels{i} = reshape(labels{i}, reshape_labels{i});
    end
end

end

function data = prepare_data(data, data_param)
num_samples = size(data,2);

if ~isfield(data_param, 'data_split_div')
    data_split_div = size(data,1);
else
    data_split_div = data_param.data_split_div;
    assert(sum(data_split_div) == size(data,1));
end

data = mat2cell(data, data_split_div, size(data,2))';

if isfield(data_param, 'reshape_data')
    reshape_data = data_param.reshape_data;
    num_data_types = length(data);
    assert(length(data) == length(reshape_data));
    for i = 1:num_data_types
        assert(length(reshape_data{i}) == 4 || length(reshape_data{i}) == 3);
        if length(reshape_data{i}) == 3
            reshape_data{i}(4) = num_samples;
        end
        data{i} = reshape(data{i}, reshape_data{i});
    end
end
end

function [mAR, AR] = compute_mAR_of_bbox_reg(bbox_init, all_bbox_pred, all_bbox_targets, all_mask, num_classes)

AR = zeros(1,num_classes);
fprintf('\nAveRecall: ')
for c = 1:num_classes
    bbox_pred    = all_bbox_pred(:,    (c-1)*4 + (1:4));
    bbox_targets = all_bbox_targets(:, (c-1)*4 + (1:4));
    mask         = all_mask(:,         (c-1)*4 + (1:4));
    is_not_gt    = ~all(bbox_init == bbox_targets, 2);
    bbox_pred    = bbox_pred(   is_not_gt, :);
    bbox_targets = bbox_targets(is_not_gt, :);
    mask         = mask(        is_not_gt, :);
    mask         = all(mask,2);
    bbox_pred    = bbox_pred(   mask,:);
    bbox_targets = bbox_targets(mask,:);
    AR(c)        = compute_ave_recall_of_bbox( bbox_pred, bbox_targets );
    fprintf('[%d]=%.2f ',c, AR(c));
end

fprintf('\n')
mAR = mean(AR);

end

author: gidariss

repo: mrcnn-object-detection


text:
%% Simplicial Complex in 3D
%
% We describe general idea of the data structure generated for three
% dimensional triangular grids, suppose $N,NE,NF,NT$ represents the number
% of vertices, edges, faces and elements resprectively.

%% Construct Data Structure
% elem is $NT\times 4$ matrix with $NT$ be the number of the tetrahedrons.
% elem(i,:) contain the global index of the vertex of the i-th tetrahedron.
% we sort the vertices of elem such
% that:elem(i,1)<elem(i,2)<elem(i,3)<elem(i,4). In this way, the volume is
% not always positive, because we sacrifice the righthand orientation of tet.
% We also sort vertices of face such that:face(j,1)<face(j,2)<face(j,3),and
% the vertices of edge such that: edge(k,1)<edge(k,2).
%
% The local ordering of face in a tetrahedron is locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3].
% The local ordering of edge in a face is locEdge = [2 3;1 3;1 2].
% Now the ordering of local indices is always consistent with that of global
% indices.
%
% In the following example, elem2face, a $NT\times 4$ matrix, is the
% elementwise pointer from elem to face indicies. elem(i,:) stores the
% global indicies of faces for tetrahedron i. elem2edge, a $NT\times 6$ matrix, is the
% elementwise pointer from elem to edge indicies. elem(j,:) stores the
% global indicies of edges for tetrahedron j.

%% Example
% Generate

elem = [1 4 5 8;1 4 5 7];
node = [1,0,0;1,1,1;1,-1,-1; 0,1,0; 0,0,0;1,1,-1; 0,0,1;-1,-1,-1];
[elem2edge,~,edge] = dof3edge(elem);
[elem2face,~,face] = dof3RT0(elem);
figure(1);clf;
showmesh3(node,elem);
view(-20,20);
findnode3(node,[1,4,5,7,8]);
findelem3(node,elem);
display(edge); display(elem2edge);
display(face); display(elem2face);

%% The directions.
%% Direction of face.
% While we sort the elem, face and edge in the increasing order of
% vertices. We determine the direction of a face [i j k] by the righthand
% rule. The direction of face is not always outword normal direction of
% a certain tetrahedron. We denote the direction as +1 if the direction
% of a face is the same with the normal direction in a certain elem, and -1
% otherwise. Then for locface of a tetrahedron, the directions should be
% sign(volume)*[+1 -1 +1 -1].
%% Direction of edge.
% Similarly, we determine the direction of a edge [i j] is from the node with
% smaller global index to bigger one(that is from i to j). The direction of
% edge is not always consistent with the counter clockwise direction in a
% face. When the edge direction is the same with the local counter
% clockwise direction, we set the direction as +1, otherwise -1. Then for
% locedge of a face, the directions should be [+1 -1 +1].
%% Direction of the basis.
% According to our sorted elem and face, the RT0 basis on face [i j k].
% Defined as $\phi_{i,j,k} = 2(\lambda_i \bigtriangledown \lambda_j \times
% \bigtriangledown \lambda_k+
% \lambda_j \bigtriangledown \lambda_k \times \bigtriangledown
% \lambda_i+\lambda_k \bigtriangledown \lambda_i \times \bigtriangledown
% \lambda_j)$ The direction is always consistent with the face direction.
%
%% Inward normal direction.
% Dlambda(t,:,1) is the gradient of lambda of the 1st index of the
% t-th tet. The direction of Dlambda is always inward normal direction.
[Dlambda, ~] = gradbasis3(node,elem(2,:));
display(Dlambda);

repo: Aero-matlab


text:
function y = mnorm_pdf(x,mu,S)
%MNORM_PDF Multivariate-Normal probability density function (pdf).
%
%   Description
%   Y = MNORM_PDF(X,MU,SIGMA) Returns the multivariate-normal
%   pdf with mean, MU, and covariance matrix, SIGMA, at the values in X.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.
%
% Copyright (c) 1998-2005 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3,
  S = 1;
end

if nargin < 2;
  mu = 0;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

[n,m]=size(x);
xmu=x-repmat(mu,n,1);
if ~issparse(S)
    % Use Cholesky decomposition, since it is faster and
    % numerically more stable.
    L=chol(S,'lower');
    y=zeros(n,1);
    for i1=1:n
      b=L\xmu(i1,:)';
      y(i1)=-b'*b;
    end
    y=exp(.5*y-sum(log(diag(L)))-.5*m*log(2*pi));
else
    LD = ldlchol(S);
    y=zeros(n,1);
    for i1=1:n
        xmui1=xmu(i1,:)';
      y(i1)=-xmui1'*ldlsolve(LD,xmui1);
    end
    y=exp(0.5*(y-sum(log(diag(LD)))-m*log(2*pi)));
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [theta] = normalEqn(X, y)
%NORMALEQN Computes the closed-form solution to linear regression
%   NORMALEQN(X,y) computes the closed-form solution to linear
%   regression using the normal equations.

theta=pinv(X'*X)*X'*y;

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the code to compute the closed form solution
%               to linear regression and put the result in theta.
%

% ---------------------- Sample Solution ----------------------

% -------------------------------------------------------------

% ============================================================

end

author: loserChen

repo: Coursera-MachineLearning


text:
classdef nnlayers < nncheckder
  properties (TestParameter)
    conserveMemory = {false, true}
  end
  properties
    currentConserveMemory
  end

  methods (Test)
    function testLayers(test, conserveMemory)
      test.currentConserveMemory = conserveMemory ;

      % use Params for all inputs so we can choose their values now
      x = Param('value', randn(7, 7, 2, 5, test.currentDataType)) ;
      w = Param('value', randn(3, 3, 2, 3, test.currentDataType)) ;
      b = Param('value', randn(3, 1, test.currentDataType)) ;
      labels = Param('value', ones(5, 1, test.currentDataType)) ;
      Layer.workspaceNames() ;

      % test several layers and syntaxes

      do(test, vl_nnrelu(x)) ;

      do(test, vl_nnconv(x, w, b)) ;

      do(test, vl_nnconv(x, w, b, 'stride', 3, 'pad', 2)) ;

      do(test, vl_nnconvt(x, permute(w, [1 2 4 3]), b)) ;

      do(test, vl_nnconvt(x, permute(w, [1 2 4 3]), b, 'upsample', 3, 'crop', 2)) ;

      do(test, vl_nnpool(x, 2)) ;

      do(test, vl_nnpool(x, [2, 2], 'stride', 2, 'pad', 1)) ;

      do(test, vl_nnloss(x, labels, 'loss', 'classerror'), labels) ;

      % dropout is composed of 2 parts: the mask generator, and the dropout
      % mask applier. run derivative check with fixed mask.
      rate = 0.1 ;
      dropout = vl_nndropout(x, 'rate', rate) ;
      mask = dropout.inputs{2} ;

      test.verifyInstanceOf(mask, 'Layer') ;
      test.verifyEqual(mask.func, @vl_nnmask) ;

      dropout.inputs{2} = vl_nnmask(x.value, rate) ;  % make it constant
      do(test, dropout) ;

      % bnorm params are single
      if strcmp(test.currentDataType, 'single')
        % batch-norm needs special handling
        bnorm = vl_nnbnorm(x) ;

        % replicate gain/bias for all output channels (ordinarily done by
        % the solver on first update)
        bnorm.inputs{2}.value = bnorm.inputs{2}.value([1; 1]) ;  % gain
        bnorm.inputs{3}.value = bnorm.inputs{3}.value([1; 1]) ;  % bias

        % ignore the 4th parameter (moments), since it is not updated by
        % gradient descent but by a moving average
        ignore = bnorm.inputs{4} ;

        do(test, bnorm, ignore) ;
      end
    end

    function testMath(test, conserveMemory)
      test.currentConserveMemory = conserveMemory ;

      % use Params for all inputs so we can choose their values now
      a = Param('value', randn(3, 3, test.currentDataType) + 0.1 * eye(3,3)) ;  % matrix
      b = Param('value', randn(3, 3, test.currentDataType) + 0.1 * eye(3,3)) ;  % matrix
      c = Param('value', 2 * ones(1, 1, test.currentDataType)) ;  % scalar
      d = Param('value', randn(3, 1, test.currentDataType)) ;  % vector
      e = Param('value', rand(3, 3, test.currentDataType) + 1e-3 * ones(3,3)) ;  % non-negative matrix
      f = Param('value', rand(3, 3, test.currentDataType) * 2 - 1) ;  % matrix in -1..1
      Layer.workspaceNames() ;

      % test several operations

      % weighted sums
      do(test, a + b) ;
      do(test, 10 * a) ;
      do(test, a + 2 * b - c) ;  % collected arguments in a single wsum

      % matrix
      do(test, a * b) ;
      do(test, a') ;
      do(test, inv(a)) ;
      do(test, a / b) ;
      do(test, a \ b) ;

      % binary with expanded dimensions
      do(test, a .* d) ;
      do(test, a ./ d) ;
      do(test, a .^ c, [], 1e-6 * test.range, 1e-3) ;  % higher tolerance

      do(test, atan2(a, b)) ;

      % matrix ops should deal with scalars gracefully
      do(test, a * c) ;
      do(test, a / c) ;

      % unary
      do(test, sqrt(e)) ;
      do(test, sin(a)) ;
      do(test, cos(a)) ;
      do(test, tan(a)) ;
      do(test, asin(f)) ;
      do(test, acos(f)) ;
      do(test, atan(a)) ;

      % sorting is a kind of math
      do(test, sort(a)) ;
    end

    function testConv(test, conserveMemory)
      test.currentConserveMemory = conserveMemory ;

      % extra conv tests
      if strcmp(test.currentDataType, 'double'), return, end

      x = Param('value', randn(7, 7, 2, 5, test.currentDataType)) ;

      % 'size' syntax
      do(test, vl_nnconv(x, 'size', [7, 7, 2, 5])) ;

      % bias
      layer = vl_nnconv(x, 'size', [7, 7, 2, 5], 'hasBias', false) ;
      do(test, layer) ;
      test.verifyEmpty(layer.inputs{3}) ;

      % Param learning arguments
      layer = vl_nnconv(x, 'size', [7, 7, 2, 5], ...
          'learningRate', [1, 2], 'weightDecay', [3, 4]) ;
      do(test, layer) ;
      test.eq(layer.inputs{2}.learningRate, 1) ;
      test.eq(layer.inputs{3}.learningRate, 2) ;
      test.eq(layer.inputs{2}.weightDecay, 3) ;
      test.eq(layer.inputs{3}.weightDecay, 4) ;
    end
  end

  methods
    function do(test, output, varargin)
      % show layer for debugging
      display(output) ;

      % use parent class's derivative check (defined in nncheckder)
      test.checkDer(output, test.currentConserveMemory, varargin{:}) ;
    end
  end
end

repo: MFT


text:
function p3 = mpower(p1,p2)
%PREAL/MPOWER Overloaded mpower (^) operator for class preal.
% Note: preal/mpower is the same as preal/power.

global useUnitsFlag

if ~(useUnitsFlag) % If physunits is disabled...
    p3=mpower(double(p1),double(p2)); % ... treat as double.
    return
end

p1=preal(p1);
p2=preal(p2);
tol=0.001;
if isscalar(p1)
    p1=preal(p1.value*ones(size(p2)),p1.units);
elseif isscalar(p2)
    p2=preal(p2.value*ones(size(p1)),p2.units);
elseif ndims(p1)==ndims(p2)&size(p1)==size(p2)
    % do nothing
else
    error('Array dimensions must agree')
end
p3=preal(ones(size(p1)));
for k=1:numel(p1)
    if any(abs(p2(k).units)>tol)
        error('Exponent must be non-dimensional (pure) number')
    end
    p3(k).value=p1(k).value^p2(k).value;
    p3(k).units=p1(k).units*p2(k).value;
end

author: Sable

repo: mcbench-benchmarks


text:
function x = cvx_accept_concave( x )
global cvx___
if isa( x, 'cvx' ),
    t = cvx_vexity( x ) < 0;
    if any( t( : ) ),
        prob = cvx___.problems(end).self;
        if all( t ),
            src = x;
            dst = newtemp( prob, size( src ) );
            x = dst;
        else
            src = x( t );
            dst = newtemp( prob, size( src ) );
            x( t ) = dst;
        end
        newcnstr( prob, src(:), dst(:), '>=' );
    end
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function test_pull1412

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_heartrate

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/pull1412'));

%%
% this corresponds to the preprocessed dataset 006_3013065.02_rest1 from bug3433

load datappg

cfg = [];
cfg.channel = 'HR';
cfg.threshold = 0.7;
cfg.method = 'findpeaks';
% cfg.method = 'pantompkin'; this does not work very well on the PPG data
heartrate0 = ft_heartrate(cfg, data);

figure
plot(heartrate0.time{1}, heartrate0.trial{1}(1,:), '-');

%%
% this corresponds to the ECG channel from ArtifactMEG.ds as documented on
% http://www.fieldtriptoolbox.org/example/use_independent_component_analysis_ica_to_remove_ecg_artifacts/

load dataecg

cfg = [];
cfg.channel = 'ECG';
cfg.threshold = 1.2;
cfg.method = 'findpeaks';
cfg.flipsignal = 'no';
heartrate1 = ft_heartrate(cfg, data);

cfg = [];
cfg.channel = 'ECG';
cfg.method = 'pantompkin';
heartrate2 = ft_heartrate(cfg, data);

figure
plot(heartrate1.time{1}, heartrate1.trial{1}(1,:), 'b-');
hold on
plot(heartrate2.time{1}, heartrate2.trial{1}(1,:), 'rx');
legend({'findpeaks', 'pantompkin'});


text:
%----STORED STREAMING CREDENTIALS----%
my_credentials = loadplotlycredentials;
try
    my_stream_token = my_credentials.stream_ids{1};
catch
    fprintf(['\nOops - No stream_keys found! please run: >>saveplotlycredentials(',...
        ' ''username'',''api_key'',''stream_key).'' \n',...
        'Your stream key(s) can be found online at: https://plot.ly or contact chuck@plot.ly',...
        'for more information.\n\n']);
    return
end

%----SETUP-----%

p = plotlyfig('visible','off');
p.data{1}.x = [];
p.data{1}.y = [];
p.data{1}.type = 'scatter';
p.data{1}.stream.token = my_stream_token;
p.data{1}.stream.maxpoints = 30;
p.PlotOptions.Strip = false;
p.PlotOptions.FileName = 'stream_test';
p.PlotOptions.FileOpt = 'overwrite';

%----PLOTLY-----%

p.plotly;

%----CREATE A PLOTLY STREAM OBJECT----%

ps = plotlystream(my_stream_token);

%----OPEN THE STREAM----%

ps.open();

%----WRITE TO THE STREAM----%

for i = 1:2000
    mydata.x = i;
    mydata.y = rand;
    ps.write(mydata);
    %take a breath
    pause(0.05);
end

%----CLOSE THE STREAM----%
ps.close;

author: goGPS-Project


text:
function rx = rxLoadMrVistaAlignment(rx,mrSessPath)
%
% rx = rxLoadMrVistaAlignment([rx],[mrSessPath])
%
% Load a mrVista alignment from a mrSESSION.mat file
% into mrRx.
%
% ras 03/05.
if notDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if notDefined('mrSessPath')
    mrSessPath = fullfile(pwd,'mrSESSION.mat');
end

if ~exist(mrSessPath,'file')
    msg = sprintf('%s not found.',mrSessPath);
    myErrorDlg(msg);
end

h = msgbox('Loading mrVista Alignment...');

load(mrSessPath);

newXform = mrSESSION.alignment;

% flip to (x,y,z) instead of (y,x,z):
newXform(:,[1 2]) = newXform(:,[2 1]);
newXform([1 2],:) = newXform([2 1],:);

rx = rxSetXform(rx,newXform,0);

rxStore(rx,'mrVista Alignment');

close(h);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%     Steffen Urban email: steffen.urban@kit.edu
%     Copyright (C) 2014  Steffen Urban
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

% 04.03.2014 by Steffen Urban
% this is a modified file from
% Davide Scaramuzzas Toolbox OcamCalib
% filename: get_checkerboard_corners.m
% Code was changed to force subpixel corners

function [callBack, Xp_abs, Yp_abs] = get_checkerboard_cornersUrban(kk, use_corner_find, calib_data)

fclose('all');
if exist('autoCornerFinder/cToMatlab/cornerInfo.txt','file')
    delete('autoCornerFinder/cToMatlab/cornerInfo.txt');
end
if exist('autoCornerFinder/cToMatlab/cornersX.txt','file')
    delete('autoCornerFinder/cToMatlab/cornersX.txt');
end
if exist('autoCornerFinder/cToMatlab/cornersY.txt','file')
    delete('autoCornerFinder/cToMatlab/cornersY.txt');
end
if exist('autoCornerFinder/cToMatlab/error.txt','file')
    delete('autoCornerFinder/cToMatlab/error.txt');
end

%INITIALIZATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\nProcessing image %s...',calib_data.L{kk});

I = imread(calib_data.L{kk});

Xp_abs = 0;
Yp_abs = 0;

%Tell the automatic corner extractor, which image file to process
fid = fopen('./autoCornerFinder/pictures.txt','w');
fprintf(fid,'../%s',calib_data.L{kk});
fclose(fid);

%Call the automatic corner extraction algorithm
%Width and height in the algorithm are defined differently than
%in this toolbox: Its the number of inernal corners instead of
%internal quadrangles.
%-m specifies the number of corners the algorithm has to find, before it
%terminates
%If callback = -1   ->An error occured, automatic corner finding is aborted
%            =  0   ->Not enough corners have been found, add some manually
%            =  1   ->Enough corners have been found for calibration

%Visualization turned OFF
cd autoCornerFinder;

callString = (['FindCorners.exe -w ' num2str(calib_data.n_sq_x+1) ' -h ' num2str(calib_data.n_sq_y+1) ' -m ' num2str((calib_data.n_sq_x+1) * (calib_data.n_sq_y+1)) ' pictures.txt']);

if ~ispc %if not Windows
    callString = ['./' callString];
end
%Visualization turned ON
%callString = (['cd autoCornerFinder & FindCornersVisual.exe -w ' num2str(n_sq_x+1) ' -h ' num2str(n_sq_y+1) ' -m ' num2str((n_sq_x+1) * (n_sq_y+1)) ' pictures.txt']);

%Visualization turned ON and Saving of the images turned ON
%WARNING: Does somehow not work under Windows 2000...
%callString = (['cd autoCornerFinder & FindCornersVisualSave.exe -w ' num2str(n_sq_x+1) ' -h ' num2str(n_sq_y+1) ' -m ' num2str((n_sq_x+1) * (n_sq_y+1)) ' pictures.txt']);

callBack   = system(callString);
cd ..
%Do error checking
if callBack == 0
    %Display the error message
    fprintf(1,'Image omitted -- Not all corners found\n');
    return;
end

if callBack ~= 1
    %Display the error message
    filename = 'autoCornerFinder/cToMatlab/error.txt';
    fid = fopen(filename, 'r');
    line = fgetl(fid);
    fclose(fid);
    fprintf(1,'Image omitted -- During corner finding an error occured: %s \n',line);
    return;
end

%Open the corner size information file
filename = 'autoCornerFinder/cToMatlab/cornerInfo.txt';
fid = fopen(filename, 'r');
cornerInfo = fscanf(fid, '%g %g', [1 2]);
fclose(fid);

%Open the files with the found corners
filename = 'autoCornerFinder/cToMatlab/cornersX.txt';
fid = fopen(filename, 'r');
cornersX = fscanf(fid, '%g %g', [cornerInfo(2) cornerInfo(1)]);
cornersX = cornersX';
fclose(fid);

filename = 'autoCornerFinder/cToMatlab/cornersY.txt';
fid = fopen(filename, 'r');
cornersY = fscanf(fid, '%g %g', [cornerInfo(2) cornerInfo(1)]);
cornersY = cornersY';
fclose(fid);

%VARIABLE DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numCorners = (calib_data.n_sq_x+1)*(calib_data.n_sq_y+1);
numCornerThreshold = (calib_data.n_sq_x+1)*(calib_data.n_sq_y+1);
numOfFoundCorners = 0;
cornerNumber = 0;
cornersNumbering = -1*ones(cornerInfo(1),cornerInfo(2));
startingCorner = [];
nextFoundCorner = [-1,-1];
deltaCols = 0;
startingCornerSet = false;
%Index of the corner with the smallest index
min_i = [];
min_j = [];
%Used for image zoom in
[size1, size2] = size(I);
min_x = max(size1, size2) + 1;
max_x = 0;
min_y = max(size1, size2) + 1;
max_y = 0;

%INPUT CORNER VALUE AND MATRIX SIZE ADAPTATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%If there are entire rows of zeros at the end of the file, this means that
%the cornerfinder is not entirely shure about the location of the board.
%Then eventually additional corners are needed in the front.
[iSave, dummy] = find(cornersX >= 0);
iSave = max(iSave);
%Add at least one row! This is because the columns are maybe ambiguous,
%and then we don't know whether they belong to the current row or to the
%next one...
if (cornerInfo(1) - iSave >= 0)
    cornersX = [-1* ones(cornerInfo(1) - iSave + 1, cornerInfo(2)); cornersX];
    cornersY = [-1* ones(cornerInfo(1) - iSave + 1, cornerInfo(2)); cornersY];
end
[rows, cols] = size(cornersX);

%Add one pixel to every non "-1" value, since Matlab starts numbering
%at one, whereas c++ starts at zero.
flagStart = true;
for i = 1:1:rows
    for j = 1:1:cols
        %Define the starting corner as the first found corner
        %which has a neighbor corner which was also found
        if j ~= cols && flagStart == true
            if cornersX(i,j) >= 0 && cornersX(i,j+1) >= 0
                startingCorner = [i,j];
                flagStart = false;
            end
        end
        if cornersX(i,j) >= 0
            cornersX(i,j) = cornersX(i,j) + 1;
            %Count the number of found corners
            numOfFoundCorners = numOfFoundCorners + 1;
            %Determine the minimal and maximal cordinates of all
            %found corners. ->Needed further down
            if cornersX(i,j) > max_x
                max_x = cornersX(i,j);
            end
            if cornersX(i,j) < min_x
                min_x = cornersX(i,j);
            end
        end
        if cornersY(i,j) >= 0
            cornersY(i,j) = cornersY(i,j) + 1;
            if cornersY(i,j) > max_y
                max_y = cornersY(i,j);
            end
            if cornersY(i,j) < min_y
                min_y = cornersY(i,j);
            end
        end
    end
end

%PREPARATIONS FOR PROPER PLOT ZOOM-IN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
min_x = min_x - (max_x - min_x)*(1 - numOfFoundCorners/numCorners + 0.2);
max_x = max_x + (max_x - min_x)*(1 - numOfFoundCorners/numCorners + 0.2);
min_y = min_y - (max_y - min_y)*(1 - numOfFoundCorners/numCorners + 0.2);
max_y = max_y + (max_y - min_y)*(1 - numOfFoundCorners/numCorners + 0.2);

min_x = max(min_x,0);
max_x = min(max_x,size2);
min_y = max(min_y,0);
max_y = min(max_y,size1);

% %ORIENTATION AMBIGUITY?
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Decide whether the position and orientation of the checkerboard
% %can be unambiguously determined
% if min((n_sq_y+1), (n_sq_x+1)) ~= min(cornerInfo(1),cornerInfo(2))
%     if flagStart == true
%         %DISPLAY AN ERROR AND RETURN
%     end
%     %There is some ambiguity, get rid of it
%     figure(2);
%     imagesc(I);
%     colormap(gray);
%     set(2,'color',[1 1 1]);
%     title({['Image ' num2str(kk)]});
%     h = get(gca, 'title');
%     set(h, 'FontWeight', 'bold')
%     axis([min_x max_x min_y max_y]);
%
%
%     figure(2); hold on;
%     i = startingCorner(1)
%     j = startingCorner(2)
%     %Plot the starting corner and its neighbor
%     plot( cornersX(i,j),cornersY(i,j),'+','color','red','linewidth',2);
%     plot( cornersX(i,j+1),cornersY(i,j+1),'+','color','red','linewidth',2);
%     text( cornersX(i,j)+3,cornersY(i,j)+3,num2str(0) )
%     text( cornersX(i,j+1)+3,cornersY(i,j+1)+3,num2str(1) )
%     set(findobj('type', 'text'), 'color', 'red');
%     hold off;
%
%     %Get user input on behalf of the board orientation
%     numCornersDirZeroOne = input(['The automatic corner finder was not able to decide how the pattern is oriented. Please indicate the number of corners in direction of corners [0 -> 1]: ']);
%     %We look in row direction
%     deltaCols = cornerInfo(2) - numCornersDirZeroOne;
% end

%DRAW AND NUMBER FOUND CORNERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Draw the found corners onto the image and number them
%Define the first encountered found corner as number "zero"
PlotXxiX = [];
PlotXxiY = [];
PlotCornersX = [];
PlotCornersY = [];
PlotCornerNumber = [];

for i = 1:1:rows
    for j = 1:1:(cols-deltaCols)
        if cornersX(i,j) ~= -1
            %Save for plotting later
            PlotCornersX = [PlotCornersX,cornersX(i,j)];
            PlotCornersY = [PlotCornersY,cornersY(i,j)];
            PlotCornerNumber = [PlotCornerNumber, cornerNumber];
            if use_corner_find
                %Apply the corner finder
%                 [xxi] = cornerfinder([cornersX(i,j);cornersY(i,j)],I,winty,wintx);
                %% ================
                %  added(changed code steffen urban
                wintx = ceil(calib_data.ocam_model.height/100);
                winty = ceil(calib_data.ocam_model.height/100);
                [xxi] = cornerfinder([cornersX(i,j);cornersY(i,j)], single(I), winty, wintx);

                cornersX(i,j) = xxi(1);
                cornersY(i,j) = xxi(2);

                %% ================
                %Save for plotting later
                PlotXxiX = [PlotXxiX,xxi(1)];
                PlotXxiY = [PlotXxiY,xxi(2)];
            end
            if cornerNumber == 0
                %Change starting corner to this one
                %Needed further down
                startingCorner = [i,j];
                startingCornerSet = true;
            end
        end
        if startingCornerSet == true
            cornerNumber = cornerNumber + 1;
        end
    end
end

if 0
    figure(2);
    imagesc(I);
    colormap(gray);
    set(2,'color',[1 1 1]);
    title({['Image ' num2str(kk)]; [num2str(numOfFoundCorners) ' / ' num2str(numCorners) ' corner have been found.']; ['Press ENTER to continue.']});
    h = get(gca, 'title');
    set(h, 'FontWeight', 'bold')
    axis([min_x max_x min_y max_y])

    %Plot the original corners
    figure(2); hold on;
    plot( PlotCornersX,PlotCornersY,'+','color','red','linewidth',2);
    if use_corner_find
        %Plot the "corner finder enhanced" corners
        plot( PlotXxiX,PlotXxiY,'+','color','blue','linewidth',2);
    end

    text( PlotCornersX'+3,PlotCornersY'+3,num2str(PlotCornerNumber') )
    set(findobj('type', 'text'), 'color', 'red');
    drawnow;
    hold off
end
%If no negative corners were designated, then "min_i" and "min_j"
%are still empty
if isempty(min_i)
    min_i = startingCorner(1);
    min_j = startingCorner(2);
end

%SAVE CORNER INFORMATION IN 2 VECTORS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Save all corners in two arrays for further processing
%by other functions
x = [];
y = [];

%Start with the smallest found corner and then append the larger ones
iteration = 0;
while true
    %Iterators
    i = min_i + floor(iteration/(cols-deltaCols));
    j = mod( (min_j - 1 + iteration),cols-deltaCols ) + 1;
    iteration = iteration + 1;

    x = [x,cornersX(i,j)];
    y = [y,cornersY(i,j)];

    if(iteration >= numCorners | i > rows)
        break
    end
end

%     %DOES NOT WORK RELIABLY!
%     %3. In case of a n x m board, where n is even and m is odd (or vice
%     %versa) there still exists a 180 degree ambiguity. This could be get
%     %rid off here.
%     %We define the starting corner as the corner where "a black checker is
%     %outernmost". Only proceed if one dimension is odd and the other even!
%     if (mod(n_cor_min,2) ~= 0 & mod(n_cor_max,2) == 0) | (mod(n_cor_min,2) == 0 & mod(n_cor_max,2) ~= 0)
%         %Determine the intensity at the given locations in the image
%         [dummy,lengthl] = size(x);
%         intens1 = I( floor((y(1,1)+y(1,n_cor_max+2))/2), floor((x(1,1)+x(1,n_cor_max+2))/2) )
%         intens2 = I( floor((y(1,lengthl)+y(1,lengthl-n_cor_max-1))/2), floor((x(1,lengthl)+x(1,lengthl-n_cor_max-2))/2) )
%         I(10,10)
%         if intens1 > intens2
%             %We need to reverse the numbering
%             xTemp = x;
%             yTemp = y;
%             for i = 1:1:lengthl
%                 x(i) = xTemp(lengthl+1 - i);
%                 y(i) = yTemp(lengthl+1 - i);
%             end
%         end
%     end

%REPOSITION CORNERS (IF NEEDED)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Allow the user (if needed) to reposition any of the corners
[dummy,sizex] = size(x);
iNumber = 1:1:sizex;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check, whether the numbering increases along the longer or the shorter pattern dimention:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_cor_min = min(calib_data.n_sq_x+1, calib_data.n_sq_y+1);
n_cor_max = max(calib_data.n_sq_x+1, calib_data.n_sq_y+1);

hmin = zeros(1,n_cor_max*n_cor_min-1);
hmin(1:n_cor_min:end) = 1;
hmax = zeros(1,n_cor_max*n_cor_min-1);
hmax(1:n_cor_max:end) = 1;

dxy = sqrt(diff(x).^2+diff(y).^2);

pmin = conv(hmin,dxy);
pmin = max(pmin(:));
pmax=conv(hmax,dxy);
pmax = max(pmax(:));
if pmin(1)>pmax(1)
    inc_dir = n_cor_min;
    oth_dir = n_cor_max;
else
    inc_dir = n_cor_max;
    oth_dir = n_cor_min;
end

% Rearrange numbering from starting point
xTemp = x;
yTemp = y;
area = [];
for i = 1:length(x)-1

    xborder = [xTemp(1,1:inc_dir-1),xTemp(1,inc_dir:inc_dir:end-inc_dir),xTemp(1,end:-1:end-inc_dir+2),xTemp(1,end-inc_dir+1:-inc_dir:1)];
    yborder = [yTemp(1,1:inc_dir-1),yTemp(1,inc_dir:inc_dir:end-inc_dir),yTemp(1,end:-1:end-inc_dir+2),yTemp(1,end-inc_dir+1:-inc_dir:1)];
    area(i) = abs(trapz(xborder,yborder));
    xTemp = [xTemp(1,2:end),xTemp(1,1)]; %Shift points one step forward
    yTemp = [yTemp(1,2:end),yTemp(1,1)]; %Shift points one step forward
end
shift = find( area == max(area) ) - 1;
if shift > 0
    x = [x(1,1+shift:end) , x(1,1:shift)];
    y = [y(1,1+shift:end) , y(1,1:shift)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
    figure(2);
    imagesc(I);
    colormap(gray);
    hold on;
    %Plot the corners
    plot( x,y,'+','color','red','linewidth',2);
    %Plot the corresponding number
    text( x'+3,y'+3,num2str(iNumber') )
    set(findobj('type', 'text'), 'color', 'red');
    axis([min_x max_x min_y max_y]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ASSIGN STARTING CORNER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This algorithm was first designed to be used with square patterns only,
%where the starting corner is meaningless, since every orientation +n*90
%degrees is structurally equivalent.
%With the extension to non-square checker boards, this is no longer the
%case and, depending on the specific board, 2 or 4 orientations can be
%distinguished.
%1. Check, whether we are dealing with a non-square board:
if calib_data.n_sq_x ~= calib_data.n_sq_y
    %2. Check, whether the numbering increases along the longer pattern
    %dimention:
%     n_cor_min = min(n_sq_x+1, n_sq_y+1);
%     n_cor_max = max(n_sq_x+1, n_sq_y+1);
%     dist1 = (x(1,n_cor_min)-x(1,n_cor_min+1))^2 + (y(1,n_cor_min)-y(1,n_cor_min+1))^2;
%     dist2 = (x(1,n_cor_max)-x(1,n_cor_max+1))^2 + (y(1,n_cor_max)-y(1,n_cor_max+1))^2;

    n_cor_x = calib_data.n_sq_x+1;
    n_cor_y = calib_data.n_sq_y+1;
    dist1 = (x(1,n_cor_x)-x(1,n_cor_x+1))^2 + (y(1,n_cor_x)-y(1,n_cor_x+1))^2;
    dist2 = (x(1,n_cor_y)-x(1,n_cor_y+1))^2 + (y(1,n_cor_y)-y(1,n_cor_y+1))^2;

    if dist1 > dist2
       %We have it wrongly numbered, renumber
       xTemp = x;
       yTemp = y;
       [dummy,lengthl] = size(x);
       iterMult = n_cor_x;
       iterOffset = 0;
       for i = 1:1:lengthl
           j = mod(i-1,n_cor_y)+1;
           x(i) = xTemp(j*iterMult-iterOffset);
           y(i) = yTemp(j*iterMult-iterOffset);
           if j*iterMult > n_cor_x*(n_cor_y-1)
                iterOffset = iterOffset + 1;
           end
       end
    end
end

Yp_abs = x';
Xp_abs = y';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1
    figure(2);
    imagesc(I);
    colormap(gray);
    title(['Image ' num2str(kk)]);
    h = get(gca, 'title');
    set(h, 'FontWeight', 'bold')
    hold on;
    %Plot the corners
    plot( x,y,'+','color','red','linewidth',2);
    %Plot the corresponding number
    text( x'+3,y'+3,num2str(iNumber') )
    set(findobj('type', 'text'), 'color', 'red');
    axis([min_x max_x min_y max_y]);
    draw_axes(Xp_abs, Yp_abs,calib_data.n_sq_y);
    drawnow;
    hold off;
%     close(2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DELETE FILES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%We delete the interface files between Matlab and c++, in order to prevent
%reloading old data in case of some errors.
delete('autoCornerFinder/cToMatlab/cornerInfo.txt');
delete('autoCornerFinder/cToMatlab/cornersX.txt');
delete('autoCornerFinder/cToMatlab/cornersY.txt');
delete('autoCornerFinder/cToMatlab/error.txt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Display finished message
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'Done\n');

author: urbste

repo: ImprovedOcamCalib


text:
function el = export_graphic_layer_module_field(args)
%"export_graphic_layer_module_field"
%   Returns a properly populated and formatted instance of tags in
%   graphic layer sequence.
%
%   For speed, tag must be a decimal representation of the 8 digit
%   hexidecimal DICOM tag desired, ie instead of '00100010', pass
%   hex2dec('00100010');
%
%   Arguments are passed in a structure, arg:
%       arg.tag         = decimal tag of field to fill
%       arg.data        = CERR structure(s) to fill from
%       arg.template    = an empty template of the module created by the
%                         function build_module_template.m
%
%   This function requires arg.data is a single doseS.
%
%APA 07/26/2019
%
%Usage:
%   dcmobj = export_rt_dose_module_field(args)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
template    = args.template;

switch tag
    case 7340128   % 0070,0060
        templateEl  = template.getValue(tag);
        fHandle = @export_graphic_layer_sequence;

        %New null sequence
        tmp = org.dcm4che3.data.Attributes;
        el = tmp.newSequence(tag, 0);

        i = 1; % only one layer
        dcmobj = export_sequence(fHandle, templateEl, gspsS);
        %dcmobj = export_sequence(fHandle, tag, {structS(i), i});
        el.add(i-1, dcmobj);

        %get attribute to return
        el = el.getParent();

    otherwise
        warning(['No methods exist to populate DICOM graphic_layer module field ' dec2hex(tag,8) '.']);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [hpb] = MW2hpb(MW)
% Convert power from megawatts to boiler horsepower.
% Chad A. Greene 2012
hpb = MW*101.941995005;


text:
function []=Rotor(be,th,ze)% be=[0 10 2];
% th=[0 12 8];
% ze=[0 0 0];

load NACA632_615.txt;
xy=NACA632_615;

x=xy(:,1);
y=xy(:,2);

rR=10;
R=100;
Spar_Loc=[2*3.5;0;2*0.3];

x=x/5-Spar_Loc(1);
y=y/5-Spar_Loc(3);
x=[x;x(1)];
y=[y;y(1)];

deg=pi/180;
b=be*deg;
t=th*deg;
z=ze*deg;

beta=inline('b1+b2*cos(psi)+b3*sin(psi)','b1','b2','b3','psi');
zeta=inline('z1+z2*cos(psi)+z3*sin(psi)','z1','z2','z3','psi');
theta=inline('t1+t2*cos(psi)+t3*sin(psi)','t1','t2','t3','psi');

psi=linspace(0,2*pi,50);
n=10;
r=linspace(0,100,n);

[rr ssi]=meshgrid(r,psi);

bbb=beta(b(1),b(2),b(3),ssi);
zzz=zeta(z(1),z(2),z(3),ssi);

xx=rr.*cos(bbb).*cos(ssi-zzz);
yy=rr.*cos(bbb).*sin(ssi-zzz);
zz=rr.*sin(bbb);

zmin=-15;
zmax=max(max(zz));

Rotor_Cone=surf(xx,yy,zz,repmat(7,size(zz)),'FaceAlpha',0);   % handle !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set(Rotor_Cone,'Visible','off')
hold on
HP=fill( 100*cos(psi),100*sin(psi),'w' );                       % handle !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set(HP,'FaceAlpha',0)
set(HP,'Visible','off')

TPP=fill3(xx(:,n),yy(:,n),zz(:,n),'g','FaceAlpha',0.5);         % handle !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set(TPP,'Visible','off')

NFP=fill3( 100*cos(psi),100*sin(psi),zeros(size(psi)),'m' ,'FaceAlpha',0.5);
rotate(NFP,[1 0 0],th(2),[0 0 0])
rotate(NFP,[0 1 0],th(3),[0 0 0])
set(NFP,'Visible','off')

HP_Cross(1)=plot3([-100 100],[0 0],[0 0]);
HP_Cross(2)=plot3([0 0],[-100 100],[0 0]);
set(HP_Cross,'Color','y','Linewidth',2)
set(HP_Cross,'Visible','off')

% ############################################## Plot limits ####################################################################
xlim([-100 100])
ylim([-100 100])
zlim([zmin zmax])

[Rotor_Shaft(1) Rotor_Shaft(2) Rotor_Shaft(3)]=Cylinder( [0 0 0],2.5,zmin,'z',15,'closed' );
set(Rotor_Shaft,'FaceColor','w')
set(Rotor_Shaft,'Visible','off')

[XAXIS(1) XAXIS(2) XAXIS(3)]=arrow3d( [0 0 0.5]',0.5,90,8,'x');
[YAXIS(1) YAXIS(2) YAXIS(3)]=arrow3d( [0 0 0.5]',0.5,90,8,'y');
[ZAXIS(1) ZAXIS(2) ZAXIS(3)]=arrow3d( [0 0 0.5]',0.5,40,4,'z');

XAXIS(4)=text(110,0,0,'X','Color','r','FontSize',12);
YAXIS(4)=text(0,110,0,'Y','Color','r','FontSize',12);
ZAXIS(4)=text(0,0,zmax+10,'Z','Color','r','FontSize',12);

set(XAXIS(1),'FaceColor','k','EdgeAlpha',0)
set(XAXIS(2),'FaceColor','r')
set(XAXIS(3),'FaceColor','k')

set(YAXIS(1),'FaceColor','k','EdgeAlpha',0)
set(YAXIS(2),'FaceColor','r')
set(YAXIS(3),'FaceColor','k')

set(ZAXIS(1),'FaceColor','k','EdgeAlpha',0)
set(ZAXIS(2),'FaceColor','r')
set(ZAXIS(3),'FaceColor','k')

set(XAXIS,'Visible','off')
set(YAXIS,'Visible','off')
set(ZAXIS,'Visible','off')

PSI=0:10:360;

Targetx=100*cos(PSI*pi/180);
Targety=100*sin(PSI*pi/180);
looktarget=plot3(Targetx',zeros(size(Targetx')),Targety',Targetx',Targety',zeros(size(Targetx')),zeros(size(Targetx')),Targetx',Targety');
set(looktarget,'Visible','off')

BETA(1,:)=beta(b(1),b(2),b(3),PSI*deg);
THETA(1,:)=theta(t(1),t(2),t(3),PSI*deg);
ZETA(1,:)=zeta(z(1),z(2),z(3),PSI*deg);

tr=[t(1) t(2)-b(3) b(2)+t(3)];
br=[b(1) 0 0];
zr=z;

BETA(2,:)=beta(br(1),br(2),br(3),PSI*deg);
THETA(2,:)=theta(tr(1),tr(2),tr(3),PSI*deg);
ZETA(2,:)=zeta(zr(1),zr(2),zr(3),PSI*deg);

br=[b(1) b(2)+t(3) b(3)-t(2)];
tr=[t(1) 0 0];
zr=z;

BETA(3,:)=beta(br(1),br(2),br(3),PSI*deg);
THETA(3,:)=theta(tr(1),tr(2),tr(3),PSI*deg);
ZETA(3,:)=zeta(zr(1),zr(2),zr(3),PSI*deg);

for i=1:37;
    Blade_Disk(i,:)=Make_Blade(x,y,rR,R,PSI(i),THETA(1,i)/deg,BETA(1,i)/deg,ZETA(1,i)/deg);
    set(Blade_Disk(i,:),'Visible','off')

    a_HP(i,1)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0,['HP: \psi = ' num2str(PSI(i)) ],'Color','k','FontSize',12);
    a_HP(i,2)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0-10,['HP: \beta = ' num2str(BETA(1,i)/deg) ],'Color','r','FontSize',12);
    a_HP(i,3)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0+10,['HP: \theta = ' num2str(THETA(1,i)/deg) ],'Color','m','FontSize',12);
    a_HP(i,4)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0+20,['HP: \zeta = ' num2str(ZETA(1,i)/deg) ],'Color','b','FontSize',12);
    set(a_HP(i,:),'Visible','off')

    a_TPP(i,1)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0,['TPP: \psi = ' num2str(PSI(i))],'Color','k','FontSize',12);
    a_TPP(i,2)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0-10,['TPP: \beta = ' num2str(BETA(2,i)/deg) ],'Color','r','FontSize',12);
    a_TPP(i,3)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0+10,['TPP: \theta = ' num2str(THETA(2,i)/deg) ],'Color','m','FontSize',12);
    a_TPP(i,4)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0+20,['TPP: \zeta = ' num2str(ZETA(2,i)/deg) ],'Color','b','FontSize',12);
    set(a_TPP(i,:),'Visible','off')

    a_NFP(i,1)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0,['NFPP: \psi = ' num2str(PSI(i)) ],'Color','k','FontSize',12);
    a_NFP(i,2)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0-10,['NFPP: \beta = ' num2str(BETA(3,i)/deg) ],'Color','r','FontSize',12);
    a_NFP(i,3)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0+10,['NFPP: \theta = ' num2str(THETA(3,i)/deg) ],'Color','m','FontSize',12);
    a_NFP(i,4)=text(90*cos(PSI(i)*deg),90*sin(PSI(i)*deg),0+20,['NFPP: \zeta = ' num2str(ZETA(3,i)/deg) ],'Color','b','FontSize',12);
    set(a_NFP(i,:),'Visible','off')
end

axis equal
axis off

daspect([1 1 1])
daspect('manual')
camlookat(looktarget)
camzoom(1.5)

save para Rotor_Shaft HP TPP Rotor_Cone  HP_Cross XAXIS YAXIS ZAXIS Blade_Disk NFP a_HP a_TPP a_NFP

save looktargetfile looktarget


text:
function out = UpDown(in)

UD = {'outOfPlane','intoPlane'};

if nargin == 0
  out = UD;
elseif isnumeric(in)
  out = UD{in};
elseif ischar(in)
  ind = strcmpi(in,UD);
  out = find(ind);
else
  error('Input must be a string or a number!');
end


text:
function x = p42_start ( n )

%*****************************************************************************80
%
%% P42_START returns a starting point for optimization for problem 42.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = [ 0.0, 1.0, 2.0 ]';

  return
end


text:
function [in2] = cm22in2(cm2)
% Convert area from square centimeters to square inches.
% Chad A. Greene 2012
in2 = cm2*0.1550003100006 ;


text:
function [uout, tout] = spinsphere(varargin)
%SPINSPHERE  Solve stiff PDEs on the sphere, double Fourier sphere method and
%implicit-explicit schemes.
%
%   UOUT = SPINSPHERE(PDECHAR) solves the PDE specified by the string PDECHAR,
%   and plays a movie of the solution. Possible strings include 'AC' and 'GL'
%   for the Allen-Cahn and Ginzburg-Landau equations. Other PDEs are available,
%   see Remark 1 and Examples 1-4. The output UOUT is a SPHEREFUN corresponding
%   to the solution at the final time (a CHEBMATRIX for systems of equations,
%   each row representing one variable).
%
%   UOUT = SPINSPHERE(S, N, DT) solves the PDE specified by the SPINOPSPHERE S
%   with N grid points in each direction (longitude/latitude) and time-step DT,
%   and plays a movie of the solution. See HELP/SPINOPSPHERE and Example 5.
%
%   UOUT = SPINSPHERE(S, N, DT, PREF) allows one to use the preferences
%   specified by the SPINPREFSPHERE object PREF. See HELP/SPINPREFSPHERE and
%   Example 6.
%
%   [UOUT, TOUT] = SPINSPHERE(...) also returns the times chunks TOUT at which
%   UOUT was computed.
%
%   Users of SPINSPHERE will quickly find they want to vary aspects of the
%   plotting. The fully general syntax for this involves using preferences
%   specified by a SPINPREFSPHERE object PREF. See HELP/SPINPREFSPHERE and
%   Example 6. However for many purposes it is most convenient to use the syntax
%
%   UOUT = SPINSPHERE(..., 'PREF1', VALUE1, 'PREF2', VALUE2, ...)
%
%   For example:
%
%   UOUT = SPINSPHERE(..., 'Clim', [a b]) changes colorbar limits to [a b]
%   UOUT = SPINSPHERE(..., 'colormap', 'jet') changes the colormap to 'jet'
%   UOUT = SPINSPHERE(..., 'dataplot', 'abs') plots absolute value
%   UOUT = SPINSPHERE(..., 'grid', 'on') for lagitude/longitude circles
%   UOUT = SPINSPHERE(..., 'iterplot', 4) plots only every 4th time step
%   UOUT = SPINSPHERE(..., 'Nplot', 256) plays a movie at 256x256 resolution
%   UOUT = SPINSPHERE(..., 'plot', 'off') for no movie
%   UOUT = SPINSPHERE(..., 'view', [a b]) changes the view angle to [a b]
%
% Remark 1: List of PDEs (case-insensitive)
%
%    - 'AC' for the Allen-Cahn equation,
%    - 'GL' for the Ginzburg-Landau equation,
%    - 'GM' for the Gierer-Meinhardt equations,
%    - 'NLS' for the focusing nonlinear Schroedinger equation.
%
% Example 1: Allen-Cahn equation (metastable solutions)
%
%        u = spinsphere('AC');
%
%    solves the Allen-Cahn equation
%
%        u_t = 1e-2*laplacian(u) + u - u^3
%
%    on the sphere from t=0 to t=60, with initial condition
%
%        u0(x, y, z) = cos(cosh(5*x*z) - 10*y).
%
% Example 2: Ginzburg-Landau equation (spiral waves)
%
%        u = spinsphere('GL');
%
%    solves the Ginzburg-Landau equation
%
%        u_t = 1e-3*laplacian(u) + u - (1+1.5i)*u*|u|^2,
%
%    on the sphere from t=0 to t=100 with a RANDNFUNSPHERE initial condition.
%    The movie shows the real part of u.
%
% Example 3: Gierer-Meinhardt equations (pattern formation - spots)
%
%        u = spinsphere('GM);
%
%    solves the Gierer-Meinhardt equations,
%
%       u_t = 1e-2*laplacian(u) + u^2/v - u,
%       v_t = 1e-1*laplacian(v) + u^2 - v,
%
%    on the sphere from t=0 to t=80, with initial condition
%
%       u0(x,y,z) = 1 + .1*(cos(20*x) + cos(20*y) + cos(20*z)),
%       v0(x,y,z) = 1 - .1*(cos(20*x) + cos(20*y) + cos(20*z)).
%
% Example 4: Nonlinear Schroedinger equation (spherical harmonic & breather)
%
%        u = spinsphere('NLS');
%
%    solves the focusing nonlinear Schroedinger equation
%
%        u_t = 1i*laplacian(u) + 1i*u|u|^2,
%
%    on the sphere from t=0 to t=3, with initial condition
%
%     u0(lam, th) = .1*(2*B^2./(2 - sqrt(2)*sqrt(2-B^2)*cos(A*B*th)) - 1)*A
%                  + Y_8^6(lam, th), with A=1 and B=1.
%
%    The movie shows the absolute value of u.
%
% Example 5: PDE specified by a SPINOPSPHERE
%
%       tspan = [0 100];
%       S = spinopsphere(tspan);
%       S.lin = @(u) 1e-3*lap(u);
%       S.nonlin = @(u) u - (1 + 1.5i)*u.*(abs(u).^2);
%       S.init = randnfunsphere(.1);
%       S.init = S.init/norm(S.init, inf);
%       u = spinsphere(S, 128, 1e-1);
%
%   is equivalent to u = spinsphere('GL');
%
% Example 6: Using preferences
%
%       pref = spinprefsphere('Clim', [-1 1]);
%       S = spinopsphere('AC');
%       u = spinsphere(S, 128, 1e-1, pref);
%   or simply,
%       u = spinsphere(S, 128, 1e-1, 'Clim', [-1 1]);
%
%   solves the Allen-Cahn equation using N=128 grid points in each direction
%   and a time-step dt=1e-1, and sets the limits of the colorbar to [-1 1].
%
% See also SPINOPSPHERE, SPINPREFSPHERE, IMEX.

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

% We are going to parse the inputs and call SOLVEPDE in the following ways,
%
%       SPINOPERATOR.SOLVEPDE(S, N, dt)
%  or
%       SPINOPERATOR.SOLVEPDE(S, N, dt, pref)
%
% where S is a SPINOPSPHERE object, N is the number of grid points in each
% direction, DT is the time-step and PREF is a SPINPREFSPHERE object.

% CASE 1. U = SPINSPHERE('GL'):
if ( nargin == 1 )

    try spinopsphere(varargin{1});
    catch
        error('Unrecognized PDE. See HELP/SPINSPHERE for the list of PDEs.')
    end
    [S, N, dt, pref] = parseInputs(varargin{1});
    varargin{1} = S;
    varargin{2} = N;
    varargin{3} = dt;
    varargin{4} = pref;

% CASE 2. U = SPINSPHERE('GL', 'PREF1', VALUE1) or U = SPINSPHERE(S, N, DT):
elseif ( nargin == 3 )

    % CASE 2.1. U = SPINSPHERE('GL', 'PREF1', VALUE1):
    if ( isa(varargin{1}, 'char') == 1 && isa(varargin{2}, 'char') == 1 )
        [S, N, dt, pref] = parseInputs(varargin{1});
        pref.(varargin{2}) = varargin{3};
        varargin{1} = S;
        varargin{2} = N;
        varargin{3} = dt;
        varargin{4} = pref;

    % CASE 2.2. U = SPINSPHERE(S, N, DT):
    else
        % Nothing to do here.
    end

% CASE 3. U = SPINSPHERE(S, N, DT, PREF)
elseif ( nargin == 4 )
    % Nothing to do here.

% CASE 4.
elseif ( nargin >= 5 )

    % CASE 4.1. U = SPINSPHERE('GL', 'PREF1', VALUE1, 'PREF2', VALUE2, ...)
    if ( isa(varargin{1}, 'char') == 1 && isa(varargin{2}, 'char') == 1 )
        [S, N, dt, pref] = parseInputs(varargin{1});
        j = 2;
        while j < nargin
            pref.(varargin{j}) = varargin{j+1};
            varargin{j} = [];
            varargin{j+1} = [];
            j = j + 2;
        end
        varargin{1} = S;
        varargin{2} = N;
        varargin{3} = dt;
        varargin{4} = pref;
        varargin = varargin(~cellfun(@isempty, varargin));

    % CASE 4.2. U = SPINSPHERE(S, N, DT, 'PREF1', VALUE1, 'PREF2', VALUE2, ...)
    else
        pref = spinprefsphere();
        j = 4;
        while j < nargin
            pref.(varargin{j}) = varargin{j+1};
            varargin{j} = [];
            varargin{j+1} = [];
            j = j + 2;
        end
        varargin{4} = pref;
        varargin = varargin(~cellfun(@isempty, varargin));
    end

end

% SPINSPHERE is a wrapper for SOLVPDE:
[uout, tout] = spinoperator.solvepde(varargin{:});

end

function [S, N, dt, pref] = parseInputs(pdechar)
%PARSEINPUTS   Parse the inputs.

pref = spinprefsphere();
S = spinopsphere(pdechar);
if ( strcmpi(pdechar, 'AC') == 1 )
    dt = 1e-1;
    N = 128;
    pref.Clim = [-1 1];
    pref.iterplot = 2;
    pref.Nplot = 256;
elseif ( strcmpi(pdechar, 'GL') == 1 )
    dt = 1e-1;
    N = 128;
    pref.Clim = [-1 1];
    pref.iterplot = 2;
    pref.Nplot = 256;
elseif ( strcmpi(pdechar, 'GM') == 1 )
    dt = 2e-1;
    N = 64;
    pref.Clim = [0 3 0.5 2];
    pref.iterplot = 4;
    pref.Nplot = 128;
elseif ( strcmpi(pdechar, 'NLS') == 1 )
    dt = 1e-2;
    N = 128;
    pref.colormap = 'jet';
    pref.dataplot = 'abs';
    pref.Clim = [0 1];
    pref.iterplot = 1;
    pref.Nplot = 256;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function output = dissimilarInsert(dataStruct, newElement, index, nullData)
%"dissimilarInsert"
%   Add a dissimilar element to a datastructure, in position index.  If the
%   element has more fields than the structure, the new fields are added to
%   the structure and filled with nullData.  If the structure has more fields
%   than the element, the reverse occurs. Any data in that position is
%   deleted, as it would be in a normal indexing operation. If index is blank
%   or not entered it is assumed to be the end, resulting in append.
%
%   If nulldata is undefined [] is used.
%
%JRA 10/7/03
%JOD, 10 Nov 03, define nullData if needed.
%JRA, 10/1/04, Modifications to make compatible with ML6.1.
%JRA, 11/4/04, BugFix: Crashed if identical but misordered fieldnames.
%              Also added some more comments, and fixed initialization.
%
%Usage:
%   function output = dissimilarInsert(dataStruct, newElement, index, nullData)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%nullData defaults to [].
if nargin < 4
    nullData = [];
end

%Insert index defaults to appending.
if ~exist('index','var') || (exist('index','var') && isempty(index))
    index = length(dataStruct)+1;
end

%Get all fieldnames, find common fields.
structNames = fieldnames(dataStruct);
newNames = fieldnames(newElement);
commonNames = intersect(structNames,newNames);

%If same number of fields in both...
if length(commonNames) == length(structNames) && length(commonNames) == length(newNames)
    differences = find(~strcmpi(structNames, newNames));
    %If no differences...
    if isempty(differences)
        %then insert normally.
        dataStruct(index) = newElement;
    %If there were differences...
    else
        %Fieldnames were identical but misordered, insert with reordering.
        for i=1:length(newNames)
            dataStruct = setfield(dataStruct, {index}, newNames{i}, getfield(newElement, newNames{i}));
        end
        %warning('Element''s fieldnames were misordered, but the same.')
    end
    output = dataStruct;
    return;
end

%If more fields in new element...
if length(commonNames) < length(newNames)
    newFieldsForStruct = setdiff(newNames, commonNames);
    for i = 1:length(newFieldsForStruct)
        %dataStruct = setfield(dataStruct, {}, newFieldsForStruct{i}, []);
        for j=1:length(dataStruct)
            dataStruct = setfield(dataStruct, {j}, newFieldsForStruct{i}, nullData);
        end
    end
end

%If more fields in data structure...
if length(commonNames) < length(structNames)
    newFieldsForElement = setdiff(structNames, commonNames);
    for i=1:length(newFieldsForElement)
        newElement = setfield(newElement, newFieldsForElement{i}, nullData);
    end
end

%Now actually perform the insert.
finalNames = fieldnames(newElement);
siz = size(dataStruct);
if prod(siz) == 0
    indV = 1:2:2*length(finalNames);
    argV = {};
    [argV{indV}] = deal(finalNames{:});
    argV{end+1} = [];
    dataStruct = struct(argV{:});
end

for i=1:length(finalNames)
    dataStruct = setfield(dataStruct, {index}, finalNames{i}, getfield(newElement, finalNames{i}));
end
output = dataStruct;
return;


text:
function [J grad] = nnCostFunction(nn_params, ...
                                   input_layer_size, ...
                                   hidden_layer_size, ...
                                   num_labels, ...
                                   X, y, lambda)
%NNCOSTFUNCTION Implements the neural network cost function for a two layer
%neural network which performs classification
%   [J grad] = NNCOSTFUNCTON(nn_params, hidden_layer_size, num_labels, ...
%   X, y, lambda) computes the cost and gradient of the neural network. The
%   parameters for the neural network are "unrolled" into the vector
%   nn_params and need to be converted back into the weight matrices.
%
%   The returned parameter grad should be a "unrolled" vector of the
%   partial derivatives of the neural network.
%

% Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
% for our 2 layer neural network
Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...
                 hidden_layer_size, (input_layer_size + 1));

Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...
                 num_labels, (hidden_layer_size + 1));

% Setup some useful variables

m = size(X, 1);
X = [ones(size(X,1),1), X];

% You need to return the following variables correctly
J = 0;
Theta1_grad = zeros(size(Theta1));
Theta2_grad = zeros(size(Theta2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the code by working through the
%               following parts.
%
% Part 1: Feedforward the neural network and return the cost in the
%         variable J. After implementing Part 1, you can verify that your
%         cost function computation is correct by verifying the cost
%         computed in ex4.m
%
% Part 2: Implement the backpropagation algorithm to compute the gradients
%         Theta1_grad and Theta2_grad. You should return the partial derivatives of
%         the cost function with respect to Theta1 and Theta2 in Theta1_grad and
%         Theta2_grad, respectively. After implementing Part 2, you can check
%         that your implementation is correct by running checkNNGradients
%
%         Note: The vector y passed into the function is a vector of labels
%               containing values from 1..K. You need to map this vector into a
%               binary vector of 1's and 0's to be used with the neural network
%               cost function.
%
%         Hint: We recommend implementing backpropagation using a for-loop
%               over the training examples if you are implementing it for the
%               first time.
%
% Part 3: Implement regularization with the cost function and gradients.
%
%         Hint: You can implement this around the code for
%               backpropagation. That is, you can compute the gradients for
%               the regularization separately and then add them to Theta1_grad
%               and Theta2_grad from Part 2.
%

Z2 = X*Theta1';
A2 = sigmoid(Z2);

A2 = [ones(size(A2,1),1), A2];
Z3 = A2*Theta2';

A3 = sigmoid(Z3);

% Y = zeros(m, max(y(:)));
% Y(y,:)

Y = dummyvar(y);

J = sum(sum(-1.*Y.*log(A3)-(1-Y).*log(1-A3),2))/m;

reg = (sum(sum(Theta1(:,2:end).^2)) + sum(sum(Theta2(:,2:end).^2)))*lambda/(2*m);
J = J+reg;

delta3 = A3 - Y;

Theta2_grad = Theta2_grad + (delta3'*A2)/m;

delta2 = delta3*Theta2(:,2:end).*sigmoidGradient(Z2);

Theta1_grad = Theta1_grad+(delta2'*X)/m;

Theta1_grad(:,2:end) = Theta1_grad(:,2:end) + lambda*Theta1(:,2:end)/m;
Theta2_grad(:,2:end) = Theta2_grad(:,2:end) + lambda*Theta2(:,2:end)/m;
% -------------------------------------------------------------

% =========================================================================

% Unroll gradients
grad = [Theta1_grad(:) ; Theta2_grad(:)];

end

repo: machine-learning-coursera