text:
function model = vargplvmInitDynamics(model,optionsDyn)
% VARGPLVMINITDYNAMICS Initialize the dynamics of a var-GPLVM model.
% FORMAT
% ARG optionsDyn : the VARGPLVM structure with the options used to
% create the vargplvm model..
%
%
% COPYRIGHT : Michalis K. Titsias, 2011
% COPYRIGHT : Neil D. Lawrence, 2011
% COPYRIGHT : Andreas C. Damianou, 2011
%
% SEEALSO : vargplvmOptionsDyn, vargplvmAddDynamics

% VARGPLVM

if ~isfield(model, 'dynamics') || isempty(model.dynamics)
    error(['model does not have field dynamics']);
end

%model.dynamics.kern.comp{1}.inverseWidth = 200./(((max(model.dynamics.t)-min(model.dynamics.t))).^2);
%params = vargplvmExtractParam(model);
%model = vargplvmExpandParam(model, params);

% Initialize barmu: first initialize X (mu) and then (see lines below)
% initialize barmu. The first 'if' statement checks if optionsDyn.initX is
% the initial X itself (e.g. calculated before in the demo) and there is no
% need to recalculate it.
if isfield(optionsDyn,'initX') && ~isstr(optionsDyn.initX)
    X = optionsDyn.initX;
else
    if isfield(optionsDyn,'initX')
        initFunc = str2func([optionsDyn.initX 'Embed']);
    else
        initFunc=str2func('ppcaEmbed');
    end

    % If the model is in "D greater than N" mode, then we want initializations
    % to be performed with the original model.m
    if isfield(model, 'DgtN') && model.DgtN
        X = initFunc(model.mOrig, model.q);
    else
        X = initFunc(model.m, model.q);
    end
end

%---------------------- TEMP -------------------------------%
% The following are several tries to make the init. of the means better.
% It does that by "squeezing" the ends of the means to values closer to
% zero.

% X(1,:)=zeros(size(X(1,:)));
% X(end,:)=zeros(size(X(1,:)));

% X = 0.1*randn(size(X));
% X(1,:)=zeros(size(X(1,:)));
% X(end,:)=zeros(size(X(1,:)));

% m = floor((size(X,1))/2);
% p=max(max(abs(X))); p=p/10; %p = p^(1/m);
% l = logspace(p,0,m);
% l2 = logspace(p,0,m);
% l2 = l2(end:-1:1);
% l = [l l2]';
% mask = repmat(l,1,size(X,2));
% X = X.*(1./mask);
% X(1,:) = zeros(size(X(1,:)));
% X(end,:) = zeros(size(X(1,:)));
% X(end-1,:) = zeros(size(X(1,:)));

%save 'TEMPX.mat' 'X' %%TEMP

if optionsDyn.regularizeMeans
    model.dynamics.regularizeMeans = 1;
    m = floor((size(X,1))/7);  % 7
    p=ones(1,model.q);
    mm=max(X).^(1/2); %1/2 % the smallest this ratio(the exponent) is, the largest the effect
    for i=0:m-4
        X(end-i,:) = X(end-i,:) -  X(end-i,:).*p;
        p = p./mm;
    end
    p=ones(1,model.q);
    for i=1:m
        X(i,:) = X(i,:) -  X(i,:).*p;
        p = p./mm;
    end
end

% The initial covariances we aim for, via calibrating the reparametrised
% covariances and measuring their median
if ~isfield(optionsDyn, 'initCovarMedian') || isempty(optionsDyn.initCovarMedian)
    initCovarMedian = 0.18;
else
    initCovarMedian = optionsDyn.initCovarMedian;
end
% The lowest acceptable value of the above
if ~isfield(optionsDyn, 'initCovarMedianLowest') || isempty(optionsDyn.initCovarMedianLowest)
    initCovarMedianLowest = initCovarMedian/1.8;
else
    initCovarMedianLowest = optionsDyn.initCovarMedianLowest;
end
% The highest acceptable value of the above
if ~isfield(optionsDyn, 'initCovarMedianHighest') || isempty(optionsDyn.initCovarMedianHighest)
    initCovarMedianHighest = initCovarMedian/0.3;
else
    initCovarMedianHighest = optionsDyn.initCovarMedianHighest;
end

%model.X = X;

% m = floor((size(X,1))/2);
% l = logspace(10,1,m);
% l2 = logspace(10,1,m);
% l2 = l2(end:-1:1);
% mask = repmat([l l2]',1,size(X,2));
% X2 = X-X.*mask;

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

vX = var(X);
noise = 0;%.01;
for q=1:model.q
    Lkt = chol(model.dynamics.Kt + noise*vX(q)*eye(model.N))';
    % barmu = inv(Kt + s*I)*X, so that  mu = Kt*barmu =  Kt*inv(Kt +
    % s*I)*X, which is jsut the GP prediction, is temporally smoothed version
    % of the PCA latent variables X (for the data Y)
    model.dynamics.vardist.means(:,q) = Lkt'\(Lkt\X(:,q));
end

if isfield(optionsDyn, 'vardistCovars') && ~isempty(optionsDyn.vardistCovars)
    if length(optionsDyn.vardistCovars) ==1
        model.dynamics.vardist.covars = 0.1*ones(model.N,model.q) + 0.001*randn(model.N,model.q);
        model.dynamics.vardist.covars(model.dynamics.vardist.covars<0.05) = 0.05;
        model.dynamics.vardist.covars = optionsDyn.vardistCovars * ones(size(model.dynamics.vardist.covars));
    elseif size(optionsDyn.vardistCovars) == size(model.dynamics.vardist.covars)
        model.dynamics.vardist.covars = optionsDyn.vardistCovars;
    end
else % Try to automatically find a satisfactory value for the initial reparametrized
     % covariances (lambda) so that the median of the real covariances is
     % as close as possible to an ideal value, eg 0.18
    fprintf(1, '# Finding an initial value for the reparametrized covariance, aiming at median(S)=%.2f...',initCovarMedian)
    % Try an initial calibration in the default search region of the
    % parameters
    [bestVdc, bestMed] = calibrateLambda(model, 0.05, 0.05, 3, initCovarMedian);
    if bestMed < initCovarMedianLowest || bestMed > initCovarMedianHighest
        % The above search failed. Try some other intervals for the search
        % parameter.
        fprintf(1, '.')
        [bestVdcNew, bestMedNew] = calibrateLambda(model, 0.0008 , 0.05, 0.05, initCovarMedian);
        if abs(bestMedNew - initCovarMedian) < abs(bestMed - initCovarMedian)
            bestMed = bestMedNew;
            bestVdc = bestVdcNew;
        end
    end
    if bestMed < initCovarMedianLowest || bestMed > initCovarMedianHighest
        % Repeat as above for a different interval.
        fprintf(1, '.')
        [bestVdcNew, bestMedNew] = calibrateLambda(model, 3, 0.05, 6.5, initCovarMedian);
        if abs(bestMedNew - initCovarMedian) < abs(bestMed - initCovarMedian)
            bestMed = bestMedNew;
            bestVdc = bestVdcNew;
        end
    end
    fprintf(1, '\n')
    model.dynamics.vardist.covars = 0.1*ones(model.N,model.q) + 0.001*randn(model.N,model.q);
    model.dynamics.vardist.covars(model.dynamics.vardist.covars<0.05) = 0.05;
    model.dynamics.vardist.covars = bestVdc * ones(size(model.dynamics.vardist.covars));
    model = vargplvmDynamicsUpdateStats(model);
    fprintf('# Automatically calibrated lambda so that median(covars)=%.2f\n',median(model.vardist.covars(:)))
end

% smaller lengthscales for the base model
%model.kern.comp{1}.inputScales = 5./(((max(X)-min(X))).^2);
params = vargplvmExtractParam(model);
model = vargplvmExpandParam(model, params);

if isfield(optionsDyn,'X_u')
    model.X_u = optionsDyn.X_u;
else
    % inducing point need to initilize based on model.vardist.means
    if model.k <= model.N % model.N = size(mode.vardist.means,1)
        if ~isfield(optionsDyn, 'labels')
            ind = randperm(model.N);
            ind = ind(1:model.k);
            model.X_u = model.X(ind, :);
        else
            % in the case that class labels are supplied, make sure that inducing inputs
            % from all classes are chosen
            [idcs, nSmpls] = class_samples( optionsDyn.labels, model.k );

            count = 1;
            midx = [];
            for inds = idcs
                ind   = inds{:};
                ind   = ind(randperm(numel(ind)));
                idx  = ind(1:nSmpls(count));

                % test that there is no overlap between index sets
                assert(isempty(intersect(midx, idx)));
                midx = [midx, idx];

                count = count+1;
            end
            model.X_u = model.X(midx,:);
        end
    else
        samplingInd=0; %% TEMP
        if samplingInd
            % !!! We could also try to sample all inducing points (more uniform
            % solution)
            % This only works if k<= 2*N
            model.X_u=zeros(model.k, model.q);
            ind = randperm(model.N);
            %ind = ind(1:model.N);
            model.X_u(1:model.N,:) = model.X(ind, :);

            % The remaining k-N points are sampled from the (k-N) first
            % distributions of the variational distribution (this could be done
            % randomly as well).
            dif=model.k-model.N;
            model.X_u(model.N+1:model.N+dif,:)=model.vardist.means(1:dif,:) + rand(size(model.vardist.means(1:dif,:))).*sqrt(model.vardist.covars(1:dif,:));  % Sampling from a Gaussian.
        else
            % !!! The following is not a good idea because identical initial
            % ind.points are difficult to optimise... better add some
            % additional noise.
            model.X_u=zeros(model.k, model.q);
            for i=1:model.k
                %ind=randi([1 size(model.vardist.means,1)]);
                % Some versions do not have randi... do it with rendperm
                % instead:
                ind=randperm(size(model.vardist.means,1));
                ind=ind(1);
                model.X_u(i,:) = model.vardist.means(ind,:)+0.001*randn(1,model.q);
            end
        end
    end
end

if isfield(optionsDyn, 'testReoptimise')
    model.dynamics.reoptimise = optionsDyn.testReoptimise;
end

model.dynamics.learnVariance = optionsDyn.learnVariance; % DEFAULT: 0

if isfield(optionsDyn, 'constrainType') && ~isempty(optionsDyn.constrainType)
    model.dynamics.constrainType = optionsDyn.constrainType;
end

params = vargplvmExtractParam(model);
model = vargplvmExpandParam(model, params);

med = median(model.vardist.covars(:));
minC = min(model.vardist.covars(:));
maxC = max(model.vardist.covars(:));
if med < initCovarMedianLowest || med > initCovarMedianHighest
   warning('!!! [Min Median Max] value of variational covariances is [%1.2f %1.2f %1.2f].\n', minC, med, maxC);
else
    fprintf('# [Min Median Max] value of variational covariances is [%1.2f %1.2f %1.2f].\n', minC, med, maxC);
end
end

function [bestVdc, bestMed] = calibrateLambda(model, from, step, to, initCovarMedian)
    if nargin < 2 || isempty(from)
        from = 0.05;
    end
    if nargin < 3 || isempty(step)
        step = 0.05;
    end
    if nargin < 4 || isempty(to)
        to = 3;
    end
    if nargin < 5 || isempty(initCovarMedian)
        initCovarMedian = 0.18;
    end

    bestVdc = 0.05; bestDiff = 1000;
    for vdc = [from:step:to]
        model.dynamics.vardist.covars = 0.1*ones(model.N,model.q) + 0.001*randn(model.N,model.q);
        model.dynamics.vardist.covars(model.dynamics.vardist.covars<0.05) = 0.05;
        model.dynamics.vardist.covars = vdc * ones(size(model.dynamics.vardist.covars));
        modelNew = vargplvmDynamicsUpdateStats(model);
        med = median(modelNew.vardist.covars(:));
        curDiff = abs(med - initCovarMedian);
        if (curDiff) < bestDiff % "Ideal" median is initCovarMedian
            bestVdc = vdc;
            bestDiff = curDiff;
            bestMed = med;
        end
    end
end

author: SheffieldML

repo: vargplvm


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = RegCoeff_M(S, M, Ns, g, df, B, Nb, Nr)
% Calculates the regression coefficiients

v = g(:,end);   % start for backward induction

f = zeros(Nb,Nr-1);

% backward induction and regression from t_{Nr-1} up to t_1
for i = Nr-1:-1:1
        index = find(g(:,i+1) > 0); % all ITM paths
        s = S(index,i+1);           % values of S at given time points
        m = M(index,i+1);           % values of M at given time points
        v = v * df(i+1);            % option value at t_i

        Acell = B(s,m);             % evaluate basis function in cell array B
        A = cell2mat(Acell{:,:});   % convert to matrix

        f(:,i) = (A'*A)\(A'*v(index));  % determine coefficients
        c = A*f(:,i);                   % continuation value

        exercise = g(index,i+1) >= c;    % early exercise
        v(index(exercise)) = g(index(exercise),i+1);
end

y = f;

end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = rts_mainGUI(varargin)
% RTS_MAINGUI MATLAB code for rts_mainGUI.fig
%      RTS_MAINGUI, by itself, creates a new RTS_MAINGUI or raises the existing
%      singleton*.
%
%      H = RTS_MAINGUI returns the handle to a new RTS_MAINGUI or the handle to
%      the existing singleton*.
%
%      RTS_MAINGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in RTS_MAINGUI.M with the given input arguments.
%
%      RTS_MAINGUI('Property','Value',...) creates a new RTS_MAINGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before rts_mainGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to rts_mainGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help rts_mainGUI

% Last Modified by GUIDE v2.5 09-Sep-2014 18:16:51

% Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @rts_mainGUI_OpeningFcn, ...
                       'gui_OutputFcn',  @rts_mainGUI_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT
end

% --- Executes just before rts_mainGUI is made visible.
function rts_mainGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to rts_mainGUI (see VARARGIN)

% Choose default command line output for rts_mainGUI
handles.output = hObject;

%%% REPLACE THIS WITH A CALL TO GET CEL BODY DATA FROM MAIN KSPTOT GUI %%%
celBodyData = varargin{1};
setappdata(handles.rtsMainGUI,'CelBodyData',celBodyData);
setappdata(handles.rtsMainGUI,'ChildGUIs',{});

rtsHost = varargin{2};
setappdata(handles.rtsMainGUI,'RHost',rtsHost);

mainGUIFig = varargin{3};
setappdata(handles.rtsMainGUI,'mainGUIFig',mainGUIFig);

rtsOpts = struct();
rtsOpts.tmTransRate = '1';
rtsOpts.procExecMode = 'fixedRate';
rtsOpts.processFreq = '0.5';
setappdata(handles.rtsMainGUI,'RTSOptions',rtsOpts);

% Initialize RTS Main Panel
initRTSMainPanel(handles.rtsMainGUI);

% Update handles structure
if(ishandle(hObject))
    guidata(hObject, handles);
end

% UIWAIT makes rts_mainGUI wait for user response (see UIRESUME)
% uiwait(handles.rtsMainGUI);
end

function initRTSMainPanel(hRtsMainGUI)
    rts_StartTMStream(hRtsMainGUI);
end

% --- Outputs from this function are returned to the command line.
function varargout = rts_mainGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    varargout{1} = hObject;
end

% --- Executes on button press in rts_FlightPageButton.
function rts_FlightPageButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_FlightPageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
end

% --- Executes on button press in rts_OrbOpsPageButton.
function rts_OrbOpsPageButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_OrbOpsPageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA).
    vesselGUID = getappdata(handles.rtsMainGUI,'VesselGUID');
    h = rts_OrbOpsConsole(handles.rtsMainGUI, vesselGUID);
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    ChildGUIs{end+1} = h;
    setappdata(handles.rtsMainGUI,'ChildGUIs',ChildGUIs);
end

% --- Executes on button press in rts_SCDynPageButton.
function rts_SCDynPageButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_SCDynPageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    vesselGUID = getappdata(handles.rtsMainGUI,'VesselGUID');
    h = rts_DynConsole(handles.rtsMainGUI, vesselGUID);
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    ChildGUIs{end+1} = h;
    setappdata(handles.rtsMainGUI,'ChildGUIs',ChildGUIs);
end

% --- Executes on button press in rts_PwrThrmlPageButton.
function rts_PwrThrmlPageButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_PwrThrmlPageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    vesselGUID = getappdata(handles.rtsMainGUI,'VesselGUID');
    h = rts_PwrThermConsole(handles.rtsMainGUI, vesselGUID);
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    ChildGUIs{end+1} = h;
    setappdata(handles.rtsMainGUI,'ChildGUIs',ChildGUIs);
end

% --- Executes on button press in rts_PropPageButton.
function rts_PropPageButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_PropPageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    vesselGUID = getappdata(handles.rtsMainGUI,'VesselGUID');
    h = rts_PropConsole(handles.rtsMainGUI, vesselGUID);
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    ChildGUIs{end+1} = h;
    setappdata(handles.rtsMainGUI,'ChildGUIs',ChildGUIs);
end

% --- Executes on button press in rts_ClockUtilButton.
function rts_ClockUtilButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_ClockUtilButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
end

% --- Executes when user attempts to close rtsMainGUI.
function rtsMainGUI_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to rtsMainGUI (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
    rts_EndTMStream(handles.rtsMainGUI, true, false, false);
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    for(i=1:length(ChildGUIs))
        if(ishandle(ChildGUIs{i}))
            delete(ChildGUIs{i});
        end
    end
    delete(hObject);
end

% --------------------------------------------------------------------
function connect2KSPTOTConnect_Callback(hObject, eventdata, handles)
% hObject    handle to connect2KSPTOTConnect (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    rts_StartTMStream(handles.rtsMainGUI);
end

% --------------------------------------------------------------------
function disconnectFromKSPTOTConnect_Callback(hObject, eventdata, handles)
% hObject    handle to disconnectFromKSPTOTConnect (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    rts_EndTMStream(handles.rtsMainGUI, false, false, false);
end

% --- Executes on button press in rts_TMMonUtilButton.
function rts_TMMonUtilButton_Callback(hObject, eventdata, handles)
% hObject    handle to rts_TMMonUtilButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    vesselGUID = getappdata(handles.rtsMainGUI,'VesselGUID');
    h = rts_TmMonitorConsole(handles.rtsMainGUI, vesselGUID);
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    ChildGUIs{end+1} = h;
    setappdata(handles.rtsMainGUI,'ChildGUIs',ChildGUIs);
end

% --------------------------------------------------------------------
function exitGUI_Callback(hObject, eventdata, handles)
% hObject    handle to exitGUI (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    ChildGUIs = getappdata(handles.rtsMainGUI,'ChildGUIs');
    for(i=1:length(ChildGUIs))
        if(ishandle(ChildGUIs{i}))
            close(ChildGUIs{i});
        end
    end
    close(handles.rtsMainGUI);
end

% --------------------------------------------------------------------
function connectionMenu_Callback(hObject, eventdata, handles)
% hObject    handle to connectionMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
end

% --------------------------------------------------------------------
function connect2NewVessel_Callback(hObject, eventdata, handles)
% hObject    handle to connect2NewVessel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    rts_EndTMStream(handles.rtsMainGUI, true, false, true);
    rts_StartTMStream(handles.rtsMainGUI);
end

% --------------------------------------------------------------------
function rtsOptions_Callback(hObject, eventdata, handles)
% hObject    handle to rtsOptions (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    rts_kspTotConnectOptions(handles.rtsMainGUI);
end

% --------------------------------------------------------------------
function terminateTMStream_Callback(hObject, eventdata, handles)
% hObject    handle to terminateTMStream (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    rts_EndTMStream(handles.rtsMainGUI, true, false, true);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%BUILDOPTICALFLOWPYRAMID  Constructs the image pyramid which can be passed to cv.calcOpticalFlowPyrLK
%
%     pyramid = cv.buildOpticalFlowPyramid(img)
%     [pyramid,maxLvl] = cv.buildOpticalFlowPyramid(img)
%     pyramid = cv.buildOpticalFlowPyramid(img, 'OptionName',optionValue, ...)
%
% ## Input
% * __img__ 8-bit input image.
%
% ## Output
% * __pyramid__ output pyramid. A cell-array of matrices (pyramid levels),
%   interleaved with corresponding gradients if `WithDerivatives` flag is
%   true. Each gradient is a matrix of type `int16` of same size as pyramid at
%   that level with 2-channels.
% * __maxLvl__ number of levels in constructed pyramid. Can be less than
%   `MaxLevel`.
%
% ## Options
% * __WinSize__ window size of optical flow algorithm. Must be not less than
%   `WinSize` argument of cv.calcOpticalFlowPyrLK. It is needed to calculate
%   required padding for pyramid levels. default [21,21]
% * __MaxLevel__ 0-based maximal pyramid level number. default 3
% * __WithDerivatives__ set to precompute gradients for the every pyramid
%   level. If pyramid is constructed without the gradients then
%   cv.calcOpticalFlowPyrLK will calculate them internally. default true
% * __PyrBorder__ the border mode for pyramid layers. See cv.copyMakeBorder.
%   default 'Reflect101'
% * __DerivBorder__ the border mode for gradients. See cv.copyMakeBorder.
%   default 'Constant'
% * __TryReuseInputImage__ put ROI of input image into the pyramid if
%   possible. You can pass false to force data copying. default true
%
% See also: cv.calcOpticalFlowPyrLK, cv.buildPyramid
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function x = colvec(x)
% Reshape the input into a column vector
% This function is useful when a function returns a row vector or matrix
% and you want to immediately reshape it in a functional way. Suppose
% f(a,b,c) returns a row vector, matlab will not let you write f(a,b,c)(:)
% - you would have to first store the result. With this function you can
% write colvec(f(a,b,c)) and be assured that the result is a column vector.

% This file is from pmtk3.googlecode.com

x = x(:);
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function varargout = CompressionPreferences(varargin)
% COMPRESSIONPREFERENCES MATLAB code for CompressionPreferences.fig
%      COMPRESSIONPREFERENCES, by itself, creates a new COMPRESSIONPREFERENCES or raises the existing
%      singleton*.
%
%      H = COMPRESSIONPREFERENCES returns the handle to a new COMPRESSIONPREFERENCES or the handle to
%      the existing singleton*.
%
%      COMPRESSIONPREFERENCES('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in COMPRESSIONPREFERENCES.M with the given input arguments.
%
%      COMPRESSIONPREFERENCES('Property','Value',...) creates a new COMPRESSIONPREFERENCES or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before CompressionPreferences_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to CompressionPreferences_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help CompressionPreferences

% Last Modified by GUIDE v2.5 24-Nov-2011 13:22:56

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @CompressionPreferences_OpeningFcn, ...
                   'gui_OutputFcn',  @CompressionPreferences_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% CompressionPreferences(jlabelGuiData)
function CompressionPreferences_OpeningFcn(hObject, eventdata, handles, varargin) %#ok<*INUSL>
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to CompressionPreferences (see VARARGIN)

assert(isa(varargin{1},'JLabelGUIData'));
handles.guidata = varargin{1};

handles = Initialize(handles);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes CompressionPreferences wait for user response (see UIRESUME)
uiwait(handles.figure_CompressionPreferences);

function handles = Initialize(handles)

% set last values
% useVideoWriter
% outavi_compression -> compression_format: profile (VideoWriter) or
% compression (avifile)
% outavi_fps -> fps
% outavi_quality -> quality: Quality (VideoWriter + Motion JPEG AVI) or
% Compression Ratio and LosslessCompression (VideoWriter + !Motion JPEG AVI
% or quality (avifile)

jlgd = handles.guidata;
if exist('VideoWriter','file')>0
  handles.useVideoWriter = jlgd.useVideoWriter;
else
  handles.useVideoWriter = true;
end
handles.compression_format = jlgd.outavi_compression;
handles.fps = jlgd.outavi_fps;
handles.quality = jlgd.outavi_quality;

handles = SetCompressionFormats(handles);

handles = UpdateGUIPrompts(handles);

set(handles.radiobutton_useVideoWriter,'Value',handles.useVideoWriter);
set(handles.radiobutton_useAviFile,'Value',~handles.useVideoWriter);
if ~jlgd.mat_lt_8p4,
  set(handles.radiobutton_useVideoWriter,'enable','off')
  set(handles.radiobutton_useAviFile,'enable','off');
end
%i = find(strcmpi(handles.compression_format,handles.compression_formats),1);
%set(handles.popupmenu_compression_profile,'Value',i);
set(handles.edit_fps,'String',num2str(handles.fps));
set(handles.edit_quality,'String',num2str(handles.quality));

function handles = UpdateGUIPrompts(handles)

% useVideoWriter
% outavi_compression -> compression_format: profile (VideoWriter) or
% compression (avifile)
% outavi_fps -> fps
% outavi_quality -> quality: Quality (VideoWriter + Motion JPEG AVI) or
% Compression Ratio and LosslessCompression (VideoWriter + ~Motion JPEG AVI
% or quality (avifile)

if ~exist('VideoWriter','file'),
  if handles.useVideoWriter,
    warning('This should never happen: useVideoWriter == true and VideoWriter does not exist');
    handles.useVideoWriter = false;
  end
  set(handles.uipanel_advanced,'Visible','off');
else
  set(handles.uipanel_advanced,'Visible','on');
end

if handles.useVideoWriter && strcmpi(handles.compression_format,'Motion JPEG AVI'),
  set(handles.text_quality,'String','Compression Quality');
else
  set(handles.text_quality,'String','Compression Ratio');
end

if (handles.useVideoWriter && ...
    ismember(lower(handles.compression_format),{'archival','none'})) || ...
    (~handles.useVideoWriter && strcmpi(handles.compression_format,'None')),
  set(handles.text_quality,'Enable','off');
  set(handles.edit_quality,'Enable','off');
else
  set(handles.text_quality,'Enable','on');
  set(handles.edit_quality,'Enable','on');
end

function handles = SetCompressionFormats(handles)

if handles.useVideoWriter,
  profiles = VideoWriter.getProfiles();
  handles.compression_formats = {profiles.Name};
  isvalid = profiles.isvalid;
  if ~all(isvalid),
    handles.compresion_formats(~isvalid) = [];
  end
  i = find(strcmpi(handles.compression_formats,'Uncompressed AVI'),1);
  if ~isempty(i),
    handles.compression_formats{i} = 'None';
  end
else
  if ispc,
    compressors = mmcompinfo;
    handles.compression_formats = [{'None'},{compressors.videoCompressors.Name}];
  else
    handles.compression_formats = {'None'};
  end
end

i = find(strcmpi(handles.compression_formats,handles.compression_format),1);
if isempty(i),
  i = 1;
  handles.compression_format = handles.compression_formats{1};
end

set(handles.popupmenu_compression_profile,'String',handles.compression_formats,'Value',i);

% --- Outputs from this function are returned to the command line.
function varargout = CompressionPreferences_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

varargout = {};
%delete(handles.figure_CompressionPreferences);

% --- Executes on button press in radiobutton_useVideoWriter.
function radiobutton_useVideoWriter_Callback(hObject, eventdata, handles) %#ok<*DEFNU>
% hObject    handle to radiobutton_useVideoWriter (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_useVideoWriter
handles.useVideoWriter = get(handles.radiobutton_useVideoWriter,'Value');
set(handles.radiobutton_useAviFile,'Value',~handles.useVideoWriter);
% translatable compression format
handles.compression_format = 'none';
handles.quality = 100;
handles = SetCompressionFormats(handles);
handles = UpdateGUIPrompts(handles);
guidata(hObject,handles);

% --- Executes on button press in radiobutton_useAviFile.
function radiobutton_useAviFile_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton_useAviFile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton_useAviFile
handles.useVideoWriter = get(handles.radiobutton_useAviFile,'Value') == 0;
set(handles.radiobutton_useVideoWriter,'Value',handles.useVideoWriter);
% translatable compression format
handles.compression_format = 'none';
handles.quality = 100;
handles = SetCompressionFormats(handles);
handles = UpdateGUIPrompts(handles);
guidata(hObject,handles);

% --- Executes on selection change in popupmenu_compression_profile.
function popupmenu_compression_profile_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu_compression_profile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu_compression_profile contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu_compression_profile
old_compression_format = handles.compression_format;
handles.compression_format = handles.compression_formats{get(handles.popupmenu_compression_profile,'Value')};
% translate quality
if handles.useVideoWriter && ...
    ismember(lower(handles.compression_format),{'archival','none'}),
  handles.quality = 1;
  set(handles.edit_quality,'String',num2str(handles.quality));
elseif ~handles.useVideoWriter && strcmpi(handles.compression_format,'None'),
  handles.quality = 100;
  set(handles.edit_quality,'String',num2str(handles.quality));
elseif strcmpi(old_compression_format,'Motion JPEG AVI') && ...
    ~strcmpi(handles.compression_format,'Motion JPEG AVI'),
  % >99 -> 1, <=99 -> 10
  if handles.quality >= 99,
    handles.quality = 1;
  else
    handles.quality = 10;
  end
  set(handles.edit_quality,'String',num2str(handles.quality));
elseif ~strcmpi(old_compression_format,'Motion JPEG AVI') && ...
    strcmpi(handles.compression_format,'Motion JPEG AVI'),
  % <= 1.5 -> 100, >1.5 -> 75
  if handles.quality <= 1.5,
    handles.quality = 100;
  else
    handles.quality = 75;
  end
  set(handles.edit_quality,'String',num2str(handles.quality));
end

handles = UpdateGUIPrompts(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function popupmenu_compression_profile_CreateFcn(hObject, eventdata, handles) %#ok<*INUSD>
% hObject    handle to popupmenu_compression_profile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function edit_quality_Callback(hObject, eventdata, handles)
% hObject    handle to edit_quality (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_quality as text
%        str2double(get(hObject,'String')) returns contents of edit_quality as a double

v = str2double(get(hObject,'String'));
if isnan(v),
  uiwait(warndlg(sprintf('%s must be a number',get(handles.text_quality,'String'))));
  set(hObject,'String',num2str(handles.quality));
  return;
end
if handles.useVideoWriter && ~strcmpi(handles.compression_format,'Motion JPEG AVI'),
  if v < 1,
    uiwait(warndlg('Compression ratio should be >= 1'));
    set(hObject,'String',num2str(handles.quality));
    return;
  end
else
  if v <= 0 || v > 100,
    uiwait(warndlg('Compression quality should be > 0 and <= 100'));
    set(hObject,'String',num2str(handles.quality));
    return;
  end
end
handles.quality = v;
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_quality_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_quality (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton_done.
function pushbutton_done_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_done (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% store current configuration
jlgd = handles.guidata;
jlgd.outavi_fps = handles.fps;
jlgd.outavi_compression = handles.compression_format;
jlgd.useVideoWriter = handles.useVideoWriter;
jlgd.outavi_quality = handles.quality;

delete(handles.figure_CompressionPreferences);

% --- Executes during object creation, after setting all properties.
function pushbutton_done_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pushbutton_done (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --- Executes on button press in pushbutton_cancel.
function pushbutton_cancel_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

delete(handles.figure_CompressionPreferences);

function edit_fps_Callback(hObject, eventdata, handles)
% hObject    handle to edit_fps (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_fps as text
%        str2double(get(hObject,'String')) returns contents of edit_fps as a double

v = str2double(get(hObject,'String'));
if isnan(v) || v <= 0 || round(v) ~= v,
  uiwait(warndlg('Frames per second must be a whole number greater than zero'));
  set(hObject,'String',num2str(handles.fps));
  return;
end
handles.fps = v;
guidata(hObject,handles);

author: kristinbranson

repo: JAABA


text:
function combo_test17 ( )

%*****************************************************************************80
%
%% COMBO_TEST17 tests MOUNTAIN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST17\n' );
  fprintf ( 1, '  MOUNTAIN computes mountain numbers.\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '   Y    MXY\n' );
  fprintf ( 1, ' \n' );

  for y = 0 : n
    fprintf ( 1, '  %2d   ', y );
    for x = 0 : 2 * n
      mxy = mountain ( n, x, y );
      fprintf ( 1, '%4d', mxy );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

%% Syntax
% |[E,VE,Fq1,Fq2]=patchThick(F,V,dirSet,layerThickness,numSteps);|

%% Description
% Use |patchThick| to thicken a quadrilateral mesh to create hexahedral
% elements.

%%

clear; close all; clc;

%% Examples

%%
% PLOT SETTINGS
fontSize=20;

%% Example: Using |patchThick| to thicken quads into hexahedra

%%
% Creating an example polygon
ns=15;
t=linspace(0,pi,ns);
x=cos(t);
y=sin(t);
z=zeros(size(x));
Vc=flipud([x(:) y(:) z(:)]);

%%
% Extruding polygon to a quadrilateral surface
cPar.depth=2;
cPar.patchType='quad';
cPar.dir=0;
cPar.closeLoopOpt=0;
cPar.numSteps=8;
[F,V]=polyExtrude(Vc,cPar);

%%
% Thickening quadrilaterial elements to hexahedral elements
layerThickness=0.5;
numSteps=3;
[E,VE]=patchThick(F,V,1,layerThickness,numSteps);

%Use element2patch to get patch data
FE=element2patch(E);

%%
% Visualize mesh

cFigure; hold on;
title('Hexahedral mesh');
gpatch(FE,VE,'bw','k',1);
axisGeom;
camlight headlight;
drawnow;

%% Example: Using |patchThick| to thicken triangles into pentahedra

%%
% Create test data set
[F,V,C]=hemiSphereMesh(1,1,0);
V(:,3)=-V(:,3);
F=fliplr(F);

%%
% Thickening quadrilaterial surface to hexahedral elements
layerThickness=0.4;
numSteps=2;
[E,VE,Fq1,Fq2]=patchThick(F,V,1,layerThickness,numSteps);

%Use element2patch to get patch data
FE=element2patch(E,[],'penta6');

%%
% Visualize mesh

cFigure; hold on;
title('Pentahedral mesh');
gpatch(FE,VE,'bw','k',1);
% patchNormPlot(FE,VE);
axisGeom;
camlight headlight;
drawnow;

%% Example: Using |patchThick| to thicken a mixed mesh of quads and triangles

%%
% Create test data set
[F,V,C]=hemiSphereMesh(1,1,0);
V(:,3)=-V(:,3);
F=fliplr(F);

%%
% Convert to mixed mesh
optionStruct.maxAngleDeviation=90*(pi/180);
optionStruct.selectionMethod='best';
optionStruct.triangleConvert=0;
optionStruct.fourConnectConvert=0;
[F,V]=tri2quadGroupSplit(F,V,optionStruct);

%%
% Thickening quadrilaterial surface to hexahedral elements
layerThickness=0.4;
numSteps=2;
[E,VE,Fq1,Fq2]=patchThick(F,V,1,layerThickness,numSteps);

%Use element2patch to get patch data
FE=E;
for q=1:1:numel(E)
    if size(E{q},2)==6
        elementType='penta6';
    elseif size(E{q},2)==8
        elementType='hex8';
    end
    FE{q}=element2patch(E{q},[],elementType);
end

%%
% Visualize mesh

cFigure; hold on;
title('Mixed hexahedral-pentahedral mesh');
for q=1:1:numel(FE)
    gpatch(FE{q},VE,'bw','k',1);
end
% patchNormPlot(FE,VE);
axisGeom;
camlight headlight;
drawnow;

%%
%
% <<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:
% DESCRIPTION:
%       subscript to scale source terms to the correct units
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 15th February 2012
%       last update - 20th January 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% get the dimension size
N = kgrid.dim;

% check for non-uniform grid and give error for source terms that haven't
% yet been implemented
if (nonuniform_grid && ( uy_source || uz_source || transducer_source))
    disp('WARNING: source scaling not implemented for non-uniform grids with given source condition');
    return
end

% Scale the input pressure by 1/c^2 (to convert to units of density), then
% by 1/N (to split the input across the split density field). If the
% pressure is injected as a mass source, also scale the pressure by
% 2*dt*c/dx to account for the time step and convert to units of
% [kg/(m^3 s)]
if p_source
    if strcmp(source.p_mode, 'dirichlet')
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous
            % sound speed
            source.p = source.p ./ (N*c^2);
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for p_index = 1:length(source.p(:,1))
                source.p(p_index, :) = source.p(p_index, :) ./ (N*c(p_source_index(p_index))^2);
            end
        end
    else
        if nonuniform_grid

            % create empty matrix
            grid_point_sep = zeros(size(kgrid.x));

            % compute averaged grid point seperation map, the interior
            % points are calculated using the average distance to all
            % connected grid points (the edge values are not calculated
            % assuming there are no source points in the PML)
            switch kgrid.dim
                case 1
                    grid_point_sep(2:end-1) = kgrid.x_size*(kgrid.xn(3:end, 1) - kgrid.xn(1:end-2, 1))/2;
                case 2
                    grid_point_sep(2:end-1, 2:end-1) = ...
                        kgrid.x_size*(kgrid.xn(3:end, 2:end-1) - kgrid.xn(1:end-2, 2:end-1))/4 + ...
                        kgrid.y_size*(kgrid.yn(2:end-1, 3:end) - kgrid.yn(2:end-1, 1:end-2))/4;
                case 3
                    grid_point_sep(2:end-1, 2:end-1, 2:end-1) = ...
                        kgrid.x_size*(kgrid.xn(3:end, 2:end-1, 2:end-1) - kgrid.xn(1:end-2, 2:end-1, 2:end-1))/6 + ...
                        kgrid.y_size*(kgrid.yn(2:end-1, 3:end, 2:end-1) - kgrid.yn(2:end-1, 1:end-2, 2:end-1))/6 + ...
                        kgrid.z_size*(kgrid.zn(2:end-1, 2:end-1, 3:end) - kgrid.zn(2:end-1, 2:end-1, 1:end-2))/6;
            end

            % compute and apply scale parameter
            for p_index = 1:length(source.p(:,1))
                if numel(c) == 1
                    % compute the scale parameter based on the homogeneous sound speed
                    source.p(p_index, :) = source.p(p_index, :) .* (2.*dt./(N*c*grid_point_sep(p_source_index(p_index))));
                else
                    % compute the scale parameter based on the sound speed at that position
                    source.p(p_index, :) = source.p(p_index, :) .* (2.*dt./(N*c(p_source_index(p_index)).*grid_point_sep(p_source_index(p_index))));
                end
            end

            % clear unused variables
            clear grid_point_sep;

        else
            if numel(c) == 1
                % compute the scale parameter based on the homogeneous
                % sound speed
                source.p = source.p .* (2*dt./(N*c*kgrid.dx));
            else
                % compute the scale parameter seperately for each source
                % position based on the sound speed at that position
                for p_index = 1:length(source.p(:,1))
                    source.p(p_index, :) = source.p(p_index, :) .* (2.*dt./(N*c(p_source_index(p_index)).*kgrid.dx));
                end
            end
        end
    end
end

% scale the stress source by 1/N to divide amoungst the split field
% components, and if source.s_mode is not set to 'dirichlet', also scale by
% 2*dt*c/dx to account for the time step and convert to units of
% [kg/(m^3 s)]
if sxx_source
    if strcmp(source.s_mode, 'dirichlet') || p0_source
        source.sxx = source.sxx ./ N;
    else
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound
            % speed
            source.sxx = source.sxx .* (2*dt*c./(N*kgrid.dx));
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for s_index = 1:length(source.sxx(:,1))
                source.sxx(s_index, :) = source.sxx(s_index, :) .* (2.*dt.*c(s_source_index(s_index))./(N*kgrid.dx));
            end
        end
    end
end
if syy_source
    if strcmp(source.s_mode, 'dirichlet') || p0_source
        source.syy = source.syy ./ N;
    else
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound
            % speed
            source.syy = source.syy .* (2*dt*c./(N*kgrid.dx));
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for s_index = 1:length(source.syy(:,1))
                source.syy(s_index, :) = source.syy(s_index, :) .* (2.*dt.*c(s_source_index(s_index))./(N*kgrid.dx));
            end
        end
    end
end
if szz_source
    if strcmp(source.s_mode, 'dirichlet') || p0_source
        source.szz = source.szz ./ N;
    else
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound
            % speed
            source.szz = source.szz .* (2*dt*c./(N*kgrid.dx));
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for s_index = 1:length(source.szz(:,1))
                source.szz(s_index, :) = source.szz(s_index, :) .* (2.*dt.*c(s_source_index(s_index))./(N*kgrid.dx));
            end
        end
    end
end
if sxy_source
    if strcmp(source.s_mode, 'dirichlet')
        source.sxy = source.sxy ./ N;
    else
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound
            % speed
            source.sxy = source.sxy .* (2*dt*c./(N*kgrid.dx));
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for s_index = 1:length(source.sxy(:,1))
                source.sxy(s_index, :) = source.sxy(s_index, :) .* (2.*dt.*c(s_source_index(s_index))./(N*kgrid.dx));
            end
        end
    end
end
if sxz_source
    if strcmp(source.s_mode, 'dirichlet')
        source.sxz = source.sxz ./ N;
    else
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound
            % speed
            source.sxz = source.sxz .* (2*dt*c./(N*kgrid.dx));
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for s_index = 1:length(source.sxz(:,1))
                source.sxz(s_index, :) = source.sxz(s_index, :) .* (2.*dt.*c(s_source_index(s_index))./(N*kgrid.dx));
            end
        end
    end
end
if syz_source
    if strcmp(source.s_mode, 'dirichlet')
        source.syz = source.syz ./ N;
    else
        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound
            % speed
            source.syz = source.syz .* (2*dt*c./(N*kgrid.dx));
        else
            % compute the scale parameter seperately for each source
            % position based on the sound speed at that position
            for s_index = 1:length(source.syz(:,1))
                source.syz(s_index, :) = source.syz(s_index, :) .* (2.*dt.*c(s_source_index(s_index))./(N*kgrid.dx));
            end
        end
    end
end

% if source.u_mode is not set to 'dirichlet', scale the x-direction
% velocity source terms by 2*dt*c/dx to account for the time step and
% convert to units of [m/s^2]
if ux_source && ~strcmp(source.u_mode, 'dirichlet')

    if nonuniform_grid

        % create empty matrix
        grid_point_sep = zeros(size(kgrid.x));

        % compute averaged grid point seperation map, the interior
        % points are calculated using the average distance to all
        % connected grid points (the edge values are not calculated
        % assuming there are no source points in the PML)
        grid_point_sep(2:end-1, :, :) = kgrid.x_size*(kgrid.xn(3:end, :, :) - kgrid.xn(1:end-2, :, :))/2;

        % compute and apply scale parameter
        for u_index = 1:length(source.ux(:,1))
            if numel(c) == 1
                % compute the scale parameter based on the homogeneous sound speed
                source.ux(u_index, :) = source.ux(u_index, :) .* (2*c*dt./(grid_point_sep(u_source_index(u_index))));
            else
                % compute the scale parameter based on the sound speed at that position
                source.ux(u_index, :) = source.ux(u_index, :) .* (2*c(u_source_index(u_index))*dt./(grid_point_sep(u_source_index(u_index))));
            end
        end

        % clear unused variables
        clear grid_point_sep;

    else

        if numel(c) == 1
            % compute the scale parameter based on the homogeneous sound speed
            source.ux = source.ux .* (2*c*dt./kgrid.dx);
        else
            % compute the scale parameter seperately for each source position
            % based on the sound speed at that position
            for u_index = 1:length(source.ux(:,1))
                source.ux(u_index, :) = source.ux(u_index, :) .* (2.*c(u_source_index(u_index)).*dt./kgrid.dx);
            end
        end
    end
end

% if source.u_mode is not set to 'dirichlet', scale the y-direction
% velocity source terms by 2*dt*c/dy to account for the time step and
% convert to units of [m/s^2]
if uy_source && ~strcmp(source.u_mode, 'dirichlet')
    if numel(c) == 1
        % compute the scale parameter based on the homogeneous sound speed
        source.uy = source.uy .* (2*c*dt./kgrid.dy);
    else
        % compute the scale parameter seperately for each source position
        % based on the sound speed at that position
        for u_index = 1:length(source.uy(:,1))
            source.uy(u_index, :) = source.uy(u_index, :) .* (2.*c(u_source_index(u_index)).*dt./kgrid.dy);
        end
    end
end

% if source.u_mode is not set to 'dirichlet', scale the z-direction
% velocity source terms by 2*dt*c/dz to account for the time step and
% convert to units of [m/s^2]
if uz_source && ~strcmp(source.u_mode, 'dirichlet')
    if numel(c) == 1
        % compute the scale parameter based on the homogeneous sound speed
        source.uz = source.uz .* (2*c*dt./kgrid.dz);
    else
        % compute the scale parameter seperately for each source position
        % based on the sound speed at that position
        for u_index = 1:length(source.uz(:,1))
            source.uz(u_index, :) = source.uz(u_index, :) .* (2.*c(u_source_index(u_index)).*dt./kgrid.dz);
        end
    end
end

% scale the transducer source term by 2*dt*c/dx to account for the time
% step and convert to units of [m/s^2]
if transducer_source
    if numel(c) == 1
        transducer_input_signal = transducer_input_signal .* (2*c*dt./kgrid.dx);
    else
        % compute the scale parameter based on the average sound speed at the
        % transducer positions (only one input signal is used to drive the
        % transducer)
        transducer_input_signal = transducer_input_signal.* (2*(mean(c(u_source_index)))*dt./kgrid.dx);
    end
end

% clear subscript variables
clear N;

repo: Aero-matlab


text:
%% 1D Heat Transfer LBM D1Q2
% by Manuel Diaz
clc; clear; %close all;

%% Domain Grid
L = 30;    % Length
n = 30;    % Nodes
dx = L/n; dt = 1;
x  = 1:dx:L;

%% Initial Condition
f1 = zeros(1,n);
f2 = zeros(1,n);
rho= zeros(1,n); % initial value of the dependent variable rho = T(x,t)
feq1= zeros(1,n);
feq2= zeros(1,n);

%% Constants
csq   = (dx^2)/(dt^2);
alpha = 0.25;
omega = 1/(alpha/(dt*csq)+0.5);
w     = [1/2,1/2];
cx    = [  1, -1];
cy    = [  0,  0];
links = [  1,  2];
tEnd  = 200; % time steps
tPlot = 10;
frames= tEnd/tPlot;
count = 0;

%% Movie Parameters
tPlot = 10; frames= tEnd/tPlot; count = 0;
figure(1)
colordef white

%% Boundary Conditions
twall = 1;

%% Main Loop
for cycle = 1:tEnd;
    % Collision process
    rho = f1+f2;

    % even that for this case k1=k2=1/2m then feq1=feq2=feq
    feq1 = 0.5*rho;
    feq2 = 0.5*rho;

    f1 = (1-omega) * f1 + omega * feq1;
    f2 = (1-omega) * f2 + omega * feq2;
    % Streaming process
    f1 = stream2d( f1 , [cy(1),cx(1)]);
    f2 = stream2d( f2 , [cy(2),cx(2)]);

%   % Original streaming process in A.A.Mohamad
%     for i = 2:n-1;  % Double streaming in one loop
%         f1(n-i+1) = f1(n-i);  % Streaming
%         f2( i-1 ) = f2( i );  % Streaming
%     end

    % Boundary condition
    f1(1) = twall-f2(1);    %Dirichlet BC
    f1(n) = f1(n-1);        %Neumann BC
    f2(n) = f2(n-1);        %Neumann BC

    % Visualization every tPlot
    if mod(cycle,tPlot) == 1
        count = count + 1;
        plot(x,rho,'.');
        title '1D Heat Equation using LBM D1Q2'
        xlabel 'x cell'
        ylabel 'Temperature'
        M(count)=getframe;
    end

    % Animated gif file
    if mod(cycle,tPlot) == 1
        F = getframe;
        if count == 1
            [im,map] = rgb2ind(F.cdata,256,'nodither');
            im(1,1,1,tEnd/tPlot) = 0;
        end
        im(:,:,1,count) = rgb2ind(F.cdata,map,'nodither');
    end
end

%% Make pretty figures
%plot(x,rho,'.');

%% Make Movie
movie(M,3,10); % movie(M,n,fps)

%% Export to Gif
imwrite(im,map,'heat_eq_d1q2.gif','DelayTime',0,'LoopCount',3)


text:
% SCRIPT TEST FOR THE KINEMATIC PROBLEM FOR SERIAL ROBOTS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
%
% qinv =
%
%          0         0         0         0    3.1416    3.1416    3.1416    3.1416
%    -0.0412   -0.0412    1.6597    1.6597   -1.5798   -1.5798   -0.3550   -0.3550
%     0.0720    0.0720    3.0696    3.0696   -0.3825   -0.3825   -2.7591   -2.7591
%     0.0000   -3.1416    3.1416   -0.0000    0.0000    3.1416    0.0000   -3.1416
%    -0.0308    0.0308   -1.5539    1.5539   -1.1793    1.1793   -0.0275    0.0275
%    -0.0000    3.1416   -3.1416         0    3.1416   -0.0000    3.1416   -0.0000

function test_kinematics_irb140

close all

fprintf('\nTHE DEMO PRESENTS THE DIRECT AND INVERSE KINEMATIC PROBLEM')

%load robot parameters. You can try different robots%
robot=load_robot('ABB', 'IRB140');
%adjust 3D view as desired
adjust_view(robot)

% reach this!
T =[0.0000   -0.0000    1.0000    0.4000;
   -0.0000    1.0000    0.0000   -0.0000;
   -1.0000   -0.0000    0.0000    0.6000;
   0         0         0    1.0000];

%Call the inversekinematic for this robot. All the possible solutions are
%stored at qinv. At least, one of the possible solutions should match q
qinv = inversekinematic(robot, T)

test_joints(robot, qinv)

test_solutions(robot, qinv, T);

test_movement(robot, qinv, T);

function test_movement(robot, qinv, T)

q0 = [0, 0, 0, 0, 0, 0];
n = 5
for i=1:size(qinv,2)
   qq = [];
   qi = qinv(:,i);
   for j=1:6
        v = linspace(q0(j), qi(j), n);
        qq = [qq; v];
   end
   animate(robot, qq)
end

function test_solutions(robot, qinv, T)
n_solutions = 8;

fprintf('\nNOW WE CAN REPRESENT THE DIFFERENT SOLUTIONS TO ACHIEVE THE SAME POSITION AND ORIENTATION\n')
fprintf('\nNote that some solutions may not be feasible since some joints may be out of range.\n')
correct=zeros(1,n_solutions);
%check that all of them are possible solutions!
for i=1:size(qinv,2)

    Ti = directkinematic(robot, qinv(:,i)) %Ti is constant for the different solutions

    % Note that all the solutions may not be feasible. Some of the joints may
    % be out of range. You can test this situation with test_joints
    test_joints(robot, qinv(:,i));

    %now draw the robot to see the solution
    drawrobot3d(robot, qinv(:,i))

    pause(0.5);

    k=sum(sum((T-Ti).^2));
    if k < 0.01 % a simple threshold to find differences in the solution
        correct(1,i)= 1;
    else
        correct(1,i)= 0; %uncorrect solution
        fprintf('\nERROR: One of the solutions seems to be uncorrect. Sum of errors: %f', i, k);
    end
end

fprintf('\n************** RESULTS **************')

%Display a message if any of the solutions is not correct
if sum(correct)==n_solutions
    fprintf('\nTEST 1--> OK: Every solution in qinv yields the same position/orientation T');
else
    fprintf('\nTEST 1--> ERROR: One or more of the solutions seem to be uncorrect.');
end

fprintf('\n************** ****** **************\n')

repo: ARTE


text:
function [league] = m2league(m)
% Convert length from meters to leagues.
% Chad A. Greene 2012
league = m*0.0002071237307458;


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%

% this wrapper around MatLab internal map is meant for
% polymorphic function signature compatability with other
% classdef and JAVA class in CytoGenie AutoGate
classdef BasicMap < Map
    properties
        pu;
        parentCmpForPopup;
        needToAskForUmap=true;
        noDbscan=false;
        urlMap;
    end

    properties(SetAccess=private)
        appFolder;
        contentFolder;
        toolBarSize=0;
        toolBarFactor=0;
        highDef=false;
        supStart='';
        supEnd='';
        subStart='';
        subEnd='';
        smallStart='';
        smallEnd='';
        h3Start='
';
        h3End='
';
        h2Start='
';
        h2End='
';
        h1Start='
';
        h1End='
';
        whereMsgOrQuestion='center';
    end

    methods(Static)

        function this=Global(closeIfTrueOrMap)
            persistent singleton;
            if nargin>0 && islogical(closeIfTrueOrMap) && closeIfTrueOrMap
                clear singleton;
                singleton=[];
                disp('Resetting global BasicMap');
                this=[];
            else
                if nargin==1
                    try
                        priorMap=closeIfTrueOrMap;
                        %test method compatibility
                        priorMap.size
                        prop='IsBasicMap';
                        if ~priorMap.has(prop)
                            priorMap.set(prop, 'false')
                            priorMap.get(prop, 'true')
                            priorMap.remove(prop);
                        else
                            priorMap.get(prop, 'true')
                        end
                        singleton=priorMap;
                        this=singleton;
                        return;
                    catch ex
                        ex.getReport
                    end
                end
                if isempty(singleton)
                    singleton=BasicMap;
                    singleton.highDef=Gui.hasPcHighDefinitionAProblem(2000, 2500, false);
                    BasicMap.SetHighDef(singleton, singleton.highDef);
                end
                this=singleton;
            end
        end

        function path=Path
            path=BasicMap.Global.contentFolder;
        end

        function obj=SetHighDef(obj, hasHighDef)
            factor=0;
            NORMAL_FONT_SIZE=12;
            SMALL_FONT_SIZE=2;
            H3_FONT_SIZE=3;
            H2_FONT_SIZE=3.5;
            H1_FONT_SIZE=4;
            if hasHighDef
                obj.highDef=true;
                factor=javax.swing.UIManager.getFont('Label.font').getSize...
                    /NORMAL_FONT_SIZE;
            else
                if ismac
                    %factor=1.6;
                end
            end
            if factor>0
                obj.toolBarFactor=factor;
                obj.toolBarSize=floor(16*factor);
                smallSize=floor(SMALL_FONT_SIZE*factor);
                if ispc
                    smallSize=smallSize+1;
                end
                obj.smallStart=[''];
                obj.smallEnd='';
                obj.subStart=obj.smallStart;
                obj.supStart=obj.smallStart;
                obj.subEnd=obj.smallEnd;
                obj.supEnd=obj.smallEnd;
                h1Size=floor(H1_FONT_SIZE *factor);
                if ispc
                    h1Size=h1Size+1;
                end
                obj.h1Start=[''];
                obj.h1End='
';

                h2Size=floor(H2_FONT_SIZE *factor);
                if ispc
                    h2Size=h2Size+1;
                end
                obj.h2Start=[''];
                obj.h2End='
';

                h3Size=floor(H3_FONT_SIZE *factor);
                if ispc
                    h3Size=h3Size+1;
                end
                obj.h3Start='
';
                obj.h3End='
';
            else
                obj.toolBarSize=0;
                obj.toolBarFactor=0;
                obj.highDef=false;

                obj.smallStart='';
                obj.smallEnd='';
                obj.h3Start='
';
                obj.h3End='
';
                obj.h2Start='
';
                obj.h2End='
';
                obj.h1Start='
';
                obj.h1End='
';

            end
        end

        function nums=GetNumbers(props, name)
            nums=props.get(name, []);
            if ~isempty(nums)
                nums=str2num(nums);
            else
                nums=[];
            end
        end

    end

    methods
        function this=BasicMap(keysOrFileName, values)
            if nargin<2
                values={};
                if nargin<1
                    keysOrFileName=[];
                end
            end
            this=this@Map(keysOrFileName, values);
            this.contentFolder=fileparts(mfilename('fullpath'));
            this.appFolder=fullfile(File.Home, '.run_umap');
            File.mkDir(this.appFolder);
            this.urlMap=Map;
            %'https://1drv.ms/u/s!AkbNI8Wap-7_jNJYg4RNDTKR4mkYOg?e=pfwWfO'
        end

        function closeToolTip(this)
        end

        function retry=reportProblem(this, exception, retryId)
            retry=false;
            exception.getReport;
        end

        function [name, found]=getMatchName(this, name, showEppImg, b1, b2)
            if nargin<4
                if this.highDef
                    b1='';
                    b2='';
                else
                    b1='';
                    b2='';
                end
                if nargin<3
                    if nargin<2
                        showEppImg=true;
                    end
                end
            end
            idx=strfind(name, ',\bf');
            if ~isempty(idx)
                sIdx=4;
                idx=idx(1);
            else
                idx=strfind(name, '\bf');
                sIdx=3;
            end
            if ~isempty(idx)
                matchName=strtrim(name(idx(1)+sIdx:end));
                eppName=name(1:idx(1)-1);
                if showEppImg
                    eppName =[ eppName ' ' this.eppImg];
                end
                if idx==1
                    name=[matchName eppName];
                else
                    name=[matchName ' ' b1 eppName b2];
                end
                found=true;
            else
                found=false;
            end
        end

    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function test_failed=test_dwilt2

test_failed=0;

disp(' ===============  TEST_DWILT2 ================');

% Run some fixed test to test the interface.
% This is not a thourough tester.

% --- test 1 ----------------

L=64;
M=8;
Lf=63;
W=3;

f=tester_rand(Lf,Lf,W);

g=pgauss(L,1);
gd=wildual(g,M);

[c,Ls]=dwilt2(f,g,M);
r=idwilt2(c,gd,Ls);

res=r-f;

nres=norm(res(:));

[test_failed,fail]=ltfatdiditfail(nres,test_failed);
% failed='';
% if nres>10e-10
  % failed='FAILED';
  % test_failed=test_failed+1;
% end;

s=sprintf('DWILT2 Lf:%3i L:%3i %0.5g %s',Lf,L,nres,fail);
disp(s)

% --- test 2 -------------------
L=256;
M1=16;
M2=32;
W=1;

f=tester_rand(L,L,1);

g=pgauss(L,1);

gd1=wildual(g,M1);
gd2=wildual(g,M2);

c=dwilt2(f,g,[M1,M2]);
c2=ref_dwilt2(f,g,g,M1,M2);

rc=c-c2;
nres=norm(rc(:));
[test_failed,fail]=ltfatdiditfail(nres,test_failed);
% failed='';
% if nres>10e-10
  % failed='FAILED';
  % test_failed=test_failed+1;
% end;

s=sprintf('DWILT2 REF M1:%3i M2:%3i %0.5g %s',M1,M2,nres,fail);
disp(s)

r=idwilt2(c,gd1,gd2);

res=r-f;

nres=norm(res(:));
[test_failed,fail]=ltfatdiditfail(nres,test_failed);
% failed='';
% if nres>10e-10
  % failed='FAILED';
  % test_failed=test_failed+1;
% end;

s=sprintf('DWILT2 INV M1:%3i M2:%3i %0.5g %s',M1,M2,nres,fail);
disp(s)

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function K = matching( k, d1, d2, ind1, ind2, par)

%
%  Matching type similarity function (actually not a kernel) as described in the paper:
%  "Recognition with local features: the kernel recipe"
%  by C. Wallraven, B.Caputo and A.Graf
%
%  takes 2 parameters: { minorKernelName, minorKernelParameter}
%
%  implemented minor kernel names are:
%   'linear', 'rbf' ....  (standard kernels)
%   'normCC' ...........  normalized cross-correlation kernel
%   'spider' ...........  uses a spider kernel object (passed in minorKernelParameter),
%                         (is slower due to object-handling overhead)
%

  global VERBOSITY

  x1 = get_x( d1, ind1);
  x2 = get_x( d2, ind2);
  sz1 = size( x1, 1);
  sz2 = size( x2, 1);
  K = zeros( sz2, sz1);
  nTot = sz1*sz2;

  for ii = 1:sz1
    for jj = 1:sz2
      K( jj, ii) = Fkernel( double( x1{ ii}), double( x2{ jj}), par);
    end
    if VERBOSITY Fprogress( (sz2*(ii-1)+jj)/nTot); end
  end

  K = 0.5*( K + K');

function k = Fkernel( a, b, par)

  switch par{ 1}
   case 'linear'
    M = b*a';

   case 'rbf'
    M = rbf( a, b, par{ 2});

   case 'normCC'
    M = normCC( a, b, par{ 2});

   case 'spider'
    M = calc( par{ 2}, data( a), data( b));

   otherwise
    error( 'Kernel name unknown or not implemented: %s', par{ 1});
  end

  k = mean( max( M));

function K = rbf( a, b, sigma)
% exponential kernel
  sigma2 = sigma * sigma;
  K = b*a';
  K = K + K; % *2
  K1 = sum( a.^2, 2);
  K2 = sum( b.^2, 2);
  K = ones( length( K2), 1)*K1' + K2*ones( 1, length( K1)) - K;
  K = exp( -K ./ ( 2*sigma2));

function K = normCC( dat1, dat2, rho)
% normalized cross-correlation

  dat1 = dat1 - repmat( mean( dat1), size( dat1, 1), 1);
  dat2 = dat2 - repmat( mean( dat2), size( dat2, 1), 1);
  dat1 = dat1 ./ repmat( sqrt( diag( dat1*dat1')), 1, size( dat1, 2));
  dat2 = dat2 ./ repmat( sqrt( diag( dat2*dat2')), 1, size( dat2, 2));

  K = exp( -rho*( 1 - dat1*dat2'));


text:
function exact = p08_exact ( )

%*****************************************************************************80
%
%% P08_EXACT returns the estimated integral for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 0.86697298733991103757;

  return
end


text:
function [p, npix] = histroi(f, c, r)
%HISTROI Computes the histogram of an ROI in an image.
%   [P, NPIX] = HISTROI(F, C, R) computes the histogram, P, of a
%   polygonal region of interest (ROI) in image F.  The polygonal
%   region is defined by the column and row coordinates of its
%   vertices, which are specified (sequentially) in vectors C and R,
%   respectively. All pixels of F must be >= 0. Parameter NPIX is the
%   number of pixels in the polygonal region.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/09/05 16:14:35 $

% Generate the binary mask image.
B = roipoly(f, c, r);

% Compute the histogram of the pixels in the ROI.
p = imhist(f(B));

% Obtain the number of pixels in the ROI if requested in the output.
if nargout > 1
   npix = sum(B(:));
end

repo: weiyanmin


text:
function adv = adversarial_perturbation(x,l,Df_base,f_out,opts)
NUM_LABELS = 10;
OS = 0.02;
Q = 2;
MAX_ITER = 100;
if(nargin==5)
    if isfield(opts,'labels_limit') NUM_LABELS = opts.labels_limit;end;
    if isfield(opts,'overshoot') OS = opts.overshoot;end;
    if isfield(opts,'norm_p')
        Q = opts.norm_p/(opts.norm_p-1);
        if opts.norm_p==Inf
            Q = 1;
        end
    end
    if isfield(opts,'max_iter') MAX_ITER = opts.max_iter;end;
end

Df = @(y,idx) Df_base(y,l,idx);

ff = f_out(x,0);
ff = ff-ff(l);
[~,I] = sort(ff,'descend');
labels = I(2:NUM_LABELS);

r = x*0;
x_u = x;

itr = 0;
while(f_out(x+(1+OS)*r,1)==l && itr<MAX_ITER)
    itr = itr + 1;

    ff = f_out(x_u,0);
    ff = ff-ff(l);

    idx = [l labels];
    ddf = Df(x_u,idx);

    dr = project_boundary_polyhedron(ddf,ff(idx),Q);

    x_u = x_u+dr;
    r = r + dr;
end

adv.r = (1+OS)*r;
adv.new_label = f_out(x+(1+OS)*r,1);
adv.itr = itr;
end

function dir = project_boundary_polyhedron(Df,f,Q)
res = abs(f)./arrayfun(@(idx) norm(Df(:,idx),Q), 1:size(Df,2));
[~,ii]=min(res);
if isinf(Q)
    dir = res(ii).*(abs(Df(:,ii))>=max(abs(Df(:,ii)))).*sign(Df(:,ii));
elseif(Q==1)
    dir = res(ii).*sign(Df(:,ii));
else
    dir = res(ii)*(abs(Df(:,ii))/norm(Df(:,ii),Q)).^(Q-1).*sign(Df(:,ii));
end
end

repo: DeepFool


text:
  classdef viewPC < matlab.System & matlab.system.mixin.CustomIcon
    % Untitled Add summary here
    %
    % This template includes the minimum set of functions required
    % to define a System object with discrete state.

    % Public, tunable properties
    properties

    end

    properties(DiscreteState)
        xyz;
        p1;
        p2;
        p3;
    end

    % Pre-computed constants
    properties(Access = private)

    end

    methods(Access = protected)
        function stepImpl(obj,xyzpoint)
            % Perform one-time calculations, such as computing constants
            obj.p1 = xyzpoint(:,1);
            obj.p2 = xyzpoint(:,2);
            obj.p3 = xyzpoint(:,3);

            obj.xyz = [obj.p1,obj.p2,obj.p3];
            pcshow(obj.xyz)
        end

        function num = getNumInputsImpl(~)
            num = 1;
        end
        function num = getNumOutputsImpl(~)
            num = 0;
        end
        function resetImpl(obj)

        end
    end
  end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function x = project_consensus(varargin)
% PROJECT_CONSENSUS    Bring a set of points into consensus.
%
%   project_consensus(v1,v2,...,vn) returns the elementwise
%   average of v1, ..., vn, i.e., the projection of the vi
%   onto the consensus set.

    N = length(varargin);
    x = zeros(size(varargin{1}));
    for i = 1:N
        x = x + varargin{i};
    end
    x = x./N;
end

author: cvxgrp

repo: proximal

save_path: github-repos/MATLAB/cvxgrp-proximal


text:
function results = TC_ADAL( data, params )

tic;
Y = data.X;
b = data.b;
Omega = data.Omega;

mu = params.mu0;
N = length( size(Y) );
lambda = params.lambda;

Lamb = cell( 1, N );
Xs = cell(1,N);
for i = 1:N
    Lamb{i} = data.Lamb{i};
    Xs{i} = data.X;
end

for iter = 1:params.max_iter
    % solve for X_{N+1} = Y
    Yprev = Y;
    R = tenzeros( size(Y) );
    R(Omega) = lambda*b;
    R = R + ten_sum_all( Lamb ) + ten_sum_all( Xs )/mu;
    Y = R * mu/N;
    Y(Omega) = R(Omega) / (lambda+N/mu);

    % solve for X_i's
%     Xs = cell(1,N);
    for i = 1:N
        Xs{i} = tensor_shrinkage( Y-mu*Lamb{i}, mu, i );
        Lamb{i} = Lamb{i} - (Y-Xs{i})/mu;
    end

    % compute optimality stats
    pinf = 0;
    for i = 1:N
        pinf = pinf + norm( double(tenmat( Xs{i}-Y, i )), 'fro' )^2;
    end
    ynorm = norm(double(tenmat(Y,1)),'fro');
    pinf = sqrt(pinf / N) / ynorm;
%     pinf = pinf / ynorm;
    dinf = norm( double(tenmat(Y - Yprev, 1 )), 'fro' ) / ynorm;
%     dinf = dinf / ynorm;

    % print
    fprintf('Iter: %d,   pinf: %3.2e,   dinf: %3.2e\n', iter, pinf, dinf);

    if pinf < params.opt_tol && dinf < params.opt_tol
        break;
    end

    % update mu
%     if rem( iter, 10 ) == 0
%         mu = mu / 5;
%     end
end

diffNorm = norm( double( tenmat( Y - data.X, 1 )), 'fro' );
trueNorm = norm( double( tenmat( data.X, 1 )), 'fro' );
rel_err = diffNorm / trueNorm;

% save results
results.X = Y;
results.iter = iter;
results.cpu = toc;
results.rel_err = rel_err;

end

author: andrewssobral

repo: lrslibrary


text:
function cost = terminalcosts(k, x, mpcModel, varargin)
    cost = 0.0;
end

author: juchengquan


text:
function cinds = APPgetSpIndsOld(map, imsegs);
% Get pixel indices for each superpixel

nc = max(map(:));
cinds(nc) = struct('inds', [], 'count', 0);

t1 = cputime;
for c = 1:nc
    spind = find(map==c);
    nind = sum(imsegs.npixels(spind));
    cinds(c).inds = zeros(nind, 1);
    cinds(c).count = 0;
end

t2 = cputime;
nseg = imsegs.nseg;
for k = 1:nseg
    ind = find(imsegs.segimage==k);
    c = map(k);
    cinds(c).inds(cinds(c).count+1:cinds(c).count+length(ind)) = ind;
    cinds(c).count = cinds(c).count + length(ind);
end

for c = 1:length(cinds)
    cinds(c).inds = cinds(c).inds(1:cinds(c).count);
end
%segimage = imsegs.segimage;
%for k = 1:numel(segimage)
%    c = map(segimage(k));
%    cinds(c).count = cinds(c).count + 1;
%    cinds(c).inds(cinds(c).count) = k;
%end

author: zouchuhang

repo: LayoutNet


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured_SYM';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';
showBC = true;

HJiter0 = 1;
e2 = 100;
N_holes = [12 5 5];
R_holes = 0.2;
phase_holes = [0 0 0];
nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 0.7;
optimality_initial = 5e-2;
constr_initial = 5e-2;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function trigamma_values_test ( )

%*****************************************************************************80
%
%% TRIGAMMA_VALUES_TEST demonstrates the use of TRIGAMMA_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIGAMMA_VALUES_TEST:\n' );
  fprintf ( 1, '  TRIGAMMA_VALUES stores values of\n' );
  fprintf ( 1, '  the TriGamma function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X               FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = trigamma_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function [fDeltaBValue, fBValueLearning, fBValueObserved] = kj_calcbvalues(mCatalog, fSplitTime, bLearningPeriod, fLearningPeriod, bObservedPeriod, fObservedPeriod, nMinimumNumber, nCalculateMC)
% function [fDeltaBValue, fBValueLearning, fBValueObserved]
%   = kj_calcbvalues(mCatalog, fSplitTime, bLearningPeriod, fLearningPeriod,
%                    bObservedPeriod, fObservedPeriod, nMinimumNumber, nCalculateMC)
% ----------------------------------------------------------------------------------
% Calculation of the b-values in the learning and observation (forecast) period
%
% Input parameters:
%   mCatalog            Earthquake catalog
%   fSplitTime          Time at which the catalog will be split
%   bLearningPeriod     Fix duration of the learning period (1 = fix, 0 = use catalog from start)
%   fLearningPeriod     Duration of the learning period
%   bObservedPeriod     Fix duration of the observation period (1 = fix, 0 = use catalog till end)
%   fObservedPeriod     Duration of the observation period
%   nMinimumNumber      Minimum number of earthquakes in the catalog for calculating the output values
%   nCalculateMC        Method to determine the magnitude of completeness (see also: help calc_Mc)
%
% Output parameters:
%   fDeltaBValue        Difference of both of the b-values (fBValueObserved - fBValueLearning)
%   fBValueLearning     b-value of the learning period
%   fBValueObserved     b-value of the observation (forecast) period
%
% Danijel Schorlemmer
% April 23, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Init output variables
fDeltaBValue = nan;
fBValueLearning = nan;
fBValueObserved = nan;

try
  % Create the catalogs for the learning and observation periods
  [mLearningCatalog, mObservedCatalog] = ex_SplitCatalog(mCatalog, fSplitTime, bLearningPeriod, fLearningPeriod, bObservedPeriod, fObservedPeriod);

  % Determine magnitude of completeness
  fMc = calc_Mc(mLearningCatalog, nCalculateMC);

  % Calculate the b-value of the learning period
  vSelection = mLearningCatalog(:,6) >= fMc;
  mLearningCatalog = mLearningCatalog(vSelection,:);
  if length(mLearningCatalog(:,1)) > nMinimumNumber
    [vDummy, fBValueLearning, vDummy, vDummy] =  bmemag(mLearningCatalog);
  end

  % Calculate the b-value of the observation period
  vSelection = mObservedCatalog(:,6) >= fMc;
  mObservedCatalog = mObservedCatalog(vSelection,:);
  if length(mObservedCatalog(:,1)) > nMinimumNumber
    [vDummy, fBValueObserved, vDummy, vDummy] =  bmemag(mObservedCatalog);
  end

  % Calculate the difference
  if (~isnan(fBValueLearning)) && (~isnan(fBValueObserved))
    fDeltaBValue = fBValueObserved - fBValueLearning;
  end
catch
end

author: CelsoReyes


text:
function gescatter(filename,lon,lat,c,varargin)
% GESCATTER - create a scatter plot in Google Earth
%
%   GESCATTER(FILENAME,LON,LAT,C) - creates a .kml file that
%   displays colored circles at the locations specified by the
%   vectors LON and LAT similar to ML's builtin function, SCATTER.
%   The color of the circles is scaled relative to the
%   values provided in third input, C.
%
%   OPTIONS AND SYNTAX - Optional inputs are entered as
%   property/value pairs.  Valid optional properties are:
%
%   GESCATTER(...,'colormap','hot') - uses Matlabs 'hot'
%   colormap instead of the default (jet). Also accepts function
%   handles (@hot), or custom colormaps (m x 3 matrices).
%
%   GESCATTER(...,'clims',[low high]) - limit the color
%   values to a specified values (similar to CAXIS on a ML
%   figure). Clims should be supplied as a 2-element array.
%
%   GESCATTER(...,'time',timevector) - assigns a time to each
%   point. The length of the timevector array should be the same
%   as LAT, LON, and C.
%
%   GESCATTER(...,'scale',size) - scales the size of the dots in the
%   Google Earth file.  Default value is 0.4.
%
%   EXAMPLE
%
%   %generate some data
%   x=(0:0.05:6*pi);
%   lon = -122.170087 + cos(x)*0.01;
%   lat = 37.455697 + x*0.001;
%
%   %color the points according to their latitude
%   gescatter('foo.kml',lon,lat,lat)
%
% SEE ALSO scatter

% A. Stevens @ USGS 3/04/2009
% astevens@usgs.gov

%default values
clims=[min(c) max(c)];
cmap=fliplr(jet);
t=[];
scale=0.5;

%parse inputs and do some error-checking
if nargin>0
    [m,n]=size(varargin);
    opts={'clims','time','scale','colormap'};
    for i=1:n;
        indi=strcmpi(varargin{i},opts);
        ind=find(indi==1);
        if isempty(ind)~=1
            switch ind
                case 1
                    clims=varargin{i+1};
                    if numel(clims)~=2
                        error('Clims should be a two-element array.')
                    end

                case 2
                    t=varargin{i+1};
                    if any(isnan(t))
                        error('Time vector should not contain NaNs.')
                    end
                    if ~isnumeric(t)
                        error('Time should be entered in ML Datenum format.')
                    end
                case 3
                    scale=varargin{i+1};
                case 4
                    cmap=varargin{i+1};
                    %check size of numeric colormap input
                    if isa(cmap,'numeric')
                        [m,n]=size(cmap);
                        if n~=3
                            error('Custom colormap must have 3 columns.')
                        end
                        cmap=fliplr(cmap);
                    else
                        %if standard colormap is supplied
                        if isa(cmap,'function_handle')
                            cmap= func2str(cmap);
                        end
                        cmap=fliplr(feval(cmap));
                    end

            end
        end
    end
end

[pathstr,namer] = fileparts(filename);

%get rid on nans
gind=(isfinite(lon) & isfinite(lat) & isfinite(c));
lon=lon(gind);
lat=lat(gind);
c=c(gind);

%figure out the rgb colors of each value
cvals=[-inf;linspace(clims(1),clims(2),...
    length(cmap(:,1))-2)';inf];
[n,bin]=histc(c,cvals);
colors=cmap(bin,:);

%convert to GE's hex format
rgb=cellfun(@(x)(dec2hex(floor(x.*255),2)),...
    num2cell(colors),'uni',0);

%write the GE file
header=['<?xml version="1.0" encoding="UTF-8"?>',...
    '<kml xmlns="http://www.opengis.net/kml/2.2">',...
    '<Document><name>',namer,'</name>'];
footer='</Document></kml>';

h = waitbar(0,'Creating file, Please wait...');
set(h,'name','Creating Google Earth file')

fid = fopen(filename, 'wt');
fprintf(fid, '%s \n',header);

for i=1:length(lon)

    %create a style to hide each point in one document
    fprintf(fid,'%s \n','');

    %define the point style
    fprintf(fid,'%s \n','');

    %add the placemark
    fprintf(fid, '%s \n','<Placemark>');
    fprintf(fid,'%s \n','');

    %create a simple description for each point
    fprintf(fid, '%s \n','<description><![CDATA[

');
    fprintf(fid, '%s \n',['
Filename: ',namer,'
']);
    fprintf(fid, '%s \n',['
Value: ',sprintf('%.1f',c(i)),'
']);
    if ~isempty(t)
        fprintf(fid, '%s \n',['Time (GMT): ',datestr(t(i)),'
']);
    end
    fprintf(fid, '%s \n',']]></description>');

    fprintf(fid,'%s \n','<Point>');
    fprintf(fid,'%s','<coordinates>');
    fprintf(fid, ' %.6f, %.6f, %.2f', [lon(i) lat(i) c(i)]);
    fprintf(fid,'%s \n','</coordinates>');
    fprintf(fid,'%s \n','</Point>');

    if ~isempty(t)
        fprintf(fid,'%s \n','<TimeSpan>');
        fprintf(fid,'%s \n',['<begin>',datestr(t(1),29),...
            'T',datestr(t(1),13),'Z</begin>']);
        fprintf(fid,'%s \n',['<end>',datestr(t(end),29),...
            'T',datestr(t(end),13),'Z</end>']);
        fprintf(fid,'%s \n','</TimeSpan>');
    end

    fprintf(fid, '%s \n','</Placemark>');

    waitbar(i/length(lon),h,sprintf('%d%% complete...',...
        round((i/length(lon))*100)));

end

fprintf(fid, '%s \n','');
fprintf(fid, '%s \n',footer);

close(h);
fclose(fid);

repo: mcbench-benchmarks


text:
function out = cfg_run_named_file(job)

% Return selected files as separate output arguments.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_named_file.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

out.files = job.files;
cind = 1;
for k = 1:numel(job.files)
    out.index{k} = cind:(cind+numel(job.files{k})-1);
    cind = cind+numel(job.files{k});
end;

author: spm


text:
function g = restrict(f, s)
%RESTRICT Restrict an UNBNDFUN to a subinterval.
%   RESCTRICT(F, S) returns a FUN object that is restricted to the subinterval
%   [S(1), S(2)] of the domain of F.
%
%   If length(S) > 2, i.e., S = [S1, S2, ..., Sn], then RESCTRICT(F, S) returns
%   an array of FUN objects, where the cells contain F restricted to each of
%   the subintervals defined by S. If there is only one FUN to be returned,
%   that is, if length(S) == 2, then the FUN object g is returned. This
%   facilitates the use of the result by other functions, e.g. plot etc.

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

% Deal with empty case:
if ( isempty(f) )
    g = f;
    return
end

% Check if the argument s is actually a subinterval:
if ( s(1) < f.domain(1) || s(end) > f.domain(2) || any(diff(s) <= 0) )
    error('CHEBFUN:UNBNDFUN:restrict:badInterval', 'Not a valid interval.')
elseif ( numel(s) == 2 && all(s == f.domain) )
    % Nothing to do here!
    g = f;
    return
end

% Number of subdomains:
numSubDom = numel(s) - 1;

% Preallocate the output cell:
g = cell(1, numSubDom);

% Any exponents?
exps = [];

if ( issing(f) )

    % Grab the exponents:
    exps = get(f, 'exponents');

    % If there is a non-trivial exponent:
    if ( any(exps) )

        % Interior exponents:
        interiorExps = zeros(1, numSubDom-1);

        % Insert the interior exponents:
        exps = [exps(1) interiorExps exps(2)];

        % Negate the exponents for infinite endpoint, since the exponents stored
        % in SINGFUN are the negated values of those supplied to the UNBNDFUN
        % constructor:

        if ( s(1) == -Inf )
            exps(1) = -exps(1);
        elseif ( f.domain(1) == -Inf )
            exps(1) = 0;
        end

        if ( s(end) == Inf )
            exps(end) = -exps(end);
        elseif ( f.domain(end) == Inf )
            exps(end) = 0;
        end

    end

end

% Grab the vscale:
vscale = get(f, 'vscale');

% Loop over each subdomain:
for k = 1:numSubDom

    % Initialize preferences:
    pref = chebfunpref();

    % Pass the information about exponents through preference:
    if ( ~isempty(exps) && any(exps(k:k+1)) )
        data.exponents = exps(k:k+1);

        if ( (k == 1) || (k == numSubDom) )
            pref.techPrefs.extrapolate = 1;
        end
    else
        data.exponents = [];
    end

    data.domain = s(k:k+1);
    data.vscale = vscale;
    g{k} = classicfun.constructor(@(x) feval(f, x), data, pref);
end

% When there is only one cell, return the UNBNDFUN instead:
if ( numSubDom == 1 )
    g = g{1};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun