text:
function om = display_soln(om, set_type, name, idx)
%DISPLAY_SOLN  Display solution values.
%   OM.DISPLAY_SOLN()
%   OM.DISPLAY_SOLN(SET_TYPE)
%   OM.DISPLAY_SOLN(SET_TYPE, NAME)
%   OM.DISPLAY_SOLN(SET_TYPE, NAME, IDX)
%
%   Displays the model solution, including values, bounds and shadow
%   prices for variables and linear constraints, values and shadow
%   prices for nonlinear constraints, and individual cost components.
%
%   Results are displayed for each SET_TYPE or specified SET_TYPE and
%   for each named/indexed set or a specified NAME/IDX.
%
%   Inputs:
%       SET_TYPE - one of the following, specifying the type of set:
%           'var' - variables
%           'lin' - linear constraints
%           'nle' - nonlinear equality constraints
%           'nli' - nonlinear inequality constraints
%           'nlc' - nonlinear costs
%           'qdc' - quadratic costs
%         or
%           a cell array of one or more of the above
%         or
%           '' or 'all' - indicating to display all
%       NAME - (optional) char array specifying the name of the set
%       IDX  - (optional) cell array specifying the indices of the set
%
%   Examples:
%       om.display_soln('var');
%       om.display_soln({'nle', 'nli'});
%       om.display_soln('var', 'P');
%       om.display_soln('lin', 'lin_con_1');
%       om.display_soln('nle', 'nle_con_b', {2,3});
%
%   See also GET_SOLN, PARSE_SOLN.

%   MP-Opt-Model
%   Copyright (c) 2023, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% input arg handling
if nargin < 4
    idx = [];
    if nargin < 3
        name = [];
        if nargin < 2
            set_type = 'all';
        end
    end
end

fid = 1;

mu_thresh = 1e-7;
num_inf = 1e10;

%% print header
if om.is_solved()
    s = om.soln;
    if strcmp(set_type, 'all')
        set_types = fieldnames(om.set_types);   %% all set types
    elseif ~iscell(set_type)
        set_types = {set_type}; %% make set_type a cell array of char arrays
    else
        set_types = set_type;
    end

    for ss = 1:length(set_types)
        st = set_types{ss};
        om_st = om.(st);

        if om.(st).N
            if isempty(name)            %% all indices for set type
                idxs = (1:om_st.N);
            elseif isempty(idx)         %% all indices for set type & name
                idxs = [];
                for k = 1:length(om_st.order)
                    if strcmp(om_st.order(k).name, name)
                        i1 = om_st.idx.i1.(name)(om_st.order(k).idx{:});
                        iN = om_st.idx.iN.(name)(om_st.order(k).idx{:});
                        idxs = [idxs (i1:iN)];
                    end
                end
            else                        %% indices for set type, name, idx
                idxs = (om_st.idx.i1.(name)(idx{:}):om_st.idx.iN.(name)(idx{:}));
            end

            %% set type headers
            fprintf(fid, '=====  %s  =====\n', om.set_types.(st));
            %% common columns
            hdr1 = {'  idx    description                ', ...
                    '------- ----------------------------' };
            %% type specific columns
            switch st
            case {'var', 'lin'}
                hdr2 = {'   mu_lb     lb       val      ub      mu_ub', ...
                        ' -------- -------- -------- -------- --------' };
                switch st
                case 'var'
                    [v0, vl, vu] = om.params_var();
                    v = s.x;
                    mu_l = s.lambda.lower;
                    mu_u = s.lambda.upper;
                case 'lin'
                    [A, vl, vu] = om.params_lin_constraint();
                    v = A * s.x;
                    mu_l = s.lambda.mu_l;
                    mu_u = s.lambda.mu_u;
                end
            case 'nle'
                hdr2 = {'    val    lambda', ...
                        ' -------- --------' };
                lam = s.lambda.eqnonlin;
                v = om.eval_nln_constraint(s.x, 1);
            case 'nli'
                hdr2 = {'    val      ub      mu_ub', ...
                        ' -------- -------- --------' };
                mu_u = s.lambda.ineqnonlin;
                v = om.eval_nln_constraint(s.x, 0);
            case 'qdc'
                hdr2 = {'   cost  =  quad    linear  constant  average', ...
                        ' -------- -------- -------- -------- --------' };
                c = [];
                c_k = [];
                c_lin = [];
                c_avg = [];
                vv = om.get_idx();
                for k = 1:length(om_st.order)
                    n = om_st.order(k).name;
                    i = om_st.order(k).idx;
                    [QQ, cc, kk, vs] = om.params_quad_cost(n, i);
                    xx = om.varsets_x(s.x, vs, 'vector');
                    c_total = om.eval_quad_cost(s.x, n, i);
                    len = length(c_total);
                    if len == 1
                        c_constant = kk;
                        c_linear = cc' * xx;
                        c_average = c_total / sum(xx);
                    else
                        c_linear = cc .* xx;
                        c_average = c_total ./ xx;
                        if isscalar(kk)
                            c_constant = ones(len, 1)*kk/len;
                        else
                            c_constant = kk;
                        end
                    end
                    if sum(sum(QQ)) == 0 && sum(kk) == 0 && len == length(cc)
                        c_average = cc;
                    end
                    c = [c; c_total];
                    c_lin = [c_lin; c_linear];
                    c_k = [c_k; c_constant];
                    c_avg = [c_avg; c_average];
                end
                c_quad = c - c_lin - c_k;
            case 'nlc'
                c = [];
                for k = 1:length(om_st.order)
                    n = om_st.order(k).name;
                    i = om_st.order(k).idx;
                    c_total = om.eval_nln_cost(s.x, n, i);
                    c = [c; c_total];
                end
                hdr2 = {'   cost', ...
                        ' --------' };
            otherwise
                hdr2 = {'', ...
                        '' };
            end

            %% print headers
            for h = 1:length(hdr1)
                fprintf(fid, '%s\n', [hdr1{h} hdr2{h}]);
            end

            %% print data
            none = '- ';
            % none = '';
            for k = 1:length(idxs)
                ii = sprintf('%d', idxs(k));
                fmt = sprintf('%%-%ds', length(ii)+ceil((7-length(ii))/2));
                fprintf(fid, '%7s %-28s', sprintf(fmt, ii), om.describe_idx(st, idxs(k)));

                switch st
                case {'var', 'lin'}
                    if abs(mu_l(idxs(k))) < mu_thresh
                        mu_lb = sprintf(none);
                    else
                        mu_lb = sprintf_num(8, mu_l(idxs(k)));
                    end
                    if abs(mu_u(idxs(k))) < mu_thresh
                        mu_ub = sprintf(none);
                    else
                        mu_ub = sprintf_num(8, mu_u(idxs(k)));
                    end
                    if vl(idxs(k)) < -num_inf
                        lb = sprintf(none);
                    else
                        lb = sprintf_num(8, vl(idxs(k)));
                    end
                    if vu(idxs(k)) > num_inf
                        ub = sprintf(none);
                    else
                        ub = sprintf_num(8, vu(idxs(k)));
                    end
                    fprintf(fid, '%9s%9s%9s%9s%9s\n', ...
                        mu_lb, lb, sprintf_num(8, v(idxs(k))), ub, mu_ub);
                case 'nle'
                    if abs(lam(idxs(k))) < mu_thresh
                        mu_ub = sprintf(none);
                    else
                        mu_ub = sprintf_num(8, lam(idxs(k)));
                    end
                    fprintf(fid, '%9s%9s\n', sprintf_num(8, v(idxs(k))), ...
                        mu_ub);
                case 'nli'
                    if abs(mu_u(idxs(k))) < mu_thresh
                        mu_ub = sprintf(none);
                    else
                        mu_ub = sprintf_num(8, mu_u(idxs(k)));
                    end
                    fprintf(fid, '%9s%9s%9s\n', sprintf_num(8, v(idxs(k))), ...
                        '0', mu_ub);
                case 'qdc'
                    if isnan(c_avg(idxs(k)))
                        cc_avg = '- ';
                    else
                        cc_avg = sprintf_num(8, c_avg(idxs(k)));
                    end
                    fprintf(fid, '%9s%9s%9s%9s%9s\n', ...
                        sprintf_num(8, c(idxs(k))), ...
                        sprintf_num(8, c_quad(idxs(k))), ...
                        sprintf_num(8, c_lin(idxs(k))), ...
                        sprintf_num(8, c_k(idxs(k))), ...
                        cc_avg);
                case 'nlc'
                    fprintf(fid, '%9s\n', sprintf_num(8, c(idxs(k))));
                otherwise
                    fprintf(fid, '\n');
                end
            end     %% loop over idxs

            %% section footer
            switch st
            case {'var', 'lin'}
                fprintf(fid, '%s\n', [hdr1{2} hdr2{2}]);
                fprintf(fid, '%7s %-28s%9s%9s%9s%9s%9s\n', '', 'Min', ...
                    sprintf_num(8, min(mu_l)), sprintf_num(8, min(vl)), ...
                    sprintf_num(8, min(v)), ...
                    sprintf_num(8, min(vu)), sprintf_num(8, min(mu_u)));
                fprintf(fid, '%7s %-28s%9s%9s%9s%9s%9s\n', '', 'Max', ...
                    sprintf_num(8, max(mu_l)), sprintf_num(8, max(vl)), ...
                    sprintf_num(8, max(v)), ...
                    sprintf_num(8, max(vu)), sprintf_num(8, max(mu_u)));
                fprintf(fid, '\n');
            case 'nle'
                fprintf(fid, '%s\n', [hdr1{2} hdr2{2}]);
                fprintf(fid, '%7s %-28s%9s%9s\n', '', 'Min', ...
                    sprintf_num(8, min(v)), ...
                    sprintf_num(8, min(lam)));
                fprintf(fid, '%7s %-28s%9s%9s\n', '', 'Max', ...
                    sprintf_num(8, max(v)), ...
                    sprintf_num(8, max(lam)));
                fprintf(fid, '\n');
            case 'nli'
                fprintf(fid, '%s\n', [hdr1{2} hdr2{2}]);
                fprintf(fid, '%7s %-28s%9s%9s%9s\n', '', 'Max', ...
                    sprintf_num(8, max(v)), '0', ...
                    sprintf_num(8, max(mu_u)));
                fprintf(fid, '\n');
            case 'qdc'
                fprintf(fid, '%s\n', [hdr1{2} hdr2{2}]);
                fprintf(fid, '%7s %-28s%9s%9s%9s%9s%9s\n', '', ...
                    'Sum of Displayed Costs', ...
                    sprintf_num(8, sum(c(idxs))), ...
                    sprintf_num(8, sum(c_quad(idxs))), ...
                    sprintf_num(8, sum(c_lin(idxs))), ...
                    sprintf_num(8, sum(c_k(idxs))), '');
                fprintf(fid, '\n');
            case 'nlc'
                if length(c) > 1
                    fprintf(fid, '%s\n', [hdr1{2} hdr2{2}]);
                    fprintf(fid, '%7s %-28s%9s%9s%9s%9s%9s\n', '', ...
                        'Sum of Displayed Costs', ...
                        sprintf_num(8, sum(c(idxs))));
                end
                fprintf(fid, '\n');
            otherwise
                fprintf(fid, '\n');
            end
        end
    end             %% loop over set types
else
    fprintf(fid, 'Not a solved model.\n');
end

function str = sprintf_num(width, val)
val = full(val);
fmt = sprintf('%%%dg', width);
str = sprintf(fmt, val);
if length(str) > width
    fmt = sprintf('%%%d.*g', width);
    for p = width-2:-1:0
        str = sprintf(fmt, p, val);
        if length(str) <= width
            break;
        end
    end
end
assert(length(str) <= width)

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function plotkde(X, sigma2)
% Plot 2d kernel density.
% Written by Mo Chen (sth4nth@gmail.com).
if nargin < 2
    sigma2 = 1e-1;
end
level = 64;
n = 256;

X = standardize(X);

spread(X);
x_range = xlim;
y_range = ylim;

x = linspace(x_range(1),x_range(2), n);
y = linspace(y_range(2),y_range(1), n);

[a,b] = meshgrid(x,y);

z = exp(logkdepdf([a(:)';b(:)'],X,sigma2));

z = z-min(z);
z = floor(z/max(z)*(level-1));

figure;
image(reshape(z,n,n));
colormap(jet(level));
set(gca, 'XTick', [1 256]);
set(gca, 'XTickLabel', [min(x) max(x)]);
set(gca, 'YTick', [1 256]);
set(gca, 'YTickLabel', [min(y) max(y)]);
axis off

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function D = driving_function_mono_nfchoa(x0,xs,src,f,conf)
%DRIVING_FUNCTION_MONO_NFCHOA driving signal for NFC-HOA
%
%   Usage: D = driving_function_mono_nfchoa(x0,xs,src,f,conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [nx7]
%       xs          - position of virtual source or direction of plane
%                     wave / m [1x3] or [1x6]
%                     OR position and orientation of a line source [1x6]
%       src         - source type of the virtual source
%                         'pw' - plane wave (xs is the direction of the
%                                plane wave in this case)
%                         'ls' - line source
%                         'ps' - point source
%                         'fs' - focused source
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   See also: plot_sound_field, sound_field_mono_nfchoa,
%             driving_function_imp_nfchoa

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargsecondarysource(x0);
isargxs(xs);
isargpositivescalar(f);
isargchar(src);
isargstruct(conf);

%% ===== Computation ====================================================

% Secondary source positions
x0 = x0(:,1:3);

% Get maximum order of spherical harmonics
N = nfchoa_order(size(x0,1),conf);

% Source position/direction/orientation
xs = repmat(xs,[size(x0,1) 1]);

% Get driving signals
if strcmp('pw',src)
    % === Plane wave =====================================================
    % Direction of plane wave
    nk = bsxfun(@rdivide,xs,vector_norm(xs,2));
    % Driving signal
    D = driving_function_mono_nfchoa_pw(x0,nk,f,N,conf);

elseif strcmp('ps',src)
    % === Point source ===================================================
    % Driving Signal
    D = driving_function_mono_nfchoa_ps(x0,xs,f,N,conf);

elseif strcmp('ls',src)
    % === Line source ====================================================
    % Driving signal
    D = driving_function_mono_nfchoa_ls(x0,xs,f,N,conf);

elseif strcmp('fs',src)
    % === Focused source =================================================
    % Driving Signal
    D = driving_function_mono_nfchoa_fs(x0,xs,f,N,conf);

else
    error('%s: %s is not a known source type.',upper(mfilename),src);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [date, doy, dow] = gps2date(gps_week, gps_sow)

% SYNTAX:
%   [date, doy, dow] = gps2date(gps_week, gps_sow);
%
% INPUT:
%   gps_week = GPS week
%   gps_sow  = GPS seconds of week
%
% OUTPUT:
%   date = date [year month day hour min sec]
%   doy  = day of year
%   dow  = day of week
%
% DESCRIPTION:
%   Conversion from GPS time to calendar date and day of year (DOY).

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

gps_start_datenum = 723186; %This is datenum([1980,1,6,0,0,0])

gps_dow = fix(gps_sow/86400);                             %day of week
date = datevec(gps_start_datenum + 7*gps_week + gps_dow); %calendar date up to days
gps_sod = gps_sow - gps_dow*86400;                        %seconds of day
date(:,4) = floor(gps_sod/3600);                          %hours
date(:,5) = floor(gps_sod/60 - date(:,4)*60);             %minutes
date(:,6) = gps_sod - date(:,4)*3600 - date(:,5)*60;      %seconds

%day of year (DOY)
if (nargout > 1)
    doy = date2doy(datenum(date));
    doy = floor(doy);
end

%day of week
if (nargout > 2)
    dow = gps_dow;
end

author: goGPS-Project


text:
function [ y ] = vl_nneuclideanloss( x,r,vis,p )
% Euclidean Loss Layer
if nargin<4
    %forward
    %delta = x - r ;
    %y = sum(delta(:).^2) ;

    %vis=ones(1,1,21,1,'single');

    delta = (x - r).^2 ;
    temp(1,1,:,:)=vis;
    temp(1,2,:,:)=vis;
    delta = delta.*temp;
    y = sum(delta(:));

else
    %backward
    %y = 2 * p * (x - r) ;

    %vis=ones(1,1,21,1,'single');
    dx = 2 * (x - r) ;
    dx(:,1,:,:)=dx(:,1,:,:).*vis;
    dx(:,2,:,:)=dx(:,2,:,:).*vis;
    y = p*dx;

end

end

author: anilbas


text:
function Derive_EoM()
% Derive_EoM()
%
% This function generates the equations of motion for what I will call the
% "Retractable Double Pendulum" model of walking. It uses the Matlab
% symbolic toolbox to generate the equations of motion, and then
% automatically writes them to a file.
%
% The model consists of a point mass at each foot and at the hip. The feet
% are connected to this hip by an extensible, actuated, leg. Each foot has
% an 'ankle actuator' to provide a control torque, and there is another
% torque actuator at the hip, connecting the two legs.
%
% Written by Matthew Kelly
% November 29, 2013
% Cornell University
%
% See also WRITE_CONTINUOUSDYNAMICS
clc; clear; commandwindow;
addpath ../Shared
Directory = '../computerGeneratedCode';  %Write all code in this directory
disp('Running Derive_EoM...')
disp(' -> Defining Model');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                             Model                                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% This model of walking is comprised of three point masses: one at each
% foot and one at the hip. The equations are derived for three phases of
% motion: Flight, Single Stance, and Double Stance. The lets in this
% walking model are massless and have variable length. There are five
% controls available to the system: force actuator in each leg, foot torque
% when foot is in contact with the ground, and hip torque. There are six
% degrees of freedom in this model: position of foot one (x1,y1), position
% of the hip (x0,y0), and position of foot two (x2, y2).

% The position of the hip
x0 = sym('x0','real'); % Horizontal position of the hip
y0 = sym('y0','real'); % Vertical position of the hip

% The position of foot one
x1 = sym('x1','real'); % Horizontal position of foot one
y1 = sym('y1','real'); % Vertical position of foot one

% The position of foot two
x2 = sym('x2','real'); % Horizontal position of foot two
y2 = sym('y2','real'); % Vertical position of foot two

% Each leg has a small mass at the foot and a large mass at the hip. Since
% the two hip masses are coincident, they are treated as a single mass.
m1 = sym('m1','real');  % Mass of foot one
m2 = sym('m2','real');  % Mass of foot two
M = sym('M','real');  % Mass of the hip of the robot

% The system experiences a constant gravitational acceleration:
g = sym('g','real');

% Axial Force along each leg. This force is considered to be an actuator
% input to the system. Compression is positive.
F1 = sym('F1','real'); % Force in the stance leg
F2 = sym('F2','real'); % Force in the swing leg

% Constraint force normal to each leg. Connects ankle torques to system
N1 = sym('N1','real'); % Constraint force at foot one
N2 = sym('N2','real'); % constraint force between legs

% There is a torque motor connecting the two legs, and an ankle motor on
% the stance leg.
T1 = sym('T1','real'); % Ankle Torque, acting on leg one
T2 = sym('T2','real'); % Ankle Torque, acting on leg two
Thip = sym('Thip','real'); % Hip torque. From leg one acting on leg two

% Ground contact forces at each foot
H1 = sym('H1','real'); % Horizontal contact force at foot one
V1 = sym('V1','real'); % Vertical contact force at foot one
H2 = sym('H2','real'); % Horizontal contact force at foot two
V2 = sym('V2','real'); % Vertical contact force at foot wto

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        State Derivatives                                %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% The first time derivative of each state. These are considered known.
% The position of the hip
dx0 = sym('dx0','real'); % Horizontal velocity of the hip
dy0 = sym('dy0','real'); % Vertical velocity of the hip
dx1 = sym('dx1','real'); % Horizontal velocity of foot one
dy1 = sym('dy1','real'); % Vertical velocity of foot one
dx2 = sym('dx2','real'); % Horizontal velocity of foot two
dy2 = sym('dy2','real'); % Vertical velocity of foot two

% The second time derivative of each state. Goal is to find these.
ddx0 = sym('ddx0','real'); % Horizontal acceleration of the hip
ddy0 = sym('ddy0','real'); % Vertical acceleration of the hip
ddx1 = sym('ddx1','real'); % Horizontal acceleration of foot one
ddy1 = sym('ddy1','real'); % Vertical acceleration of foot one
ddx2 = sym('ddx2','real'); % Horizontal acceleration of foot two
ddy2 = sym('ddy2','real'); % Vertical acceleration of foot two

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                 Coordinate System & Kinematics                          %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Inertial reference frame:
i = [1;0;0];  %Positive horizontal axis
j = [0;1;0];  %Positive vertical axis
k = [0;0;1];  %Positive lateral axis

%Relative change in each ordinate along the legs
X1 = (x1-x0); dX1 = (dx1-dx0);
Y1 = (y1-y0); dY1 = (dy1-dy0);
X2 = (x2-x0); dX2 = (dx2-dx0);
Y2 = (y2-y0); dY2 = (dy2-dy0);

%length of both legs:
% This is a commonly used expression that is costly to calculate, so I
% numerically calculate it as an intermediate step.
L1 = sym('L1','real');
L2 = sym('L2','real');
dL1 = sym('dL1','real');
dL2 = sym('dL2','real');

%Angles of both legs, as measured in the k direction from the -j axis:
% th1 = sym('th1','real');
% th2 = sym('th2','real');
dth1 = sym('dth1','real');
dth2 = sym('dth2','real');

% Matlab symbolic toolbox was not being friendly, so I did a small amount
% of the math out by hand. I believe the following to be true:
%
%  L^2 = x^2 + y^2
% (d/dt)(L) = (x*dx + y*dy)/L
% (d/dt)(atan2(x,y)) = (dx*y - x*dy)/L^2
%

Kinematics.L1 = sqrt(X1^2 + Y1^2);
Kinematics.L2 = sqrt(X2^2 + Y2^2);

Kinematics.dL1 = (X1*dX1 + Y1*dY1)/L1;
Kinematics.dL2 = (X2*dX2 + Y2*dY2)/L2;

Kinematics.th1 = atan2(X1,Y1);
Kinematics.th2 = atan2(X2,Y2);
Kinematics.dth1 = (dX1*Y1 - X1*dY1)/L1^2;
Kinematics.dth2 = (dX2*Y2 - X2*dY2)/L2^2;

% Unit vectors pointing from the hip to foot one (a1) and it's normal (b1)
a1 = (X1*i + Y1*j)/L1;
b1 = (-Y1*i + X1*j)/L1;

% Unit vectors pointing from the hip to foot two (a2) and it's normal (b2)
a2 = (X2*i + Y2*j)/L2;
b2 = (-Y2*i + X2*j)/L2;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Position Vectors                                 %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Point 0 = Hip
% Point 1 = Foot One
% Point 2 = Foot Two

r0 = x0*i + y0*j;         % Position of Hip(Absolute)
r1 = x1*i + y1*j;         % Position of Foot One (Absolute)
r2 = x2*i + y2*j;         % Position of Foot Two (Absolute)

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                      Position Derivatives                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

dr0 = dx0*i + dy0*j;         % Velocity of Hip(Absolute)
dr1 = dx1*i + dy1*j;         % Velocity of Foot One (Absolute)
dr2 = dx2*i + dy2*j;         % Velocity of Foot Two (Absolute)

ddr0 = ddx0*i + ddy0*j;         % Acceleration of Hip(Absolute)
ddr1 = ddx1*i + ddy1*j;         % Acceleration of Foot One (Absolute)
ddr2 = ddx2*i + ddy2*j;         % Acceleration of Foot Two (Absolute)

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                Linear Momentum Balance on Foot One                      %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

sum_of_forces = simplify(V1*j + H1*i + N1*b1 - m1*g*j + F1*a1);
linear_momentum_rate = m1*ddr1;

LMB_F1 = simplify(sum_of_forces - linear_momentum_rate);
LMB_F1_i = dot(LMB_F1,i);
LMB_F1_j = dot(LMB_F1,j);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                Linear Momentum Balance on Foot Two                      %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

sum_of_forces = simplify(V2*j + H2*i + N2*b2 - m2*g*j + F2*a2);
linear_momentum_rate = m2*ddr2;

LMB_F2 = simplify(sum_of_forces - linear_momentum_rate);
LMB_F2_i = dot(LMB_F2,i);
LMB_F2_j = dot(LMB_F2,j);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                Linear Momentum Balance on Hip                           %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

sum_of_forces = simplify(-F1*a1 - N1*b1 - F2*a2 - N2*b2 - M*g*j);
linear_momentum_rate = M*ddr0;

LMB_H = simplify(sum_of_forces - linear_momentum_rate);
LMB_H_i = dot(LMB_H,i);
LMB_H_j = dot(LMB_H,j);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%               Angular Momentum Balance on Leg One                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

sum_of_torques = (T1 + Thip - N1*L1)*k;
angular_momentum_rate = 0*k;   %Leg One has no mass

AMB_L1_k = dot(simplify(sum_of_torques - angular_momentum_rate),k);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%               Angular Momentum Balance on Leg One                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

sum_of_torques = (T2 - Thip - N2*L2)*k;
angular_momentum_rate = 0*k;   %Leg Two has no mass

AMB_L2_k = dot(simplify(sum_of_torques - angular_momentum_rate),k);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%      Collect Implicit Equations of Motion and Phase constraints         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% All phases of motion share these equations. An additional four equations
% are required to solve the system, which come from the 'definition' of
% each phase of motion.

Physics = [         LMB_F1_i;
                    LMB_F1_j;
                    LMB_F2_i;
                    LMB_F2_j;
                    LMB_H_i;
                    LMB_H_j;
                    AMB_L1_k;
                    AMB_L2_k    ];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Solve Flight Dynamics                            %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% During flight, the contact forces are all equal to zero. Treat as knowns.
% %             H1==0;
% %             V1==0;
% %             H2==0;
% %             V2==0;

Unknowns = [    ddx0;    ddy0;
                ddx1;    ddy1;
                ddx2;    ddy2;
                N1;     N2];

disp(' -> Solving Flight Dynamics')
Soln = jacobSolve(Physics,Unknowns);

Dyn.Flight = Soln;
Dyn.Flight.H1 = sym('0');
Dyn.Flight.V1 = sym('0');
Dyn.Flight.H2 = sym('0');
Dyn.Flight.V2 = sym('0');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                 Solve Single Stance One Dynamics                        %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% During single stance one, contact forces at Foot Two are equal to zero.
% Additionally, we know that the velocity and acceleration of Foot One are
% both also equal to zero.
% %             ddx1==0;
% %             ddy1==0;
% %             H2==0;
% %             V2==0;

Unknowns = [    ddx0;   ddy0;
                H1;     V1
                ddx2;   ddy2;
                N1;     N2];

disp(' -> Solving Single Stance One Dynamics')
Soln = jacobSolve(Physics,Unknowns);

Dyn.SingleOne = Soln;
Dyn.SingleOne.ddx1 = sym('0');
Dyn.SingleOne.ddy1 = sym('0');
Dyn.SingleOne.H2 = sym('0');
Dyn.SingleOne.V2 = sym('0');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                 Solve Single Stance Two Dynamics                        %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% During single stance two, contact forces at Foot One are equal to zero.
% Additionally, we know that the velocity and acceleration of Foot Two are
% both also equal to zero.
%
%             ddx2 == 0;
%             ddy2 == 0;
%             H1==0;
%             V1==0;

Unknowns = [    ddx0;   ddy0;
                ddx1;   ddy1;
                H2;     V2;
                N1;     N2];

disp(' -> Solving Single Stance Two Dynamics')
Soln = jacobSolve(Physics,Unknowns);

Dyn.SingleTwo = Soln;
Dyn.SingleTwo.ddx2 = sym('0');
Dyn.SingleTwo.ddy2 = sym('0');
Dyn.SingleTwo.H1 = sym('0');
Dyn.SingleTwo.V1 = sym('0');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                   Solve Double Stance Dynamics                          %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% During double stance, both feet remain stationary.
% Additionally, we know that the velocity and acceleration of Foot One are
% both also equal to zero.
%               ddx1 == 0;
%               ddy1 == 0;
%               ddx2 == 0;
%               ddy2 == 0;

% In Double Stance, Both feet remain stationary. Add constraint equation.
% Note that the constraint is on acceleration, NOT velocity. This means
% that these equations will only yield the desired solutions if the initial
% velocity of Foot Two == 0;

Unknowns = [    ddx0;   ddy0
                N1;     N2;
                H1;     V1;
                H2;     V2];

disp(' -> Solving Double Stance Dynamics')
Soln = jacobSolve(Physics,Unknowns);

Dyn.Double = Soln;
Dyn.Double.ddx1 = sym('0');
Dyn.Double.ddy1 = sym('0');
Dyn.Double.ddx2 = sym('0');
Dyn.Double.ddy2 = sym('0');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%              Write Continuous Dynamics Function Files                   %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
disp(' -> Writing Continuous Dynamics Functions')

States = cell(12,2);
States(1,:) = {'x0','(m) Hip horizontal position'};
States(2,:) = {'y0','(m) Hip vertical position'};
States(3,:) = {'x1','(m) Foot One horizontal position'};
States(4,:) = {'y1','(m) Foot One vertical position'};
States(5,:) = {'x2','(m) Foot Two horizontal position'};
States(6,:) = {'y2','(m) Foot Two vertical position'};

States(7,:) = {'dx0','(m/s) Hip horizontal velocity'};
States(8,:) = {'dy0','(m/s) Hip vertical velocity'};
States(9,:) = {'dx1','(m/s) Foot One horizontal velocity'};
States(10,:) = {'dy1','(m/s) Foot One vertical velocity'};
States(11,:) = {'dx2','(m/s) Foot Two horizontal velocity'};
States(12,:) = {'dy2','(m/s) Foot Two vertical velocity'};

Contacts = cell(4,2);
Contacts(1,:) = {'H1','(N) Foot One, horizontal contact force'};
Contacts(2,:) = {'V1','(N) Foot One, vertical contact force'};
Contacts(3,:) = {'H2','(N) Foot Two, horizontal contact force'};
Contacts(4,:) = {'V2','(N) Foot Two, vertical contact force'};

Actuators = cell(5,2);
Actuators(1,:) = {'F1','(N) Compresive axial force in Leg One'};
Actuators(2,:) = {'F2','(N) Compresive axial force in Leg Two'};
Actuators(3,:) = {'T1','(Nm) External torque applied to Leg One'};
Actuators(4,:) = {'T2','(Nm) External torque applied to Leg Two'};
Actuators(5,:) = {'Thip','(Nm) Torque acting on Leg Two from Leg One'};

Parameters = cell(4,2);
Parameters(1,:) = {'m1','(kg) Foot One mass'};
Parameters(2,:) = {'m2','(kg) Foot Two mass'};
Parameters(3,:) = {'M','(kg) Hip mass'};
Parameters(4,:) = {'g','(m/s^2) Gravity'};

CommonExpressions = cell(2,2);
CommonExpressions(1,:) = {'L1s','L1.^2'};
CommonExpressions(2,:) = {'L2s','L2.^2'};

FileData = cell(3,4);
FileData(1,:) = {'Flight','dynamics_flight',...
    'Dymanics Model: retractable double pendulum biped',...
    'Motion Phase: Flight'};
FileData(2,:) = {'SingleOne','dynamics_singleStanceOne',...
    'Dymanics Model: retractable double pendulum biped',...
    'Motion Phase: Single Stance One'};
FileData(3,:) = {'SingleTwo','dynamics_singleStanceTwo',...
    'Dymanics Model: retractable double pendulum biped',...
    'Motion Phase: Single Stance Two'};
FileData(4,:) = {'Double','dynamics_doubleStance',...
    'Dymanics Model: retractable double pendulum biped',...
    'Motion Phase: Double Stance'};

FileWritingSetup = Make_Struct(Dyn,States,Contacts,Actuators,...
    Kinematics,CommonExpressions,Parameters,FileData,Directory);

Write_ContinuousDynamics(FileWritingSetup);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                          Find System Energy                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

Energy.Potential.m1 = simplify(m1*g*dot(r1,j));
Energy.Potential.m2 = simplify(m2*g*dot(r2,j));
Energy.Potential.M = simplify(M*g*dot(r0,j));

Energy.Kinetic.m1 = simplify(0.5*m1*norm(dr1).^2);
Energy.Kinetic.m2 = simplify(0.5*m2*norm(dr2).^2);
Energy.Kinetic.M = simplify(0.5*M*norm(dr0)^2);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Write Energy Function                           %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
disp(' -> Writing Kinematics Function')

FileWritingSetup = Make_Struct(Energy,States,Parameters,Directory);

Write_Energy(FileWritingSetup);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                          Write Power Function                           %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
disp(' -> Writing Actuator Power Function')

Power.legOne = F1*dL1;
Power.legTwo = F2*dL2;
Power.ankleOne = T1*dth1;
Power.ankleTwo = T2*dth2;
Power.hip = Thip*(dth2-dth1);

FileWritingSetup = Make_Struct(Power,Kinematics,...
    States,Actuators,Directory);

Write_ActuatorPower(FileWritingSetup);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                          Write Kinematics Function                           %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
disp(' -> Writing Kinematics Function')

FileWritingSetup = Make_Struct(Kinematics, States, Directory);

Write_Kinematics(FileWritingSetup);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                             Impact Equations                            %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Assume that any transitions between phases can be modeled by external
% impluses being applied to either Foot One or Foot Two.
%
% Each of the point masses in the model is connected to the others through
% force and torque actuators. Assume that at any instant in time there is a
% finite force or torque being provided by each of those actuators. In the
% limit as the duration of the impact goes to zero, the impulse transfered
% across those actuators also goes to zero.
%
% Thus, I conclude that the collision maps can be applied by solving the
% kinematic constraints before and after the collision.
%
% Since the masses in this system are essientally decoupled from eachother
% as far as impacts are concerned, the impact map can be simplified to the
% following statement: "The state is unchanged, except for the velocity of
% the mass that strikes the ground, which should go to zero".

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Write Impact Map Function Files                      %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
disp(' -> Writing Phase Map Function')

FileWritingSetup = ...
    Make_Struct(States, Directory);

Write_PhaseMap(FileWritingSetup);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Write Conversion Function                            %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
disp(' -> Writing Conversion Function')

FileWritingSetup = Make_Struct(States,Actuators,Contacts,Directory);

Write_Convert(FileWritingSetup);

disp('DONE!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    SUB-FUNCTIONS                                  %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Soln = jacobSolve(Equations,Unknowns)
%
% FUNCTION:
%   Solve a nonlinear system of equations by assuming that it is linear in
%   the unknown variables (which is true for these mechanics problems)
%
% ARGUMENTS:
%   Equations = [Nx1] vector of symbolic expressions that are equal to zero
%   Unknowns = [Nx1] vector of symbolic variables to solve Equations for
%
% RETURNS:
%   Soln = a struct with a field for each unknown
%
% The matlab solve command seems to have a problem with solving large
% systems of non-linear equations. In the case of classical mechanics
% problems, it turns out that these systems are not too hard to solve
% because they are actually linear in the accelerations and constraint
% forces. Assuming that this is true, then you can transform the equations
% into a linear system by taking partial derivatives. Once this step is
% done, then matlab does a great job of solving the linear system.
%
% MATH:
%   Equations = 0;                  % By Definition
%   Equations = A*x + b;            % Assume: form, A independant* of x
%   A = jacobian(Equations wrt x);  %
%   b = Equations - A*x;            %
%   0 = A*x + b;                    %
%   x = -A\b;                       % Solved!
%

A = simplify(jacobian(Equations,Unknowns));
b = simplify(Equations - A*Unknowns);
x = simplify(-A\b);

for i=1:length(Unknowns)
   Soln.(char(Unknowns(i))) = x(i);
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function sesync_data = convert2sesync(measurements)

sesync_data.R = {measurements.R};
sesync_data.t = {measurements.t};
sesync_data.kappa = {measurements.kappa};
sesync_data.tau = {measurements.tau};
edge_start = cell2mat( {measurements.i} );
edge_end = cell2mat( {measurements.j} );
sesync_data.edges = [edge_start', edge_end'];

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function volr = volresize_t2(vol, k, varargin)
    % volresize(vol, k, {interp_type}
    interp_type = 1;
    if nargin >= 3
        interp_type = varargin{1};
    end
    if numel(k) == 1
        k = k * [1,1,1];
    end

    newsz = round(k.*size(vol));

    tmp = zeros(max(newsz, size(vol)));
    tmp(1:size(vol,1), 1:size(vol,2), 1:size(vol,3)) = vol;
    [n1, n2, n3] = ndgrid(1:size(tmp,1), 1:size(tmp,2), 1:size(tmp,3));

    k=1./k;
    d = 1 - k;
    T = cat(4, k(1) * n1 - n1 + d(1), k(2) * n2 - n2 + d(2), k(3) * n3 - n3  + d(3));

    volr = imdeform3(tmp, T, interp_type);

    volr = volr(1:newsz(1), 1:newsz(2), 1:newsz(3));
end

repo: pTVreg


text:
% File Name: checkValidity.m
% Author: Heung-Il Suk
% Cite: H.-I. Suk and S.-W. Lee, "A Novel Bayesian Framework for Discriminative
% Feature Extraction in Brain-Computer Interfaces," IEEE Trans. on PAMI,
% 2012. (Accepted)

function sample = opt_checkValidity( sample )
% We enforce the samples to be within min_freq and max_freq

min_freq = 0.05;
max_freq = 40;

if sample(1) <= min_freq
    sample(1) = min_freq;
end

if sample(2) <= min_freq
    sample(2) = min_freq;
end

if sample(1) > max_freq
    sample(1) = max_freq;
end

if sample(2) > max_freq
    sample(2) = max_freq;
end

if sample(1) > sample(2)
    temp = sample(2);
    sample(2) = sample(1);
    sample(1) = temp;
end

if (sample(1) - sample(2)) < 1
    sample(1) = sample(1);
    sample(2) = sample(2) + 0.5;
end

author: PatternRecognition

repo: OpenBMI


text:
function [Rz] = R_z(psi)

% Input sanity check
if nargin < 1
    psi = sym('psi_t','real');
end
assert(numel(psi) == 1,'Angle must be a scalar [1x1].');
assert(isa(psi,'sym') || isnumeric(psi),'Angle must be a symbolic or numeric value.');

% Rotate about the z-axis (yaw)
Rz = [cos(psi) -sin(psi) 0;
      sin(psi)  cos(psi) 0;
             0         0 1];
end

repo: OpenMAS


text:
function [src,sfi] = zSpecularFreeImage(src)
Lambda = 0.6;
camDark = 20;
r = src.rgb(:,:,1); g = src.rgb(:,:,2); b = src.rgb(:,:,3);
src.i(intersect(intersect(find(r<camDark),find(g<camDark)),...
    find(b<camDark))) = z.CAMERA_DARK;
c = z.MaxChroma(src.rgb);
dI = (z.Max(src.rgb).*(3*c-1))./(c*(3*Lambda-1));
sI = (z.Total(src.rgb)-dI)/3;
drgb = min(255,max(0,src.rgb-sI));
sfi.rgb = drgb;
end

author: vitorsr

repo: SIHR

save_path: github-repos/MATLAB/vitorsr-SIHR


text:
function value = i4_is_odd ( i )

%*****************************************************************************80
%
%% I4_IS_ODD returns TRUE if I is odd.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the integer to be tested.
%
%    Output, logical VALUE, is TRUE if I is odd.
%
  value = ( mod ( i+1, 2 ) == 0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = witps(varargin)
% VL_WITPS  Inverse thin-plate spline warping
%   [X1,X2]=VL_WITPS(XP1,XP2,Y,Yp) computes the inverse thin-plate spline
%   (TPS) warp of the points XP1,XP2.
%
%   Remark::
%     The inverse of a thin-plate spline in general is NOT a
%     thin-plate spline and some splines do not have an inverse.  This
%     function uses Gauss-Newton to compute a set of points (X1,X2)
%     such that [XP1,XP2]=VL_WTPS(X1,X2,Y,Yp).
%
%   See also: VL_WTPS(), VL_HELP().
[varargout{1:nargout}] = vl_witps(varargin{:});

author: yihui-he

repo: panorama

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


text:
function out = cwpt2_interface(in, direction, filter_type, trans_type, btree)
% function out = cwpt2_interface(in, direction, filter_type, trans_type, btree)
%
% Forward 2D Continuous Wavelet Transform.
% Inverse 2D Continuous Wavelet Transform.
% Shifted Inverse 2D Continuous Wavelet Transform.
%
% Forward 2D Continuous Wavelet Packet Transform.
% Inverse 2D Continuous Wavelet Packet Transform.
%
% This is a wrapper for cwpt2, cwpt2i, c3wpt, c3wpti.
% It allows to select arbitrary filters to be used for the transform.
% It allows classical (quad) or 3-way wavelet elementary transforms.
%
% direction =   'forward' | 'inverse'
% filter_type = '7-9' | 'spline'
% trans_type =  'quad' | 'tri'
% btree =       a) for wavelet packet transforms, the corresponding btree
%               b) for forward wavelet transform, the number of scales
%               c) for inverse wavelet transform, must bet set to 0
%               d) for shifted inverse wavelet transform, must be set to 1
%
% A forward transform given a M x N image provides a M x N x S 3-D matrix.
%
% ================
% For Wavelet Packet Transform, type 'help cwpt2' for more information about the
% order of subbands.
%
% ================
% For Wavelet Transform, S = ( 1 + 3 * nb_scales ) and the order of subbands is:
%
% HL, LH, HH (scale 1),
% HL, LH, HH (scale 2),
% ...
% LL, HL, LH, HH (scale nb_scales),
%
% where the first filtering direction is the first memory raster direction of data.
% For instance, LH means for Matlab images low-filtering along columns (first dimension)
% and high-filtering along rows (second dimension).
%
% ================
% For 3-way Wavelet Transform, S = ( 1 + 2 * nb_scales ) and the order of subbands is:
%
% H0, 0H (scale 1),
% H0, 0H (scale 2),
% ...
% LL, H0, 0H (scale nb_scales),
%
% ================
% A shifted inverse (classical or 3-way) wavelet transform considers a transform of
% scale N as a transform of scale (N+1)
%

if ~isequal(class(in),'double')
    warning('Input value should be of type ''double''.');
    in = double(in);
end;

switch lower(filter_type)
case '7-9'
    % analysis
    alf = [ 0.03782845550699547033,-0.02384946501937996663,-0.1106244044184234859, 0.3774028556126538536, 0.8526986790094033264, 0.3774028556126538536,-0.1106244044184233194,-0.02384946501937999092, 0.03782845550699546339];
    ahf = [ 0.06453888262893847649,-0.04068941760955852721,-0.4180922732222120963, 0.7884856164056645023,-0.4180922732222124294,-0.04068941760955835374, 0.06453888262893842098];
    % reconstruction
    rlf = [-0.03226944131446923825,-0.02034470880477926361, 0.2090461366111060482, 0.3942428082028322511, 0.2090461366111062147,-0.02034470880477917687,-0.03226944131446921049];
    % ahf / 2; * -1 ^ odd
    rhf = [ 0.01891422775349773516, 0.01192473250968998331,-0.05531220220921174296,-0.1887014278063269268, 0.4263493395047016632,-0.1887014278063269268,-0.0553122022092116597, 0.01192473250968999546, 0.01891422775349773169];
    % alf / 2; * -1 ^ even
case 'spline'
    h = MakeONFilter_copy('Battle',3);
    g = ((-1 * ones(1,41)).^(0:40)).*h;

    alf = h;
    ahf = g;
    rlf = h./2;
    rhf = g./2;
otherwise
    error('Unknown filter type');
end;

switch lower(direction)
case 'forward'

    if length(btree)==1                         % wavelet; btree = nb_scales
        btree = cwpt2_btree(btree, 2);          % 2 = wavelet
    else                                        % packets
        if isequal(trans_type,'tri')
            warning('3-way transform has not been tested for WP transform');
        end;
    end;

    switch lower(trans_type)
    case 'quad'
        out = cwpt2(in, btree, alf, floor(length(alf)/2), ahf, floor(length(ahf)/2));
    case 'tri'
        out = c3wpt(in, btree, alf, floor(length(alf)/2), ahf, floor(length(ahf)/2));
    otherwise
        error('Unknown transform type');
    end;

case 'inverse'

    if length(size(in))~=3
        error('The input argument must be a M x N x S 3-D matrix.');
    end;

    if length(btree)==1                         % wavelet; btree = scale_shift
        if btree~=0 && btree~=1
            error('btree parameter must be 0 or 1 for inverse wavelet transform')
        end;

        nb_scales = size(in,3) - 1;
        if isequal(trans_type,'quad')
            if mod(nb_scales,3)~=0
                error('Invalid data: the size along the dimension 3 must be ( 1 + 3 * nb_scales ).');
            end;
            nb_scales = btree + (nb_scales / 3);
        else
            if mod(nb_scales,2)~=0
                error('Invalid data: the size along the dimension 3 must be ( 1 + 2 * nb_scales ).');
            end;
            nb_scales = btree + (nb_scales / 2);
        end;

        btree = cwpt2_btree(nb_scales, 2);      % 2 = wavelet
    else                                        % packets
        if isequal(trans_type,'tri')
            warning('3-way transform has not been tested for WP transform');
        end;
    end;

    switch lower(trans_type)
    case 'quad'
        out = cwpt2i(in, btree, rlf, floor(length(rlf)/2), rhf, floor(length(rhf)/2));
    case 'tri'
        out = c3wpti(in, btree, alf, floor(length(alf)/2), ahf, floor(length(ahf)/2), rlf, floor(length(rlf)/2), rhf, floor(length(rhf)/2));
    otherwise
        error('Unknown transform type');
    end;

otherwise
    error('Incorrect transform direction');
end;

% if the output argument is a string, this is an error message
if isequal(class(out),'char')
    error(out);
end;

% 2D Continuous Wavelet Packet Transform package
% (c) 2002-2005 Let It Wave, all rights reserved

author: gpeyre

repo: matlab-toolboxes


text:
function [x0] = unpackStereoParams(camK, camKc, prjK, prjKc, R, T)
% Unpack parameters for camera and projector calibration using lsqnonlin
% This function unpacks vectorized camera and projector parameters from
% matrices.
% See also: lsqnonlin, Calibration.packStereoParams

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%% camera

% intrinsics
fx0_c = camK(1,1);
fy0_c = camK(2,2);
cx0_c = camK(1,3);
cy0_c = camK(2,3);

%% projector
% intrinsics
fx0_p = prjK(1,1);
fy0_p = prjK(2,2);
cx0_p = prjK(1,3);
cy0_p = prjK(2,3);
% sk0_p = prjK(1,2);

if(nargin == 2)
    x0 = [fx0_c, fy0_c, cx0_c, cy0_c, fx0_p, fy0_p, cx0_p, cy0_p]; % ignore skew
else
    if(nargin > 4)
        % rotation matrix to rotation vector
        % rvec = rotationMatrixToVector(R);

        rvec = cv.Rodrigues(R)';
        x0 = [fx0_c, fy0_c, cx0_c, cy0_c, camKc, fx0_p, fy0_p, cx0_p, cy0_p, prjKc, rvec, T']; % ignore skew
    else
        x0 = [fx0_c, fy0_c, cx0_c, cy0_c, camKc, fx0_p, fy0_p, cx0_p, cy0_p, prjKc]; % ignore skew
    end
end

end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function varargout = linop_test( op, cmode, maxits )

%LINOP_TEST Performs an adjoint test on a linear operator.
%    LINOP_TEST( OP ) attempts to verify that a linear operator OP obeys
%    the inner product test: <A*x,y> = <x,A'*y> for all x, y. OP must be a
%    TFOCS linear operator with hard-coded size information; that is,
%    OP([],0) must return valid size info.
%
%    When called with a single argument, LINOP_TEST creates real test
%    vectors for X and Y. To test complex operators, use the two-argument
%    version LINOP_TEST( OP, cmode ), where:
%        cmode = 'R2R': real input, real output
%        cmode = 'R2C': real input, complex output
%        cmode = 'R2CC': real input, conjugate-symmetric complex output
%        cmode = 'C2R': complex input, real output
%        cmode = 'CC2R': conjugate-symmetric complex input, real output
%        cmode = 'C2C': complex input, complex output
%
%    The conjugate-symmetric options follow the symmetry conventions of Matlab's FFT.
%
%    LINOP_TEST( OP, CMODE, MAXITS ) performs MAXITS iterations of the
%    test loop. MAXITS=25 is the default.
%
%   myNorm = LINOP_TEST(...) returns an estimate of the myNorm of
%       the linear operator.

error(nargchk(1,3,nargin));
if isnumeric( op ),
    if nargin < 2 || isempty( cmode )
        disp('warning: for matrix inputs, this assumes the cmode is ''C2C'' unless you specify otherwise');
        cmode = 'C2C';
    end
    op = linop_matrix( op, cmode );
end
y_conjSymmetric     = false;
x_conjSymmetric     = false;
if nargin < 2 || isempty( cmode ),
    x_real = true;
    y_real = true;
else
    switch upper( cmode ),
        case 'R2R', x_real = true; y_real = true;
        case 'R2C', x_real = true; y_real = false;
        case 'R2CC', x_real = true; y_real = false;  y_conjSymmetric    = true;
        case 'C2R', x_real = false; y_real = true;
        case 'CC2R', x_real = false; y_real = true;  x_conjSymmetric    = true;
        case 'C2C', x_real = false; y_real = false;
        case 'CC2CC', x_real = false; y_real = false;
            x_conjSymmetric    = true;
            y_conjSymmetric    = true; % this is probably never going to happen though
        otherwise, error( 'Invalid cmode: %s', cmode );
    end
end
if nargin < 3 || isempty(maxits),
    maxits = 25;
end
sz = op([],0);
if ~iscell(sz)
    sz = { [sz(2),1], [sz(1),1] };
elseif isempty(sz{1})
    disp('warning: could not detect the size; this often happens when using a scalar input, which represents scaling');
    if isempty( sz{2} )
        disp('  Proceeding under the assumption that the domain is 1D');
        sz{1} = [1,1];
        sz{2} = [1,1];
    else
        disp('  Proceeding under the assumption that the domain equals the range');
        sz{1}   = sz{2};
    end
elseif isempty(sz{2})
    disp('warning: could not detect the size; this often happens when using a scalar input, which represents scaling');
    disp('  Proceeding under the assumption that the domain equals the range');
    sz{2}   = sz{1};
end
nf = 0;
na = 0;
errs = zeros(1,maxits+1);
nxe = 0; nye = 0;
for k = 1 : maxits,

    %
    % The adjoint test
    %

    if x_real,
        x = randn(sz{1});
    else
        x = randn(sz{1})+1j*randn(sz{1});
        if x_conjSymmetric
            x = make_conj_symmetrix( x );
        end
    end

    if y_real,
        y = randn(sz{2});
    else
        y = randn(sz{2})+1j*randn(sz{2});
        if y_conjSymmetric
            y = make_conj_symmetrix( y );
        end
    end

    nx = myNorm(x);
    Ax = op(x,1);
    nf = max( nf, myNorm(Ax)/nx );
    Ax_y = tfocs_dot( Ax, y );

    ny = myNorm(y);
    Ay = op(y,2);
    na = max( na, myNorm(Ay) / ny );
    Ay_x = tfocs_dot( x, Ay );

    errs(k) = abs(Ax_y-Ay_x)/(nx*ny);

    %
    % The myNorm iteration
    %

    if nxe == 0,
        if x_real,
            xx = randn(sz{1});
        else
            xx = randn(sz{1}) + 1j*randn(sz{1});
            if x_conjSymmetric
                xx = make_conj_symmetrix( xx );
            end
        end
        nxe = myNorm(xx);
    end
    yy = op(xx/nxe,1);
    nye = max(realmin,myNorm(yy));
    xx = op(yy/nye,2);
    nxe = myNorm(xx);

end

%
% Use the estimated singular vectors for a final adjoint est
%

if nxe > 0,
    Ax_y = tfocs_dot( op(xx,1), yy );
    Ay_x = tfocs_dot( op(yy,2), xx );
    errs(end) = abs(Ax_y-Ay_x) / (nxe*nye);
end

%
% Display the output
%

nmax = max(nye,nxe);
myNorm_err = abs(nye-nxe) / nmax;
peak_err = max(errs) / nmax;
mean_err = mean(errs) / nmax;
rc = { 'complex', 'real', 'complex symmetric' };
fprintf( 'TFOCS linear operator test:\n' );
fprintf( '   Input size:  [' ); fprintf( ' %d', sz{1} ); fprintf( ' ], %s\n', rc{x_real+1+2*x_conjSymmetric} );
fprintf( '   Output size: [' ); fprintf( ' %d', sz{2} ); fprintf( ' ], %s\n', rc{y_real+1+2*y_conjSymmetric} );
fprintf( 'After %d iterations:\n', maxits  );
fprintf( '    myNorm estimates (forward/adjoint/error): %g/%g/%g\n', nye, nxe, myNorm_err );
fprintf( '       Gains: forward %g, adjoint %g\n', nf, na );
fprintf( '    Inner product error:\n' );
fprintf( '       Mean (absolute/relative): %g/%g\n', mean(errs), mean_err );
fprintf( '       Peak (absolute/relative): %g/%g\n', max(errs), peak_err );
fprintf( '       (inner product errors should 1e-10 or smaller)\n');

good = true;
if myNorm_err/max( nye, nxe ) > 1e-4
    fprintf('  Detected mismatch in forward/adjoint norm estimates. This is potentially a bad sign. Check your implementation\n');
    good = false;
end
if mean_err > 1e-8
    fprintf('  The mean error (relative) is high. This is a bad sign. Check your implementation\n');
    good = false;
end
if good
    fprintf('  Allowing for some roundoff error, there are no obvious errors. This is good.\n');
end

if nargout > 0
    varargout{1} = mean([nye,nxe]);
end

% Improvement as suggested by Graham Coleman
% Allows for 3D arrays.
% This also changes default behavior of 2D arrays
% to now use the Frobenius norm instead of spectral norm.
% This is wise, since it's a much quicker computation.
function y = myNorm(x)
y = norm( x(:) );

function y = make_conj_symmetrix( y )
ny = size( y, 1 );
y(1,:)  = real(y(1,:));         % DC component is 0
if round(ny/2) == ny/2  % even
    y(ny/2+1,:)     = real(y(ny/2+1,:));    % Nyquist component is 0
    y( ny:-1:(ny/2+2) )     = conj( y(2:ny/2) );
else                    % odd
    y( ny:-1:((ny+1)/2+1) ) = conj( y(2:((ny+1)/2)) );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function mustBeInRange(prop, lower, upper)

mustBeGreaterThanOrEqual(prop, lower)
mustBeLessThanOrEqual(prop, upper)

end

author: european-central-bank

repo: BEAR-toolbox


text:
function retval=statglmeval(action,fn,varargin)
%STATGLMEVAL  Evaluate or test link function in proper environment
%   STATGLMEVAL('eval',FN,ARGS,...) evaluates the function FN in an
%   environment in which certain functions such as LOGIT and
%   D_LOGIT are defined.  This allows the function FN to be either
%   a user-defined function or one of the pre-defined functions
%   provided in the GLMFIT function, without contaminating the name
%   space with those pre-defined functions.
%
%   STATGLMEVAL('testlink',FN) test for the existence of FN as a
%   function handle, inline function, or text string containing the
%   name of a function in an M-file.

%   Author:  Tom Lane, 3-6-2000
%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 1.3 $  $Date: 2002/02/04 19:25:46 $

switch(action)
 case 'eval'
   retval = feval(fn,varargin{:});

 case 'testlink'
   c = class(fn);
   fnclass = class(@logit);
   if (isequal(c,fnclass) | isequal(c,'inline') | ...
       (isequal(c,'char') & ~isempty(which(fn))))
      retval = 1;
   else
      retval = 0;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following functions are related to the link function that
% links the distribution parameter mu with the linear combination
% eta of predictor variables.  For example, the logit link defines
%
%        eta = log(mu/(1-mu))

%%%% functions for identity link
function a=identity(b)
a=b;

function a=d_identity(b)
a = ones(size(b));

function b=i_identity(a)
b=a;

%%%% functions for logit link
function a=logit(p)
a = log(p ./ (1-p));

function a=d_logit(p)
a = 1 ./ max(eps, (p .* (1-p)));

function p=i_logit(a)
p = 1 ./ (1 + exp(-a));

%%%% functions for probit link
function a=probit(p)
a = norminv(p);

function a=d_probit(p)
a = 1 ./ max(eps, normpdf(norminv(p)));

function p=i_probit(a)
p = normcdf(a);

%%%% functions for complementary log log link
function a=comploglog(p)
a = log(-log(1-max(eps,p)));

function a=d_comploglog(p)
a = 1 ./ -(max(eps,1-p) .* log(1-max(eps,p)));

function p=i_comploglog(a)
p = 1 - exp(-exp(a));

%%%% functions for log log link
function a=logloglink(p)
a = log(-log(max(eps,p)));

function a=d_logloglink(p)
a = 1 ./ (max(eps, p) .* log(max(eps,p)));

function p=i_logloglink(a)
p = exp(-exp(a));

%%%% functions for log link
function a=d_log(b);
a = 1 ./ max(eps,b);

function b=i_log(a);
b = exp(a);

%%%% functions for reciprocal link
function a=reciprocal(b)
a = 1 ./ max(eps, b);

function a=d_reciprocal(b);
a = -1 ./ max(eps,b).^2;

function b=i_reciprocal(a);
b = 1 ./ max(eps, a);

%%%% functions for power link
function a=power(b,p)
if (p==0)
   a = log(max(eps,b));
else
   a = max(eps,b) .^ p;
end

function a=d_power(b,p);
if (p==0)
   a = 1 ./ max(eps,b);
else
   a = p * max(eps,b).^(p-1);
end

function b=i_power(a,p);
if (p==0)
   b = exp(a);
else
   b = max(eps,a) .^ (1/p);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following functions define the variance

function a=normalvariance(b)
a = ones(size(b));

function a=poissonvariance(b)
a = b;

function a=binomialvariance(p,N)
a = p .* (1-p) ./ N;

function a=gammavariance(b)
a = b.^2;

function a=inversegaussianvariance(b)
a = b.^3;

author: zouchuhang

repo: LayoutNet


text:
function zc = zerocross(x,wintype,winamp,winlen)
%ENERGY   Short-time energy computation.
%   y = ZEROCROSS(X,WINTYPE,WINAMP,WINLEN) computes the short-time enery of
%   the sequence X.
%
%   WINTYPE defines the window type. RECTWIN, HAMMING, HANNING, and
%   BLACKAMN are the possible choices. WINAMP sets the amplitude of the
%   window and the length of the window is WINLEN.
%
%   See also RECTWIN, HAMMING, HANNING, BARTLETT, BLACKMAN.
%
%   Author: Nabin Sharma
%   Date: 2009/03/15

error(nargchk(1,4,nargin,'struct'));

% generate x[n] and x[n-1]
x1 = x;
x2 = [0, x(1:end-1)];

% generate the first difference
firstDiff = sgn(x1)-sgn(x2);

% magnitude only
absFirstDiff = abs(firstDiff);

% lowpass filtering with window
zc = winconv(absFirstDiff,wintype,winamp,winlen);

author: Sable

repo: mcbench-benchmarks


text:
function cleanimg=deislands2d(img,sizelim)
%
% cleanimg=deislands2d(img,sizelim)
%
% remove isolated islands on a 2D image below speicified size limit
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
% 	img: a 2D binary image
% 	sizelim: a integer as the maximum pixel size of a isolated region
%
% output:
% 	cleanimg: a binary image after removing islands below sizelim
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

img=squeeze(img);
maxisland=-1;
if(nargin==2) maxisland=sizelim; end

islands={};

cleanimg=zeros(size(img));
if(sum(img(:)))
    img=imclose(img, strel('disk',3));
    islands=bwislands(img);
end

if(length(islands))
    % remove small islands of the foreground
    maxblock=-1;
    maxblockid=-1;
    if(maxisland<0)
      for i=1:length(islands)
        if(length(islands{i})>maxblock)
            maxblockid=i;
            maxblock=length(islands{i});
        end
      end
      if(maxblock>0)
          cleanimg(islands{maxblockid})=1;
      end
    else
      for i=1:length(islands)
        if(length(islands{i})>maxisland)
            cleanimg(islands{i})=1;
        end
      end
    end

    % remote small islands of the background

    cleanimg=imfill(cleanimg,'holes');
end

author: fangq


text:
function [S] = spm_mci_random (mcmc,R,v,M,U,Y)
% Random effects estimation
% FORMAT [S] = spm_mci_random (mcmc,R,v,M,U,Y)
%
% mcmc  Sampling parameters
% R     Priors on random effects (R.pE, R.pC)
% v     Fixed effects
% M     Model Structure (single subject)
% U     Inputs (single subject)
% Y     Data (single subject)
%
% S     Samples, [maxits x M.n]
%
% Uses Langevin Monte Carlo
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_random.m 6697 2016-01-27 14:57:28Z spm $

% Data precision

try, verbose=mcmc.verbose; catch, verbose=0; end
try, maxits=mcmc.maxits; catch, maxits=64; end
try, plot_int=mcmc.plot_int; catch, plot_int=1; end
try, h=mcmc.h; catch, h=0.5; end

% Assign init/flow/out params as fixed/random effects
try, assign=mcmc.assign; catch, assign=[]; end

% Compute eigen-parameterisation
M = spm_mci_minit (M);

try, init=mcmc.init; catch, init=R.pE; end
% Initial param in eigenspace
xinit = init;

% Read data points and time indices
try, ind=Y.ind; catch, ind=1:M.N; end
Nt=length(ind);
y=Y.y;

% Sample matrix
Nx=size(R.pE,1);
x = zeros(maxits,Nx);
x(1,:) = xinit';

ipC=pinv(R.pC);
logdet_RCp=spm_logdet(R.pC);

if verbose, figure; end

% Tune h by monitoring acceptance rate
tune_h=1;
acc_block=64;
acc_low=0.3;
acc_high=0.7;
total_acc_target=64; % Number of accepted samples to get
acc=zeros(maxits,1);

i=1;
while (i <= maxits),

    if verbose
        if mod(i,plot_int) == 0 && i > 2
            spm_mci_progress (x,E,i);
        end
    end

    if mod(i,acc_block)==0 && tune_h
        % Change step size h ?
        Nacc=sum(acc(i-acc_block+1:i-1));
        prop_acc=Nacc/acc_block;
        if prop_acc < acc_low
            if verbose, disp('Decreasing step size ...'); end
            h=h/2;
        elseif prop_acc > acc_high
            if verbose, disp('Increasing step size ...'); end
            h=h*2;
        end
    end

    % Proposal (first proposal always accepted)
    if i==1
        pos=x(1,:)';
        curr=[];
    else
        pos=spm_normrnd(curr.mu,curr.Cp,1);
    end

    % Quantities re proposal
    prop.pos=pos;
    if isfield(M,'IS')
        % Other model types
        [j,iCpY,st,L] = spm_mci_joint_grad (pos,M,U,Y);
    else
        % Differential equation models
        [dLdp,iCpY,st] = spm_mci_grad_curve (assign,pos,v,M,U,Y,'random');
        % Gradient of log prior
        ep = pos-R.pE;
        dlogprior=-ep'*ipC;
        j=dLdp+dlogprior;

        [p_init,p_flow] = spm_mci_init_flow (assign,pos,v,M);
        log_like = spm_mci_like_ind (p_flow,p_init,M,U,Y);
        log_prior = -0.5*ep'*ipC*ep + logdet_RCp-0.5*M.Np*log(2*pi);
        L=log_like+log_prior;
    end

    if st==-1
        error('Integration problem in spm_mci_random.m');
    end
    prop.L = L;

    % Posterior covariance under local linear approximation
    prop.Cp = h*inv(iCpY+ipC);
    prop.mu = pos+0.5*h*prop.Cp*j(:);

    prop.iCp = inv(prop.Cp);
    prop.logdetCp = spm_logdet(prop.Cp);

    % Accept proposal ?
    [curr,accepted,bayes_fb(i),dL(i)] = spm_mci_mh_update(curr,prop,verbose);
    acc(i)=accepted;
    E(i) = -curr.L;
    if i > 1
        dEdit(i-1)=100*(E(i)-E(i-1))/E(i-1);
    end
    x(i,:)=curr.pos;

    i=i+1;
end

S=x(1:i-1,:)';

author: spm


text:
function printTrainBatch(data, params)
  if params.isBi
    printSent(2, data.srcInput(1, :), params.srcVocab, '  src input 1: ');
    fprintf(2, '  src mask: %s\n', num2str(data.srcMask(1, :)));
  end

  printSent(2, data.tgtInput(1, :), params.tgtVocab, '  tgt input 1: ');
  printSent(2, data.tgtOutput(1, :), params.tgtVocab, '  tgt output 1: ');

  fprintf(2, '  tgt mask: %s\n', num2str(data.tgtMask(1, :)));
end

author: lmthang


text:
function gpcf = gpcf_linearLogistic(varargin)
%GPCF_LINEARLOGISTIC  Create a covariance function corresponding to
%                     logistic mean function
%
%  Description
%    GPCF = GPCF_LINEARLOGISTIC('PARAM1',VALUE1,'PARAM2,VALUE2,...) creates
%    a covariance function structure corresponding to logistic mean
%    function in which the named parameters have the specified values. Any
%    unspecified parameters are set to default values.
%
%    GPCF = GPCF_LINEARLOGISTIC(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    The logistic functional form is given by
%         h(x) = w* (logitinv(a.*x + b) - 0.5);
%    By giving a zero mean Gaussian prior for weight,
%         w ~ N(0, coeffSigma2)
%    the prior for h(x) is
%         h(x) ~ N(0, H(x)*H(x)'*coeffSigma2) )
%    where H(x) = [h(x(1)), ... , h(x(n))]' and, hence,
%    H(x)*H(x)'*coeffSigma2) is the covariance function related to the
%    logistic mean function.
%
%    Parameters for linearLogistic (dot product) covariance function
%      a                 - regression coefficient of linear part [1].
%                          For identifiability a is restricted to positive
%                          values by log transformation
%      b                 - intercept of the linear part [0].
%                          b can be positive or negative.
%      a_prior           - prior for a [prior_gaussian('s2',10)]
%      b_prior           - prior for b [prior_gaussian('s2',10)]
%      coeffSigma2       - prior variance for regressor coefficients [10]
%                          This can be either scalar corresponding
%                          to a common prior variance or vector
%                          defining own prior variance for each
%                          coefficient.
%      coeffSigma2_prior - prior structure for coeffSigma2 [prior_logunif]
%      selectedVariables - vector defining which inputs are used [all]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%  Example:
%   a=0.2;
%   b = -10;
%   x = linspace(0,100,100)';
%   y = 3.*(logitinv(a.*x + b) - 0.5) + 0.1*randn(100,1);
%   cf = gpcf_linearLogistic('a', a, 'b', b, 'selectedVariables', 1)  ;
%   gp = gp_set('cf', cf);
%   Ef = gp_pred(gp,x,y,x);
%   figure,plot(x,y,'.'), hold on, plot(x,Ef,'k')
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*, MEAN_*
%
% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2008-2010 Jaakko Riihimäki
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GPCF_LINEARLOGISTIC';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('coeffSigma2',10, @(x) isvector(x) && all(x>0));
  ip.addParamValue('a',1, @(x) isvector(x) && all(x>0));
  ip.addParamValue('b',0, @(x) isvector(x) );
  ip.addParamValue('coeffSigma2_prior',prior_logunif, @(x) isstruct(x) || isempty(x));
  ip.addParamValue('a_prior',prior_gaussian('s2',10), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('b_prior',prior_gaussian('s2',10), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('selectedVariables',[], @(x) isvector(x) && all(x>0));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    init=true;
    gpcf.type = 'gpcf_linearLogistic';
  else
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_linearLogistic')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  % Initialize parameter
  if init || ~ismember('coeffSigma2',ip.UsingDefaults)
    gpcf.coeffSigma2=ip.Results.coeffSigma2;
  end
  if init || ~ismember('a',ip.UsingDefaults)
    gpcf.a=ip.Results.a;
  end
  if init || ~ismember('b',ip.UsingDefaults)
    gpcf.b=ip.Results.b;
  end

  % Initialize prior structure
  if init
    gpcf.p=[];
  end
  if init || ~ismember('coeffSigma2_prior',ip.UsingDefaults)
    gpcf.p.coeffSigma2=ip.Results.coeffSigma2_prior;
  end
  if init || ~ismember('a_prior',ip.UsingDefaults)
      gpcf.p.a=ip.Results.a_prior;
  end
  if init || ~ismember('b_prior',ip.UsingDefaults)
      gpcf.p.b=ip.Results.b_prior;
  end

  if ~ismember('selectedVariables',ip.UsingDefaults)
    selectedVariables=ip.Results.selectedVariables;
    if ~isempty(selectedVariables)
      gpcf.selectedVariables = selectedVariables;
    end
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_linearLogistic_pak;
    gpcf.fh.unpak = @gpcf_linearLogistic_unpak;
    gpcf.fh.lp = @gpcf_linearLogistic_lp;
    gpcf.fh.lpg = @gpcf_linearLogistic_lpg;
    gpcf.fh.cfg = @gpcf_linearLogistic_cfg;
    gpcf.fh.cfdg = @gpcf_linearLogistic_cfdg;
    gpcf.fh.cfdg2 = @gpcf_linearLogistic_cfdg2;
    gpcf.fh.ginput = @gpcf_linearLogistic_ginput;
    gpcf.fh.ginput2 = @gpcf_linearLogistic_ginput2;
    gpcf.fh.ginput3 = @gpcf_linearLogistic_ginput3;
    gpcf.fh.ginput4 = @gpcf_linearLogistic_ginput4;
    gpcf.fh.cov = @gpcf_linearLogistic_cov;
    gpcf.fh.trcov  = @gpcf_linearLogistic_trcov;
    gpcf.fh.trvar  = @gpcf_linearLogistic_trvar;
    gpcf.fh.recappend = @gpcf_linearLogistic_recappend;
  end

end

function [w, s, h] = gpcf_linearLogistic_pak(gpcf, w)
%GPCF_GPCF_LINEARLOGISTIC_PAK  Combine GP covariance function parameters into one vector
%
%  Description
%    W = GPCF_GPCF_LINEARLOGISTIC_PAK(GPCF) takes a covariance function
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for
%    example in energy and gradient computations.
%
%       w = [ log(gpcf.coeffSigma2)
%             (hyperparameters of gpcf.coeffSigma2)
%              log(gpcf.a)
%             (hyperparameters of gpcf.a)
%              gpcf.b
%             (hyperparameters of gpcf.b)]'
%
%  See also
%    GPCF_GPCF_LINEARLOGISTIC_UNPAK

  w = []; s = {}; h =[];
  if ~isempty(gpcf.p.coeffSigma2)
    w = log(gpcf.coeffSigma2);
    if numel(gpcf.coeffSigma2)>1
      s = [s; sprintf('log(linearLogistic.coeffSigma2 x %d)',numel(gpcf.coeffSigma2))];
    else
      s = [s; 'log(linearLogistic.coeffSigma2)'];
    end
    h = [h ones(1, numel(gpcf.coeffSigma2))];
    % Hyperparameters of coeffSigma2
    [wh, sh, hh] = gpcf.p.coeffSigma2.fh.pak(gpcf.p.coeffSigma2);
    sh=strcat(repmat('prior-', size(sh,1),1),sh);
    w = [w wh];
    s = [s; sh];
    h = [h 1+hh];
  end

  if ~isempty(gpcf.p.a)
      w = [w log(gpcf.a)];
      if numel(gpcf.a)>1
          s = [s; sprintf('log(linearLogistic.a x %d)',numel(gpcf.a))];
      else
          s = [s; 'log(linearLogistic.a)'];
      end
      h = [h ones(1, numel(gpcf.a))];
      % Hyperparameters of a
      [wh, sh, hh] = gpcf.p.a.fh.pak(gpcf.p.a);
      sh=strcat(repmat('prior-', size(sh,1),1),sh);
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
  end

  if ~isempty(gpcf.p.b)
      w = [w gpcf.b];
      if numel(gpcf.b)>1
          s = [s; sprintf('linearLogistic.b x %d',numel(gpcf.b))];
      else
          s = [s; 'linearLogistic.b'];
      end
      h = [h ones(1, numel(gpcf.b))];
      % Hyperparameters of b
      [wh, sh, hh] = gpcf.p.b.fh.pak(gpcf.p.b);
      sh=strcat(repmat('prior-', size(sh,1),1),sh);
      w = [w wh];
      s = [s; sh];
      h = [h 1+hh];
  end

end

function [gpcf, w] = gpcf_linearLogistic_unpak(gpcf, w)
%GPCF_GPCF_LINEARLOGISTIC_UNPAK  Sets the covariance function parameters
%                   into the structure
%
%  Description
%    [GPCF, W] = GPCF_GPCF_LINEARLOGISTIC_UNPAK(GPCF, W) takes a covariance
%    function structure GPCF and a hyper-parameter vector W, and
%    returns a covariance function structure identical to the
%    input, except that the covariance hyper-parameters have been
%    set to the values in W. Deletes the values set to GPCF from
%    W and returns the modified W. This is a mandatory subfunction
%    used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = [ log(gpcf.coeffSigma2)
%             (hyperparameters of gpcf.coeffSigma2)
%              log(gpcf.a)
%             (hyperparameters of gpcf.a)
%              log(gpcf.b)
%             (hyperparameters of gpcf.b)]'
%
%  See also
%   GPCF_GPCF_LINEARLOGISTIC_PAK

  gpp=gpcf.p;

  if ~isempty(gpp.coeffSigma2)
    i2=length(gpcf.coeffSigma2);
    i1=1;
    gpcf.coeffSigma2 = exp(w(i1:i2));
    w = w(i2+1:end);

    % Hyperparameters of coeffSigma2
    [p, w] = gpcf.p.coeffSigma2.fh.unpak(gpcf.p.coeffSigma2, w);
    gpcf.p.coeffSigma2 = p;
  end

  if ~isempty(gpp.a)
      i2=length(gpcf.a);
      i1=1;
      gpcf.a = exp(w(i1:i2));
      w = w(i2+1:end);

      % Hyperparameters of a
      [p, w] = gpcf.p.a.fh.unpak(gpcf.p.a, w);
      gpcf.p.a = p;
  end
  if ~isempty(gpp.b)
      i2=length(gpcf.b);
      i1=1;
      gpcf.b = w(i1:i2);
      w = w(i2+1:end);

      % Hyperparameters of b
      [p, w] = gpcf.p.b.fh.unpak(gpcf.p.b, w);
      gpcf.p.b = p;
  end
end

function lp = gpcf_linearLogistic_lp(gpcf)
%GPCF_GPCF_LINEARLOGISTIC_LP  Evaluate the log prior of covariance function
%                             parameters
%
%  Description
%    LP = GPCF_GPCF_LINEARLOGISTIC_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%   GPCF_GPCF_LINEARLOGISTIC_PAK, GPCF_GPCF_LINEARLOGISTIC_UNPAK,
%   GPCF_GPCF_LINEARLOGISTIC_LPG, GP_E

% Evaluate the prior contribution to the error. The parameters that
% are sampled are from space W = log(w) where w is all the "real" samples.
% On the other hand errors are evaluated in the W-space so we need take
% into account also the  Jacobian of transformation W -> w = exp(W).
% See Gelman et al. (2013), Bayesian Data Analysis, third edition, p. 21.
  lp = 0;
  gpp=gpcf.p;

  if ~isempty(gpp.coeffSigma2)
    lp = gpp.coeffSigma2.fh.lp(gpcf.coeffSigma2, gpp.coeffSigma2) + sum(log(gpcf.coeffSigma2));
  end
  if ~isempty(gpp.a)
    lp = lp + gpp.a.fh.lp(gpcf.a, gpp.a) + sum(log(gpcf.a));
  end
  if ~isempty(gpp.b)
    lp = lp + gpp.b.fh.lp(gpcf.b, gpp.b);
  end
end

function lpg = gpcf_linearLogistic_lpg(gpcf)
%GPCF_GPCF_LINEARLOGISTIC_LPG  Evaluate gradient of the log prior with respect
%                 to the parameters.
%
%  Description
%    LPG = GPCF_GPCF_LINEARLOGISTIC_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_GPCF_LINEARLOGISTIC_PAK, GPCF_GPCF_LINEARLOGISTIC_UNPAK,
%    GPCF_GPCF_LINEARLOGISTIC_LP, GP_G

  lpg = [];
  gpp=gpcf.p;

  if ~isempty(gpcf.p.coeffSigma2)
      lll=length(gpcf.coeffSigma2);
      lpgs = gpp.coeffSigma2.fh.lpg(gpcf.coeffSigma2, gpp.coeffSigma2);
      lpg = [lpg lpgs(1:lll).*gpcf.coeffSigma2+1 lpgs(lll+1:end)];
  end
  if ~isempty(gpcf.p.a)
      lll=length(gpcf.a);
      lpgs = gpp.a.fh.lpg(gpcf.a, gpp.a);
      lpg = [lpg lpgs(1:lll).*gpcf.a+1 lpgs(lll+1:end)];
  end
  if ~isempty(gpcf.p.b)
      lll=length(gpcf.b);
      lpgs = gpp.b.fh.lpg(gpcf.b, gpp.b);
      lpg = [lpg lpgs(1:lll) lpgs(lll+1:end)];
  end

end

function DKff = gpcf_linearLogistic_cfg(gpcf, x, x2, mask, i1)
%GPCF_GPCF_LINEARLOGISTIC_CFG  Evaluate gradient of covariance function
%                 with respect to the parameters
%
%  Description
%    DKff = GPCF_GPCF_LINEARLOGISTIC_CFG(GPCF, X) takes a covariance
%    function structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_GPCF_LINEARLOGISTIC_CFG(GPCF, X, X2) takes a covariance
%    function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_GPCF_LINEARLOGISTIC_CFG(GPCF, X, [], MASK) takes a
%    covariance function structure GPCF, a matrix X of input vectors and
%    returns DKff, the diagonal of gradients of covariance matrix
%    Kff = k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_GPCF_LINEARLOGISTIC_CFG(GPCF,X,X2,MASK,i) takes a
%    covariance function structure GPCF, a matrix X of input vectors and
%    returns DKff, the gradient of covariance matrix Kff =
%    k(X,X2), or k(X,X) if X2 is empty, with respect to ith
%    hyperparameter. This subfunction is needed when using
%    memory save option in gp_set.
%
%  See also
%   GPCF_GPCF_LINEARLOGISTIC_PAK, GPCF_GPCF_LINEARLOGISTIC_UNPAK,
%   GPCF_GPCF_LINEARLOGISTIC_LP, GP_G

  [n, m] =size(x);

  DKff = {};

  if nargin==5
    % Use memory save option
    savememory=1;
    if i1==0
      % Return number of hyperparameters
      DKff=0;
      if ~isempty(gpcf.p.coeffSigma2)
        DKff=length(gpcf.coeffSigma2);
      end
      if ~isempty(gpcf.p.a)
          DKff=DKff+length(gpcf.a);
      end
      if ~isempty(gpcf.p.b)
          DKff=DKff+length(gpcf.b);
      end
      return
    end
  else
    savememory=0;
  end

  % Evaluate: DKff{1} = d Kff / d coeffSigma2
  % NOTE! Here we have already taken into account that the parameters are
  % transformed through log() and thus dK/dlog(p) = p * dK/dp

  h = logitinv(gpcf.a.*x + gpcf.b) - 0.5;

  % evaluate the gradient for training covariance
  if nargin == 2 || (isempty(x2) && isempty(mask))

    if isfield(gpcf, 'selectedVariables')
      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*h(:,gpcf.selectedVariables)*(h(:,gpcf.selectedVariables)');
        else
          if ~savememory
            i1=1:length(gpcf.coeffSigma2);
          end
          for ii1=i1
            DD = gpcf.coeffSigma2(ii1)*h(:,gpcf.selectedVariables(ii1))*(h(:,gpcf.selectedVariables(ii1))');
            DD(abs(DD)<=eps) = 0;
            DKff{ii1}= (DD+DD')./2;
          end
        end
      end
      if ~isempty(gpcf.p.a)
          ii1= length(DKff) +1 ;
          hh = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b).^2.*...
              gpcf.a.*x(:,gpcf.selectedVariables).*exp(-gpcf.a.*x(:,gpcf.selectedVariables)-gpcf.b);
          DKff{ii1} = gpcf.coeffSigma2.* (hh*h(:,gpcf.selectedVariables)' + h(:,gpcf.selectedVariables)*hh');
      end
      if ~isempty(gpcf.p.b)
          ii1= length(DKff) +1 ;
          hh = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b).^2.*exp(-gpcf.a.*x(:,gpcf.selectedVariables)-gpcf.b);
          DKff{ii1} = gpcf.coeffSigma2* (hh*h(:,gpcf.selectedVariables)' + h(:,gpcf.selectedVariables)*hh');
      end
    else
      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*h*(h');
        else
          if isa(gpcf.coeffSigma2,'single')
            epsi=eps('single');
          else
            epsi=eps;
          end
          if ~savememory
            i1=1:length(gpcf.coeffSigma2);
          end
          DKff=cell(1,length(i1));
          for ii1=i1
            DD = gpcf.coeffSigma2(ii1)*h(:,ii1)*(h(:,ii1)');
            DD(abs(DD)<=epsi) = 0;
            DKff{ii1}= (DD+DD')./2;
          end
        end
      end
      if ~isempty(gpcf.p.a)
          ii1= length(DKff) +1 ;
          hh = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b).^2.*gpcf.a.*x.*exp(gpcf.a.*x+gpcf.b);
          DKff{ii1} = gpcf.coeffSigma2* (hh*h' + h*hh');
      end
      if ~isempty(gpcf.p.b)
          ii1= length(DKff) +1 ;
          hh = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b).^2.*exp(gpcf.a.*x+gpcf.b);
          DKff{ii1} = gpcf.coeffSigma2* (hh*h' + h*hh');
      end
    end

    % Evaluate the gradient of non-symmetric covariance (e.g. K_fu)
  elseif nargin == 3 || isempty(mask)
    if size(x,2) ~= size(x2,2)
      error('gpcf_linearLogistic -> _ghyper: The number of columns in x and x2 has to be the same. ')
    end
    error('gpcf_linearLogistic -> _ghyper: "nargin == 3 || isempty(mask)" not implemented')
    if isfield(gpcf, 'selectedVariables')
      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*x(:,gpcf.selectedVariables)*(x2(:,gpcf.selectedVariables)');
        else
          if ~savememory
            i1=1:length(gpcf.coeffSigma2);
          end
          for ii1=i1
            DKff{ii1}=gpcf.coeffSigma2(ii1)*x(:,gpcf.selectedVariables(ii1))*(x2(:,gpcf.selectedVariables(ii1))');
          end
        end
      end
    else
      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*x*(x2');
        else
          if ~savememory
            i1=1:m;
          end
          for ii1=i1
            DKff{ii1}=gpcf.coeffSigma2(ii1)*x(:,ii1)*(x2(:,ii1)');
          end
        end
      end
    end
    % Evaluate: DKff{1}    = d mask(Kff,I) / d coeffSigma2
    %           DKff{2...} = d mask(Kff,I) / d coeffSigma2
  elseif nargin == 4 || nargin == 5
    error('gpcf_linearLogistic -> _ghyper: "nargin == 4 || nargin == 5" not implemented')
    if isfield(gpcf, 'selectedVariables')
      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*sum(x(:,gpcf.selectedVariables).^2,2); % d mask(Kff,I) / d coeffSigma2
        else
          if ~savememory
            i1=1:length(gpcf.coeffSigma2);
          end
          for ii1=i1
            DKff{ii1}=gpcf.coeffSigma2(ii1)*(x(:,gpcf.selectedVariables(ii1)).^2); % d mask(Kff,I) / d coeffSigma2
          end
        end
      end
    else
      if ~isempty(gpcf.p.coeffSigma2)
        if length(gpcf.coeffSigma2) == 1
          DKff{1}=gpcf.coeffSigma2*sum(x.^2,2); % d mask(Kff,I) / d coeffSigma2
        else
          if ~savememory
            i1=1:m;
          end
          for ii1=i1
            DKff{ii1}=gpcf.coeffSigma2(ii1)*(x(:,ii1).^2); % d mask(Kff,I) / d coeffSigma2
          end
        end
      end
    end
  end
  if savememory
    DKff=DKff{i1};
  end
end

function C = gpcf_linearLogistic_cov(gpcf, x1, x2, varargin)
%GP_GPCF_LINEARLOGISTIC_COV  Evaluate covariance matrix between two input
%                            vectors
%
%  Description
%    C = GP_GPCF_LINEARLOGISTIC_COV(GP, TX, X) takes in covariance function of
%    a Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%  See also
%    GPCF_GPCF_LINEARLOGISTIC_TRCOV, GPCF_GPCF_LINEARLOGISTIC_TRVAR,
%    GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  if isfield(gpcf, 'selectedVariables')
      h1 = logitinv(gpcf.a.*x1(:,gpcf.selectedVariables) + gpcf.b) - 0.5;
      h2 = logitinv(gpcf.a.*x2(:,gpcf.selectedVariables) + gpcf.b) - 0.5;
      C = h1*diag(gpcf.coeffSigma2)*(h2');
  else
      h1 = logitinv(gpcf.a.*x1 + gpcf.b) - 0.5;
      h2 = logitinv(gpcf.a.*x2 + gpcf.b) - 0.5;
      C = h1*diag(gpcf.coeffSigma2)*(h2');
  end
  C(abs(C)<=eps) = 0;
end

function C = gpcf_linearLogistic_trcov(gpcf, x)
%GP_GPCF_LINEARLOGISTIC_TRCOV  Evaluate training covariance matrix of
%                              inputs
%
%  Description
%    C = GP_GPCF_LINEARLOGISTIC_TRCOV(GP, TX) takes in covariance function
%    of a Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This
%    is a mandatory subfunction used for example in prediction and
%    energy computations.
%
%  See also
%    GPCF_GPCF_LINEARLOGISTIC_COV, GPCF_GPCF_LINEARLOGISTIC_TRVAR, GP_COV,
%    GP_TRCOV

  if isfield(gpcf, 'selectedVariables')
      h = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b) - 0.5;
      C = h*diag(gpcf.coeffSigma2)*(h');
  else
      h = logitinv(gpcf.a.*x + gpcf.b) - 0.5;
      C = h*diag(gpcf.coeffSigma2)*(h');
  end
  C(abs(C)<=eps) = 0;
  C = (C+C')./2;

end

function C = gpcf_linearLogistic_trvar(gpcf, x)
%GP_GPCF_LINEARLOGISTIC_TRVAR  Evaluate training variance vector
%
%  Description
%    C = GP_GPCF_LINEARLOGISTIC_TRVAR(GPCF, TX) takes in covariance
%    function of a Gaussian process GPCF and matrix TX that contains
%    training inputs. Returns variance vector C. Every element i
%    of C contains variance of input i in TX. This is a mandatory
%    subfunction used for example in prediction and energy computations.
%
%
%  See also
%    GPCF_GPCF_LINEARLOGISTIC_COV, GP_COV, GP_TRCOV

  if length(gpcf.coeffSigma2) == 1
    if isfield(gpcf, 'selectedVariables')
        h = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b) - 0.5;
      C=gpcf.coeffSigma2.*sum(h.^2,2);
    else
        h = logitinv(gpcf.a.*x + gpcf.b) - 0.5;
      C=gpcf.coeffSigma2.*sum(h.^2,2);
    end
  else
    if isfield(gpcf, 'selectedVariables')
        h = logitinv(gpcf.a.*x(:,gpcf.selectedVariables) + gpcf.b) - 0.5;
      C=sum(repmat(gpcf.coeffSigma2, size(x,1), 1).*h.^2,2);
    else
        h = logitinv(gpcf.a.*x + gpcf.b) - 0.5;
      C=sum(repmat(gpcf.coeffSigma2, size(h,1), 1).*h.^2,2);
    end
  end
  C(abs(C)<eps)=0;

end

function reccf = gpcf_linearLogistic_recappend(reccf, ri, gpcf)
%RECAPPEND Record append
%
%  Description
%    RECCF = GPCF_GPCF_LINEARLOGISTIC_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the parameters. Returns RECCF which contains all
%    the old samples and the current samples from GPCF. This
%    subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_linearLogistic';

    % Initialize parameters
    reccf.coeffSigma2= [];
    reccf.a= [];
    reccf.b= [];

    % Set the function handles
    reccf.fh.pak = @gpcf_linearLogistic_pak;
    reccf.fh.unpak = @gpcf_linearLogistic_unpak;
    reccf.fh.lp = @gpcf_linearLogistic_lp;
    reccf.fh.lpg = @gpcf_linearLogistic_lpg;
    reccf.fh.cfg = @gpcf_linearLogistic_cfg;
    reccf.fh.cfdg = @gpcf_linearLogistic_cfdg;
    reccf.fh.cfdg2 = @gpcf_linearLogistic_cfdg2;
    reccf.fh.ginput = @gpcf_linearLogistic_ginput;
    reccf.fh.ginput2 = @gpcf_linearLogistic_ginput2;
    reccf.fh.ginput3 = @gpcf_linearLogistic_ginput3;
    reccf.fh.ginput4 = @gpcf_linearLogistic_ginput4;
    reccf.fh.cov = @gpcf_linearLogistic_cov;
    reccf.fh.trcov  = @gpcf_linearLogistic_trcov;
    reccf.fh.trvar  = @gpcf_linearLogistic_trvar;
    reccf.fh.recappend = @gpcf_linearLogistic_recappend;
    reccf.p=[];
    reccf.p.coeffSigma2=[];
    if ~isempty(ri.p.coeffSigma2)
      reccf.p.coeffSigma2 = ri.p.coeffSigma2;
    end
    if ~isempty(ri.p.a)
      reccf.p.a = ri.p.a;
    end
    if ~isempty(ri.p.b)
      reccf.p.b = ri.p.b;
    end

  else
    % Append to the record
    gpp = gpcf.p;

    % record coeffSigma2
    reccf.coeffSigma2(ri,:)=gpcf.coeffSigma2;
    if isfield(gpp,'coeffSigma2') && ~isempty(gpp.coeffSigma2)
      reccf.p.coeffSigma2 = gpp.coeffSigma2.fh.recappend(reccf.p.coeffSigma2, ri, gpcf.p.coeffSigma2);
    end

    reccf.a(ri,:)=gpcf.a;
    if isfield(gpp,'a') && ~isempty(gpp.a)
      reccf.p.a = gpp.a.fh.recappend(reccf.p.a, ri, gpcf.p.a);
    end

    reccf.b(ri,:)=gpcf.b;
    if isfield(gpp,'b') && ~isempty(gpp.b)
      reccf.p.b = gpp.b.fh.recappend(reccf.p.b, ri, gpcf.p.b);
    end

    if isfield(gpcf, 'selectedVariables')
      reccf.selectedVariables = gpcf.selectedVariables;
    end
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function x = InvMeyerPartition(y, L, deg)
% InvMeyerPartition: Inverse the scale partitioning
%  Usage:
%     y = MeyerPartition(xhat, L, deg);
%  Inputs:
%    y      Vector of length 2*n; signal separated into disjoint scales
%    L      Coarsest Scale
%    deg    Degree of the polynomial
%  Outputs:
%    x     Vector of length n
%  Description
%    Performs the inverse of MeyerPartition; i.e., reconstruct
%    an object from its different scale contributions.
%
% By Emmanuel candes, 2003-2004

        n = length(y)/2;
	J = log2(n);
	x = zeros(1,n);

%
%  Unfold Partition at Coarse Level.
%
        [index, window] = CoarseMeyerWindow(L-1,deg);
	l_index = reverse(n/2 + 1 - index);
	r_index = n/2 + index;
	x(l_index) = y(1:2^L).* reverse(window);
	x(r_index) = y((2^L+1):2^(L+1)).* window;

%
%  Loop to Unfold Partition for  j = L - 1, ..., J - 3.
%
	for j = L-1:(J-3),
	  dyadic_points = [2^j 2^(j+1)];
	  [index, window] = DetailMeyerWindow(dyadic_points,deg);
	  yy = y((2^(j+2)+1):(2^(j+3)));
	  m = length(yy)/2;
	  l_index = reverse(n/2 + 1 - index);
	  x(l_index) = x(l_index) + yy(1:m).* reverse(window);
	  r_index = n/2 + index;
	  x(r_index) = x(r_index) + yy((m+1):(2*m)).* window;
	end

%
%  Finest Subband (for j = J - 2).
%

        j = J - 2;
        [index, window] = FineMeyerWindow(j,deg);
	yy = y((2^(j+2)+1):(2^(j+3)));
	m = length(yy)/2;
	l_index = reverse(n/2 + 1 - index);
        x(l_index) = x(l_index) + yy(1:m).* reverse(window);
        r_index = n/2 + index;
	x(r_index) = x(r_index) + yy((m+1):(2*m)).* window;

author: thomaskuestner


text:
function y = callback_fft(x,dir,options)

% callback_fft - callback for sparsity with FFT
%
%   y = callback_fft(x,dir,options);
%
%   Works in 1D and 2D. Orthogonal transforms.
%
%   Copyright (c) 2008 Gabriel Peyre

options.null = 0;

%% Detect dimension
if size(x,1)==1 || size(x,2)==1
    ndims = 1;
else
    ndims = 2;
end
ndims = getoptions(options, 'ndims', ndims);
isreal = getoptions(options, 'isreal', +1);
remove_high_freq = getoptions(options, 'remove_high_freq', 0);

%% Transform
if ndims==1
    %% 1D
    n = length(x);
    if dir==1
        y = ifft(x) * sqrt(n);
        if isreal
            y = real(y);
        end
    else
        y = fft(x) / sqrt(n);
    end
else
    n = size(x,1); p = size(x,2);
    %% 2D
    if dir==1
        if remove_high_freq>0
            x = remove_freq(x,remove_high_freq,dir);
        end
        y = ifft2(x) * sqrt(n*p);
        if isreal
            y = real(y);
        end
    else
        y = fft2(x) / sqrt(n*p);
        if remove_high_freq>0
            y = remove_freq(y,remove_high_freq,dir);
        end
    end
end

%%
function y = remove_freq(y,s,dir)

%y = fftshift(y);
% y(1:s,:) = 0; y(:,1:s) = 0;
% y(end-s+2:end,:) = 0; y(:,end-s+2:end) = 0;
% y = fftshift(y);

% remove only corners
y([end-s+2:end 1:s], end/2-s:end/2+s) = 0;
y(end/2-s:end/2+s, [end-s+2:end 1:s]) = 0;
y(end/2-s:end/2+s, end/2-s:end/2+s) = 0;

return;

global Qshuf;
global Ishuf;
if size(Qshuf)~=3*(2*s-1)^2
    A = zeros(size(y));
    n = size(y,1);
    s1 = n/2-s+1:n/2+s-1;
    s2 = [n-s+2:n 1:s];
    A(s2, s1) = 1;
    A(s1, s2) = 1;
    A(s1,s1) = 1;
    Ishuf = find(A==1);
    rand('state', 123456);
    [Qshuf,R] = qr(rand(length(Ishuf)));
end
if dir==1
    y(Ishuf) = Qshuf*y(Ishuf);
else
    y(Ishuf) = Qshuf'*y(Ishuf);
end

return;


text:
function nodes2 = grNeighborNodes(edges, node)
%GRNEIGHBORNODES Find adjacent nodes of a given node
%
%   Deprecated, use 'grAdjacentNodes' instead.
%
%   NEIGHS = grNeighborNodes(EDGES, NODE)
%   EDGES: the complete edges list
%   NODE: index of the node
%   NEIGHS: the nodes adjacent to the given node.
%
%   NODE can also be a vector of node indices, in this case the result is
%   the set of neighbors of any input node.
%
%   see also
%   grAdjacentNodes

%   HISTORY
%   10/02/2004 documentation
%   13/07/2004 faster algorithm
%   03/10/2007 can specify several input nodes

warning('MatGeom:deprecated', ...
    'function grNeighborNodes is obsolete, use grAdjacentNodes instead');

[i, j] = find(ismember(edges, node)); %#ok<NASGU>
nodes2 = edges(i,1:2);
nodes2 = unique(nodes2(:));
nodes2 = sort(nodes2(~ismember(nodes2, node)));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function pass = test_real( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e4*pref.techPrefs.chebfuneps;

% Example 1
f = real(ballfun(@(x,y,z)x+1i*y));
exact = ballfun(@(x,y,z)x);
pass(1) = norm( f - exact ) < tol;

% Example 2
f = real(ballfun(@(x,y,z)sin(z)+1i*cos(y)));
exact = ballfun(@(x,y,z)sin(z));
pass(2) = norm( f - exact ) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun