text:
classdef ColorGraphFuxin < AbstractGraph
%COLORGRAPHFUXIN is a variant of ColorGraph, which generates graphs with
% seeds and unaries using color information. This class is derived from
% AbstractGraph. To undestand how unaries are generated, its worth seeing
% how AbstractGraph.prepare_graphs() calls create_unaries() in this class
% to compute unary values.
%
% @authors:     Fuxin Li
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    properties
        precompute
    end

    methods
        function obj = ColorGraphFuxin(seg_obj, segm_index)
            obj.init(seg_obj, segm_index);

            obj.precompute.bg_unary_values = [];
            obj.precompute.clr_unary_hash = {};
            obj.precompute.clr_unary_values = {};
        end

        function create_unaries(obj, graph_sub_method, seed_sets)
            switch graph_sub_method
                case 'internal'
                    obj.generate_color_unaries('internal', seed_sets);
                case 'external'
                    obj.generate_color_unaries('external', seed_sets);
                case 'external2'
                    obj.generate_color_unaries('external2',seed_sets);
                case 'subframe'
                    error('ColorGraph:compute_graph', ...
                          '''subframe'' not implemented yet');
                otherwise
                    error('ColorGraph:compute_graph', ...
                          '''%s'' is an invalid graph method', ...
                          graph_sub_method);
            end
        end

        function generate_color_unaries(obj, internal_external, seed_sets)
            inside_frame = ~obj.seg_obj.sp_frame_set;

            funcs = create_unary_aux_funcs(obj, internal_external);

            [bg_unary_values, clr_unary_vals] = ...
                precompute_return_data(obj, funcs, seed_sets, ...
                                       inside_frame, internal_external);

            % find the offset according to the number of unary graphs
            % already added to graph_unaries_all
            ofst = size(obj.graph_unaries_all.nonlambda_s,2);

            % initialize the matrices for holding all the unaries (see the
            % comment for this structure in AbstractGraph.prepare_graphs())
            obj.graph_unaries_all.nonlambda_s = ...
                [obj.graph_unaries_all.nonlambda_s, zeros(size(seed_sets))];
            obj.graph_unaries_all.nonlambda_t = ...
                [obj.graph_unaries_all.nonlambda_t, zeros(size(seed_sets))];
            obj.graph_unaries_all.lambda_s = ...
                [obj.graph_unaries_all.lambda_s, zeros(size(seed_sets))];
            obj.graph_unaries_all.lambda_t = ...
                [obj.graph_unaries_all.lambda_t, zeros(size(seed_sets))];

            % iterate over all seed locations
            for i = 1:size(seed_sets, 2)
                if strcmp(internal_external, 'internal')
                    [nonlambda_s, lambda_t] = internal_color(obj, ...
                                                   seed_sets(:,i), ...
                                                   clr_unary_vals(:,i));

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                elseif strcmp(internal_external,'external')
                    [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                        external_color(obj, seed_sets(:,i), ...
                                       inside_frame, bg_unary_values, ...
                                       clr_unary_vals(:,i), ...
                                       funcs.parametric_func, ...
                                       funcs.frame_cost_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.nonlambda_t(:,ofst+i) = ...
                                                               nonlambda_t;
                    obj.graph_unaries_all.lambda_s(:,ofst+i) = lambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                else
                    [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                        external_color2(obj, seed_sets(:,i), ...
                                       inside_frame, bg_unary_values, ...
                                       clr_unary_vals(:,i), ...
                                       funcs.parametric_func);

                    obj.graph_unaries_all.nonlambda_s(:,ofst+i) = ...
                                                               nonlambda_s;
                    obj.graph_unaries_all.nonlambda_t(:,ofst+i) = ...
                                                               nonlambda_t;
                    obj.graph_unaries_all.lambda_s(:,ofst+i) = lambda_s;
                    obj.graph_unaries_all.lambda_t(:,ofst+i) = lambda_t;
                end
            end
        end

        function [nonlambda_s, lambda_t] = internal_color(obj, ...
                                                fg_seeds, unary_values)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(fg_seeds),1);
            lambda_t = zeros(length(fg_seeds),1);

            bground_pixels = ~fg_seeds;
            nonlambda_s(fg_seeds) = inf;
            nonlambda_s(bground_pixels) = unary_values(bground_pixels)*2;
            lambda_t(bground_pixels) = 1;
        end

        function [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                external_color(obj, fg_seeds, inside_frame, ...
                               bg_unary_values, fg_unary_values, ...
                               parametric_func, frame_cost_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(fg_seeds),1);
            nonlambda_t = zeros(length(fg_seeds),1);
            lambda_s = zeros(length(fg_seeds),1);
            lambda_t = zeros(length(fg_seeds),1);

            non_seed_pixels = inside_frame & ~fg_seeds;
            lambda_weights = parametric_func(non_seed_pixels);
            bground_pixels = ~(inside_frame | obj.seg_obj.sp_frame_set);

            % internal seed - topological constraint
            nonlambda_s(fg_seeds) = inf;
            % external seed - topological constraint
            nonlambda_t(bground_pixels) = inf;

            % not we dont need to add values because non_seed_pixels is a
            % subset of inside_frame and bground_pixels is outside
            % inside_frame
            nonlambda_t(non_seed_pixels) = bg_unary_values(non_seed_pixels) * 2;
            % unary potential reflecting similarity to fg_seeds pixels
            nonlambda_s(non_seed_pixels) = fg_unary_values(non_seed_pixels) * 2;
            % parametric potential
            lambda_s(non_seed_pixels) = lambda_weights;

            frame_cost = frame_cost_func(fg_seeds);

%           frame_cost = obj.frame_bg_bias(fg_seeds, size(obj.I));
            %frame_cost = 1;
            % always have some frame soft bias to belong in the background
            nonlambda_t(obj.seg_obj.sp_frame_set) = ...
                nonlambda_t(obj.seg_obj.sp_frame_set) + ...
                frame_cost*obj.SEED_FRAME_WEIGHT;
            % parametric potential
            lambda_t(non_seed_pixels) = lambda_weights;
        end

        function [nonlambda_s, nonlambda_t, lambda_s, lambda_t] = ...
                external_color2(obj, fg_seeds, inside_frame, ...
                               bg_unary_values, fg_unary_values, ...
                               parametric_func)
            % see the comment for these output matrices in
            % AbstractGraph.prepare_graphs()
            nonlambda_s = zeros(length(fg_seeds),1);
            nonlambda_t = zeros(length(fg_seeds),1);
            lambda_s = zeros(length(fg_seeds),1);
            lambda_t = zeros(length(fg_seeds),1);

            non_seed_pixels = ~fg_seeds;
            lambda_weights = parametric_func(non_seed_pixels);
            bground_pixels = ~(inside_frame | obj.seg_obj.sp_frame_set);

            % internal seed - topological constraint
            nonlambda_s(fg_seeds) = inf;
            % external seed - topological constraint
            nonlambda_t(bground_pixels) = inf;

            % not we dont need to add values because non_seed_pixels is a
            % subset of inside_frame and bground_pixels is outside
            % inside_frame
            nonlambda_t(non_seed_pixels) = bg_unary_values(non_seed_pixels) * 2;
            % unary potential reflecting similarity to fg_seeds pixels
            nonlambda_s(non_seed_pixels) = fg_unary_values(non_seed_pixels) * 2;
            % parametric potential
            lambda_s(non_seed_pixels) = lambda_weights;

%             frame_cost = frame_cost_func(fg_seeds);
%
% %           frame_cost = obj.frame_bg_bias(fg_seeds, size(obj.I));
%             %frame_cost = 1;
%             % always have some frame soft bias to belong in the background
%             nonlambda_t(obj.seg_obj.sp_frame_set) = ...
%                 nonlambda_t(obj.seg_obj.sp_frame_set) + ...
%                 frame_cost*obj.SEED_FRAME_WEIGHT;
            % parametric potential
            lambda_t(non_seed_pixels) = lambda_weights;
        end

        function [bg_unary_values, clr_unary_vals] = ...
                precompute_return_data(obj, funcs, seed_sets, ...
                                       inside_frame, unary_desc)

            % precompute the non-changing bg color unary values
            if isempty(obj.precompute.bg_unary_values)
                % unary potential reflecting similarity to bground pixels
                bndry_for_estimation = funcs.bndry_est_func(inside_frame);
                obj.precompute.bg_unary_values = ...
                    funcs.color_unary_func(bndry_for_estimation);
            end
            bg_unary_values = obj.precompute.bg_unary_values;

            seed_hash = num2str(sum(getByteStreamFromArray(seed_sets)));
%            seed_hash = DataHash(seed_sets, ...
%                                 obj.seg_obj.other_params.hash_params);
            % check if this seed type has not already been computed
            precomputed_unary_idx = find(strcmp(seed_hash, ...
                                         obj.precompute.clr_unary_hash));
            % precompute the fg color unary values
            if isempty(precomputed_unary_idx)
                clr_unary_values = zeros(size(seed_sets));
                % iterate over all seeds - compute unary values for each sp
                for i = 1:size(seed_sets, 2)
                    clr_unary_values(:,i) = ...
                        funcs.color_unary_func(seed_sets(:,i));
                end
                obj.precompute.clr_unary_values{end+1} = clr_unary_values;
                obj.precompute.clr_unary_hash{end+1} = seed_hash;
                precomputed_unary_idx = ...
                    length(obj.precompute.clr_unary_values);
            end
            clr_unary_vals = ...
                obj.precompute.clr_unary_values{precomputed_unary_idx};

            % if using the color segmenter for 'internal' type, normalize
            % by superpixel size
            if strcmpi(unary_desc, 'internal')
                clr_unary_vals = bsxfun(@rdivide, clr_unary_vals, ...
                                        obj.seg_obj.sp_data.sp_seg_szs);
            end
        end
    end
end

author: Cloud-CV

repo: object-proposals


text:
function [L, V] = model_sort(m, i, L, V)
% Topological sort of the nonterminal symbols in m's grammar.
%   [L, V] = model_sort(m, i, L, V)
%
% Return values
%   L   Symbols visited in post order
%   V   (internal use) Symbol visitation status
%
% Arguments
%   m   Object model
%   i   (internal use) Current symbol
%   L   (internal use) Symbols visited thus far in post order
%   V   (internal use) Symbol visitation status thus far

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% initialize depth-first search at start symbol
if nargin < 2
  i = m.start;
  L = [];
  V = zeros(m.numsymbols, 1);
end

% check for cycle containing symbol i
if V(i) == 1
  error('Cycle detected in grammar!');
end

% mark symbol i as pre-visit
V(i) = 1;
for r = 1:length(m.rules{i})
  for s = m.rules{i}(r).rhs
    % recurse if s is a nonterminal and not already visited
    if m.symbols(s).type == 'N' && V(s) < 2
      [L, V] = model_sort(m, s, L, V);
    end
  end
end
% mark symbol i as post-visit
V(i) = 2;
L = [L i];

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function G = gsp_graph(W, coords, limits)
%GSP_GRAPH  Create a graph given weighted adjacency matrix
%   Usage:  G = gsp_graph(W);
%           G = gsp_graph(W, coords);
%           G = gsp_graph(W, coords, limits);
%
%   Input parameters:
%         W     : (n by n) Weighted adjacency matrix
%         coords: (n by 2) or (n by 3) Coordinates of the points (optional)
%         limits: limits for the coordinates (optional)
%   Output parameters:
%         G     : Graph structure.
%
%   'gsp_graph(W, coords, limits)' initializes a graph structure with W as
%   weight matrix.
%
%   Example:::
%
%          W = rand(10);
%          W = W - diag(diag(W));
%          W = (W + W')/2;
%          G = gsp_graph(W);
%

% Author: Nathanael Perraudin
% Date: 16 March 2014

gsp_check_weights(W);

G.W = W;

% Create coordinates
if nargin > 1
    G.coords = coords;
end
if nargin > 2
    G.plotting.limits = limits;
end

G.type = 'from weight';

G = gsp_graph_default_parameters(G);

end

repo: gspbox


text:
function [rotationAngle,X,Y,finalImage,errors,finalOriginalImage] = ...
            alignTwoImages(image1,image2,angleGuess,spacing,...
                    fractionalPixelAccuracy,noRotation,originalImage)
%alignTwoImages rotationally and translationally aligns an image with a
%background image
%
%   Input variables:
%
%       image1 -> background image
%       image2 -> image to be aligned
%       angleGuess -> initial guess to eliminate 180 degree degeneracy
%       spacing -> angular spacing in Radon transform
%       fractionalPixelAccuracy -> accuracy of translational alignment
%       noRotation -> true if only translational alignment used
%       originalImage -> full version of image2 (optional)
%
%
%   Output variables:
%
%       rotationAngle -> rotational alignment angle
%       X, Y -> translational alignment values (in pixels)
%       finalImage -> aligned version of image2
%       errors -> errors in alignment
%       finalOriginalImage -> aligned version of originalImage (optional)
%
%
% (C) Gordon J. Berman, 2014
%     Princeton University

    if nargin < 3 || isempty(angleGuess) == 1
        angleGuess = 0;
    else
        angleGuess = mod(angleGuess,360);
    end

    angleGuess = angleGuess*pi/180;

    if nargin < 4 || isempty(spacing) == 1
        spacing = .5;
    end
    N = 180/spacing;

    if nargin < 5 || isempty(fractionalPixelAccuracy) == 1
        fractionalPixelAccuracy = .25;
    end

    if nargin < 6 || isempty(noRotation)
        noRotation = false;
    end

    if nargin < 7
        originalImage = [];
        finalOriginalImage = [];
    end

    errors = zeros(2,1);

    s = size(image1);

    if ~noRotation

        thetas = linspace(0, 180-spacing, N);

        %Find fft of the Radon transform
        F1 = abs(fft(radon(image1, thetas)));
        F2 = abs(fft(radon(image2, thetas)));

        %Find the index of the correlation peak
        correlation = sum(fft2(F1) .* fft2(F2));
        peaks = real(ifft(correlation));
        peakIndex = find(peaks==max(peaks));

        if length(peakIndex) > 1
            peakIndex = peakIndex(1);
        end

        %Find rotation angle via quadratic interpolation
        if (peakIndex~=1) && (peakIndex ~= N)
            p=polyfit(thetas((peakIndex-1):(peakIndex+1)),peaks((peakIndex-1):(peakIndex+1)),2);
            rotationAngle = -.5*p(2)/p(1);
            errors(1) = polyval(p,rotationAngle);
        else
            if peakIndex == 1
                p = polyfit([thetas(end)-180,thetas(1),thetas(2)],peaks([N,1,2]),2);
                rotationAngle = -.5*p(2)/p(1);
                errors(1) = polyval(p,rotationAngle);
                if rotationAngle < 0
                    rotationAngle = 180 + rotationAngle;
                end
            else
                p = polyfit([thetas(end-1),thetas(end),180+thetas(1)],peaks([N-1,N,1]),2);
                rotationAngle = -.5*p(2)/p(1);
                errors(1) = polyval(p,rotationAngle);
                if rotationAngle >= 180
                    rotationAngle = rotationAngle - 180;
                end
            end
        end

        %Check to see if rotation angle is in the correct direction
        rA = rotationAngle*pi/180;
        test = dot([cos(rA),sin(rA)],[cos(angleGuess),sin(angleGuess)]);
        if test < 0
            rotationAngle = mod(rotationAngle-180,360);
        end
        rotationAngle = mod(rotationAngle,360);
        toRotate = mod(-rotationAngle,360);

        %Rotate Image & Crop to original Size
        rotatedImage = imrotate(image2,toRotate,'crop');

    else

        rotationAngle = mod(angleGuess,360);
        toRotate = mod(-rotationAngle,360);
        rotatedImage = imrotate(image2,toRotate,'crop');

    end

    % Take 2D FFT of each image
    F1 = fft2(image1);
    F2 = fft2(rotatedImage);

    shifts = dftregistration(F1,F2,round(1/fractionalPixelAccuracy));
    X = shifts(4);
    Y = shifts(3);

    errors(2) = shifts(1);

    T = maketform('affine',[1 0 0 ;0 1 0;X Y 1]);
    if nargout > 3
        finalImage = imtransform(rotatedImage,T,'XData',[1 s(2)],'YData',[1 s(1)]);
    end

    if isempty(originalImage) == 0
        rotatedImage2 = imrotate(originalImage,toRotate,'crop');
        finalOriginalImage = imtransform(rotatedImage2,T,'XData',[1 s(2)],'YData',[1 s(1)]);
    end

author: gordonberman

repo: MotionMapper


text:
function msh = mrmInitHostWindow(msh)
% Initialize the host and window for a mesh.
%
%   msh = mrmInitHostWindow(msh)
%
% The msh structure contains slots for the windowID and the local host. If
% the currently identified windowID is open and the server is running on
% the host, nothing happens.
%
% If the current windowID is not open, then a new one is opened and the
% windowID field is updated.
%
% If the server is not running and the window is not opened, both are
% started and assigned.
%
% Examples:
%   msh = viewGet(VOLUME{1},'mesh');
%   msh = mrmInitHostWindow(msh);
%
% BW (c) Stanford VISTASOFT Team, many moons ago

if ieNotDefined('msh'), error('You must define a mesh.'); end

windowID = meshGet(msh,'id');
host = meshGet(msh,'host');

% If the server is not started, start it with windowID 1
if ~mrmCheckServer,  windowID = mrmStart(1); end

% Open the window.
msh = meshSet(msh,'windowID',mrmOpenWindow(host,windowID));

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Demo for 2D neural field model.
% This demo inverts a linear 2D neural field model.
% NB: in this instance, the neural field is deterministic.

clear variables
close all

%---- Choose basic settings for simulations ----%

% dimensions
n_t             = 2e2;                      % number of time points
ns              = 5;                        % size of the square grid
deltat          = 5e-2;                     % time discretization
f_fname         = @f_2DneuralField;         % Evolution function
g_fname         = @g_Id;                    % Observation function

% Sinusoidal input to the approximate center of the neural field
u               = zeros(ns^2,n_t);
ohmega          = 0.91;
centre          = 0.5*ns*(ns-1);
u(centre-1:centre+1,2)     = 10;
u(centre-1:centre+1,:)     = repmat(1*sin(ohmega*deltat*[1:n_t]),3,1);

% parameters of the simulation
alpha           = Inf;%1e6;
sigma           = 1e3;%Inf;
theta           = [0.15;1;1];
phi             = [1;1;0.1];

% Build options structure for temporal integration of SDE
I = speye(ns,ns);
E = sparse(2:ns,1:ns-1,1,ns,ns);
D = E+E'-2*I;
inF.L           = kron(D,I)+kron(I,D);
inF.deltat      = deltat;
inG.ind         = 1:ns^2*2; % this makes states observable (not temporal derivatives)
options.inF     = inF;
options.inG     = inG;
% options.u0      = 0*ones(ns^2,1);

% Build priors for model inversion
priors.muX0 = 0*ones(ns^2*2,1);
priors.SigmaX0 = 0*speye(ns^2*2);%1e-1*speye(ns^2*2);
priors.muTheta = [.5;.5;.5];
priors.SigmaTheta = 1e0*speye(size(theta,1));
priors.muPhi = [0.1;0.9;0];
priors.SigmaPhi = 1e0*speye(size(phi,1));
priors.SigmaPhi(2,2) = 1e-4;
priors.a_alpha = Inf;%1e3;
priors.b_alpha = 0;%1e1;
priors.a_sigma = 1e3;
priors.b_sigma = 1e1;

% Build options and dim structures for model inversion
options.priors      = priors;
dim.n_theta         = size(theta,1);
dim.n_phi           = size(phi,1);
dim.n               = (ns^2)*2;

% Build time series of hidden states and observations
[y,x,x0,eta,e] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options,priors.muX0);

% display time series of hidden states and observations
displaySimulations(y,x,eta,e);
% disp('--paused--')
% pause

% Call inversion routine
% [posterior,out] = VBA_onlineWrapper(y,u,f_fname,g_fname,dim,options);
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% Display results
displayResults(posterior,out,y,x,x0,theta,phi,alpha,sigma);

% Make predictions
try
    options = out.options;
    [xs,ys,xhat,vx,yhat,vy] = VBA_comparePredictions(...
        n_t,theta,phi,zeros(size(u)),alpha,sigma,options,posterior,dim);
catch
    disp('------!!Unable to form predictions!!------')
end

author: MBB-team

repo: VBA-toolbox


text:
classdef LaunchVehicleBasicElectricalBatteryState < AbstractLaunchVehicleEpsStorageState
    %LaunchVehicleBasicElectricalBatteryState Summary of this class goes here
    %   Detailed explanation goes here

    properties
        stageState = LaunchVehicleStageState.empty(1,0)% LaunchVehicleStageState
        battery% LaunchVehicleBasicElectricalBattery

        stateOfCharge = 0%(1,1) double = 0; %mT
        active(1,1) logical = true;
    end

    methods
        function obj = LaunchVehicleBasicElectricalBatteryState(stageState, battery)
            obj.stageState = stageState;
            obj.battery = battery;
        end

        function epsStorageComponent = getEpsStorageComponent(obj)
            epsStorageComponent = [obj.battery];
        end

        function active = getActiveState(obj)
            active = [obj.active];
        end

        function setActiveState(obj,active)
            obj.active = active;
        end

        function stateOfCharge = getStateOfCharge(obj)
            stateOfCharge = obj.stateOfCharge;
        end

        function setStateOfCharge(obj, stateOfCharge)
            obj.stateOfCharge = stateOfCharge;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function imResult = blendMode_Glow(A, B, offsetW, offsetH)
%% Glow blending mode: This mode is a variation of reflect mode (base and
%   blend color swapped).
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_Glow));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

C = blendMode_Reflect(B, A);

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end

author: Sable

repo: mcbench-benchmarks


text:
function a = conex2_inverse ( alpha )

%*****************************************************************************80
%
%% CONEX2_INVERSE returns the inverse of the CONEX2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar defining A.
%    A common value is 100.0.  ALPHA must not be zero.
%
%    Output, real A(3,3), the matrix.
%
  a = zeros ( 3, 3 );

  if ( alpha == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CONEX2_INVERSE - Fatal error!\n' );
    fprintf ( 1, '  The input value of ALPHA was zero.\n' );
    error ( 'CONEX2_INVERSE - Fatal error!' );
  end

  a(1,1) = 1.0;
  a(1,2) = ( 1.0 - alpha * alpha ) / alpha;
  a(1,3) = ( 1.0 + alpha * alpha ) / alpha^2;

  a(2,1) = 0.0;
  a(2,2) = alpha;
  a(2,3) = 1.0;

  a(3,1) = 0.0;
  a(3,2) = 0.0;
  a(3,3) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function line_cvt_lloyd_test04 ( )

%*****************************************************************************80
%
%% LINE_CVT_LLOYD_TEST04 tests the constrained computation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_CVT_LLOYD_TEST04:\n' );
  fprintf ( 1, '  Test the constrained computation.\n' );
  fprintf ( 1, '  SORT the initial points before use.\n' );

  n = 25;
  a = 0.0;
  b = 1.0;
  it_num = 200;
  x = a + ( b - a ) * rand ( n, 1 );
  x = sort ( x );
  header = 'test04';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use %d points in the interval [%g,%g]\n', n, a, b );
  fprintf ( 1, '  Take %d iterations.\n', it_num );
  fprintf ( 1, '  Call this calculation "%s"\n', header );
  fprintf ( 1, '  Expect a uniform spacing of %g\n', ( b - a ) / n );

  x = line_ccvt_lloyd ( n, a, b, it_num, header, x );

  return
end


text:
function img=ssimReadImg(imgName,parms)
% Read an image, and converts it into appropriate
% channels depending on the parameters

img=imread(imgName);
if parms.nChannels==3
  if(size(img,3)==1)
    img(:,:,2)=img(:,:,1);img(:,:,3)=img(:,:,1);
  end;
elseif parms.nChannels==1
  if(size(img,3)==3)
    img=rgb2gray(img);
  end;
else
  error('Invalid number of image channels\n');
end

img=double(img);

author: adikhosla

repo: feature-extraction


text:
function s=get_name(a,detailed)
 s=get_name(a.algorithm);
 eval_name

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef SomeOptimalSuperEllipseExponentVsStress < handle

    properties (Access = private)
        samplePoints
        qOpt
        phiMin
        phiMax
        nPhi
        rhoV
        xiV
        rho
        xi
        outputPath
        phiDV
        phiV
        iIn
        iOut
        qMin
        qMax
        qMean
    end

    methods (Access = public)

        function obj = SomeOptimalSuperEllipseExponentVsStress()
            obj.init();
            for iTest = 1:length(obj.xiV)
                obj.rho   = obj.rhoV(iTest);
                obj.xi   = obj.xiV(iTest);
                obj.phiV  = obj.createPhi();
                obj.phiDV = obj.createPhiInInterval();
                obj.createSamplePoints();
                obj.computeOptimalSuperEllipseExponent();
                obj.computeMeanSuperEllipseExponent();
                obj.plotQoptAndGaussianVsPhi(iTest);
            end
        end

    end

    methods (Access = private)

        function init(obj)
            obj.outputPath = '/home/alex/Dropbox/PaperStress/';
            obj.phiMin = 0;
            obj.phiMax = pi;
            obj.nPhi = 3;%50;
            obj.rhoV  = [0.9,0.9,0.5,0.5];
            obj.xiV  = pi/2 - [0.1083,0.557,0.88974,1.0984];
        end

        function itIs = isInInterval(obj,phi)
            itIs = phi>= obj.phiMin & phi <= obj.phiMax;
        end

        function phiD = createPhiInInterval(obj)
            phi = obj.phiV;
            iI  = obj.isInInterval(phi);
            iO  = ~iI;
            phiD(iI) = phi(iI);
            phiD(iO) = sign(phi(iO)).*(phi(iO)-pi/2);
        end

        function phi = createPhi(obj)
            phi = linspace(obj.phiMin,obj.phiMax,obj.nPhi);
        end

        function createSamplePoints(obj)
            s.type = 'FromFixedRhoAndTxi';
            s.rho0 = obj.rho;
            s.txi  = obj.xi;
            s.phi  = obj.phiV;
            sample = SamplePointsCreatorForOptimalExponentComputer.create(s);
            sample.compute();
            obj.samplePoints = sample;
        end

        function computeOptimalSuperEllipseExponent(obj)
            s.samplePoints = obj.samplePoints;
            rhoT = strrep(num2str(obj.rho),'.','_');
            txiT = strrep(num2str(round(obj.xi,3)),'.','_');
            fN = ['AveragingSuperEllipseRho',rhoT,'Txi',txiT];
            s.fileName = fN;
            exponentComputer = OptimalExponentComputer(s);
            exponentComputer.compute();
            obj.qOpt = exponentComputer.qOpt;
            obj.qMax = exponentComputer.qMax;
            obj.qMin = exponentComputer.qMin;
        end

        function computeMeanSuperEllipseExponent(obj)
            s.phiV = obj.phiV;
            qMeanC    = SuperEllipseMeanAndDesvExponentComputer(s);
            obj.qMean = qMeanC.computeMean(obj.qOpt,obj.xi);
        end

        function plotQoptAndGaussianVsPhi(obj,itxi)
            s.phiV = obj.phiV;
            s.phiMin = obj.phiMin;
            s.phiMax = obj.phiMax;
            s.qOpt = obj.qOpt;
            s.qMin = obj.qMin;
            s.qMax = obj.qMax;
            s.qMean = obj.qMean;
            s.xi = obj.xi;
            s.rho = obj.rho;
            p = OptimalSuperEllipseExponentVsGaussianPlotter(s);
            p.plot(itxi);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%%% Setup
% restoredefaultpath;
close, clear, clc;

%%% Path
addpath(fullfile(pwd,'utils'));

%%% Params
params.sequences_path = fullfile(pwd,'sequences');
params.sequences_mat = fullfile(pwd,'sequences_mat');
params.sequences_name = get_sequences_name(params);
params.sequences_format = 'PNG';
params.sequences_ext = strcat('*.',params.sequences_format);
params.debug = 0;
%params.current_sequence = 1;
clc;

%% Frame selection
displog('Performing frame selection');
for seq = 1:size(params.sequences_name,2)
    %seq = 1;
    displog(['Processing sequence: ' params.sequences_name(seq).name]);
    params.current_sequence = seq;
    perform_frame_selection(params);
    %break;
end
clear seq;

%% Matrix Completion
clc;
params.algs_path = 'algs_mc';
params.algs_name = get_algs_name(params);
params.results_path = 'results_mc';
displog('--- Matrix Completion ---');
for alg = 1:size(params.algs_name,2)
  %alg = 1;
  displog(['Current algorithm: ' params.algs_name(alg).name]);
  params.current_algorithm = alg;
  for seq = 1:size(params.sequences_name,2)
    %seq = 1;
    displog(['Processing sequence: ' params.sequences_name(seq).name]);
    params.current_sequence = seq;
    perform_matrix_completion(params);
    %break;
  end
  %break;
end
clear alg seq;

%% Tensor Completion
clc;
params.algs_path = 'algs_tc';
params.algs_name = get_algs_name(params);
params.results_path = 'results_tc';
displog('--- Tensor Completion ---');
for alg = 1:size(params.algs_name,2)
  %alg = 1;
  displog(['Current algorithm: ' params.algs_name(alg).name]);
  params.current_algorithm = alg;
  for seq = 1:size(params.sequences_name,2)
    %seq = 1;
    displog(['Processing sequence: ' params.sequences_name(seq).name]);
    params.current_sequence = seq;
    perform_tensor_completion(params);
    %break;
  end
  %break;
end
clear alg seq;

author: andrewssobral


text:
function input = prepare_data_for_input(root, list_pos, box_pos, pts_pos, list_neg, box_neg, ignore_neg, param)
    input = cell(length(param.gpu_id), 1);
    input(1:length(param.pos_gpu)) = prepare_pos_data_for_input( root, list_pos, box_pos, pts_pos, param );
    input(length(param.pos_gpu)+1:end) = prepare_neg_data_for_input( root, list_neg, box_neg, ignore_neg, param );
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function realimag(c)
%REALIMAG     Display real and imaginary part of interval hessians separately
%
%   realimag(c)
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 02/11/06     S.M. Rump  SparseInfNanFlag removed
%

  loose = strcmp(get(0,'FormatSpacing'),'loose');

  name = inputname(1);
  if isempty(name)                    % happens for display(hessianinit(random))
    name = 'ans';
  end

  if isreal(c.x)
    display(c,name)
  else
    if loose, disp(' '); end
    display(real(c),['real(' name ')'])
    if loose, disp(' '); end

    if loose, disp(' '); end
    display(imag(c),['imag(' name ')'])
    if loose, disp(' '); end
  end

repo: OpenMAS


text:
% Some Utility Functions used in the Library
%
% Description
% Including hertz to angular coversion,
% frame_by_frame_calculation used for computer conputation,
% and inner() to calculate inner product of two matrix.
%
% Methods
% angular = hertz_to_angular(hertz, samp_rate)
% hertz = angular_to_hertz(angle, samp_rate)
% coeffs = frame_by_frame_calculation(computer, signal, chunk_size)
% y = inner(a,b)
%
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef Util < handle

    methods (Static)
        function angular = hertz_to_angular(hertz, samp_rate)
            % Convert cycles/sec to radians/sec
            %
            % Input
            % samp_rate : [Hz] sampling rate in Hz
            %
            % Output
            % hertz : [radian] sampling rate in radian
            %
            % Example
            % >> samp_rate = 16
            % >> Angular = hertz_to_angular(samp_rate);
            %
            angular = hertz * 2 * pi / samp_rate;
        end

        function hertz = angular_to_hertz(angle, samp_rate)
            % Convert radians/sec to cycles/sec
            %
            % Input
            % samp_rate : [radian] sampling rate in radian
            %
            % Output
            % hertz : [Hz] sampling rate in Hz
            %
            % Example
            % >> samp_rate = 16
            % >> Hertz = angular_to_hertz(samp_rate);
            %
            hertz = angle * samp_rate / (2 * pi);
        end

        function coeffs = frame_by_frame_calculation(computer, signal, chunk_size)
            % Compute feature representation of entire signal iteratively
            %
            % Description
            % This function constructs a feature matrix of a signal through
            % successive calls to `computer.compute_chunk`. Its return value
            % should be identical to that of calling
            % `computer.compute_full(signal)`, but is possibly much slower.
            % `computer.compute_full` should be favoured.
            %
            % Inputs
            % computer   : Computer Object
            % signal     : array-like,  A 1D float array of the entire signal
            % chunk_size : int
            %              Length of the signal buffer to process at a given time
            %
            % Outputs
            % coeffs :  A 2D float array of shape ``(num_frames, num_coeffs)``.
            % ``num_frames`` is nonnegative (possibly 0). Contains some number
            % of feature vectors, ordered in time over axis 0.
            %
            % Example
            % >> bank = GaborFilterBank(MelScaling());
            % >> frame_length_ms = 25;
            % >> frame_shift_ms = 10;
            % >> frame_style = 'causal';
            % >> include_energy = true;
            % >> pad_to_nearest_power_of_two = true;
            % >> use_log = true;
            % >> use_power = true;
            % >> computer = ShortTimeFourierTransformFrameComputer(...
            %                                         bank, ...
            %                                         frame_length_ms, ...
            %                                         frame_shift_ms, ...
            %                                         frame_style, ...
            %                                         include_energy, ...
            %                                         pad_to_nearest_power_of_two, ...
            %                                         use_log, ...
            %                                         use_power);
            % feats_framewise = Util.frame_by_frame_calculation(computer, buff);
            %

            if nargin == 2
                chunk_size = 2 ^ 10;
            end
            if (computer.started == true)
                error('Already started computing frames');
            end
            coeffs = [];

            while length(signal) > 0
                coeff = computer.compute_chunk(signal(1:min(chunk_size, length(signal))));
                coeffs = [coeffs; ...
                            coeff];
                if chunk_size <= length(signal)
                    signal = signal(chunk_size+1:end);
                else
                    signal = [];
                end
            end

            chunk_finalize = computer.finalize();

            coeffs = [coeffs; chunk_finalize];
            coeffs = horzcat(coeffs);
        end

        function res = inner(a,b)
            % Compute the inner product of two vectors a and b.
            %
            % Inputs
            % a : vector
            % b : vector
            %   If a and b are nonscalar, their last dimensions must match.
            %
            % Outputs:
            % res : The value of the inner product of a and b.
            %       res.shape = a.shape[:-1] + b.shape[:-1]
            %
            % Example
            % >> y = inner(a,b) or inner(a,b)
            %

            c = 0;
            n = length(a);
            for k=1:n
                c=c+a(k)*b(k);
            end
            res = c;
        end

    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function varargout=var(varargin)
%VAR (overloaded)
%
% V = var(x)

x = varargin{1};

if nargin > 1 | min(size(x))>1
    error('SDPVAR/VAR only supports simple 1-D variance'),
end

switch length(x)
    case 1
        varargout{1} = x;
    otherwise
        x = reshape(x,length(x),1);
        m = sum(x)/length(x);
        varargout{1} = ((x-m)'*(x-m)) / (length(x) - 1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function out = Solve_Optimal(M,b1,b2,lambda)

[m, n]=size(b1);

smallNum = 0.000001;

a = abs(imfilter(b2,ones(7)/(sum(sum(ones(7)))),'replicate'));
A = 1./(a + smallNum);

A = A(:);
A = sparse(1:1:m*n,1:1:m*n,A(1:1:m*n)',m*n,m*n);

II=sparse(1:1:m*n,1:1:m*n,ones(1,m*n),m*n,m*n);

denominator=2*II+lambda*(A+A');

numerator=2*M(:)+lambda*(A+A')*b1(:);

out=denominator\numerator;
out=reshape(out,m,n);

end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [j, Grad] = costfunction(Theta, X, Y)

n = length(Y);

j = 0;

j = (1/n)*(-Y'* log(Sigmoid(X * Theta)) - (1 - Y)'* log(1 - Sigmoid(X * Theta))); % See the sigmoid function.

Grad = zeros(size(Theta));

Grad_ = (1 / n) * ((Sigmoid(X * Theta) - Y)' * X);

Grad = Grad_';

end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function c=ref_dgt_3(f,g,a,M)
%REF_DGT_3  DGT algorithm 3

L=size(g,1);
N=L/a;
b=L/M;

[c,h_a,h_m]=gcd(-a,M);
p=a/c;
q=M/c;
d=N/q;

w=zeros(M,N);

if 0

  % This version uses the definition.

  F=zeros(c,d,p,q);
  G=zeros(c,d,p,q);

  for r=0:c-1
    for s=0:d-1
      for k=0:p-1
	for l=0:q-1
	  for st=0:d-1
	    F(r+1,s+1,k+1,l+1)=F(r+1,s+1,k+1,l+1)+f(mod(r+k*M+st*p*M-l*h_a*a,L)+1)*exp(-2*pi*i*s*st/d);
	    G(r+1,s+1,k+1,l+1)=G(r+1,s+1,k+1,l+1)+g(mod(r+k*M-l*a+st*p*M,L)+1)*exp(-2*pi*i*s*st/d);
	  end;
	end;
      end;
    end;
  end;

 for r=0:c-1
    for l=0:q-1
      for u=0:q-1
	for s=0:d-1
	  for v=0:d-1
	    for k=0:p-1
	      w(r+l*c+1,mod(u+s*q-l*h_a,N)+1)=w(r+l*c+1,mod(u+s*q-l*h_a,N)+1)+...
		  1/d*F(r+1,v+1,k+1,l+1)*conj(G(r+1,v+1,k+1,u+1))*exp(2*pi*i*v*s/d);
	    end;
	  end;
	end;
      end;
    end;
  end;

else

  % This version uses matrix-vector products and ffts

  F=zeros(c,d,p,q);
  G=zeros(c,d,p,q);
  C=zeros(c,d,q,q);

  % Set up the matrices
  for r=0:c-1
    for s=0:d-1
      for k=0:p-1
	for l=0:q-1
	  F(r+1,s+1,k+1,l+1)=f(mod(r+k*M+s*p*M-l*h_a*a,L)+1);
	  G(r+1,s+1,k+1,l+1)=sqrt(M*d)*g(mod(r+k*M-l*a+s*p*M,L)+1);
	end;
      end;
    end;
  end;

  % fft them
  F=dft(F,[],2);
  G=dft(G,[],2);

  % Multiply them
  for r=0:c-1
    for s=0:d-1
      GM=reshape(G(r+1,s+1,:,:),p,q);
      FM=reshape(F(r+1,s+1,:,:),p,q);
      C(r+1,s+1,:,:)=GM'*FM;
    end;
  end;

  % Inverse fft
  C=idft(C,[],2);

  % Place the result
  for r=0:c-1
    for l=0:q-1
      for u=0:q-1
	for s=0:d-1
	  w(r+l*c+1,mod(u+s*q-l*h_a,N)+1)=C(r+1,s+1,u+1,l+1);
	end;
      end;
    end;
  end;

end;

c=dft(w);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function r = GetRadar(dt)
%
%
persistent posp

if isempty(posp)
  posp = 0;
end

vel = 100  +  5*randn;
alt = 1000 + 10*randn;

pos = posp + vel*dt;

v = 0 + pos*0.05*randn;

r = sqrt(pos^2 + alt^2) + v;

posp = pos;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function [eu,EUhm] = hmg2euc(hm)

% HMG2EUC Homogeneous to Euclidean point transform.
%   UMG2EUC(HM) is the Euclidean point corresponding to the homogeneous
%   point HM.
%
%   [e, E_hm] = HMG2EUC(HM) returns the Jacobian of the transformation.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

s = size(hm,1);
n = size(hm,2);

eu = hm(1:s-1,:)./repmat(hm(s,:),s-1,1);

if nargout > 1 % Jacobians
    if n == 1
        EUhm = [1/hm(s)*eye(s-1) -eu/hm(s)]; % this takes any dimension
    else
        error('Jacobians not computed for multiple points in HM')
    end
end

return

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
for i = 60:63
	n = 2^i;
	p = [0.4; 0.6];
	x = sample_hist(p, n);
	assert(sum(x) == n);
	assert(max(abs(x/n - p)) < 1e-4);
end

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed


text:
function [a,t,tLength,k]=BSplineInterpIntMultiDim(t,tLength,a,k,intDims)
%%BSPLINEINTERPINTMULTIDIM Given a hypermatrix of multivariate b-spline
%            interpolation weights a and a set of knots t,compute the
%            modified weights and knots to interpolate the integral of the
%            a given point using b-spline interpolation. If all knots at
%            the ith boundary are repeated k(i) times, then an integral
%            over the ith dimension will start from the beginning of the
%            interpolation region. Otherwise, this can be considered an
%            indefinite integral with a particular additive constant and
%            differencing and be used to get a definite integral.
%
%INPUTS: t The maxNumKnotsXnumDims set of knots fo r the interpolation
%          function. These are actually the values in each coordinate
%          dimension given in ascending order. The full knots are implied
%          [t1,t2,etc.]=ndgrid(t(1:tLengths(1),1),t(1:tLengths(2),2),...)
%          and the dimensions can be put together into the full set of
%          numDimsXtotalNumPoints points as tTotal=[t1(:).';t2(:).';...].
%          The first and last k(curDim)-1 knots in each dimension are
%          outside of the ends of the region with data or mark the ends of
%          the region with data. The number of rows is the maximum number
%          needed for all of the dimensions. tLength says how many items
%          in each column are actually used. These values must be real.
%  tLength A numDimX1 vector where tLength(i) says the number of elements
%          in t(:,i).
%        a A hypermatrix with numDim+1 indices containing the set of
%          coefficients for the b-splines covering all of the interpolation
%          intervals for all of the sets of interpolation values. If there
%          is only one set, the final dimension is unitary meaning that
%          there are effectively only numDim dimensions. The derivatives
%          will be applied to all sets present. These values can be real or
%          complex.
%        k A numDimsX1 or 1XnumDims set of the order of the b-splines in
%          each dimension. The value k-1 is the polynomial order of the
%          approximation being performed. If the order is the same in all
%          dimensions, then a single scalar can be passed.
%  intDims This is a numDimsX1 or a 1XnumDims boolean vector indicating
%          over which dimensions the integrals should be taken.
%
%OUTPUTS: a The modified coefficient matrix. This is one larger in each
%           dimensions where an integral was taken.
%         t The modified knots.
%   tLength The modified vector indicating how many knows are present.
%         k The modified set of orders.
%
%Equation 22 in Chapter X of [1] has 1D integration. Chapter XVII describes
%the use of tensor product splines for more than one dimension. The
%generalization to more than one dimension comes from just evaluating the
%integrals across each dimension one at a time, if they are taken.
%
%EXAMPLE 1:
% f=@(x,y,z)(x.^4-2*x.^2+x).*(y.^4-2*y.^2+y).*(z.^4-3*z.^2+z);
% %The definite integral of f with respect to x and z starting at -1.5 is
% fIntXZ=@(x,y,z)((-891+240*x.^2-320*x.^3+96*x.^5).*y.*(1-2*y+y.^3).*(-477+80*z.^2-160*z.^3+32*z.^5))/76800;
%
% numDims=3;
% numPointsX=10;
% numPointsY=11;
% numPointsZ=9;
% tauLengths=[numPointsX;numPointsY;numPointsZ];
% tau=zeros(max(tauLengths),numDims);
% tau(1:tauLengths(1),1)=linspace(-1.5,1.5,numPointsX);
% tau(1:tauLengths(2),2)=linspace(-1.5,1.5,numPointsY);
% tau(1:tauLengths(3),3)=linspace(-1.5,1.5,numPointsZ);
% %Note that using meshgrid would have put the elements in the wong order.
% %ndgrid must be used.
% [tau1,tau2,tau3]=ndgrid(tau(1:tauLengths(1),1),tau(1:tauLengths(2),2),tau(1:tauLengths(3),3));
% y=f(tau1,tau2,tau3);
%
% k=[5;5;5];
% [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,y(:),k);
% intDims=[1;0;1];
% [a,t,tLength,k]=BSplineInterpIntMultiDim(t,tLength,a,k,intDims);
%
% numPointsY=21;
% numPointsZ=19;
% pointsX=linspace(-1.5,1.5,numPointsX);
% pointsY=linspace(-1.5,1.5,numPointsY);
% pointsZ=linspace(-1.5,1.5,numPointsZ);
% [X,Y,Z]=ndgrid(pointsX,pointsY,pointsZ);
% x=[X(:).';Y(:).';Z(:).'];
%
% zIntTrue=fIntXZ(X,Y,Z);
% z=BSplineInterpValMultiDim(x,t,tLength,a,k);
% max(abs(z(:)-zIntTrue(:)))
%
%EXAMPLE 2:
%This is similar to example 1, except two sets of values are fitted,
%integrated, and interpolated at once.
% f1=@(x,y,z)(x.^4-2*x.^2+x).*(y.^4-2*y.^2+y).*(z.^4-3*z.^2+z);
% f2=@(x,y,z)(x.^3-2*x.^2+x).*(y.^4-2*y.^2+y).*(2*z.^3-2*z.^2+1);
% %The definite integral of f with respect to x and z starting at -1.5 is
% fIntXZ1=@(x,y,z)((-891+240*x.^2-320*x.^3+96*x.^5).*y.*(1-2*y+y.^3).*(-477+80*z.^2-160*z.^3+32*z.^5))/76800;
% fIntXZ2=@(x,y,z)(((-891+16*x.^2.*(6+x.*(-8+3*x))).*y.*(1-2*y+y.^3).*(-315+16*z.*(6+z.^2.*(-4+3*z))))/18432);
%
% numDims=3;
% numPointsX=10;
% numPointsY=11;
% numPointsZ=9;
% tauLengths=[numPointsX;numPointsY;numPointsZ];
% tau=zeros(max(tauLengths),numDims);
% tau(1:tauLengths(1),1)=linspace(-1.5,1.5,numPointsX);
% tau(1:tauLengths(2),2)=linspace(-1.5,1.5,numPointsY);
% tau(1:tauLengths(3),3)=linspace(-1.5,1.5,numPointsZ);
% %Note that using meshgrid would have put the elements in the wong order.
% %ndgrid must be used.
% [tau1,tau2,tau3]=ndgrid(tau(1:tauLengths(1),1),tau(1:tauLengths(2),2),tau(1:tauLengths(3),3));
% y1=f1(tau1,tau2,tau3);
% y2=f2(tau1,tau2,tau3);
% y=[y1(:),y2(:)];
%
% k=[5;5;5];
% [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,y,k);
% intDims=[1;0;1];
% [a,t,tLength,k]=BSplineInterpIntMultiDim(t,tLength,a,k,intDims);
%
% numPointsY=21;
% numPointsZ=19;
% pointsX=linspace(-1.5,1.5,numPointsX);
% pointsY=linspace(-1.5,1.5,numPointsY);
% pointsZ=linspace(-1.5,1.5,numPointsZ);
% [X,Y,Z]=ndgrid(pointsX,pointsY,pointsZ);
% x=[X(:).';Y(:).';Z(:).'];
%
% zIntTrue1=fIntXZ1(X,Y,Z);
% zIntTrue2=fIntXZ2(X,Y,Z);
% z=BSplineInterpValMultiDim(x,t,tLength,a,k);
% max(abs(z(:,1)-zIntTrue1(:)))
% max(abs(z(:,2)-zIntTrue2(:)))
%
%EXAMPLE 3:
%This is an example of a bivariate complex function whose integral over the
%first dimension is interpolated.numPoints=80;
% xMin=0;
% xMax=3;
% yMin=0;
% yMax=3;
% x=linspace(xMin,xMax,numPoints);
% y=linspace(yMin,yMax,numPoints);
% [X,Y]=ndgrid(x,y);
% f=@(x,y)(sin(5*y)+1j*cos(10*x));
% %The integral from 0 in the x dimension.
% fIntx=@(x,y)((1/10)*1j*sin(10*x)+x.*sin(5*y));
%
% fxy=f(X(:),Y(:));
% k=[5;5];
% tau=[x(:),y(:)];
% tauLengths=[numPoints,numPoints];
% [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,fxy,k);
% intDims=[1;0];
% [a,t,tLength,k]=BSplineInterpIntMultiDim(t,tLength,a,k,intDims);
%
% %Interpolate the derivative
% numPoints=100;
% x=linspace(xMin,xMax,numPoints);
% y=linspace(xMin,xMax,numPoints);
% [X,Y]=ndgrid(x,y);
%
% fxy=fIntx(X(:),Y(:));
% pts=[X(:).';Y(:).'];
% fxyInterp=reshape(BSplineInterpValMultiDim(pts,t,tLength,a,k),[numPoints,numPoints]);
% fxy=reshape(fxy,[numPoints,numPoints]);
%
% figure(1)
% clf
% hold on
% surface(X,Y,real(fxy),'EdgeColor','none')
% title('Real Integral Values')
% colorbar()
%
% figure(2)
% clf
% hold on
% surface(X,Y,imag(fxy),'EdgeColor','none')
% title('Imaginary Integral Values')
% colorbar()
%
% figure(3)
% clf
% hold on
% surface(X,Y,abs(real(fxyInterp)-real(fxy)),'EdgeColor','none')
% title('Real Interpolation Error')
% colorbar()
%
% figure(4)
% clf
% hold on
% surface(X,Y,abs(imag(fxyInterp)-imag(fxy)),'EdgeColor','none')
% title('Imaginary Interpolation Error')
% colorbar()
%
%REFERENCES:
%[1] C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag,
%    1978.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDims=length(tLength);
numA=size(a);

if(length(numA)==numDims)
    numSets=1;
    numA=[numA,1];
else
    numSets=numA(end);
end

tNew=zeros(max(tLength)+2,numDims);
for curDim=1:numDims
    if(intDims(curDim)==0)
        tNew(1:tLength(curDim),curDim)=t(1:tLength(curDim),curDim);
        continue;
    end

    temp=zeros(1,numDims);
    temp(curDim)=1;
    aNew=zeros(numA+[temp,0]);

    %Next, we go through all tuples of values for the dimensions other than
    %this one. In each instance, we must take the integral of the a terms
    %that are present. These indices select all of the dimensions that are
    %not the current one across which derivatives are being taken.
    for curSet=1:numSets
        idxList=[1:(curDim-1),(curDim+1):numDims];

        maxVals=numA(idxList)-1;
        numACur=numA(curDim);
        curTuple=getNextTuple(numDims-1);
        idxCell=cell(1,numDims+1);
        idxCell{numDims+1}=curSet;
        tCur=t(1:tLength(curDim),curDim);
        kCur=k(curDim);

        outputIdx=1:(numACur+1);
        while(~isempty(curTuple))
            idxCell{curDim}=1:numA(curDim);%Mark the free dimension.

            %We have to select the elements in a to work on based on the
            %current tuple.
            for curIdx=1:(numDims-1)
                idxCell{idxList(curIdx)}=curTuple(curIdx)+1;
            end

            aCur=reshape(a(idxCell{:}),[numACur,1]);

            [~,aCur]=BSplineInterpInt(tCur,aCur,kCur);
            idxCell{curDim}=outputIdx;
            aNew(idxCell{:})=aCur;

            curTuple=getNextTuple(curTuple,maxVals);
        end
    end
    numA(curDim)=numA(curDim)+1;
    a=aNew;
    k(curDim)=k(curDim)+1;
    tNew(1:(tLength(curDim)+2),curDim)=[t(1,curDim);t(1:tLength(curDim),curDim);t(tLength(curDim),curDim)];
    tLength(curDim)=tLength(curDim)+2;
end
a=aNew;
t=tNew;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [fMinLimit, fMaxLimit, nTickStep] = ex_CalcFigureLimits(fMin, fMax)

fDiff = fMax - fMin;
nDigits = ceil(log10(fDiff));
nTickStep = 10^(nDigits-1);

bMinPositive = (fMin >= 0);
if ~bMinPositive
  fMin = (-1) * fMin;
end
bMaxPositive = (fMax >= 0);
if ~bMaxPositive
  fMax = (-1) * fMax;
end

if (fMin == 0)
  fMinLimit = 0;
else
  if bMinPositive
    fMinLimit = (fMin - 10^(nDigits-1));
    fMinLimit = ceil(fMinLimit * 10^(-nDigits+1))/10^(-nDigits+1);
  else
    fMinLimit = (fMin + 10^(nDigits-1));
    fMinLimit = floor(fMinLimit * 10^(-nDigits+1))/10^(-nDigits+1);
  end
  if ~bMinPositive
    fMinLimit = (-1) * fMinLimit;
  end
end
if (fMax == 0)
  fMaxLimit = 0;
else
  fMaxLimit = (fMax + 10^(nDigits-1));
  fMaxLimit = floor(fMaxLimit * 10^(-nDigits+1))/10^(-nDigits+1);
  if ~bMaxPositive
    fMaxLimit = (-1) * fMaxLimit;
  end
end

author: CelsoReyes


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, sp, u, gnum] = ...
              mp_solve_laplace_2d (problem_data, method_data)

warning ('geopdes:obsolete', 'Function MP_SOLVE_LAPLACE_2D is obsolete. Using MP_SOLVE_LAPLACE instead')

[geometry, msh, sp, u, gnum] = mp_solve_laplace (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function maxfig(varargin)
% maxfig(fig,extend,mode)  MAXimize a FIGure to fill the screen.
%   Subsequent calls will restore the maximized figure to its previous state.
%   Examples:
%       maxfig, maxfig(3,1+i), maxfig([],-1), maxfig('normal')
%   "fig" (first numeric input) is the handle of the figure to maximize
%       (default is current figure).
%   "extend" (second numeric input) is a complex number specifying how
%       figures are to be created on the extended desktop.
%       Examples: +1/-1 is one screen to the right/left,
%                 0 is the primary screen (default),
%                 +i is one screen up,
%                 -1-0.5i is half a screen down and to the left, etc...
%   "mode" is a string specifying whether to fully maximize the figure
%       (borders extend beyond the monitor, menu and tool bars are removed)
%       or set to normal mode with borders, menu & toolbars displayed.
%       Value for maximize mode is any of {m,max,y,yes,t,true}(default).
%       Value for normal mode is any of {n,normal,no,f,false}.

%   Copyright 2006 Mirtech, Inc.
%   created 08/23/2006  by Mirko Hrovat on Matlab Ver. 7.2
%   Mirtech, Inc.   email: mhrovat@email.com
% More notes:
% Will work with docked figures.
% Inspired by MAXFIGSIZE (ID 3036) by Duane Hanselman.
% Extended monitors must have the same resolution as the primary monitor.
% The figure properties, TAG and USERDATA, are used to keep track of the
% maximized figure.

% --------- Key Constants -----------
% To allow for user or OS borders set the appropriate value in pixels.
% For example one might want to set the bottomborder for the windows
% taskbar (try a value of 40).
topborder       =0;
bottomborder    =0;
rightborder     =0;
leftborder      =0;
% Border properties for figures - Adjust if needed.
% determined by the difference between OuterPosition and Position values.
% GETFIGDIM can return these values for your system.
figureborder    =5;
titleheight     =26;        %#ok
% ----------------------------------

fullflag=1;     % set default state for maximization choice
extend=[];
% just in case, check all figures for 'MaxFig' tag
set(0,'ShowHiddenHandles','on')
fig= findobj('Tag','MaxFig');
set(0,'ShowHiddenHandles','off')
% If "MaxFig" figure is found then reset size to original values
if ~isempty(fig),
    set(fig,'Tag','')
    userdata=get(fig,'UserData');
    set(fig,'Units',userdata{1},'OuterPosition',userdata{2});
    set(fig,'MenuBar',userdata{3},'ToolBar',userdata{4},...
        'WindowStyle',userdata{5});
    drawnow
end
if nargin==0,
    if ~isempty(fig),
        return
    % return if no arguments, maxfig is only resetting the figure.
    end
else
    count=0;
    for n=1:length(varargin),
        arg=varargin{n};
        if isnumeric(arg),      % note "empty" is also considered to be numeric!
            count=count+1;
            switch count
                case 1
                    fig=arg;
                case 2
                    extend=arg;
            end  % switch count
        elseif ischar(arg),
            if any(strcmpi(arg,{'y','yes','t','true','m','max'})),
                fullflag=1;
            elseif any(strcmpi(arg,{'n','no','f','false','normal'})),
                fullflag=0;
            end
        end
    end  % for
end  % if nargin
if isempty(fig),        fig=gcf;    end
if isempty(extend),     extend=0;   end

% Check WindowStyle property
winstyle=get(fig,'WindowStyle');
if strcmpi(winstyle,'modal'),
    error('  Cannot maximize this figure.')
end

% Get screen dimensions. Note that multiple monitors are assumed to have
% the same resolution.
set(0,'Units','pixels');
scnsize=get(0,'ScreenSize');
scnwidth=scnsize(3)-leftborder-rightborder;
scnheight=scnsize(4)-topborder-bottomborder;

% Calculate offsets for multiple monitors.
offx=scnsize(3)*real(extend)+leftborder;
offy=scnsize(4)*imag(extend)+bottomborder;

% Save properties for figure in UserData field
set(fig,'Tag','MaxFig')
menbar=get(fig,'MenuBar');
tbar=get(fig,'ToolBar');
set(fig,'UserData',{get(fig,'Units'), get(fig,'OuterPosition'),...
    menbar,tbar,winstyle})

if fullflag
    tbar='none';
    menbar='none';
    offx=offx-figureborder;
    offy=offy-figureborder;
    scnwidth=scnwidth+2*figureborder;
    scnheight=scnheight+2*figureborder;
%    scnheight=scnheight+2*figureborder+titleheight;
% The above seems to create problems as results don't match settings.
%   The values seem to wrap around. So it looks like there is no simple
%   way to get rid of the title bar in Matlab.
end

pos=[offx+1,offy+1,scnwidth,scnheight];
% Move, size figure and modify properties as needed.
if strcmpi(winstyle,'docked'),
    set(fig,'WindowStyle','normal')
    drawnow
end
set(fig,'Units','pixels', ...
        'OuterPosition',pos,...
        'MenuBar',menbar,...
        'ToolBar',tbar);
figure(fig)                 % bring figure to front.
end

author: Sable

repo: mcbench-benchmarks


text:
function [A,B,C,D]=era(h,n,N,Ts,def);

% Eigensystem Realization Algorithm (ERA)
%
% Author: Samuel da Silva - UNICAMP
% e-mail: samsilva@fem.unicamp.br
% Date: 2006/10/20
%
% [A,B,C,D]=era(h,n,N,Ts,def);
%
% Inputs:
%    h: discrete-time impulse response
%    n: order of the system
%    N: number of samples to assembly the Hankel matrix
%    Ts: sample time
%    def: if = 1: the output will be the discrete-time state-space model
%         if = 2: the output will be the continuous-time state-space model
%
% Otputs:
%    [A,B,C,D]: state-space model
%
% Note: For now, it works to SISO systems and it is necessary the control toolbox
%
% References: Juang, J. N. and Phan, M. Q. "Identification and Control of
% Mechanical Systems", Cambridge University Press, 2001

% Hankel matrix
H0 = hankel(h(2:N+1));            % k = 0
H1 = hankel(h(3:N+2));            % k = 1;

% Factorization of the Hankel matrix by use of SVD
[R,Sigma,S] = svd(H0);
% R and S are orthonormal and Sigma is a rectangular matrix

Sigman = Sigma(1:n,1:n);

Wo = R(:,1:n)*Sigman^0.5;           % observability matrix
Co = Sigman^.5*S(:,1:n)';           % controllability matrix

% The identified system matrix are:
A = Sigman^-.5*R(:,1:n)'*H1*S(:,1:n)*Sigman^-.5;            % dynamic matrix
B = Co(:,1);                    % input matrix
C = Wo(1,:);                    % output matrix
D = h(1);                       % direct-transmission matrix

sysdisc = ss(A,B,C,D,Ts);       % discrete-time system

if def == 2
    syscont = d2c(sysdisc,'zoh');       % Conversion of discrete LTI models to continuous time
    [A,B,C,D]=ssdata(syscont);          % continuous system
end

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


text:
function callback_tracker(mysrc,~,fn)
    % provide message showing where something is called from
    cb = mysrc.Callback;
    if iscell(cb)
        cb=char(cb{1});
    end
    if isa(cb,'function_handle')
        cb=func2str(cb);
    end
    fprintf('in: %s called from %s\n',cb,fn);
end


text:
function [pE,gE,pC,gC] = spm_ind_priors(A,B,C,Nm,Nf)
% prior moments for a neural-mass model of induced responses
% FORMAT [pE,gE,pC,gC] = spm_ind_priors(A,B,C,dipfit,Nu,Nf)
% A{2},B{m},C  - binary constraints on extrinsic connections
% Nm           - number of frequency modes used
% Nf           - number of frequency modes explained
%
% pE - prior expectation - f(x,u,P,M)
% gE - prior expectation - g(x,u,G,M)
%
% connectivity parameters
%--------------------------------------------------------------------------
%    pE.A    - trial-invariant
%    pE.B{m} - trial-dependent
%    pE.C    - stimulus-stimulus dependent
%
% stimulus and noise parameters
%--------------------------------------------------------------------------
%    pE.R    - onset and dispersion
%
% pC - prior covariances: cov(spm_vec(pE))
%__________________________________________________________________________
%
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_ind_priors.m 5900 2014-02-27 21:54:51Z karl $

% orders
%--------------------------------------------------------------------------
if nargin < 5; Nf = Nm; end
n    = size(C,1);                                 % number of sources
nu   = size(C,2);                                 % number of inputs
D    = sparse(1:Nm,1:Nm,1,Nf,Nf);

% set extrinsic connectivity - linear and nonlinear (cross-frequency)
%--------------------------------------------------------------------------
E.A  = kron(speye(Nf,Nf), -speye(n,n));
V.A  = kron(D,A{1}) + kron(D*(1 - speye(Nf,Nf))*D,A{2});

% input-dependent
%--------------------------------------------------------------------------
for i = 1:length(B)
    E.B{i} = sparse(n*Nf,n*Nf);
    V.B{i} = kron(D*ones(Nf,Nf)*D,B{i}) & V.A;
end

% exogenous inputs
%--------------------------------------------------------------------------
E.C  = kron(D*ones(Nf,1),C - C);
V.C  = kron(D*ones(Nf,1),C);

% set stimulus parameters: magnitude, onset and dispersion
%--------------------------------------------------------------------------
E.R  = kron(ones(nu,1),[0 0]);
V.R  = kron(ones(nu,1),[1/16 1/16]);

% prior moments
%--------------------------------------------------------------------------
pE   = E;
gE   = [];
pC   = diag(sparse(spm_vec(V)));
gC   = [];

author: spm


text:
function [f, d] = getSIFTFeatures(image, edgeThresh)

%convert images to greyscale
if (size(image, 3) == 3)
    Im = single(rgb2gray(image));
else
    Im = single(image);
end

% get features and descriptors
[f, d] = vl_sift(Im, 'EdgeThresh', edgeThresh);

end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
classdef nnslice < nntest
  methods (Test)

    function basic(test)
        sz = [3,3,5,4] ;
        x = test.randn(sz) ;
        dim = 4 ;
        slicePoints = 1:dim - 1 ; % slice along fourth dim
        y = vl_nnslice(x, dim, slicePoints, []) ;

        % check derivatives with numerical approximation
        dzdy = cellfun(@(x) test.randn(size(x)), y, 'Uni', 0) ;
        dzdx = vl_nnslice(x, dim, slicePoints, dzdy, 'inputSizes', {sz}) ;
        dzdy_ = cat(dim, dzdy{:}) ;
        dzdx_ = dzdx{1} ;
        test.der(@(x) forward_wrapper(x, dim, slicePoints), x, dzdy_, dzdx_, 1e-3*test.range) ;
    end
  end
end

% -----------------------------------------------------------------
function y = forward_wrapper(x, dim, slicePoints)
% -----------------------------------------------------------------
  y = vl_nnslice(x, dim, slicePoints, []) ;
  y = cat(dim, y{:}) ;
end

repo: MFT


text:
function R = imnoise2(type,M,N,a,b)
% Imnoise2 to creates a random array with the specified PDF
% R = IMNOISE2(TYPE,M,N,A,B) generates an array, R, of size M by N whose
% elements are random numbers of the specified TYPE with parameters A and
% B. IF ony TYPE is included in the input argument list, a single random
% number of the specified TYPE and and default parameters show below is
% generated. If only TYPE, M and N are provided, the default parameters
% shown below are used. If M=N=1, IMNOISE2 generates a single random number
% of the specified TYPE and parameters A and B.
%
% Valid values for TYPE and parameters A and B are:
% 'uniform'        Uniform random numbers in the interval (A,B).
%                  The default values are (0,1).
%
% 'guassian'       Gaussian random numbers with mean A and standard
%                  deviation B.The default values are A = 0, B = 1.
%
% 'salt & pepper'  Salt and pepper numbers  of amplitude 1 with probability
%                  Pa  = A, and amplitude 0 with probability Pb = B. The
%                  default values are Pa = Pb = A = B = 0.05. Note that the
%                  noise has values 0(with probability Pa = A) and 1(with
%                  probability Pb = B), so scaling is necessary if values
%                  other than 0 and 1 are required. The noise matrix R is
%                  assigned three values. If R(x,y) = 0, the noise at (x,y)
%                  is pepper(black). If R(x,y ) = 1, the noise at (x,y) is
%                  salt(white). If R(x,y) = 0.5, there is no noise assigned
%                  to coordinates (x,y).
%
% 'lognormal'      Lognormal numbers with offset A and shape parameter B.
%                  The defaults are A = 1 and B = 0.25.
%
% 'rayleigh'       Rayleigh noise with parameters A and B. The default
%                  values are A = 0 and B = 1.
%
% 'exponential'    Exponential random numbers with parameter A. The default
%                  value is A = 1.
%
% 'erlang'         Erlang(gamma) random numbers with parameters A and B. B
%                  must be a positive integer. The defaults are A = 2 and B
%                  = 5. Erlang random numbers are approximated as the sum
%                  of B exponential random numbers.
%
% set default values.
if nargin == 1
    a = 0; b = 1;
    M = 1; N = 1;
elseif nargin == 3
    a = 0; b = 1;
end
%as we need only small letters as the type so...
switch lower(type)
    case 'uniform'
        R = a + (b-a)*rand(M,N);
    case 'gaussian'
        R = a + b*randn(M,N);
    case 'salt & pepper'
        if nargin <= 3
        a = 0.05; b = 0.05;
        end
% check to make sure that Pa + Pb is not > 1.
        if (a + b) > 1
            error('The sum of the Pa and Pb cannot exeed 1.')
        end
        R(1:M,1:N) = 0.5;
% Generate an M by N array of uniformly distributed random numbers in the
% range (0,1). Then, Pa*(M*N) of them will have values <= a. The
% coordinates of these points we call 0 (pepper noise). Similarly, Pb*(M*N)
% points will have values in the range > a & <= (a+b). These we call
% (salt noise).
        X = rand(M,N);
        c = find(X<=a);
        R(c) = 1;
        u = a + b;
        c = find(X > a & X <= u);
        R(c) = 1;
    case 'lognormal'
        if nargin<=3
            a = 1; b = 0.25;
        end
        R = a*esp(b*randn(M,N));
    case 'Rayleigh'
        R = a + (-b*log(1-rand(M,N)))^0.5;
    case 'exponential'
        if nargin <= 3
            a = 1;
        end
        if a <= 0
            error('the value of a must b positive for exponential operation')
        end
        k = -1/a
        R = k*log(1 - rand(M,N));
    case 'erlang'
        if nargin <= 3
            a = 2; b = 5;
        end
        if (b ~= round(b)| b <= 0)
            error('Parameter b should b a negative value for erlang')
        end
        k = -1/a;
        R = zeros(M,N);
        for j = 1:b
            R = R + k*log(1 - rand(M,N));
        end
    otherwise
        error('Unknown distribution type.')
end


text:
function varargout = transformVector(varargin)
%TRANSFORMVECTOR Transform a vector with an affine transform.
%
%   VECT2 = transformVector(VECT1, TRANS);
%   where VECT1 has the form [xv yv], and TRANS is a [2*2], [2*3] or [3*3]
%   matrix, returns the vector transformed with affine transform TRANS.
%
%   Format of TRANS can be one of :
%   [a b]   ,   [a b c] , or [a b c]
%   [d e]       [d e f]      [d e f]
%                            [0 0 1]
%
%   VECT2 = transformVector(VECT1, TRANS);
%   Also works when PTA is a [N*2] array of double. In this case, VECT2 has
%   the same size as VECT1.
%
%   [vx2 vy2] = transformVector(vx1, vy1, TRANS);
%   Also works when vx1 and vy1 are arrays the same size. The function
%   transform each couple of (vx1, vy1), and return the result in
%   (vx2, vy2), which is the same size as (vx1 vy1).
%
%
%   See also
%   vectors2d, transforms2d, rotateVector, transformPoint

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2007-03-12
% Copyright 2007-2022 INRA - TPV URPOI - BIA IMASTE

if length(varargin)==2
    var = varargin{1};
    vx = var(:,1);
    vy = var(:,2);
    trans = varargin{2};
elseif length(varargin)==3
    vx = varargin{1};
    vy = varargin{2};
    trans = varargin{3};
else
    error('wrong number of arguments in "transformVector"');
end

% compute new position of vector
vx2 = vx*trans(1,1) + vy*trans(1,2);
vy2 = vx*trans(2,1) + vy*trans(2,2);

% format output
if nargout==0 || nargout==1
    varargout{1} = [vx2 vy2];
elseif nargout==2
    varargout{1} = vx2;
    varargout{2} = vy2;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ff,f]=v_lpccc2ff(cc,np,nc,c0)
%V_LPCCC2FF Convert complex cepstrum to complex spectrum FF=(CC,NP,NC)
%
%  Inputs: cc(nf,n)     Complex ceptral coefficients excluding c(0), one frame per row
%          np           Size of output spectrum is np+1 [n]
%                       Alternatively, a vector of output frequencies in the range 0 to 0.5
%          nc           Number of cepstral coefficients to use [np or, if np is a vector, n]
%                       Set nc=-1 to use n coefficients
%          c0(nf,1)     Cepstral coefficient c(0) [0]
%
% Outputs: ff(nf,np+2)  Complex spectrum from DC to Nyquist
%          f(1,np+2)    Normalized frequencies (0 to 0.5)
%
% The "complex cepstral coefficients", cc(n), are the inverse discrete-time Fourier transform
% of the log of the complex-valued spectrum. The cc(n) are real-valued and, for n<0, cc(n)=0.
% The "real cepstral coeffcients", rc(n), are the inverse discrete-time Fourier transform
% of the log of the magnitude spectrum; rc(0)=cc(0) and rc(n)=0.5*cc(n) for n~=0.
% For highest speed, choose np+1 to be a power of 2.

%      Copyright (C) Mike Brookes 2014
%      Version: $Id: v_lpccc2ff.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,mc]=size(cc);
if nargin<2 || ~numel(np)
    if nargout
        np=mc;
    else
        np=128;
    end
end
if nargin>=3 && numel(nc)==1 && nc==-1 nc=mc; end
if nargin<4 || ~numel(c0) c0=zeros(nf,1); end
if numel(np)>1 || np(1)<1
    if nargin<3 || ~numel(nc) nc=mc; end
    f=np(:)';
    if nc==mc
        ff=exp([c0 cc]*exp(-2i*pi*(0:mc)'*f));
    else
        ff=exp([c0 lpccc2cc(cc,nc)]*exp(-2i*pi*(0:nc)'*f));
    end
else
    if nargin<3 || ~numel(nc) nc=np; end
    if nc==mc
        ff=exp(v_rfft([c0 cc].',2*np).');
    else
        ff=exp(v_rfft([c0 v_lpccc2cc(cc,nc)].',2*np).');
    end
    f=linspace(0,0.5,np+1);
end
if ~nargout
    subplot(2,1,2);
    plot(f,unwrap(angle(ff.')));
    xlabel('Normalized frequency f/f_s');
    ylabel('Phase (rad)');
    subplot(2,1,1);
    plot(f,db(abs(ff.')));
    xlabel('Normalized frequency f/f_s');
    ylabel('Gain (dB)');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function DecayMatrix = getDecayMatrix(EchoTimes,T2vals)
    DecayMatrix = zeros(length(EchoTimes),length(T2vals));
    for j = 1:length(T2vals)
        DecayMatrix(:,j) = exp(-EchoTimes/T2vals(j))';
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef CEC2010_F13 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- 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
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{13};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 500;
            obj.upper    = zeros(1,obj.D) + 500;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = mean(-Z.*sin(sqrt(abs(Z))),2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = -50 + mean(Z.^2,2)/100;
            PopCon(:,2) = 50*mean(sin(pi*Z/50),2);
            PopCon(:,3) = 75 - 50*(sum((Z.^2/4000),2)-prod((cos(Z./repmat(sqrt(1:size(Z,2)),size(Z,1),1))),2)+1);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ x, w ] = nco_set ( n )

%*****************************************************************************80
%
%% NCO_SET sets abscissas and weights for open Newton-Cotes quadrature.
%
%  Discussion:
%
%    The open Newton-Cotes rules use equally spaced abscissas, and
%    hence may be used with equally spaced data.
%
%    The rules are called "open" because they do not include the interval
%    endpoints.
%
%    Most of the rules involve negative weights.  These can produce loss
%    of accuracy due to the subtraction of large, nearly equal quantities.
%
%    The integral:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Abramowitz and Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%    Daniel Zwillinger, editor,
%    Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input, integer N, the order.
%    N must be between 1 and 7, and 9.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  if ( n == 1 )

    w(1) = 2.0;

  elseif ( n == 2 )

    w(1) = 1.0;
    w(2) = 1.0;

  elseif ( n == 3 )

    d = 3.0;

    w(1) =   4.0 / d;
    w(2) = - 2.0 / d;
    w(3) =   4.0 / d;

  elseif ( n == 4 )

    d = 12.0;

    w(1) = 11.0 / d;
    w(2) =  1.0 / d;
    w(3) =  1.0 / d;
    w(4) = 11.0 / d;

  elseif ( n == 5 )

    d = 10.0;

    w(1) =   11.0 / d;
    w(2) = - 14.0 / d;
    w(3) =   26.0 / d;
    w(4) = - 14.0 / d;
    w(5) =   11.0 / d;

  elseif ( n == 6 )

    d = 1440.0;

    w(1) =  1222.0 / d;
    w(2) = - 906.0 / d;
    w(3) =  1124.0 / d;
    w(4) =  1124.0 / d;
    w(5) = - 906.0 / d;
    w(6) =  1222.0 / d;

  elseif ( n == 7 )

    d = 945.0;

    w(1) =    920.0 / d;
    w(2) = - 1908.0 / d;
    w(3) =   4392.0 / d;
    w(4) = - 4918.0 / d;
    w(5) =   4392.0 / d;
    w(6) = - 1908.0 / d;
    w(7) =    920.0 / d;

  elseif ( n == 9 )

    d = 4536.0;

    w(1) =    4045.0 / d;
    w(2) = - 11690.0 / d;
    w(3) =   33340.0 / d;
    w(4) = - 55070.0 / d;
    w(5) =   67822.0 / d;
    w(6) = - 55070.0 / d;
    w(7) =   33340.0 / d;
    w(8) = - 11690.0 / d;
    w(9) =    4045.0 / d;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'NCO_SET - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of N = %d\n', n );
    fprintf ( 1, '  Legal values are 1 to 7, and 9.\n' );
    error ( 'NCO_SET - Fatal error!' );

  end
%
%  Set the abscissas.
%
  for i = 1 : n
    x(i) = ( 2 * i - n - 1 ) / ( n + 1 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = subsref(S2G,s)
%overloads subsref

switch s(1).type
  case '()'

    v = vector3d(S2G);
    [varargout{1:nargout}] = subsref(v,s);

  case '.'

    [varargout{1:nargout}] = builtin('subsref',S2G,s);

end
end

author: mtex-toolbox

repo: mtex

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


text:
function [tpl, str] = parse(tpl,target,handle,append)
%TEMPLATE/PARSE Fill in replacement fields with the class properties
%  [TPL, STR] = PARSE(TPL,TARGET,HANDLE) fills in the replacement field
%  HANDLE using previously defined variables of template TPL and store
%  it in field TARGET. HANDLE can also be a cell array of field names.
%  Output is also provided in output STR (content of TARGET).
%  [TPL, STR] = PARSE(TPL,TARGET,HANDLE,APPEND) allows to specify if
%  TARGET field is reseted before being filled or if new content is
%  appended to the previous one.

%  Copyright (C) 2003 Guillaume Flandin <Guillaume@artefact.tk>
%  $Revision: 1.0 $Date: 2003/05/05 22:19:51 $

narginchk(3,4);
if nargin == 3
	append = 0;
end

if iscellstr(handle)
	for i=1:length(handle)
		[tpl, str] = subst(tpl,handle{i});
		tpl = set(tpl,'var',target,str);
	end
elseif ischar(handle)
	[tpl, str] = subst(tpl,handle);
	if append
		tpl = set(tpl,'var',target,[get(tpl,'var',target) str]);
	else
		tpl = set(tpl,'var',target,str);
	end
else
	error('[Template] Badly formed handle.');
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function handle = plot_ranking(trackers, accuracy, robustness, varargin)
% plot_ranking Generate an A-R ranking plot
%
% The function generates an A-R ranking plot for a sequence of measurements.
%
% Input:
% - trackers (cell): A cell array of tracker structures.
% - accuracy (matrix): Accuracy ranks.
% - robustness (matrix): Robustness ranks.
% - varargin[Title] (string): A title of the plot.
% - varargin[Limit] (double): A manually set maximum rank.
% - varargin[Visible] (boolean): Is the figure visible on the display.
% - varargin[Width] (double): Figure width hint.
% - varargin[Height] (double): Figure height hint.
% - varargin[Handle] (handle): Plot on existing figure handle.
% - varargin[Legend] (boolean): Render plot legend.
%
% Output:
% - handle (handle): A figure handle.
%

    plot_title = [];
    visible = false;
    plot_limit = numel(trackers);
    width = [];
    height = [];

    handle = [];

    show_legend = true;

    for i = 1:2:length(varargin)
        switch lower(varargin{i})
            case 'title'
                plot_title = varargin{i+1};
            case 'limit'
                plot_limit = varargin{i+1};
            case 'visible'
                visible = varargin{i+1};
            case 'width'
                width = varargin{i+1};
            case 'height'
                height = varargin{i+1};
            case 'handle'
                handle = varargin{i+1};
            case 'legend'
                show_legend = varargin{i+1};
            otherwise
                error(['Unknown switch ', varargin{i},'!']) ;
        end
    end

    if isempty(handle)
        if ~visible
            handle = figure('Visible', 'off');
        else
            handle = figure();
        end
    else
        figure(handle);
    end;

    if isempty(width)
        width = iff(show_legend, 6, 4);
    end

    if isempty(height)
        height = 4;
    end

    hold on; box on; grid on;
    title(plot_title, 'interpreter','none');

    available = true(length(trackers), 1);

    for t = 1:length(trackers)

        if isnan(accuracy(t))
            available(t) = 0;
            continue;
        end;

        plot(robustness(t), accuracy(t), trackers{t}.style.symbol, 'Color', ...
            trackers{t}.style.color, 'MarkerSize', 10, 'LineWidth', trackers{t}.style.width);

    end;
    plot_labels = cellfun(@(tracker) tracker.label, trackers, 'UniformOutput', 0);
    if show_legend
        legend(plot_labels(available), 'Location', 'NorthWestOutside', 'interpreter', 'none', 'FontSize', 9);
    end;
    xlabel('Robustness rank'); set(gca, 'XDir', 'Reverse');
    ylabel('Accuracy rank'); set(gca, 'YDir', 'Reverse');
    xlim([0.9, plot_limit + 0.1]);
    ylim([0.9, plot_limit + 0.1]);

    set(handle, 'PaperUnits', 'inches', 'PaperSize', [width, height], 'PaperPosition', [0, 0, width, height]);

    hold off;

author: votchallenge

repo: toolkit-legacy


text:
function [normal,normalf] = compute_normal(vertex,face)

% compute_normal - compute the normal of a triangulation
%
%   [normal,normalf] = compute_normal(vertex,face);
%
%   normal(i,:) is the normal at vertex i.
%   normalf(j,:) is the normal at face j.
%
%   Copyright (c) 2004 Gabriel Peyr�

[vertex,face] = check_face_vertex(vertex,face);

nface = size(face,2);
nvert = size(vertex,2);
normal = zeros(3, nvert);

% unit normals to the faces
normalf = crossp( vertex(:,face(2,:))-vertex(:,face(1,:)), ...
                  vertex(:,face(3,:))-vertex(:,face(1,:)) );
d = sqrt( sum(normalf.^2,1) ); d(d<eps)=1;
normalf = normalf ./ repmat( d, 3,1 );

% unit normal to the vertex
normal = ComputeVertexNormal(face, normalf, nvert);
% normal = zeros(3,nvert);
% for i=1:nface
%     f = face(:,i);
%     for j=1:3
%         normal(:,f(j)) = normal(:,f(j)) + normalf(:,i);
%     end
% end
% normalize
d = sqrt( sum(normal.^2,1) ); d(d<eps)=1;
normal = normal ./ repmat( d, 3,1 );

% enforce that the normal are outward
v = vertex - repmat(mean(vertex,1), 3,1);
s = sum( v.*normal, 2 );
if sum(s>0)<sum(s<0)
    % flip
    normal = -normal;
    normalf = -normalf;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = crossp(x,y)
% x and y are (m,3) dimensional
z = x;
z(1,:) = x(2,:).*y(3,:) - x(3,:).*y(2,:);
z(2,:) = x(3,:).*y(1,:) - x(1,:).*y(3,:);
z(3,:) = x(1,:).*y(2,:) - x(2,:).*y(1,:);


text:
function [inputs, label] = ComputeTestExamples(curImgsLDR, curExpo, curLabel)

global param;

cropSize = param.cropSizeTraining;
border = param.border;

%%% prepare input features
[inputs, label] = PrepareInputFeatures(curImgsLDR, curExpo, curLabel, true);

% %%% crop boundaries
% inputs = CropImg(inputs, cropSize-border);
% label = CropImg(label, cropSize-border);

author: qingsenyangit

repo: AHDRNet


text:
function root_dir=mocov_util_get_root_path_containing(needle,child_dir)
% finds the nearest parent directory that contains a particular directory
%
% root_dir=mocov_util_get_root_path_containing(needle,child_dir)
%
% Inputs:
%   needle              subdirectory to look for
%   child_dir           starting directory
%
% Output:
%   rootdir             directory that contains needle and is a parent
%                       directory of child_dir. If multiple such
%                       directories exist, then the one with the longest
%                       absolute path is taken
%
% Notes:
% - this function is used by get_coverage_json to find the .git directory
%   of a repository
%
    if nargin<2
        child_dir=pwd();
    end

    root_dir=helper_get_root_dir_containing(needle,child_dir);

function root_path=helper_get_root_dir_containing(needle,child_dir)
    while true
        if has_in_dir(needle,child_dir)
            root_path=child_dir;
            return;
        end

        parent_dir=fileparts(child_dir);
        if isequal(child_dir,parent_dir)
            error('Unable to find file or directory %s',needle);
        end

        child_dir=parent_dir;
    end

function tf=has_in_dir(needle, child_dir)
    d=dir(child_dir);

    names={d.name};

    tf=(sum(strcmp(needle,names))>=1);


text:
function trials = er_assignColors(trials);
%
% trials = er_assignColors(trials);
%
% Dialog to assign colors to different conditions for
% event-related analyses.
%
% written by ras, 09/2005
if ieNotDefined('trials')
    trials = er_concatParfiles(getSelectedInplane);
end

% set up an input dialog to get colors %%%%%%%%%%%%%%%%%%%
dlg(1).fieldName = 'preset';
dlg(1).style = 'popup';
dlg(1).string = 'Use a preset color map?';
dlg(1).list = {'No, manually define below...'};
dlg(1).value = 1;
presets = mrvColorMaps;
for i = 1:length(presets), dlg(1).list{i+1} = presets{i}; end

for i = 1:length(trials.condNums)
    dlg(i+1).fieldName = sprintf('Cond%i',i);
    dlg(i+1).style = 'edit';
    dlg(i+1).string = sprintf('%s color', trials.condNames{i});
    dlg(i+1).value = num2str(trials.condColors{i});
end

resp = generalDialog(dlg, 'Set Cond Colors...');

% exit quietly if user cancels
if isempty(resp), return; end

preset = cellfind(dlg(1).list, resp.preset);
if preset>1
    % use preset value
    map = mrvColorMaps(preset-1, length(trials.condNums));
    for i = 1:size(map,1)
        trials.condColors{i} = map(i,:);
    end
else
    % parse the rest of the responses for manual definitions
    for i = 1:length(trials.condNums)
        field = sprintf('Cond%i',i);
        vals = sscanf(resp.(field),'%f %f %f');
        trials.condColors{i} = vals';
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function C=CallPrice(P, K, r, t, s)

d_1=log(P/K)+(r+s.*s/2)*t;
d_2=d_1-s*sqrt(t);

C=P.*normcdf(d_1)-K*exp(-r*t).*normcdf(d_2);


text:
function c = catalan ( n )

%*****************************************************************************80
%
%% CATALAN computes the Catalan numbers, from C(0) to C(N).
%
%  First values:
%
%     C(0)     1
%     C(1)     1
%     C(2)     2
%     C(3)     5
%     C(4)    14
%     C(5)    42
%     C(6)   132
%     C(7)   429
%     C(8)  1430
%     C(9)  4862
%    C(10) 16796
%
%  Formula:
%
%    C(N) = (2*N)! / ( (N+1) * (N!) * (N!) )
%         = 1 / (N+1) * COMB ( 2N, N )
%         = 1 / (2N+1) * COMB ( 2N+1, N+1).
%
%  Recursion:
%
%    C(N) = 2 * (2*N-1) * C(N-1) / (N+1)
%    C(N) = sum ( 1 <= I <= N-1 ) C(I) * C(N-I)
%
%  Discussion:
%
%    The Catalan number C(N) counts:
%
%    1) the number of binary trees on N vertices;
%    2) the number of ordered trees on N+1 vertices;
%    3) the number of full binary trees on 2N+1 vertices;
%    4) the number of well formed sequences of 2N parentheses;
%    5) number of ways 2N ballots can be counted, in order,
%       with N positive and N negative, so that the running sum
%       is never negative;
%    6) the number of standard tableaus in a 2 by N rectangular Ferrers diagram;
%    7) the number of monotone functions from [1..N} to [1..N} which
%       satisfy f(i) <= i for all i;
%    8) the number of ways to triangulate a polygon with N+2 vertices.
%
%  Example:
%
%    N = 3
%
%    ()()()
%    ()(())
%    (()())
%    (())()
%    ((()))
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dennis Stanton and Dennis White,
%    Constructive Combinatorics,
%    Springer Verlag, New York, 1986.
%
%  Parameters:
%
%    Input, integer N, the number of Catalan numbers desired.
%
%    Output, integer C(1:N+1), the Catalan numbers from C(0) to C(N).
%
  if ( n < 0 )
    c = [];
    return;
  end

  c(1) = 1;
%
%  The extra parentheses ensure that the integer division is
%  done AFTER the integer multiplication.
%
  for i = 1 : n
    c(i+1) = ( c(i) * 2 * ( 2 * i - 1 ) ) / ( i + 1 );
  end

  return
end


text:
%DECODER PORTION

function synth_speech = f_DECODER (aCoeff, pitch_plot, voiced, gain);

%re-calculating frame_length for this decoder,
frame_length=1;
for i=2:length(gain)
    if gain(i) == 0,
        frame_length = frame_length + 1;
    else break;
    end
end

%decoding starts here,
for b=1 : frame_length : (length(gain)),    %length(gain) should be very close
                                            %(i.e less than a frame_length error) to length(x)

    %FRAME IS VOICED OR UNVOICED
    if voiced(b) == 1,   %voiced frame
            pitch_plot_b = pitch_plot(b);
        syn_y1 = f_SYN_V (aCoeff, gain, frame_length, pitch_plot_b, b);
    else syn_y1 = f_SYN_UV (aCoeff, gain, frame_length, b); %unvoiced frame
    end

    synth_speech(b:b+frame_length-1) = syn_y1;
end


text:
function hermite_polynomial_test15 ( )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST15 tests HE_POLYNOMIAL_COEFFICIENTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2012
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST15\n' );
  fprintf ( 1, '  HE_POLYNOMIAL_COEFFICIENTS determines the probabilist''s Hermite \n' );
  fprintf ( 1, '  polynomial coefficients.\n' );

  c = he_polynomial_coefficients ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  He(%d) = \n', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( c(i+1,j+1) == 0.0 )

      elseif ( j == 0 )
        fprintf ( 1, '  %f\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %f * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '  %f * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  return
end


text:
function lp = wishpdfln(X,a,B,inverse)
%WISHPDFLN    Logarithm of Wishart probability density function.
%  See WISHPDF for argument description.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

if nargin < 3
  B = [];
end
if nargin < 4
  inverse = 0;
end

if inverse
  X = inv(X);
end
if isempty(B)
  XB = X;
  logDetB = 0;
else
  XB = X*B;
  logDetB = logdet(B);
end
d = rows(X);
d2 = (d+1)/2;
if inverse
  d2 = -d2;
end
logDetXB = (a-d2)*logdet(XB);
lp = logDetXB - trace(XB) + d2*logDetB - gammaln(a,d);

author: tminka


text:
function gB = plus(gB,xy)
% shift grains in x/y direction
%
% Syntax
%
%   % shift in x direction
%   gB = gB + [100,0]
%
% Input
%  gB - @grainBoundary
%  xy - x and y coordinates of the shift
%
% Output
%  gB - @grainBoundary

if isa(xy,'grainBoundary'), [xy,gB] = deal(gB,xy); end
if isa(xy,'vector3d'), xy = [xy.x,xy.y]; end

gB.V = gB.V + repmat(xy,size(gB.V,1),1);


text:
function mse = calc_mse(V, W, H)

    [m, n] = size(V);

    f_val = nmf_cost(V, W, H, []);

    mse = 2 * f_val / (m*n);

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function treepack_test07 ( )

%*****************************************************************************80
%
%% TREEPACK_TEST07 tests TREE_ARC_DIAM.
%
%  Discussion:
%
%    The tree is:
%
%    2---3---6---8---1---9
%       /       / \
%      7       5   4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
  nnode = 9;

  inode = [ 2, 3, 3, 6, 8, 8, 8, 1 ];
  jnode = [ 3, 7, 6, 8, 4, 5, 1, 9 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TREEPACK_TEST07\n' );
  fprintf ( 1, '  TREE_ARC_DIAM computes the diameter of a tree.\n' );

  graph_arc_print ( nnode-1, inode, jnode, '  The edge list of the tree:' );

  [ diam, label, nnode1, nnode2 ] = tree_arc_diam ( nnode, inode, jnode );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  This tree has a diameter of %d\n', diam );
  fprintf ( 1, '  between nodes %d and %d\n', nnode1, nnode2 );

  i4vec_print ( nnode, label, '  Nodes and labels:' );

  return
end


text:
function varargout = removeMeshVertices(vertices, faces, indsToRemove)
%REMOVEMESHVERTICES Remove vertices and associated faces from a mesh.
%
%   [V2, F2] = removeMeshVertices(VERTS, FACES, VERTINDS)
%   Removes the vertices specified by the vertex indices VERTINDS, and
%   remove the faces containing one of the removed vertices.
%
%   Example
%     % remove some vertices from a soccerball polyhedron
%     [v, f] = createSoccerBall;
%     plane = createPlane([.6 0 0], [1 0 0]);
%     indAbove = find(~isBelowPlane(v, plane));
%     [v2, f2] = removeMeshVertices(v, f, indAbove);
%     drawMesh(v2, f2);
%     axis equal; hold on;
%
%   See also
%     meshes3d, trimMesh

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2016-02-03, using Matlab 8.6.0.267246 (R2015b)
% Copyright 2016-2022 INRA - Cepia Software Platform

% parse inputs
if nargin == 2
    indsToRemove = faces;
    [vertices, faces] = parseMeshData(vertices);
end

% create array of indices to keep
nVertices = size(vertices, 1);
newInds = (1:nVertices)';
newInds(indsToRemove) = [];

% create new vertex array
vertices2 = vertices(newInds, :);

% compute map from old indices to new indices
oldNewMap = zeros(nVertices, 1);
for iIndex = 1:size(newInds, 1)
   oldNewMap(newInds(iIndex)) = iIndex;
end

% change labels of vertices referenced by faces
if isnumeric(faces)
    faces2 = oldNewMap(faces);
    if size(faces2,2)==1; faces2=faces2'; end
    % keep only faces with valid vertices
    faces2 = faces2(sum(faces2 == 0, 2) == 0, :);
elseif iscell(faces)
    faces2 = cell(1, length(faces));
    for iFace = 1:length(faces)
         newFace = oldNewMap(faces{iFace});
         % add the new face only if all vertices are valid
         if ~any(newFace == 0)
             faces2{iFace} = newFace;
         end
    end

    % remove empty faces
    faces2 = faces2(~cellfun(@isempty, faces2));
end

% format output arguments
varargout = formatMeshOutput(nargout, vertices2, faces2);


text:
function print_newline()
	fprintf('\n');
end

author: layumi

repo: Image-Text-Embedding