text:
classdef MimLinkedSeriesSidePanel < GemListBoxWithTitle
    % MimLinkedSeriesSidePanel. Part of the gui for the TD MIM Toolkit.
    %
    %     This class is used internally within the TD MIM Toolkit to help
    %     build the user interface.
    %
    %     MimLinkedSeriesSidePanel is part of the side panel and contains the sliding list
    %     box showing linked series
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        PatientDatabase
        GuiCallback
        LinkedRecorder
        GroupPatientsWithSameName
    end

    methods
        function obj = MimLinkedSeriesSidePanel(parent, patient_database, linked_recorder, group_patients_with_same_name, gui_callback)
            obj = obj@GemListBoxWithTitle(parent, 'LINKED SERIES', 'Import images', 'Delete images');

            obj.PatientDatabase = patient_database;
            obj.GuiCallback = gui_callback;
            obj.LinkedRecorder = linked_recorder;
            obj.GroupPatientsWithSameName = group_patients_with_same_name;
        end

        function RepopulateSidePanel(obj, patient_id, series_uid)
            obj.AddSeriesToListBox(patient_id, series_uid);
            if obj.ComponentHasBeenCreated
                obj.Resize(obj.Position);
            end
        end

        function UpdateSidePanel(obj, series_uid)
            obj.ListBox.SelectItem(series_uid, true);
        end

        function SelectSeries(obj, series_uid, selected)
            obj.ListBox.SelectItem(series_uid, selected);
        end
    end

    methods (Access = protected)
        function AddButtonClicked(obj, ~, event_data)
            obj.GuiCallback.ImportMultipleFiles;
        end

        function DeleteButtonClicked(obj, ~, event_data)
            series_uid = obj.ListBox.GetListBox.SelectedTag;
            if ~isempty(series_uid)
                parent_figure = obj.GetParentFigure;
                parent_figure.ShowWaitCursor;
                obj.GuiCallback.UnlinkDataset(series_uid);
                parent_figure.HideWaitCursor;
            end
        end
    end

    methods (Access = private)

        function AddSeriesToListBox(obj, patient_id, series_uid)

            % Get uids for every series associated with this patient
            datasets = obj.PatientDatabase.GetAllSeriesForThisPatient(obj.GuiCallback.GetCurrentProject, patient_id, obj.GroupPatientsWithSameName);
            all_uids = CoreContainerUtilities.GetFieldValuesFromSet(datasets, 'SeriesUid');
            obj.ListBox.ClearItems;

            link_map = obj.LinkedRecorder.LinkMap;
            linked_uids_list = {};
            linked_name_list = {};
            for uid = all_uids
                if link_map.isKey(uid{1})
                    link_record = link_map(uid{1});
                    linked_uids_list{end + 1} = uid{1};
                    linked_name_list{end + 1} = MimLinkedSeriesSidePanel.GuessPrimaryName(link_record);

                    linked_uids_list = [linked_uids_list, link_record.LinkMap.keys];
                    linked_name_list = [linked_name_list, link_record.LinkMap.values];
                end
            end

            if ~isempty(linked_uids_list)
                % Extract only the series which are part of the linked list map
                [uid_found_map, series_index] = ismember(linked_uids_list, all_uids);
                linked_series = datasets(series_index(uid_found_map));

                for series_index = 1 : length(linked_series)
                    series = linked_series{series_index};
                    link_name_text = linked_name_list{series_index};
                    series_item = MimSidePanelLinkedSeriesDescription(obj.ListBox.GetListBox, series.Modality, series.StudyName, series.Name, series.Date, series.Time, series.NumberOfImages, patient_id, series.SeriesUid, link_name_text, obj.GuiCallback);

                    obj.ListBox.AddItem(series_item);
                end

                obj.ListBox.SelectItem(series_uid, true);
            end
        end
    end

    methods (Static, Access = private)
        function name = GuessPrimaryName(link_record)
            % The primary dataset in a set of links does not have a name, but it is useful
            % to assign one as a visual clue in the GUI
            name = 'Primary';
            names = link_record.LinkMap.values;
            if ismember('MR', names) && ismember('CT', names)
                name = 'XE';
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function dq = av2dq(AV)

% AV2DQ   transforms an angular vector into a rotation dual quaternion
%
%       DQ = AV2DQ(AV) returns the rotation dual quaternion DQ from its
%       angular vector representation AV [deg].
%       - AV is the angular vector representation of a rotation of angle
%          theta and unitary vecotr n. AV = theta*n = theta*[nx ny nz]. It
%          is expressed in [deg]. AV is a 3-vector or an 3*N array (where
%          column i is the angular vector representation of rotation i).
%       - DQ is a rotation dual quaternion. It is an 8*N array (each column
%         is a rotation dual quaternion).
%
% See also DQ2AV

sAV = size(AV);
if sAV == [1 3]
    AV = AV.';
    sAV = size(AV);
end

% wrong format
if sAV(1) ~= 3
    error('DualQuaternion:AV2dquat:wrongsize',...
        '%d rows in the AV array. It should be 3.',sAV(1));
end

n = sAV(2);
theta = sym(sqrt(sum(AV.^2))); %default is column sum
axis = sym([ones(1,n); zeros(2,n)]);
axis(:,n) = AV(:,n)./repmat(theta(n),3,1);

dq = rot2dq(theta,axis);

author: Sable

repo: mcbench-benchmarks


text:
function [ n_data, x, fx ] = normal_01_cdf_values ( n_data )

%*****************************************************************************80
%
%% NORMAL_01_CDF_VALUES returns some values of the Normal 01 CDF.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`ContinuousDistributions`"]
%      dist = NormalDistribution [ 0, 1 ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 17;

  fx_vec = [ ...
     0.5000000000000000E+00, ...
     0.5398278372770290E+00, ...
     0.5792597094391030E+00, ...
     0.6179114221889526E+00, ...
     0.6554217416103242E+00, ...
     0.6914624612740131E+00, ...
     0.7257468822499270E+00, ...
     0.7580363477769270E+00, ...
     0.7881446014166033E+00, ...
     0.8159398746532405E+00, ...
     0.8413447460685429E+00, ...
     0.9331927987311419E+00, ...
     0.9772498680518208E+00, ...
     0.9937903346742239E+00, ...
     0.9986501019683699E+00, ...
     0.9997673709209645E+00, ...
     0.9999683287581669E+00 ];

  x_vec = [ ...
     0.0000000000000000E+00, ...
     0.1000000000000000E+00, ...
     0.2000000000000000E+00, ...
     0.3000000000000000E+00, ...
     0.4000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.6000000000000000E+00, ...
     0.7000000000000000E+00, ...
     0.8000000000000000E+00, ...
     0.9000000000000000E+00, ...
     0.1000000000000000E+01, ...
     0.1500000000000000E+01, ...
     0.2000000000000000E+01, ...
     0.2500000000000000E+01, ...
     0.3000000000000000E+01, ...
     0.3500000000000000E+01, ...
     0.4000000000000000E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef CEC2013_F2 < PROBLEM
% <single> <real> <large>
% Shifted Rastrigin's function

%------------------------------- Reference --------------------------------
% X. Li, K. Tang, M. N. Omidvar, Z. Yang, and K. Qin, Benchmark functions
% for the CEC'2013 special session and competition on large-scale global
% optimization, RMIT University, Australia, 2013.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        Xopt;	% Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2013.mat'),'Data');
            obj.Xopt = Data{2};
            obj.M    = 1;
            obj.D    = 1000;
            obj.lower    = zeros(1,obj.D) - 5;
            obj.upper    = zeros(1,obj.D) + 5;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = Rastrigin(Tdiag(Tasy(Tosz(PopDec-repmat(obj.Xopt,size(PopDec,1),1)),0.2),10));
        end
    end
end

function F = Rastrigin(X)
    F = sum(X.^2-10*cos(2*pi*X)+10,2);
end

function Z = Tosz(X)
    X1 = zeros(size(X));
    X1(X~=0) = log(abs(X(X~=0)));
    C1 = zeros(size(X)) + 5.5;
    C1(X>0) = 10;
    C2 = zeros(size(X)) + 3.1;
    C2(X>0) = 7.9;
    Z = sign(X).*exp(X1+0.049*(sin(C1.*X1)+sin(C2.*X1)));
end

function Z = Tasy(X,beta)
    Z = X.^(1+repmat(beta*linspace(0,1,size(X,2)),size(X,1),1).*sqrt(X));
    Z(X<=0) = X(X<=0);
end

function Z = Tdiag(X,alpha)
    Z = X.*repmat(sqrt(alpha).^linspace(0,1,size(X,2)),size(X,1),1);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% APM Report Problem Details
%
% y = apm_details(server,app,x,lam)
%
% This function reports the details of a problem
%   located on the APM server with the following arguments:
%
%   server = address of server
%      app = application name
%        x = values of all variables
%      lam = Lagrange multipliers
%        y = structure with problem information
%
% Reports model details such as values, multipliers, residuals,
%   first derivatives, and second derivatives
%
%      y.nvar  = number of variables
%    y.var.lb  = variable lower bounds
%    y.var.val = variable values
%    y.var.ub  = variable upper bounds
%        y.obj = objective value
%   y.obj_grad = objective gradient
%       y.neqn = number of equations
%     y.eqn.lb = equation lower bounds
%    y.eqn.res = equation residuals
%     y.eqn.ub = equation upper bounds
%        y.jac = jacobian (1st derivatives)
%        y.lam = lagrange multipliers
%    y.hes_obj = second derivative of the objective
% y.hes_eqn{i} = second derivatives of equation {i}
%        y.hes = full Hessian
%
function y = apm_details(server,app,x,lam)

persistent sol_prev sol_count

% initialize sol_count
if (isempty(sol_count)),
    sol_count = 0;
end

% tolerance for checking if the values are the same
%   as the previous call
tol = 1e-10;

% optionally load solution vector of variables
if (nargin>=3),
    % create 'warm' as a column vector
    if (size(x,1)==1),
        warm = x';
    else
        warm = x;
    end

    % see if the values are the same as a previous call
    warm_sum = sum(warm);
    for i = 1:sol_count,
        % check summation of warm first
        if (abs(sol_prev(i).y.sumx-warm_sum) <= tol),
            % check individual elements next
            if (abs(sol_prev(i).y.var.val-warm) <= tol),
                % return previous solution if values haven't changed
                y = sol_prev(i).y;
                return
            end
        end
    end

    % load warm.t0 to the server
    fid = fopen('warm.t0','w');
    fprintf(fid,'%20.12e\n',warm);
    fclose(fid);

    t0_load(server,app,'warm.t0');
    delete('warm.t0');
end

% optionally load lagrange multipliers to server
if (nargin>=4),
    % create 'lam' as a column vector
    if (size(lam,1)==1),
        lam = lam';
    end

    fid = fopen('lam.t0','w');
    fprintf(fid,'%20.12e\n',lam);
    fclose(fid);

    t0_load(server,app,'lam.t0');
    delete('lam.t0');
end

% compute details on server
output = apm(server,app,'solve');

% retrieve variables and bounds
apm_get(server,app,'apm_var.txt');
load apm_var.txt
delete('apm_var.txt');
y.nvar = size(apm_var,1); % get number of variables
y.var.lb = apm_var(:,1);
y.var.val = apm_var(:,2);
y.var.ub = apm_var(:,3);

% retrieve objective function value
apm_get(server,app,'apm_obj.txt');
load apm_obj.txt
delete('apm_obj.txt');
y.obj = apm_obj;

% retrieve objective gradient
apm_get(server,app,'apm_obj_grad.txt');
load apm_obj_grad.txt
delete('apm_obj_grad.txt');
y.obj_grad = apm_obj_grad;

% retrieve equation residuals and bounds
apm_get(server,app,'apm_eqn.txt');
load apm_eqn.txt
delete('apm_eqn.txt');
y.neqn = size(apm_eqn,1); % get number of equations
y.eqn.lb = apm_eqn(:,1);
y.eqn.res = apm_eqn(:,2);
y.eqn.ub = apm_eqn(:,3);

% retrieve jacobian
apm_get(server,app,'apm_jac.txt');
load apm_jac.txt
delete('apm_jac.txt');
jac = apm_jac;
y.jac = sparse(jac(:,1),jac(:,2),jac(:,3),y.neqn,y.nvar);

% retrieve lagrange multipliers
apm_get(server,app,'apm_lam.txt');
load apm_lam.txt
delete('apm_lam.txt');
y.lam = apm_lam;

% retrieve hessian of the objective only
apm_get(server,app,'apm_hes_obj.txt');
load apm_hes_obj.txt
delete('apm_hes_obj.txt');
hs = apm_hes_obj;
y.hes_obj = sparse(hs(:,1),hs(:,2),hs(:,3),y.nvar,y.nvar);

% retrieve hessian of the equations only
apm_get(server,app,'apm_hes_eqn.txt');
load apm_hes_eqn.txt
delete('apm_hes_eqn.txt');
hs = apm_hes_eqn;
nhs = size(apm_hes_eqn,1);
i1(1:y.neqn)=0;
i2(1:y.neqn)=0;
% equations listed in sequential order
for i = 1:nhs,
    if(i1(hs(i,1))==0),
        i1(hs(i,1)) = i;
    end
    i2(hs(i,1)) = i;
end
for i = 1:y.neqn,
    if (i1(i)==0),
        y.hes_eqn{i} = sparse(y.nvar,y.nvar);
    else
        y.hes_eqn{i} = sparse(hs(i1(i):i2(i),2),hs(i1(i):i2(i),3),hs(i1(i):i2(i),4),y.nvar,y.nvar);
    end
end

% construct hessian from obj, lam, and eqn portions
y.hes = y.hes_obj;
for i = 1:y.neqn,
    y.hes = y.hes + y.lam(i) * y.hes_eqn{i};
end

% store values
sol_count = sol_count + 1;
sol_prev(sol_count).y = y;
sol_prev(sol_count).y.sumx = sum(y.var.val);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function K = lfmXlfmKernCompute(lfmKern1, lfmKern2, t1, t2)

% LFMXLFMKERNCOMPUTE Compute a cross kernel between two LFM kernels.
% FORMAT
% DESC computes cross kernel terms between two LFM kernels for
% the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel.
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two LFM kernels for
% the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel.
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : lfmKernParamInit, lfmKernCompute, lfmKernParamInit
%
% COPYRIGHT : David Luengo, 2007, 2008, Mauricio Alvarez, 2008
%
% MODIFICATIONS : Neil D. Lawrence, 2007, 2008,
%

% KERN

if nargin < 4
    t2 = t1;
end
if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end
if lfmKern1.inverseWidth ~= lfmKern2.inverseWidth
    error('Kernels cannot be cross combined if they have different inverse widths.')
end

% Get length scale out.
sigma2 = 2/lfmKern1.inverseWidth;
sigma = sqrt(sigma2);

% Parameters of the kernel
alpha(1) = lfmKern1.damper./(2*lfmKern1.mass);
alpha(2) = lfmKern2.damper./(2*lfmKern2.mass);
omega(1) = sqrt(lfmKern1.spring./lfmKern1.mass - alpha(1)*alpha(1));
omega(2) = sqrt(lfmKern2.spring./lfmKern2.mass - alpha(2)*alpha(2));

% Creation of the time matrices

if isreal(omega)
    % Precomputations to increase speed
    gamma1 = alpha(1) + j*omega(1);
    gamma2 = alpha(2) + j*omega(2);
    preGamma(1) = gamma1 + gamma2;
    preGamma(2) = conj(gamma1) + gamma2;
    preConst = 1./preGamma;
    preExp1 = exp(-gamma1*t1);
    preExp2 = exp(-gamma2*t2);
    % Actual computation of the kernel
    sK = real(lfmComputeH3(gamma1, gamma2, sigma2, t1,t2,preConst, 0, 1) + ...
        lfmComputeH3(gamma2, gamma1, sigma2, t2,t1,preConst(2) - preConst(1), 0, 0).' + ...
        lfmComputeH4(gamma1, gamma2, sigma2, t1, preGamma, preExp2, 0, 1  ) + ...
        lfmComputeH4(gamma2, gamma1, sigma2, t2, preGamma, preExp1,0, 0 ).');
    if lfmKern1.isNormalised
        K0 = (lfmKern1.sensitivity*lfmKern2.sensitivity)/(4*sqrt(2)*lfmKern1.mass*lfmKern2.mass*prod(omega));
    else
        K0 = (sigma*sqrt(pi)*lfmKern1.sensitivity*lfmKern2.sensitivity)/(4*lfmKern1.mass*lfmKern2.mass*prod(omega));
    end
    K = K0*sK;
else
    % Precomputations to increase the speed
    preExp1 = zeros(length(t1),2);
    preExp2 = zeros(length(t2),2);
    gamma1_p = alpha(1) + j*omega(1);
    gamma1_m = alpha(1) - j*omega(1);
    gamma2_p = alpha(2) + j*omega(2);
    gamma2_m = alpha(2) - j*omega(2);
    preGamma(1) = gamma1_p + gamma2_p;
    preGamma(2) = gamma1_p + gamma2_m;
    preGamma(3) = gamma1_m + gamma2_p;
    preGamma(4) = gamma1_m + gamma2_m;
    preConst = 1./preGamma;
    preFactors(1) = preConst(2) - preConst(1);
    preFactors(2) = preConst(3) - preConst(4);
    preFactors(3) = preConst(3) - preConst(1);
    preFactors(4) = preConst(2) - preConst(4);
    preExp1(:,1) = exp(-gamma1_p*t1);
    preExp1(:,2) = exp(-gamma1_m*t1);
    preExp2(:,1) = exp(-gamma2_p*t2);
    preExp2(:,2) = exp(-gamma2_m*t2);
    % Actual computation of the kernel
    sK = (  lfmComputeH3(gamma1_p, gamma1_m, sigma2, t1,t2,preFactors([1 2]), 1) + ...
        lfmComputeH3(gamma2_p, gamma2_m, sigma2, t2,t1,preFactors([3 4]), 1).' + ...
        lfmComputeH4(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExp2, 1 ) + ...
        lfmComputeH4(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExp1, 1 ).');
    if lfmKern1.isNormalised
        K0 = (lfmKern1.sensitivity*lfmKern2.sensitivity)/(8*sqrt(2)*lfmKern1.mass*lfmKern2.mass*prod(omega));
    else
        K0 = (sigma*sqrt(pi)*lfmKern1.sensitivity*lfmKern2.sensitivity)/(8*lfmKern1.mass*lfmKern2.mass*prod(omega));
    end
    K = K0*sK;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%    The following is an implementation of the guided filter (GF) based
%    context enhancement (GFCE) through fusion of infrared and visible
%    images.
%    Ref: Zhiqiang Zhou et al. "Fusion of infrared and visible images for night-vision context
%    enhancement", Applied Optics, 55(23), 2016
%
%    Some of the test images were obtained at
%      http://www.imagefusion.org
%      http://www.ece.lehigh.edu/SPCRL/IF/image_fusion.htm
%
%    Zhiqiang Zhou, Beijing Institute of Technology
%    Apr. 2016
for iname =1:50

% close all;
nLevel = 4;
name = num2str(iname);
%  path_Vis = '.\image\b01_1.tif';      path_IR = '.\image\b01_2.tif';
path_Vis = ['..\..\road\ir\',name,'.jpg']; path_IR = ['..\..\road\vi\',name,'.jpg'];
% path_Vis = '.\image\Trees4906_Vis.jpg'; path_IR = '.\image\Trees4906_IR.jpg';
% path_Vis = '.\image\Octec_Vis.jpg';     path_IR = '.\image\Octec_IR.jpg';
% path_Vis = '.\image\Road_Vis.jpg';      path_IR = '.\image\Road_IR.jpg';
% path_Vis = '.\image\Kayak_Vis.jpg';     path_IR = '.\image\Kayak_IR.jpg';
% path_Vis = '.\image\Steamboat_Vis.jpg'; path_IR = '.\image\Steamboat_IR.jpg';
% path_Vis = '.\image\Trees4917_Vis.jpg'; path_IR = '.\image\Trees4917_IR.jpg';
% path_Vis = '.\image\Dune_Vis.jpg';      path_IR = '.\image\Dune_IR.jpg';

[img1, img2, para.name] = PickName(path_Vis, path_IR);
paraShow1.fig = 'Visible image';
paraShow2.fig = 'Infrared image';

%% ---------- Visibility enhancement for visible image--------------
img1E = Ehn_GF(img1);
img1 = img1E;
%% ---------- Infrared image normalization--------------
mi = min(img2(:));
ma = max(img2(:));
img2 = (img2-mi)/(ma-mi)*255;
%% ---------- Automatic parameter selection --------------
Rs = Relative_PS(img2, img1);
if Rs<0.8
    lambda = 100
else
    if Rs>1.6
        lambda = 2000
    else
        lambda = 2500*Rs - 1900
    end
end

%% ---------- Hybrid multiscale decomposition based on guided filter--------------
sigma = 2;  k = 2;
r0 = 2;     eps0 = 0.1;
l = 2;

M1 = cell(1, nLevel+1);
M1L = cell(1, nLevel+1);
M1{1} = img1/255;
M1L{1} = M1{1};
M1D = cell(1, nLevel+1);
M1E = cell(1, nLevel+1);
sigma0 = sigma;
r = r0;
eps = eps0;
for ii = 2:nLevel+1,

%     % ***using fast guided filter, which has the potential to achieve real-time performance when codes are fully optimized
%     % ***NOTE: large subsampling ratio may cause problem for fusion of some source images.
%     s = max(1, r/2); % subsampling ratio
%     M1{ii} = fastguidedfilter_md(M1{ii-1}, M1{ii-1}, r, 100^2, s);
%     M1L{ii} = fastguidedfilter_md(M1L{ii-1}, M1L{ii-1}, r, eps^2, s);

    M1{ii} = guidedfilter(M1{ii-1}, M1{ii-1}, r, 100^2);
    M1L{ii} = guidedfilter(M1L{ii-1}, M1L{ii-1}, r, eps^2);

    M1D{ii} = M1{ii-1} - M1L{ii};
    M1E{ii} = M1L{ii} - M1{ii};

    sigma0 = k*sigma0;
    r = k*r;
    eps = eps/l;
end

M2 = cell(1, nLevel+1);
M2L = cell(1, nLevel+1);
M2{1} = img2/255;
M2L{1} = M2{1};
M2D = cell(1, nLevel+1);
M2E = cell(1, nLevel+1);
sigma0 = sigma;
r = r0;
eps = eps0;
for ii = 2:nLevel+1,
%     s = max(1, r/2);
%     M2{ii} = fastguidedfilter_md(M2{ii-1}, M2{ii-1}, r, 100^2, s);
%     M2L{ii} = fastguidedfilter_md(M2L{ii-1}, M2L{ii-1}, r, eps^2, s);
    M2{ii} = guidedfilter(M2{ii-1}, M2{ii-1}, r, 100^2);
    M2L{ii} = guidedfilter(M2L{ii-1}, M2L{ii-1}, r, eps^2);

    M2D{ii} = M2{ii-1} - M2L{ii};
    M2E{ii} = M2L{ii} - M2{ii};

    sigma0 = k*sigma0;
    r = k*r;
    eps = eps/l;
end

%% ---------- Fusion --------------

for j = nLevel+1:-1:3
D2 = abs(M2E{j});
D1 = abs(M1E{j});
R = max(D2-D1, 0);
Rmax = max(R(:));
P = R/Rmax;

Cj = atan(lambda*P)/atan(lambda);

sigma_b = 2*sigma0;
if j == nLevel+1
    w = floor(3*sigma_b);
    h = fspecial('gaussian', [2*w+1, 2*w+1], sigma_b);
    lambda0 = lambda;
    Cb = atan(lambda0*P)/atan(lambda0);
    Cb = imfilter(Cb, h, 'symmetric');
    MB = Cb.*M2{nLevel+1} + (1-Cb).*M1{nLevel+1};
end

sigma_c = 1;
w = floor(3*sigma_c);
h = fspecial('gaussian', [2*w+1, 2*w+1], sigma_c);
Cj = imfilter(Cj, h, 'symmetric');

md = Cj.*M2E{j}+ (1-Cj).*M1E{j};
MB = MB + md;
md = Cj.*M2D{j}+ (1-Cj).*M1D{j};
MB = MB + md;
end

sigma_t = 1;
w = floor(3*sigma_t);
h = fspecial('gaussian', [2*w+1, 2*w+1], sigma_t);
C11 = double(abs(M1E{2}) < abs(M2E{2}));
C11 = imfilter(C11, h, 'symmetric');
md = C11.*M2E{2}+ (1-C11).*M1E{2};
MB = MB + md;
C10 = double(abs(M1D{2}) < abs(M2D{2}));
md = C10.*M2D{2}+ (1-C10).*M1D{2};
MB = MB + md;
FI = min(round(MB*275), 255);
FI = max(FI, 0);

paraShow.fig = 'Result';
ShowImageGrad(FI, paraShow,iname);

end

author: thfylsty


text:
% VL_DEMO_SIFT_MATCH  Demo: SIFT: basic matching

pfx = fullfile(vl_root,'figures','demo') ;
randn('state',0) ;
rand('state',0) ;
figure(1) ; clf ;

% --------------------------------------------------------------------
%                                                    Create image pair
% --------------------------------------------------------------------

Ia = imread(fullfile(vl_root,'data','a.jpg')) ;
Ib = imread(fullfile(vl_root,'data','b.jpg')) ;

% --------------------------------------------------------------------
%                                           Extract features and match
% --------------------------------------------------------------------

[fa,da] = vl_sift(im2single(rgb2gray(Ia))) ;
[fb,db] = vl_sift(im2single(rgb2gray(Ib))) ;

[matches, scores] = vl_ubcmatch(da,db) ;

[drop, perm] = sort(scores, 'descend') ;
matches = matches(:, perm) ;
scores  = scores(perm) ;

figure(1) ; clf ;
imagesc(cat(2, Ia, Ib)) ;
axis image off ;
vl_demo_print('sift_match_1', 1) ;

figure(2) ; clf ;
imagesc(cat(2, Ia, Ib)) ;

xa = fa(1,matches(1,:)) ;
xb = fb(1,matches(2,:)) + size(Ia,2) ;
ya = fa(2,matches(1,:)) ;
yb = fb(2,matches(2,:)) ;

hold on ;
h = line([xa ; xb], [ya ; yb]) ;
set(h,'linewidth', 1, 'color', 'b') ;

vl_plotframe(fa(:,matches(1,:))) ;
fb(1,:) = fb(1,:) + size(Ia,2) ;
vl_plotframe(fb(:,matches(2,:))) ;
axis image off ;

vl_demo_print('sift_match_2', 1) ;

author: jianxiongxiao

repo: ProfXkit


text:
function kmlStruct = kml2struct(kmlFile)
% kmlStruct = kml2struct(kmlFile)
%
% Import a .kml file as a vector array of shapefile structs, with Geometry, Name,
% Description, Lon, Lat, and BoundaryBox fields.  Structs may contain a mix
% of points, lines, and polygons.
%
% .kml files with folder structure will not be presented as such, but will
% appear as a single vector array of structs.
%
%

[FID msg] = fopen(kmlFile,'rt');

if FID<0
    error(msg)
end

txt = fread(FID,'uint8=>char')';
fclose(FID);

expr = '<Placemark.+?>.+?</Placemark>';

objectStrings = regexp(txt,expr,'match');

Nos = length(objectStrings);

for ii = 1:Nos
    % Find Object Name Field
    bucket = regexp(objectStrings{ii},'<name.*?>.+?</name>','match');
    if isempty(bucket)
        name = 'undefined';
    else
        % Clip off flags
        name = regexprep(bucket{1},'<name.*?>\s*','');
        name = regexprep(name,'\s*</name>','');
        if strcmp(name(1:2),'<!')
            name = regexprep(name,'<![CDATA[\s*','');
            name = regexprep(name,']]>\s*','');
        end
    end

    % Find Object Description Field
    bucket = regexp(objectStrings{ii},'<description.*?>.+?</description>','match');
    if isempty(bucket)
        desc = '';
    else
        % Clip off flags
        desc = regexprep(bucket{1},'<description.*?>\s*','');
        desc = regexprep(desc,'\s*</description>','');

        % Clip off other flags and rearrange into rows
        if strcmp(desc(1:2),'<!')
            desc = regexprep(desc,'<![CDATA[\s*','');
            desc = regexprep(desc,']]>\s*','');
            % Splitting into rows
            I = find(desc == '<');
            if ~isempty(I)
                desc_c = cell(length(I)+1,1);
                desc_c{1} = desc(1:I(1)-1);
                for i = 2:length(I)
                    desc_c{i} = desc(I(i-1)+4:I(i)-1);
                end
                desc_c{end} = desc(I(end)+4:end);
            end
        end
    end

    % Find style
    bucket = regexp(objectStrings{ii},'<styleUrl.*?>.+?</styleUrl>','match');
    if isempty(bucket)
        StyleUrl = 'undefined';
    else
        % Clip off flags
        StyleUrl = regexprep(bucket{1},'<styleUrl.*?>\s*','');
        StyleUrl = regexprep(StyleUrl,'\s*</styleUrl>','');
    end

    geom = 0;
    % Identify Object Type
    if ~isempty(regexp(objectStrings{ii},'<Point', 'once'))
        geom = 1;
    elseif ~isempty(regexp(objectStrings{ii},'<LineString', 'once'))
        geom = 2;
    elseif ~isempty(regexp(objectStrings{ii},'<Polygon', 'once'))
        geom = 3;
    end

    switch geom
        case 1
            geometry = 'Point';
        case 2
            geometry = 'Line';
        case 3
            geometry = 'Polygon';
        otherwise
            geometry = '';
    end

    % Find Coordinate Field
    bucket = regexp(objectStrings{ii},'<coordinates.*?>.+?</coordinates>','match');
    % Clip off flags
    coordStr = regexprep(bucket{1},'<coordinates.*?>(\s+)*','');
    coordStr = regexprep(coordStr,'(\s+)*</coordinates>','');
    % Split coordinate string by commas or white spaces, and convert string
    % to doubles
    coordMat = str2double(regexp(coordStr,'[,\s]+','split'));
    % Rearrange coordinates to form an x-by-3 matrix

    [m,n] = size(coordMat);
    if m*n >= 3
        coordMat = reshape(coordMat,3,m*n/3)';
    end

    % define polygon in clockwise direction, and terminate
    if exist('poly2ccw')
        [Lat, Lon] = poly2ccw(coordMat(:,2),coordMat(:,1));
    else
        Lat = coordMat(:,2); Lon = coordMat(:,1);
    end
    if geom==3
        Lon = [Lon;NaN];
        Lat = [Lat;NaN];
    end

    % Create structure
    kmlStruct(ii).Geometry = geometry;
    kmlStruct(ii).Name = name;
    if exist('desc_c','var')
        kmlStruct(ii).Description = desc_c;
    else
        kmlStruct(ii).Description = desc;
    end
    kmlStruct(ii).Lon = Lon;
    kmlStruct(ii).Lat = Lat;
    kmlStruct(ii).BoundingBox = [[min(Lon) min(Lat);max(Lon) max(Lat)]];
    kmlStruct(ii).StyleUrl = StyleUrl;
end

author: CHLNDDEV


text:
function str = showdetail(item)

% function str = showdetail(item)
% Display details for a cfg_choice and all of its options.
%
% 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: showdetail.m 4867 2012-08-30 13:04:51Z volkmar $

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

str = showdetail(item.cfg_item);
str{end+1} = 'Class  : cfg_choice';
str{end+1} = ['Set one of the options listed as a struct with a ' ...
              'single field:'];
% Display detailed help for each cfg_choice value item
tags = tagnames(item, true);
for k = 1:numel(tags)
    str{end+1} = sprintf('.%s', tags{k});
end;

author: spm


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function varargout = MU_Load_DICOM_File(varargin)
% MU_LOAD_DICOM_FILE M-file for MU_Load_DICOM_File.fig
%      MU_LOAD_DICOM_FILE, by itself, creates a new MU_LOAD_DICOM_FILE or raises the existing
%      singleton*.
%
%      H = MU_LOAD_DICOM_FILE returns the handle to a new MU_LOAD_DICOM_FILE or the handle to
%      the existing singleton*.
%
%      MU_LOAD_DICOM_FILE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MU_LOAD_DICOM_FILE.M with the given input arguments.
%
%      MU_LOAD_DICOM_FILE('Property','Value',...) creates a new MU_LOAD_DICOM_FILE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MU_Load_DICOM_File_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to MU_Load_DICOM_File_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 MU_Load_DICOM_File

% Last Modified by GUIDE v2.5 18-Sep-2013 17:17:55

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @MU_Load_DICOM_File_OpeningFcn, ...
                   'gui_OutputFcn',  @MU_Load_DICOM_File_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

% --- Executes just before MU_Load_DICOM_File is made visible.
function MU_Load_DICOM_File_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% varargin   command line arguments to MU_Load_DICOM_File (see VARARGIN)

%MatrixUser logo
MU_marks(handles.Matrix_preview_axes,'MatrixUser');

set(handles.File_list,'String',varargin{4});
set(handles.File_list,'Max',max(size(varargin{4})));
handles.pathname=varargin{5};

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes MU_Load_DICOM_File wait for user response (see UIRESUME)
% uiwait(handles.MU_load_DICOM_file);

% --- Outputs from this function are returned to the command line.
function varargout = MU_Load_DICOM_File_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);

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

% --- Executes on selection change in File_list.
function File_list_Callback(hObject, eventdata, handles)

% Hints: contents = get(hObject,'String') returns File_list contents as cell array
%        contents{get(hObject,'Value')} returns selected item from File_list

contents=get(handles.File_list,'String');
if ischar(contents)
        filename=[handles.pathname contents];
else
        selected_files=contents(get(handles.File_list,'Value'));
        filename=[handles.pathname selected_files{1}];
end
try
    info=dicominfo(filename);
    dicom_image=dicomread(filename);
catch me
    warndlg([char(39) filename char(39) ' is not in DICOM format.']);
    return;
end

% dicom_image=double(dicom_image);
set(gcf,'CurrentAxes',handles.Matrix_preview_axes);
imagesc(dicom_image);
colormap gray;
axis off;          % Remove axis ticks and numbers
axis image;        % Set aspect ratio to obtain square pixels

cstr=MU_readStruct(info);
set(handles.DICOM_header_list,'String',cstr);

% --- Executes during object creation, after setting all properties.
function File_list_CreateFcn(hObject, eventdata, handles)

% Hint: listbox 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 File_check_button.
function File_check_button_Callback(hObject, eventdata, handles)

contents=get(handles.File_list,'String');
if iscell(contents)
    selected_files=contents(get(handles.File_list,'Value'));
    set(handles.Selected_file_list,'String',selected_files);
    set(handles.Selected_file_list,'Max',max(size(selected_files)));
else
    selected_files=contents;
    set(handles.Selected_file_list,'String',selected_files);
    set(handles.Selected_file_list,'Max',1);
end

% --- Executes on button press in File_uncheck_button.
function File_uncheck_button_Callback(hObject, eventdata, handles)

if ~isempty(get(handles.Selected_file_list,'String'))
    contents=get(handles.Selected_file_list,'String');
    if iscell(contents)
        ind=get(handles.Selected_file_list,'Value');
        contents(ind)=[];
        set(handles.Selected_file_list,'String',contents);
    else
        set(handles.Selected_file_list,'String',[]);
    end
end
set(handles.Selected_file_list,'Value',1);

% --- Executes on selection change in Selected_file_list.
function Selected_file_list_Callback(hObject, eventdata, handles)

% Hints: contents = get(hObject,'String') returns Selected_file_list contents as cell array
%        contents{get(hObject,'Value')} returns selected item from Selected_file_list

% --- Executes during object creation, after setting all properties.
function Selected_file_list_CreateFcn(hObject, eventdata, handles)

% Hint: listbox 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 DICOM_convert_button.
function DICOM_convert_button_Callback(hObject, eventdata, handles)

contents=get(handles.Selected_file_list,'String');

set(handles.DICOM_convert_button,'Enable','off');
set(handles.File_uncheck_button,'Enable','off');
set(handles.File_check_button,'Enable','off');
set(handles.Load_matrix_button,'Enable','off');
set(handles.Edit_matrix_name,'Enable','off');
[created_matrix,created_matrix_name,errflag]=MU_DICOM2Mat(contents,handles.pathname,handles);
set(handles.DICOM_convert_button,'Enable','on');
set(handles.File_uncheck_button,'Enable','on');
set(handles.File_check_button,'Enable','on');
set(handles.Load_matrix_button,'Enable','on');
set(handles.Edit_matrix_name,'Enable','on');

if errflag==1
    return;
elseif errflag==0
    handles.created_matrices.(created_matrix_name)=created_matrix;
    set(handles.Created_matrix_list,'String',fieldnames(handles.created_matrices));

end
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function Created_matrix_list_CreateFcn(hObject, eventdata, handles)

% Hint: listbox 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 during object creation, after setting all properties.
function Edit_matrix_name_CreateFcn(hObject, eventdata, handles)

% 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 Load_matrix_button.
function Load_matrix_button_Callback(hObject, eventdata, handles)

if ~isfield(handles,'created_matrices')
    errordlg('No matrix was created!');
    return;
end
if isempty(handles.created_matrices)
    errordlg('No matrix was created!');
    return;
end

created_matrix_name=fieldnames(handles.created_matrices);
for i=1:max(size(created_matrix_name))
    MU_load_matrix(created_matrix_name{i}, handles.created_matrices.(created_matrix_name{i}),0);
end

MU_load_DICOM_file_CloseRequestFcn(hObject, eventdata, handles);

% --- Executes when user attempts to close MU_load_DICOM_file.
function MU_load_DICOM_file_CloseRequestFcn(hObject, eventdata, handles)

delete(handles.MU_load_DICOM_file);

% Hint: delete(hObject) closes the figure
% delete(hObject);

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

% Hint: listbox 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_matrix_name_Callback(hObject, eventdata, handles)
% hObject    handle to Edit_matrix_name (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_matrix_name as text
%        str2double(get(hObject,'String')) returns contents of Edit_matrix_name as a double

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [y,z]=amen_sum(X, c, tol, varargin)
%Approximate the linear combination of TT tensors via the AMEn iteration
%   [y,z]=amen_sum(X, v, tol, varargin)
%   Attempts to approximate the y(j) = sum_i X{i}*c(i,j)
%   with accuracy TOL using the AMEn+ALS iteration.
%
%   Options are provided in form
%   'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2 and so
%   on. The parameters are set to default (in brackets in the following)
%   The list of option names and default values are:
%       o y0 - initial approximation to y [rand rank-2]
%       o nswp - maximal number of sweeps [20]
%       o verb - verbosity level, 0-silent, 1-sweep info, 2-block info [1]
%       o kickrank - compression rank of the error,
%         i.e. enrichment size [3]
%       o init_qr - perform QR of the input (save some time in ts, etc) [true]
%       o renorm - Orthog. and truncation methods: direct (svd,qr) or gram
%         (apply svd to the gram matrix, faster for m>>n) [direct]
%       o fkick - Perform solution enrichment during forward sweeps [false]
%         (rather questionable yet; false makes error higher, but "better
%         structured": it does not explode in e.g. subsequent matvecs)
%       o z0 - initial approximation to the error Ax-y [rand rank-kickrank]
%       o can - whether the input is in CP format [false]
%       o multrank - shall we try to grow ranks multiplicatively by adding
%       the previous iterand [false]
%
%
%********
%   For description of adaptive ALS please see
%   Sergey V. Dolgov, Dmitry V. Savostyanov,
%   Alternating minimal energy methods for linear systems in higher dimensions.
%   Part I: SPD systems, http://arxiv.org/abs/1301.6068,
%   Part II: Faster algorithm and application to nonsymmetric systems, http://arxiv.org/abs/1304.1222
%
%   Use {sergey.v.dolgov, dmitry.savostyanov}@gmail.com for feedback
%********
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

% Special case: amen_sum called from another amen package
if (~isempty(varargin))
    v1 = varargin{1};
    if (isa(v1, 'cell'))
        varargin=v1;
    end;
end;

% Default valuse
nswp = 20;
kickrank = 4;
kickrank2 = 0;
verb = 1;
y = [];
z = [];
init_qr = true;
renorm = 'direct';
% renorm = 'gram';
fkick = false;
rmax = Inf;
multrank = false;

can = false; % Whether the input is the canonical format

% Read parameters from input sequence
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'y0'
            y=varargin{i+1};
        case 'z0'
            z=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case 'kickrank2'
            kickrank2=varargin{i+1};
        case 'init_qr'
            init_qr = varargin{i+1};
        case 'renorm'
            renorm = varargin{i+1};
        case 'fkick'
            fkick = varargin{i+1};
        case 'rmax'
            rmax = varargin{i+1};
        case 'can'
            can = varargin{i+1};
        case 'multrank'
            multrank = varargin{i+1};

        otherwise
            warning('Unrecognized option: %s\n',varargin{i});
    end
end

% Convert input to canonical or TT cell array
N = size(c,1);
M = size(c,2);
Xin = X;

if (can)
    N = 1;
    % X itself is already given in the required format
    d = numel(X);
    R = size(X{1}, 2);
    nrms = ones(d,1);
    n = zeros(d,1);
    for i=1:d
        n(i) = size(X{i}, 1);
    end;
    vectype = 1;
else
    R = 1;
    for j=1:N
        if (isa(Xin{j}, 'tt_tensor'))
            if (j==1)
                d = Xin{j}.d;
                n = Xin{j}.n;
                X = cell(d,N);

                % A storage for norms.
                % Z may be kept normalized, but for y we will have y_real = y*prod(nrms).
                % The same for all x. But note that we keep only the _scale_ for all x vecs.
                nrms = ones(d,1);
            end;
            X(:,j) = core2cell(Xin{j});
            vectype = 1; % tt_tensor
        else
            if (j==1)
                d = numel(Xin{j});
                n = zeros(d,1);
                for i=1:d
                    n(i) = size(Xin{j}{i}, 2);
                end;
                X = cell(d,N);
            end;
            X(:,j) = Xin{j};
            vectype = 0; % cell
        end;
    end;
end;

% Initial guess
if (isempty(y))
    init_qr = false;
    [y,ry] = gen_rand(n,d,2);
else
    if (isa(y, 'tt_tensor'))
        y = core2cell(y);
    end;
    ry = ones(d+1,1);
    for i=1:d
        ry(i+1) = size(y{i}, 3);
    end;
end;

% Enrichment vector
if (kickrank+kickrank2>0)
    if (isempty(z))
        [z,rz] = gen_rand(n,d,kickrank+kickrank2);
        init_qr_z = false;
    else
        init_qr_z = true;
        if (isa(z, 'tt_tensor'))
            z = core2cell(z);
        end;
        rz = ones(d+1,1);
        for i=1:d
            rz(i+1) = size(z{i}, 3);
        end;
    end;

    phizx = cell(d+1,N);
    for j=1:N
        phizx{1,j}=ones(1,R); phizx{d+1,j}=ones(R,1);
    end;
    phizy = cell(d+1,1); phizy{1}=1; phizy{d+1}=1;
end;

phiyx = cell(d+1,N);
for j=1:N
    phiyx{1,j}=ones(1,R); phiyx{d+1,j}=ones(R,1);
end;

% Initial ort
for i=1:d-1
    if (init_qr)
        cr = reshape(y{i}, ry(i)*n(i), ry(i+1));
        if (strcmp(renorm, 'gram'))&&(ry(i)*n(i)>5*ry(i+1))
            [cr,~,R]=svdgram(cr);
        else
            [cr,R]=qr(cr, 0);
        end;
        nrmr = norm(R, 'fro');
        if (nrmr>0)
            R = R/nrmr;
        end;
        cr2 = reshape(y{i+1}, ry(i+1), n(i+1)*ry(i+2));
        cr2 = R*cr2;
        ry(i+1) = size(cr, 2);
        y{i} = reshape(cr, ry(i), n(i), ry(i+1));
        y{i+1} = reshape(cr2, ry(i+1), n(i+1), ry(i+2));
    end;
    % We will extract the norms of Y'X
    [phiyx(i+1,:),nrms(i)] = compute_next_Phi(phiyx(i,:), y{i}, X(i,:), 'lr', can);

    if (kickrank+kickrank2>0)
        if (init_qr_z)
            cr = reshape(z{i}, rz(i)*n(i), rz(i+1));
            if (strcmp(renorm, 'gram'))&&(rz(i)*n(i)>5*rz(i+1))
                [cr,~,R]=svdgram(cr);
            else
                [cr,R]=qr(cr, 0);
            end;
            nrmr = norm(R, 'fro');
            if (nrmr>0)
                R = R/nrmr;
            end;
            cr2 = reshape(z{i+1}, rz(i+1), n(i+1)*rz(i+2));
            cr2 = R*cr2;
            rz(i+1) = size(cr, 2);
            z{i} = reshape(cr, rz(i), n(i), rz(i+1));
            z{i+1} = reshape(cr2, rz(i+1), n(i+1), rz(i+2));
        end;

        % But we have to renorm Z'X/|Y'X| to keep Z in the same
        % scale. Hope it will not make |Z'X| too small
        phizx(i+1,:) = compute_next_Phi(phizx(i,:), z{i}, X(i,:), 'lr', can, nrms(i));
        % Z'Y is actually an orthogonal matrix and does not require renorm
        phizy(i+1) = compute_next_Phi(phizy(i), z{i}, y(i), 'lr');
    end;
end;

i = d;
dir = -1;
swp = 1;
max_dx = 0;
y{d} = reshape(y{d}, ry(d)*n(d), ry(d+1));
y{d} = y{d}(:,1:min(ry(d+1), M));
y{d} = [y{d}, zeros(ry(d)*n(d), M-ry(d+1))];
y{d} = reshape(y{d}, ry(d), n(d), 1, M);
ry(d+1) = 1;

while (swp<=nswp)
    % Project the sum
    cry = proj_sum(phiyx(i,:), X(i,:), phiyx(i+1,:), c);
    % It will survive in the terminal block. All we need.
    nrms(i) = norm(cry, 'fro');
    % The main goal is to keep y{i} of norm 1
    if (nrms(i)>0)
        cry = cry/nrms(i);
    else
        nrms(i)=1;
    end;

    % Check stopping criteria
    y{i} = reshape(y{i}, ry(i)*n(i)*ry(i+1), M);
    dx = norm(cry-y{i}, 'fro')/norm(cry, 'fro');
    max_dx = max(max_dx, dx);

    % Truncation and enrichment
    if ((dir>0)&&(i<d))
        cry = reshape(cry, ry(i)*n(i), ry(i+1)*M);
        if (M==1)&&(multrank)
            % Try to accelerate the rank growth
            cry = [cry, reshape(y{i}, ry(i)*n(i), ry(i+1))];
        end;
        if (strcmp(renorm, 'gram'))
            [u,~,v]=svdgram(cry, tol/sqrt(d));
            v = v.';
            if (size(u,2)>rmax)
                u = u(:,1:rmax);
                v = v(:,1:rmax);
            end;
            r = size(u,2);
        else
            [u,s,v]=svd(cry, 'econ');
            s = diag(s);
            r = my_chop2(s, tol*norm(s)/sqrt(d));
            r = min(r,rmax);
            u = u(:,1:r);
            v = conj(v(:,1:r))*diag(s(1:r));
        end;
        if (M==1)&&(multrank)
            % Remove auxiliary previous solution
            v = reshape(v, ry(i+1), 2, r);
            v = v(:,1,:);
            v = reshape(v, ry(i+1), r);
        end;

        % Prepare enrichment, if needed
        if (kickrank+kickrank2>0)
            cry = u*v.';
            cry = reshape(cry, ry(i)*n(i)*ry(i+1), M);
            % For updating z
            crys = cry.';
            crys = reshape(crys, M*ry(i)*n(i), ry(i+1));
            crys = crys*phizy{i+1};
            crys = reshape(crys, M, ry(i)*n(i)*rz(i+1));
            crys = crys.';
            cryz = reshape(crys, ry(i), n(i)*rz(i+1)*M);
            cryz = phizy{i}*cryz;
            cryz = reshape(cryz, rz(i)*n(i)*rz(i+1), M);

            crz = proj_sum(phizx(i,:), X(i,:), phizx(i+1,:), c);
            crz = crz/nrms(i) - cryz;
            nrmz = norm(crz, 'fro');
            crz = reshape(crz, rz(i)*n(i), rz(i+1)*M);
            [crz,~,~]=svd(crz, 'econ');
            crz = crz(:, 1:min(size(crz,2), kickrank));
            if (kickrank2>0)
                crz = [crz, randn(rz(i)*n(i), kickrank2)];
                [crz,~]=qr(crz,0);
            end;
            % For adding into solution
            if (fkick)
                crs = proj_sum(phiyx(i,:), X(i,:), phizx(i+1,:), c);
                crs = crs/nrms(i) - crys;
                crs = reshape(crs, ry(i)*n(i), rz(i+1)*M);
                [crs,~,~]=svd(crs, 'econ');
                crs = crs(:, 1:min(size(crs,2), kickrank));
                u = [u,crs];
                if (strcmp(renorm, 'gram'))&&(ry(i)*n(i)>5*(ry(i+1)+rz(i+1)))
                    [u,~,R]=svdgram(u);
                else
                    [u,R]=qr(u, 0);
                end;
                v = [v, zeros(ry(i+1)*M, size(crs,2))];
                v = v*R.';
                r = size(u, 2);
            end;
        end;
        y{i} = reshape(u, ry(i), n(i), r);

        cr2 = reshape(y{i+1}, ry(i+1), n(i+1)*ry(i+2));
        v = reshape(v, ry(i+1), M*r);
        cr2 = v.'*cr2;
        cr2 = reshape(cr2, M, r*n(i+1)*ry(i+2));
        y{i+1} = reshape(cr2.', r, n(i+1), ry(i+2), M);

        ry(i+1) = r;

        [phiyx(i+1,:), nrms(i)] = compute_next_Phi(phiyx(i,:), y{i}, X(i,:), 'lr', can);

        if (kickrank+kickrank2>0)
            rz(i+1) = size(crz, 2);
            z{i} = reshape(crz, rz(i), n(i), rz(i+1));
            % z{i+1} will be recomputed from scratch in the next step
            phizx(i+1,:) = compute_next_Phi(phizx(i,:), z{i}, X(i,:), 'lr', can, nrms(i));
            phizy(i+1) = compute_next_Phi(phizy(i), z{i}, y(i), 'lr');
        end;
    elseif ((dir<0)&&(i>1))
        cry = reshape(cry.', M*ry(i), n(i)*ry(i+1));
        if (M==1)&&(multrank)
            % Try to accelerate the rank growth
            cry = [cry; reshape(y{i}, ry(i), n(i)*ry(i+1))];
        end;
        [u,s,v]=svd(cry, 'econ');
        s = diag(s);
        r = my_chop2(s, tol*norm(s)/sqrt(d));
        r = min(r,rmax);
        u = u(:,1:r)*diag(s(1:r));
        v = conj(v(:,1:r));

        if (M==1)&&(multrank)
            % Remove auxiliary previous solution
            u = reshape(u, ry(i), 2, r);
            u = u(:,1,:);
            u = reshape(u, ry(i), r);
        end;

        % Prepare enrichment, if needed
        if (kickrank+kickrank2>0)
            cry = u*v.';
            cry = reshape(cry, M, ry(i)*n(i)*ry(i+1));
            % For updating z
            crys = cry.';
            crys = reshape(crys, ry(i), n(i)*ry(i+1)*M);
            crys = phizy{i}*crys;
            crys = reshape(crys, rz(i)*n(i)*ry(i+1), M);
            cryz = reshape(crys.', M*rz(i)*n(i), ry(i+1));
            cryz = cryz*phizy{i+1};
            cryz = reshape(cryz, M, rz(i)*n(i)*rz(i+1));
            cryz = cryz.';
            crz = proj_sum(phizx(i,:), X(i,:), phizx(i+1,:), c);
            crz = crz/nrms(i) - cryz;
            nrmz = norm(crz, 'fro');
            crz = reshape(crz.', M*rz(i), n(i)*rz(i+1));
            [~,~,crz]=svd(crz, 'econ');
            crz = crz(:, 1:min(size(crz,2), kickrank));
            if (kickrank2>0)
                crz = [crz, randn(n(i)*rz(i+1), kickrank2)];
                [crz,~]=qr(crz,0);
            end;
            crz = crz';
            % To add into solution
            crs = proj_sum(phizx(i,:), X(i,:), phiyx(i+1,:), c);
            crs = crs/nrms(i) - crys;
            crs = reshape(crs.', M*rz(i), n(i)*ry(i+1));
            [~,~,crs]=svd(crs, 'econ');
            crs = crs(:, 1:min(size(crs,2), kickrank));
            v = [v,conj(crs)];
            if (strcmp(renorm, 'gram'))&&(ry(i+1)*n(i)>5*(ry(i)+rz(i)))
                [v,~,R]=svdgram(v);
            else
                [v,R]=qr(v, 0);
            end;
            u = [u, zeros(M*ry(i), size(crs,2))];
            u = u*R.';
            r = size(v, 2);
        end;
        y{i} = reshape(v.', r, n(i), ry(i+1));

        cr2 = reshape(y{i-1}, ry(i-1)*n(i-1), ry(i));
        u = reshape(u, M, ry(i)*r);
        u = u.';
        u = reshape(u, ry(i), r*M);
        cr2 = cr2*u;
        cr2 = reshape(cr2, ry(i-1), n(i-1), r, M);
        y{i-1} = cr2;

        ry(i) = r;

        [phiyx(i,:), nrms(i)] = compute_next_Phi(phiyx(i+1,:), y{i}, X(i,:), 'rl', can);

        if (kickrank+kickrank2>0)
            rz(i) = size(crz, 1);
            z{i} = reshape(crz, rz(i), n(i), rz(i+1));
            % z{i+1} will be recomputed from scratch in the next step
            phizx(i,:) = compute_next_Phi(phizx(i+1,:), z{i}, X(i,:), 'rl', can, nrms(i));
            phizy(i) = compute_next_Phi(phizy(i+1), z{i}, y(i), 'rl');
        end;
    else
        y{i} = reshape(y{i}, ry(i), n(i), ry(i+1), M);
    end;

    if (verb>1)
        fprintf('amen-sum: swp=[%d,%d], dx=%3.3e, r=%d, |y|=%3.3e, |z|=%3.3e\n', swp, i, dx, r, norm(cry, 'fro'), nrmz);
    end;

    % Stopping or reversing
    if ((dir>0)&&(i==d))||((dir<0)&&(i==1))
        if (verb>0)
            fprintf('amen-sum: swp=%d{%d}, max_dx=%3.3e, max_r=%d\n', swp, (1-dir)/2, max_dx, max(ry));
        end;
        if ((max_dx<tol)||(swp==nswp))&&(dir>0)
            break;
        end;
        max_dx = 0;
        if (dir>0); swp = swp+1; end;
        dir = -dir;
    else
        i = i+dir;
    end;
end;

y{d} = reshape(cry, ry(d), n(d), M);

% Distribute norms equally...
nrms = exp(sum(log(nrms))/d);
% ... and plug them into y
for i=1:d
    y{i} = y{i}*nrms;
end;

if (vectype==1)
    y = cell2core(tt_tensor, y);
    z = cell2core(tt_tensor, z);
end;

end

function [Phi,nrm] = compute_next_Phi(Phi_prev, x, Y, direction, can, extnrm)
% Performs the recurrent Phi (or Psi) matrix computation
% Phi = Phi_prev * (x'y).
% If direction is 'lr', computes Psi
% if direction is 'rl', computes Phi
% A can be empty, then only x'y is computed.

% Phi1:  rx1, ry1
% Phi2:  ry2, rx2

if (nargin<5)||(isempty(can))
    can = false;
end;
if (nargin<6)
    extnrm = [];
end;

rx1 = size(x,1); n = size(x,2); rx2 = size(x,3);
N = numel(Y);
if (can) % We are working with the canonical format
    Phi = cell(1,1);
    R = size(Y{1},2);
    % Phi1: rx1, R
    % Phi2: R, rx2

    if (strcmp(direction, 'lr'))
        %lr: Phi1
        x = reshape(x, rx1, n*rx2);
        Phi{1} = x'*Phi_prev{1};
        Phi{1} = reshape(Phi{1}, n*rx2, R);
        Y = repmat(Y{1}, rx2, 1); % equalize the sizes for Hadamard product
        Phi{1} = Phi{1}.*Y;
        Phi{1} = reshape(Phi{1}, n, rx2*R);
        Phi{1} = sum(Phi{1}, 1);
        Phi{1} = reshape(Phi{1}, rx2, R);
        if (nargout>1)
            % Extract the scale to prevent overload
            nrm = norm(Phi{1}, 'fro');
            if (nrm>0)
                Phi{1} = Phi{1}/nrm;
            else
                nrm=1;
            end;
        elseif (~isempty(extnrm))
            % Override the normalization by the external one
            Phi{1} = Phi{1}/extnrm;
        end;
    else
        %rl: Phi2
        x = reshape(x, rx1*n, rx2);
        Phi{1} = Phi_prev{1}*x';
        Phi{1} = reshape(Phi{1}, R, rx1*n);
        Y = Y{1}.';
        Y = repmat(Y,rx1,1); % equalize the sizes for Hadamard product
        Y = reshape(Y, R, rx1*n);
        Phi{1} = Phi{1}.*Y;
        Phi{1} = reshape(Phi{1}, R*rx1, n);
        Phi{1} = sum(Phi{1}, 2);
        Phi{1} = reshape(Phi{1}, R, rx1);
        if (nargout>1)
            % Extract the scale to prevent overload
            nrm = norm(Phi{1}, 'fro');
            if (nrm>0)
                Phi{1} = Phi{1}/nrm;
            else
                nrm=1;
            end;
        elseif (~isempty(extnrm))
            % Override the normalization by the external one
            Phi{1} = Phi{1}/extnrm;
        end;
    end;
else % a set of TT-tensors
    Phi = cell(1,N);
    if (strcmp(direction, 'lr'))
        %lr: Phi1
        x = reshape(x, rx1, n*rx2);
        if (nargout>1)
            nrm = 0;
        end;
        for i=1:N
            y = Y{i};
            ry1 = size(y,1); ry2 = size(y,3);
            Phi{i} = x'*Phi_prev{i};
            Phi{i} = reshape(Phi{i}, n, rx2*ry1);
            Phi{i} = Phi{i}.';
            Phi{i} = reshape(Phi{i}, rx2, ry1*n);
            y = reshape(y, ry1*n, ry2);
            Phi{i} = Phi{i}*y;
            Phi{i} = reshape(Phi{i}, rx2, ry2);
            if (nargout>1)
                nrm = max(nrm, norm(Phi{i}, 'fro'));
            end;
        end;
        if (nargout>1)
            % Extract the scale to prevent overload
            if (nrm>0)
                for i=1:N
                    Phi{i} = Phi{i}/nrm;
                end;
            else
                nrm=1;
            end;
        elseif (~isempty(extnrm))
            % Override the normalization
            for i=1:N
                Phi{i} = Phi{i}/extnrm;
            end;
        end;
    else
        %rl: Phi2
        x = reshape(x, rx1, n*rx2);
        if (nargout>1)
            nrm = 0;
        end;
        for i=1:N
            y = Y{i};
            ry1 = size(y,1); ry2 = size(y,3);
            y = reshape(y, ry1*n, ry2);
            Phi{i} = y*Phi_prev{i};
            Phi{i} = reshape(Phi{i}, ry1, n*rx2);
            Phi{i} = Phi{i}*x';
            Phi{i} = reshape(Phi{i}, ry1, rx1);
            if (nargout>1)
                nrm = max(nrm, norm(Phi{i}, 'fro'));
            end;
        end;
        if (nargout>1)
            % Extract the scale to prevent overload
            if (nrm>0)
                for i=1:N
                    Phi{i} = Phi{i}/nrm;
                end;
            else
                nrm=1;
            end;
        elseif (~isempty(extnrm))
            % Override the normalization
            for i=1:N
                Phi{i} = Phi{i}/extnrm;
            end;
        end;
    end;
end;

end

function [x,r]=gen_rand(n,d,r)
% Generate an orthogonal random vector
if (numel(r)==1)
    r = [1; r*ones(d-1,1); 1];
end;
x = cell(d,1);
for i=1:d
    cr = randn(r(i)*n(i), r(i+1));
    [cr,~]=qr(cr,0);
    r(i+1) = size(cr,2);
    x{i} = reshape(cr, r(i), n(i), r(i+1));
end;
end

function [y]=proj_sum(Phi1, X, Phi2, c)
N = size(c,1);
M = size(c,2);

ry1 = size(Phi1{1},1);
ry2 = size(Phi2{1},2);

if (numel(X)==1) % Canonical format
    X = X{1};
    n = size(X,1);
    Phi1 = repmat(Phi1{1}, n, 1); % size ry1*n, R
    X = reshape(X,  1, n*N);
    X = repmat(X, ry1, 1);
    X = reshape(X, ry1*n, N);
    y = X.*Phi1;
    c = repmat(c, ry2, 1); % size N*ry2, M
    c = reshape(c, N, ry2*M);
    Phi2 = repmat(Phi2{1}, 1, M); % size N, ry2*M
    Phi2 = Phi2.*c;
    y = y*Phi2;
    y = reshape(y, ry1*n*ry2, M);
    if (issparse(y))
        y = full(y);
    end;
else
    n = size(X{1}, 2);
    y = zeros(ry1*n*ry2, M);
    for j=1:N
        rx1 = size(X{j},1); rx2 = size(X{j},3);
        cry = reshape(X{j}, rx1, n*rx2);
        cry = Phi1{j}*cry;
        cry = reshape(cry, ry1*n, rx2);
        cry = cry*Phi2{j};
        cry = reshape(cry, ry1*n*ry2, 1);
        cry = cry*c(j,:);
        y = y+cry;
    end;
end;
end

author: oseledets

repo: TT-Toolbox


text:
%YUVCLOSE	Close a YUV4MPEG file
%
%	yuvclose(yuv)
%
%	Close an open YUV file.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function yuvclose(yuv)

	fclose(yuv.fp);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function variance = quasigeometric_variance ( a, b )

%*****************************************************************************80
%
%% QUASIGEOMETRIC_VARIANCE returns the variance of the Quasigeometric PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the probability of 0 successes.
%    0.0 <= A <= 1.0.
%
%    Input, real B, the depreciation constant.
%    0.0 <= B < 1.0.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = ( 1.0 - a ) * ( a + b ) / ( 1.0 - b ) / ( 1.0 - b );

  return
end


text:
function fstats = lme_F(stats,C)
% fstats = lme_F(stats,C)
%
% Inference for the fixed-effects in the linear mixed-effects model(Depends
% on the Statistics toolbox).
%
% Input
% stats: Structure obtained from any of the model fitting functions:
% lme_fit_EM, lme_fit_FS and lme_fit_NR.
% C: Contrast matrix.
%
% Output
% fstats.F: F-Statistic.
% fstats.pval: P-value of the F-Statistic.
% fstats.sgn: Sign of the contrast.
% fstats.df: Degrees of freedom of the F-Statistic.
%
% $Revision: 1.4 $  $Date: 2016/04/08 19:36:05 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2016/04/08 19:36:05 $
%    $Revision: 1.4 $
% References: Bernal-Rusiel J.L., Greve D.N., Reuter M., Fischl B., Sabuncu
% M.R., 2012. Statistical Analysis of Longitudinal Neuroimage Data with Linear
% Mixed Effects Models, NeuroImage, doi:10.1016/j.neuroimage.2012.10.065.
%
if nargin < 2
    error('Too few inputs');
end;
X = stats.X;
if size(C,2) ~= size(X,2)
     error(['The number of colums in C must be equal to the number of ' ...
                                         'colums in the design matrix X']);
end;
p = size(X,2);
ni = stats.ni;
Zcols = stats.Zcols;
Z = X(:,Zcols);
q = length(Zcols);
nth = q*(q+1)/2+1;
W = stats.W;
V = stats.SIGMA;
CBhat = stats.CovBhat;
L = chol(stats.Dhat);
phi = sqrt(stats.phisqhat);
% Computation of Rijs
Rthth = zeros(nth,nth,p,p);
jk = 0;
for k = 1:q
    for j = 1:k
        jk = jk+1;

        uv = 0;
        for v = 1:q
            for u = 1:v
                uv = uv+1;

                posi = 1; SumR = 0;
                for i = 1:length(ni)
                    posf = posi+ni(i)-1;

                    Xi = X(posi:posf,:); Zi = Z(posi:posf,:); Wi = W(posi:posf,1:ni(i));
                    Ekj = zeros(q,q); Ekj(k,j) = 1; Euv = zeros(q,q); Euv(u,v) = 1;
                    Ai = Zi*Ekj*Euv*Zi';
                    Ri = Xi'*Wi*(Ai+Ai')*Wi*Xi;
                    SumR = SumR+Ri;

                    posi = posf+1;
                end
                Rthth(jk,uv,:,:) = SumR;
            end
        end
    end
end
%Computation of Pis,Qijs and the expected information matrix EI.
[EI,Pth,Qthth] = lme_EI(X,Zcols,W,CBhat,V,L,phi,ni);
invEI = EI\eye(nth);
%Estimation of the bias in the covariance matrix and computation of the
%F-statistic and the degrees of freedom of the test.
Bias = 0;
OM = C'*(C*CBhat*C')^-1*C;
A1 = 0; A2 = 0;
Term1 = OM*CBhat;
Term2 = CBhat*OM*CBhat;
for k=1:nth
    Pk = squeeze(Pth(k,:,:));
    for j=1:nth
        Qkj = squeeze(Qthth(k,j,:,:)); Rkj = squeeze(Rthth(k,j,:,:));
        Pj = squeeze(Pth(j,:,:));
        Bias = Bias + invEI(k,j)*(Qkj-Pk*CBhat*Pj-0.25*Rkj);
        A1 = A1+invEI(k,j)*trace(Term1*Pk*CBhat)*trace(Term1*Pj*CBhat);
        A2 = A2+invEI(k,j)*trace(Term1*Pk*Term2*Pj*CBhat);
    end;
end;
szC = size(C,1);
Bdf = (A1+6*A2)/(2*szC);
g = ((szC+1)*A1-(szC+4)*A2)/((szC+2)*A2);
d = 3*szC+2*(1-g);
c1 = g/d;
c2 = (szC-g)/d;
c3 = (szC+2-g)/d;
EF = (1-A2/szC)^-1;
VF = (2/szC)*(1+c1*Bdf)/((1-c2*Bdf)^2*(1-c3*Bdf));
ro = VF/(2*EF^2);
m = 4+(szC+2)/(szC*ro-1);
l = m/(EF*(m-2));
Bias = CBhat*Bias*CBhat;
CovBhat = CBhat + 2*Bias;
Bhat = stats.Bhat;
F = l*Bhat'*C'*(C*CovBhat*C')^-1*C*Bhat/szC;
if F<0
    F = 0;
end;
fstats.F = F;
fstats.pval = max([1-fcdf(F,szC,m), 1e-30]);
fstats.sgn = sign(C*Bhat);
fstats.df = [szC m];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function i4vec_histogram_test ( )

%*****************************************************************************80
%
%% I4VEC_HISTOGRAM_TEST tests I4VEC_HISTOGRAM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_HISTOGRAM_TEST\n' );
  fprintf ( 1, '  I4VEC_HISTOGRAM histograms an integer vector.\n' );

  [ a, seed ] = i4vec_uniform_ab ( n, 0, 25, seed );

  histo_num = 20;

  histo_gram = i4vec_histogram ( n, a, histo_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Histogram of data from 0 to %d\n', histo_num );
  fprintf ( 1, '\n' );

  for i = 0 : histo_num
    if ( 0 < histo_gram(i+1) )
      fprintf ( 1, '  %6d  %6d\n', i, histo_gram(i+1) );
    end
  end

  return
end


text:
function Qb2n = ch_m2q(Cb2n)
% 姿态阵转四元数
%
% Input: Cb2n
% Output: Qb2n
%
    C11 = Cb2n(1,1); C12 = Cb2n(1,2); C13 = Cb2n(1,3);
    C21 = Cb2n(2,1); C22 = Cb2n(2,2); C23 = Cb2n(2,3);
    C31 = Cb2n(3,1); C32 = Cb2n(3,2); C33 = Cb2n(3,3);
    if C11>=C22+C33
        q1 = 0.5*sqrt(1+C11-C22-C33);
        q0 = (C32-C23)/(4*q1); q2 = (C12+C21)/(4*q1); q3 = (C13+C31)/(4*q1);
    elseif C22>=C11+C33
        q2 = 0.5*sqrt(1-C11+C22-C33);
        q0 = (C13-C31)/(4*q2); q1 = (C12+C21)/(4*q2); q3 = (C23+C32)/(4*q2);
    elseif C33>=C11+C22
        q3 = 0.5*sqrt(1-C11-C22+C33);
        q0 = (C21-C12)/(4*q3); q1 = (C13+C31)/(4*q3); q2 = (C23+C32)/(4*q3);
    else
        q0 = 0.5*sqrt(1+C11+C22+C33);
        q1 = (C32-C23)/(4*q0); q2 = (C13-C31)/(4*q0); q3 = (C21-C12)/(4*q0);
    end
    Qb2n = [q0; q1; q2; q3];

author: yandld


text:
function [lambda,v] = eig3(M,a12,a13,a22,a23,a33,varargin)
% eigenvalue and vectors of a symmetric 3x3 matrix
%
% Syntax
%
%   lambda = eig3(M)
%
%   lambda = eig3(a11,a12,a13,a22,a23,a33)
%
%
%   [v,lambda] = eig3(a11,a12,a13,a22,a23,a33)
%
% Input
%  M - array of symmetric 3x3 matrix
%  a11, a12,a13,a22,a23,a33 - vector of matrix elements
%
% Output
%  lambda - eigen values
%  v - eigen vectors
%

% get input
if nargin == 1
  a11 = M(1,1,:); a12 = M(1,2,:); a13 = M(1,3,:);
  a22 = M(2,2,:); a23 = M(2,3,:); a33 = M(3,3,:);
else
  a11 = M;
end

s = size(a11);

% input should be column vectors
a11 = a11(:).'; a12 = a12(:).'; a13 = a13(:).';
a22 = a22(:).'; a23 = a23(:).'; a33 = a33(:).';

% Given a real symmetric 3x3 matrix A, compute the eigenvalues
p1 = a12.^2 +a13.^2 + a23.^2;

q = (a11 + a22 + a33)/3;
p2 = (a11 - q).^2 + (a22 - q).^2 + (a33 - q).^2 + 2 * p1;
p = sqrt(p2 / 6);

%r = det(A-q*Id) / 2 / p^3;
r = (a11-q) .* ( (a22-q) .* (a33-q) - a23.^2) + ...
  a12 .* (a13 .* a23 - a12 .* (a33-q)) + ...
  a13 .* (a12 .* a23 - (a22-q) .* a13);
r = r / 2 ./ p.^3;

% In exact arithmetic for a symmetric matrix  -1 <= r <= 1
% but computation error can leave it slightly outside this range.
phi = acos(r) / 3;
phi(r <= -1) = pi / 3;
phi(r >= 1) = 0;

% the eigenvalues satisfy eig3 <= eig2 <= eig1
lambda = zeros(3,numel(a11));
lambda(1,:) = q + 2 * p .* cos(phi + (2*pi/3));
lambda(3,:) = q + 2 * p .* cos(phi);
lambda(2,:) = 3 * q - lambda(1,:) - lambda(3,:);     % since trace(A) = eig1 + eig2 + eig3

if nargout > 1

  % this is only required for matlab versions prior to 2015
  b11 = repmat(a11,3,1); b12 = repmat(a12,3,1); b13 = repmat(a13,3,1);
  b22 = repmat(a22,3,1); b23 = repmat(a23,3,1); b33 = repmat(a33,3,1);

  v = vector3d(b12 .* b23 - b13 .* (b22 - lambda),...
    b13 .* b12 - (b11 - lambda) .* b23,...
    (b11 - lambda) .* (b22 - lambda) - b12 .* b12,'antipodal');

  v = v.normalize;

  % fallback for special cases
  % TODO: this should be done better
  id = ~(abs(det(v(1,:),v(2,:),v(3,:))) > (1 - 1e-5)) | isnan(lambda(1,:));
  x = v.x; y = v.y; z = v.z;
  for k = find(id)

    [V,l] = eig([a11(k) a12(k) a13(k);...
      a12(k) a22(k) a23(k); ...
      a13(k) a23(k) a33(k)]);

    x(:,k) = V(1,:); y(:,k) = V(2,:); z(:,k) = V(3,:);
    lambda(:,k) = diag(l);

  end
  v.x = x; v.y = y; v.z = z;

  %a1 = vector3d(a11-lambda(1),a12,a13);
  %a2 = vector3d(a12,a22-lambda(1),23);
  %a3 = vector3d(a13,a23,33-lambda(1));

  %cross(a2,a3)

  % return only the largest eigen vector
  if check_option(varargin,'largest'), v = reshape(v(3,:),s); end

  % for some reason Matlab eig function changes to order outputs if called
  % with two arguments - so we should do the same
  [lambda,v] = deal(v,lambda);

end

end

function test

% generate random symmetric 3x3 matrixes
N = 10
a = rand(6,N);

tic
[V1,lambda1] = eig3(a(1,:),a(2,:),a(3,:),a(4,:),a(5,:),a(6,:));
toc

tic
lambda2 = zeros(3,N);
V2 = vector3d.zeros(3,N)
for i = 1:N
  A = [a(1,i),a(2,i),a(3,i);a(2,i),a(4,i),a(5,i);a(3,i),a(5,i),a(6,i)];
  [V,lambda2(:,i)] = eig(A,'vector');
  V2(1,i) = vector3d(V(:,1));
  V2(2,i) = vector3d(V(:,2));
  V2(3,i) = vector3d(V(:,3));
end
toc

norm(lambda1 - flipud(lambda2)) ./ N
max(angle(V1(:),V2(:)))

end

author: mtex-toolbox

repo: mtex

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


text:
function a = i4mat_perm ( n, a, p )

%*****************************************************************************80
%
%% I4MAT_PERM permutes the rows and columns of a square I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Nijenhuis, H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer A(N,N), the matrix to be permuted.
%
%    Input, integer P(N), the permutation.  P(I) is the new number of
%    row and column I.
%
%    Output, integer A(N,N), the permuted matrix.
%
  ierror = perm1_check ( n, p );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_PERM - Fatal error!\n' );
    fprintf ( 1, '  PERM1_CHECK says permutation is illegal.\n' );
    error ( 'I4MAT_PERM - Fatal error!' );
  end

  p = perm_cycle ( n, p, is, nc, 1 );

  for i = 1 : n

    i1 = -p(i);

    if ( 0 < i1 )

      lc = 0;

      while ( 1 )

        i1 = p(i1);
        lc = lc + 1;

        if ( i1 <= 0 )
          break
        end

      end

      i1 = i;

      for j = 1 : n

        if ( p(j) <= 0 )

          j2 = j;
          k = lc;

          while ( 1 )

            j1 = j2;
            it = a(i1,j1);

            while ( 1 )

              i1 = abs ( p(i1) );
              j1 = abs ( p(j1) );

              temp = it;
              it = a(i1,j1);
              a(i1,j1) = it;

              if ( j1 ~= j2 )
                continue
              end

              k = k - 1;

              if ( i1 == i )
                break
              end

            end

            j2 = abs ( p(j2) );

            if ( k == 0 )
              break
            end

          end

        end

      end

    end

  end
%
%  Restore the positive signs of the data.
%
  p(1:n) = abs ( p(1:n) );

  return
end


text:
function Z = xor(X,Y)
%XOR Logical EXCLUSIVE OR for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Z = tenfun(@xor,X,Y);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% This is something I always wondered...
% Manuel Diaz, NTU, 2014.01.21

%
A = magic(6);

%
idx = 1:6;
idx_new = repmat(idx,[6,1]);

%
x = linspace(0,1,6);
v = sin(2*pi*x);

% operation
b = A.*v(idx_new);
disp(b);

% conclusion:
% matlab indexes = c++ pointers ;)

repo: Aero-matlab


text:
function out = callerFile(field)
%CALLERFILE infor of current function's caller file
%USAGE
% s = callerFile
% s
% .name
% .file
% .line
% .path
%
% c = callerFile(field)
%
%

% carry with some debug info
% st(1): funcName st(2): currently running function

ST = dbstack('-completenames', 1);
L = length(ST);

if L < 2
    out.file = [];
    out.name = [];
    out.line = [];
    out.path = cd;
else
    out = ST(2);
    out.path = fileparts(out.file);
end

if nargin ~= 0
	out = out.(field);
end

% carry with some debug info
% st(1): funcName st(2): currently running function

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
function s1 = stirling1 ( n, m )

%*****************************************************************************80
%
%% STIRLING1 computes the Stirling numbers of the first kind.
%
%  Discussion:
%
%    The absolute value of the Stirling number S1(N,M) gives the number
%    of permutations on N objects having exactly M cycles, while the
%    sign of the Stirling number records the sign (odd or even) of
%    the permutations.  For example, there are six permutations on 3 objects:
%
%      A B C   3 cycles (A) (B) (C)
%      A C B   2 cycles (A) (BC)
%      B A C   2 cycles (AB) (C)
%      B C A   1 cycle  (ABC)
%      C A B   1 cycle  (ABC)
%      C B A   2 cycles (AC) (B)
%
%    There are
%
%      2 permutations with 1 cycle, and S1(3,1) = 2
%      3 permutations with 2 cycles, and S1(3,2) = -3,
%      1 permutation with 3 cycles, and S1(3,3) = 1.
%
%    Since there are N! permutations of N objects, the sum of the absolute
%    values of the Stirling numbers in a given row,
%
%      sum ( 1 <= I <= N ) abs ( S1(N,I) ) = N!
%
%  First terms:
%
%    N/M:  1     2      3     4     5    6    7    8
%
%    1     1     0      0     0     0    0    0    0
%    2    -1     1      0     0     0    0    0    0
%    3     2    -3      1     0     0    0    0    0
%    4    -6    11     -6     1     0    0    0    0
%    5    24   -50     35   -10     1    0    0    0
%    6  -120   274   -225    85   -15    1    0    0
%    7   720 -1764   1624  -735   175  -21    1    0
%    8 -5040 13068 -13132  6769 -1960  322  -28    1
%
%  Recursion:
%
%    S1(N,1) = (-1)^(N-1) * (N-1)! for all N.
%    S1(I,I) = 1 for all I.
%    S1(I,J) = 0 if I < J.
%
%    S1(N,M) = S1(N-1,M-1) - (N-1) * S1(N-1,M)
%
%  Properties:
%
%    sum ( 1 <= K <= M ) S2(I,K) * S1(K,J) = Delta(I,J)
%
%    X_N = sum ( 0 <= K <= N ) S1(N,K) X^K
%    where X_N is the falling factorial function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows of the table.
%
%    Input, integer M, the number of columns of the table.
%
%    Output, integer S1(N,M), the Stirling numbers of the first kind.
%
  if ( n <= 0 )
    s1 = [];
    return
  end

  if ( m <= 0 )
    s1 = [];
    return
  end

  s1(1,1) = 1;
  s1(1,2:m) = 0;

  for i = 2 : n

    s1(i,1) = - ( i - 1 ) * s1(i-1,1);

    for j = 2 : m
      s1(i,j) = s1(i-1,j-1) - ( i - 1 ) * s1(i-1,j);
    end

  end

  return
end


text:
function [weight] = sp_kernel(data, weight_o, gamma)
% This program output the final portfolio the SP strategy on data
%
% function [weight] = sp_kernel(data, weight_o, gamma)
%
% weight: final portfolio, used for next rebalance
%
% data: market sequence vectors
% weight_o: last portfolio, also can be last price relative adjusted.
% gamma: switching parameter
%
% Example: [weight] = sp_kernel(data, weight_o, 0.25)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

weight = sp_expert(data, weight_o, gamma);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function imP = ImToPolar (imR, rMin, rMax, M, N)
% IMTOPOLAR converts rectangular image to polar form. The output image is
% an MxN image with M points along the r axis and N points along the theta
% axis. The origin of the image is assumed to be at the center of the given
% image. The image is assumed to be grayscale.
% Bilinear interpolation is used to interpolate between points not exactly
% in the image.
%
% rMin and rMax should be between 0 and 1 and rMin < rMax. r = 0 is the
% center of the image and r = 1 is half the width or height of the image.
%
% V0.1 7 Dec 2007 (Created), Prakash Manandhar pmanandhar@umassd.edu

[Mr Nr] = size(imR); % size of rectangular image
Om = (Mr+1)/2; % co-ordinates of the center of the image
On = (Nr+1)/2;
sx = (Mr-1)/2; % scale factors
sy = (Nr-1)/2;

imP  = zeros(M,  N);

delR = (rMax - rMin)/(M-1);
delT = 2*pi/N;

% loop in radius and
for ri = 1:M
for ti = 1:N
    r = rMin + (ri - 1)*delR;
    t = (ti - 1)*delT;
    x = r*cos(t);
    y = r*sin(t);
    xR = x*sx + Om;
    yR = y*sy + On;
    imP (ri, ti) = interpolate (imR, xR, yR);
end
end

function v = interpolate (imR, xR, yR)
    xf = floor(xR);
    xc = ceil(xR);
    yf = floor(yR);
    yc = ceil(yR);
    if xf == xc & yc == yf
        v = imR (xc, yc);
    elseif xf == xc
        v = imR (xf, yf) + (yR - yf)*(imR (xf, yc) - imR (xf, yf));
    elseif yf == yc
        v = imR (xf, yf) + (xR - xf)*(imR (xc, yf) - imR (xf, yf));
    else
       A = [ xf yf xf*yf 1
             xf yc xf*yc 1
             xc yf xc*yf 1
             xc yc xc*yc 1 ];
       r = [ imR(xf, yf)
             imR(xf, yc)
             imR(xc, yf)
             imR(xc, yc) ];
       a = A\double(r);
       w = [xR yR xR*yR 1];
       v = w*a;
    end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function VOV = generateVOV(obj,theta)
% Generate the virtual objective vectors

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Tomoaki Takagi

    obj    = normalize(obj,'range');
    [N,M]  = size(obj);
    [W,N2] = UniformPoint(2e4,M,'ILD');
    VOV    = zeros(N2,M);
    flag   = false(N2,1);

    normW   = sqrt(sum(W.^2,2));
    normObj = sqrt(sum(obj.^2,2));
    for i = 1 : N2
        CosineVOV = sum(obj.*repmat(W(i,:),N,1),2)./normW(i,:)./normObj;
        d2 = normObj.*sqrt(1-CosineVOV.^2);
        [mind2,I] = min(d2);
        d1 = normObj(I)*CosineVOV(I);
        r = d1/norm(W(i,:));

        VOV(i,:) = W(i,:).*r;
        flag(i)  = mind2 < theta;
    end
    VOV = VOV(flag,:);
end

author: BIMK


text:
% KDTREE_BUILD_DEMO: illustrate the functionalities of kdtree_delete
clc, clear, close all;

mex kdtree_build.cpp
mex kdtree_delete.cpp
disp('compiled.');

p = rand( 10, 3 );
tree = kdtree_build( p );
disp(sprintf('tree structure pointer: %d',tree ));
disp('deleting tree');

kdtree_delete(tree);

disp('tree deleted correctly');

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function vl_demo_kdtree_ann
% VL_DEMO_KDTREE
%   Demonstrates the use of a kd-tree for approximate nearest neighbor
%   (ANN) queries.

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

randn('state',0) ;
rand('state',0) ;

% Generate some 2D data and a query point
X = rand(2, 100) ;
Q = rand(2,1) ;

% Buld a kd-tree
kdtree = vl_kdtreebuild(X) ;

% Query with increasing accuracy
maxNumComparisonRange = [1 10 20 30] ;
for t = [1 2 3 4]
  figure(t) ; clf ;
  vl_plotframe(X, 'ro') ;
  hold on ;
  xl = [.2, .8] ;
  yl = [.1, .7] ;
  xlim(xl) ;
  ylim(yl) ;

  %  vl_demo_kdtree_plot(kdtree, 1, xl, yl) ;

  [i, d] = vl_kdtreequery (kdtree, X, Q, ...
                           'NumNeighbors', 10, ...
                           'MaxComparisons', maxNumComparisonRange(t), ...
                           'Verbose') ;

  vl_plotframe(Q,'b*','markersize',10) ;
  for k=1:length(i)
    if i(k) == 0, continue ; end
    vl_plotframe([Q ; sqrt(d(k))],'b-','linewidth',1) ;
    vl_plotframe(X(:, i(k)), 'bx','markersize',15) ;
  end
  title(sprintf('10 ANNs with at most %d comparisions', maxNumComparisonRange(t))) ;

  axis square ;
  set(gca,'xtick',[],'ytick',[]) ;
  vl_demo_print(t, sprintf('kdtree_ann_%d', t)) ;
end

author: Cloud-CV

repo: object-proposals


text:
function output = relaxMtFitFuncLs(x, m, W_B, W_F, T2_B, R1, S0, t_m, t_s, t_r);
% function output = relaxMtFitFuncLs(x, m, W_B, W_F, T2_B, R1, S0, t_m, t_s, t_r);
% inputs:
%       x = [k, f] where k is the cross-relaxation rate constant defined
%           for the transition from free pool (F) to bound pool (B), f is the
%           fraction of bound spins expressed in terms of concentrations as f =
%           [B]/([B]+[F])
%       m: MT(:,ii) ??
%
%       W_B: Effective saturation rate of the bound pool
%            W_B(ii) = pi*(w1rms^2)*lorentzian (delta(ii), T2_B);% eqtn (2)
%       W_F: 1/R1_F .* Effective saturation rate of the free pool
%            W_F = (w1rms./(2*pi*delta)).^2/.055; --> eqtn(3)./R1_F
%            w1rms = 2400; % omega-1 RMS --> ???? where does this number
%            come from???
%       T2_B: T2 relaxation time of Bound pool, units= seconds
%           T2_B = 11e-6; "average-brain" [YY (2004), pg 411, column2, paragraph 1]
%       R1: observed relaxation rate (measured in the independent
%           experiment) [YY(2004), pg 411, column 1, parag 2]
%           nz = T1>0;
%           R1 = zeros(size(T1));
%           R1(nz) = 1./T1(nz);
%       S0: Synthetic reference image computed by equation (6) using PD(protein Density) and
%           R1 maps.
%       t_m: duration of an off-resonance RF pulse
%           t_m = 8e-3; %bese 8e-3
%       t_s: delay time BEFORE an exitation RF pulse
%           t_s = 5e-3; %bese 5e-3
%       t_r: delay time AFTER an exitation RF pulse
%           t_r = 19e-3; %bese 19e-3
%   Output:
%       output: ??? m/S0 - M_z(:,1)./m_norm(1);
%
%    Example:
%       output=relaxMtFitFunc(x, MT(:,ii), W_B, W_F, T2_B, R1(ii),S0(ii), t_m, t_s, t_r)
%

% C : diagonal matrix = diag(cos(alpha),1) corresponding to instant
% rotation of the magnetization Mz_F by an excitation pulse with a flip
% angle alpha
% cos(10*pi/180) = 0.984807753012
C = [ 0.984807753012 0;
            0        1 ];

k = x(1); % the k-parameter is passed through x
f = x(2); % the f-parameter is passed through x

if(f<=0.01 || f>=0.5 || k<=0.1 || k>=5)
   output = Inf;
   return;
end

%T2_B = 11e-6; %bese 11e-6

R1_B = 1; %bese 1
% Compute this term just one to save some cycles
kf = k*(1-f)/f;
R1_F = R1 - k*(R1_B - R1)/(R1_B - R1 + kf); % eqtn (4)
%W_F = R1_F*(w1rms./(2*pi*delta)).^2/.055; %bese .055
W_F = R1_F*W_F;

R = [ (-R1_F - k)     (kf)    ;
           (k)    (-R1_B - kf) ];

A = R1_F*R1_B + R1_F*kf + R1_B*k;
E_s = expm(R*t_s); % relaxation during delays before (t_s) an exitation RF pulse
E_r = expm(R*t_r); % relaxation during delays after (t_r) an exitation RF pulse

I = eye(2); %identity matrix

for ii = 1:length(W_B);
  W = [-W_F(ii), 0;
       0,   -W_B(ii)];
  E_m = expm((R + W)*t_m);% off-resonance saturation by an RF pulse w/ duration t_m
  D = A + (R1_F + k)*W_B(ii) + (R1_B + kf)*W_F(ii) + W_B(ii)*W_F(ii);
  M_eq = [1-f f]';
  M_ss = 1/D*[(1-f)*(A + R1_F*W_B(ii)) f*(A + R1_B*W_F(ii))]';
  term_1 = inv(I - E_s*E_m*E_r*C);
  %if (rcond(term_1)<1e-4), output = +Inf; return; end;
  term_2 = (E_s*E_m*(I-E_r) + I-E_s)*M_eq;
  term_3 = E_s*(I-E_m)*M_ss;
  M_z(ii,:) = term_1*(term_2 + term_3);
end

% Ova mozebi e nepotrebno, zasto m_norm treba da se zameni so S0
E_m = expm((R)*t_m);
D = A;
M_eq = [1-f f]';
M_ss = 1/D*[(1-f)*(A) f*(A)]';
term_1 = inv(I - E_s*E_m*E_r*C);
term_2 = (E_s*E_m*(I - E_r) + I-E_s)*M_eq;
term_3 = E_s*(I-E_m)*M_ss; %na pocetok treba E_s?
m_norm = term_1*(term_2 + term_3);
result = m/S0 - M_z(:,1)./m_norm(1);
output = sum(result.^2);
%output = result;

return;


text:
function [Qout,fcount] = quadgui_adaptive(F,a,b,tol,varargin)
%QUADGUI  Demonstrate numerical evaluation of a definite integral.
%   Q = QUADGUI(F,A,B) shows the steps in approximating the integral
%   of F(x) from A to B by adaptive extrapolated Simpson's quadrature.
%
%   The shaded area shows the integral over the current subinterval.
%   The color switches to green when the desired accuracy is obtained.
%
%   Q = QUADGUI(F,A,B,tol) uses the given tolerance instead of 1.e-4.
%
%   The first argument, F, is a function handle or an anonymous function
%   that defines F(x).
%   Arguments beyond the first four, Q = QUADGUI(F,a,b,tol,p1,p2,...),
%   are passed on to the integrand, F(x,p1,p2,..).
%
%   [Q,fcount] = QUADGUI(F,...) also counts the number of evaluations
%   of F(x).
%
%   Examples:
%        F                      a     b       tol(optional)
%      humps(x)                 0     1
%      humps(x)                 0     1       1.e-6
%      humps(x)                -1     2
%      sin(x)                   0     pi      1.e-8
%      cos(x)                   0     9*pi/2  1.e-6
%      sqrt(x)                  0     1       1.e-8
%      -sqrt(x)*log(x)          eps   1       1.e-8
%      1/(3*x-1)                0     1
%      t^(8/3)*(1-t)^(10/3)     0     1       1.e-8
%      tan(sin(x))-sin(tan(x))  0     pi
%
%      quadgui(@(x)F,a,b)
%
%   See also QUADTX, QUAD, QUADL, DBLQUAD.

%   Copyright 2014 Cleve Moler
%   Copyright 2014 The MathWorks, Inc.

shg
clf reset
set(gcf,'menubar','none','numbertitle','off','name','Quad gui')

% Default tolerance
if nargin < 4 | isempty(tol)
   tol = 1.e-4;
end

% Default function and interval.
if nargin < 3
   F = @humps;
   a = 0;
   b = 1;
end

% Initialization
c = (a + b)/2;
fa = F(a,varargin{:});
fc = F(c,varargin{:});
fb = F(b,varargin{:});

% Scale the plot
h = b - a;
x = [a c b];
y = [fa fc fb];
maxy = max(y);
miny = min(y);
for k = 1:63
   v = F(a+k*h/64,varargin{:});
   maxy = real(max(maxy,v));
   miny = real(min(miny,v));
end
set(gcf,'userdata',0)
hold on
p(1) = fill(a,fa,'k');
p(2) = fill(b,fb,'k');
p(3) = plot(x,y,'.','markersize',16);
hold off
s = (maxy - miny)/20;
axis([a b miny-s maxy+s])
q(1) = uicontrol('string','step', ...
   'units','normal','pos',[.65 .02 .08 .06], ...
   'callback','set(gcf,''userdata'',1)');
q(2) = uicontrol('string','auto', ...
   'units','normal','pos',[.75 .02 .08 .06], ...
   'callback','set(gcf,''userdata'',2)');
q(3) = uicontrol('string','quit', ...
   'units','normal','pos',[.85 .02 .08 .06], ...
   'callback','set(gcf,''userdata'',3)');

% Recursive call
[Q,k] = quadguistep(F, a, b, tol, fa, fc, fb, varargin{:});
fcount = k + 3;

% Finish
title(sprintf('Q = %8.4f, fcount = %4.0f',Q,fcount))
delete(p(1:2));
delete(q(1:2));
set(q(3),'string','close','callback','close(gcf)')
if nargout > 0, Qout = Q; end

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

function [Q,fcount] = quadguistep(F,a,b,tol,fa,fc,fb,varargin)

% Recursive subfunction used by quadtx.

h = b - a;
c = (a + b)/2;
d = (a + c)/2;
e = (c + b)/2;
fd = F(d,varargin{:});
fe = F(e,varargin{:});
Q1 = h/6 * (fa + 4*fc + fb);
Q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);

u1 = a:h/64:c;
v1 = polyinterp([a d c],[fa fd fc],u1);
u1 = [a u1 c];
v1 = [0 v1 0];
u2 = c:h/64:b;
v2 = polyinterp([c e b],[fc fe fb],u2);
u2 = [c u2 b];
v2 = [0 v2 0];
if (abs(Q2 - Q1) <= tol)
   color = [0 2/3 0];
else
   color = [.6 .6 .6];
end
p = flipud(get(gca,'child'));
x = [get(p(3),'xdata') d e];
y = [get(p(3),'ydata') fd fe];
set(p(1),'xdata',u1,'ydata',v1,'facecolor',color)
set(p(2),'xdata',u2,'ydata',v2,'facecolor',color)
set(p(3),'xdata',x,'ydata',y)
set(gca,'xtick',sort(x),'xticklabel',[]);
title(num2str(length(x)))
pause(.25)
while get(gcf,'userdata') == 0
   pause(.25)
end
if get(gcf,'userdata') == 1
   set(gcf,'userdata',0)
end

if (abs(Q2 - Q1) <= tol) | (get(gcf,'userdata') == 3)
   Q  = Q2 + (Q2 - Q1)/15;
   fcount = 2;
else
   [Qa,ka] = quadguistep(F, a, c, tol, fa, fd, fc, varargin{:});
   [Qb,kb] = quadguistep(F, c, b, tol, fc, fe, fb, varargin{:});
   Q  = Qa + Qb;
   fcount = ka + kb + 2;
end

repo: NumericalAnalysis


text:
function p = pshift ( p, x0, y0 )

%*****************************************************************************80
%
%% PSHIFT shifts a set of points by a given increment.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input/output, real P(NP,2), a set of points to be shifted.
%
%    Input, real X0, Y0, the X and Y increments to be added to each point.
%
  p(:,1) = p(:,1) + x0;
  p(:,2) = p(:,2) + y0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = isequal( f, g )
%ISEQUAL Equality test for SEPARABLEAPPROX.
%
% BOL = ISEQUAL(F, G) returns 0 or 1. If returns 1 then F and G are the same
% SEPARABLEAPPROX, up to relative machine precision. If returns 0 then F and G are
% not the same up to relative machine precision.

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

if ( isempty( f ) )
    if ( isempty( g ) )
        out = true;
    else
        out = false;
    end
    return
end

% Get the low rank representation for f.
fcols = f.cols;
frows = f.rows;
fpiv = f.pivotValues;

% Get the low rank representation for g.
gcols = g.cols;
grows = g.rows;
gpiv = g.pivotValues;

% Test every part:
out = ( isequal(fcols, gcols) & isequal(frows, grows) & isequal(fpiv, gpiv) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


function [err,yy,Value]=energyclassify(L,x,y,xTest,yTest,Kg,varargin);
% function [err,yy,Value]=energyclassify(L,xTr,yTr,xTe,yTe,Kg,varargin);
%
% INPUT:
%  L	:   transformation matrix (learned by LMNN)
%  xTr	:   training vectors (each column is an instance)
%  yTr	:   training labels  (row vector!!)
%  xTe  :   test vectors
%  yTe  :   test labels
%  Kg	:   number of nearest neighbors
%
% Good luck!
%
% copyright by Kilian Q. Weinberger, 2006

% checks
D=length(L);
x=x(1:D,:);
xTest=xTest(1:D,:);
if(size(x,1)>length(L)) error('x and L must have matching dimensions!\n');end;

% set parameters
pars.alpha=1e-09;
pars.tempid=0;
pars.save=0;
pars.speed=10;
pars.skip=0;
pars.factor=1;
pars.correction=15;
pars.prod=0;
pars.thresh=1e-16;
pars.ifraction=1;
pars.scale=0;
pars.obj=0;
pars.union=1;
pars.margin=0;
pars.tabularasa=Inf;
pars.blocksize=500;
pars=extractpars(varargin,pars);

pars

tempname=sprintf('temp%i.mat',pars.tempid);
% Initializationip
[D,N]=size(x);
[gen,NN]=getGenLS(x,y,Kg,pars);

if(pars.scale)
 fprintf('Scaling input vectors!\n');
 sc=sqrt(mean(sum( ((x-x(:,NN(end,:)))).^2)));
 x=x./sc;
 xTest=xTest./sc;
end;

Lx=L*x;
Lx2=sum(Lx.^2);
LxT=L*xTest;

for inn=1:Kg
 Ni(inn,:)=sum((Lx-Lx(:,NN(inn,:))).^2)+1;
end;

MM=min(y);
y=y-MM+1;
un=unique(y);
Value=zeros(length(un),length(yTest));

B=pars.blocksize;
if(size(x,2)>50000) B=250;end;
NTe=size(xTest,2);
for n=1:B:NTe
  fprintf('%2.2f%%: ',n/NTe*100);
  nn=n:n+min(B-1,NTe-n);
  DD=distance(Lx,LxT(:,nn));
 for i=1:length(un)
 % Main Loopfor iter=1:maxiter
  testlabel=un(i);
  fprintf('%i.',testlabel+MM-1);

  enemy=find(y~=testlabel);
  friend=find(y==testlabel);

  Df=mink(DD(friend,:),Kg);
  Value(i,nn)=sumiflessv2(DD,Ni(:,enemy),enemy)+sumiflessh2(DD,Df,enemy)+sum(Df);
%  Value(i,nn)=sumiflessh2(DD,Df+pars.margin,enemy)+sum(Df);
 end;
 fprintf('\n');
end;

 fprintf('\n');
 [temp,yy]=min(Value);

 yy=un(yy)+MM-1;
err=sum(yy~=yTest)./length(yTest);
fprintf('Energy error:%2.2f%%\n',err*100);

function [gen,NN]=getGenLS(x,y,Kg,pars);
fprintf('Computing nearest neighbors ...\n');
[D,N]=size(x);
if(pars.skip) load('.LSKGnn.mat');
else
un=unique(y);
Gnn=zeros(Kg,N);
for c=un
fprintf('%i nearest genuine neighbors for class %i:',Kg,c);
i=find(y==c);
nn=LSKnn(x(:,i),x(:,i),2:Kg+1);
Gnn(:,i)=i(nn);
fprintf('\n');
end;

end;
NN=Gnn;
gen1=vec(Gnn(1:Kg,:)')';
gen2=vec(repmat(1:N,Kg,1)')';

gen=[gen1;gen2];

if(pars.save)
save('.LSKGnn.mat','Gnn');
end;

function NN=LSKnn(X1,X2,ks,pars);
B=2000;
[D,N]=size(X2);
NN=zeros(length(ks),N);
DD=zeros(length(ks),N);

for i=1:B:N
  BB=min(B,N-i);
  fprintf('.');
  Dist=distance(X1,X2(:,i:i+BB));
  fprintf('.');
%  [dist,nn]=sort(Dist);
  [dist,nn]=mink(Dist,max(ks));
  clear('Dist');
  fprintf('.');
%  keyboard;
  NN(:,i:i+BB)=nn(ks,:);
  clear('nn','dist');
  fprintf('(%i%%) ',round((i+BB)/N*100));
end;

function v=vec(M);
% vectorizes a matrix

v=M(:);


text:
function out = cfg_run_save_vars(job)

% Save input variables to .mat file - either as a struct array, or as
% individual variables.
%
% 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_save_vars.m 6929 2016-11-14 13:07:31Z guillaume $

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

[p, n, e] = fileparts(job.name);
out.file{1} = fullfile(job.outdir{1}, [n '.mat']);
for k = 1:numel(job.vars)
    svar.(job.vars(k).vname) = job.vars(k).vcont; %#ok<STRNU>
end;
if job.saveasstruct
    save(out.file{1}, 'svar','-v6');
else
    save(out.file{1}, '-struct','svar','-v6');
end;


text:
function [KB] = PB2KB(PB)
% Convert computery things from petabytes to kilobytes.
% Chad A. Greene 2012
KB = PB*1099511627776;

author: Sable

repo: mcbench-benchmarks


text:
function [signal, lineNoiseOut] = cleanLineNoise(signal, lineNoiseIn)
% Remove sharp spectral peaks from signal using Sleppian filters
%
% Usage:
% signal = cleanLineNoise(signal)
% [signal, lineNoiseOut] = hcleanLineNoise(signal, lineNoiseIn)
%
% Parameters:
%    signal          Structure with .data and .srate fields
%    lineNoiseIn     Input structure with fields described below
%
% Structure parameters (lineNoiseIn):
%    fPassBand       Frequency band used (default [0, Fs/2] = entire band)
%    Fs 	            Sampling frequency
%    fScanBandWidth  +/- bandwidth centered on each f0 to scan for significant
%                       lines (TM)
%    lineFrequencies Line frequencies to be removed (default
%                       [60, 120, 180, 240, 300])
%    lineNoiseChannels  Channels to remove line noise from (default
%                       size(data, 1))
%    maximumIterations   Maximum times to iterate removal (default = 10)
%    p               Significance level cutoff (default = 0.01)
%    pad             FFT padding factor ( -1 corresponds to no padding,
%                       0 corresponds to padding to next highest power of 2
%                       etc.) (default is 0)
%    pnts
%    tapers          Precomputed tapers from dpss
%    taperBandWidth  Taper bandwidth (default 2 Hz)
%    taperWindowSize Taper sliding window length (default 4 sec)
%    taperWindowStep Sliding window step size (default 4 sec = no overlap)
%    tau             Window overlap smoothing factor (default 100)
%
% This function is based on code originally written by Tim Mullen in a
% package called tmullen-cleanline which is based on the chronux_2
% libraries.
%

lineNoiseOut = lineNoiseIn;
%% Remove line frequencies that are greater than Nyquist frequencies
tooLarge = lineNoiseOut.lineFrequencies >= lineNoiseOut.Fs/2;
if any(tooLarge)
    warning('cleanLineNoise:LineFrequenciesTooLarge', ...
        'Eliminating frequencies greater than half the sampling rate');
    lineNoiseOut.lineFrequencies(tooLarge) = [];
    lineNoiseOut.lineFrequencies = squeeze(lineNoiseOut.lineFrequencies);
end

%% Set up multi-taper parameters
hbw = lineNoiseOut.taperBandWidth/2;   % half-bandwidth
lineNoiseOut.taperTemplate = [hbw, lineNoiseOut.taperWindowSize, 1];
Nwin = round(lineNoiseOut.Fs*lineNoiseOut.taperWindowSize); % number of samples in window
lineNoiseOut.tapers = checkTapers(lineNoiseOut.taperTemplate, Nwin, lineNoiseOut.Fs);

%% Perform the calculation for each channel separately
signal.data = double(signal.data);
chans = lineNoiseOut.lineNoiseChannels;
data = signal.data(chans, :);
parfor ch = 1:size(data, 1)
    data(ch, :) = removeLinesMovingWindow(squeeze(data(ch, :)), lineNoiseOut);
end
signal.data(chans, :) = data;
clear data;

author: VisLab

repo: EEG-Clean-Tools