text:
function PsiDraw = samplePsi_DK(YData,YsData,B,H,V,G,CD,meanTS,kappa,dataValues,KFSmatrices)
% This function takes a draw from the conditional posterior distribution of
% Psi (local mean).
% B=Bdraw;
% H=Hdraw;
% V=Vdraw;
% G=Gdraw;
% CD=CDdraw;
% kappa=priorValues.kappa
%% Initialize
yData=YData;
zData=YsData;

[T,n] = size(yData);
Mz = size(zData,2); %number of survey forecasts
p=size(B,1)/n;     %lags
na=n*(p+1);        %number of variables in state space

%% Matrices of state space representation
Zmatrix=KFSmatrices.Zmatrix; %(matrix for measurement equation)
Tmatrix=KFSmatrices.Tmatrix; %(matrix for state equation)
Tmatrix(n+1:2*n,n+1:end) = B'; %put the VAR coefficients in
Rmatrix=KFSmatrices.R;        %scaling in front of errors for transition equation
Lmatrix=KFSmatrices.L;        %scaling in front of errors for measurement equation

%% Data vector and KF initialisation
ZYcol = [zData';yData']; %concatenate surveydata and var data

psiInit    = mean(ZYcol(Mz+1:end,1:p),2); %rowwhise mean
y_PsiStart = reshape(ZYcol(Mz+1:end,p:-1:1),n*p,1)-repmat(psiInit,p,1); %initial values in deviation from mean
aStart=[meanTS;  y_PsiStart]; %start values for transition equation (RHS) with mean followed by lags of yt in deviation from mean

Pstart=kappa*eye(na); %variance for the initial state

aplus = nan(na,T-p); %matrix storing all the period specific for transition equation
QQ = nan(na,na,T-p);
LL = nan(Mz+n,Mz+n,T-p);
yplus = nan(n+Mz,T-p);
aplus(:,1) = aStart+mvnrnd(zeros(na,1),Pstart,1)'; %initialiye RHS of transition equation equation
yplus(:,1) = Zmatrix*aplus(:,1)+[sqrt(G(p+1,:)').*randn(Mz,1);zeros(n,1)]; %initialize the RHS of measurement equation
QQ(:,:,1) = Rmatrix*blkdiag(diag(V(p+1,:)),H(:,:,p+1))*Rmatrix'; %initial VCV of transition equation
LL(:,:,1) = Lmatrix*diag(G(p+1,:))*Lmatrix';                     %initial VCV of measurement equation

%recursively simulate the matrices by drawing from multivariate normal
%using time period specific VCV where
%H is the VCV of the VAR residuals,
%V is the VCV of the state transition residuals and
%G is the VCV of the measurement equation for the survey data

%mean correction simulation smoother
for t = 2:T-p
    %simulate the aplus matrix (constant in state transition equation) by drawing from the period specific distribution of the transition equation and VAR residuals
    aplus(:,t) =  Tmatrix*aplus(:,t-1)+ [sqrt(V(t+p,:)').*randn(n,1);  mvnrnd(zeros(n,1),H(:,:,t+p),1)';zeros(na-2*n,1)];
    %simulate the yplus matrix (containing survey data and var data) by drawing from the period specific
    %distribution of the residuals from the measurement equation
    yplus(:,t) = Zmatrix*aplus(:,t)+[sqrt(G(t+p,:)').*randn(Mz,1);zeros(n,1)];
    %period specific VCV of transition equation
    QQ(:,:,t) = Rmatrix*blkdiag(diag(V(t+p,:)),H(:,:,t+p))*Rmatrix';
   %period specific VCV of measurement equation
    LL(:,:,t) = Lmatrix*diag(G(t+p,:))*Lmatrix';
end

ystar = ZYcol(:,p+1:T)-yplus; %
ahatstar = bear.runKF_DK(ystar, Tmatrix, Zmatrix, QQ, LL, zeros(na,1), Pstart);
atilda = ahatstar+aplus;
PsiDraw = [nan(p,n);atilda(1:n,:)'];

end

author: european-central-bank

repo: BEAR-toolbox


text:
%% Script: Illustrate use of FanChart function

clear all; clc;

% Generate artifical history and forecasts
T=18; % number of historical observations
H=12; % number of forecast periods
fore=randn(H,100);
for h=1:H
   fore(h,:)=(fore(h,:).*h)./200;
end
hist=.1.*randn(T,1);

% Call fan chart function
FanChart(fore,hist,'Jun2009',6);

author: Sable

repo: mcbench-benchmarks


text:
function [lambda_perc, alpha_rad] = calcWheelSlips(omega_rad, VehicleState, DeltaWheel_rad, tw_front_m, tw_rear_m, l_front_m, l_rear_m, tyreradius_front_m, tyreradius_rear_m, v_min)

% initialize variables
alpha_rad = zeros(4, 1);
v_ref = zeros(4, 1);
v_wheel = zeros(4, 1);

vxCG_mps = VehicleState(1);
vyCG_mps = VehicleState(2);
dPsi_radps = VehicleState(3);

% calculate chassis velocities at wheel centers projected to the ground
vx_fr = vxCG_mps + (dPsi_radps*tw_front_m*0.5);
vx_fl = vxCG_mps - (dPsi_radps*tw_front_m*0.5);
vx_rr = vxCG_mps + (dPsi_radps*tw_rear_m*0.5);
vx_rl = vxCG_mps - (dPsi_radps*tw_rear_m*0.5);
vy_fr = vyCG_mps + (dPsi_radps*l_front_m);
vy_fl = vyCG_mps + (dPsi_radps*l_front_m);
vy_rr = vyCG_mps - (dPsi_radps*l_rear_m);
vy_rl = vyCG_mps - (dPsi_radps*l_rear_m);
% calculate rotation matrix for front tires to transform to tire coordinates
R = [cos(-DeltaWheel_rad), -sin(-DeltaWheel_rad);...
  sin(-DeltaWheel_rad), cos(-DeltaWheel_rad)];
% transform to tire coordinates
vxT_fl = R(1, :)*[vx_fl; vy_fl];
vyT_fl = R(2, :)*[vx_fl; vy_fl];
vxT_fr = R(1, :)*[vx_fr; vy_fr];
vyT_fr = R(2, :)*[vx_fr; vy_fr];
% calculate tire slip angles
if(vxCG_mps > v_min)
  alpha_rad(1) = atan2(-vyT_fl, vxT_fl);
  alpha_rad(2) = atan2(-vyT_fr, vxT_fr);
  alpha_rad(3) = atan2(-vy_rl, vx_rl);
  alpha_rad(4) = atan2(-vy_rr, vx_rr);

    % calculate absolute reference velocities of tire
    v_ref(1) = vxT_fl;
    v_ref(2) = vxT_fr;
    v_ref(3) = vx_rl;
    v_ref(4) = vx_rr;
    % calculate wheel over ground velocities
    v_wheel(1:2) = tyreradius_front_m.*omega_rad(1:2);
    v_wheel(3:4) = tyreradius_rear_m.*omega_rad(3:4);
    % calculate tireslips
    lambda_perc = 100*(v_wheel-v_ref)./(max(v_min, v_ref));
else
  alpha_rad = zeros(4,1);
  lambda_perc = zeros(4, 1);
end

author: TUMFTM


text:
function motionCompPlayMeanMovie(slice,view)
%
%    gb 04/30/05
%
%    motionCompPlayMeanMovie(view,slice)
%
% Plays the movie of the mean maps of the choosen slice, or creates the
% movie if it does not exist.
%
% The argument view is optional but an error in thrown if it is not
% passed in and the movie does not exist
%
if ~ieNotDefined('view')
    cd(dataDir(view))
end

if ieNotDefined('slice')
    slice = 18;
end

if ~exist(['movie_' num2str(slice) '.avi'])
    try
        motionCompMeanMovie(view,slice)
    catch
        error('The movie has not been created. Please insert a view argument');
    end
end

aviobj = aviread(['movie_' num2str(slice)]);
mplay(aviobj);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef BOWImgDescriptorExtractor < handle
    %BOWIMGDESCRIPTOREXTRACTOR  Class to compute an image descriptor using the bag of visual words
    %
    % Such a computation consists of the following steps:
    %
    % 1. Compute descriptors for a given image and its keypoints set.
    % 2. Find the nearest visual words from the vocabulary for each keypoint
    %    descriptor.
    % 3. Compute the bag-of-words image descriptor as is a normalized
    %    histogram of vocabulary words encountered in the image. The i-th bin
    %    of the histogram is a frequency of i-th word of the vocabulary in
    %    the given image.
    %
    % ## Example
    %
    %     % create bag of visual words
    %     trainer = cv.BOWKMeansTrainer(K);
    %     dictionary = trainer.cluster(train_descs);
    %
    %     % Compute histogram of visual word occurrences of an image
    %     extractor = cv.BOWImgDescriptorExtractor('SIFT','BruteForce');
    %     extractor.Vocabulary = dictionary;
    %     descs = extractor.compute(im, keypoints);
    %
    % See also: cv.BOWImgDescriptorExtractor.BOWImgDescriptorExtractor,
    %  cv.BOWKMeansTrainer, bagOfFeatures, trainImageCategoryClassifier,
    %  imageCategoryClassifier, indexImages, retrieveImages
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Visual vocabulary
        %
        % Vocabulary (can be trained using cv.BOWKMeansTrainer). Each row of
        % the vocabulary is a visual word (cluster center).
        Vocabulary
    end

    methods
        function this = BOWImgDescriptorExtractor(dextractor, dmatcher)
            %BOWIMGDESCRIPTOREXTRACTOR  The constructor
            %
            %     extractor = cv.BOWImgDescriptorExtractor(dextractor, dmatcher)
            %     extractor = cv.BOWImgDescriptorExtractor({dextractor, 'key',val,...}, {dmatcher, 'key',val,...})
            %
            % ## Input
            % * __dextractor__ Descriptor extractor that is used to compute
            %   descriptors for an input image and its keypoints. It can be
            %   specified by a string containing the type of descriptor
            %   extractor, such as 'SIFT' or 'SURF'. See
            %   cv.DescriptorExtractor.DescriptorExtractor for possible types.
            % * __dmatcher__ Descriptor matcher that is used to find the
            %   nearest word of the trained vocabulary for each keypoint
            %   descriptor of the image. It can be spacified by a string
            %   specifying the type of descriptor extractor, such as
            %   'BruteForce' or 'BruteForce-L1'. See
            %   cv.DescriptorMatcher.DescriptorMatcher for possible types.
            %   default 'BruteForce'
            %
            % In the first variant, it creates descriptor extractor/matcher of
            % the given types using default parameters (by calling the default
            % constructors).
            %
            % In the second variant, it creates descriptor extractor/matcher
            % of the given types using the specified options. Each algorithm
            % type takes optional arguments. Each of the extractor/matcher are
            % specified by a cell-array that starts with the type name
            % followed by option arguments, as in:
            % `{'Type', 'OptionName',optionValue, ...}`. Refer to the
            % individual extractor/matcher functions to see a list of possible
            % options of each algorithm.
            %
            % ## Example
            %
            %     % first variant
            %     extractor = cv.BOWImgDescriptorExtractor('ORB', 'BruteForce');
            %
            %     % second variant
            %     extractor = cv.BOWImgDescriptorExtractor(...
            %         {'FastFeatureDetector', 'Threshold',10}, ...
            %         {'BFMatcher', 'NormType','L2'});
            %
            % See also: cv.DescriptorExtractor, cv.DescriptorMatcher
            %
            if nargin < 2, dmatcher = 'BruteForce'; end
            this.id = BOWImgDescriptorExtractor_(0, 'new', dextractor, dmatcher);
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     extractor.delete()
            %
            % See also: cv.BOWImgDescriptorExtractor
            %
            if isempty(this.id), return; end
            BOWImgDescriptorExtractor_(this.id, 'delete');
        end

        function sz = descriptorSize(this)
            %DESCRIPTORSIZE  Returns image discriptor size
            %
            %     sz = extractor.descriptorSize()
            %
            % ## Output
            % * __sz__ Returns an image discriptor size if the vocabulary is
            %   set. Otherwise, it returns 0.
            %
            % This is basically `size(Vocabulary,1)` (i.e number of clusters).
            %
            % See also: cv.BOWImgDescriptorExtractor.descriptorType
            %
            sz = BOWImgDescriptorExtractor_(this.id, 'descriptorSize');
        end

        function dtype = descriptorType(this)
            %DESCRIPTORTYPE  Returns image descriptor type
            %
            %     dtype = extractor.descriptorType()
            %
            % ## Output
            % * __dtype__ Returns an image descriptor type, one of numeric
            %   MATLAB class names.
            %
            % Always `single` for BOWImgDescriptorExtractor.
            %
            % See also: cv.BOWImgDescriptorExtractor.descriptorSize,
            %  cv.DescriptorExtractor.descriptorType
            %
            dtype = BOWImgDescriptorExtractor_(this.id, 'descriptorType');
        end

        function [bow,idx,kptDescs] = compute(this, img, keypoints)
            %COMPUTE  Computes an image descriptor using the set visual vocabulary
            %
            %     [bow,idx,kptDescs] = extractor.compute(img, keypoints)
            %
            % ## Input
            % * __img__ Image, for which the descriptor is computed.
            % * __keypoints__ Keypoints detected in the input image. It is a
            %   struct array that is returned by cv.FeatureDetector.detect.
            %
            % ## Output
            % * __bow__ Computed output image descriptor. A vector of the same
            %   length as the vocabulary dimension.
            % * __idx__ Indices of keypoints that belong to the cluster. A
            %   cell array of integer vectors. This means that `idx{i}` are
            %   keypoint indices that belong to the i-th cluster (word of
            %   vocabulary).
            % * __kptDescs__ Descriptors of the image keypoints, as returned
            %   by cv.DescriptorExtractor.compute.
            %
            % See also: cv.BOWImgDescriptorExtractor.compute2,
            %  cv.BOWImgDescriptorExtractor.compute1
            %
            if nargout > 2
                [bow,idx,kptDescs] = BOWImgDescriptorExtractor_(this.id, 'compute', img, keypoints);
            elseif nargout > 1
                [bow,idx] = BOWImgDescriptorExtractor_(this.id, 'compute', img, keypoints);
            else
                bow = BOWImgDescriptorExtractor_(this.id, 'compute', img, keypoints);
            end
        end

        function [bow,idx] = compute1(this, kptDescs)
            %COMPUTE1  Computes an image descriptor using keypoint descriptors
            %
            %     [bow,idx] = extractor.compute1(kptDescs)
            %
            % ## Input
            % * __kptDescs__ Computed descriptors to match with vocabulary. It
            %   is a numeric matrix that is returned by
            %   cv.DescriptorExtractor.compute.
            %
            % ## Output
            % * __bow__ Computed output image descriptor. A vector of the same
            %   length as the vocabulary dimension.
            % * __idx__ Indices of keypoints that belong to the cluster. A
            %   cell array of integer vectors. This means that `idx{i}` are
            %   keypoint indices that belong to the i-th cluster (word of
            %   vocabulary).
            %
            % See also: cv.BOWImgDescriptorExtractor.compute
            %
            if nargout > 1
                [bow,idx] = BOWImgDescriptorExtractor_(this.id, 'compute1', kptDescs);
            else
                bow = BOWImgDescriptorExtractor_(this.id, 'compute1', kptDescs);
            end
        end

        function bow = compute2(this, img, keypoints)
            %COMPUTE2  Computes an image descriptor using the set visual vocabulary
            %
            %     bow = extrctor.compute2(img, keypoints)
            %
            % ## Input
            % * __img__ Image, for which the descriptor is computed.
            % * __keypoints__ Keypoints detected in the input image. It is a
            %   struct array that is returned by cv.FeatureDetector.detect.
            %
            % ## Output
            % * __bow__ Computed output image descriptor. A vector of the same
            %   length as the vocabulary dimension.
            %
            % See also: cv.BOWImgDescriptorExtractor.compute
            %
            bow = BOWImgDescriptorExtractor_(this.id, 'compute2', img, keypoints);
        end
    end

    %% Getters/Setters
    methods
        function value = get.Vocabulary(this)
            value = BOWImgDescriptorExtractor_(this.id, 'get', 'Vocabulary');
        end
        function set.Vocabulary(this, value)
            BOWImgDescriptorExtractor_(this.id, 'set', 'Vocabulary', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


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

%%
clear; close all; clc;

%% Syntax
% |[meshStruct]=tetMeshBox(boxDim,pointSpacing);|

%% Description
% This function generates a mesh structure containing element and node data
% for a tetrahedral element meshed box. The box dimensions in each of the 3
% directions are based on the boxDim input (1x3 vector). The number of
% elements in each direction is based on the pointSpacing input.

%% Examples

%%
% Plot settings

fontSize=20;
faceAlpha1=0.8;

%% CREATING A MESHED BOX
boxDim=[5 6 7]; % Box dimenstions
pointSpacing=1;

[meshStruct]=tetMeshBox(boxDim,pointSpacing);

%%
% Acces output fields
E=meshStruct.elements;
V=meshStruct.nodes;
F=meshStruct.faces;
Fb=meshStruct.facesBoundary;
faceBoundaryMarker=meshStruct.boundaryMarker;

%%
% Plotting model
cFigure;
title('Box boundaries faces','FontSize',fontSize);
hold on;

gpatch(Fb,V,faceBoundaryMarker,'k',faceAlpha1);
% patchNormPlot(Fb,V);

axisGeom(gca,fontSize);
colormap(gjet(6)); icolorbar;
drawnow;

%%
% Visualizing model internal mesh with |meshView|

meshView(meshStruct);

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function triangulation_test21 ( )

%*****************************************************************************80
%
%% TEST21 tests TRIANGULATION_ORDER3_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  node_num = 9;
  triangle_num = 12;

  node_xy(1:2,1:node_num) = [ ...
       0.0, 0.0; ...
       0.0, 1.0; ...
       0.2, 0.5; ...
       0.3, 0.6; ...
       0.4, 0.5; ...
       0.6, 0.4; ...
       0.6, 0.5; ...
       1.0, 0.0; ...
       1.0, 1.0 ]';
  triangle_node(1:3,1:triangle_num) = [ ...
       2, 1, 3; ...
       3, 1, 6; ...
       2, 3, 4; ...
       4, 3, 5; ...
       7, 4, 5; ...
       5, 3, 6; ...
       7, 5, 6; ...
       9, 4, 7; ...
       6, 1, 8; ...
       7, 6, 8; ...
       7, 8, 9; ...
       2, 4, 9 ]';
  triangle_neighbor(1:3,1:triangle_num) = [ ...
       -28,   2,  3; ...
         1,   9,  6; ...
         1,   4, 12; ...
         3,   6,  5; ...
         8,   4,  7; ...
         4,   2,  7; ...
         5,   6, 10; ...
        12,   5, 11; ...
         2, -34, 10; ...
         7,   9, 11; ...
        10, -38,  8; ...
         3,   8, -3 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST21\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_PRINT prints out a triangulation.\n' );

  triangulation_order3_print ( node_num, triangle_num, node_xy, ...
    triangle_node, triangle_neighbor );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

function DecodedMarginalsToChars(decodedMarginals)
    chars = 'abcdefghijklmnopqrstuvwxyz';
    fprintf('%c', chars(decodedMarginals));
    fprintf('\n');
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function y = nanmax(x,dim)
% nanmax - Max ignoring NaNs.
%
% Synopsis:
%   y = nanmax(x)
%   y = nanmax(x,dim)
%
% Arguments:
%   x: Matrix or vector
%   dim: Dimension along which sum operates. Default: First non-singleton
%       dimension.
%
% Returns:
%   y: max along the chosen dimension, treating NaNs as missing values.
%
% Description:
%   For vectors, nanmax(X) is the max of the non-NaN elements in
%   X. For matrices, nanmax(X) is a row vector containing the max
%   of the non-NaN elements in each column of X. For N-D arrays,
%   nanmax(X) operates along the first non-singleton dimension.
%
%   nanmax(X,dim) determaxes the max along the dimension dim.
%
% Examples:
%   nanmax([1 2 NaN]) returns 2.
%   nanmax([1 2 NaN], 1) returns [1 2 NaN].
%   nanmax([1 2 NaN], 2) returns 2.
%
% See also: nansum,nanmean,nanstd
%

% Author(s): Benjamax Blankertz, Anton Schwaighofer, Aug 2005

if nargin<2,
  % Operate along the first non-singleton dimension
  dim = max(find(size(x)~=1));
  if isempty(dim),
    dim = 1;
  end
end
% Replace NaNs with zeros.
nans = isnan(x);
x(nans) = -inf;

% Protect against an entire column of NaNs
y = max(x, [], dim);
allNaNs = all(nans, dim);
y(allNaNs) = NaN;

author: bbci


text:
function [G]=gsp_low_stretch_tree(k)
%GSP_LOW_STRETCH_TREE  Initialize a low stretch tree
%   Usage:  G = gsp_low_stretch_tree(k);
%           G = gsp_low_stretch_tree();
%
%   Input parameters:
%         k     : 2^k points on each side of the grid of vertices. (default 6)
%   Output parameters:
%         G     : Graph structure.
%
%   'gsp_create_low_stretch_tree(k)' initializes a graph structure containing
%   the weighted adjacency matrix (G.W), the number of vertices (G.N), the
%   plotting coordinates (G.coords), the plotting coordinate limits
%   (G.limits), and the root of a low stretch tree on a grid of points.
%   There are $2^k$ points on each side of the grid, and therefore $2^{2k}$
%   total vertices. The edge weights are all equal to 1.
%
%   Example:::
%
%          G = gsp_low_stretch_tree(3);
%          paramplot.show_edges = 1;
%          gsp_plot_graph(G,paramplot);
%
%   See also: gsp_graph

% Author : David I Shuman, Nathanael Perraudin

if nargin < 1
   k = 6;
end

ii=[2 3 1 1 4 3];
jj=[1 1 2 3 3 4];

XCoords=[1, 2, 1, 2];
YCoords=[1, 1, 2, 2];

for p=2:k
    % create the indicies of the weight matrice
    ii_new=[ii,ii+4^(p-1),ii+2*4^(p-1),ii+3*4^(p-1)];
    ii_new_middle=[4^(p-1),4^(p-1),4^(p-1)+(4^p+2)/3,5/3*4^(p-1)+1/3,4^(p-1)+(4^p+2)/3,3*4^(p-1)+1];
    ii=[ii_new,ii_new_middle];

    jj_new=[jj,jj+4^(p-1),jj+2*4^(p-1),jj+3*4^(p-1)];
    jj_new_middle=[5/3*4^(p-1)+1/3,4^(p-1)+(4^p+2)/3,3*4^(p-1)+1,4^(p-1),4^(p-1),4^(p-1)+(4^p+2)/3];
    jj=[jj_new,jj_new_middle];

    % Create Coords
    YCoords=repmat(YCoords,1,2);
    YCoords_new=[YCoords,YCoords+2^(p-1)];
    YCoords=YCoords_new;

    XCoords_new=[XCoords,XCoords+2^(p-1)];
    XCoords=repmat(XCoords_new,1,2);
end

G.W=sparse(ii,jj,ones(size(ii)));
G.coords=[XCoords',YCoords'];

G.limits=[0,2^k+1,0,2^k+1];
G.N=(2^k)^2;
G.root=4^(k-1);

G.plotting.edge_width=1.25;
G.plotting.vertex_size=75;

G.type = 'low strech tree';

G = gsp_graph_default_parameters(G);

end

repo: gspbox


text:
function[x] = IMSVD(Y,U)
% inverse MSVD (IMSVD)
% Inputs-> Y: MSVD coefficients & U: unitary matrix (U in SVD)
% output-> x: image (spaitial domain)

[m,n] = size(Y.LL);
mn = m*n;
T = zeros(4,mn);
T(1,:) = reshape(Y.LL,1,mn);
T(2,:) = reshape(Y.LH,1,mn);
T(3,:) = reshape(Y.HL,1,mn);
T(4,:) = reshape(Y.HH,1,mn);
A = U * T;
x = zeros(m*2,n*2);
for j = 1:n
    for i = 1:m
        x((i-1)*2+(1:2), (j-1)*2+(1:2)) = reshape(A(:,i+(j-1)*m),2,2);
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function X=diag(X,k)
%DIAG (overloaded)

n = X.dim(1);
m = X.dim(2);

if nargin == 1
    k = 0;
end

if min([n m])==1
    Y = X;
	diagX = diag(X.basis(:,1),k);
	Y.basis = diagX(:);
	for i = 1:length(X.lmi_variables)
		diagX = diag(X.basis(:,i+1),k);
		Y.basis(:,i+1) = diagX(:);
    end
	Y.lmi_variables = X.lmi_variables;
	Y.dim(1) = size(diagX,1);
	Y.dim(2) = size(diagX,2);
    X=Y;
else
    index = diag(reshape(1:n*m,n,m),k);
    X.basis = X.basis(index,:);
	X.dim(1) = length(index);
	X.dim(2) = 1;
    X = clean(X);
end
% Reset info about conic terms
if isa(X,'sdpvar')
    X.conicinfo = [0 0];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [ T, trnorm, U ] = tensor_shrinkage( X, tau, mode )
% X and T are tensors
global tmode;
tmode = mode;

Xmat = tenmat(X,mode);
[Tmat, trnorm, U] = matrix_shrinkage( double(Xmat), tau );
% [ U, S, V ] = svd( double(Xmat), 'econ' );
% s = diag(S);
% s = s - tau;
% s( s < 0 ) = 0;
% Tmat = scale_matrix( U, s, 1 ) * V';
T = tensor( tenmat( Tmat, Xmat.rdims, Xmat.cdims, Xmat.tsize ) );

end

author: andrewssobral

repo: lrslibrary


text:
function p = root()
    %ROOT  Obtain mexopencv root path
    %
    persistent p_cache
    if isempty(p_cache)
        p_cache = fileparts(fileparts(mfilename('fullpath')));
    end
    p = p_cache;
end


text:
function [gamma,phase,dta,alpha,gamma_xmesh,gamma_ymesh,gamma_zmesh]=dicomrt_GAMMAcal2DVol(evalm,refm,dose_xmesh,dose_ymesh,dose_zmesh,resf,range,dta_criteria,dd_criteria,voi,voiselect,voilookup)
% dicomrt_GAMMAcal2DVol(evalm,refm,dose_xmesh,dose_ymesh,dose_zmesh,resf,range,dta_criteria,dd_criteria,voi,voiselect,voilookup)
%
% Calculates 2D gamma distribution for a 3D dataset using a 2D algorithm.
% Calculates 2D gamma maps on every slice of a user selected voi and stack 2D maps onto a 3D array.
%
% evalm is the 2D matrix to evaluate.
% refm is the reference 2D matrix.
%   Both eval and ref can be a TPS generated dataset or a MC generated dataset.
%   Doses are normalised to the prescribed target dose whenever is possible.
%   If no normalization dose is provided within the data it is assumed that matrices are already normalised.
% dose_xmesh, dose_ymesh, are the coordinates of the center of the pixel for eval and ref.
% slice is the number of the slice quantities should be calculated for.
% resf is the "resolution factor". Each voxel is divided resf times to allow dose to be interpolated
%   and quantities calculated. The higher resf the more accurate the calculation is, the slower the
%   function will be.
% range is the "search range". Range is the number of pixel about (j,i) that is considered for calculation.
%   If a dta or a dd match is not found within the search range gamma and all the other quantities in that point
%   will be NaN (Not a Number).
% dta_criteria is the Distance-to-agreement criteria in cm (e.g. 0.3).
% dd_criteria:
%   1) IF matrices ARE NOT NORMALIZED before calling the function
%      dta_criteria is the percentage dose difference (e.g. =0.03);
%   2) IF matrices ARE already NORMALIZED before calling the function
%      dd_criteria must be given accordingly to the normalization applyed
%      (e.g. =0.03 for 3% over dose norm =1Gy, or =1.98 for 3% over dose norm=66Gy).
% voi and voiselect are the vois' cell array and the # of the voi to be used for the gamma calculation respectively.
%   They have to be specified together. Both matrices will be masked and reduced in size accordingly with the selected
%   voi's dimensions. This reduce calculation time especially for the 3D algorithm.
%
% NOTE: the use of tyhis function is equivalent to looping calls to dicomrt_GAMMAcal2D, storing
% results in a 3D matrix. Useful and not mandatory parameter that this function returns is zmesh, which
% is needed by dicomrt_gvhcal.
%
% Example:
%
% [gamma1_2,phase1_2,dta1_2,alpha1_2,xmesh,ymesh,zmesh]=dicomrt_GAMMAcal2DVol(image_eval,image_ref, ...
%                                                             xmesh_red,ymesh_red,40,2,4,0.3,0.03);
%
% returns the calculated gamma function for image_eval vs image_ref at slice # 40 in gamma1_2.
% Gamma is calculated with 3%-3mm DD-DTA criteria.
% The phase is returned in phase1_2, the dta matrix in dta1_2.
% The function return also the direction cosines in alpha1_2. These are the
% direction cosines on the gamma vector in the two dimensional space xy. The direction
% cosines may provide useful information in detecting systematic spatial displacements between eval and ref.
% Coordinates of the xyz location where gamma is defined are returned in xmesh, ymesh, and zmesh.
%
% The concept of gamma function was developed by Low et al Med. Phys. 25 656-661.
%
% See also dicomrt_GAMMAcal3DP, dicomrt_loaddose, dicomrt_loadmcdose, dicomrt_gvhcal
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(11,12,nargin))

if exist('voilookup')~=1
    voilookup=1; % assume VOI number 1 is patient outline
end

% Check case and set-up some parameters and variables
[evalm_temp,evalm_type,evalm_label,evalm_PatientPosition]=dicomrt_checkinput(evalm);
evalm=dicomrt_varfilter(evalm_temp);
[refm_temp,refm_type,refm_label,refm_PatientPosition]=dicomrt_checkinput(refm);
refm=dicomrt_varfilter(refm_temp);
[voi_temp]=dicomrt_checkinput(voi);
voi=dicomrt_varfilter(voi_temp);
voitype=dicomrt_checkvoitype(voi_temp);

% Retrieve
%voi_start=dicomrt_findslice(voi_temp,voiselect,1,voilookup);
%voi_stop=voi_start+size(voi{voiselect,2},1)-1;
%bias=voi_start-1;

gamma=[];
phase=[];
dta=[];
alpha=[];
gamma_zmesh=[];

% Mask the matrix using the selected VOI
voiZ=dicomrt_makevertical(dicomrt_getvoiz(voi_temp,voiselect));

for i=1:length(voiZ)
    [locate_slice]=dicomrt_findpointVECT(dose_zmesh,voiZ(i),evalm_PatientPosition);
    disp(['Working on ',voi{voiselect,1}, ' slice: ', num2str(i)]);
    [gamma(:,:,i),phase(:,:,i),dta(:,:,i),alpha(:,:,i),gamma_xmesh,gamma_ymesh,zmesh_temp]=...
        dicomrt_GAMMAcal2D(evalm_temp,refm_temp,dose_xmesh,dose_ymesh,dose_zmesh,...
        locate_slice,resf,range,dta_criteria,dd_criteria,voi_temp,voiselect);
    gamma_zmesh=[gamma_zmesh zmesh_temp];
end

gamma=dicomrt_restorevarformat(evalm_temp,gamma);

% Label info and update time of creation
gamma{1,1}{1}.RTPlanLabel=[gamma{1,1}{1}.RTPlanLabel,'-GAMMA'];
gamma{1,1}{1}.RTPlanDate=date;
time=fix(clock);
creationtime=[num2str(time(4)),':',num2str(time(5))];
gamma{1,1}{1}.RTPlanTime=creationtime;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [predict_label,accuracy] = jdrf(x_train,y_train,x_test,y_test,n_trees)
%%jdrf:packaged TreeBagger from Matlab,just for easier use.
%%==============================================================================
%%input:
%%------x_train,...,y_test  :   training and testing sets.        [required]
%%------n_trees             :   number of trees,default is 10.    [not required]

%%output:
%%------predict_label       :   predicted label vector for test case.
%%------accuracy            :   accuracy
%%==============================================================================
	predict_label = [];
    accuracy = 0;
    switch nargin
        case 4
            tree_para = 10;
        case 5
            tree_para = n_trees;
        otherwise
            fprintf('++++++Fatal error!Please check the input!');
            return
	end
	factor = TreeBagger(tree_para,x_train,y_train);
    [predict_labels,scores] = predict(factor,x_test);
    predictY = str2num(char(predict_labels));
    accuracy = mean(predictY == y_test);
    predict_label = predictY;
end

author: jindongwang

repo: activityrecognition


text:
function layer = genNetworkLSTMEncoderDecoder2(para)
if isfield(para, 'LastActivation4MSE')==0
    para.LastActivation4MSE = 'linear';
end

layer{1}.name = 'Input';        % this is an input layer
layer{end}.inputIdx = 1;    % specifies the index of GCC in Visible_tr
layer{end}.dim = [1 1]*double(para.inputDim);             % [input dim; output dim];

if isfield(para, 'useFbank') && para.useFbank
    MelWindow = mel_window_FE(40, 256, 16000)';
    MelWindow(:,end+1) = 0;
    layer{end+1}.name = 'Affine';
    layer{end}.prev = -1;
    layer{end}.W = MelWindow;
    layer{end}.b = zeros(40,1);
    layer{end}.dim = [40 layer{end-1}.dim(1)];
    layer{end}.update = 0;

    layer{end+1}.name = 'log';
    layer{end}.prev = -1;
    layer{end}.const = 1e-7;
    layer{end}.dim = [1 1]*layer{end-1}.dim(1);
end

layer{end+1}.name = 'delta';
layer{end}.prev = -1;
layer{end}.dim = [3 1]*layer{end-1}.dim(1);

layer{end+1}.name = 'Affine';
layer{end}.prev = -1;
layer{end}.W = [];
layer{end}.b = [];
layer{end}.dim = [1 1]*layer{end-1}.dim(1);
layer{end}.update = 0;

layerLSTM = genNetworkLSTM(para);
for i=1:length(layerLSTM)
    layerLSTM{i} = rmfield(layerLSTM{i}, 'next');
end

layerLSTM = layerLSTM(2:end);

layerLSTM{1}.prev = -1;
layerLSTM{1}.dim(2) = layer{end}.dim(1);

layer = [layer layerLSTM];

layer{end-2}.name = 'multi_softmax';
layer{end-2}.TaskVocabSizes = para.outputDim/para.ngram_len;
if isfield(para, 'labelDelay')
    shiftLayer.name = 'frame_shift';
    shiftLayer.prev = -1;
    shiftLayer.delay = para.labelDelay;

end

layer{end}.name = 'multi_cross_entropy';
layer{end}.TaskVocabSizes = para.outputDim/para.ngram_len;
layer = FinishLayer(layer);
end

author: singaxiong

repo: SignalGraph


text:
%  rouwenhorst approximates and AR(1) process with a Markov chain
%
%  ::
%
%    [Thetai,y]=rouwenhorst(mu,rho,sigma)
%    [Thetai,y]=rouwenhorst(mu,rho,sigma,N)
%    [Thetai,y]=rouwenhorst(mu,rho,sigma,N,p)
%    [Thetai,y]=rouwenhorst(mu,rho,sigma,N,p,q)
%
%  Args:
%
%     mu (numeric): unconditonal mean of the ar(1) process
%     rho [numeric): autoregressive coefficient
%     sigma (numeric): standard deviation of the shock
%     N (numeric | {2}): number of states of the markov chain
%     p (numeric | {.5*(1+rho)}): first parameter of the procedure
%     q (numeric | {.5*(1+rho)}): second parameter of the procedure
%
%  Returns:
%     :
%
%     - **Thetai** [numeric] : NxN transition matrix
%     - **y** [numeric] : vector of nodes or states
%
%  Note:
%
%     - The process is assumed to be of the form
%       x_t-mu=rho*(x_{t-1}-mu)+sigma*error_term where error_term ~ N(0,1)
%
%  Reference:
%
%     - Karen A. Kopecky, Richard M. H. Suen (2010): "Finite state
%       Markov-chain approximations to highly persistent processes". Review of
%       Economic Dynamics 13, pp 701-714
%
%

author: jmaih


text:
function f = transpose( f )
% .'   SEPARABLEAPPROX transpose.
%    F.' is the non-conjugate transpose of a F.
%
% See also CTRANSPOSE.

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

% Empty check:
if ( isempty( f ) )
    return
end

% Swap columns and rows.
temp = f.cols;
f.cols = f.rows;
f.rows = temp;

% Adjust the domain
f.domain = [f.domain(3:4) f.domain(1:2)];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef StepBase < matlab.hwmgr.internal.hwsetup.TemplateBase

  properties
    InfoText1
    InfoText2
    InfoText3
  end

  methods
    function obj = StepBase(varargin)
      obj@matlab.hwmgr.internal.hwsetup.TemplateBase(varargin{:});

      parentContainer = obj.ContentPanel;

      obj.InfoText1 = matlab.hwmgr.internal.hwsetup.Label.getInstance(parentContainer);
      obj.InfoText1.Position = [20 290 400 80];
      obj.InfoText1.Text = '';
      obj.InfoText1.Color = matlab.hwmgr.internal.hwsetup.util.Color.WHITE;

      obj.InfoText2 = matlab.hwmgr.internal.hwsetup.Label.getInstance(parentContainer);
      obj.InfoText2.Position = [20 230 400 50];
      obj.InfoText2.Text = '';
      obj.InfoText2.Color = matlab.hwmgr.internal.hwsetup.util.Color.WHITE;

      obj.InfoText3 = matlab.hwmgr.internal.hwsetup.Label.getInstance(parentContainer);
      obj.InfoText3.Position = [20 170 400 50];
      obj.InfoText3.Text = '';
      obj.InfoText3.Color = matlab.hwmgr.internal.hwsetup.util.Color.WHITE;

      obj.HelpText.AboutSelection = '';
      obj.HelpText.WhatToConsider = '';

      if strcmp(obj.Workflow.FirstScreenID, class(obj))
        obj.BackButton.Visible = 'off';
        obj.BackButton.Enable = 'off';
      end
    end

    function restoreScreen(obj)
      obj.enableScreen();
    end
  end

%   methods (Access = protected, Abstract)
%     status = installDriverImpl(obj)
%   end

end

author: analogdevicesinc


text:
%% Epipolar Geometry
%
% In this sample:
%
% * We will learn about the basics of multiview geometry
% * We will see what is epipole, epipolar lines, epipolar constraint etc.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/da/de9/tutorial_py_epipolar_geometry.html>
%

%% Theory
%
% When we take an image using pin-hole camera, we loose an important
% information, i.e depth of the image, or how far is each point in the image
% from the camera because it is a 3D-to-2D conversion. So it is an important
% question whether we can find the depth information using these cameras. And
% the answer is to use more than one camera. Our eyes works in similar way
% where we use two cameras (two eyes) which is called stereo vision. So let's
% see what OpenCV provides in this field.
%
% (_Learning OpenCV_ by Gary Bradsky has a lot of information in this field.)
%
% Before going to depth images, let's first understand some basic concepts in
% multiview geometry. In this section we will deal with epipolar geometry. See
% the image below which shows a basic setup with two cameras taking the image
% of same scene.
%
% <<https://docs.opencv.org/3.2.0/epipolar.jpg>>
%
% If we are using only the left camera, we can't find the 3D point
% corresponding to the point $x$ in image because every point on the line
% $OX$ projects to the same point on the image plane. But consider the right
% image also. Now different points on the line $OX$ projects to different
% points ($x'$) in right plane. So with these two images, we can triangulate
% the correct 3D point. This is the whole idea.
%
% The projection of the different points on $OX$ form a line on right plane
% (line $l'$). We call it *epiline* corresponding to the point $x$. It means,
% to find the point $x$ on the right image, search along this epiline. It
% should be somewhere on this line. (Think of it this way, to find the
% matching point in other image, you need not search the whole image, just
% search along the epiline. So it provides better performance and accuracy).
% This is called *Epipolar Constraint*. Similarly all points will have its
% corresponding epilines in the other image. The plane $XOO'$ is called
% *Epipolar Plane*.
%
% $O$ and $O'$ are the camera centers. From the setup given above, you can see
% that projection of right camera $O'$ is seen on the left image at the point,
% $e$. It is called the *epipole*. Epipole is the point of intersection of
% line through camera centers and the image planes. Similarly $e'$ is the
% epipole of the left camera. In some cases, you won't be able to locate the
% epipole in the image, they may be outside the image (which means, one camera
% doesn't see the other).
%
% All the epilines pass through its epipole. So to find the location of
% epipole, we can find many epilines and find their intersection point.
%
% So in this session, we focus on finding epipolar lines and epipoles. But to
% find them, we need two more ingredients, *Fundamental Matrix (F)* and
% *Essential Matrix (E)*. Essential Matrix contains the information about
% translation and rotation, which describe the location of the second camera
% relative to the first in global coordinates. See the image below
% (Image courtesy: _Learning OpenCV_ by Gary Bradsky):
%
% <<https://docs.opencv.org/3.2.0/essential_matrix.jpg>>
%
% But we prefer measurements to be done in pixel coordinates, right?
% Fundamental Matrix contains the same information as Essential Matrix in
% addition to the information about the intrinsics of both cameras so that we
% can relate the two cameras in pixel coordinates. (If we are using rectified
% images and normalize the point by dividing by the focal lengths, $F=E$). In
% simple words, Fundamental Matrix F, maps a point in one image to a line
% (epiline) in the other image. This is calculated from matching points from
% both the images. A minimum of 8 such points are required to find the
% fundamental matrix (while using 8-point algorithm). More points are
% preferred and use RANSAC to get a more robust result.
%

%% Code
%
% First we need to find as many possible matches between two images to find
% the fundamental matrix. For this, we use SIFT descriptors with FLANN based
% matcher and ratio test.
%
% Next we find the Fundamental Matrix from the list of best matches from both
% the images.
%
% Then we find the epilines. Epilines corresponding to the points in first
% image are drawn on second image. So mentioning of correct images are
% important here. We get an array of lines. So we define a new function to
% draw these lines on the images.
%
% Below is the result we get:
%
% <<https://docs.opencv.org/3.2.0/epiresult.jpg>>
%
% You can see in the left image that all epilines are converging at a point
% outside the image at right side. That meeting point is the epipole.
%
% For better results, images with good resolution and many non-planar points
% should be used.
%
% Notes:
%
% * One important topic is the forward movement of camera. Then epipoles will
%   be seen at the same locations in both with epilines emerging from a fixed
%   point. See
%   <http://answers.opencv.org/question/17912/location-of-epipole/ this>.
% * Fundamental Matrix estimation is sensitive to quality of matches, outliers
%   etc. It becomes worse when all selected matches lie on the same plane. See
%   <http://answers.opencv.org/question/18125/epilines-not-correct/ this>.
%

function epipolar_geometry_demo()
    %%
    % a pair of stereo images (grayscale)
    img1 = cv.imread(fullfile(mexopencv.root(),'test','books_left.jpg'), ...
        'Grayscale',true);  % query image
    img2 = cv.imread(fullfile(mexopencv.root(),'test','books_right.jpg'), ...
        'Grayscale',true); % train image

    %%
    % detect keypoints and calculate descriptors using SIFT
    obj = cv.SIFT('ConstrastThreshold',0.03);
    [kp1,desc1] = obj.detectAndCompute(img1);
    [kp2,desc2] = obj.detectAndCompute(img2);

    %%
    % match descriptors using FLANN
    matcher = cv.DescriptorMatcher('FlannBasedMatcher', ...
        'Index',{'KDTree', 'Trees',5}, 'Search',{'Checks',50});
    m = matcher.knnMatch(desc1, desc2, 2);

    %%
    % keep only "good" matches (ratio test as per Lowe's paper)
    m = cat(1, m{:});    % N-by-2 array of structs
    idx_good = ([m(:,1).distance] < 0.8*[m(:,2).distance]);
    m = m(idx_good,1);

    %%
    % extract keypoints from filtered matches
    pts1 = cat(1, kp1([m.queryIdx]+1).pt);
    pts2 = cat(1, kp2([m.trainIdx]+1).pt);
    if true
        pts1 = int32(pts1);
        pts2 = int32(pts2);
    end

    %%
    % find Fundamental matrix
    [F,mask] = cv.findFundamentalMat(pts1, pts2, 'Method','LMedS');
    mask = logical(mask);

    %%
    % select only inlier points
    pts1 = pts1(mask,:);
    pts2 = pts2(mask,:);

    %%
    % random colors to draw matching points and lines
    clrs = randi([0 255], [size(pts1,1) 3], 'uint8');
    clrs(:,4) = 0;

    %%
    % find epilines corresponding to points in right image (second image)
    % and draw its lines on left image
    lines1 = cv.computeCorrespondEpilines(pts2, F, 'WhichImage',2);
    [img11,~] = drawlines(img1, img2, lines1, pts1, pts2, clrs);

    %%
    % find epilines corresponding to points in left image (first image)
    % and draw its lines on right image
    lines2 = cv.computeCorrespondEpilines(pts1, F, 'WhichImage',1);
    [img22,~] = drawlines(img2, img1, lines2, pts2, pts1, clrs);

    %%
    % show result
    if ~mexopencv.isOctave() && mexopencv.require('images')
        imshowpair(img11, img22, 'montage')
        title('Left/Right')
    else
        figure('Position',[200 200 1200 400])
        subplot(121), imshow(img11), title('Left')
        subplot(122), imshow(img22), title('Right')
    end
end

function [img1,img2] = drawlines(img1, img2, lines, pts1, pts2, clrs)
    %DRAWLINES  Draw epilines and points on images
    %
    %     [img1,img2] = drawlines(img1, img2, lines, pts1, pts2, clrs)
    %
    % ## Input
    % * __img1__ first image
    % * __img2__ second image
    % * __lines__ epilines corresponding to `pts2` in `img2`
    % * __pts1__ points in `img1`
    % * __pts2__ points in `img2`
    % * __clrs__ color of each line and matching points
    %
    % ## Output
    % * __img1__ image with drawn points and epilines for the points in `img2`
    % * __img2__ image with drawn points
    %
    % Epilines corresponding to the points in 1st image is drawn on 2nd image.
    %
    % See also: estimateFundamentalMatrix, epipolarLine, lineToBorderPoints
    %

    % convert to RGB
    img1 = cv.cvtColor(img1, 'GRAY2RGB');
    img2 = cv.cvtColor(img2, 'GRAY2RGB');

    % epilines
    w = size(img1,2);
    N = size(lines,1);
    p1 = int32([zeros(N,1), -lines(:,3)./lines(:,2)]);
    p2 = int32([ones(N,1)*w, -(lines(:,3)+lines(:,1)*w)./lines(:,2)]);
    img1 = cv.line(img1, p1, p2, 'Colors',clrs, 'LineType','AA');

    % matching points
    img1 = cv.circle(img1, pts1, 5, 'Colors',clrs, 'Thickness','Filled');
    img2 = cv.circle(img2, pts2, 5, 'Colors',clrs, 'Thickness','Filled');
end


text:
function varargout = cvx_values( varargin )
varargout = cvx_value( varargin );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [sensor, crash] = checkCrash(car, obstacle_polygon, obstacle_circle)

    %% ----- Initialization

    A = get(car(1),'Vertices');     %Get vertices of the vehicle
    dist_crash = 1.5;               %Sensor distance (1.5 = 15cm) for crash
    obs_crash = false;              %Obstacle crash
    sens_crash = false;             %Crash sensor values are below 1.5.
    status = true;                  %Help variable

    %% ----- Check collision between the vehicle and obstacles

    n_polygons = length(obstacle_polygon);
    n_circles = length(obstacle_circle);
    for i = 1:n_polygons

        B = get(obstacle_polygon(i),'Vertices');

        obs_crash = obstacleCrash(A, B, 'polygon');

        if(obs_crash)
            status = false;
            break;
        end

    end

    if(status)

        for i = 1:n_circles

            C = get(obstacle_circle(i),'Vertices');

            obs_crash = obstacleCrash(A, C, 'circle');

            if(obs_crash)
                break;
            end
        end
    end

    %% ----- Check collision with sensors

    n_sensors = length(car(2:end));
    sensor = zeros(1, n_sensors);

    % Loop through the sensors
    for i = 1:n_sensors

        % Get vertices of the sensors
        D = get(car(i+1),'Vertices');
        sensor_temp = 0;

        % For each sensor, check collision with all the obstacles
        for j = 1:n_circles

            C = get(obstacle_circle(j), 'Vertices');

            % Store all the sensor values in a vector
            [sensor_temp(j), crash] = sensorValues(D, C, 'circle', dist_crash);

            % Check collision
            if(crash)
                sens_crash = true;
                break;
            end
        end

        if( ~sens_crash )

            for k = (j+1):(j + n_polygons)

                B = get(obstacle_polygon(k-j), 'Vertices');

                [sensor_temp(k), crash] = sensorValues(D, B, 'polygon', dist_crash);

                if(crash)

                    sens_crash = true;
                end
            end
        end

        % store the lowest sensor value
        sensor(i) = min(sensor_temp);

    end

    %sensor = round(10.*sensor);
    sensor = 10.*sensor;                % Scale the sensor values with 10
    crash = [obs_crash, sens_crash];    % Store the results in a vector
end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
clear

[solversToUse] = prepareTest('requireOneSolverOf',{'gurobi','ibm_cplex'});

printLevel=1;
solverOK = changeCobraSolver('ibm_cplex','LP');
modelToLoad='Recon3DModel_301';
%modelToLoad='ecoli_core_model';

switch modelToLoad
    case 'ecoli_core_model'
        load('ecoli_core_model')
    case 'Recon3DModel_301'
        load('Recon3DModel_301.mat')

    case 'Recon3D'
        load('Recon3D_301.mat')
        model = Recon3D;
    case 'Recon2betaModel'
        %graphStoich/data/modelCollection/121114_Recon2betaModel.mat
        load 121114_Recon2betaModel.mat
        model=modelRecon2beta121114;
    case 'KEGGMatrix'
        load ~/work/modeling/projects/graphStoich/data/modelCollectionBig/KEGGMatrix.mat
        model=KEGG;
end

if 1
    %[SConsistentMetBool, SConsistentRxnBool, SInConsistentMetBool, SInConsistentRxnBool, unknownSConsistencyMetBool, unknownSConsistencyRxnBool, model, stoichConsistModel] = findStoichConsistentSubset(model, massBalanceCheck, printLevel, fileName, epsilon)
    [SConsistentMetBool, SConsistentRxnBool, SInConsistentMetBool, SInConsistentRxnBool, unknownSConsistencyMetBool, unknownSConsistencyRxnBool, model, stoichConsistModel] = findStoichConsistentSubset(model,0,1);
    if strcmp(modelToLoad,'Recon3DModel_301')
        % load('Recon3DModel_301.mat') should be fully stoichiometrically
        % consistent
        assert(nnz(SConsistentMetBool)==5835)
        assert(nnz(SConsistentRxnBool)==8791)
        assert(nnz(SInConsistentMetBool)==0)
        assert(nnz(SInConsistentRxnBool)==1809)
        assert(nnz(unknownSConsistencyMetBool)==0)
        assert(nnz(unknownSConsistencyRxnBool)==0)
    end
else

    %finds the exchange reactions
    model=findSExRxnInd(model);

    [nMet,nIntRxn]=size(model.S(:,model.SIntRxnBool));
    m=sparse(nMet,15);
    i=1;

    if 1
        if 1
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='gurobi5';
            method.param.Method=1;
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 1
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='gurobi5';
            method.param.Method=2;
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 1
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='gurobi5';
            method.param.Method=3;
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 1
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='gurobi5';
            method.param.Method=4;
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 0
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='mosek_linprog';
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 0
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='mosek';
            method.param.MSK_IPAR_OPTIMIZER='MSK_OPTIMIZER_PRIMAL_SIMPLEX';
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 0
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='mosek';
            method.param.MSK_IPAR_OPTIMIZER='MSK_OPTIMIZER_DUAL_SIMPLEX';
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
            i=i+1;
        end

        if 0 %still debugging this solve with Erling @ Mosek
            %by default, the check for stoichiometric consistency omits the columns
            %of S corresponding to exchange reactions
            clear method
            method.interface='solveCobraLP';
            method.solver='mosek';
            method.param.MSK_IPAR_OPTIMIZER='MSK_OPTIMIZER_INTPNT';
            printLevelZ=3;
            [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevelZ,method);
            i=i+1;
        end
    end

    if 1
        clear method
        method.interface='solveCobraLP';
        method.solver='ibm_cplex';
        [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel+1,method);
        i=i+1;
    end

    if 0
        clear method
        method.interface='solveCobraLP';
        method.solver='pdco';
        [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel+1,method);
        i=i+1;
    end

    if 0
        clear method
        method.interface='cvx';
        method.solver='gurobi';
        [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
        i=i+1;
    end

    if 0
        clear method
        method.interface='cvx';
        method.solver='mosek';
        [inform(i),m(:,i),models{i}]=checkStoichiometricConsistency(model,printLevel,method);
        i=i+1;
    end

    m=full(m);

end

author: opencobra

repo: cobratoolbox


text:
classdef AD9364Tests < HardwareTests

    properties
        uri = 'ip:192.168.2.1';
        author = 'ADI';
    end

    methods(TestClassSetup)
        % Check hardware connected
        function CheckForHardware(testCase)
            Device = @()adi.AD9364.Rx;
            testCase.CheckDevice('ip',Device,testCase.uri(4:end),false);
        end
    end

    methods (Static)
        function estFrequency(data,fs)
            nSamp = length(data);
            FFTRxData  = fftshift(10*log10(abs(fft(data))));
            df = fs/nSamp;  freqRangeRx = (-fs/2:df:fs/2-df).'/1000;
            plot(freqRangeRx, FFTRxData);
        end
    end

    methods (Test)

        function testAD9364Rx(testCase)
            % Test Rx DMA data output
            rx = adi.AD9364.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            [out, valid] = rx();
            rx.release();
            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
        end

        function testAD9364TooManyChannels(testCase)
            rx = adi.AD9364.Rx();
            testCase.verifyError(@SetchanelCount4,?MException);
            function SetchanelCount4
                rx.channelCount = 4;
            end
        end

        function testAD9364RxWithTxDDS(testCase)
            % Test DDS output
            tx = adi.AD9364.Tx('uri',testCase.uri);
            tx.DataSource = 'DDS';
            toneFreq = 5e5;
            tx.DDSFrequencies = repmat(toneFreq,2,2);
            tx.AttenuationChannel0 = -10;
            tx();
            pause(1);
            rx = adi.AD9364.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            rx.kernelBuffersCount = 1;
            for k=1:10
                valid = false;
                while ~valid
                    [out, valid] = rx();
                end
            end
            rx.release();

%             plot(real(out));
%             testCase.estFrequency(out,rx.SamplingRate);
            freqEst = meanfreq(double(real(out)),rx.SamplingRate);

            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(freqEst,toneFreq,'RelTol',0.01,...
                'Frequency of DDS tone unexpected')

        end

        function testAD9364RxWithTxData(testCase)
            % Test Tx DMA data output
            amplitude = 2^15; frequency = 0.12e6;
            swv1 = dsp.SineWave(amplitude, frequency);
            swv1.ComplexOutput = true;
            swv1.SamplesPerFrame = 1e4*10;
            swv1.SampleRate = 3e6;
            y = swv1();

            tx = adi.AD9364.Tx('uri',testCase.uri);
            tx.DataSource = 'DMA';
            tx.EnableCyclicBuffers = true;
            tx.AttenuationChannel0 = -10;
            tx(y);
            rx = adi.AD9364.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            rx.kernelBuffersCount = 1;
            for k=1:10
                valid = false;
                while ~valid
                    [out, valid] = rx();
                end
            end
            rx.release();

%             plot(real(out));
%             testCase.estFrequency(out,rx.SamplingRate);
            freqEst = meanfreq(double(real(out)),rx.SamplingRate);

            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(freqEst,frequency,'RelTol',0.01,...
                'Frequency of ML tone unexpected')
        end
    end

end


text:
% bug_3s.m
% explore an imperfect adjoint that used to exist in 3s

if 0 % rotmex bug, fixed by AY on 2006-3-8:
	im(rotmex(ones(200, 'single'), 40.97561, int32(3)))
end

if ~isvar('G')
	nx = 12;
	ny = 12;
	nz = 8;
	na = 6;

%	mask = []; % for 3s object!
	f.sys_type = sprintf('3s@1,1,1,360,0,%d%s@%s@%s@-%d,%d,%d', ...
			1, ',none', '-', '-', nx, nz, na)

	G = Gtomo3(f.sys_type, [], nx, ny, nz);

	ig = image_geom('nx', nx, 'ny', ny, 'nz', nz, 'dx', 1);
	ig.mask = G.arg.mask;
prompt
end

if ~isvar('A')
	[A Aa] = test_adjoint(G);
	im(131, A), im(132, Aa'), im(133, Aa'-A)
end

jj = find(sum(abs(Aa'-A)) > 0.1);
if isempty(jj), printm 'no bug now!', return, end
jj = jj(4)

a1 = reshape(A(:,jj), nx, nz, na, []);
a2 = reshape(Aa(jj,:)', nx, nz, na, []);

im pl 3 1
im(1, a1, 'A ej'), cbar
im(2, a2, 'A'''' ej'), cbar
im(3, a2-a1, 'diff'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [x,y]=als_prec(mat,rhs,y,niter)
%Computes an approximate rank-1 solution to the linear system
%   [Y]=ALS_PREC(MAT,RHS,X) Computes (approximate) rank-1 solution to the
%   problems MAT*P = RHS, where A is a low-Kronecker rank matrix
%   P=X x Y, RHS is a low-Kronecker rank matrix
%   B is given as a two-2d cell array, the same is for X
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

a=mat{1}; %The first cell
b=mat{2};
R=size(a,3); %The rank
n=size(a,1);
m=size(b,1);
a1=permute(a,[1,3,2]); a1=reshape(a1,[n*R,n]);
b1=permute(b,[1,3,2]); b1=reshape(b1,[m*R,m]);
rhs1=rhs;
if ( nargin < 5 )
   niter=10;
end
u=rhs{1};
v=rhs{2};
rr=size(u,3);
u1=reshape(u,[n*n,rr]);
v1=reshape(v,[m*m,rr]);

for i=1:niter
    %fprintf('i=%d \n',i);
    %at=a1*x;
    %bt=b1*y;
    %Iterate over x
    bt=b1*y;
    %bt is n*R*n, scalar product of
    bt1=reshape(bt,[m,R,m]); bt1=permute(bt1,[1,3,2]); bt1=reshape(bt1,[m*m,R]);
    p=bt1'*bt1; %p matrix is ready
    %Compute the local matrix sum_{a,b} p_{ab} A^{\top}_b A_a,
    %it is A(k,i,b)*A(k,j,a)*p(a,b)
    a2=reshape(a,[n*n,R]); a2=a2*p; %a2 is (k,j,b), sum over (k,b) with (k,i,b)
    a3=reshape(a,[n,n,R]); a3=permute(a3,[2,1,3]); a3=reshape(a3,[n,n*R]);
    a2=reshape(a2,[n,n,R]); a2=permute(a2,[2,1,3]); a2=reshape(a2,[n,n*R]);
    loc_mat=a3*a2'; %Should be valid even for complex numbers
    %For the right hand side
    %Compute the q matrix <b,v> = b is m x m x R, V is m x m x rr
    %the result is R x rr
    q=bt1'*v1; %q is R x rr
    %Now compute right-hand side as
    %q(R,rr)*U(i,k,rr)*A(j,k,R)

    u2=u1*q'; % u2 is (i,k,R)
    u2=reshape(u2,[n,n*R]);
    a2=reshape(a,[n,n*R]);
    rhs=u2*a2';
    x= loc_mat \ rhs;
    %cond(loc_mat)
    %Iterate over y
    at=a1*x;
    %bt is n*R*n, scalar product of
    at1=reshape(at,[n,R,n]); at1=permute(at1,[1,3,2]); at1=reshape(at1,[n*n,R]);
    p=at1'*at1; %p matrix is ready
    %Compute the local matrix sum_{a,b} p_{ab} A^{\top}_b A_a,
    %or the summation: \sum_{a,b,k} p(a,b) A(i,k,a)*A(j,k,b)
    b2=reshape(b,[m*m,R]); b2=b2*p; %a2 is (k,j,b), sum over (k,b) with (k,i,b)
    b3=reshape(b,[m,m,R]); b3=permute(b3,[2,1,3]); b3=reshape(b3,[m,m*R]);
    b2=reshape(b2,[m,m,R]); b2=permute(b2,[2,1,3]); b2=reshape(b2,[m,m*R]);
    loc_mat=b3*b2'; %Should be valid even for complex numbers
    %For the right hand side
    %Compute the q matrix <b,v> = b is m x m x R, V is m x m x rr
    %the result is R x rr
    q=at1'*u1; %q is R x rr
    %Now compute right-hand side as
    %q(R,rr)*U(i,k,rr)*A(j,k,R)

    v2=v1*q'; % u2 is (i,k,R)
    v2=reshape(v2,[m,m*R]);
    b2=reshape(b,[m,m*R]);
    rhs=v2*b2';
    y= loc_mat \ rhs;
    %    norm(tt_matrix(mat)*kron(tt_matrix(x,1e-10),tt_matrix(y,1e-10))-tt_matrix(rhs1))
    %keyboard;
    %norm(tt_matrix(mat)*kron(tt_matrix(x,1e-10),tt_matrix(y,1e-10))-tt_mat
    %rix(rhs1))
end
return

author: oseledets

repo: TT-Toolbox


text:
function tidepl()
    %  tidpl  plots a time depth plot of the seismicity
    %  Stefan Wiemer 5/95
    %
    %TODO delete this, replaced by TimeDepthPlotter.plot(catalog)
    global maepi dep1 dep2 dep3 newt2 ms6 ty1 ty2 ty3 a wex wey
    report_this_filefun(mfilename('fullpath'));

    newcat = a;
    xt2  = [ ];
    meand = [ ];
    er = [];
    ind = 0;

    % Find out of figure already exists
    %
    figNumber = findobj(0,'Tag','time_depth_figure');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(figNumber)

        figNumber = figure_w_normalized_uicontrolunits(...
            'Name','Time Depth',...
            'visible','off',...
            'NumberTitle','off', ...
            'MenuBar','none', ...
            'NextPlot','new', ...
            'Tag','time_depth_figure',...
            'Units','Pixel',  'Position',[wex wey 550 400'])
        hold on
        matdraw

        % Make the menu to change symbol size and type
        %
        symbolmenu = uimenu('Label',' Symbol ');
        SizeMenu = uimenu(symbolmenu,'Label',' Symbol Size ');
        TypeMenu = uimenu(symbolmenu,'Label',' Symbol Type ');
        uimenu(SizeMenu,'Label','3','Callback','ms6 =3;eval(cal6)');
        uimenu(SizeMenu,'Label','6','Callback','ms6 =6;eval(cal6)');
        uimenu(SizeMenu,'Label','9','Callback','ms6 =9;eval(cal6)');
        uimenu(SizeMenu,'Label','12','Callback','ms6 =12;eval(cal6)');
        uimenu(SizeMenu,'Label','14','Callback','ms6 =14;eval(cal6)');
        uimenu(SizeMenu,'Label','18','Callback','ms6 =18;eval(cal6)');
        uimenu(SizeMenu,'Label','24','Callback','ms6 =24;eval(cal6)');

        uimenu(TypeMenu,'Label','dot',...
            'Callback','ty1=''.'';ty2=''.'';ty3=''.'';eval(cal6)');
        uimenu(TypeMenu,'Label','red+ blue o green x',...
            'Callback','ty1=''+'';ty2=''o'';ty3=''x'';eval(cal6)');
        uimenu(TypeMenu,'Label','o','Callback',...
            'ty1=''o'';ty2=''o'';ty3=''o'';eval(cal6)');
        uimenu(TypeMenu,'Label','x','Callback',...
            'ty1=''x'';ty2=''x'';ty3=''x'';eval(cal6)');
        uimenu(TypeMenu,'Label','*',...
            'Callback','ty1=''*'';ty2=''*'';ty3=''*'';eval(cal6)');
        uimenu(TypeMenu,'Label','none','Callback','set(deplo1,''visible'',''off'');set(deplo2,''visible'',''off'');set(deplo3,''visible'',''off''); ');
        cal6 = ...
            [ 'set(deplo1,''MarkerSize'',ms6,''LineStyle'',ty1,''visible'',''on'');',...
            'set(deplo2,''MarkerSize'',ms6,''LineStyle'',ty2,''visible'',''on'');',...
            'set(deplo3,''MarkerSize'',ms6,''LineStyle'',ty3,''visible'',''on'');' ];

    end  % if figure exist
    figure(figNumber)
    p5 = findobj(figNumber,'Tag','time_depth_plot');%assumes we only need one plot
    if isempty(p5)
        orient tall
        rect = [0.15, 0.15, 0.75, 0.65];
        p5 = gca;
        p5.Tag = 'time_depth_plot';
        p5.Position = rect;
    end

    depidx = newt2.Depth<=dep1;
    plotTimeDepthSection(p5, newt2, depidx, [ty1 'b'], 'shallower', 'timedepth_depths1');
    %{
    deplo1 = findobj(p5,'Tag','timedepth_depths1');
    if isempty(deplo1)
        deplo1 =plot(p5, newt2.Date(depidx), -newt2.Depth(depidx),'.b','DisplayName','shallower','Tag','timedepth_depths1');
    else
        set(deplo1,'Xdata',newt2.Date(depidx),'YData',newt2.Depth(depidx));
    end
    set(deplo1,'MarkerSize',ms6,'Marker',ty1)
    %}
    hold on

    depidx = newt2.Depth<=dep2&newt2.Depth>dep1;
    plotTimeDepthSection(p5, newt2, depidx, [ty2 'g'], 'mid-depth', 'timedepth_depths2');

    %deplo2 =plot(p5, newt2.Date(depidx),-newt2.Depth(depidx),'.g','DisplayName','mid-depth','Tag','timedepth_depths2');
    %set(deplo2,'MarkerSize',ms6,'Marker',ty2);

    depidx = dep3 >= newt2.Depth & newt2.Depth > dep2;
    plotTimeDepthSection(p5, newt2, depidx, [ty3 'r'], 'mid-depth', 'timedepth_depths3');

    %deplo3 =plot(p5, newt2.Date(depidx),-newt2.Depth(depidx),'.r', 'DisplayName','deeper','Tag','timedepth_depths3');
    %set(deplo3,'MarkerSize',ms6,'Marker',ty3)

    if ~isempty(maepi)
        if isnumeric(maepi)
            maepi = ZmapCatalog(maepi);
            warning('converted maepi to ZmapCatalog in tidepl');
        end
        pl = findobj(p5,'Tag','timedepth_bigevents');
        if isempty(pl)
            pl =   plot(p5, maepi.Date, -maepi.Depth,'xm','Tag','timedepth_bigevents');
        else
            set(pl,'XData',maepi.Date, 'YData', -maepi.Depth,'Color','m','Marker','x');
        end
        set(pl,'LineWidth',2.0)
    end

    xlabel(p5,'Time in Years','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    ylabel(p5,'Depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

    set(p5,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

    grid
    hold off
    done
end
function plotTimeDepthSection(ax, catalog, idx, linechars, displayname, tag)
    % add a part to the time-depth plot
    % reuses if possible
    global ms6
    lineh = findobj(ax,'Tag', tag);
    if isempty(lineh)
        lineh = plot(ax, catalog.Date(idx), -catalog.Depth(idx),linechars,'DisplayName',displayname,'Tag',tag);
    else

        set(lineh,'Xdata',catalog.Date(idx),'YData',-catalog.Depth(idx), 'color',linechars(2),'Marker',linechars(1));
    end
    set(lineh,'MarkerSize',ms6);%,'Marker',ty1);
end

author: CelsoReyes


text:
function title = p07_title ( )

%*****************************************************************************80
%
%% P07_TITLE returns a title for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'f(x) = Langerman2 function';

  return
end

author: johannesgerer


text:
function [X, info] = IRhtv(A, b, varargin)
% IRhtv Least squares solver with heuristic total variation penalization
%
% options  = IRhtv('defaults')
% [X,info] = IRhtv(A,b)
% [X,info] = IRhtv(A,b,K)
% [X,info] = IRhtv(A,b,options)
% [X,info] = IRhtv(A,b,K,options)
%
% This penalized restarted iteration method incorporates a heuristic TV
% penalization term (it does not produce a strict TV solution).
%
% It is assumed that x represents an n-times-n image such that N = n^2.
%
% IRhtv is a simplified driver for IRrestart, which uses an inner-outer
% iteration scheme. Semi-convergent or hybrid iterative solvers are used
% in the inner iterations, using one of the iterative methods in IRtools
% (e.g., IRhybrid_gmres). In the case of IRhtv, the TV penalization
% is updated at each outer iteration.
%
% The regularization parameter and number of inner iterations influence
% the behavior and convergence of the outer iterations.
%
% With 'defaults' as input returns the default options.  Otherwise outputs
% the iterates specified in K, using max(K) as MaxIter, and using all other
% default options.  With options as input: uses the user-specified options
% and all the other default options.
%
% Inputs:
%  A : either (a) a full or sparse matrix
%             (b) a matrix object that performs the matrix*vector operation
%             (c) user-defined function handle
%  b : right-hand side vector
%  K : (optional) integer vector that specifies which (total) iterates are
%      returned in X; the maximum number of iterations is assumed to be max(K)
%      [ positive integer | vector of positive components ]
%  options : structure with the following fields (optional)
%      x0          - initial guess for the iterations; default = zero vector
%                    [ array | {'none'} ]
%      MaxIterIn   - maximum number of inner iterations
%      MaxIterOut  - maximum number of outer iterations
%      x_true      - true solution; allows us to returns error norms with
%                    respect to x_true at each iteration
%                    [ array | {'none'} ]
%      RegParam    - a value or a method to find the regularization used in
%                    the inner iterations
%                    [non-negative scalar | {'gcv'} | 'discrep' ]
%                    This also determines which stopping rule is used for
%                    the inner iterations.
%                    If 'gcv' is chosen, the inner iteration is stopped when
%                      the GCV function minimum stabilizes or increases
%                      within a certain window of iterations (see 'stopGCV',
%                      'FlatTol' and 'MinTol').
%                     If 'discrep' is chosen, and NoiseLevel is rovided,
%                       then the discrepancy principle is used as stopping
%                       criterion (see 'NoiseLevel' and 'eta').
%      stopGCV      - stopping criterion for the inner iterations when
%                     GCV is used
%                     [ 'GCVvalues' | {'resflat'} ]
%      FlatTol      - tolerance for detecting flatness (stabilization)
%                     in the GCV function as a stopping criterion for the
%                     inner iterations
%                     [ {10^-6} | non-negative scalar ]
%      MinTol       - window of iterations - if the GCV minimum continues
%                     to increase over this window, then the inner
%                     iterations are stopped:
%                     [ {3} | positive integer ]
%      RegMatrix    - priorconditioner for the inner iterations
%                     [ {'identity'} | square nonsingular matrix |
%                     function handle ]
%      NoiseLevel   - norm of noise in rhs divided by norm of rhs
%                     (must be assigned in RegParam is 'discrep')
%                     [ {none} | nonnegative scalar ]
%      eta          - safety factor for the discrepancy principle
%                     [ {1.01} | scalar greater than (and close to) 1 ]
%      RegParam0    - first regularization parameter, used only on the
%                     very first iteration (needed if RegParam is 'discrep')
%                     [ {1} | positive scalar ]
%      stopOut      - stopping criterion for the outer iterations;
%                     [ {'xstab'} | 'Lxstab' | 'regPstab' ]
%      inSolver     - solver to be employed during the inner iterations
%                     [ {'gmres'} | 'lsqr' ]
%      adaptConstr  - approximate constraint or regularization to be
%                     incorporated
%                     [ {'tv'} | 'tvnn' ]
%      nonnegativity - may be used to also impose nonnegativity
%                      (similarly to 'tvnn')
%                      [ 'on' | {'off'} ]
%      IterBar       - shows the progress of the outer iterations
%                      [ {'on'} | 'off' ]
%      NoStopIn      - specifies whether the inner iterations should
%                      proceed after a stopping criterion has been satisfied
%                      [ 'on' | {'off'}]
%      NoStopOut     - specifies whether the outer iterations should
%                      proceed after a stopping criterion is satisfied
%                      [ 'on' | {'off'} ]
%      verbosity     - switch on or off the "verbosity" of the function
%                      [ {'on'} | 'off' ]
% Note: the options structure can be created using the function IRset.
%
% Outputs:
%   X : computed solutions, stored column-wise (at the iterations listed in K)
%   info: structure with the following fields:
%      its          - number of the last computed iteration
%      saved_iterations - iteration numbers of iterates stored in X
%      StopFlag_in  - string that describes the inner stopping condition:
%                       * Stopping criterion of the inner iterations is
%                         never satisfied
%                       * Stopping criterion is satisfied at least once
%                         during the inner iterations
%      StopFlag_out - string that describes the outer stopping condition;
%                     depending on the inputs it can be one of the following:
%                       * Outer stopping criterion is never satisfied
%                       * Diagonal weighting matrix is numerically zero
%                       * Solution stabilizes
%                       * Transformed solution stabilizes
%                       * Regularization parameter stabilizes
%      Rnrm     - relative residual norms at each iteration
%      Xnrm     - solution norms at each iteration
%      Enrm     - relative error norms (requires x_true) at each iteration
%      StopReg  - struct containing information about the solution that
%                 satisfies the stopping criterion.  Fields:
%                   It   : iteration where the stopping criterion is satisfied
%                   X    : solution satisfying the stopping criterion
%                   Enrm : the corresponding relative error (requires x_true)
%      BestReg  - struct containing information about the solution that
%                 minimizes Enrm (requires x_true). Fields:
%                   It   : iteration where the minimum is attained
%                   X    : best solution
%                   Enrm : best relative error
%      Xout     - approximate solutions at the end of each inner cycle,
%                 stored column-wise
%      itsInOut - 3-column matrix whose the columns store
%                   1. outer iteration count
%                   2. inner iteration count (i.e., for each cycle)
%                   3. total iteration count
%
% See also: IRell1, IRirn, IRrestart, IRget, IRset

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('x0', 'none', 'MaxIterIn', 30 , 'MaxIterOut', 20 , ...
    'RegParam', 'gcv', 'stopGCV', 'resflat', ...
    'resflatTol', 0.05, 'GCVflatTol', 10^-6, 'GCVminTol', 3,...
    'x_true', 'none', 'IterBar', 'on', 'NoStop', 'off', 'NoStopIn', 'off',...
    'NoStopOut', 'off', 'stopOut', 'xstab', 'stabOut', 1e-6, ...
    'thr0', 1e-10, 'NoiseLevel', 'none', 'eta', 1.01, 'RegParam0', 1,...
    'inSolver', 'gmres', 'adaptConstr', 'tv', 'nonnegativity', 'off', ...
    'verbosity', 'off');

% If input is 'defaults,' return the default options in X.
if nargin==1 && nargout <= 1 && isequal(A,'defaults')
    X = defaultopt;
    return;
end

% Check for acceptable number of optional input arguments.
switch length(varargin)
    case 0
        K = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = [];
        else
            K = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = varargin{2};
        else
            K = varargin{2}; options = varargin{1};
        end
    otherwise
        error('Too many input parameters')
end

n = length(b(:));
nosquare = 0;
test_sq = ones(n,1);
try
    test_sq = A_times_vec(A, test_sq);
    if (length(test_sq)~=n)
        nosquare = 1;
    end
catch
    nosquare = 1;
end

if isfield(options, 'inSolver') && ~isempty(options.inSolver)
    inSolver   = IRget(options, 'inSolver', [], 'fast');
    if strcmp(inSolver, 'gmres')
        if nosquare
            warning(['The matrix A is rectangular, and a solver like hybrid gmres cannot handle it. ',...
                'The solver is changed to hybrid lsqr.'])
            options.inSolver = 'lsqr';
        end
    end
end

if isempty(options)
    options = defaultopt;
end

options = IRset(defaultopt, options);
inSolver   = IRget(options, 'inSolver', [], 'fast');

if nosquare && strcmp(inSolver, 'gmres')
    options.inSolver = 'lsqr';
end

nn = IRget(options, 'nonnegativity',  [], 'fast');

if strcmp(nn, 'on')
    options.adaptConstr = 'tvnn';
end

% Call IRrestart with the specified options.
options = rmfield(options, 'nonnegativity');
[X, info] = IRrestart(A, b, K, options);

repo: IRtools


text:
function plotFourier(SO3F,varargin)
% visualize the harmonic coefficients
%
% Syntax
%   plotFourier(SO3F)
%   plotFourier(SO3F,'bandwidth',32)
%
% Input
%  SO3F - @SO3Fun
%

warning('The command plotFourier is depreciated! Please use plotSpektra instead.')

plotSpektra(SO3F,varargin)

author: mtex-toolbox

repo: mtex

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


text:
function [ind,t0,s0,t0close,s0close] = crossing(S,t,level,imeth)
% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input.

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% brueckner@sbrs.net

% check the number of input arguments
narginchk(1,4);

% check the time vector input for consistency
if nargin < 2 || isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end

% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); ind0(ind0==1)=[];

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 ); ind1(ind1==1)=[];

% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind);
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        if abs(S(ind(ii))) > eps(S(ind(ii)))
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            DELTA =  NUM / DEN;
            t0(ii) = t0(ii) - S(ind(ii)) * DELTA;
            % I'm a bad person, so I simply set the value to zero
            % instead of calculating the perfect number ;)
            s0(ii) = 0;
        end
    end
end

% Addition:
% Some people like to get the data points closest to the zero crossing,
% so we return these as well
[CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1);
ind2 = ind + (II-2); %update indices

t0close = t(ind2);
s0close = S(ind2);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function exact = p37_exact ( )

%*****************************************************************************80
%
%% P37_EXACT returns the exact integral for problem 37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  alpha = p37_param_get ( );

  exact = atan ( ( 4.0 - pi ) * 4.0^( alpha - 1.0 ) ) ...
        + atan (         pi   * 4.0^( alpha - 1.0 ) );

  return
end

repo: jburkardt-m


text:
% Mathematics Q2812691
% https://math.stackexchange.com/questions/2812691
% Least Squares with Euclidean L2 Norm Constraint
% References:
%   1.  aa
% Remarks:
%   1.  See the trick for Semi Definiteness http://ask.cvxr.com/t/5168
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     16/06/2018
%   *   First release.

%% General Parameters

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;

%% Simulation Parameters

numRows = 4; %<! Parameter 'n' in the question
numCols = 3; %<! Parameter 'k' in the question

%% Generate Data

mA = randn([numRows, numRows]);
mA = mA.' * mA;
mB = randn([numRows, numCols]);
mI = eye(numCols);

mZ = chol(mA); %<! mA = mZ.' * mZ

paramLambda = 0.1;

normConst = 1;

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable mX(numRows, numCols)
    % minimize( norm(mZ * mX, 'fro') + (paramLambda * trace(mX.' * mB)) );
    minimize( norm(mA * mX - mB, 'fro') );
    subject to
        % ((mX.' * mX) - mI) == semidefinite(numCols);
        % norm(mX.' * mX, 2) <= 1;
        [eye(numCols), mX.'; mX, eye(numRows)] == semidefinite(numCols + numRows)
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Solution by KKT Solution

vParamLambda    = [0:0.1:25];
vNormVal        = zeros([length(vParamLambda), 1]);

mAA = mA.' * mA;
mAb = mA.' * vB;
mI  = eye(size(mA, 2));

for ii = 1:length(vParamLambda)
    paramLambda = vParamLambda(ii);

    vNormVal(ii) = norm((mAA + (paramLambda * mI)) \ mAb, 2);

end

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineSeries = plot(vParamLambda, [vNormVal, normConst * ones([length(vParamLambda), 1])]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', ['Tikhonov Regularization Least Squares Solution Norm'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', '\lambda', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'L_2 Norm', ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Solution Norm'], ['Constraint Norm']});
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

vX = SolveLsNormConst(mA, vB, normConst);

objVal = sum((mA * vX - vB) .^ 2);

disp([' ']);
disp(['Projected Gradient Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(objVal)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function layer = vl_nndropout(x, varargin)
%VL_NNDROPOUT Additional options for vl_nndropout (CNN dropout)
%   Y = Layer.vl_nndropout(X) applies dropout to the data X. See help
%   vl_nndropout for more details.
%
%   This method overloads MatConvNet's vl_nndropout function for Layer
%   objects, so that instead of executing vl_nndropout, a new Layer object
%   is returned.
%
%   The dropout layer is composed of a vl_nndropout_wrapper and an
%   auxiliary mask generator layer, vl_nnmask.

% Copyright (C) 2016-2017 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  opts.rate = 0.5 ;
  opts = vl_argparse(opts, varargin) ;

  % create mask generator layer
  maskLayer = Layer(@vl_nnmask, x, opts.rate) ;

  % create dropout wrapper layer
  layer = Layer(@vl_nndropout_wrapper, x, maskLayer, Input('testMode')) ;

  % vl_nndropout_wrapper doesn't return a derivative for the mask
  layer.numInputDer = 1 ;

end

repo: MFT


text:
function [mssim, ssim_map, mcs, cs_map] = ssim_index_new(img1, img2, K, window)

if (nargin < 2 | nargin > 4)
   ssim_index = -Inf;
   ssim_map = -Inf;
   return;
end

if (size(img1) ~= size(img2))
   ssim_index = -Inf;
   ssim_map = -Inf;
   return;
end

[M N] = size(img1);

if (nargin == 2)
   if ((M < 11) | (N < 11))
	   ssim_index = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian', 11, 1.5);	%
   K(1) = 0.01;										% default settings
   K(2) = 0.03;										%
end

if (nargin == 3)
   if ((M < 11) | (N < 11))
	   ssim_index = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian', 11, 1.5);
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   ssim_index = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   ssim_index = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 4)
   [H W] = size(window);
   if ((H*W) < 4 | (H > M) | (W > N))
	   ssim_index = -Inf;
	   ssim_map = -Inf;
      return
   end
   if (length(K) == 2)
      if (K(1) < 0 | K(2) < 0)
		   ssim_index = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   ssim_index = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

C1 = (K(1)*255)^2;
C2 = (K(2)*255)^2;
window = window/sum(sum(window));

mu1   = filter2(window, img1, 'valid');
mu2   = filter2(window, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 & C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
   cs_map = (2*sigma12 + C2)./(sigma1_sq + sigma2_sq + C2);
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;

   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);

   cs_map = ones(size(mu1));
   index = denominator2 > 0;
   cs_map(index) = numerator2(index)./denominator2(index);
end

mssim = mean2(ssim_map);
mcs = mean2(cs_map);

return

author: thomas-koehler

repo: SupER


text:
function [newVal,newUnitStr,oldUnitStr] = unitConvert(oldVal,valueType,oldUnitStr,newUnitStr)
% [newVal,newUnitStr,oldUnitStr] = unitConvert(oldVal,valType,oldUnitStr)
%
% Value types:
%  'length' - recognized units: 'nm', 'um', 'mm', 'cm'
%  'time'   - recognized units: 'usec', 'msec', 'sec'
%
% Adapted from isetbio unitConvert
%
% Examples
%
% newVal = unitConvert(1500,'time','ms' ,'s')
% newVal = unitConvert(1,'length','cm' ,'mm')

% Switch on value type
switch (valueType)
    case 'length'
        % Length

        % Factor from old to default
        switch oldUnitStr
            case 'm'
                oldConversionFactor = 1;
            case 'nm'
                oldConversionFactor = 1e-9;
            case 'um'
                oldConversionFactor = 1e-6;
            case 'mm'
                oldConversionFactor = 1e-3;
            case 'cm'
                oldConversionFactor = 1e-2;
            otherwise
                error('Bad units %s passed for type %s',oldUnitStr,valueType);
        end

         % Factor from old to default
        switch newUnitStr
            case 'm'
                newConversionFactor = 1;
            case 'nm'
                newConversionFactor = 1e9;
            case 'um'
                newConversionFactor = 1e6;
            case 'mm'
                newConversionFactor = 1e3;
            case 'cm'
                newConversionFactor = 1e2;
            otherwise
                error('Bad units %s passed for type %s',oldUnitStr,valueType);
        end

    case 'time'
        % Time

        % Factor from old to default
        switch oldUnitStr
            case {'s' 'sec' 'seconds'}
                oldConversionFactor = 1;
            case {'ms' 'msec' 'milliseconds'}
                oldConversionFactor = 1e-3;
            otherwise
                error('Bad units %s passed for type %s',oldUnitStr,valueType);
        end

         % Factor from old to default
        switch newUnitStr
            case {'s' 'sec' 'seconds'}
                newConversionFactor = 1;
            case {'ms' 'msec' 'milliseconds'}
                newConversionFactor = 1e3;
            otherwise
                error('Bad units %s passed for type %s',oldUnitStr,valueType);
        end

    otherwise
        error('Unknown value type passed');
end

% Convert
newVal = oldConversionFactor*newConversionFactor*oldVal;

end

author: vistalab

repo: vistasoft


text:
function geometry_test1836 ( )

%*****************************************************************************80
%
%% GEOMETRY_TEST1836 tests SPHERE_EXP2IMP_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2011
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  p = [ ...
    4.0, 2.0, 3.0;
    1.0, 5.0, 3.0;
    1.0, 2.0, 6.0;
   -2.0, 2.0, 3.0 ]';
  pc_true = [ 1.0, 2.0, 3.0 ]';
  r_true = 3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEOMETRY_TEST1836\n' );
  fprintf ( 1, '  SPHERE_EXP2IMP_ND: explicit sphere => implicit form;\n' );

  r8mat_transpose_print ( n, n + 1, p, '  Initial form of explicit sphere:' );

  [ r, pc ] = sphere_exp2imp_nd ( n, p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed form of implicit sphere:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Imputed radius = %f\n', r );
  fprintf ( 1, '  True radius =    %f\n', r_true );

  r8vec_print ( n, pc, '  Imputed center' );

  r8vec_print ( n, pc_true, '  True center' );

  return
end


text:
function [pass, u1, u2, info1, info2] = test_scalarODEdamping(pref)
% A nonlinear CHEBOP test. This test tests a scalar ODE, where no breakpoints
% occur. It solves the problem using chebcolloc1, chebcolloc2 and ultraS
% discretizations. The problem solved requires damping for the Newton iteration
% to converge.
%
% Asgeir Birkisson, May 2014.

%% Setup
dom = [0 pi];
if ( nargin == 0 )
    pref = cheboppref;
end

N = chebop(@(x,u) .05*diff(u,2) + cos(5*x).*sin(u), dom);
N.lbc = @(u) u - 2;
N.rbc = @(u) u - 3;
rhs = 0;

pref.bvpTol = 1e-10;

% Try different discretizations:

%% Start with chebcolloc2
pref.discretization = @chebcolloc2;
pref.bvpTol = 1e-13;
[u1, info1] = solvebvp(N, rhs, pref);
err(1) = norm(N(u1));

%% Change to chebcolloc1
pref.discretization = @chebcolloc1;
[u2, info2] = solvebvp(N, rhs, pref);
err(2) = norm(N(u2));

%% Change to ultraS
pref.discretization = @ultraS;
pref.bvpTol = 1e-12;
[u3, info3] = solvebvp(N, rhs, pref);
err(3) = norm(N(u3));

%% Did we pass?
% To pass, both residuals have to be small, but we should not expect u1 and u2
% to be identical!

tol = 1e-9;
pass = err < tol;

end


text:
%This function is used for mapping the fixation record into the corresponding fixation maps.
screen_res_x = 1440;
screen_res_y = 900;

datasetFile1 = 'movie';
datasetFile = 'videodata';
gazeFile = 'gaze_record';

videoFiles = dir(fullfile('./', datasetFile));
videoNUM = length(videoFiles)-2;
rate = 30;

 for videonum = 1:700
        videofolder =  videoFiles(videonum+2).name
        options.infolder = fullfile( './', datasetFile,  videofolder, 'images' );
        % Cache all frames in memory
        [data.frames,names,video_res_y,video_res_x,nframe ]= readAllFrames( options );
        a=video_res_x/screen_res_x;
        b=(screen_res_y-video_res_y/a)/2;
        all_fixation = zeros(video_res_y,video_res_x,nframe);
        for person = 1:17
            txtloc = fullfile('./', gazeFile, sprintf('P%03d',person), [sprintf('P%03d_Trail',person), sprintf('%04d.txt',videonum)]);
            if exist(txtloc, 'file')
                [time,x_screen,y_screen,event]=textread(txtloc,'%f%f%f%s','headerlines',1);
                if size(time,1)
                    time = time-time(1);
                    event = cellfun(@(x) x(1), event);
                    for index = 1:nframe
                            eff = find( ((index-1)<rate*time/1000000)&(rate*time/1000000<index)&event=='F'); %framerate = 10;
                            x_stimulus=int32(a*x_screen(eff));
                            y_stimulus=int32(a*(y_screen(eff)-b));
                            t = x_stimulus<=0|x_stimulus>=video_res_x|y_stimulus<=0|y_stimulus>=video_res_y;
                            all_fixation(y_stimulus(~t),x_stimulus(~t),index) = 1;
                    end
                end
            end
        end
end

author: wenguanwang


text:
% TD | Tucker-ALS | Tucker Decomposition solved by Alternating Least Squares
% process_video('TD', 'Tucker-ALS', 'dataset/demo.avi', 'output/demo_Tucker-ALS.avi');

r = [size(T,1) size(T,2) 1];
A = double(T);
L = double(tucker_als(T,r));
S = A - L;

author: andrewssobral

repo: lrslibrary


text:
function addToExecutionWarnings(warnStr, eventNum, bodyID, celBodyData)
%addToExecutionWarnings Summary of this function goes here
%   Detailed explanation goes here
    global script_execution_warnings;

    if(bodyID == -1)
        bNameStr = '';
    else
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
        bName = bodyInfo.name;
        bNameStr = [' (',bName ')'];
    end

    if(isempty(eventNum) || eventNum == -1)
        eventStr = '';
    else
        eventStr = ['Event ', num2str(eventNum), bNameStr, ': '];
    end

    execWarn = script_execution_warnings;
    execWarn{end+1} = [eventStr, warnStr];
    script_execution_warnings = execWarn;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot