text:
function mAP = rfcn_test(conf, imdb, roidb, varargin)
% Image-level testing using final D model and RPN proposals
% --------------------------------------------------------
% D&T implementation
% Modified from MATLAB  R-FCN (https://github.com/daijifeng001/R-FCN/)
% and Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2017, Christoph Feichtenhofer
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

%% inputs
    ip = inputParser;
    ip.addRequired('conf',                              @isstruct);
    ip.addRequired('imdb',                              @isstruct);
    ip.addRequired('roidb',                             @isstruct);
    ip.addParamValue('net_def_file',    '', 			@isstr);
    ip.addParamValue('net_file',        '', 			@isstr);
    ip.addParamValue('cache_name',      '', 			@isstr);
    ip.addParamValue('suffix',          '',             @isstr);
    ip.addParamValue('ignore_cache',    false,          @islogical);
    ip.addParamValue('proposal_path',    '',          @isstr);
    ip.addParamValue('visualize_every',   Inf,          @isscalar);
    ip.addParamValue('test_classes',    [] 			);

    ip.parse(conf, imdb, roidb, varargin{:});
    opts = ip.Results;
    num_proposals = 300;
    if isfield(conf, 'num_proposals')
      num_proposals = conf.num_proposals;
    end
    visualize = false;

%%  set cache dir
    cache_dir = fullfile(conf.root_path, 'output', 'rfcn_cachedir', opts.cache_name, imdb.name);
    mkdir_if_missing(cache_dir);

%%  init log
    timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
    mkdir_if_missing(fullfile(cache_dir, 'log'));
    log_file = fullfile(cache_dir, 'log', ['test_', timestamp, '.txt']);
    diary(log_file);

    num_images = length(imdb.image_ids);
    num_classes = imdb.num_classes;

    try
      aboxes = cell(num_classes, 1);
      if opts.ignore_cache
          throw('');
      end
      for i = 1:num_classes
        load(fullfile(cache_dir, [imdb.classes{i} '_boxes_' imdb.name opts.suffix]));
        aboxes{i} = boxes;
      end
    catch
%%      testing
        % init caffe net
        caffe_log_file_base = fullfile(cache_dir, 'caffe_log');
        caffe.init_log(caffe_log_file_base);
        caffe_net = caffe.Net(opts.net_def_file, 'test');
        caffe_net.copy_from(opts.net_file);

        % set random seed
        prev_rng = seed_rand(conf.rng_seed);
        caffe.set_random_seed(conf.rng_seed);

        % set gpu/cpu
        if conf.use_gpu
            caffe.set_mode_gpu();
        else
            caffe.set_mode_cpu();
        end

        % determine the maximum number of rois in testing
        max_rois_num_in_gpu = 10000;

        disp('opts:');
        disp(opts);
        disp('conf:');
        disp(conf);

        %heuristic: keep an average of 160 detections per class per images prior to NMS
        max_per_set = 160 * num_images;
        % heuristic: keep at most 400 detection per class per image prior to NMS
        max_per_image = 400;
        % detection thresold for each class (this is adaptively set based on the max_per_set constraint)
        thresh = -inf * ones(num_classes, 1);
        % top_scores will hold one minheap of scores per class (used to enforce the max_per_set constraint)
        top_scores = cell(num_classes, 1);
        % all detections are collected into:
        %    all_boxes[cls][image] = N x 5 array of detections in
        %    (x1, y1, x2, y2, score)
        aboxes = cell(num_classes, 1);
        box_inds = cell(num_classes, 1);
        for i = 1:num_classes
            aboxes{i} = cell(length(imdb.image_ids), 1);
            box_inds{i} = cell(length(imdb.image_ids), 1);
        end
        count = 0;
        t_start = tic;
        for i = 1:1:num_images
            d = roidb.rois(i);
            if ~isempty(opts.test_classes),
              testThis = false;
              for ii = 1:length(opts.test_classes)
                if any(opts.test_classes(ii) == d.class), testThis = true; end
              end
              if ~testThis, continue; end
            end
            count = count + 1;
            fprintf('%s: test (%s) %d/%d ', procid(), imdb.name, count, num_images);
            th = tic;

            if ~isempty(opts.proposal_path)
              file = imdb.image_ids{i};
              file(file == '\') = '/';
              ind = find(file == '/', 1, 'first');
              videoFrame = file ;
              if strcmp(file(1:20), 'ILSVRC2015_VID_train')
                videoFrame = file(ind+1:end);
              end
              ab_fetch_file = [opts.proposal_path '/' videoFrame '.mat'];
              if ~exist(ab_fetch_file,'file')
                  fprintf('In roidb_from_proposal_manual.m --> error:file does not exist : \n %s \n', ab_fetch_file);
                  fprintf('generate this file and run script again, pausing for now...\n');
                  pause;
              end

              abst = load(ab_fetch_file);
              if isfield(abst, 'boxes')
                proposals = abst.boxes;
              else
                proposals = abst.aboxes;
              end

              if num_proposals < 1
                  proposals = proposals(proposals(:,5) > num_proposals,1:4);
              else
                  proposals = proposals(1:min(num_proposals,size(proposals,1)),1:4);
              end
              if isempty(proposals), continue; end

              d.boxes = cat(1, d.boxes, proposals);
              d.gt = cat(1, d.gt, zeros(size(proposals,1),1));
            end

            im_path = imdb.image_at(i);
            if isfield(conf, 'input_modality') && ~strcmp(conf.input_modality, 'rgb')

              img_file_u = strrep(im_path, ['Data' filesep 'VID'], ...
                ['Data' filesep 'VID' filesep 'tvl1_flow_600' filesep 'u'] ) ;
              img_file_v = strrep(im_path, ['Data' filesep 'VID'], ...
                ['Data' filesep 'VID' filesep 'tvl1_flow_600' filesep 'v'] ) ;
              try
                 nFrames = 1;
                 if isfield(conf, 'nFrames')
                   nFrames = conf.nFrames;
                 end
                [path_u, ~, ~] = fileparts(img_file_u);
                [path_v, frame, ext] = fileparts(img_file_v);
                frames = str2num(frame)-nFrames+1:str2num(frame)+nFrames-1;
                vid = fileparts(videoFrame);
                frames(frames <0) = 0; frames(frames > imdb.nFrames(vid)-2) =  imdb.nFrames(vid)-2;
                im_u = []; im_v = [];
                for frame=frames
                  img_file_u = fullfile(path_u, [sprintf('%06d',frame), ext]) ;
                  img_file_v = fullfile(path_v, [sprintf('%06d',frame), ext]) ;
                  im_u{end+1} = imread(img_file_u); im_v{end+1} = imread(img_file_v);
                end
                im_u = cat(3,im_u{:});im_v = cat(3,im_v{:});

              catch % the last frame is simply replicated in optical flow case
                continue;
              end
                switch conf.input_modality
                  case 'flow'
                    flow = single(cat(3,im_u,im_v)) - 128;
                    flow = bsxfun(@minus,flow, median(median(flow,1),2)) ;
                    mag_flow = sqrt(sum(flow.^2,3)) - 128;
                    im = (cat(3,flow,mag_flow));
                    im = imresize( im, imdb.sizes(i,:));
                  case 'grayflow'
                    imgray = imresize( rgb2gray(imread(im_path)), sz(1:2));
                    im = cat(3,im_u,im_v,imgray);
                  case 'rgbflow'

                    if isfield(conf, 'subflow_mean') && ~conf.subflow_mean
                      flow = single(cat(3,im_u,im_v));
                    else
                      flow = single(cat(3,im_u,im_v)) - 128;
                    end

                    imrgb = single(imread(im_path));

                    if size(imrgb,3) < 3
                      imrgb = imrgb(:,:,[1 1 1]);
                    end
                    sz = size(imrgb);
                    im = cat(3,imrgb,imresize( flow, sz(1:2)));

                end
              else
                  im = imread(im_path);
              end

            [boxes, scores] = rfcn_im_detect(conf, caffe_net, im, d.boxes(d.gt==0, :), max_rois_num_in_gpu);

            for j = 1:num_classes
                inds = find(scores(:, j) > thresh(j));
                if ~isempty(inds)
                    [~, ord] = sort(scores(inds, j), 'descend');
                    ord = ord(1:min(length(ord), max_per_image));
                    inds = inds(ord);
                    cls_boxes = boxes(inds, (1+(j-1)*4):((j)*4));
                    cls_scores = scores(inds, j);
                    aboxes{j}{i} = [aboxes{j}{i}; cat(2, single(cls_boxes), single(cls_scores))];
                    box_inds{j}{i} = [box_inds{j}{i}; inds];
                else
                    aboxes{j}{i} = [aboxes{j}{i}; zeros(0, 5, 'single')];
                    box_inds{j}{i} = box_inds{j}{i};
                end
            end

            fprintf(' time %.3fs\n', toc(th));

            if mod(count,opts.visualize_every) == 0
             boxes_cell = cell(length(imdb.classes), 1);

              vis_thres = 0.3;
              for k = 1:length(boxes_cell)
                  boxes_cell{k} = [boxes(:, (1+(k-1)*4):(k*4)), scores(:, k)];

                  boxes_cell{k} = boxes_cell{k}(nms(boxes_cell{k}, 0.3), :);

                  I = boxes_cell{k}(:, 5) >= vis_thres;
                  boxes_cell{k} = boxes_cell{k}(I, :);
              end
              figure(1);
              showboxes(uint8(im(:,:,1:3)), boxes_cell, imdb.classes, 'default'); drawnow;
              if size(im,3) > 3
                  figure(2);
                  showboxes(uint8(im(:,:,4:6)+128), boxes_cell, imdb.classes, 'default'); drawnow;
              end

            end

            if mod(count, 1000) == 0
                for j = 1:num_classes
                [aboxes{j}, box_inds{j}, thresh(j)] = ...
                    keep_top_k(aboxes{j}, box_inds{j}, i, max_per_set, thresh(j));
                end
                disp(thresh);
            end
        end

        for j = 1:num_classes
            [aboxes{j}, box_inds{j}, thresh(j)] = ...
                keep_top_k(aboxes{j}, box_inds{j}, i, max_per_set, thresh(j));
        end
        disp(thresh);

        for i = 1:num_classes
            top_scores{i} = sort(top_scores{i}, 'descend');
            if (length(top_scores{i}) > max_per_set)
                thresh(i) = top_scores{i}(max_per_set);
            end

            % go back through and prune out detections below the found threshold
            for j = 1:length(imdb.image_ids)
                if ~isempty(aboxes{i}{j})
                    I = find(aboxes{i}{j}(:,end) < thresh(i));
                    aboxes{i}{j}(I,:) = [];
                    box_inds{i}{j}(I,:) = [];
                end
            end

            save_file = fullfile(cache_dir, [imdb.classes{i} '_boxes_' imdb.name opts.suffix]);
            boxes = aboxes{i};
            inds = box_inds{i};
            save(save_file, 'boxes', 'inds');
            clear boxes inds;
        end
        fprintf('test all images in %f seconds.\n', toc(t_start));

        caffe.reset_all();
        rng(prev_rng);
    end

    % ------------------------------------------------------------------------
    % Peform AP evaluation
    % ------------------------------------------------------------------------

    if isequal(imdb.eval_func, @imdb_eval_voc)
        new_parpool();
        parfor model_ind = 1:num_classes
          cls = imdb.classes{model_ind};
          res(model_ind) = imdb.eval_func(cls, aboxes{model_ind}, imdb, opts.cache_name, opts.suffix);
        end
    else
      if 0
        % sanity test:
        boxes_gt = cell(num_classes,1);
        classes = {roidb.rois.class};
        boxes = {roidb.rois.boxes};
        boxes = cellfun(@(x) cat(2,x,ones(size(x,1),1)),boxes,'UniformOutput' , false);
        for c=1:num_classes
          boxes_gt{c} = cell(num_images,1);
        end
        for j=1:num_images
          if ~isempty(classes{j})
            for c=unique(classes{j})'
              boxes_gt{c}{j} = boxes{j}(classes{j}==c,:);
            end
          end
        end
        res_gt = imdb.eval_func(boxes_gt, imdb, conf, opts.suffix);
      end
      % ilsvrc
      res = imdb.eval_func(aboxes, imdb, conf, opts.suffix);
    end

    if ~isempty(res)
        fprintf('\n~~~~~~~~~~~~~~~~~~~~\n');
        fprintf('Results:\n');
        aps = [res(:).ap]' * 100;
        disp(aps);
        disp(mean(aps));
        fprintf('~~~~~~~~~~~~~~~~~~~~\n');
        mAP = mean(aps);
    else
        mAP = nan;
    end

    diary off;
end

% ------------------------------------------------------------------------
function [boxes, box_inds, thresh] = keep_top_k(boxes, box_inds, end_at, top_k, thresh)
% ------------------------------------------------------------------------
    % Keep top K
    X = cat(1, boxes{1:end_at});
    if isempty(X)
        return;
    end
    scores = sort(X(:,end), 'descend');
    thresh = scores(min(length(scores), top_k));
    for image_index = 1:end_at
        if ~isempty(boxes{image_index})
            bbox = boxes{image_index};
            keep = find(bbox(:,end) >= thresh);
            boxes{image_index} = bbox(keep,:);
            box_inds{image_index} = box_inds{image_index}(keep);
        end
    end
end

author: feichtenhofer

repo: Detect-Track


text:
function bti_test04 ( )

%*****************************************************************************80
%
%% BTI_TEST04 tests BURGERS_TIME_INVISCID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTI_TEST04\n' );
  fprintf ( 1, '  Test BURGERS_TIME_INVISCID\n' );

  method = 4;
  nx = 81;
  nt = 200;
  t_max = 2.0;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Method: 4, Lax Wendroff.\n' );
  fprintf ( 1, '  Initial condition: expansion\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_inviscid ( method, @ic_expansion, nx, nt, t_max, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 4 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers, inviscid, initial expansion, Lax Wendroff' )

  filename = 'bti_test04.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dim_num = p02_dim_num ( )

%*****************************************************************************80
%
%% P02_DIM_NUM returns the spatial dimension for problem p02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DIM_NUM, the spatial dimension of the
%    dependent variables.
%
  dim_num = 2;

  return
end


text:
classdef ParadigmRSSD2 < ParadigmDataflowSimplified
    % Advanced paradigm for time-varying oscillatory processes.
    %
    % This paradigm allows to learn the joint space/time/frequency structure in EEG, under the
    % assumption that an overcomplete ICA decomposition can produce a reasonably complete coverage
    % of potentially relevant sources. For each of the resulting (more or less independent) source
    % signals, a full time/freq decomposition is computed in the time epoch of interest, and a
    % regression model is learned which maps from the (very high-dimensional) space of these
    % features onto the desired output variable. The regression is by default logistic, and is
    % "rank-sparse" (using the trace norm regularization), meaning that the time/freq weights
    % learned for each signal component will tend to have low rank (reducing the effective
    % # of parameters) and also a sparse set of component signals will have non-zero weights.
    %
    % In addition, a prior can be imposed on the weights, which is parameterized in space (brain area), time, and
    % frequency. This makes the overall model a fairly general-purpose approach to time-varying / non-stationary
    % oscillations.
    %
    % References:
    %  [1] Ryota Tomioka and Klaus-Robert Mueller, "A regularized discriminative framework for EEG analysis with application to brain-computer interface",
    %      Neuroimage, 49 (1) pp. 415-432, 2010.
    %  [2] Ryota Tomioka & Masashi Sugiyama, "Dual Augmented Lagrangian Method for Efficient Sparse Reconstruction",
    %      IEEE Signal Procesing Letters, 16 (12) pp. 1067-1070, 2009.
    %  [3] Makeig S, Debener S, Onton J, Delorme A, "Mining event-related brain dynamics"
    %      Trends in Cognitive Science, 8(5):204-210, 2004.
    %
    % Name:
    %   Regularized Spatio-Spectral Dynamics v2
    %
    %                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                                2011-01-28

    methods
        function defaults = preprocessing_defaults(self)
            defaults = { ...
                'Resampling',128, ...
                'FIRFilter',{[0.5 1],'highpass'}, ...
                'EpochExtraction',[-2 2], ...
                'ICA','beamica', ...
                ...%'DipoleFitting',{'ConfusionRange',7}, ...
                'Projection','on', ...
                'ERSPTransform',{'WindowStep',1/15,'SpectralMap','sqrt'}, ...
                };
        end

        function defaults = machine_learning_defaults(self)
            defaults = {'dal', 2.^(4:-0.25:-3), 'scaling','none'};
            %defaults = {'logreg', 'variant',{'lars','ElasticMixing',0.5}};
        end

        function model = feature_adapt(self,varargin)
            args = arg_define(varargin, ...
                    arg_norep('signal'), ...
                    arg({'spectral_prior','SpectralPrior'},@(f)1,[],'Spectral prior. Likelihood function of frequency in Hz.','guru',true), ...
                    arg({'temporal_prior','TemporalPrior'},@(t)1,[],'Temporal prior. Likelihood function of time in s.','guru',true), ...
                    arg({'spatial_prior','SpatialPrior'},@(p)1,[],'Spatial prior. Likelihood function of MNI coordinate vector.','guru',true), ...
                    arg({'anatomical_prior','AnatomicalPrior'},{'Left Cerebrum','Right Cerebrum','Left Cerebellum','Right Cerebellum','Left Brainstem','Right Brainstem','Inter-Hemispheric'},{'-- Hemispheres --','Left Cerebrum','Right Cerebrum','Left Cerebellum','Right Cerebellum','Left Brainstem','Right Brainstem','Inter-Hemispheric','-- Lobes --','Anterior Lobe','Frontal Lobe','Frontal-Temporal Space','Limbic Lobe','Medulla','Midbrain','Occipital Lobe','Parietal Lobe','Pons','Posterior Lobe','Sub-lobar','Temporal Lobe','-- Gyri --','Angular Gyrus','Anterior Cingulate','Caudate','Cerebellar Lingual','Cerebellar Tonsil','Cingulate Gyrus','Claustrum','Culmen','Culmen of Vermis','Cuneus','Declive','Declive of Vermis','Extra-Nuclear','Fastigium','Fourth Ventricle','Fusiform Gyrus','Inferior Frontal Gyrus','Inferior Occipital Gyrus','Inferior Parietal Lobule','Inferior Semi-Lunar Lobule','Inferior Temporal Gyrus','Insula','Lateral Ventricle','Lentiform Nucleus','Lingual Gyrus','Medial Frontal Gyrus','Middle Frontal Gyrus','Middle Occipital Gyrus','Middle Temporal Gyrus','Nodule','Orbital Gyrus','Paracentral Lobule','Parahippocampal Gyrus','Postcentral Gyrus','Posterior Cingulate','Precentral Gyrus','Precuneus','Pyramis','Pyramis of Vermis','Rectal Gyrus','Subcallosal Gyrus','Sub-Gyral','Superior Frontal Gyrus','Superior Occipital Gyrus','Superior Parietal Lobule','Superior Temporal Gyrus','Supramarginal Gyrus','Thalamus','Third Ventricle','Transverse Temporal Gyrus','Tuber','Tuber of Vermis','Uncus','Uvula','Uvula of Vermis'}, 'Anatomical prior. Select anatomical structures that are likely to contain processes of interest.'),...
                    arg({'vectorize_features','VectorizeFeatures'},true,[],'Vectorize feature tensors. This is for classifiers that cannot handle matrix or tensor-shaped features.'),...
                    arg({'normalize_features','NormalizeFeatures'},true,[],'Normalize time/frequency features. If enabled, features will be normalized by a rank-1 normalization matrix (rather than pixelwise).'));
%
            model = rmfield(args,'signal');
            % determine data rescaling factors
            if args.normalize_features
                tdata = permute(args.signal.data,[1 2 4 3]);
                sdata = permute(args.signal.data,[1 4 2 3]);
                for c=size(args.signal.data,1):-1:1
                    % temporal scale vector
                    temp = 1./median(abs(bsxfun(@minus,tdata(c,:,:),median(tdata(c,:,:),3))),3);
                    % spectral scale vector
                    spec = 1./median(abs(bsxfun(@minus,sdata(c,:,:),median(sdata(c,:,:),3))),3);
                    % time/freq scaling matrix
                    model.scaling(c,:,:) = sqrt(temp')*sqrt(spec);
                end
            end
            % determine model shape
            model.shape = size(args.signal.data);
            model.shape = model.shape([2 4 1]);
        end

        function [features,shape] = feature_extract(self,signal,featuremodel)
            signal.data = permute(signal.data,[1 2 4 3]);
            % optionally apply normalization
            if featuremodel.normalize_features
                features = bsxfun(@times,signal.data,featuremodel.scaling);
            else
                features = signal.data;
            end
            features = permute(features,[2 3 1 4]);
            % determine feature shape
            siz = size(features);
            shape = siz(1:3);
            % do final vectorization if desired
            if featuremodel.vectorize_features
                features = reshape(features,[],signal.trials)'; end
        end

        function layout = dialog_layout_defaults(self)
            layout = {'SignalProcessing.Resampling.SamplingRate','SignalProcessing.FIRFilter.Frequencies', '', ...
                'SignalProcessing.ICA.DataCleaning.DataSetting', '', 'SignalProcessing.ICA.Variant', ...
                'SignalProcessing.EpochExtraction','', ...
                'SignalProcessing.ERSPTransform','', ...
                'Prediction.FeatureExtraction.SpectralPrior', 'Prediction.FeatureExtraction.TemporalPrior', ...
                'Prediction.FeatureExtraction.SpatialPrior','Prediction.FeatureExtraction.AnatomicalPrior', '', ...
                'Prediction.MachineLearning.Learner.Lambdas', 'Prediction.MachineLearning.Learner.LossFunction','Prediction.MachineLearning.Learner.Regularizer'};
        end
    end
end

% TODO: reimplement priors

%         function [featuremodel,conditioningmodel,predictivemodel] = calibrate_prediction_function(self,varargin)
%             args = arg_define(varargin, ...
%                 arg_norep('signal'), ...
%                 arg_sub({'fex','FeatureExtraction'},{},...
%                     {arg({'spectral_prior','SpectralPrior'},@(f)1,[],'Spectral prior. Likelihood function of frequency in Hz.','guru',true), ...
%                      arg({'temporal_prior','TemporalPrior'},@(t)1,[],'Temporal prior. Likelihood function of time in s.','guru',true), ...
%                      arg({'spatial_prior','SpatialPrior'},@(p)1,[],'Spatial prior. Likelihood function of MNI coordinate vector.','guru',true), ...
%                      arg({'anatomical_prior','AnatomicalPrior'},{'Left Cerebrum','Right Cerebrum','Left Cerebellum','Right Cerebellum','Left Brainstem','Right Brainstem','Inter-Hemispheric'},{'-- Hemispheres --','Left Cerebrum','Right Cerebrum','Left Cerebellum','Right Cerebellum','Left Brainstem','Right Brainstem','Inter-Hemispheric','-- Lobes --','Anterior Lobe','Frontal Lobe','Frontal-Temporal Space','Limbic Lobe','Medulla','Midbrain','Occipital Lobe','Parietal Lobe','Pons','Posterior Lobe','Sub-lobar','Temporal Lobe','-- Gyri --','Angular Gyrus','Anterior Cingulate','Caudate','Cerebellar Lingual','Cerebellar Tonsil','Cingulate Gyrus','Claustrum','Culmen','Culmen of Vermis','Cuneus','Declive','Declive of Vermis','Extra-Nuclear','Fastigium','Fourth Ventricle','Fusiform Gyrus','Inferior Frontal Gyrus','Inferior Occipital Gyrus','Inferior Parietal Lobule','Inferior Semi-Lunar Lobule','Inferior Temporal Gyrus','Insula','Lateral Ventricle','Lentiform Nucleus','Lingual Gyrus','Medial Frontal Gyrus','Middle Frontal Gyrus','Middle Occipital Gyrus','Middle Temporal Gyrus','Nodule','Orbital Gyrus','Paracentral Lobule','Parahippocampal Gyrus','Postcentral Gyrus','Posterior Cingulate','Precentral Gyrus','Precuneus','Pyramis','Pyramis of Vermis','Rectal Gyrus','Subcallosal Gyrus','Sub-Gyral','Superior Frontal Gyrus','Superior Occipital Gyrus','Superior Parietal Lobule','Superior Temporal Gyrus','Supramarginal Gyrus','Thalamus','Third Ventricle','Transverse Temporal Gyrus','Tuber','Tuber of Vermis','Uncus','Uvula','Uvula of Vermis'}, 'Anatomical prior. Select anatomical structures that are likely to contain processes of interest.'),...
%                     }, 'Parameters for the feature-adaptation function. These parameters control how features are statistically adapted and extracted from the filtered data before they are passed int othe machine learning stage.'), ...
%                 arg_sub({'cond','Conditioning'},{},@self.feature_adapt_conditioning,'Feature conditioning parameters. Allows to further process features for better usability with classifiers.'), ...
%                 arg_sub({'ml','MachineLearning'},{'Learner',self.machine_learning_defaults()},@ml_train,'Machine learning stage of the paradigm. Operates on the feature vectors that are produced by the feature-extraction stage.'));
%
%             try
%                 data = self.rssd_load_overcomplete(args.signal);
%                 % read out component dipole fits
%                 if ~isempty(data.dipfit)
%                     if isfield(data.dipfit,'multimodel')
%                         dipfits = [data.dipfit.multimodel{:}];
%                     else
%                         dipfits = data.dipfit.model;
%                     end
%                 else
%                     dipfits = [];
%                 end
%                 if ~isempty(args.fex.anatomical_prior) && ~isequal(args.fex.anatomical_prior,false) && ~isempty(dipfits)
%                     % if an anatomical prior was given, we can pre-prune the potenial ERSPs
%                     ersprange = [];
%                     for k=1:size(data.icaweights,1)
%                         matches{k} = intersect(dipfits(k).structures,args.fex.anatomical_prior);
%                         if ~isempty(matches{k})
%                             ersprange(end+1) = k; end
%                     end
%                 else
%                     % otherwise we need to compute them all
%                     ersprange = 1:size(data.icaweights,1);
%                 end
%
%                 % compute ERSPs
%                 disp('Now computing time/frequency decompositions...');
%                 [T,X,freqs,times] = evalc('self.ersp_compute(data,args,ersprange)'); %#ok<ASGLU>
%                 retain_ics = [];
%                 structures = {};
%                 probabilities = [];
%                 summed_probabilities = [];
%                 for k = ersprange
%                     % get left/right normalization vectors lhs/rhs
%                     lhs = args.fex.spectral_prior(freqs) .* sum(reshape(var(X{k},[],2),size(X{k},1),[]),2).^-0.25;
%                     rhs = args.fex.temporal_prior(times) .* sum(reshape(var(X{k},[],1),[],size(X{k},3)),2).^-0.25;
%                     % make sure that they are well-behaved
%                     lhs(~isfinite(lhs)) = 0; medlhs = median(lhs); lhs(lhs>3*medlhs) = medlhs;
%                     rhs(~isfinite(rhs)) = 0; medrhs = median(rhs); rhs(rhs>3*medrhs) = medrhs;
%                     % turn into a scaling matrix
%                     prior{k} =  diag(lhs) * ones(length(freqs),length(times)) * diag(rhs);
%                     % incorporate the spatial prior
%                     if ~isempty(dipfits)
%                         prior{k} = prior{k} * args.fex.spatial_prior(dipfits(k).posxyz); end
%                 end
%                 for k = ersprange
%                     % incorporate the anatomical prior
%                     if ~isempty(args.fex.anatomical_prior) && ~isequal(args.fex.anatomical_prior,false) && ~isempty(dipfits)
%                         [matches,idx] = intersect(dipfits(k).structures,args.fex.anatomical_prior); %#ok<ASGLU>
%                         % sum the probabilities for being in each of the accepted structures (can be > 1 as the structures are highly overlapping)
%                         prior{k} = sum(dipfits(k).probabilities(idx)) * prior{k};
%                         structures{k} = dipfits(k).structures(idx);
%                         probabilities{k} = dipfits(k).probabilities(idx);
%                         summed_probabilities(k) = sum(dipfits(k).probabilities(idx));
%                         % figure; topoplot(data.icawinv(:,k),data.chanlocs(data.icachansind),'electrodes','labels'); title([hlp_tostring(xdipfits(k).structures(idx)) ' - ' hlp_tostring(dipfits(k).probabilities(idx))]);
%                     end
%                     if ~all(all(prior{k}==0))
%                         retain_ics(end+1) = k; end
%                 end
%                 featuremodel.prior = prior;
%                 featuremodel.retain_ics = retain_ics;
%                 featuremodel.args = args;
%                 blocksizes = cellfun(@size,featuremodel.prior,'UniformOutput',false);
%
%                 % update machine learning parameters
%                 args.ml.learner.shape = vertcat(blocksizes{retain_ics}); %vertcat(blocksizes{cellfun(@prod,blocksizes)});
%
%         end
%
%         function features = feature_extract(self,signal,featuremodel)
%             try
%                 data = self.rssd_load_overcomplete(signal);
%                 % compute ERSP features
%                 [T,features] = evalc('self.ersp_compute(data,featuremodel.args,featuremodel.retain_ics)'); %#ok<ASGLU>
%                 % scale each component-ERSP by the respective prior & sparsify
%                 for c=featuremodel.retain_ics
%                     features{c} = bsxfun(@times,features{c},featuremodel.prior{c}); end
%                 % generate block-compressed version of the data
%                 features = reshape([features{featuremodel.retain_ics}],[],data.trials)';
%             catch
%                 disp('Glitch in para_rssd. halting...');
%                 keyboard;
%             end
%         end

repo: REST


text:
function [] = fix_pos(blocka,blockb,r,t)

% fix the height distance between 2 blocks
% blocka acts as the origin of coordinates
% blockb is the one move around the blocka
% t is for the down, r is for the left

Pos = get_param(blocka, 'Position');
Posb = get_param(blockb, 'Position');
B = Posb(3)-Posb(1);
H = Posb(4)-Posb(2);
pos1 = Pos(1)+r;
pos3 = pos1+B;
pos2 = Pos(4)+t;
pos4 = pos2+H;
set_param(blockb, 'Position',[pos1 pos2 pos3 pos4]);

author: Sable

repo: mcbench-benchmarks


text:
function value = r8_cosh ( x )

%*****************************************************************************80
%
%% R8_COSH evaluates the hyperbolic cosine of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the hyperbolic cosine of X.
%
  persistent ymax

  if ( isempty ( ymax ) )
    ymax = 1.0 / sqrt ( r8_mach ( 3 ) );
  end

  y = exp ( abs ( x ) );

  if ( y < ymax )
    value = 0.5 * ( y + 1.0 / y );
  else
    value = 0.5 * y;
  end

  return
end


text:
function logInfo(str, varargin)
% get caller function ID and display log msg
if nargin > 1
    cmd = 'sprintf(str';
    for i = 1:length(varargin)
        cmd = sprintf('%s, varargin{%d}', cmd, i);
    end
    str = eval([cmd, ');']);
end
[st, i] = dbstack();
caller = st(2).name;
fprintf('%s: %s\n', caller, str);
end

author: kunhe

repo: FastAP-metric-learning


text:
%changeCobraSolver('tomlab_cplex','lp');
%modelClosed = modelConsistent;
clear FF R
% add demands for all metabolites in Recon
%modelClosed = model;
%modelClosed = addDemandReaction(modelClosed,modelClosed.mets);

modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
modelexchanges4 = strmatch('EX_',modelClosed.rxns);
modelexchanges2 = strmatch('DM_',modelClosed.rxns);
modelexchanges3 = strmatch('sink_',modelClosed.rxns);
selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
modelClosed.lb(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges)))=0;
%modelClosed.ub(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;

[LeakMets,modelClosed] = fastLeakTest(modelClosed, modelClosed.rxns(modelexchanges),0);

author: opencobra

repo: cobratoolbox


text:
function this = compute_stat_images(this)
% computes statistical images for time series QC, such as
%   mean, standard deviation (sd), signal-to-noise ratio (snr)
%   coefficient of variation (1/snr) images
%   difference images between last/first (drift) and odd/even (image noise)
%   volumes
%
%   Y = MrSeries()
%   Y.compute_stat_images(inputs)
%
% This is a method of class MrSeries.
%
% IN
%   parameters.compute_stat_images
%
% OUT
%   this.mean
%   this.snr
%   this.sd
%   this.coeff_var
%   this.diffLastFirst
%   this.diffOddEven
%
% EXAMPLE
%   compute_stat_images
%
%   See also MrSeries MrImage MrImage.compute_stat_image

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-06
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

this.init_processing_step('compute_stat_images');
applicationDimension = this.parameters.compute_stat_images.applicationDimension;

% compute statistical images via MrImage method and update save-parameters
[~, nameStatImageArray] = this.get_all_image_objects('stats');

for iImage = 1:numel(nameStatImageArray)
    img = nameStatImageArray{iImage};
    parameters = this.(img).parameters;
    this.(img) = this.data.compute_stat_image(img, ...
        'applicationDimension', applicationDimension);
    this.(img).name = sprintf('%s (%s)', img, this.name);
    this.(img).parameters.save = parameters.save;
end

this.finish_processing_step('compute_stat_images', this.(img));

author: translationalneuromodeling

repo: tapas


text:
function cx = getcx(X)

cx = X.cx;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function z = cvxprob( varargin )

global cvx___
if ~iscellstr( varargin ),
    error( 'Arguments must be strings.' );
end

%
% Clear out any old problems left at this depth. These will be here due to
% an error while constructing a CVX model, or due to the user deciding to
% start over when constructing a model
%

cvx_global
st = dbstack;
depth = length( st ) - 2;
if length(st) <= 2,
    name = '';
else
    name = st(3).name;
    if cvx___.mversion < 7,
        name = name(max(find(name==filesep))+1:end-2);
    end
end
if ~isempty( cvx___.problems ),
    ndx = find( [ cvx___.problems.depth ] >= depth );
    if ~isempty( ndx ),
        temp = cvx___.problems(ndx(1));
        if temp.depth == depth && ( ~isempty(temp.objective) || ~isempty(temp.variables) || ~isempty(temp.duals) || nnz(temp.t_variable) > 1 );
            warning( 'CVX:Empty', 'A non-empty cvx problem already exists in this scope.\n   It is being overwritten.', 1 ); %#ok
        end
        cvx_pop( temp.self, 'reset' );
    end
end

%
% Grab the latest defaults to place in the new problem
%

if ~isempty( cvx___.problems ),
    nprec  = cvx___.problems( end ).precision;
    ngprec = cvx___.problems( end ).gptol;
    nrprec = cvx___.problems( end ).rat_growth;
    nsolv  = cvx___.problems( end ).solver;
    nquiet = cvx___.problems( end ).quiet;
else
    nprec  = cvx___.precision;
    ngprec = cvx___.gptol;
    nrprec = cvx___.rat_growth;
    nsolv  = cvx___.solver;
    nquiet = cvx___.quiet;
end

%
% Construct the object
%

z = class( struct( 'index_', length( cvx___.problems ) + 1 ), 'cvxprob', cvxobj );
temp = struct( ...
    'name',          name,   ...
    'complete',      true,   ...
    'sdp',           false,  ...
    'gp',            false,  ...
    'separable',     false,  ...
    'locked',        false,  ...
    'precision',     nprec,  ...
    'solver',        nsolv,  ...
    'gptol',         ngprec, ...
    'quiet',         nquiet, ...
    'cputime',       cputime, ...
    'rat_growth',    nrprec, ...
    't_variable',    logical( sparse( length( cvx___.reserved ), 1 ) ), ...
    'n_equality',    length(cvx___.equalities), ...
    'n_linform',     length(cvx___.linforms), ...
    'n_uniform',     length(cvx___.uniforms), ...
    'variables',     [],         ...
    'duals',         [],         ...
    'dvars',         [],         ...
    'direction',     '',         ...
    'geometric',     [],         ...
    'objective',     [],         ...
    'status',        'unsolved', ...
    'result',        [],         ...
    'depth',         depth, ...
    'self',          z );
temp.t_variable( 1 ) = true;

%
% Process the argument strings
%

for k = 1 : nargin,
    mode = varargin{k};
    switch lower( mode ),
        case 'quiet',
            temp.quiet = true;
        case 'set',
            temp.complete  = false;
            temp.direction = 'find';
        case 'sdp',
            if temp.gp,
                error( 'The GP and SDP modifiers cannot be used together.' );
            end
            temp.sdp = true;
        case 'gp',
            if temp.sdp,
                error( 'The GP and SDP modifiers cannot be used together.' );
            end
            temp.gp = true;
            if cvx___.expert == 0,
                cvx___.expert = -1;
            end
        case 'separable',
            temp.separable = true;
        otherwise,
            error( 'Invalid CVX problem modifier: %s', mode );
    end
end

%
% Add the problem to the stack
%

if isempty( cvx___.problems ),
    cvx___.problems = temp;
    cvx_setpath(1);
    if cvx___.profile,
        profile resume
    end
else
    cvx___.problems( end + 1 ) = temp;
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function u = bsv ( a, b, alpha, beta, nu, n, output )

%*****************************************************************************80
%
%% BSV applies Newton's method to a discretized steady viscous Burgers equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the left and right endpoints.
%
%    Input, real ALPHA, BETA, the Dirichlet boundary values at the left
%    and right.
%
%    Input, real NU, the viscosity.  Normally, 0 < NU.
%
%    Input, integer N, the number of nodes to use between A and B.
%
%    Input, logical OUTPUT, is TRUE if printout is desired.
%
%    Output, real U(N), the computed discretized solution.
%
  if ( n < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BSV - Fatal error!\n' );
    fprintf ( 1, '  N < 2.\n' );
    error ( 'BSV - Fatal error!\n' );
  end
%
%  Set some iteration parameters.
%
  newton_step_max = 50;
  newton_resid_tol = sqrt ( eps );
  newton_step_tol = sqrt ( eps );
%
%  Use equally spaced nodes from A to B.
%
  dx = ( b - a ) / ( n - 1 );
%
%  The initial guess will be the linear interpolant to the boundary conditions.
%
  u = ( linspace ( alpha, beta, n ) )';
%
%  Prepare for the Newton iteration.
%
  J = sparse ( [], [], [], n, n, 3 * n );
  newton_step = 0;

  if ( output )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Step  ||F(U)||\n' );
    fprintf ( 1, '\n' );
  end

  while ( newton_step <= newton_step_max )

    f(1,1) = u(1) - alpha;

%   for i = 2 : n - 1
%     f(i,1) = 0.5 * ( u(i+1)^2 - u(i-1)^2 ) / ( 2.0 * dx ) ...
%          - nu * ( u(i+1) - 2.0 * u(i) + u(i-1) ) / ( dx^2 );
%   end

    f(2:n-1,1) = 0.5 * ( u(3:n).^2 - u(1:n-2).^2 ) / ( 2.0 * dx ) ...
           - nu * ( u(3:n) - 2.0 * u(2:n-1) + u(1:n-2) ) / ( dx^2 );

    f(n,1) = u(n) - beta;

    f_norm = norm ( f, inf );

    if ( output )
      fprintf ( 1, '  %4d  %g\n', newton_step, f_norm );
    end

    if ( f_norm < newton_resid_tol )
      break
    end
%
%  Define the Jacobian matrix.
%
    J(1,1) = 1.0;
    for i = 2 : n - 1
      J(i,i-1) = - 2.0 * u(i-1) / ( 4.0 * dx )       - nu / dx^2;
      J(i,i)   =                                 2.0 * nu / dx^2;
      J(i,i+1) =   2.0 * u(i+1) / ( 4.0 * dx )       - nu / dx^2;
    end
    J(n,n) = 1.0;
%
%  Solve the linear system J * du = -f to get the Newton update.
%
    du = J \ f;
    du_norm = norm ( du, inf );

    u_norm = norm ( u, inf );
    if ( du_norm < newton_step_tol * ( u_norm + 1.0 ) )
      break
    end

    u = u - du;
    newton_step = newton_step + 1;

  end

  if ( newton_step_max < newton_step )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BSV - Warning!\n' );
    fprintf ( 1, '  The Newton iteration did not converge.\n' );
  end

  return
end


text:
classdef TestBlur
    %TestBlur

    properties (Constant)
        img = 255*uint8([...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 1 1 1 0 0 0 0;...
            0 0 0 1 1 1 0 0 0 0;...
            0 0 0 1 1 1 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
        ]);
    end

    methods (Static)
        function test_1
            result = cv.blur(TestBlur.img);
            validateattributes(result, {class(TestBlur.img)}, ...
                {'size',size(TestBlur.img)});
        end

        function test_2
            result = cv.blur(TestBlur.img, 'KSize', [3,7]);
            validateattributes(result, {class(TestBlur.img)}, ...
                {'size',size(TestBlur.img)});
        end

        function test_3
            result = cv.blur(TestBlur.img, 'Anchor', [0,1]);
            validateattributes(result, {class(TestBlur.img)}, ...
                {'size',size(TestBlur.img)});
        end

        function test_4
            result = cv.blur(TestBlur.img, 'BorderType', 'Constant');
            validateattributes(result, {class(TestBlur.img)}, ...
                {'size',size(TestBlur.img)});
        end

        function test_5
            img = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
            result = cv.blur(img);
            validateattributes(result, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.blur();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [traj, infStates] = tapas_hgf_binary_pu_tbt(r, p, varargin)
% Calculates the trajectories of the agent's representations under the HGF
%
% This function can be called in two ways:
%
% (1) tapas_hgf_binary_pu_tbt(r, p)
%
%     where r is the structure generated by tapas_fitModel and p is the parameter vector in native space;
%
% (2) tapas_hgf_binary_pu_tbt(r, ptrans, 'trans')
%
%     where r is the structure generated by tapas_fitModel, ptrans is the parameter vector in
%     transformed space, and 'trans' is a flag indicating this.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016-2017 Rebecca Lawson, Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform paramaters back to their native space if needed
if ~isempty(varargin) && strcmp(varargin{1},'trans');
    p = tapas_hgf_binary_pu_tbt_transp(r, p);
end

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = (length(p)+1)/5;

    if l ~= floor(l)
        error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
    end
end

% Unpack parameters
mu_0 = p(1:l);
sa_0 = p(l+1:2*l);
rho  = p(2*l+1:3*l);
ka   = p(3*l+1:4*l-1);
om   = p(4*l:5*l-2);
th   = exp(p(5*l-1));
eta0 = p(5*l);
eta1 = p(5*l+1);

% Add dummy "zeroth" trial
u = [0; r.u(:,1)];
al= [0; r.u(:,2)];

% Number of trials (including prior)
n = length(u);

% Assume that if u has more than one column, the last contains t
try
    if r.c_prc.irregular_intervals
        if size(u,2) > 1
            t = [0; r.u(:,end)];
        else
            error('tapas:hgf:InputSingleColumn', 'Input matrix must contain more than one column if irregular_intervals is set to true.');
        end
    else
        t = ones(n,1);
    end
catch
    if size(u,2) > 1
        t = [0; r.u(:,end)];
    else
        t = ones(n,1);
    end
end

% Initialize updated quantities

% Representations
mu = NaN(n,l);
pi = NaN(n,l);

% Other quantities
muhat = NaN(n,l);
pihat = NaN(n,l);
v     = NaN(n,l);
w     = NaN(n,l-1);
da    = NaN(n,l);

% Representation priors
% Note: first entries of the other quantities remain
% NaN because they are undefined and are thrown away
% at the end; their presence simply leads to consistent
% trial indices.
mu(1,1) = tapas_sgm(mu_0(1), 1);
pi(1,1) = Inf;
mu(1,2:end) = mu_0(2:end);
pi(1,2:end) = 1./sa_0(2:end);

% Pass through representation update loop
for k = 2:1:n
    if not(ismember(k-1, r.ign))

        %%%%%%%%%%%%%%%%%%%%%%
        % Effect of input u(k)
        %%%%%%%%%%%%%%%%%%%%%%

        % 2nd level prediction
        muhat(k,2) = mu(k-1,2) +t(k) *rho(2);

        % 1st level
        % ~~~~~~~~~
        % Prediction
        muhat(k,1) = tapas_sgm(ka(1) *muhat(k,2), 1);

        % Precision of prediction
        pihat(k,1) = 1/(muhat(k,1)*(1 -muhat(k,1)));

        % Mean update
        mu(k,1) = u(k);
        und1 = exp(-(u(k) -eta1)^2/(2*al(k)));
        und0 = exp(-(u(k) -eta0)^2/(2*al(k)));
        mu(k,1) = muhat(k,1) *und1 /(muhat(k,1) *und1 +(1 -muhat(k,1)) *und0);

        % Prediction error
        da(k,1) = mu(k,1) -muhat(k,1);

        % 2nd level
        % ~~~~~~~~~
        % Prediction: see above

        % Precision of prediction
        pihat(k,2) = 1/(1/pi(k-1,2) +exp(ka(2) *mu(k-1,3) +om(2)));

        % Updates
        pi(k,2) = pihat(k,2) +ka(1)^2/pihat(k,1);
        mu(k,2) = muhat(k,2) +ka(1)/pi(k,2) *da(k,1);

        % Implied posterior precision at first level
        sgmmu2 = tapas_sgm(ka(1) *mu(k,2), 1);
        pi(k,1) = pi(k,2)/(sgmmu2*(1-sgmmu2));

        % Volatility prediction error
        da(k,2) = (1/pi(k,2) +(mu(k,2) -muhat(k,2))^2) *pihat(k,2) -1;

        if l > 3
            % Pass through higher levels
            % ~~~~~~~~~~~~~~~~~~~~~~~~~~
            for j = 3:l-1
                % Prediction
                muhat(k,j) = mu(k-1,j) +t(k) *rho(j);

                % Precision of prediction
                pihat(k,j) = 1/(1/pi(k-1,j) +t(k) *exp(ka(j) *mu(k-1,j+1) +om(j)));

                % Weighting factor
                v(k,j-1) = t(k) *exp(ka(j-1) *mu(k-1,j) +om(j-1));
                w(k,j-1) = v(k,j-1) *pihat(k,j-1);

                % Updates
                pi(k,j) = pihat(k,j) +1/2 *ka(j-1)^2 *w(k,j-1) *(w(k,j-1) +(2 *w(k,j-1) -1) *da(k,j-1));

                if pi(k,j) <= 0
                    error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
                end

                mu(k,j) = muhat(k,j) +1/2 *1/pi(k,j) *ka(j-1) *w(k,j-1) *da(k,j-1);

                % Volatility prediction error
                da(k,j) = (1/pi(k,j) +(mu(k,j) -muhat(k,j))^2) *pihat(k,j) -1;
            end
        end

        % Last level
        % ~~~~~~~~~~
        % Prediction
        muhat(k,l) = mu(k-1,l) +t(k) *rho(l);

        % Precision of prediction
        pihat(k,l) = 1/(1/pi(k-1,l) +t(k) *th);

        % Weighting factor
        v(k,l)   = t(k) *th;
        v(k,l-1) = t(k) *exp(ka(l-1) *mu(k-1,l) +om(l-1));
        w(k,l-1) = v(k,l-1) *pihat(k,l-1);

        % Updates
        pi(k,l) = pihat(k,l) +1/2 *ka(l-1)^2 *w(k,l-1) *(w(k,l-1) +(2 *w(k,l-1) -1) *da(k,l-1));

        if pi(k,l) <= 0
            error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
        end

        mu(k,l) = muhat(k,l) +1/2 *1/pi(k,l) *ka(l-1) *w(k,l-1) *da(k,l-1);

        % Volatility prediction error
        da(k,l) = (1/pi(k,l) +(mu(k,l) -muhat(k,l))^2) *pihat(k,l) -1;
    else

        mu(k,:) = mu(k-1,:);
        pi(k,:) = pi(k-1,:);

        muhat(k,:) = muhat(k-1,:);
        pihat(k,:) = pihat(k-1,:);

        v(k,:)  = v(k-1,:);
        w(k,:)  = w(k-1,:);
        da(k,:) = da(k-1,:);

    end
end

% Implied learning rate at the first level
sgmmu2 = tapas_sgm(ka(1) *mu(:,2), 1);
lr1    = diff(sgmmu2)./da(2:n,1);
lr1(da(2:n,1)==0) = 0;

% Remove representation priors
mu(1,:)  = [];
pi(1,:)  = [];

% Check validity of trajectories
if any(isnan(mu(:))) || any(isnan(pi(:)))
    error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
else
    % Check for implausible jumps in trajectories
    dmu = diff(mu(:,2:end));
    dpi = diff(pi(:,2:end));
    rmdmu = repmat(sqrt(mean(dmu.^2)),length(dmu),1);
    rmdpi = repmat(sqrt(mean(dpi.^2)),length(dpi),1);

    jumpTol = 16;
    if any(abs(dmu(:)) > jumpTol*rmdmu(:)) || any(abs(dpi(:)) > jumpTol*rmdpi(:))
        error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
    end
end

% Remove other dummy initial values
muhat(1,:) = [];
pihat(1,:) = [];
v(1,:)     = [];
w(1,:)     = [];
da(1,:)    = [];

% Create result data structure
traj = struct;

traj.mu     = mu;
traj.sa     = 1./pi;

traj.muhat  = muhat;
traj.sahat  = 1./pihat;

traj.v      = v;
traj.w      = w;
traj.da     = da;

% Updates with respect to prediction
traj.ud = mu -muhat;

% Psi (precision weights on prediction errors)
psi        = NaN(n-1,l);
psi(:,2)   = 1./pi(:,2);
psi(:,3:l) = pihat(:,2:l-1)./pi(:,3:l);
traj.psi   = psi;

% Epsilons (precision-weighted prediction errors)
epsi        = NaN(n-1,l);
epsi(:,2:l) = psi(:,2:l) .*da(:,1:l-1);
traj.epsi   = epsi;

% Full learning rate (full weights on prediction errors)
wt        = NaN(n-1,l);
wt(:,1)   = lr1;
wt(:,2)   = psi(:,2);
wt(:,3:l) = 1/2 *(v(:,2:l-1) *diag(ka(2:l-1))) .*psi(:,3:l);
traj.wt   = wt;

% Create matrices for use by the observation model
infStates = NaN(n-1,l,4);
infStates(:,:,1) = traj.muhat;
infStates(:,:,2) = traj.sahat;
infStates(:,:,3) = traj.mu;
infStates(:,:,4) = traj.sa;

return;


text:
classdef ObjectiveFunction < handle

    properties (Access = public)
        value
        gradient
    end

    properties (Access = protected)
        cost
        constraint
        dualVariable
    end

    properties (Access = private)
        valueOld
    end

    methods (Access = public, Static)

        function obj = create(cParams)
            f = ObjectiveFunctionFactory();
            obj = f.create(cParams);
        end

    end

    methods (Access = protected)

        function init(obj,cParams)
            obj.cost         = cParams.cost;
            obj.constraint   = cParams.constraint;
            obj.dualVariable = cParams.dualVariable;
        end

    end

    methods (Access = public)

        function incr = computeIncrement(obj)
            v  = obj.value;
            vI = obj.valueOld;
            incr = (v - vI)/abs(vI);
        end

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

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

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [V,F,bi,bo] = annulus(s,r,varargin)
  % ANNULUS Construct a triangle mesh of a unit annulus.
  %
  % [V,F] = annulus(s,r)
  % [V,F] = annulus(s,r,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   s  number of samples on the inner ring
  %   r  radius of the inner ring
  %   Optional:
  %     'Flags'  followed by flags to pass to Triangle for meshing.
  %        {'-q30 -aX'} where X is squared boundary edge length
  %     'R'  followed by outer ring radius {1}
  % Outputs:
  %   V  #V by 2 list of mesh vertex positions
  %   F  #F by 3 list of triangle mesh indices
  %   bi  list of vertices on inner boundary
  %   bo  list of vertices on outer boundary
  %

  flags = nan;
  R = 1;
  params_to_variables = containers.Map( ...
    {'Flags','R'},{'flags','R'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  if isnan(flags)
    flags = sprintf('-q30 -a%0.17f',(2*pi*r/s)^2);
  end

  theta = linspace(0,2*pi,s+1)';theta = theta(1:end-1);
  Vr = r*[cos(theta) sin(theta)];
  Er = fliplr([1:size(Vr,1);2:size(Vr,1) 1]');
  theta = linspace(0,2*pi,ceil(R/r)*s+1)';theta = theta(1:end-1);
  VR = R*[cos(theta) sin(theta)];
  ER = [1:size(VR,1);2:size(VR,1) 1]';

  [V,F] = triangulate([Vr;VR],[Er;size(Vr,1)+ER],'Holes',[0 0],'Flags',flags);
  b = unique(outline(F));
  bi = intersect(find(normrow(V)< 0.5*(R+r)),b);
  bo = intersect(find(normrow(V)> 0.5*(R+r)),b);
end

author: alecjacobson

repo: gptoolbox


text:
function pass = test_mldivide(pref)

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

% Generate a few random points to use as test values.
seedRNG(6178);
xr = 2 * rand(100, 1) - 1;

T = restrict(chebpoly(0:3), [-1 -0.5 0 0.5 1]);
L = restrict(legpoly(0:3), [-1 0 1]);

%% Scalar A:
pass(1) = normest(2\T - .5*T) < 10*eps;

%% Numeric A:
B = T.';
A = (1:4).';
x = A\B;
x0 = feval(x, 0);
x0_true = -1/15;
pass(2) = abs(x0 - x0_true) < 10*eps;

A = eye(4);
X = A\B;
X0 = feval(X, 0);
pass(3) = norm(X0 - [1 0 -1 0].') < 10*eps;

%% Transposed A:
A = T.';
B = (1:4).';
X = A\B;
X0 = feval(X, 0);
X0_true = -2.625;
pass(4) = abs(X0 - X0_true) < 100*eps;

%% Else

X = T\L;
C = diag([1 1 .75 .625]);
C(1, 3) = .25;
C(2, 4) = .375;
pass(5) = norm(X - C, inf) < 10*eps;

%% Test a bug from #437

x = chebfun('x');
A = [];
for j = 0:6
   xj = -1 + j/4;
   A = [ A , max(0, 1-4*abs(x-xj)) ];
end
u = A\x;
expected = [...
  -0.999851249504164
  -0.750297500991670
  -0.498958746529156
  -0.253867512891710
   0.014428798095994
   0.196152320507735
   0.700961919873066];
err = norm(u - expected, inf);
pass(6) = err < 1e3*eps;

%% Test on SINGFUN:

% scalar * [1 x INF] = [1 x INF] => scalar\row SINGFUN:
f = chebfun(@(x)sin(20*x)./(x+1), 'exps', [-1 0], 'splitting', 'on');
f = f.';
g = 3\f;
op = @(x) sin(20*x)./(3*(x+1));
g_vals = feval(g, xr);
g_exact = op(xr).';
err = g_vals - g_exact;
pass(7) = norm(err, inf) < 1e5*vscale(g)*eps;

% [1 x INF] * [INF x 1] = scalar => row SINGFUN\scalar:
f = chebfun(@(x)(sin(100*x).^2+1)./((x+1).^0.25), 'exps', [-0.25 0], 'splitting', 'on');
f = f.';
g = f\3;
err = f*g - 3;
pass(8) = abs(err) < 10*vscale(g)*eps;

% [INF x 1] * SCALAR = [INF x 1] => column SINGFUN\column SINGFUN:

f = chebfun(@(x)3*(x.^2+3)./(x+1).^0.4, 'exps', [-0.4 0], 'splitting', 'on');
g = chebfun(@(x)(x.^2+3)./(x+1).^0.4, 'exps', [-0.4 0], 'splitting', 'on');
h = f\g;
err = h - 1/3;
pass(9) = norm(err, inf) < vscale(f)*eps;

% [TODO]: Revive the following test.

%% Test for function defined on unbounded domain:

% % Set the domain:
% dom = [-Inf 3*pi];
% domCheck = [-1e6 3*pi];
%
% % Generate a few random points to use as test values:
% x = diff(domCheck) * rand(100, 1) + domCheck(1);
%
% % A*X = B, where both A and B are UNBNDFUNs ==> X = A\B
%
% opA = @(x) [exp(x) x.*exp(x) (1-exp(x))./x];
% A = chebfun(opA, dom);
% opB = @(x) [(2*x+1).*exp(x) exp(x) 2*(1-exp(x))./x];
% B = chebfun(opB, dom);
% X = A\B;
% res = A*X - B;
% err = feval(res, x);
% pass(10) = norm(err(:), inf) < max([eps*get(A,'vscale') ...
%     eps*get(B,'vscale')]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%  Return orientation of the current image:
%	orient is orientation 1x3 matrix, in that:
%	Three elements represent: [x y z]
%	Element value:	1 - Left to Right; 2 - Posterior to Anterior;
%			3 - Inferior to Superior; 4 - Right to Left;
%			5 - Anterior to Posterior; 6 - Superior to Inferior;
%  e.g.:
%	Standard RAS Orientation:	[1 2 3]
%	Standard RHOS Orientation:	[2 4 3]

%  Jimmy Shen (jimmy@rotman-baycrest.on.ca), 26-APR-04
%
function orient = rri_orient_ui(varargin)

   if nargin == 0
      init;
      orient_ui_fig = gcf;
      uiwait;					% wait for user finish

      orient = getappdata(gcf, 'orient');

      if isempty(orient)
         orient = [1 2 3];
      end

      if ishandle(orient_ui_fig)
         close(gcf);
      end

      return;
   end

   action = varargin{1};

   if strcmp(action, 'done')
      click_done;
   elseif strcmp(action, 'cancel')
      uiresume;
   end

   return;						% rri_orient_ui

%----------------------------------------------------------------------
function init

   save_setting_status = 'on';
   rri_orient_pos = [];

   try
      load('pls_profile');
   catch
   end

   try
      load('rri_pos_profile');
   catch
   end

   if ~isempty(rri_orient_pos) & strcmp(save_setting_status,'on')

      pos = rri_orient_pos;

   else

      w = 0.35;
      h = 0.4;
      x = (1-w)/2;
      y = (1-h)/2;

      pos = [x y w h];

   end

   handles.figure = figure('Color',[0.8 0.8 0.8], ...
        'Units','normal', ...
        'Name', 'Convert to standard RAS orientation', ...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'Position',pos, ...
        'WindowStyle', 'normal', ...
        'ToolBar','none');

   h0 = handles.figure;
   Font.FontUnits  = 'point';
   Font.FontSize   = 12;

   margin = .1;
   line_num = 6;
   line_ht = (1 - margin*2) / line_num;

   x = margin;
   y = 1 - margin - line_ht;
   w = 1 - margin * 2;
   h = line_ht * .7;

   pos = [x y w h];

   handles.Ttit = uicontrol('parent', h0, ...
	'style','text', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'background', [0.8 0.8 0.8], ...
	'string', 'Please input orientation of the current image:');

   y = y - line_ht;
   w = .2;

   pos = [x y w h];

   handles.Tx_orient = uicontrol('parent', h0, ...
	'style','text', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'background', [0.8 0.8 0.8], ...
	'string', 'X   Axes:');

   y = y - line_ht;

   pos = [x y w h];

   handles.Ty_orient = uicontrol('parent', h0, ...
	'style','text', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'background', [0.8 0.8 0.8], ...
	'string', 'Y   Axes:');

   y = y - line_ht;

   pos = [x y w h];

   handles.Tz_orient = uicontrol('parent', h0, ...
	'style','text', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'background', [0.8 0.8 0.8], ...
	'string', 'Z   Axes:');

   choice = {	'From Left to Right', 'From Posterior to Anterior', ...
		'From Inferior to Superior', 'From Right to Left', ...
		'From Anterior to Posterior', 'From Superior to Inferior'	};

   y = 1 - margin - line_ht;
   y = y - line_ht;
   w = 1 - margin - x - w;
   x = 1 - margin - w;

   pos = [x y w h];

   handles.x_orient = uicontrol('parent', h0, ...
	'style','popupmenu', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'string', choice, ...
	'value', 1, ...
	'background', [1 1 1]);

   y = y - line_ht;

   pos = [x y w h];

   handles.y_orient = uicontrol('parent', h0, ...
	'style','popupmenu', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'string', choice, ...
	'value', 2, ...
	'background', [1 1 1]);

   y = y - line_ht;

   pos = [x y w h];

   handles.z_orient = uicontrol('parent', h0, ...
	'style','popupmenu', ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','left',...
	'string', choice, ...
	'value', 3, ...
	'background', [1 1 1]);

   x = margin;
   y = y - line_ht * 1.5;
   w = .3;

   pos = [x y w h];

   handles.done = uicontrol('parent', h0, ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','center',...
        'callback', 'rri_orient_ui(''done'');', ...
	'string', 'Done');

   x = 1 - margin - w;

   pos = [x y w h];

   handles.cancel = uicontrol('parent', h0, ...
	'unit', 'normal', ...
	Font, ...
	'Position',pos, ...
	'HorizontalAlignment','center',...
        'callback', 'rri_orient_ui(''cancel'');', ...
	'string', 'Cancel');

   setappdata(h0, 'handles', handles);
   setappdata(h0, 'orient', [1 2 3]);

   return;						% init

%----------------------------------------------------------------------
function click_done

   handles = getappdata(gcf, 'handles');

   x_orient = get(handles.x_orient, 'value');
   y_orient = get(handles.y_orient, 'value');
   z_orient = get(handles.z_orient, 'value');

   orient = [x_orient y_orient z_orient];
   test_orient = [orient, orient + 3];
   test_orient = mod(test_orient, 3);

   if length(unique(test_orient)) ~= 3
      msgbox('Please don''t choose same or opposite direction','Error','modal');
      return;
   end

   setappdata(gcf, 'orient', [x_orient y_orient z_orient]);
   uiresume;

   return;						% click_done

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function value = ball_unit_volume_3d ( )

%*****************************************************************************80
%
%% BALL_UNIT_VOLUME_3D computes the volume of a unit ball in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that
%
%      X**2 + Y**2 + Z**2 <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real BALL_UNIT_VOLUME_3D, the volume of the ball.
%
  value = ( 4.0E+00 / 3.0E+00 ) * pi;

  return
end


text:
function computeAllModelChoice_STS(pathWORK,fSetNameType,outcome,freedomMat,maxOrder,nBoot)
% -------------------------------------------------------------------------
% function computeAllModelChoice_STS(pathWORK,fSetNameType,outcome,freedomMat,maxOrder,nBoot)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection for a given feature set type
% and for all experiments with different degrees of freedom. See ref. [1]
% for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the STS WORKSPACE directory.
% - fSetNameType: String specifying the name of the type of feature set
%                 (e.g., 'PET', 'SEPARATE', 'FUSED', etc.)
% - outcome: Column vector of size [nInst X 1] specifying the outcome status
%            (1 or 0) for all instances.
% - freedomMat:  Matrix of row vectors of 1's and 0's to specify the degree
%                of freedom on texture extraction parameters for all
%                experiments. For example, for an ith experiment where
%                extraction parameters 1, 2 and 4 in paramAll are allowed
%                to vary, use freedomMat(i,:) = [1,1,0,1].
% - maxOrder: Integer specifying the maximal model order to construct.
% - nBoot: Number of bootstrap samples to use.
% -------------------------------------------------------------------------
% OUTPUTS: Mutivariable models are saved in a folder named 'MODELS' in the
%          STS WORKSPACE.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd([pathWORK,'/FSET']), pathFSET = pwd;
cd([pathWORK,'/MODELS']), pathModels = pwd;

nParamType = size(freedomMat,2);
nFreedom = size(freedomMat,1);
tStart = tic;
for i = 1:nFreedom
    cd(pathFSET)
    nameOpen=['FSET_',fSetNameType,'_'];
    for j = 1:nParamType
        nameOpen = [nameOpen,num2str(freedomMat(i,j))];
    end
    fSet = load(nameOpen); fSet = struct2cell(fSet); fSet = fSet{1};
    fprintf(['SELECTING FEATURES (MODEL ORDERS OF 1 to % u) FOR ',nameOpen,' ... '],maxOrder)
    tic
    [models] = featureSelection_STS(fSet.Data,outcome,maxOrder,nBoot,fSet.Info,'IABR');
    toc
    cd(pathModels)
    save(['MODELS',nameOpen(5:end)],'models')
end
time = toc(tStart);
fprintf('TOTAL TIME: %.2f seconds\n',time)

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
%COPYFIG Create a copy of a figure, without changing the figure
%
% Examples:
%   fh_new = copyfig(fh_old)
%
% This function will create a copy of a figure, but not change the figure,
% as copyobj sometimes does, e.g. by changing legends.
%
% IN:
%    fh_old - The handle of the figure to be copied. Default: gcf.
%
% OUT:
%    fh_new - The handle of the created figure.

% Copyright (C) Oliver Woodford 2012

% 26/02/15: If temp dir is not writable, use the dest folder for temp
%           destination files (Javier Paredes)

function fh = copyfig(fh)
% Set the default
if nargin == 0
    fh = gcf;
end
% Is there a legend?
if isempty(findall(fh, 'Type', 'axes', 'Tag', 'legend'))
    % Safe to copy using copyobj
    fh = copyobj(fh, 0);
else
    % copyobj will change the figure, so save and then load it instead
    tmp_nam = [tempname '.fig'];
    try
        % Ensure that the temp dir is writable (Javier Paredes 26/2/15)
        fid = fopen(tmp_nam,'w');
        fwrite(fid,1);
        fclose(fid);
        delete(tmp_nam);  % cleanup
    catch
        % Temp dir is not writable, so use the current folder
        [dummy,fname,fext] = fileparts(tmp_nam); %#ok<ASGLU>
        fpath = pwd;
        tmp_nam = fullfile(fpath,[fname fext]);
    end
    hgsave(fh, tmp_nam);
    fh = hgload(tmp_nam);
    delete(tmp_nam);
end
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%LINEARPOLAR  Remaps an image to polar coordinates space
%
%     dst = cv.linearPolar(src, center, maxRadius)
%     dst = cv.linearPolar(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ Source image.
% * __center__ The transformation center.
% * __maxRadius__ The radius of the bounding circle to transform. It
%   determines the inverse magnitude scale parameter too.
%
% ## Output
% * __dst__ Destination image. It will have same size and type as `src`.
%
% ## Options
% * __Interpolation__ Interpolation method, default 'Linear'. One of:
%   * __Nearest__ nearest neighbor interpolation
%   * __Linear__ bilinear interpolation
%   * __Cubic__ bicubic interpolation
%   * __Lanczos4__ Lanczos interpolation over 8x8 neighborhood
% * __FillOutliers__ flag, fills all of the destination image pixels. If some
%   of them correspond to outliers in the source image, they are set to zero.
%   default true
% * __InverseMap__ flag, inverse transformation, default false. For example,
%   polar transforms:
%   * flag is not set: Forward transformation `dst(rho,phi) = src(x,y)`
%   * flag is set: Inverse transformation `dst(x,y) = src(rho,phi)`
%
% Transform the source image using the following transformation:
%
%     dst(rho,phi) = src(x,y)
%     size(dst) <- size(src)
%
% where:
%
%     I = (dx,dy) = (x-center(2), y-center(1))
%     rho = Kx * magnitude(I)
%     phi = Ky * angle(I)_{0..360 deg}
%
% and:
%
%     Kx = size(src,2) / maxRadius
%     Ky = size(src,1) / 360
%
% Polar remaps reference:
%
% ![image](https://docs.opencv.org/3.3.1/polar_remap_doc.png)
%
% Note: To calculate magnitude and angle in degrees, cv.cartToPolar is used
% internally thus angles are measured from 0 to 360 with accuracy about 0.3
% degrees.
%
% See also: cv.logPolar, cv.remap
%


text:
function disp(sncl)
% DISP - scnl disp overloaded operator

% VERSION: 1.0 of filter objects
% AUTHOR: Celso Reyes
% LASTUPDATE: 1/30/2007
if numel(sncl) > 1;
    disp(' ');
    DispStr = sprintf('%d',size(sncl,1));
    for n = 2 : numel(size(sncl))
        DispStr = sprintf('%sx%d', DispStr, size(sncl,n));
    end
    disp(sprintf('%s %s object with fields:', DispStr, class(sncl)));
    disp('    station');
    disp('    channel');
    disp('    network');
    disp('   location');
else
    disp(['   station: ' sncl.tag.station]);
    disp(['   channel: ' sncl.tag.channel]);
    disp(['   network: ' sncl.tag.network]);
    disp(['  location: ' sncl.tag.location]);
end;

author: geoscience-community-codes

repo: GISMO


text:
% Map the numerical elements in the vector "v" onto the variables "s" which can
% be of any type. The number of numerical elements must match; on exit "v"
% should be empty. Non-numerical entries are just copied. See also unwrap.m.

function [s v] = rewrap(s, v)

if isnumeric(s)
  if numel(v) < numel(s)
    error('The vector for conversion contains too few elements')
  end
  s = reshape(v(1:numel(s)), size(s));            % numeric values are reshaped
  v = v(numel(s)+1:end);                        % remaining arguments passed on
elseif isstruct(s)
  [s p] = orderfields(s); p(p) = 1:numel(p);      % alphabetize, store ordering
  [t v] = rewrap(struct2cell(s), v);                 % convert to cell, recurse
  s = orderfields(cell2struct(t,fieldnames(s),1),p);  % conv to struct, reorder
elseif iscell(s)
  for i = 1:numel(s)             % cell array elements are handled sequentially
    [s{i} v] = rewrap(s{i}, v);
  end
end                                             % other types are not processed

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [node,elem,HB] = cubegradmesh(xh,yh,zh)
%% CUBEGRADMESH a uniform mesh of a cube
%
%
% Example
%     xh = 0:0.125:1;
%     yh = 0:0.125:1;
%     zh = 0:0.125:0.5;
%     [node,elem] = cubegradmesh(xh,yh,zh);
%     showmesh3(node,elem);
%
%   xh = (0:0.125:1).^2;
%   yh = (0:0.125:1).^2;
%   zh = 0:0.125:0.5;
%   [node,elem] = cubegradmesh(xh,yh,zh);
%   showmesh3(node,elem);
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%       8 --- 7
%      /|    /|
%     5 --- 6 |      z
%     | 4 --| 3      |  y
%     |/    |/       | /
%     1 --  2        o --- x
%
% The order of vertices is important for the uniform refinement.

[x,y,z] = ndgrid(xh,yh,zh);
node = [x(:),y(:),z(:)];
[nx,ny,nz] = size(x);
elem = zeros(6*(nx-1)*(ny-1)*(nz-1),4);
indexMap = reshape(1:nx*ny*nz,nx,ny,nz);
localIndex = zeros(8,1);
idx = 1;
for k = 1:nz-1
    for j = 1:ny-1
        for i = 1:nx-1
            localIndex(1) = indexMap(i,j,k);
            localIndex(2) = indexMap(i+1,j,k);
            localIndex(3) = indexMap(i+1,j+1,k);
            localIndex(4) = indexMap(i,j+1,k);
            localIndex(5) = indexMap(i,j,k+1);
            localIndex(6) = indexMap(i+1,j,k+1);
            localIndex(7) = indexMap(i+1,j+1,k+1);
            localIndex(8) = indexMap(i,j+1,k+1);
            elem(idx:idx+5,:) = localIndex([1 2 3 7; 1 4 3 7; 1 5 6 7;...
                                            1 5 8 7; 1 2 6 7; 1 4 8 7]);
            idx = idx + 6;
        end
    end
end
% Set this as an initial grid
N0 = size(node,1);
HB = zeros(N0,4);
HB(1:N0,1:3) = repmat((1:N0)',1,3);

repo: ifem


text:
function depth(fileName, sensorType, varName)

    %% Check if input is valid
    if ~isstring(fileName) && ~ischar(fileName)
        error("File name must a string or char array\n")
    end

    if ~strcmp(sensorType, "depth")
       error("Wrong sensor selected\n")
    end

    if ~isvarname(varName)
        error("Invalid variable name\n")
    end

    %% Create a python file
    file = fopen(strcat(fileName, '.py'), 'w');

    % Automatically genrates a python file containing the sensor call back
    % bindings
    fprintf(file, 'import numpy as np\n');
    fprintf(file, 'import cv2\n');
    fprintf(file, '\n');
    fprintf(file, 'def bindSensor(sensor):\n');
    fprintf(file, '    sensor.listen(lambda _image: do_something(_image))\n');
    fprintf(file, '\n');
    fprintf(file, 'def do_something(_image):\n');
    fprintf(file, '    global %s\n', varName);
    fprintf(file, '    data = np.frombuffer(_image.raw_data, dtype=np.dtype("uint8"))\n');
    fprintf(file, '\n');
    fprintf(file, '    # Convert to depth in meters\n');
    fprintf(file, '    data = np.reshape(data, (_image.height, _image.width, 4))\n');
    fprintf(file, '    B = data[:, :, 0].astype("float64")\n');
    fprintf(file, '    G = data[:, :, 1].astype("float64")\n');
    fprintf(file, '    R = data[:, :, 2].astype("float64")\n');
    fprintf(file, '\n');
    fprintf(file, '    data = 1000 * (R + G * 256 + B * 256 * 256) / (256 * 256 * 256 - 1)\n');
    fprintf(file, '\n');

    fprintf(file, '    # Convert the data into MATLAB cast compatible type\n');
    fprintf(file, '    %s = np.ascontiguousarray(data)\n', varName);

    fclose(file);

end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
%-------------------------------------------------------------------------------------------------------------
% This is an implementation of the TWSC algorithm for real-world image denoising
%
% Author:  Jun Xu, csjunxu@comp.polyu.edu.hk / nankaimathxujun@gmail.com
%          The Hong Kong Polytechnic University
%
% Please refer to the following paper if you find this code helps:
%
% @article{TWSC_ECCV2018,
% 	author = {Jun Xu and Lei Zhang and David Zhang},
% 	title = {A Trilateral Weighted Sparse Coding Scheme for Real-World Image Denoising},
% 	journal = {ECCV},
% 	year = {2018}
% }
%
% Please see the file License.txt for the license governing this code.
%-------------------------------------------------------------------------------------------------------------
clear;
GT_Original_image_dir = 'cc/';
GT_fpath = fullfile(GT_Original_image_dir, '*mean.png');
TT_Original_image_dir = 'cc/';
TT_fpath = fullfile(TT_Original_image_dir, '*real.png');
% GT_Original_image_dir = 'C:\Users\csjunxu\Desktop\CVPR2018 Denoising\PolyU\';
% GT_fpath = fullfile(GT_Original_image_dir, '*mean.JPG');
% TT_Original_image_dir = 'C:\Users\csjunxu\Desktop\CVPR2018 Denoising\PolyU\';
% TT_fpath = fullfile(TT_Original_image_dir, '*real.JPG');
GT_im_dir  = dir(GT_fpath);
TT_im_dir  = dir(TT_fpath);
im_num = length(TT_im_dir);

method = 'TWSC';
dataset = 'cc';
write_MAT_dir = [dataset '_Results/'];
write_sRGB_dir = [write_MAT_dir method];
if ~isdir(write_sRGB_dir)
    mkdir(write_sRGB_dir)
end

% Parameters
Par.ps = 6;        % patch size
Par.step = 3;       % the step of two neighbor patches
Par.win = 20;   % size of window around the patch
Par.Outerloop = 7;
Par.Innerloop = 2;
Par.maxIter = 10;
Par.maxrho = 100;
Par.nlspini = 70;
Par.display = 1;
Par.delta = 0;
Par.nlspgap = 0;
Par.lambda1 = 0;
Par.lambda2 = 4.9;

Par.nlspini = 70;
% set Parameters
% record all the results in each iteration
Par.PSNR = zeros(Par.Outerloop, im_num, 'double');
Par.SSIM = zeros(Par.Outerloop, im_num, 'double');
for i = 1 : im_num
    Par.nlsp = Par.nlspini;  % number of non-local patches
    Par.image = i;
    Par.nim = im2double(imread(fullfile(TT_Original_image_dir, TT_im_dir(i).name) ));
    Par.I = im2double(imread(fullfile(GT_Original_image_dir, GT_im_dir(i).name)));
    S = regexp(TT_im_dir(i).name, '\.', 'split');
    IMname = S{1};
    [h,w,ch] = size(Par.nim);
    % noise estimation
    for c = 1:ch
        Par.nSig(c) = NoiseEstimation(Par.nim(:, :, c)*255, Par.ps)/255;
    end
    % initial PSNR and SSIM
    fprintf('%s: \n', TT_im_dir(i).name);
    fprintf('The initial PSNR = %2.4f, SSIM = %2.4f. \n', csnr( Par.nim*255, Par.I*255, 0, 0 ), cal_ssim( Par.nim*255, Par.I*255, 0, 0 ));
    % denoising
    [IMout, Par]  =  TWSC_Sigma_RW(Par);
    % calculate the PSNR
    Par.PSNR(Par.Outerloop, Par.image)  =   csnr( IMout*255, Par.I*255, 0, 0 );
    Par.SSIM(Par.Outerloop, Par.image)      =  cal_ssim( IMout*255, Par.I*255, 0, 0 );
    %% output
    imwrite(IMout, [write_sRGB_dir '/' method '_' dataset '_' IMname '.png']);
    fprintf('%s : PSNR = %2.4f, SSIM = %2.4f \n',TT_im_dir(i).name, Par.PSNR(Par.Outerloop, Par.image),Par.SSIM(Par.Outerloop, Par.image)     );
end
PSNR =Par.PSNR(end,:);
SSIM = Par.SSIM(end,:);
mPSNR=mean(PSNR,2);
mSSIM=mean(SSIM,2);
matname = sprintf([write_MAT_dir method '_' dataset '.mat']);
save(matname,'PSNR','SSIM','mPSNR','mSSIM');

author: csjunxu


text:
function R = QuatToRot(q)
%QuatToRot Converts a Quaternion to Rotation matrix
%   written by Daniel Mellinger

% normalize q
q = q./sqrt(sum(q.^2));

qahat(1,2) = -q(4);
qahat(1,3) = q(3);
qahat(2,3) = -q(2);
qahat(2,1) = q(4);
qahat(3,1) = -q(3);
qahat(3,2) = q(2);

R = eye(3) + 2*qahat*qahat + 2*q(1)*qahat;

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
% The COBRAToolbox: testDeleteProtons.m
%
% Purpose:
%     - tests the deleteProtons function using the DAS network
%
% Authors:
%     - German Preciat -- August 2017
%

% save the current path
currentDir = pwd;

% initialize the test
cd(fileparts(which('testDeleteProtons')));

% Load reference data
load('refData_deleteProtons.mat')

% load the dopamine synthesis model
modelDir = getDistributedModelFolder('subDas.mat');
model = load([modelDir filesep 'subDas.mat']);
if isfield(model,'model')
    model=model.model;
end
%model = readCbModel('subDas.mat');

modelNew = deleteProtons(model);
assert(all(all(modelNew.S == modelNew0.S)), 'Reference S matrix does not match.')
assert(all(all(modelNew.S == modelNew0.S)), 'Reference metabolites do not match.')
assert(all(all(modelNew.S == modelNew0.S)), 'Reference metabolites formulas do not match.')
assert(all(all(modelNew.S == modelNew0.S)), 'Reference metabolites charges do not match.')

% change the directory
cd(currentDir)


text:
function geometry_test0205 ( )

%*****************************************************************************80
%
%% TEST0205 tests DEGREES_TO_RADIANS, RADIANS_TO_DEGREES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0205\n' );
  fprintf ( 1, '  DEGREES_TO_RADIANS converts an angle from degrees\n' );
  fprintf ( 1, '    to radians;\n' );
  fprintf ( 1, '  RADIANS_TO_DEGREES converts an angle from radians\n' );
  fprintf ( 1, '    to degrees;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Degrees     Radians     Degrees\n' );
  fprintf ( 1, '\n' );

  for i = -2 : 14
    angle_deg = 30 * i;
    angle_rad = degrees_to_radians ( angle_deg );
    angle_deg2 = radians_to_degrees ( angle_rad );
    fprintf ( 1, '  %10f  %10f  %10f\n', angle_deg, angle_rad, angle_deg2 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [node,face,elem]=meshabox(p0,p1,opt,nodesize)
%
% [node,face,elem]=meshabox(p0,p1,opt,maxvol)
%
% create the surface and tetrahedral mesh of a box geometry
%
% author: Qianqian Fang, <fangq at nmr.mgh.harvard.edu>
%
% input:
%   p0:  coordinates (x,y,z) for one end of the box diagnoal
%   p1:  coordinates (x,y,z) for the other end of the box diagnoal
%   opt: maximum volume of the tetrahedral elements
%   nodesize: 1 or a 8x1 array, size of the element near each vertex
%
% output:
%   node: node coordinates, 3 columns for x, y and z respectively
%   face: integer array with dimensions of NB x 3, each row represents
%         a surface mesh face element
%   elem: integer array with dimensions of NE x 4, each row represents
%         a tetrahedron
%
% example:
%   [node,face,elem]=meshabox([2 3 2],[6 12 15],0.1,1);
%   plotmesh(node,elem,'x>4');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<4)
   nodesize=1;
end
[node,elem,face]=surf2mesh([],[],p0,p1,1,opt,[],[],nodesize);
elem=elem(:,1:4);
face=face(:,1:3);

repo: REST


text:
function value = problem4 ( x, h, epsilon, o )

%*****************************************************************************80
%
%% PROBLEM4 evaluates problem data for problem 4.
%
%  Discussion:
%
%    case # 4: u_anl  =   x-1/4.......................................x<0.5
%                         x-1/2.......................................x>0.5
%
%              w_anl  =   0...................................x<0.5-epsilon
%                        -2/epsilon^2(-1/4(log(epsilon)-
%                         log(1/2-v)))....................0.5-epsilon<x<0.5
%                         2/epsilon^2( 1/4*(log(v-1/2)-
%                         log(epsilon)))..................0.5<x<0.5+epsilon
%                         0...................................x>0.5+epsilon
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2012
%
%  Parameters:
%
%    Input, real X(:), the evaluation points.
%
%    Input, real H, the discretization parameter.
%
%    Input, real EPSILON, the nonlocalization parameter.
%
%    Input, string O, the option:
%    'U', evaluate the exact solution.
%    'D', evaluate the derivative of the exact solution.
%    'F', evaluate the right hand side.
%    'L', evaluate the lifting function.
%
  if ( o == 'u' )
    value = (x-1/4).*(x<0.5) + (x-1/2).*(x>=0.5);
  elseif ( o == 'd' )
    value = zeros ( size ( x ) );
  elseif ( o == 'f' )
    epsilon2 = epsilon * epsilon;
    value = (x<0.5-epsilon).*0 + ...
            (x>=0.5-epsilon & x<0.5).*(-2)/epsilon2.*(-1/4* (log(epsilon)-log(1/2-x))) + ...
            (x>0.5 & x<0.5+epsilon).*2/epsilon2.*(1/4*(log(x-1/2)-log(epsilon))) + ...
            (x>=0.5+epsilon).*0;
  elseif ( o == 'l' )
    value = (x<=0).*(x-1/4) + ...
            (x>0 & x<=h).* (1/4).*(x./h-1)+ ...
            (x>h & x< 1-h).* 0 * epsilon + ...
            (x>=1-h & x<=1).* (1/2).*(1. + (x-1)./h)  + ...
            (x>1).*(x-1/2);
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PROBLEM4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized option O = "%s"\n', o );
    error ( 'PROBLEM4 - Fatal error!\n' );
  end

  return
end


text:
function [w,msg,msgobj] = my_gencoswin(varargin)
%GENCOSWIN   Returns one of the generalized cosine windows.
%   GENCOSWIN returns the generalized cosine window specified by the
%   first string argument. Its inputs can be
%     Window name    - a string, any of 'hamming', 'hann', 'blackman'.
%     N              - length of the window desired.
%     Sampling flag  - optional string, one of 'symmetric', 'periodic'.

%   Copyright 1988-2010 The MathWorks, Inc.
%   $Revision: 1.7.4.3 $  $Date: 2011/05/13 18:14:10 $

% Parse the inputs
window = varargin{1};
n = varargin{2};
msg = '';
msgobj = [];

% Check for trivial orders
% [n,w,trivialwin] = check_order(n);
% if trivialwin, return, end;

% Select the sampling option
if nargin == 2, % no sampling flag specified, use default.
    sflag = 'symmetric';
else
    sflag = lower(varargin{3});
end

% Allow partial strings for sampling options
allsflags = {'symmetric','periodic'};
sflagindex = strmatch(sflag, allsflags);
if length(sflagindex)~=1         % catch 0 or 2 matches
    msgobj = message('signal:gencoswin:BadFlag','symmetric','periodic');
    msg = getString(msgobj);
    return;
else
    sflag = allsflags{sflagindex};
end

% Evaluate the window
switch sflag
case 'periodic'
    w = sym_window(n+1,window);
    w(end) = [];
case 'symmetric'
    w = sym_window(n,window);
end

%---------------------------------------------------------------------
function w = sym_window(n,window)
%SYM_WINDOW   Symmetric generalized cosine window.
%   SYM_WINDOW Returns an exactly symmetric N point generalized cosine
%   window by evaluating the first half and then flipping the same samples
%   over the other half.

if ~rem(n,2)
    % Even length window
    half = n/2;
    w = calc_window(half,n,window);
    w = [w; w(end:-1:1)];
else
    % Odd length window
    half = (n+1)/2;
    w = calc_window(half,n,window);
    w = [w; w(end-1:-1:1)];
end

%---------------------------------------------------------------------
function w = calc_window(m,n,window)
%CALC_WINDOW   Calculate the generalized cosine window samples.
%   CALC_WINDOW Calculates and returns the first M points of an N point
%   generalized cosine window determined by the 'window' string.

x = (0:m-1)'/(n-1);

switch window
case 'hann'
    % Hann window
    % w = 0.5 * (1 - cos(2*pi*(0:m-1)'/(n-1)));
    w = 0.5 - 0.5*cos(2*pi*x);
case 'hamming'
    % Hamming window
    % w = (54 - 46*cos(2*pi*(0:m-1)'/(n-1)))/100;
    w = 0.54 - 0.46*cos(2*pi*x);
case 'blackman'
    % Blackman window
    % Force end points to zero to avoid close-to-zero negative values caused
    % by roundoff errors.
    % w = (42 - 50*cos(2*pi*(0:m-1)/(n-1)) + 8*cos(4*pi*(0:m-1)/(n-1)))'/100;
    w = 0.42 - 0.5*cos(2*pi*x) + 0.08*cos(4*pi*x);
    w(1) = 0;
case 'flattopwin'
    % Flattop window
    % Coefficients as defined in the reference [1] (see flattopwin.m)
    a0 = 0.21557895;
    a1 = 0.41663158;
    a2 = 0.277263158;
    a3 = 0.083578947;
    a4 = 0.006947368;
    w = a0 - a1*cos(2*pi*x) + a2*cos(4*pi*x) - a3*cos(6*pi*x) + ...
      a4*cos(8*pi*x);
end

% [EOF] gencoswin.m

author: singaxiong

repo: SignalGraph


text:
function [mu, b, se, LR_chi2, convState] = drxlr_logistic_regression(x,y,iterlim,convcrit)
% regression by logistic model
% discrete reponse and different types of explanatory variables
% x: column vector of explanatory variables
% y: column vector of observations
% solved by maximum likelihood estimates using iterative  weighted
% least-squared.
% written by Issam El Naqa Spring 2003
% Extracted for generalized use 2005, AJH
% LM: APA 07/13/2006, added convState as output parameter to represent the
% state of convergence. convState is a 1x3 vector with 1st element being
% converged Yes or no, second element being number of iterations for
% convergence and third element being convergence criteria condition
%
% Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
% This file is part of the Dose Response Explorer System (DREES).
%
% DREES development has been led by:  Issam El Naqa, Aditya Apte, Gita Suneja, and Joseph O. Deasy.
%
% DREES has been financially supported by the US National Institutes of Health under multiple grants.
%
% DREES is distributed under the terms of the Lesser GNU Public License.
%
%     This version of DREES 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.
%
% DREES 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 DREES.  If not, see <http://www.gnu.org/licenses/>.

% Initialize parameters
%warning('off');
x = [x ones(size(x,1),1)]; % constant factor
iter = 0;
[n p] = size(x);
seps = sqrt(eps);
eta = drxlr_rlogit((y + 0.5)/2); % initialize
b = zeros(p,1);
b0 = b+1;
% Start iterating
while(1)
    iter = iter+1;
    mu = drxlr_invlogit(eta);
    mu = max(0, min(1, mu));
    % Check stopping conditions
    if (~any(abs(b-b0) > convcrit * max(seps, abs(b0))))
        convState = [1 iter-1 max(abs(b-b0))];
        break;
    end
    if (iter>iterlim)
        warning('drxlr_logistic_regression: Number of iterations exceeded without convergence (convergence tolerance = %e', convcrit)
        convState = [2 iter-1 max(abs(b-b0))];
        break;  % iteration limit reached, exit!
    end
    deta = drxlr_d_logit(mu);
    z = eta + (y - mu) .* deta;
    vy=drxlr_binomial_variance(mu);
    w = 1 ./ max(eps, (deta .^ 2) .* vy);
    b0 = b;
    [b,R,se] = drxlr_wfit(z, x, w, p);
    eta = x * b;
end

%%% calculate likelihood ratio test
LR_chi2=drxlr_get_lrt(y,mu);
return


text:
function varargout = cutMeshByPlane(v, f, plane, varargin)
%CUTMESHBYPLANE Cut a mesh by a plane.
%
%   [ABOVE, IN, BELOW] = cutMeshByPlane(MESH, PLANE)
%   where MESH, ABOVE, IN, BELOW are structs with the fields vertices and
%   faces, and PLANE is given as a row containing initial point and 2
%   direction vectors. ABOVE, IN, BELOW contain the corresponding parts of
%   the input mesh.
%
%   [ABOVE_V, ABOVE_F, IN_V, IN_F, BELOW_V,BELOW_F] = ...
%       cutMeshByPlane(V, F, PLANE) where V is a [NVx3] array containing
%   coordinates and F is a [NFx3] array containing indices of vertices of
%   the triangular faces.
%
%   BELOW = cutMeshByPlane(V, F, PLANE, 'part', 'below') BELOW is a struct
%   with the fields vertices and faces. Other options are:
%       'part'  -   'above': Faces above the plane
%               -   'in'   : Faces in the plane
%               -   'below': Faces below the plane
%
%   [BELOW_V, BELOW_F] = cutMeshByPlane(MESH, PLANE, 'part', 'below') is
%   possible, too.
%

% ---------
% Authors: oqilipo, David Legland
% Created: 2017-07-09
% Copyright 2017

narginchk(2,5)
nargoutchk(1,6)

%% Parse inputs
% If first argument is a struct
if nargin == 2 || nargin == 4
    if ~isempty(varargin)
        varargin = [{plane}, varargin(:)'];
    end
    plane = f;
    [v, f] = parseMeshData(v);
end

p = inputParser;
addRequired(p,'plane',@isPlane)
validStrings = {'above','in','below'};
addParameter(p,'part','above',@(x) any(validatestring(x, validStrings)))
parse(p, plane, varargin{:});
part=p.Results.part;

%% Algorithm
% Logical index to the vertices below the plane
VBPl_LI = isBelowPlane(v, plane);

% Logical index to three vertices of each face
FBP_LI = VBPl_LI(f);
switch nargout
    case {1, 2}
        switch part
            case 'above'
                % Faces above the plane, all three vertices == 0 -> sum has to be 0
                above = removeMeshFaces(v, f, ~(sum(FBP_LI, 2) == 0) );
            case 'in'
                % Faces in the plane, 1 or 2 vertices == 0 -> sum can be 1 or 2
                inside = removeMeshFaces(v, f, ~((sum(FBP_LI, 2) > 0 & sum(FBP_LI, 2) < 3)));
            case 'below'
                % Faces below the plane, all three vertices == 1 -> sum has to be 3
                below = removeMeshFaces(v, f, ~(sum(FBP_LI, 2) == 3) );
        end
    case {3, 6}
        % Faces above the plane, all three vertices == 0 -> sum has to be 0
        above = removeMeshFaces(v, f, ~(sum(FBP_LI, 2) == 0) );
        % Faces in the plane, 1 or 2 vertices == 0 -> sum can be 1 or 2
        inside = removeMeshFaces(v, f, ~((sum(FBP_LI, 2) > 0 & sum(FBP_LI, 2) < 3)));
        % Faces below the plane, all three vertices == 1 -> sum has to be 3
        below = removeMeshFaces(v, f, ~(sum(FBP_LI, 2) == 3) );
    otherwise
        error('Invalid number of output arguments')
end

%% Parse outputs
switch nargout
    case 1
        switch part
            case 'above'
                varargout{1}=above;
            case 'in'
                varargout{1}=inside;
            case 'below'
                varargout{1}=below;
        end
    case 2
        switch part
            case 'above'
                varargout{1}=above.vertices;
                varargout{2}=above.faces;
            case 'in'
                varargout{1}=inside.vertices;
                varargout{2}=inside.faces;
            case 'below'
                varargout{1}=below.vertices;
                varargout{2}=below.faces;
        end
    case 3
        varargout{1}=above;
        varargout{2}=inside;
        varargout{3}=below;
    case 6
        varargout{1}=above.vertices;
        varargout{2}=above.faces;
        varargout{3}=inside.vertices;
        varargout{4}=inside.faces;
        varargout{5}=below.vertices;
        varargout{6}=below.faces;
    otherwise
        error('Invalid number of output arguments')
end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function render3DTraces(mtx,colors);
% render3DTraces(mtx,[colors]);
%
% Using plot3, render the columns in a 2D matrix mtx
% as a set of traces in a 3D graph (in current axes).
%
% colors: optional cell specifying the color to use for
% each column (string -- see help plot -- or 3-vector).
%
% 07/04 ras.
if ~exist('colors','var')
    colOrder = 'rbgcky';
    for i = 1:size(mtx,2)
        colors{i} = colOrder(mod(i,length(colOrder))+1);
    end
end

cla;
hold on;

X = 1:size(mtx,1);

for i = 1:size(mtx,2)
    Y = i*ones(1,size(mtx,1));
    htmp = plot3(X,Y,mtx(:,i)');
    set(htmp,'LineWidth',1.5,'Color',colors{i});
end

axis auto;
view(-8,60);
grid on;
hold off;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft