text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script is used for processing the raw data collectd by TI awr1843
% radar
% Author : Xiangyu Gao (xygao@uw.edu), University of Washingyton
% Input: raw I-Q radar data
% Output: range-angle (RA) image, 3D point clouds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clearvars
close all

% parameter setting
params = get_params_value();
% constant parameters
c = params.c; % Speed of light in air (m/s)
fc = params.fc; % Center frequency (Hz)
lambda = params.lambda;
Rx = params.Rx;
Tx = params.Tx;

% configuration parameters
Fs = params.Fs;
sweepSlope = params.sweepSlope;
samples = params.samples;
loop = params.loop;

Tc = params.Tc; % us
fft_Rang = params.fft_Rang;
fft_Vel = params.fft_Vel;
fft_Ang = params.fft_Ang;
num_crop = params.num_crop;
max_value = params.max_value; % normalization the maximum of data WITH 1843

% Creat grid table
rng_grid = params.rng_grid;
agl_grid = params.agl_grid;
vel_grid = params.vel_grid;

% Algorithm parameters
data_each_frame = samples*loop*Tx;
set_frame_number = 30;
frame_start = 1;
frame_end = set_frame_number;
Is_Windowed = 1;% 1==> Windowing before doing range and angle fft
Is_plot_rangeDop = 1;

% specify data name and load data as variable data_frames
seq_name = 'pms1000_30fs.mat';
seq_dir = strcat('.\template data\', seq_name);
load(seq_dir); % load data as variable data_frames

for i = frame_start:frame_end
    % read the data of each frame, and then arrange for each chirps
    data_frame = data_frames(:, (i-1)*data_each_frame+1:i*data_each_frame);
    data_chirp = [];
    for cj = 1:Tx*loop
        temp_data = data_frame(:, (cj-1)*samples+1:cj*samples);
        data_chirp(:,:,cj) = temp_data;
    end

    % separate the odd-index chirps and even-index chirps for TDM-MIMO with 2 TXs
    chirp_odd = data_chirp(:,:,1:2:end);
    chirp_even = data_chirp(:,:,2:2:end);

    % permutation with the format [samples, Rx, chirp]
    chirp_odd = permute(chirp_odd, [2,1,3]);
    chirp_even = permute(chirp_even, [2,1,3]);

    % Range FFT for odd chirps
    [Rangedata_odd] = fft_range(chirp_odd,fft_Rang,Is_Windowed);

    % Range FFT for even chirps
    [Rangedata_even] = fft_range(chirp_even,fft_Rang,Is_Windowed);

    % Doppler FFT
    Dopplerdata_odd = fft_doppler(Rangedata_odd, fft_Vel, 0);
    Dopplerdata_even = fft_doppler(Rangedata_even, fft_Vel, 0);
    Dopdata_sum = squeeze(mean(abs(Dopplerdata_odd), 2));

    % Plot range-Doppler image
    if Is_plot_rangeDop
        plot_rangeDop(Dopdata_sum,rng_grid,vel_grid);
    end

    % CFAR detector on Range-Velocity to detect targets
    % Output format: [doppler index, range index(start from index 1), ...
    % cell power]
    Pfa = 1e-4; % probability of false alarm
    [Resl_indx] = cfar_RV(Dopdata_sum, fft_Rang, num_crop, Pfa);
    detout = peakGrouping(Resl_indx);

    % doppler compensation on Rangedata_even using the max-intensity peak
    % on each range bin
    for ri = num_crop+1:fft_Rang-num_crop
        find_idx = find(detout(2, :) == ri);
        if isempty(find_idx)
            continue
        else
            % pick the first larger velocity
            pick_idx = find_idx(1);
            % phase compensation for virtual elements
            pha_comp_term = exp(-1i * pi * (detout(1,pick_idx)-fft_Vel/2-1) / fft_Vel);
            Rangedata_even(ri, :, :) = Rangedata_even(ri, :, :) * pha_comp_term;
        end
    end

    Rangedata_merge = [Rangedata_odd, Rangedata_even];

    % Angle FFT
    Angdata = fft_angle(Rangedata_merge,fft_Ang,Is_Windowed);
    Angdata_crop = Angdata(num_crop + 1:fft_Rang - num_crop, :, :);
    [Angdata_crop] = Normalize(Angdata_crop, max_value);

    % Angle estimation for detected point clouds
    Dopplerdata_merge = permute([Dopplerdata_odd, Dopplerdata_even], [2, 1, 3]);
    [Resel_agl, ~, rng_excd_list] = angle_estim_dets(detout, Dopplerdata_merge, fft_Vel, ...
        fft_Ang, Rx, Tx, num_crop);

    % Transform bin index to range/velocity/angle
    Resel_agl_deg = agl_grid(1, Resel_agl)';
    Resel_vel = vel_grid(detout(1,:), 1);
    Resel_rng = rng_grid(detout(2,:), 1);

    % save_det data format below
    % [range bin, velocity bin, angle bin, power, range(m), velocity (m/s), angle(degree)]
    save_det_data = [detout(2,:)', detout(1,:)', Resel_agl', detout(3,:)', ...
        Resel_rng, Resel_vel, Resel_agl_deg];

    % filter out the points with range_bin within the crop region
    if ~isempty(rng_excd_list)
        save_det_data(rng_excd_list, :) = [];
    end

    % Plot range-angle (RA) image
    plot_rangeAng(Angdata_crop,rng_grid(num_crop+1:fft_Rang-num_crop),agl_grid);

    % Point obtained clouds
    plot_pointclouds(save_det_data);

    break

end

author: Xiangyu-Gao


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function [X_clean, bag_data] = cleanLiDARTarget(scan_num, tag_num, bag_data, X, target_len, opt)
    N = 4; % clean up using N std for x axis
    M = 3; % clean up using M std for y, z axis
    opt = optimizeCost(opt, X, target_len, 0.001);
    X_ref = opt.H_opt * X;
%     distance = sum(X_ref(1:3,:), 1); % L1
    L_infinity = max(abs(X_ref(1:3,:)));  % L infinity
    K = find(L_infinity < (target_len/2)*1.025);
    X_ref_clean_yz = X_ref(:, K); % clean up y and z axis
    L_infinity = max(X_ref_clean_yz(1:3,:));
    K_center = find(L_infinity < target_len/4);
    X_std = std(X_ref_clean_yz(:, K_center), 1, 2);
    Q = find(abs(X_ref_clean_yz(1,:)) < N*(X_std(1))); % clean up x with 2 std
    X_ref_clean = X_ref_clean_yz(:, Q);
    X_clean = inv(opt.H_opt) * X_ref_clean;
    bag_data.lidar_target(tag_num).scan(scan_num).clean_up.std = N*(X_std(1));
    bag_data.lidar_target(tag_num).scan(scan_num).clean_up.L_infinity = L_infinity;
    bag_data.lidar_target(tag_num).scan(scan_num).clean_up.L_1 = sum(X_ref(1:3,:), 1);
%     figure(200);
%     clf('reset')
%     scatter3(X_ref(1,:), X_ref(2,:), X_ref(3,:))
%     xlabel('x')
%     ylabel('y')
%     zlabel('z')
%     axis equal
%     hold on;
%     scatter3(X_ref_clean_yz(1,:), X_ref_clean_yz(2,:), X_ref_clean_yz(3,:))
%     scatter3(X_ref_clean(1,:), X_ref_clean(2,:), X_ref_clean(3,:))
%     axis equal
%     view(90,0)
%     hold off;
end

author: UMich-BipedLab


text:
function [A,B,C,D] = vibsHmxBlockOperator(omega,U,sigma,u,cL,cT,rhoS,c0,rho0,f,tol)
%+========================================================================+
%|                                                                        |
%|                 OPENVIBS - LIBRARY FOR VIBRO-ACOUSTIC                  |
%|           openVibs is part of the GYPSILAB toolbox for Matlab          |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal, Marc Bakry (c) 2017-2019.                 |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             marc.bakry@polytechnique.edu                               |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : vibsHmxBlockOperator.m                        |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Marc Bakry                  |
%|  ( # )  |   CREATION   : 14.03.2019                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+

% Constants
mu     = rhoS*cT^2;
lambda = rhoS*(cL^2 - 2*cT^2);
w      = 2*pi*f;
k      = w/c0;

% Dimension
n = size(omega.msh.elt,2)-1;

% Green kernel function
if (n == 2)
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[H0(kr)]',k);
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]1',k);
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]2',k);
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]3',k);
    G0          = '[log(r)]';
    gradyG0     = 'grady[log(r)]';
    cteGxy      = 1i/4;
    cteG0       = -1/(2*pi);

elseif (n == 3)
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[exp(ikr)/r]',k);
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',k);
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',k);
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',k);
    G0          = '[1/r]';
    gradyG0     = 'grady[1/r]';
    cteGxy      =  1/(4*pi);
    cteG0       =  1/(4*pi);

else
    error('vibsNeumannBW.m : unavailable case.')
end

% Coupling coeff for Brackage-Werner simulation
beta = 1i*k;

%%%%%%%%%%%%%%%%%%%%%%%%%%% ELASTO-ELASTO (A11) %%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialization
A = cell(n,n);

% Static part
GG = integral(omega,grad(U),grad(U));
for i = 1:n
    for j = 1:n
        % Operator div(U):div(U)
        DD = integral(omega,grad(U,i),grad(U,j));

        % Operator e(U):e(U)
        EE = integral(omega,grad(U,j),grad(U,i));
        if (i==j)
            EE = EE + GG;
        end

        % Summation
        A{i,j} = lambda.*DD + mu.*EE;
    end
end

% Dynamic part
if (f ~= 0)
    Id = integral(omega,U,U);
    for i = 1:n
       A{i,i} = A{i,i} - (rhoS*w^2) .* Id;
    end
end

% Final form (sparse)
A = cell2mat(A);

%%%%%%%%%%%%%%%%%%%%%%%%%% ELASTO-ACOUSTIC (A12) %%%%%%%%%%%%%%%%%%%%%%%%%%

% Gaussian quadrature
[Xqud,Wx] = sigma.qud;
Wx        = spdiags(Wx,0,length(Wx),length(Wx));

% Collocation mass operator
Id = u.uqm(sigma);

% Collocation boundary operator
S = cteGxy .* integral(Xqud,sigma,Gxy,u,tol) + ...
    cteG0  .* regularize(Xqud,sigma,G0,u);

% Collocation boundary operator
D = cteGxy .* integral(Xqud,sigma,gradyGxy,ntimes(u),tol) + ...
    cteG0  .* regularize(Xqud,sigma,gradyG0,ntimes(u));

% Final operator Brackage-Werner : [1i*k*beta*S - (Id/2 + D)]
B.Mr = beta.*S - (0.5*Id + D);

% Normal trace of the volumn element matrix
nU   = ntimes(U);
nPHI = nU.uqm(sigma);

% Coupling to FEM
B.Ml = cell(n,1);
for i = 1:n
    B.Ml{i} = (nPHI{i}' * Wx);
end
B.Ml = cell2mat(B.Ml);

%%%%%%%%%%%%%%%%%%%%%%%%%% ELASTO-ACOUSTIC (A21) %%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialization
C = cell(1,n);

% Coupling FEM
for i = 1:n
    C{i} = (rho0*w^2) .* integral(sigma,u,ntimes(U,i));
end

% Final format (sparse)
C = cell2mat(C);

%%%%%%%%%%%%%%%%%%%%%%%%%% ACOUSTIC-ACOUSTIC (A22) %%%%%%%%%%%%%%%%%%%%%%%%

% Finite element mass matrix
Id = integral(sigma,u,u);

% Finite element boundary operator
H  = cteGxy .* (k^2 * integral(sigma,sigma,ntimes(u),Gxy,ntimes(u),tol) ...
    - integral(sigma,sigma,nxgrad(u),Gxy,nxgrad(u),tol));
Hr = cteG0  .* (k^2 * regularize(sigma,sigma,ntimes(u),G0,ntimes(u)) ...
    - regularize(sigma,sigma,nxgrad(u),G0,nxgrad(u)));

% Finite element boundary operator
D  = cteGxy .* integral(sigma,sigma,u,gradyGxy,ntimes(u),tol);
Dr = cteG0  .* regularize(sigma,sigma,u,gradyG0,ntimes(u));

% Final operator Brackage-Werner : - [1i*k*beta*(-Id/2 + Dt) - H]
D = - (beta.*(-0.5*Id + (D+Dr).') - (H+Hr));
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function pde = Darcydata4
%% Darcydata4 diagonal and anisotropic tensor
%
% p = x^2 + y^2
% K = [1 + 4*(x^2 + y^2)   0 ]
%   = [0   1 + 11*(x^2 + y^2)]

pde = struct('exactp',@exactp,'K',@K,'exactu',@exactu,'Dp',@gradp,...
             'f', @f,'g_D',@g_D,'g_N',@g_N);

    function s = K(pt)
       x = pt(:,1); y = pt(:,2);
%        s(:,3) = 3*x.*y;
       s(:,1) = 1 + 4*(x.^2 + y.^2);
       s(:,2) = 1 + 11*(x.^2 + y.^2);
    end
    function s = exactp(pt)
       x = pt(:,1); y = pt(:,2);
%        s = x.^3.*y + y.^4 + sin(pi*x).*cos(pi*y);
       s = x.^2 + y.^2;
    end
    function s = gradp(pt)
       x = pt(:,1); y = pt(:,2);
%        s(:,2) = x.^3 + 4*y.^3 - pi*sin(pi*x).*sin(pi*y);
%        s(:,1) = 3*x.^2.*y + pi*cos(pi*x).*cos(pi*y);
       s(:,2) = 2*y;
       s(:,1) = 2*x;
    end
    function s = exactu(pt)
       Dp = gradp(pt);  % u = K*grad(p)
       K = K(pt);
       s(:,2) = K(:,2).*Dp(:,2);
       s(:,1) = K(:,1).*Dp(:,1);
    end
    function s = f(pt)
        x = pt(:,1); y = pt(:,2);
%         s = 2+24*x.^2+14*y.^2+6*x.^2+2+22*x.^2+66*y.^2;
        s = 2+24*x.^2+16*y.^2++2+22*x.^2+66*y.^2;
        s = -s;
    end
    function s = g_D(pt)
       s = exactp(pt);
    end
    function f = g_N(p,vargin)
        if nargin > 1
            f = dot(exactu(p),vargin,2);
        else
            f = zeros(size(p,1),1);
            x = p(:,1);
            y = p(:,2);
            uprime = exactu(p);
            leftbd = (abs(x)<eps);  % n = (-1,0);
            f(leftbd) = - uprime(leftbd,1);
            rightbd = (abs(x-1)<eps); % n = (1,0);
            f(rightbd) = uprime(rightbd,1);
            topbd = (abs(y-1)<eps);   % n = (0,1)
            f(topbd) = uprime(topbd,2);
            bottombd = (abs(y)<eps);% n = (0,-1)
            f(bottombd) = - uprime(bottombd,2);
        end
    end
end

repo: ifem


text:
% Type information on a property.
% This class encapsulates type, shape and domain constraints on property
% values.
%
% Examples:
% PropertyType('denserealsingle', 'scalar')
%    creates a non-sparse single real scalar with unrestricted value
% PropertyType('denserealdouble', 'scalar', [-1,1])
%    creates a non-sparse double real scalar with value in range -1 to 1
% PropertyType('char', 'row', {'spring','summer','fall','winter'})
%    creates a character array which may be either be 'spring', 'summer',
%    'fall' or 'winter'
% PropertyType('logical', 'column', {'A','B','C'})
%    creates a set whose elements are 'A', 'B' and 'C'; the logical vector
%    is an indexer into this set, e.g. [1 1 0] maps to the set {'A','B'}
%
% See also: PropertyGridField

% Copyright 2010 Levente Hunyadi
classdef PropertyType
    properties
        % The underlying MatLab type for the property.
        PrimitiveType = 'denserealdouble';
        % The expected dimensions for the property.
        Shape = 'scalar';
        % A cell array of possible property values (or {} if unrestricted).
        Domain = [];
    end
    properties (Access = protected)
        ObjectType = '';
    end
    methods
        function self = PropertyType(type, shape, domain)
            self.PrimitiveType = type;
            self.Shape = shape;
            if nargin > 2
                self.Domain = domain;
            end
        end

        function s = char(self)
        % Compact textual representation of property type information.
            if ~isempty(self.Shape)
                s = sprintf('%s/%s', self.PrimitiveType, self.Shape);
            else
                s = self.PrimitiveType;
            end
            if ~isempty(self.Domain)
                if iscellstr(self.Domain)
                    s = [s ' ' strjoin(',', self.Domain)];
                elseif isnumeric(self.Domain)
                    s = [s ' ' mat2str(self.Domain)];
                end
            end
        end

        function tf = islogical(self)
            tf = strcmp(self.PrimitiveType, 'logical');
        end

        function tf = isinteger(self)
            tf = any(strmatch(self.PrimitiveType, {'int8','uint8','int16','uint16','int32','uint32','int64','uint64'}, 'exact'));
        end

        function tf = isfloat(self)
            tf = any(strmatch(self.PrimitiveType, {'denserealsingle','denserealdouble','densecomplexsingle','densecomplexdouble','sparserealsingle','sparserealdouble','sparsecomplexsingle','sparsecomplexdouble'}, 'exact'));
        end

        function tf = isnumeric(self)
            tf = isinteger(self) || isfloat(self);
        end

        function tf = isreal(self)
            tf = any(strmatch(self.PrimitiveType, {'denserealsingle','denserealdouble','sparserealsingle','sparserealdouble'}, 'exact'));
        end

        function tf = is2d(self)
            tf = any(strmatch(self.Shape, {'row','column','matrix'}, 'exact'));
        end

        function tf = isvector(self)
            tf = any(strmatch(self.Shape, {'row','column'}, 'exact'));
        end

        function tf = isstring(self)
            tf = strcmp(self.PrimitiveType, 'char') && strcmp(self.Shape, 'row');
        end

        function tf = isset(self)
            tf = islogical(self) && isvector(self) && ~isempty(self.Domain);
        end

        function self = set.PrimitiveType(self, type)
            validateattributes(type, {'char'}, {'nonempty','row'});
            self.PrimitiveType = validatestring(type, { ...
                'logical', ...
                'char', ...
                'int8','uint8','int16','uint16','int32','uint32','int64','uint64', ...
                'denserealsingle','denserealdouble', ...
                'densecomplexsingle','densecomplexdouble', ...
                'sparserealsingle','sparserealdouble', ...
                'sparsecomplexsingle','sparsecomplexdouble', ...
                'cellstr', ...
                'object'});
        end

        function self = set.Shape(self, shape)
            validateattributes(shape, {'char'}, {'nonempty','row'});
            self.Shape = validatestring(shape, {'scalar','row','column','matrix','empty'});  % 'empty' is a pseudo-shape for the 0-by-0 empty matrix
        end

        function self = set.Domain(self, domain)
            if isempty(domain)
                self.Domain = [];
                return;
            end
            if ~isscalar(self) && ~isstring(self) && ~(islogical(self) && isvector(self))
                error('PropertyType:InvalidArgumentValue', ...
                    'Domain can be set only for scalars, strings or logical vectors.');
            end
            if iscell(domain)
                validateattributes(domain, {'cell'}, {'vector'});
                if ~(islogical(self) && isvector(self))  % interpret as exhaustive enumeration of domain elements
                    for k = 1 : numel(domain)
                        assert(self.CanAccept(domain{k}), ...
                            'PropertyType:InvalidArgumentValue', ...
                            'Invalid value specified for property domain.');
                    end
                end
            elseif isnumeric(domain)  % interpret as numeric range for domain elements
                validateattributes(domain, {'numeric'}, {'vector'});
                assert(isnumeric(self), 'PropertyType:ArgumentTypeMismatch', ...
                    'Closed interval domain is meaningful only for numeric data.');
                assert(length(domain) == 2, 'PropertyType:InvalidArgumentValue', ...
                    'Closed interval data is expected in the form [a,b] where a and b are numbers.');
                domain = [min(domain) max(domain)];
            else
                error('PropertyType:InvalidArgumentValue', ...
                    'Closed interval specification [a,b] or explicit enumeration {a,b,c,...} is expected.');
            end
            self.Domain = domain;
        end

        function tf = CanAccept(self, value, precisionloss)
        % Whether the property can accept the given value.
        %
        % Output arguments:
        % tf:
        %    false if any of the property constaints (type, shape or
        %    domain) would be violated
            validateattributes(self, {'PropertyType'}, {'scalar'});
            if nargin < 3
                precisionloss = false;
            end
            if isobject(value)
                tf = ~isempty(self.ObjectType) && metaclass(value) <= meta.class.fromName(self.ObjectType);
            elseif (isnumeric(value) || islogical(value) || ischar(value)) && any(isnan(value(:)));
                tf = false;
            else
                type = PropertyType.AutoDiscoverType(value);
                shape = PropertyType.AutoDiscoverShape(value);
                tf = PropertyType.CanAcceptType(self.PrimitiveType, type, precisionloss) ...
                    && PropertyType.CanAcceptShape(self.Shape, shape) ...
                    && self.CanAcceptValue(value);
            end
        end

        function clazz = GetPrimitiveMatLabType(self)
            validateattributes(self, {'PropertyType'}, {'scalar'});
            switch (self.PrimitiveType)
                case {'logical','char','int8','uint8','int16','uint16','int32','uint32','int64','uint64'}
                    clazz = self.PrimitiveType;
                case {'denserealsingle','densecomplexsingle','sparserealsingle','sparsecomplexsingle'}
                    clazz = 'single';
                case {'denserealdouble','densecomplexdouble','sparserealdouble','sparsecomplexdouble'}
                    clazz = 'double';
                case 'cellstr'
                    clazz = 'cell';
            end
        end

        function javatype = GetPrimitiveJavaType(self)
            validateattributes(self, {'PropertyType'}, {'scalar'});
            switch self.PrimitiveType
                case {'denserealdouble','sparserealdouble'}  % add a double-precision floating point property
                    matlabtype = 'double';  % MatLab type that is marshalled to Java
                case {'denserealsingle','sparserealsingle'}  % add a single-precision floating point property
                    matlabtype = 'single';
                case {'int8','uint8','int16'}
                    matlabtype = 'int16';
                case {'uint16','int32'}  % add an integer property
                    matlabtype = 'int32';
                case {'uint32','int64'}
                    matlabtype = 'int64';
                case 'logical'  % add a logical property
                    matlabtype = 'logical';
                case {'densecomplexdouble','sparsecomplexdouble','densecomplexsingle','sparsecomplexsingle'}
                    matlabtype = 'cellstr';
                otherwise
                    error('MatrixEditor:ArgumentTypeMismatch', 'Type %s is not supported.', self.PrimitiveType);
            end
            javatype = javaclass(matlabtype);
        end

        function javavalue = GetPrimitiveJavaValue(self, value)
            validateattributes(self, {'PropertyType'}, {'scalar'});
            switch self.PrimitiveType
                case {'denserealdouble','sparserealdouble'}  % add a double-precision floating point property
                    javavalue = java.lang.Double(value);
                case {'denserealsingle','sparserealsingle'}  % add a single-precision floating point property
                    javavalue = java.lang.Float(value);
                case {'int8','uint8','int16'}
                    javavalue = java.lang.Short(value);
                case {'uint16','int32'}  % add an integer property
                    javavalue = java.lang.Integer(value);
                case {'uint32','int64'}
                    javavalue = java.lang.Long(value);
                case 'logical'  % add a logical property
                    javavalue = java.lang.Boolean(value);
                case {'densecomplexdouble','sparsecomplexdouble','densecomplexsingle','sparsecomplexsingle'}
                    javavalue = java.lang.String(mat2str(value));
                otherwise
                    if isnumeric(value) || islogical(value)
                        javavalue = mat2str(value);
                    else
                        javavalue = char(value);
                    end
                    %error('PropertyType:ArgumentTypeMismatch', 'Type %s is not supported.', self.PrimitiveType);
            end
        end

        function javamatrix = GetJavaVectorOfVectors(self, matrix)
            assert(isnumeric(self) || islogical(self), 'PropertyType:InvalidOperation', ...
                'Operation supported on numeric and logical matrices only.');
            javamatrix = java.util.Vector(int32(size(matrix,1)));
            for m = 1 : size(matrix,1)
                javarow = java.util.Vector(int32(size(matrix,2)));
                javamatrix.add(javarow);
                for n = 1 : size(matrix,2)
                    javarow.add(self.GetPrimitiveJavaValue(matrix(m,n)));
                end
            end
        end

        function javamatrix = GetJavaMatrix(self, matrix)
        % Get MatLab value as Java 2d array of proper type.
            assert(isnumeric(self), 'PropertyType:InvalidOperation', ...
                'Operation supported on numeric matrices only.');
            assert(numel(matrix) > 0, 'PropertyType:InvalidOperation', ...
                'Empty values are not supported.');
            javatype = self.GetPrimitiveJavaType();
            javamatrix = javaArray(char(javatype.getName()), size(matrix,1), size(matrix,2));
            for m = 1 : size(matrix,1)
                for n = 1 : size(matrix,2)
                    javamatrix(m,n) = self.GetPrimitiveJavaValue(matrix(m,n));
                end
            end
        end

        function javavector = GetJavaVector(self, vector)
        % Get MatLab value as a Java vector of proper type.
            assert(isnumeric(self), 'PropertyType:InvalidOperation', ...
                'Operation supported on numeric matrices only.');
            assert(numel(vector) > 0, 'PropertyType:InvalidOperation', ...
                'Empty values are not supported.');
            javatype = self.GetPrimitiveJavaType();
            javavector = javaArray(char(javatype.getName()), numel(vector));
            for k = 1 : numel(vector)
                javavector(k) = self.GetPrimitiveJavaValue(vector(k));
            end
        end

        function javavalue = ConvertToJava(self, value)
        % Converts a MatLab value into the appropriate Java object.
            switch self.Shape
                case 'scalar'
                    javavalue = self.GetPrimitiveJavaValue(value);
                case {'row','column'}
                    switch self.PrimitiveType
                        case 'char'  % add a string property
                            if strcmp(self.Shape, 'row')
                                javavalue = value;
                            else
                                javavalue = mat2str(value);
                            end
                        case 'cellstr'
                            javavalue = java.lang.String(strjoin(sprintf('\n'), value));
                        case 'logical'
                            if ~isempty(self.Domain)
                                javavalue = javaStringArray(self.Domain(value));  % value is an indicator vector
                            else
                                javavalue = mat2str(value);
                            end
                        otherwise
                            if isnumeric(value) || islogical(value)
                                javavalue = mat2str(value);
                            else
                                javavalue = char(value);
                            end
                    end
                case 'empty'
                    javavalue = [];
                otherwise
                    if isnumeric(value) || islogical(value)
                        javavalue = mat2str(value);
                    else
                        javavalue = char(value);
                    end
            end
        end

        function [value,stringconversion] = ConvertFromJava(self, javavalue)
        % Converts a Java object into the appropriate MatLab value.
            stringconversion = false;
            switch self.Shape
                case 'scalar'
                    switch self.PrimitiveType
                        case 'denserealdouble'
                            value = full(double(javavalue));
                        case 'sparserealdouble'
                            value = sparse(double(javavalue));
                        case 'denserealsingle'
                            value = full(single(javavalue));
                        case 'sparserealsingle'
                            value = sparse(single(javavalue));
                        case {'int8','uint8','int16','uint16','int32','uint32','int64','logical'}
                            value = cast(javavalue, self.PrimitiveType);
                        case {'densecomplexdouble','sparsecomplexdouble','densecomplexsingle','sparsecomplexsingle'}
                            value = self.ConvertFromString(javavalue);
                            stringconversion = true;
                    end
                case {'row','column'}
                    switch self.PrimitiveType
                        case 'char'  % add a string property
                            if strcmp(self.Shape, 'row')
                                value = char(javavalue);
                            else
                                value = self.ConvertFromString(javavalue);
                                stringconversion = true;
                                return;
                            end
                        case 'cellstr'
                            value = strsplit(javavalue);
                        case 'logical'
                            if ~isempty(self.Domain)
                                value = strsetmatch(cell(javavalue), self.Domain);
                            else
                                value = self.ConvertFromString(javavalue);
                                stringconversion = true;
                                return;
                            end
                        case {...
                                'denserealdouble','sparserealdouble','denserealsingle','sparserealsingle',...
                                'densecomplexdouble','sparsecomplexdouble','densecomplexsingle','sparsecomplexsingle',...
                                'int8','uint8','int16','uint16','int32','uint32','int64','uint64'}
                            value = self.ConvertFromString(javavalue);
                            stringconversion = true;
                    end
                    switch self.Shape
                        case 'row'
                            value = reshape(value, 1, numel(value));
                        case 'column'
                            value = reshape(value, numel(value), 1);
                    end
                case 'matrix'
                    switch self.PrimitiveType
                        case {...
                                'denserealdouble','sparserealdouble','denserealsingle','sparserealsingle',...
                                'densecomplexdouble','sparsecomplexdouble','densecomplexsingle','sparsecomplexsingle',...
                                'char','int8','uint8','int16','uint16','int32','uint32','int64','uint64','logical'}
                            value = self.ConvertFromString(javavalue);
                            stringconversion = true;
                    end
            end
        end

        function value = ConvertFromString(self, text)
        % Converts matrix string representation to its numeric value.
        %
        % Input arguments:
        % text:
        %    a textual representation of a numeric scalar, vector or matrix
        %
        % Output arguments:
        % value:
        %    a numeric value satisfying all type constraints, or NaN if the
        %    string cannot be converted or the value cannot be coerced into
        %    the appropriate type
            text = char(text);
            switch self.Shape
                case 'scalar'
                    value = str2double(text);
                    status = ~isnan(value);
                otherwise
                    [value,status] = str2num(text); %#ok<ST2NM>
            end
            if all(int32(value(:)) == value(:))  % all numbers are integers
                value = int32(value);            % cast to integer
            end
            if ~status || ~self.CanAccept(value, true)
                value = NaN;
                return;
            end
            value = self.ConvertFromMatLab(value);
        end

        function value = ConvertFromMatLab(self, value)
        % Casts a MatLab value to the prescribed property type.
            if ~isreal(value) && isnumeric(self) && isreal(self)
                error('PropertyType:ArgumentTypeMismatch', 'Complex value cannot be coerced into a real type.');
            end
            if iscellstr(value) && strcmp(self.PrimitiveType, 'cellstr') || isobject(value) && strcmp(self.PrimitiveType, 'object')
                return
            end
            switch self.PrimitiveType
                case {'denserealdouble','densecomplexdouble'}
                    value = full(double(value));
                case {'sparserealdouble','sparsecomplexdouble'}
                    value = sparse(double(value));
                case {'denserealsingle','densecomplexsingle'}
                    value = full(single(value));
                case {'sparserealsingle','sparsecomplexsingle'}
                    value = sparse(single(value));
                case {'int8','uint8','int16','uint16','int32','uint32','int64','logical'}
                    value = cast(value, self.PrimitiveType);
                case {'char'}
                    value = char(value);
                otherwise
                    error('PropertyType:ArgumentTypeMismatch', 'Cannot coerce type %s into type %s.', class(value), self.PrimitiveType);
            end
        end
    end
    methods (Access = private)
        function tf = CanAcceptValue(self, value)
        % Whether the property constraints allow the value to be accepted.
            if islogical(self) && isvector(self)  % interpret logical vector as a set
                tf = isempty(self.Domain) || numel(value) == numel(self.Domain);  % size of logical vector corresponds to size of universe
            else  % test whether domain contains value
                tf = PropertyType.IsInDomain(self.Domain, value);
            end
        end
    end
    methods (Static)
        function obj = AutoDiscover(value)
        % Constructs a property type instance based on a value.
        % The property type is chosen to be the most specific possible that
        % fits the given value.
            if isobject(value)
                obj = PropertyType('object', PropertyType.AutoDiscoverShape(value));
                obj.ObjectType = class(value);
            else
                obj = PropertyType(PropertyType.AutoDiscoverType(value), PropertyType.AutoDiscoverShape(value));
            end
        end

        function type = AutoDiscoverType(value)
            clazz = class(value);
            switch clazz
                case {'logical','char','int8','uint8','int16','uint16','int32','uint32','int64','uint64'}
                    type = clazz;
                case {'single','double'}
                    if issparse(value)
                        sparsity = 'sparse';
                    else
                        sparsity = 'dense';
                    end
                    if isreal(value)
                        complexity = 'real';
                    else
                        complexity = 'complex';
                    end
                    type = [ sparsity complexity clazz ];
                case 'cell'
                    if iscellstr(value)
                        type = 'cellstr';
                    else
                        error('PropertyType:InvalidArgumentValue', ...
                            'Cell arrays other than cell array of strings are not supported.');
                    end
                otherwise
                    error('PropertyType:InvalidArgumentValue', ...
                        'Argument type "%s" is not supported.', class(value));
            end
        end

        function shape = AutoDiscoverShape(value)
            if ndims(value) > 2
                error('PropertyType:InvalidArgumentValue', ...
                    'Dimensions higher than 2 are not supported.');
            end
            if size(value,1) == 1 && size(value,2) == 1
                shape = 'scalar';
            elseif size(value,1) == 1
                shape = 'row';
            elseif size(value,2) == 1
                shape = 'column';
            elseif size(value,1) == 0 && size(value,2) == 0
                shape = 'empty';  % no dimensions
            else
                shape = 'matrix';
            end
        end

        function tf = CanAcceptType(generaltype, specifictype, precisionloss)
        % Type assignment check.
        %
        % Output arguments:
        % tf:
        %    true if a variable of the specific type could be assigned to a
        %    variable of the general type
            validateattributes(generaltype, {'char'}, {'nonempty','row'});
            validateattributes(specifictype, {'char'}, {'nonempty','row'});
            if nargin < 3
                precisionloss = false;
            end
            if strcmp(generaltype, specifictype)
                tf = true;
                return;
            end
            switch generaltype
                case {'densecomplexdouble','sparsecomplexdouble'}
                    type = {...
                        'int8','uint8','int16','uint16','int32','uint32','int64','uint64', ...
                        'denserealsingle','denserealdouble', ...
                        'densecomplexsingle','densecomplexdouble', ...
                        'sparserealsingle','sparserealdouble', ...
                        'sparsecomplexsingle','sparsecomplexdouble'};
                case {'denserealdouble','sparserealdouble'}
                    type = {...
                        'int8','uint8','int16','uint16','int32','uint32','int64','uint64', ...
                        'denserealsingle','denserealdouble', ...
                        'sparserealsingle','sparserealdouble'};
                case {'densecomplexsingle','sparsecomplexsingle'}
                    type = {...
                        'int8','uint8','int16','uint16','int32','uint32','int64','uint64', ...
                        'denserealsingle','densecomplexsingle','sparserealsingle','sparsecomplexsingle'};
                    if precisionloss  % permit loss of precision
                        type = [type {'denserealdouble','densecomplexdouble','sparserealdouble','sparsecomplexdouble'}];
                    end
                case {'denserealsingle','sparserealsingle'}
                    type = {...
                        'int8','uint8','int16','uint16','int32','uint32','int64','uint64', ...
                        'denserealsingle','sparserealsingle'};
                    if precisionloss  % permit loss of precision
                        type = [type {'denserealdouble','sparserealdouble'}];
                    end
                case {'int64','uint64'}
                    type = {'int8','uint8','int16','uint16','int32','uint32'};
                case {'int32','uint32'}
                    type = {'int8','uint8','int16','uint16'};
                case {'int16','uint16'}
                    type = {'int8','uint8'};
                case {'int8','uint8'}
                    type = {};
                case 'char'
                    type = {};
                case 'logical'
                    type = {'int8','uint8','int16','uint16','int32','uint32'};
                case 'cellstr'
                    type = {};
                otherwise
                    type = {};
            end
            tf = any(strmatch(specifictype, type, 'exact'));
        end

        function tf = CanAcceptShape(generalshape, specificshape)
        % Shape assignment check.
        %
        % Output arguments:
        % tf:
        %    true if a variable of the specific shape could be assigned to
        %    a variable of the general shape
            validateattributes(generalshape, {'char'}, {'nonempty','row'});
            validateattributes(specificshape, {'char'}, {'nonempty','row'});
            if strcmp(generalshape, specificshape)
                tf = true;
                return;
            end
            switch generalshape
                case 'matrix'
                    shape = {'scalar','row','column','empty'};
                case 'row'
                    shape = {'scalar','empty'};
                case 'column'
                    shape = {'scalar','empty'};
                case 'scalar'
                    shape = {};
                case 'empty'
                    shape = {};
                otherwise
                    shape = {};
            end
            tf = any(strmatch(specificshape, shape, 'exact'));
        end

        function tf = IsInDomain(domain, value)
            if isempty(domain)
                tf = true;
            elseif iscellstr(domain)
                tf = any(strcmp(value, domain));
            elseif iscell(domain)
                tf = any(cellfun(@(v) v==value, domain));
            elseif isnumeric(domain) && length(domain) == 2
                tf = value >= min(domain) && value <= max(domain);
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function ranks = tied_ranks(x,tol)
% like tiedrank(), but faster and without the bells & whistles.
%
%                               Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                               2013-07-22

if ~exist('tol','var')
    tol = 10*eps; end
ranks = 1:length(x);
[sorted,order] = sort(x);
tie_ranges = reshape(find(diff([false abs(diff(sorted(:)'))<tol false])),2,[])';
for r=1:size(tie_ranges,1)
    range = tie_ranges(r,1):tie_ranges(r,2);
    ranks(range) = sum(ranks(range))/(range(end)-range(1)+1);
end
ranks(order) = ranks;
ranks = reshape(ranks,size(x));

repo: REST


text:
function scimat = scimat_squeeze(scimat, todouble)
% SCIMAT_SQUEEZE  Remove dummy dimension from SCIRUNNRRD struct so that it
% becomes a SCIMAT struct.
%
% SCIMAT = scimat_squeeze(SCIRUNNRRD)
%
%   When the application Seg3D saves images to Matlab format, it saves them
%   to a SCIRUNNRRD struct, which is basically SCIMAT with a dummy
%   dimension.
%
%   SCIMAT is a struct that contains a 3D volume. We use SCIMAT structs
%   widely in Gerardus, because that way we have the image data and
%   metainformation (e.g. voxel size) together in the same
%   variable. For details on SCIMAT structs, see "help scimat".
%
%   This function removes the dummy dimension.
%
% SCIMAT = scimat_squeeze(..., TODOUBLE)
%
%   TODOUBLE is a flag to convert the image data from uint8 to double. This
%   is useful for image processing, but the data will be 8 times larger. By
%   default, TODOUBLE=false.
%
% See also: scimat_load.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010-2014 University of Oxford
% Version: 0.2.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 2);
nargoutchk(0, 1);

% defaults
if (nargin < 2 || isempty(todouble))
    todouble = false;
end

% remove dummy dimension
scimat.data = squeeze(scimat.data);
if (length(scimat.axis) > 3)
    scimat.axis = scimat.axis(2:end);
end

% convert data to double
if (todouble)
    scimat.data = double(scimat.data);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function value = r4_aie ( x )

%*****************************************************************************80
%
%% R4_AIE evaluates the exponential scaled Airy function Ai of an R4 argument.
%
%  Discussion:
%
%    If X <= 0
%      R4_AIE ( X ) = R4_AI ( X )
%    else
%      R4_AIE ( X ) = R4_AI ( X ) * exp ( 2/3 X^(3/2) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the Airy function Ai of X.
%
  persistent aifcs
  persistent aigcs
  persistent aipcs
  persistent naif
  persistent naig
  persistent naip
  persistent x32sml
  persistent x3sml
  persistent xbig

  if ( isempty ( naif ) )
    aifcs = [ ...
     -0.03797135849666999750E+00, ...
      0.05919188853726363857E+00, ...
      0.00098629280577279975E+00, ...
      0.00000684884381907656E+00, ...
      0.00000002594202596219E+00, ...
      0.00000000006176612774E+00, ...
      0.00000000000010092454E+00, ...
      0.00000000000000012014E+00, ...
      0.00000000000000000010E+00 ]';
    aigcs = [ ...
      0.01815236558116127E+00, ...
      0.02157256316601076E+00, ...
      0.00025678356987483E+00, ...
      0.00000142652141197E+00, ...
      0.00000000457211492E+00, ...
      0.00000000000952517E+00, ...
      0.00000000000001392E+00, ...
      0.00000000000000001E+00 ]';
    aipcs = [ ...
     -0.0187519297793868E+00, ...
     -0.0091443848250055E+00, ...
      0.0009010457337825E+00, ...
     -0.0001394184127221E+00, ...
      0.0000273815815785E+00, ...
     -0.0000062750421119E+00, ...
      0.0000016064844184E+00, ...
     -0.0000004476392158E+00, ...
      0.0000001334635874E+00, ...
     -0.0000000420735334E+00, ...
      0.0000000139021990E+00, ...
     -0.0000000047831848E+00, ...
      0.0000000017047897E+00, ...
     -0.0000000006268389E+00, ...
      0.0000000002369824E+00, ...
     -0.0000000000918641E+00, ...
      0.0000000000364278E+00, ...
     -0.0000000000147475E+00, ...
      0.0000000000060851E+00, ...
     -0.0000000000025552E+00, ...
      0.0000000000010906E+00, ...
     -0.0000000000004725E+00, ...
      0.0000000000002076E+00, ...
     -0.0000000000000924E+00, ...
      0.0000000000000417E+00, ...
     -0.0000000000000190E+00, ...
      0.0000000000000087E+00, ...
     -0.0000000000000040E+00, ...
      0.0000000000000019E+00, ...
     -0.0000000000000009E+00, ...
      0.0000000000000004E+00, ...
     -0.0000000000000002E+00, ...
      0.0000000000000001E+00, ...
     -0.0000000000000000E+00 ]';
    eta = 0.1 * r4_mach ( 3 );
    naif = r4_inits ( aifcs, 9, eta );
    naig = r4_inits ( aigcs, 8, eta );
    naip = r4_inits ( aipcs, 34, eta );
    x3sml = eta^0.3333;
    x32sml = 1.3104 * x3sml * x3sml;
    xbig = r4_mach ( 2 )^0.6666;
  end

  if ( x < - 1.0 )
    [ xm, theta ] = r4_aimp ( x );
    value = xm * cos ( theta );
  elseif ( abs ( x ) <= x32sml )
    z = 0.0;
    value = 0.375 + ( r4_csevl ( z, aifcs, naif ) ...
      - x * ( 0.25 + r4_csevl ( z, aigcs, naig ) ) );
  elseif ( abs ( x ) <= x3sml )
    z = 0.0;
    value = 0.375 + ( r4_csevl ( z, aifcs, naif ) ...
      - x * ( 0.25 + r4_csevl ( z, aigcs, naig ) ) );
    value = value * exp ( 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x <= 1.0 )
    z = x * x * x;
    value = 0.375 + ( r4_csevl ( z, aifcs, naif ) ...
      - x * ( 0.25 + r4_csevl ( z, aigcs, naig ) ) );
    value = value * exp ( 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x < xbig )
    sqrtx = sqrt ( x );
    z = 2.0 / ( x * sqrtx ) - 1.0;
    value = ( 0.28125 + r4_csevl ( z, aipcs, naip ) ) / sqrt ( sqrtx );
  else
    sqrtx = sqrt ( x );
    z = - 1.0;
    value = ( 0.28125 + r4_csevl ( z, aipcs, naip ) ) / sqrt ( sqrtx );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [success,msg,fps] = ...
  ReadFPS_Ctrax(varargin)

success = false;
msg = ''; %#ok<NASGU>

[fps,intrxfile,annfile] = myparse(varargin,...
  'fps',1,'intrxfile','','annfile','');

% try to read from the trx file first
if isempty(intrxfile),
  msg = 'Input trxfile not yet set';
  return;
end

try
  [trx,~,success1,timestamps] = load_tracks(intrxfile,'',...
    'dosave',false,'annname',annfile,'verbose',false);
  if ~success1,
    msg = sprintf('Could not load tracks from trxfile %s',intrxfile);
    return;
  end
  readfps = false;
  if isfield(trx,'fps'),
    fps = [trx.fps];
    if nnz(~isnan(fps)) > 0,
      fps = nanmedian(fps);
      readfps = true;
    end
  end
  if ~readfps,
    if nnz(~isnan(timestamps)) == 0,
      msg = sprintf('No timestamps read in from trxfile %s',intrxfile);
      return;
    end
    fps = 1/nanmedian(diff(timestamps));
  end
catch ME,
  msg = sprintf('Could not load from trxfile %s: %s',intrxfile,getReport(ME));
  return;
end

success = true;
msg = 'Read fps from trxfile';

author: kristinbranson

repo: JAABA


text:
% clear,clc
%
% close all
global bound rng
l1=1;l2=1;l3=0.5;

hh=findobj(gcf,'tag','x1');
xs=str2num(get(hh,'string'));
hh=findobj(gcf,'tag','y1');
ys=str2num(get(hh,'string'));
hh=findobj(gcf,'tag','phi');
phis=str2num(get(hh,'string'));
hh=findobj(gcf,'tag','x2');
xf=str2num(get(hh,'string'));
hh=findobj(gcf,'tag','y2');
yf=str2num(get(hh,'string'));

phi=phis*pi/180;
[k,confi]=invkini(xs,ys,phi);
if k > 1
    warndlg('non valid initial conditions..... change the initial conditions','!! error !!')
    clear
    return
end
if sqrt(xf^2+yf^2)==2.5
   warndlg('there is no redundancy..... change final position','!! error !!')
   clear
   return
end
para=[confi xf yf];
%GA parameters
pops=40;
crossprop=0.8;
mutprop=0.05;
maxgen=80;
bound=[ -pi     pi;% qm1
        -pi     pi;% qm2
        -pi     pi;% qm3
        -pi     pi;% phif
       -pi/4   pi/4;% vqm1
       -pi/4   pi/4;% vqm2
       -pi/4   pi/4;% vqm3
        0       8;% t1
        0       8];%t2

%initialization
numvar=size(bound,1);
rng=(bound(:,2)-bound(:,1))';
pop=zeros(pops,numvar);
%pop = initial population
pop(:,1:numvar)=(ones(pops,1)*rng).*(rand(pops,numvar))+...
    (ones(pops,1)*bound(:,1)');
tic
for it=1:maxgen
    [fpop,z]=fitnesstra3f(pop,para);
    pop(:,4)=z';
    [cs,inds]=max(fpop);bchrom=pop(inds,:);
    % tournament selection
    toursize=5;
    players=ceil(pops*rand(pops,toursize));
    scores=fpop(players);
    [a,m]=max(scores');
    pind=zeros(1,pops);
    for ii=1:pops
        pind(ii)=players(ii,m(ii));
        parent(ii,:)=pop(pind(ii),:);
    end
    %arithmatic crossover
    offs=cross_singlepoint(parent,crossprop);
%     offs=cross1(parent,crossprop);
    %mutate1=uniform mutation.
    moffs=mutate1(offs,mutprop);
    pop=moffs;
    [mm,z]=fitnesstra3f(pop,para);
    pop(:,4)=z';
    maxf(it)=max(mm);
    [bfit,bind]=max(mm);
    bsol=pop(bind,:);    % best solution.
    rec=recor(bsol,para);
    trec(it)=rec(1);
    qrec(it)=rec(2);
    drec(it)=rec(3);
    torrec(it)=rec(4);
    if mm(inds) < cs
    pop(inds,:)=bchrom;
    end
    %*************************plot results

    qs1=para(1);qs2=para(2);qs3=para(3);
    qm1=bsol(1);qm2=bsol(2);qm3=bsol(3);
    [w,conff]=invkin3(para(4),para(5),bsol(4));
    qg1=conff(1);qg2=conff(2);qg3=conff(3);

    [xs1,ys1]=pol2cart(qs1,l1);
    [xs2,ys2]=pol2cart(qs2+qs1,l2);
    [xs3,ys3]=pol2cart(qs3+qs2+qs1,l3);

    xxs1=linspace(0,xs1);
    yys1=linspace(0,ys1);
    xxs2=linspace(0,xs2);xxs2=xxs2+xs1;
    yys2=linspace(0,ys2);yys2=yys2+ys1;
    xxs3=linspace(0,xs3);xxs3=xxs3+xs2+xs1;
    yys3=linspace(0,ys3);yys3=yys3+ys2+ys1;

    [xm1,ym1]=pol2cart(qm1,l1);
    [xm2,ym2]=pol2cart(qm2+qm1,l2);
    [xm3,ym3]=pol2cart(qm3+qm2+qm1,l3);

    xxm1=linspace(0,xm1);
    yym1=linspace(0,ym1);
    xxm2=linspace(0,xm2);xxm2=xxm2+xm1;
    yym2=linspace(0,ym2);yym2=yym2+ym1;
    xxm3=linspace(0,xm3);xxm3=xxm3+xm2+xm1;
    yym3=linspace(0,ym3);yym3=yym3+ym2+ym1;

    [xg1,yg1]=pol2cart(qg1,l1);
    [xg2,yg2]=pol2cart(qg2+qg1,l2);
    [xg3,yg3]=pol2cart(qg3+qg2+qg1,l3);

    xxg1=linspace(0,xg1);
    yyg1=linspace(0,yg1);
    xxg2=linspace(0,xg2);xxg2=xxg2+xg1;
    yyg2=linspace(0,yg2);yyg2=yyg2+yg1;
    xxg3=linspace(0,xg3);xxg3=xxg3+xg2+xg1;
    yyg3=linspace(0,yg3);yyg3=yyg3+yg2+yg1;
    xt=[xxs1;xxs2;xxs3;xxm1;xxm2;xxm3;xxg1;xxg2;xxg3];
    yt=[yys1;yys2;yys3;yym1;yym2;yym3;yyg1;yyg2;yyg3];
    cond=[confi,conff];
    chrom=[bsol(1:3),bsol(5:end)];
    kk=trajt3(cond,chrom);
    pq=kk(1:3,[5,7,10,13,15,17,20,23,25,27,30,33,35]);
    pcart=forkin3(kk(1,:),kk(2,:),kk(3,:));
    px=pcart(1,:);py=pcart(2,:);
    if it==maxgen
      [xxt,yyt]=angls2links(pq);
      figure,plot(xxt',yyt')
      xlabel('x(m)')
      ylabel('y(m)')
      hold on
    end
    plot(px,py)
    hold on
    axis([-2.7 2.7 -2.7 2.7])
    text(0,2.6,['gen. no.',num2str(it)])
    plot(xt',yt')
    hold off
    pause(0)
end
toc
e=[1:maxgen];
figure,plot(e,1./maxf)
xlabel('generation')
ylabel('min. fitness')

tt=torque3(kk);

t1=bsol(8);ti1=linspace(0,t1,20);
ti2=bsol(9);ti2=linspace(t1,ti2+t1,20);
time=[ti1,ti2];

figure,plot(time,kk(1,:),'r--',time,kk(2,:),'g--+',time,kk(3,:),'b-*',time(20),kk(1,20),'ko',time(20),kk(2,20),'ko',time(20),kk(3,20),'ko')
h = legend('joint 1','joint 2','joint 3',2);
xlabel('Time(s)')
ylabel('joint angle(rad)')
figure,plot(time,kk(4,:),'r--',time,kk(5,:),'g--+',time,kk(6,:),'b-*',time(20),kk(4,20),'ko',time(20),kk(5,20),'ko',time(20),kk(6,20),'ko')
h = legend('joint 1','joint 2','joint 3',2);
xlabel('Time(s)')
ylabel('joint vilocity(rad/s)')
figure,plot(time,kk(7,:),'r--',time,kk(8,:),'g--+',time,kk(9,:),'b-*',time(20),kk(7,20),'ko',time(20),kk(8,20),'ko',time(20),kk(9,20),'ko')
h = legend('joint 1','joint 2','joint 3',2);
xlabel('Time(s)')
ylabel('joint acceleration(rad/s^2)')

figure,plot(time,tt(1,:),'r--',time,tt(2,:),'g--+',time,tt(3,:),'b-*',time(20),tt(1,20),'ko',time(20),tt(2,20),'ko',time(20),tt(3,20),'ko')
h = legend('joint 1','joint 2','joint 3',2);
xlabel('Time(s)')
ylabel('joint tourque(N.m)')

figure,plot(e,trec)
xlabel('generation')
ylabel('consumed time for point to point motion(s)')
figure,plot(e,qrec)
xlabel('generation')
ylabel('total joint distance(rad)')
figure,plot(e,drec)
xlabel('generation')
ylabel('total cartesian trajectory length(m)')
figure,plot(e,torrec)
xlabel('generation')
ylabel('total excessive torque(N.m)')
bsol'

author: Sable

repo: mcbench-benchmarks


text:
function [roiInd coords] = roiGetAllIndices(vw)
% Get the indices of all open ROIs
% [roiInd coords] = roiGetAllCoords(vw)

if notDefined('vw'), vw = getCurView; end

roistodisplay = viewGet(vw, 'ROIs to Display');
ROIs = viewGet(vw, 'ROIs');

for ii = 1:length(roistodisplay)
    thisroi = roistodisplay(ii);

    if ii == 1,  coords = ROIs(thisroi).coords;
    else         coords = [coords ROIs(thisroi).coords]; end
end

% eliminate redundant voxels
if ~exist('coords', 'var') || isempty(coords)
    roiInd = []; coords = []; return;
end
coords = intersectCols(coords, coords);

% coords2index
[roiInd coords] = roiIndices(vw, coords);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [beta,betastar] = igls_orig(y, x, Vy, type)
% function [betastar] = igls(y, x, Vy)
%
% Variance Component Estimation using IGLS/RIGLS
%
% y = d + cx + epsilon  where  epsilon ~ N(0,sigma*Vy)
%
% d ~ N(0, sigma_d) and c ~ N(0, sigma_c)
%
% Calculate d, c, sigma, sigma_d and sigma_c using Maximum Likelihood
% methods (IGLS) and Restricted Maximum Likelihood methods (RIGLS).
%
% y - matrix T x subjects
% x - matrix T x subjects
% Vy - matrix T x T x subjects
% type = 'i' IGLS
% type = 'r' RIGLS
%
% By Martin Lindquist, April 2007
%
% Example:
%
% len = 200; sub =20;
% x = zeros(len,sub);
% x(11:20,:) = 2;
% x = x+ normrnd(0,0.1,len,sub);
% c = normrnd(0.5,0.1,sub,1);
% d = normrnd(3,0.2,sub,1);
% y=x;
% for i=1:sub, y(:,i) = d(i) + c(i).*x(:,i) + normrnd(0,0.5,len,1); end;
% Vy = zeros(len,len,sub);
% for i=1:sub, Vy(:,:,i) = eye(len,len).*0.5; end;
%
% [beta, betastar] = igls(y, x, Vy,'i')
%
c1= clock;
[T, sub] = size(y);             % Length of y vector and Number of subjects

one = zeros(T,1)+1;             % Vector of ones
null = zeros(T,1);              % Vector of zeros

epsilon = 0.001;
num_iter = 5;

len = sub*T;                                    % Total number of observations
z = reshape(y,len,1);                           % Concatenated data
D = [zeros(len,1)+1 reshape(x,len,1)];          % Design matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Step 1: Find the OLS solution
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

beta = pinv(D)*z;                           % Beta values
resid = z - D*beta;                         % Residuals

ystar = [];
if (type == 'i'),           % IGLS
     for i=1:sub,
         tmp = vech(resid(((i-1)*T+1):(i*T))*resid(((i-1)*T+1):(i*T))');                 % Find vech of estimated covariance
         ystar = [ystar; tmp];
     end;
elseif (type == 'r'),       % RIGLS
    for i=1:sub,
        Dtmp = D((((i-1)*T+1):(i*T)),:);
        rtmp = resid(((i-1)*T+1):(i*T));
        rig = rtmp*rtmp' + Dtmp*inv(Dtmp'*Dtmp)*Dtmp';
        tmp = vech(rig);                              % Find vech of estimated covariance
        ystar = [ystar; tmp];
    end;
end;

clear tmp rtmp rig Dtmp
G = Create_Design_Eq2(y, x, Vy);                      % Create design matrix for variance estimation
betastar = pinv(G)*ystar;                             % Estimate variance components

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Step 2: Iterate
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cnt = 0;
betastar_old = betastar+10;

iSigma = zeros(len,len);

while(cnt < num_iter | sum((betastar-betastar_old).^2)> epsilon),

    num = size(G,1)/sub;

    for i=1:sub,
        iSigma(((i-1)*T+1):(i*T),((i-1)*T+1):(i*T)) = ivech(G(((i-1)*num+1):(i*num),:)*betastar);
    end;

    beta = inv(D'*iSigma*D)*D'*iSigma*z;          % Beta values
    resid = z - D*beta;                           % Residuals

    ystar = [];

    if (type == 'i'),           % IGLS
         for i=1:sub,
             tmp = vech(resid(((i-1)*T+1):(i*T))*resid(((i-1)*T+1):(i*T))');                 % Find vech of estimated covariance
             ystar = [ystar; tmp];
         end;
    elseif (type == 'r'),       % RIGLS
         for i=1:sub,
            Dtmp = D((((i-1)*T+1):(i*T)),:);
            rtmp = resid(((i-1)*T+1):(i*T));
            rig = rtmp*rtmp' + Dtmp*inv(Dtmp'*iSigma((((i-1)*T+1):(i*T)),(((i-1)*T+1):(i*T)))*Dtmp)*Dtmp';
            tmp = vech(rig);                                                                 % Find vech of estimated covariance
            ystar = [ystar; tmp];
         end;
    end;

    clear tmp rtmp rig Dtmp

    betastar_old = betastar;
    betastar = pinv(G)*ystar;                                         % Estimate variance components

    cnt = cnt+1
end;

c2 = clock;
c2 - c1
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Subfunctions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [G] = Create_Design_Eq2(y, x, Vy)
% function [G] = Create_Design(y, x, Vy)
%
% Create Design matrix for estimation of variance componets
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T, sub] = size(x);
len = T*(T+1)/2;
ONE = zeros(len,1)+1;

G = [];
H = [];

for i=1:sub,

    XX = x(:,i)*x(:,i)';

    Gtmp = zeros(len,3);
    Gtmp(:,1) = ONE;
    Gtmp(:,2) = vech(XX);
    Gtmp(:,3) = vech(Vy(:,:,i));
    G = [G; Gtmp];

end;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function V = vech(Mat)
% function V = vech(Mat)
%
% Calculate vech for the matrix Mat
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

V = Mat(logical(tril(ones(size(Mat)))));

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Mat = ivech(V)
% function Mat = vech(V)
%
% Calculate the "inverse" of the vech function
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

len = length(V);
dim = -0.5 + sqrt(0.25 + 2*len);
Mat = zeros(dim,dim);
ind=1;

for i=1:dim
    for j=i:dim
        Mat(j,i)=V(ind);
        ind=ind+1;
    end
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ddbs = select_outcomes(ddbs)
%function ddbs = select_outcomes(ddbs)
%
%This function populates outcomes from Excel file into the ddbs structure
%
%APA, 12/24/2009
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

try
    %Get Excel File
    [fName,pName] = uigetfile('.xls','Select Excel file containing data');
    if isnumeric(fName)
        return;
    end
    xlFile = fullfile(pName,fName);
    [num,txt,raw] = xlsread(xlFile);

    columnNameRawC = {};
    indNaN = [];
    for i = 1:length(raw(1,:))
        if ~isnan(raw{1,i})
            columnNameRawC{end+1} = raw{1,i};
        else
            indNaN = [indNaN i];
        end
    end

    raw(:,indNaN) = [];

    %Get FileNames to match with DVH extraction
    [selIndFileName, OK] = listdlg('ListString',columnNameRawC,'PromptString','Select Column containing CERR fileNames');
    if OK && (length(selIndFileName) == 1)
        rowIndicesV = [];
        for j=1:size(raw,1)-1
            if ~isnan(raw{j+1,selIndFileName})
                excel_fileNameC{j} = raw{j+1,selIndFileName};
                rowIndicesV = [rowIndicesV j];
            end
        end
    else
        error('Please select ONE column and try again')
    end

    ddbs_fileNameC = {ddbs(:).fileName};
    ddbs_fileNameC = strtok(ddbs_fileNameC,'.');

    [jnk,indexMapV] = ismember(ddbs_fileNameC,excel_fileNameC);

    %replace spaces with underscores
    for i=1:size(raw,2)
        raw{1,i} = repSpaceHyp(raw{1,i});
    end

%     %removes columns with non-numeric data
%     indToRemove  = [];
%     for i=1:size(raw,2)
%         if ~isnumeric([raw{2:size(raw,1),i}]) || all(isnan([raw{2:size(raw,1),i}]))
%             indToRemove = [indToRemove i];
%         end
%     end
%     raw(:,indToRemove) = [];

    %get variables from the selection list
    [selMetricsIndV, OK] = listdlg('ListString',raw(1,:),'PromptString','Select Metrics');

    %get outcome from the selection list
    [selOutcomesIndV, OK] = listdlg('ListString',raw(1,:),'SelectionMode','single','PromptString','Select Outcome');

    %set variable fields
    for i=1:length(selMetricsIndV)
        for j = 1:length(indexMapV)
            if isnumeric(raw{indexMapV(j)+1,selMetricsIndV(i)})
                ddbs(j).(raw{1,selMetricsIndV(i)}) = raw{indexMapV(j)+1,selMetricsIndV(i)};
            else
                ddbs(j).(raw{1,selMetricsIndV(i)}) = NaN;
            end
        end
    end

    %set outcomes field
    for i=1:length(selOutcomesIndV)
        for j = 1:length(indexMapV)
            if isnumeric(raw{indexMapV(j)+1,selOutcomesIndV(i)})
                ddbs(j).outcome.(raw{1,selOutcomesIndV(i)}) = raw{indexMapV(j)+1,selOutcomesIndV(i)};
            else
                ddbs(j).outcome.(raw{1,selOutcomesIndV(i)}) = NaN;
            end
        end
    end

catch
    errordlg('There was an errror reading the DVH data. Please check the data format and try again.','DVH Import error','modal')
end

return;

% --------- supporting sub-functions
function str = repSpaceHyp(str)
indSpace = strfind(str,' ');
indDot = strfind(str,'.');
str(indDot) = [];
indOpenParan = strfind(str,'(');
indCloseParan = strfind(str,')');
indToReplace = [indSpace indOpenParan indCloseParan];
str(indToReplace) = '_';
indGreaterThan = strfind(str,'>');
indLessThan = strfind(str,'<');
str(indGreaterThan) = 'G';
str(indLessThan) = 'L';
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% NMFLAB for Signal Processing written by A. Cichocki and R. Zdunek
% in cooperation with other members of Laboratory for Advanced Brain Signal
% Processing, BSI, RIKEN, Saitama, JAPAN

function [x,grad_x,k] = nmf_proj_grad(b,A,x,tol,no_iter,beta_dec,loss_fun,beta)

% Projected Gradient Algorithm based on the Armijo rule
%
% [x, grad_x, k]=nmf_proj_grad(b,A,x,tol,no_iter,beta_dec,loss_fun,beta) finds
% such x that solves the equation Ax = b.
%
% INPUTS:
% b - data matrix of dimension [m by T]
% A - fixed matrix of dimension [m by r]
% x - initial guess
% tol - tolerance
% no_iter - maximum number of iterations
% beta_dec - multiplier of learning rate in descent gradient direction
% loss_fun - index for the loss function (1--4)
% beta - parameter "alpha" in the Amari alpha-divergence
%
% OUTPUTS:
% x - estimated matrix of dimension [r by T]
% grad_x - gradient of x
% k - number of performed iterations
%
% #########################################################################

if beta == 0
   loss_fun = 2;
end

sigma = 0.01; alpha = 1; kinc = 1;

for k=1:no_iter

 switch loss_fun

      case 1 % Frobenius

           res = b - A*x;
           grad_x = -A'*res;
           F = .5*norm(res,'fro')^2;
           lower_bound = 0;

      case 2 % KL

           bx = A*x;
           grad_x = A'*(1 - b./(bx + eps));
           F = sum(sum(b.*log(b./(bx + eps)) + bx - b));
           lower_bound = 0.1;

      case 3 % Dual KL

           bx = A*x;
           grad_x = A'*log(bx./b);
           F = sum(sum(bx.*log(bx./b) + b - bx));
           lower_bound = 1E-4;

      case 4 % Amari alpha divergence

           bx = A*x;
           grad_x = (1/beta)*A'*(1 - (b./bx).^beta);
           F = sum(sum(b.*((b./bx).^(beta - 1) - 1)/(beta^2 - beta) + (bx - b)/beta));
           lower_bound = 1E-2;

  end  % switch

  tol_grad = norm(grad_x(grad_x < lower_bound | x > lower_bound));
  if tol_grad < tol,
     break
  end

%  alpha = 1;
  s = 0;

% search step size
while (alpha > 1E-15) & (alpha < 1E8)

  s = s + 1;

  xn = max(x - alpha*grad_x,lower_bound);
  delta = xn-x;

   switch loss_fun

      case 1 % Frobenius

           Fn = .5*norm(b - A*xn,'fro')^2;

      case 2 % KL

           bx = A*xn;
           Fn = sum(sum(b.*log(b./(bx + eps)) + bx - b));

      case 3 % Dual KL

           bx = A*xn;
           Fn = sum(sum(bx.*log(bx./b) + b - bx));

      case 4 % Amari alpha divergence

           bx = A*xn;
           Fn = sum(sum(b.*((b./bx).^(beta - 1) - 1)/(beta^2 - beta) + (bx - b)/beta));

  end  % switch

  cond_alpha = ((Fn - F) <= sigma*grad_x(:)'*delta(:));

  if cond_alpha | (xn == x)
   % Increase
     if ~kinc & (s > 1), x = xn; break; end
     alpha = alpha/beta_dec;
     kinc = 1; xp = xn;
  else
   % Decrease
     if kinc & (s > 1), x = xp; alpha = alpha*beta_dec; break; end
     alpha = alpha*beta_dec;
     kinc = 0;
  end

end % while

end % for k
%

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function UNew = fluidNeumann3D(varargin);
% fluidNeumann3D: solve fluid registraion in 3D with Neumann
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[DU,F,mu,lambda,PixSize,M,N,P,RegularizerFactor,HX,HY,HZ] = parse_inputs(varargin{:});

% multiply F by adjoint of Navier-Lame equations
FNew = adjointNL(F/RegularizerFactor,mu,lambda,0,M,N,P);

% compute sine transform of new force field
FS = discreteCosineTransform(FNew,M,N,P);

% construct images of coordinates scaled by pi/(N or M or P)
[a,b,c] = ndgrid(pi*(0:(M-1))/(M-1),pi*(0:(N-1))/(N-1),pi*(0:(P-1))/(P-1));

% construct LHS factor
LHSfactor = mu.*(lambda+2*mu).*(2*cos(a) + 2*cos(b) + 2*cos(c) - 6).^2;

% if gamma is zero, set origin term to 1, as DC term does not matter
LHSfactor(1,1,1) = 1;

% solve for FFT of U
VS = cat(4,FS(:,:,:,1)./LHSfactor,FS(:,:,:,2)./LHSfactor,FS(:,:,:,3)./LHSfactor);

% perform inverse DST
V = discreteCosineTransform(VS,M,N,P);

% now perform Euler integration to construct new displacements
UNew = zeros(M,N,P,3);
UNew(:,:,:,1) = (1 - imfilter(V(:,:,:,1),HX,'replicate','same')).*V(:,:,:,1) - ...
    imfilter(V(:,:,:,2),HY,'replicate','same').*V(:,:,:,2) - ...
    imfilter(V(:,:,:,3),HZ,'replicate','same').*V(:,:,:,3);
UNew(:,:,:,2) = -imfilter(V(:,:,:,1),HY,'replicate','same').*V(:,:,:,1) + ...
    (1 - imfilter(V(:,:,:,2),HY,'replicate','same')).*V(:,:,:,2) - ...
    imfilter(V(:,:,:,3),HY,'replicate','same').*V(:,:,:,3);
UNew(:,:,:,3) = -imfilter(V(:,:,:,1),HZ,'replicate','same').*V(:,:,:,1) - ...
    imfilter(V(:,:,:,2),HZ,'replicate','same').*V(:,:,:,2) + ...
    (1 - imfilter(V(:,:,:,3),HZ,'replicate','same')).*V(:,:,:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FS = discreteCosineTransform(F,M,N,P);
% compute discrete cosine transform of 3-D vector field

% initialize resulting array
FS = F;

% first perform cosine transform down columns
len = 2*M-2; ind = 1:M;
for p=1:P
    for n=1:N
        s = fft(FS(:,n,p,:),len,1);
        FS(:,n,p,:) = real(s(ind,:,:,:));
    end
end
FS = sqrt(2/(M-1))*FS;

% next perform cosine transform across rows
len = 2*N-2; ind = 1:N;
for p=1:P
    for m=1:M
        s = fft(FS(m,:,p,:),len,2);
        FS(m,:,p,:) = real(s(:,ind,:,:));
    end
end
FS = sqrt(2/(N-1))*FS;

% finally perform cosine transform across pages
len = 2*P-2; ind = 1:P;
for n=1:N
    for m=1:M
        s = fft(FS(m,n,:,:),len,3);
        FS(m,n,:,:) = real(s(:,:,ind,:));
    end
end
FS = sqrt(2/(P-1))*FS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FNew = adjointNL(F,mu,lambda,gamma,M,N,P);
% multiply vector field F by adjoint Navier-Lame equations

% initialize FNew
FNew = zeros(M,N,P,3);

% construct filter that implements 3-D Laplacian
L = (lambda+2*mu)*cat(3,[0 0 0;0 1 0;0 0 0],[0 1 0;1 -6 1;0 1 0],[0 0 0;0 1 0;0 0 0]);

% we will need to use L to form two different filters
% L1 = -(lambda+2*mu)*L; L1(2,2,2) = gamma + L1(2,2,2);
% L2 = -mu*L; L2(2,2,2) = gamma + L2(2,2,2);

% construct grad div filters
GD11 = (lambda+mu)*cat(3,zeros(3,3),[0 1 0;0 -2 0;0 1 0],zeros(3,3));
GD22 = ipermute(GD11,[2 1 3]);
GD33 = ipermute(GD11,[3 2 1]);
GD23 = zeros(3,3,3);
GD23(2,1,1) = 1; GD23(2,3,3) = 1; GD23(2,1,3) = -1; GD23(2,3,1) = -1;
GD23 = GD23*(lambda+mu)/4;
GD12 = ipermute(GD23,[3 1 2]);
GD13 = ipermute(GD23,[2 3 1]);

% perform filtering
FNew(:,:,:,1) = imfilter(F(:,:,:,1),L-GD11,'replicate') + ...
    imfilter(F(:,:,:,2),-GD12,'replicate') + ...
    imfilter(F(:,:,:,3),-GD13,'replicate');
FNew(:,:,:,2) = imfilter(F(:,:,:,1),-GD12,'replicate') + ...
    imfilter(F(:,:,:,2),L-GD22,'replicate') + ...
    imfilter(F(:,:,:,3),-GD23,'replicate');
FNew(:,:,:,3) = imfilter(F(:,:,:,1),-GD13,'replicate') + ...
    imfilter(F(:,:,:,2),-GD23,'replicate') + ...
    imfilter(F(:,:,:,3),L-GD33,'replicate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DU,F,mu,lambda,PixSize,M,N,P,RegularizerFactor,HX,HY,HZ] = parse_inputs(varargin);

% get displacement field and check size
F = varargin{2};
PixSize = varargin{4}(1:3);
M = varargin{5};
N = varargin{6};
P = varargin{7};
mu = varargin{8};
lambda = varargin{9};
RegularizerFactor = varargin{10};
DU = varargin{11};
HX = varargin{12};
HY = varargin{13};
HZ = varargin{14};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
function variable = saveRemoteVariable(variable)
%"saveRemoteVariable"
%   Takes a remote variable and, if it was loaded, saves the data to the
%   file specified in remotePath, and clears the isLoaded field.  Returns
%   the modified variable.
%
% JRA 10/15/04
%
%Usage:
%   function variable = saveRemoteVariable(variable)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if isstruct(variable) & isfield(variable, 'remotePath') & isfield(variable, 'isLoaded') & isfield(variable, 'data') & isfield(variable, 'info') & strcmpi(variable.info, 'Remotely stored data.  Use getRemoteVariable to access.')'

    if variable.isLoaded
        switch upper(variable.storageType)
            case 'ZIP'
                fullTemp = fullfile(tempdir, variable.remotePath);
                CERRRemoteVariable = variable.data;

                %Save in ML6 style.
                saveOpt = getSaveInfo;
                if ~isempty(saveOpt);
                    save(fullTemp, 'CERRRemoteVariable', saveOpt);
                else
                    save(fullTemp, 'CERRRemoteVariable');
                end

                variable.data = [];
                variable.isLoaded = 0;
            case {'LOCAL', 'NETWORK'}
                error('Local and Network storage not implemented yet.');
        end
    else
        return;
    end
else
    error('Invalid remote data struct.')
end


text:
% @TENMAT
%
% Files
%   ctranspose - Complex conjugate transpose for tenmat.
%   disp       - Command window display of a matricized tensor (tenmat).
%   display    - Command window display of a tenmat.
%   double     - Convert tenmat to double array.
%   end        - Last index of indexing expression for tenmat.
%   minus      - Binary subtraction (-) for tenmat.
%   mtimes     - Multiplies two tenmat objects.
%   norm       - Frobenius norm of a tenmat.
%   plus       - Binary addition (+) for tenmat.
%   size       - Size of tenmat.
%   subsasgn   - Subscripted assignment for tenmat.
%   subsref    - Subscripted reference for tenmat.
%   tenmat     - Create a matricized tensor.
%   tsize      - Tensor size of tenmat.
%   uminus     - Unary minus (-) for tenmat.
%   uplus      - Unary plus (+) for tenmat.

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function value = ch_to_rot13 ( ch )

%*****************************************************************************80
%
%% CH_TO_ROT13 converts a character to its ROT13 equivalent.
%
%  Discussion:
%
%    Two applications of CH_TO_ROT13 to a character will return the original.!
%
%    As a further scrambling, digits are similarly rotated using
%    a "ROT5" scheme.
%
%  Example:
%
%    Input:  Output:
%
%    a       n
%    C       P
%    J       W
%    1       6
%    5       0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CH, the character to be converted.
%
%    Output, character VALUE, the ROT13 equivalent of the character.
%

%
%  [0:4] -> [5:9]
%
  if ( 48 <= ch && ch <= 52 )
    value = ch + 5;
%
%  [5:9] -> [0:4]
%
  elseif ( 53 <= ch && ch <= 57 )
    value = ch - 5;
%
%  [A:M] -> [N:Z]
%
  elseif ( 65 <= ch && ch <= 77 )
    value = ch + 13;
%
%  [N:Z] -> [A:M]
%
  elseif ( 78 <= ch && ch <= 90 )
    value = ch - 13;
%
%  [a:m] -> [n:z]
%
  elseif ( 97 <= ch && ch <= 109 )
    value = ch + 13;
%
%  [n:z] -> [a:m]
%
  elseif ( 110 <= ch && ch <= 122 )
    value = ch - 13;
  else
    value = ch;
  end

  value = char ( value );

  return
end


text:
function N = count_network_parameters(net)
% -------------------------------------------------------------------------
%   Description:
%       Count the total number of network parameters
%
%   Input:
%       - net   : network
%
%   Output:
%       - N     : #parameters
%
%   Citation:
%       Fast and Accurate Image Super-Resolution with Deep Laplacian Pyramid Networks
%       Wei-Sheng Lai, Jia-Bin Huang, Narendra Ahuja, and Ming-Hsuan Yang
%       arXiv, 2017
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

    N = 0;
    for i = 1:length(net.params)

        N = N + numel(net.params(i).value);

    end

end

repo: LapSRN


text:
classdef TensorVoigtConverterFactory < handle

    properties (Access = protected)
        input
    end

    methods (Access = protected)

        function itIs = isStrainTensor(obj)
            itIs = strcmp(obj.input.getFieldName(),'strain');
        end

        function itIs = isStressTensor(obj)
            itIs = strcmp(obj.input.getFieldName(),'stress');
        end

        function itIs = isComplianceTensor(obj)
            itIs = strcmp(obj.input.getFieldName,'compliance');
        end

        function itIs = isStiffnessTensor(obj)
            itIs = strcmp(obj.input.getFieldName,'stiffness');
        end

        function itIs = is3D(obj)
            itIs = strcmp(obj.input.getElasticityCase,'3D');
        end

        function itIs = isPlaneStress(obj)
            itIs = strcmp(obj.input.getElasticityCase,'planeStress');
        end

    end

    methods (Access = private, Static)

        function showError()
            error('Not admitted object to transform')
        end
    end

end


text:
function test_bug629

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_header ft_read_data ft_read_event read_mff_header read_mff_data

% note that read_mff_event does not exist, handling of the events is coded in ft_read_event

datadir = dccnpath('/home/common/matlab/fieldtrip/data/test/bug629');

cd(datadir);
dataset = 'pilot05_test 20110120 1433.mff';

hdr   = ft_read_header(dataset, 'headerformat', 'egi_mff_v1');
event = ft_read_event(dataset, 'eventformat', 'egi_mff_v1');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [feat, x, y, wid, hgt] = extract_sift(img, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~exist('c', 'var'))
  c = conf();
end

feature = 'sift';
p = c.feature_config.(feature);

if(size(img, 3)==3)
  img = rgb2gray(img);
end

[hgt wid] = size(img);
grid_spacing = p.grid_spacing;
x = cell(length(p.patch_sizes), 1);
y = cell(length(p.patch_sizes), 1);
feat = cell(length(p.patch_sizes), 1);

for i=1:length(p.patch_sizes)
  patch_size = p.patch_sizes(i);
  [x{i}, y{i}, gridX, gridY] = create_grid(hgt, wid, grid_spacing, patch_size);
  feat{i} = sp_find_sift_grid(img, gridX, gridY, patch_size, 0.8);
  feat{i} = sp_normalize_sift(feat{i});
end

x = cell2mat(x);
y = cell2mat(y);
feat = cell2mat(feat);

author: adikhosla

repo: feature-extraction


text:
function cvx_erase( p )

global cvx___
pstr = cvx___.problems( p );
% n_variables, n_equalities, n_cones
cp = pstr.checkpoint + 1;
nf = cp(1);
if nf <= 2,
    cvx___.classes     = int8(3);
    cvx___.cones       = zeros(0,1);
    cvx___.exponential = zeros(0,1);
elseif length( cvx___.classes ) >= nf,
    cvx___.classes( nf : end, : ) = [];
    if ~isempty( cvx___.exponential ),
        cvx___.exponential( nf : end, : ) = [];
        cvx___.logarithm( nf : end, : ) = [];
        if ~any( cvx___.exponential ),
            cvx___.exponential = zeros(0,1);
            cvx___.logarithm = zeros(0,1);
        end
    end
end
ne = cp(2);
if ne <= 1,
    cvx___.equalities = {};
    cvx___.inequality = false(0,1);
    cvx___.n_equality = 0;
elseif length( cvx___.equalities ) >= ne,
    cvx___.n_equality = cvx___.n_equality - ...
        sum( cellfun( @(x)size(x,2), cvx___.equalities( ne : end ) ) );
    cvx___.equalities( ne : end ) = [];
    cvx___.inequality( ne : end ) = [];
end
nc = cp(3);
if nc <= 1,
    cvx___.cones = [];
else
    cvx___.cones( nc : end ) = [];
end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% test of hufmann coding by concatening several token

% probability of having 0
t = .12;
n = 4096*2;
x = (rand(n,1)>t)+1;

% entropy lower bound
p = [t 1-t];
e =  -sum(p.*log2(p));

% create a new vector by lifting
q = 3;
n1 = ceil(n/q)*q;
x1 = x;
x1(end+1:n1) = 1;
x1 = reshape(x1,[q n1/q]);
[Y,X] = meshgrid(1:n1/q,0:q-1);
x1 = sum( (x1-1) .* (2.^X), 1 )' + 1;

% generate probability table
P = p(:); p = p(:);
for i=1:q-1
    Pold = P;
    P = [];
    for i=1:length(p)
        P = [P; Pold*p(i)];
    end
end

% compute the tree
T = compute_hufftree(P);
% do the coding
y = perform_huffcoding(x1,T,+1);
% average number of bits
e1 = length(y)/length(x);

disp(['Entropy=' num2str(e) ', Huffman(block size ' num2str(q) ')=' num2str(e1)]);

author: gpeyre

repo: matlab-toolboxes


text:
function  plotmap_all(pathres, sep, offset, zvec, scalecont)
% cd ~project/data/qdots/S53
% sep=[10 20 50];
% offset=[100 1000];
cd (pathres)
% sufix = {'avg1', 'avg2', 'avg10', 'avg50' 'avg0'};
sufix = {'avg1', 'avg10', 'avg50' 'avg0'};
for is=1:size(sep,2)
    limvec1=[1:15];
    limvec2=[1:15];
    if sep(is)>=0.5
        limvec1=[1:13];
        limvec2=[1:13];
    end

    for io=1:size(offset,2)
        for ii=1:size(sufix,2)
            namebase=['S44_sep_' num2str(100*sep(is)) 'offset_' num2str(offset(io))];
            nameload=[namebase '/' namebase '_DdviMap_' sufix{ii} '.mat'];
            load (nameload)
%             plotmap(res.X1,res.X2,log(res.dh),[namebase '/dh_' sufix{ii}],zvec,limvec1,limvec2,scalecont,'on')
            minval=min(res.mhd(:));
            plotmap(res.X1,res.X2,log(res.mhd-minval+0.01),[namebase '/dh_' sufix{ii}],1,limvec1,limvec2,scalecont,'on');
        end
    end
end


text:
%%  Introduction to *wavelet_factory_2d*
% *wavelet_factory_2d* is computing operators(except the modulus) and
% filters required to compute the next layer of a scattering network.
%
%% Usage
% [Wop, filters] = wavelet_factory_2d(size_in, filt_opt, scat_opt), documentation is given in
% <matlab:doc('wavelet_factory_2d') wavelet_factory_2d>
%
%% Description
% Given a size image, some filters options and scattering options, this
% function comput the linear operators necessar to compute the next
% coefficients of scattering.

x = mandrill;

% Create $ U[\empty]x $
[Wop, filters] = wavelet_factory_2d(size(x));

% Then one can apply Wop as in *scat*. Please reference to its
% documentation.

%% Options
% filt_opt has the same fields as in *morlet_filter_bank_2d*.
%
% scat_opt has the same fields as in *wavelet_layer_2d*.
%
% See their documentation for more details.

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function rts_retrieveOrbitAndStore(~, ~,tcpipClient)
%rts_retrieveOrbitAndStore Summary of this function goes here
%   Detailed explanation goes here
    rts_retrieveDataAndStore(tcpipClient, 19, 'd');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot