text:
clear all;
cd('./MLP');
disp('Testing training a multilayer perceptron.')
Main_MNIST_MLP_RMSPROP();
cd ..

clear all;

cd('./ReinforcementLearning');
disp('Testing training a Q-network.')
Main_Cart_Pole_Q_Network
cd ..

clear all;

cd('./ReinforcementLearning');
disp('Testing training a policy network.')
Main_Cart_Pole_Policy_Network
cd ..

clear all;
cd('./RNN');
disp('Testing training a Quasi-RNN.')
Main_Char_RNN();
cd ..

clear all;
cd('./CNN');
disp('Testing using a pretrained ImageNet convolutional neural network model.')
Main_CNN_ImageNet_minimal();
cd ..

clear all;
cd('./CNN');
disp('Testing training a new convolutional neural network using SGD.')
disp('An Nvidia GPU is required by default.')
disp('If you do not have it please set use_gpu=0')
disp('Neural Network Toolbox is required by default.')
disp('If you do not have it please set opts.use_nntoolbox=0')
Main_CIFAR_CNN_SGD();
cd ..

clear all;
disp('Second-order SGD test 1.')
cd('./SGD2');
disp('Testing training a new convolutional neural network using SGD2.')
disp('An Nvidia GPU is required by default.')
disp('If you do not have it please set use_gpu=0')
disp('Neural Network Toolbox is required by default.')
disp('If you do not have it please set opts.use_nntoolbox=0')
Main_CIFAR_CNN_SGD2();
cd ..

clear all;
cd('./SGD2');
disp('Second-order SGD test 2.')
disp('SGD2 initialization tests.')
MNIST_Deep_MLP_Init_Comparison();
cd ..

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text: any(cellfun(@exist,{'fast_median', 'fast_median_ip', 'nth_element', 'nth_element_ip'})==3)

repo: REST


text:
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% This code is also for our paper in Energy&Fuels
% I'll have to work on some details
% It works fine for now
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
%% define the geometry
Nx = 50; % number of cells in x direction
Ny = 30; % number of cells in y direction
W = 300; % [m] length of the domain in x direction
H = 30; % [m] length of the domain in y direction
m = createMesh1D(Nx, W); % creates a 1D mesh
%% define the initial and boundary conditions for concentrations
c_inj_ion = 1.0; % injection concentration
c_init_ion = 0.0;
BCc = createBC(m);
BCc.left.a = 0.0; BCc.left.b = 1.0; BCc.left.c = c_inj_ion; % left boundary
c_ion_old = createCellVariable(m, c_init_ion, BCc);
c_ion = c_ion_old;
[M_bc_c, RHS_bc_c]=boundaryCondition(BCc);

%% define the diffusion adsorption domains for each cell
% partameters
rho_s = 2700; % kg/m3
a_s = 2000; % m2/kg
a = 1e9; % m^2/m^3
k_lang = 1e-7; % Langmuir adsorption coefficient
betta = 1.0; % Langmuir adsorption coefficient 2
% Define the domain and create a mesh structure
L = 1e-6; % [m]  % domain length
N = 15; % number of cells
m_diff = createMesh1D(N, L);
% Create the boundary condition structure
BC = cell(Nx, 1);
for i = 1:Nx
    BC{i} = createBC(m_diff); % all Neumann boundary condition structure
    BC{i}.left.a = 0; BC{i}.left.b=1; BC{i}.left.c=0; % left boundary
    BC{i}.right.a = 0; BC{i}.right.b=1; BC{i}.right.c=0; % right boundary
end
D_val = 1e-13; % m^2/s effective diffusivity
D = createCellVariable(m_diff, D_val);
D_face = harmonicMean(D);
Mdiff = diffusionTerm(D_face);
alfa = createCellVariable(m_diff, 1);
% initial condition
c_init = 0;
c_old = cell(Nx);
for i = 1:Nx
    c_old{i} = createCellVariable(m_diff, c_init, BC{i}); % initial values
end
c = c_old; % assign the old value of the cells to the current values
%% define the physical parametrs
% define the oil-wet rel-perm
krw0 = 0.5;
kro0 = 0.76;
nw = 2.4;
no = 2.0;
sor=0.15;
swc=0.09;
sws=@(sw)((sw>swc).*(sw<1-sor).*(sw-swc)/(1-sor-swc)+(sw>=1-sor).*ones(size(sw)));
kro=@(sw)((sw>=swc).*kro0.*(1-sws(sw)).^no+(sw<swc).*(1+(kro0-1)/swc*sw));
krw=@(sw)((sw<=1-sor).*krw0.*sws(sw).^nw+(sw>1-sor).*(-(1-krw0)/sor.*(1.0-sw)+1.0));
dkrwdsw=@(sw)((sw<=1-sor).*nw.*krw0.*(1/(1-sor-swc)).*sws(sw).^(nw-1)+(sw>1-sor)*((1-krw0)/sor));
dkrodsw=@(sw)((sw>=swc).*(-kro0*no*(1-sws(sw)).^(no-1))/(-swc-sor+1)+(sw<swc).*((kro0-1)/swc));

% define water wet relperms
krw0_ww = 0.3;
kro0_ww = 0.96;
nw_ww = 2.4;
no_ww = 2.0;
sor_ww =0.05;
swc_ww =0.09;
sws_ww=@(sw)((sw>swc_ww).*(sw<1-sor_ww).*(sw-swc_ww)/(1-sor_ww-swc_ww)+(sw>=1-sor_ww).*ones(size(sw)));
kro_ww=@(sw)((sw>=swc_ww).*kro0_ww.*(1-sws_ww(sw)).^no_ww+(sw<swc_ww).*(1+(kro0_ww-1)/swc_ww*sw));
krw_ww=@(sw)((sw<=1-sor_ww).*krw0_ww.*sws_ww(sw).^nw_ww+(sw>1-sor_ww).*(-(1-krw0_ww)/sor_ww.*(1.0-sw)+1.0));
dkrwdsw_ww=@(sw)((sw<=1-sor_ww).*nw.*krw0_ww.*(1/(1-sor-swc_ww)).*sws_ww(sw).^(nw_ww-1)+(sw>1-sor_ww)*((1-krw0_ww)/sor_ww));
dkrodsw_ww=@(sw)((sw>=swc_ww).*(-kro0_ww*no*(1-sws_ww(sw)).^(no_ww-1))/(-swc_ww-sor_ww+1)+(sw<swc_ww).*((kro0_ww-1)/swc_ww));

%% wettability modifier
SF=createCellVariable(m, 0.0); % 1 is water wet, 0 is oil wet
SF_face = arithmeticMean(SF);

p0 = 100e5; % [bar] pressure
pin = 150e5; % [bar] injection pressure at the left boundary
u_in= 1.0/(24*3600); % [m/s] equal to 1 m/day
sw0 = swc+0.1; % initial water saturation
sw_in = 1;
mu_oil = 2e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 2e-12; % [m^2] average reservoir permeability
phi0 = 0.2; % average porosity
clx=1.2;
cly=0.2;
V_dp=0.7; % Dykstra-Parsons coef.
perm_val= k0; %field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
%% Define the boundaries
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% left boundary pressure gradient
BCp.left.a(:)=(krw(sw_in)*lw.xvalue(1,:)+kro(sw_in)*lo.xvalue(1,:)); BCp.left.b(:)=0; BCp.left.c(:)=-u_in;
% change the right boandary to constant pressure (Dirichlet)
% BCp.left.a(:)=0; BCp.left.b(:)=1; BCp.left.c(:)=pin;
BCp.right.a(:)=0; BCp.right.b(:)=1; BCp.right.c(:)=p0;
% change the left boundary to constant saturation (Dirichlet)
BCs.left.a(:)=0; BCs.left.b(:)=1; BCs.left.c(:)=1;
%% define the time step and solver properties
% dt = 1000; % [s] time step
dt=(W/Nx)/u_in/10; % [s]
t_end = 400*dt; % [s] final time
eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-5; % saturation accuracy
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
oil_init=domainInt(1-sw_old);

p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
p = p_old;
uw = -gradientTerm(p_old); % an estimation of the water velocity
%% start the main loop
% generate intial pressure profile (necessary to initialize the fully
% implicit solver)
rec_fact=0;
t_day=0;
dp_alwd= 100.0; % Pa
dsw_alwd= 0.05;
t = 0;
% fprintf(1, 'progress (%%):  ');
while (t<t_end)
    error_p = 1e5;
    error_sw = 1e5;
    % Implicit loop
    loop_count=0;
    while ((error_p>eps_p) || (error_sw>eps_sw))
        loop_count=loop_count+1;
        if loop_count>10
            break
        end
        % calculate parameters
        pgrad = gradientTerm(p);
        sw_face = upwindMean(sw, -pgrad); % average value of water saturation
        labdao = ((1-SF_face).*funceval(kro, sw_face)+SF_face.*funceval(kro_ww, sw_face)).*lo;
        labdaw = ((1-SF_face).*funceval(krw, sw_face)+SF_face.*funceval(krw_ww, sw_face)).*lw;
        dlabdaodsw = ((1-SF_face).*funceval(dkrodsw, sw_face)+SF_face.*funceval(dkrodsw_ww, sw_face)).*lo;
        dlabdawdsw = ((1-SF_face).*funceval(dkrwdsw, sw_face)+SF_face.*funceval(dkrwdsw_ww, sw_face)).*lw;
        labda = labdao+labdaw;
        dlabdadsw = dlabdaodsw+dlabdawdsw;
        % compute [Jacobian] matrices
        Mdiffp1 = diffusionTerm(-labda);
        Mdiffp2 = diffusionTerm(-labdaw);
        Mconvsw1 = convectionUpwindTerm(-dlabdadsw.*pgrad);
        Mconvsw2 = convectionUpwindTerm(-dlabdawdsw.*pgrad);
        [Mtranssw2, RHStrans2] = transientTerm(sw_old, dt, phi);
        % Compute RHS values
        RHS1 = divergenceTerm(-dlabdadsw.*sw_face.*pgrad);
        RHS2 = divergenceTerm(-dlabdawdsw.*sw_face.*pgrad);
        % include boundary conditions
        [Mbcp, RHSbcp] = boundaryCondition(BCp);
        [Mbcsw, RHSbcsw] = boundaryCondition(BCs);
        % Couple the equations; BC goes into the block on the main diagonal
        M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw];
        RHS = [RHS1+RHSbcp; RHS2+RHStrans2+RHSbcsw];
        % solve the linear system of equations
        x = M\RHS;
        % x = agmg(M, RHS, [], 1e-10, 500, [], [p.value(:); sw.value(:)]);
        % separate the variables from the solution
        p_new = reshapeCell(m,full(x(1:(Nx+2))));
        sw_new = reshapeCell(m,full(x((Nx+2)+1:end)));
        % calculate error values
        error_p = max(abs((p_new(:)-p.value(:))./p_new(:)));
        error_sw = max(abs(sw_new(:)-sw.value(:)));
        % assign new values of p and sw
        p.value = p_new;
        sw.value = sw_new;
    end
    if loop_count>10
      p=p_old;
      sw=sw_old;
      dt=dt/5;
      continue
    end
    % solve the ions flow in the domain
    uw = -labdaw.*pgrad; % water velocity
    Mconv = convectionUpwindTerm(uw);
    dswdt=(sw-sw_old)/dt;
    Msc=linearSourceTerm(phi.*dswdt);
    for j = 1:3
        [M_trans, RHS_trans] = transientTerm(c_ion_old, dt, sw.*phi+a_s*rho_s*k_lang.*(1-phi)./(1+betta*c_ion));
        M = M_trans+Mconv+M_bc_c+Msc;
        RHS = RHS_trans+RHS_bc_c;
        c_ion = solvePDE(m,M, RHS);
    end
    c_ion_old = c_ion;

    % solve the ions diffusion in the water films
    for i = 1:Nx
        BC{i}.left.a = 0; BC{i}.left.b=1; BC{i}.left.c=c_ion.value(i+1); % left boundary
        BC{i}.right.a = 0; BC{i}.right.b=1; BC{i}.right.c=c_ion.value(i+1); % right boundary
        [Mbc, RHSbc] = boundaryCondition(BC{i});
        for j=1:3
            [M_trans, RHS_trans] = transientTerm(c_old{i}, dt, 1.0+a*k_lang./(1+betta*c{i}));
            M = M_trans-Mdiff+Mbc;
            RHS = RHS_trans+RHSbc;
            c{i} = solvePDE(m_diff,M, RHS);
        end

        ind_80 = find(c{i}.value(2:end-1)<0.8*c_inj_ion,1, 'first');
        if isempty(ind_80)
            SF.value(i+1) = 1.0;
        else
            SF.value(i+1) = 2*(ind_80-1)/(N-2);
        end
        c_old{i} = c{i};
    end

    SF_face = arithmeticMean(SF);

    dsw=max(abs(sw_new(:)-sw_old.value(:))./sw_new(:));
    t=t+dt;
    fprintf(1,'\b\b\b\b\b\b\b\b\b\b\b\b\b\bProgress: %d %%',floor(t/t_end*100));
    dt=min([dt*(dsw_alwd/dsw), 1.05*dt, t_end-t]);
    p_old = p;
    sw_old = sw;
    rec_fact=[rec_fact (oil_init-domainInt(1-sw))/oil_init];
    t_day=[t_day t];
    figure(1);visualizeCells(sw); drawnow;
    figure(2); visualizeCells(c_ion);

%     figure(2);plot(t_day/3600/24, rec_fact)
%     xlabel('time [day]');
%     ylabel('recovery factor');
%     title([num2str(t/3600/24) ' day']); drawnow;
end
fprintf('\n')

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


function [rec,prec,ap] = VOCevaldet_imglistin(VOCopts,id,cls,imglist,draw, ovthresh)
tic;
if(~exist('ovthresh', 'var'))
	ovthresh=0.5;
end
% load test set

cp=sprintf(VOCopts.annocachepath,'newval');
if exist(cp,'file')
    fprintf('%s: pr: loading ground truth\n',cls);
    load(cp,'gtids','recs');
else
    [gtids,t]=textread(sprintf(VOCopts.imgsetpath,VOCopts.testset),'%s %d');
	gtids=gtids(ismember(gtids, imglist));
    for i=1:length(gtids)
        % display progress
        if toc>1
            fprintf('%s: pr: load: %d/%d\n',cls,i,length(gtids));
            drawnow;
            tic;
        end

        % read annotation
        recs(i)=PASreadrecord(sprintf(VOCopts.annopath,gtids{i}));
    end
    save(cp,'gtids','recs');
end

fprintf('%s: pr: evaluating detections\n',cls);

% hash image ids
%hash=VOChash_init(gtids);
hash = containers.Map;
for i = 1:length(gtids)
  hash(gtids{i}) = i;
end
% extract ground truth objects

npos=0;
gt(length(gtids))=struct('BB',[],'diff',[],'det',[]);
for i=1:length(gtids)
    % extract objects of class
    clsinds=strmatch(cls,{recs(i).objects(:).class},'exact');
    gt(i).BB=cat(1,recs(i).objects(clsinds).bbox)';
    gt(i).diff=[recs(i).objects(clsinds).difficult];
    gt(i).det=false(length(clsinds),1);
    npos=npos+sum(~gt(i).diff);
end

% load results
[ids,confidence,b1,b2,b3,b4]=textread(sprintf(VOCopts.detrespath,id,cls),'%s %f %f %f %f %f');
BB=[b1 b2 b3 b4]';

% sort detections by decreasing confidence
[sc,si]=sort(-confidence);
ids=ids(si);
BB=BB(:,si);

% assign detections to ground truth objects
nd=length(confidence);
tp=zeros(nd,1);
fp=zeros(nd,1);
tic;
for d=1:nd
    % display progress
    if toc>1
        fprintf('%s: pr: compute: %d/%d\n',cls,d,nd);
        drawnow;
        tic;
    end

    % find ground truth image
    %i=VOChash_lookup(hash,ids{d});
    i = hash(ids{d});
	if isempty(i)
        error('unrecognized image "%s"',ids{d});
    elseif length(i)>1
        error('multiple image "%s"',ids{d});
    end

    % assign detection to ground truth object if any
    bb=BB(:,d);
    ovmax=-inf;
    for j=1:size(gt(i).BB,2)
        bbgt=gt(i).BB(:,j);
        bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
        iw=bi(3)-bi(1)+1;
        ih=bi(4)-bi(2)+1;
        if iw>0 & ih>0
            % compute overlap as area of intersection / area of union
            ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
               (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
               iw*ih;
            ov=iw*ih/ua;
            if ov>ovmax
                ovmax=ov;
                jmax=j;
            end
        end
    end
    % assign detection as true positive/don't care/false positive
    if ovmax>=ovthresh%VOCopts.minoverlap
        if ~gt(i).diff(jmax)
            if ~gt(i).det(jmax)
                tp(d)=1;            % true positive
		gt(i).det(jmax)=true;
            else
                fp(d)=1;            % false positive (multiple detection)
            end
        end
    else
        fp(d)=1;                    % false positive
    end
end

% compute precision/recall
fp=cumsum(fp);
tp=cumsum(tp);
rec=tp/npos;
prec=tp./(fp+tp);

ap=VOCap(rec,prec);

if draw
    % plot precision/recall
    plot(rec,prec,'-');
    grid;
    xlabel 'recall'
    ylabel 'precision'
    title(sprintf('class: %s, subset: %s, AP = %.3f',cls,VOCopts.testset,ap));
end


text:
% function w = phase_cwt(Wx, dWx, opt)
%
% Calculate the phase transform at each (scale,time) pair:
%   w(a,b) = Im( (1/2pi) * d/db [ Wx(a,b) ] / Wx(a,b) )
% Uses direct differentiation by calculating dWx/db in frequency
% domain (the secondary output of cwt_fw, see help cwt_fw)
%
% This is the analytic implementation of Eq. (7) of [1].
%
% 1. G. Thakur, E. Brevdo, N.-S. Fučkar, and H.-T. Wu,
% "The Synchrosqueezing algorithm for time-varying spectral analysis: robustness
%  properties and new paleoclimate applications," Signal Processing, 93:1079-1094, 2013.
%
% % 2. I. Daubechies, J. Lu, H.T. Wu, "Synchrosqueezed Wavelet Transforms: an
% empricial mode decomposition-like tool", Applied and Computational Harmonic Analysis
% 30(2):243-261, 2011.
%
% Input:
%  Wx: wavelet transform of x (see help cwt_fw)
%  dWx: samples of time derivative of wavelet transform of x (see help cwt_fw)
%  opt: options struct,
%    opt.gamma: wavelet threshold (default: sqrt(machine epsilon))
%
% Output:
%  w: phase transform, size(w) = size(Wx)
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo, Gaurav Thakur
%---------------------------------------------------------------------------------
function w = phase_cwt(Wx, dWx, opt)
    if nargin<3, opt = struct(); end
    % epsilon, gamma from [1]
    if ~isfield(opt, 'gamma'); opt.gamma = sqrt(eps); end

    % Calculate phase transform for each ai, normalize by (2*pi)
	if strcmpi(opt.dtype,'phase')
		u = unwrap(angle(Wx)).';
		w = [diff(u);u(end,:)-u(1,:)].'/(2*pi);
	else
		w = abs(imag(dWx./Wx/(2*pi)));
	end
    w(abs(Wx)<opt.gamma) = Inf;
end

author: ebrevdo

repo: synchrosqueezing


text:
function [model,delta] = getPredictors(stimList, HRF, varargin)
% Build predictors and delta functions, given a condition function or delta
% function and either a convolution matrix or vector.
%
% :Usage:
% ::
%
%     [model,delta] = getPredictors(stimList, HRF, varargin)
%
% IMPORTANT: YOU MUST ADD THE INTERCEPT YOURSELF!
%
% :Inputs:
%
%   **stimList:**
%        condition function OR delta function (1/0 indicator)
%
%   **HRF:**
%        1) hemodynamic response function
%        2) Basis set (columns)
%        3) or convolution matrix (columns
%           are HRF), defined as:
%           HRF = tril(toeplitz(hrf));
%
%           multiple column vectors for HRF are treated as basis functions!
%
%   **varargin** for downsampling:
%        'dsrate': takes every nth element of the design matrix
%
%        'dslen': the target number (length) you want to downsample to
%
% :Other Optional Inputs:
%
%   **'force_delta':**
%        getPredictors tries to determine if the input stimList
%        is a condition function with integers or a delta function with
%        indicators, but this can fail in some cases.  Use this to force
%        it to treat as a delta function.
%
%   1. a col. vector of stimulus conditions OR a delta function matrix
%
%   2. an HRF vector sampled at the frequency of the stimulus vector, OR
%        a convolution matrix H (empty for default)
%
%   3. Optional: downsampling factor for final design (i.e., TR)
%
%   4. Optional: parametric modulator keyword and modulator values
%
%      'parametric_singleregressor' : Parametrically modulate onsets by
%               modulator values, using single regressor with modulated amplitude
%               Enter a cell array with modulator values for each event
%               type, with a column vector (empty cell for no modulation)
%
%       'parametric_standard' : Parametrically modulate onsets by
%               modulator values, using two regressors per event type - One
%               to model the average response, and one for the
%               mean-centered modulator values
%
% :Outputs:
%
%   1. a n x 2 matrix of regressors (cols) for each condition
%   2. a n x k delta matrix with onsets
%
% :Example: TR = 2, 16 samples per second in hi-res delta dhr
% ::
%
%    X = getPredictors(dhr, hrf, 'dsrate', res*TR);
%    X = getPredictors(dhr, hrf, 'dslen', len/(res*TR));
%
% stimList can be condition function e.g., [1 3 2 4 3 2 1]' or
% delta matrix (n x k), n samples and k conditions, e.g., [1 0 0 0 1 0 1]'
%
% :Resampling: the default N in matlab resample has built-in antialiasing,
% but may not be good for fmri designs!  The appropriate downsampling
% is expected to be res*TR (res is units of samples/s), but we use 0
% because the model will depend on the analysis method used, and this is
% the most veridical approach.  With N = 0, every ith sample is used, where
% i is the downsampling factor you input.  Popular choices are 16*TR (for
% onsets2delta.m), using the SPM default res of 16.
% Delta is NOT resampled.
%
% :Example: TR = 2, 16 samples per second in hi-res delta dhr
% ::
%
%    [tmp,d] = downsample_delta(dhr,16*2); X=getPredictors(d,hrf);
%
% ..
%    Notes:
%    Tor Wager, last modified 2/22/04 to center predictors
%    modified to optionally take H convolution matrix as input
%    in place of HRF vector.  See convmtx.m (Buracas,mseq toolbox)
%
%    Modified 10/2011 by Tor to add parametric modulators
%
%    3/15/12: Tor updated to consider any stimlist with > 30 conditions to be
%    a delta matrix, so that one can input continuous neural response
%    functions and treat them as delta matrices for direct convolution.
%
%    9/3/12: Wani updated this function for two reasons. 1) to fix an error
%    when one condition has two or more parametric modulators. 2) to fix an
%    error due to RT with decimal points. For 2), make two options for downsampling
%   'dsrate' and 'dslen'.
%
%    1/2020: Tor modified to increase tolerance/flexibilty for fractional
%    TRs, add better error reporting for length mismatch
% ..

model = [];
issquare = all(size(HRF) == size(HRF,1));   % if square mtx, assume convolution mtx

% make sure HRF is column vector, if single vector
if size(HRF,1) == 1, HRF = HRF'; end

% parametric modulators
doing_parametric_standard = any(strcmp(varargin, 'parametric_standard'));
doing_parametric_singleregressor = any(strcmp(varargin, 'parametric_singleregressor'));

if doing_parametric_standard || doing_parametric_singleregressor
    wh = find(strcmp(varargin, 'parametric_standard') | strcmp(varargin, 'parametric_singleregressor'));
    pm_vals = varargin{wh(1) + 1};

    if ~iscell(pm_vals)
        error('pm_vals must be cell array.')
    end
end

use_downsample_rate = any(strcmp(varargin, 'dsrate'));
use_downsample_length = any(strcmp(varargin, 'dslen'));
if use_downsample_rate, wh = find(strcmp(varargin, 'dsrate')); rate_downsample = varargin{wh(1) +1}; end
if use_downsample_length, wh = find(strcmp(varargin, 'dslen')); length_downsample = varargin{wh(1) +1}; end

if isdeltamtx(stimList, varargin{:})  % delta matrix
    % -------------------------------------------------------------------------------------------------
    %    * If delta function
    % -------------------------------------------------------------------------------------------------

    delta = stimList;

    if ~issquare
        for i = 1:size(delta,2)

            for j = 1:size(HRF,2)

                if doing_parametric_standard || doing_parametric_singleregressor
                    if length(pm_vals) < i, error('PM values cell array is too short!'); end

                    model = [model pmconv(double(delta(:,i)), HRF(:,j), pm_vals{i}, doing_parametric_standard, doing_parametric_singleregressor)];
                else
                    model(:,end+1) = conv(double(delta(:,i)), HRF(:,j)); % Changed for Matlab 7.9+ compatibility - Thanks, Liane
                end
            end
        end
    end

else
    % -------------------------------------------------------------------------------------------------
    %    * If condition function
    % -------------------------------------------------------------------------------------------------

    for i = 1:max(stimList(:,1)) % condition function

        %delta(:,i) = (stimList == i);
        delta(:,i) = cast((stimList == i), 'single'); % Changed for Matlab 7.9+ compatibility -- thanks, Liane Montana and Bruce McCandliss

        if doing_parametric_standard || doing_parametric_singleregressor
            if length(pm_vals) < i, error('PM values cell array is too short!'); end

            model = [model pmconv(double(delta(:,i)), HRF(:,j), pm_vals{i}, doing_parametric_standard, doing_parametric_singleregressor)];
        end

        if ~issquare
            for j = 1:size(HRF,2)
                model(:,end+1) = conv(double(delta(:,i)), HRF(:,j)); % Changed for Matlab 7.9+ compatibility - Thanks, Liane
            end
        end

    end
end

% -------------------------------------------------------------------------------------------------
%    * If conv. matrix
% -------------------------------------------------------------------------------------------------

if issquare % convolution matrix
    model = HRF * delta;
end

model = model(1:size(stimList,1),:);      	% eliminate extra values

% downsample, if necessary
if use_downsample_rate || use_downsample_length % Wani modified this
% if ~isempty(varargin)

    % dsrate = varargin{1}; % Wani
    [n, k] = size(model);
    % nt = n ./ dsrate; % Wani
    if use_downsample_rate, nt = n ./ rate_downsample; end % Wani modified this line
    if use_downsample_length, nt = length_downsample; rate_downsample = n/length_downsample; end % Wani added this line
    t = (1:n)';             % minor change to save time, tor: 10/12/10

    if abs(nt-round(nt)) > .1 % Tor increased tolerance to 0.1 units (usually TR/16)
        fprintf('X samples: %3.0f, Xout (image) samples: %3.2f, DSrate: %3.2f\n', n, nt, rate_downsample);
        error('Length of stimList is not evenly divisible by downsampling factor.');
    end

    modeli = zeros(round(nt) , k);

    for i = 1:size(model, 2)

        xi = 1:rate_downsample:n;       % downsample rate

        modeli(:, i) = interp1(t, model(:, i), xi, 'linear', 'extrap'); % tor: 10/12/10 : minor change: add extrap
    end

    model = modeli;
    %model = model(1:varargin{1}:end,:); % equivalent to resample(X,1,varargin{1},0)
end

% do not do this if you want to do nonlinear saturation!
% not necessary, i think.
%model = model - repmat(mean(model),size(model,1),1);    % center predictors

end

function isdelta = isdeltamtx(stimList, varargin)

isdelta = 0;

if ~isempty(varargin) && any(strcmp(varargin, 'force_delta'))
    isdelta = 1;
    return
end

if islogical(stimList) || all( sum(stimList == 0 | stimList == 1) == size(stimList, 1) )  % all zeros or ones
    isdelta = 1;

elseif min(size(stimList)) > 1 % multiple columns
    isdelta = 1;

elseif length(unique(stimList)) > 30 % custom neural response function?

end

end % subfunction

function [model, delta] = pmconv(delta, HRF, pm_vals, doing_parametric_standard, doing_parametric_singleregressor)

if doing_parametric_standard && doing_parametric_singleregressor
    error('Cannot do both standard and single-regressor parametric modulator.')
end

% HRF can be basis set, but SINGLE basis function is entered
% this is for one regressor/event type, with multiple basis images

% must identify FIRST events, i.e., preceded by 0
% for epochs... and must get back epoch duration, too.
%wh = find(delta);
wh = find(delta & [1; diff(delta) > 0 ]);  % tor modified 6/5/16 to generalize to non-integer input vals

% wh_end values, for durs
wh_end = find(~delta & [0; diff(delta) < 0]);
if delta(end), wh_end(end+1) = length(delta); end % fix if last event goes to end

% transpose if horiz
if diff(size(pm_vals)) > 0, pm_vals = pm_vals'; end

if ~isempty(pm_vals) && length(wh) ~= length(pm_vals)
    disp('Number of events in parametric modulator does not match number of events in delta function.')
    fprintf('%3.0f events in delta function, %3.0f events in param. modulator\n', sum(delta), length(pm_vals));

    if length(pm_vals) > length(wh)
        disp('Truncating pm_vals. Could be an onset occurred after the end of the run?');
        pm_vals = pm_vals(1:length(wh));
    end

    disp('Bad input?  Stopping in debugger so you can check. Type dbcont or dbquit to continue.')
    keyboard
end

% get modulated delta function

if ~isempty(pm_vals) && doing_parametric_standard
    % create a second delta function that is modulated
    pm_vals = scale(pm_vals, 1);
    delta = repmat(delta, 1, size(pm_vals,2)+1); % Wani modified this line to fix an error when pm_vals has two or more columns.
    % delta(:, end + 1) = delta;

elseif ~isempty(pm_vals) && doing_parametric_singleregressor
    % do nothing; we will just use pm vals

elseif isempty(pm_vals)
else
    error('This should not happen');
end

% This is where they get added.
% tor adjusted aug 2012 to account for epoch case.
if ~isempty(pm_vals)
    for i = 1:length(wh)
        delta(wh(i):wh_end(i), end) = pm_vals(i);
    end
end

% Convolve all preds for this trial type
lendelta = size(delta, 1);

for i = 1:size(delta, 2)

    %for j = 1:size(HRF, 2)
    tmp = conv(double(delta(:,i)), HRF); % Changed for Matlab 7.9+ compatibility - Thanks, Liane

    X{i} = tmp(1:lendelta);
    %end

end

model = cat(2, X{:});

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [v,beta,s] = house (x)
%HOUSE find a Householder reflection.
% real or complex case.
% Example:
%   [v,beta,s] = house (x)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

v = x ;
s = norm (x) ;
if (s == 0)
    beta = 0 ;
    v (1) = 1 ;
else
    if (x (1) ~= 0)
        s = sign (x (1)) * s ;
    end
    v (1) = v (1) + s ;
    beta = 1 / real (conj (s) * v (1)) ;
end
s = - s ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%GRIDFACE Generate grid cell face numbering.
%
%   [ F, FV ] = GRIDFACE( C, B ) Generates an array F with global face
%   numbers and FV with grid point indices per face. Input data is an
%   array C with cell connectivities, and optionally the boundary
%   struct B if upper triangular (OpenFOAM) ordering should be used
%   (default normal ordering).
%
%       Input       Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       c           (n_v,n_c)              Cell connectivity pointer array where the
%                                          entries in each column correspond to the
%                                          cell vertex numbers specifying the cells.
%       b           struct                 Boundary struct to use upper triangular ordering
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       f           (n_f,n_c)              Array with global face numbers for all cells.
%       fv          (n_f,3/4)              Array with indices to vertices for each face.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function g = gammaEval(j, k, LR, N, dim, nVars)
%GAMMAEVAL   Evaluate a gamma-function.
%   g = GAMMAEVAL(L, LR, N, DIM, NVARS) evaluates the gamma-function (J, K) with
%   the contour LR, N grid points, in dimension DIM and with NVARS variables.
%
% See also EXPINTEG/GAMMAFUN.

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

% Get a function handle to the phi function of index L:
g = expinteg.gammaFun(j, k);

% Evaluate it with a contour integral:
g = mean(feval(g, LR), 2);

% Reshape it when nVars>1 or/and dim>1:
g = reshape(g, nVars*N, N^(dim>1), N^(dim>2));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function L = lt(p1,p2)
%PREAL/lt Overloaded LT (<) operator for class preal.

global useUnitsFlag

if ~(useUnitsFlag) % If physunits is disabled...
    L=lt(double(p1),double(p2)); % ... treat as double.
    return
end

p1=preal(p1);
p2=preal(p2);
tol=0.001;
if isscalar(p1)
    p1=preal(p1.value*ones(size(p2)),p1.units);
elseif isscalar(p2)
    p2=preal(p2.value*ones(size(p1)),p2.units);
elseif ndims(p1)==ndims(p2)&size(p1)==size(p2)
    % do nothing
else
    error('Array dimensions must agree')
end

L=false(size(p1));
for k=1:numel(L)
    if any(abs(p1(k).units-p2(k).units)>tol)
        error('Attempt to compare dimensionally inconsistent preals')
    end
    L(k)=lt(p1(k).value,p2(k).value);
end

author: Sable

repo: mcbench-benchmarks


text:
function [stlStruct] = import_STL_bin(fileName)

% function [stlStruct] = import_STL_bin(fileName)
% ------------------------------------------------------------------------
%
% This function reads an STL file in binary format into stlStruct
%
%
% Based on code originally written by: Doron Harlev, Eric C. Johnson
% (11-Dec-2008), and Francis Esmonde-White (May 2010).
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2016/02/24 Added to GIBBON and altered syntax and behaviour
%------------------------------------------------------------------------

%% Access file

fid=fopen(fileName, 'r'); %Open the file, assumes STL Binary format.
if fid == -1
    error('File could not be opened, check name or path.')
end

%% Read title
solidNameBin=fread(fid,80,'uchar=>schar'); % Read file title
solidName = char(solidNameBin');

%% Get number of faces
numFaces=fread(fid,1,'int32'); % Read number of Faces

%% Read remainder

T = fread(fid,inf,'uint8=>uint8'); % read the remaining values
fclose(fid); %Close file

%% Process import

% Each facet is 50 bytes
%  - Three single precision values specifying the face normal vector
%  - Three single precision values specifying the first vertex (XYZ)
%  - Three single precision values specifying the second vertex (XYZ)
%  - Three single precision values specifying the third vertex (XYZ)
%  - Two color bytes (possibly zeroed)

% 3 dimensions x 4 bytes x 4 vertices = 48 bytes for triangle vertices
% 2 bytes = color (if color is specified)

trilist = 1:48;

ind = reshape(repmat(50*(0:(numFaces-1)),[48,1]),[1,48*numFaces])+repmat(trilist,[1,numFaces]);
coordNormData = reshape(typecast(T(ind),'single'),[3,4,numFaces]);

N=squeeze(coordNormData(:,1,:))';
N=double(N);

V=coordNormData(:,2:4,:);
V=reshape(V,[3,3*numFaces]);
V=double(V)';

F=reshape(1:3*numFaces,[3,numFaces])';

c0 = typecast(T(49:50),'uint16');
if (bitget(c0(1),16)==1)
    trilist = 49:50;
    ind = reshape(repmat(50*(0:(numFaces-1)),[2,1]),[1,2*numFaces])+repmat(trilist,[1,numFaces]);
    c0 = reshape(typecast(T(ind),'uint16'),[1,numFaces]);

    r=bitshift(bitand(2^16-1, c0),-10);
    g=bitshift(bitand(2^11-1, c0),-5);
    b=bitand(2^6-1, c0);
    C=[r; g; b]';
else
    C = zeros(numFaces,3);
end

%% Arrange output in structure array
stlStruct.solidNames{1}=solidName;
stlStruct.solidVertices{1}=V;
stlStruct.solidFaces{1}=F;
stlStruct.solidNormals{1}=N;
stlStruct.solidColors{1}=C;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function F = imag( F )
%IMAG   Imaginary part of a DISKFUNV.
%   IMAG(F) returns the imaginary part of a DISKFUNV.
%
%   Since all DISKFUNV objects are real-valued, the result of this
%   function is always the original DISKFUNV F.
%
% See also DISKFUNV/CONJ, DISKFUNV/REAL.

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

% Empty check:
if ( isempty( F ) )
    return
end

% Take imaginary part of each component:
F.components = cellfun( @imag, F.components, 'UniformOutput', false );

end


text:
classdef GenericQuatInterpSteeringModel < AbstractSteeringModel
    %GenericQuatInterpSteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        t0(1,1) double = 0;
        tDur(1,1) double = 0;

        alpha0(1,1) double = 0;
        beta0(1,1) double = 0;
        gamma0(1,1) double = 0;

        alpha1(1,1) double = 0;
        beta1(1,1) double = 0;
        gamma1(1,1) double = 0;

        initQuat(1,4) double = [1,0,0,0];
        finalQuat(1,4) double = [1,0,0,0];

        angleContinuity(1,1) double = false;

        refFrame AbstractReferenceFrame
        controlFrame AbstractControlFrame
    end

    methods
        function dcm = getBody2InertialDcmAtTime(obj, ut, rVect, vVect, bodyInfo)
            if(obj.tDur > 0)
                if(ut < obj.t0)
                    f = 0;
                elseif(ut > (obj.t0 + obj.tDur))
                    f = 1;
                else
                    f = ((ut - obj.t0)/(obj.tDur));

                    if(f > 1)
                        f = 1;
                    elseif(f < 0)
                        f = 0;
                    end
                end

                q = quatinterp(obj.initQuat, obj.finalQuat, f);
            else
                q = obj.initQuat;
            end

%             elemSet = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());
%             if(not(isempty(obj.refFrame.getOriginBody())))
%                 elemSet = elemSet.convertToFrame(obj.refFrame);
%             end

            [alphaAng, betaAng, gammaAng] = quat2angle(q, 'ZYX');

            dcm = obj.controlFrame.computeDcmToInertialFrame(ut, rVect, vVect, bodyInfo, gammaAng, betaAng, alphaAng, obj.refFrame);
            dcm = real(dcm);
        end

        function t0 = getT0(obj)
            t0 = obj.t0;
        end

        function setT0(obj, newT0)
            obj.t0 = newT0;
        end

        function angleCont = getContinuityTerms(obj)
            angleCont = obj.angleContinuity;
        end

        function setContinuityTerms(obj, angleContinuity)
            obj.angleContinuity = angleContinuity;
        end

        function tDur = getDuration(obj)
            tDur = obj.tDur;
        end

        function setDuration(obj, tDur)
            obj.tDur = tDur;
        end

        function [gamma0, beta0, alpha0] = getInitAngles(obj)
            gamma0 = obj.gamma0;
            beta0 = obj.beta0;
            alpha0 = obj.alpha0;
        end

        function [gamma1, beta1, alpha1] = getFinalAngles(obj)
            gamma1 = obj.gamma1;
            beta1 = obj.beta1;
            alpha1 = obj.alpha1;
        end

        function setInitAngles(obj, gamma0, beta0, alpha0)
            obj.gamma0 = gamma0;
            obj.beta0 = beta0;
            obj.alpha0 = alpha0;
        end

        function setFinalAngles(obj, gamma1, beta1, alpha1)
            obj.gamma1 = gamma1;
            obj.beta1 = beta1;
            obj.alpha1 = alpha1;
        end

        function setConstsFromDcmAndContinuitySettings(obj, dcm, ut, rVect, vVect, bodyInfo)
            if(obj.angleContinuity)
                elemSet = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());

%                 if(not(isempty(obj.refFrame.getOriginBody())))
%                     elemSet = elemSet.convertToFrame(obj.refFrame);
%                 end

                [gammaAngle, betaAngle, alphaAngle] = obj.controlFrame.getAnglesFromInertialBodyAxes(LaunchVehicleAttitudeState(dcm), elemSet.time, elemSet.rVect(:), elemSet.vVect(:), elemSet.frame.getOriginBody(), obj.refFrame);

                obj.alpha0 = alphaAngle;
                obj.beta0 = betaAngle;
                obj.gamma0 = gammaAngle;
            end

            obj.setQuatsByAngles(); %need this here because I need to initialize the quaternions correctly from the angles at action init
        end

        function setInitialAttitudeFromState(obj, stateLogEntry, tOffsetDelta)
            %nothing
        end

        function [angle1Name, angle2Name, angle3Name] = getAngleNames(obj)
%             angle1Name = 'Gamma';
%             angle2Name = 'Beta';
%             angle3Name = 'Alpha';

            angleNames = obj.controlFrame.getControlFrameEnum().angleNames;
            angle1Name = angleNames{1};
            angle2Name = angleNames{2};
            angle3Name = angleNames{3};
        end

        function newSteeringModel = deepCopy(obj)
            newSteeringModel = GenericQuatInterpSteeringModel(obj.t0, obj.tDur, obj.alpha0, obj.beta0, obj.gamma0, obj.alpha1, obj.beta1, obj.gamma1);

            newSteeringModel.angleContinuity = obj.angleContinuity;

            newSteeringModel.refFrame = obj.refFrame;
            newSteeringModel.controlFrame = obj.controlFrame;
        end

        function optVar = getNewOptVar(obj)
            optVar = SetGenericQuatInterpSteeringModelActionOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesRefFrame(~)
            tf = true;
        end

        function refFrame = getRefFrame(obj)
            refFrame = obj.refFrame;
        end

        function setRefFrame(obj, refFrame)
            obj.refFrame = refFrame;
        end

        function tf = usesControlFrame(~)
            tf = true;
        end

        function cFrame = getControlFrame(obj)
            cFrame = obj.controlFrame;
        end

        function setControlFrame(obj, cFrame)
            obj.controlFrame = cFrame;
        end

        function tf = requiresReInitAfterSoIChange(obj)
            tf = false;
        end

        function enum = getSteeringModelTypeEnum(~)
            enum = SteerModelTypeEnum.QuaterionInterp;
        end

        function [addActionTf, steeringModel] = openEditSteeringModelUI(obj, lv, useContinuity)
            output = AppDesignerGUIOutput({false, obj});
            lvd_EditActionSetQuatInterpSteeringModelGUI_App(obj, lv, useContinuity, output);
            addActionTf = output.output{1};
            steeringModel = output.output{2};
        end
    end

    methods(Access=private)
        function obj = GenericQuatInterpSteeringModel(t0, tDur, alpha0, beta0, gamma0, alpha1, beta1, gamma1)
            obj.t0 = t0;
            obj.tDur = tDur;

            obj.alpha0 = alpha0;
            obj.beta0 = beta0;
            obj.gamma0 = gamma0;

            obj.alpha1 = alpha1;
            obj.beta1 = beta1;
            obj.gamma1 = gamma1;

            obj.setQuatsByAngles();
        end

        function setQuatsByAngles(obj)
            obj.initQuat = angle2quat(obj.alpha0, obj.beta0, obj.gamma0, 'ZYX');
            obj.finalQuat = angle2quat(obj.alpha1, obj.beta1, obj.gamma1, 'ZYX');
        end
    end

    methods(Static)
        function model = getDefaultSteeringModel()
            t0 = 0;
            tDur = 1;

            alpha0 = 0;
            beta0 = 0;
            gamma0 = 0;

            alpha1 = 0;
            beta1 = 0;
            gamma1 = 0;

            model = GenericQuatInterpSteeringModel(t0, tDur, alpha0, beta0, gamma0, alpha1, beta1, gamma1);
        end

        function typeStr = getTypeNameStr()
            typeStr = SteeringModelEnum.GenericQuatInterp.nameStr;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  [T]= direct_kinematics_3RRR_numerical(ROBOT, THETA, threshold, max_iter)
%
%  Compute the direct kinematic problem in position for a 3RRR arm.
%  Returns the matrix T with position/orientation of the end effector.
%
%   T=[cos(Phi) -sin(Phi) 0 xA;
%      sin(Phi) cos(Phi) 0  yA;
%        0         0     1   0;
%        0          0    0   1];
%
%   The found solution of joint coordinates is drawn.
%
%  Input parameters:
%  ROBOT: robot structure (parallel).
%  THETA: The input vector of active joint coordinates THETA = [theta1
%           theta2 theta3];
%  Threshold: Target sum of squares. When the sum of squares of the
%           function values is below this value the method stops (default 0.01).
%  Max_iter: Maximum number of iteration steps in the algorithm (default 50).
%
%  WARNING: The geometric parameters have been fixed in this function, they
%  should be extracted from the robot structure.
%
%
% THE FUNCTION SHOWS HOW TO SOLVE THE DIRECT KINEMATICS OF THE 3RRR USING A
% NUMERICAL METHOD. THIS METHOD MAY BE APPLIED TO ANY PARALELL ROBOT FOR WHICH YOU
% MAY NOT FIND A CLOSED SOLUTION. OF COURSE, BE AWARE THAT A CLOSED
% SOLUTION EXISTS FOR THE 3RRR PARALLEL ROBOT, HOWEVER, FOR SIMPLICITY, WE
% APPLY THE METHOD TO A WELL-KNOWN, SIMPLE ROBOT.
%
% The solution uses a vanilla Newton-Raphson method, also known as a
% Levenberg-Marquardt method with Lambda=1.
%
%   Please, for more information see:
%   http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
%
% Example:
%   >>robot=load_robot('example','3RRR');
%   >>T=directkinematics_3RRR_numerical(robot, [pi/2 pi/4 -pi/2])
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   20/01/2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

function T=direct_kinematics_3RRR_numerical(robot, theta, randomize, threshold, max_iter)

global xQ yQ xR yR b1 b2 b3 a1 a2 a3 h
close all
L=3; % basic mechanism length
%position of the fixed points xQ, yQ, xR, yR
xQ=L; yQ=0; xR=L/2; yR=L;

% Geometry
b1=1; b2=1; b3=1; a1=1; a2=1;a3=1; % m
%end effector is a triangle of side h
h=0.5; %m

%check convergence parameters
if ~exist('threshold','var')
    threshold=0.01;
end

if ~exist('max_iter','var')
    max_iter=50;
end

%randomize solutions
if ~exist('randomize','var')
    randomize=0;
end

%input parameters for the direct kinematics
th1=theta(1);%1.9948;
th2=theta(2);%3.2869;
th3=theta(3);%-1.3998;

theta = [th1 th2 th3];

%initial beta, different solutions may be found, depending on the
%starting values for the parameters beta
% beta = [xA yA Phi phi1 phi2 phi3]
beta=[0.0 0.0 -pi/4 0 0 0]';

%randomize solution
if randomize
    beta = beta + randn(1,6)';
end

%the vector y is the vector of values that the Gamma equations should have
%in this case, they should be all null
% Gamma=[Gamma1, Gamma2, Gamma3..., Gamma6];
y=[0 0 0 0 0 0]';

fs=[];

for i=1:max_iter,
    f=compute_gamma_values(beta,theta);
    J=compute_jacobian_beta(beta,theta);

    delta=inv(J'*J)*J'*(y-f);

    %update beta
    beta = beta + delta;

    sumofsqes=(y-f)'*(y-f);
    fs=[fs sumofsqes];

    if sumofsqes < threshold
        fprintf('\ndirect_kinematics_3RRR::Solution found in %d iterations',i);
        fprintf('\nPlease note that only ONE possible solution is returned');
        break;
    end
end

%check whether the algorithm converged in less than max_iter steps
if i==max_iter
  fprintf('ERROR:: direct_kinematics_3RRR::Could not converge in %d iterations',i);
end

fprintf('FINAL VALUE of Gamma function');
f

fprintf('FINAL VALUE of sumof squares');
sumofsqes

%build solution from the beta parameters.
T=eye(4);
T(1,4)=beta(1);
T(2,4)=beta(2);
T(1,1)=cos(beta(3));
T(2,1)=sin(beta(3));
T(1,2)=-sin(beta(3));
T(2,2)=cos(beta(3));

%draw the given solution
%q=[th1 phi1 th2 phi2 th3 phi3]
q=[th1 beta(4) th2 beta(5) th3 beta(6)];
drawrobot3d(robot, q), pause(2);

figure, plot(fs), title('Sum of squares vs. iteration step')

%compute the gamma functions, each function fi corresponds to a loop closing
%equation in the 3RRR mechanism, given the current estimates of beta
%the theta (theta1, theta2, theta3) may be considered here as constants,
%since they are known values for the direct kinematic solution
function f=compute_gamma_values(beta, theta)

global xQ yQ xR yR b1 b2 b3 a1 a2 a3 h

%just assign the values to some variables, so that the equations are easier
%to read
xA=beta(1);
yA=beta(2);
Phi=beta(3);
phi1=beta(4);
phi2=beta(5);
phi3=beta(6);

%known values of the direct kinematic problems
th1=theta(1);
th2=theta(2);
th3=theta(3);

%Loop closing equations for the 3RRR mechanism
f1=eval('xA-a1*cos(th1)-b1*cos(th1+phi1)');
f2=eval('yA-a1*sin(th1)-b1*sin(th1+phi1)');
f3=eval('xA-xQ-a2*cos(th2)-b2*cos(th2+phi2)+h*cos(Phi)');
f4=eval('yA-yQ-a2*sin(th2)-b2*sin(th2+phi2)+h*sin(Phi)');
f5=eval('xA-xR-a3*cos(th3)-b3*cos(th3+phi3)+h*cos(Phi+pi/3)');
f6=eval('yA-yR-a3*sin(th3)-b3*sin(th3+phi3)+h*sin(Phi+pi/3)');

%form vector with values,
%note that the f functions are the loop closing equations that should be
%equal to zero
f=[f1 f2 f3 f4 f5 f6]';

%Compute Jacobian with respect to beta parameters
%xA yA Phi, ph1, ph2, ph3
function J=compute_jacobian_beta(beta, theta)

global xQ yQ xR yR b1 b2 b3 a1 a2 a3 h

%just assign the values to some variables, so that the equations are easier
%to read
%xA=beta(1);
%yA=beta(2);
Phi=beta(3);
ph1=beta(4);
ph2=beta(5);
ph3=beta(6);

%known values of the direct kinematic problems
th1=theta(1);
th2=theta(2);
th3=theta(3);

%Compute Jacobian
%Gamma1 with respect to xA yA Phi, ph1, ph2, ph3
J11=1;
J12=0;
J13=0;
J14=b1*sin(th1+ph1);
J15=0;
J16=0;

%Gamma2 with respect to xA yA Phi, ph1, ph2, ph3
J21=0;
J22=1;
J23=0;
J24=-b1*cos(th1+ph1);
J25=0;
J26=0;

%Gamma3 with respect to xA yA Phi, ph1, ph2, ph3
J31=1;
J32=0;
J33=-h*sin(Phi);
J34=0;
J35=b2*sin(th2+ph2);
J36=0;

%Gamma4 with respect to xA yA Phi, ph1, ph2, ph3
J41=0;
J42=1;
J43=h*cos(Phi);
J44=0;
J45=-b2*cos(th2+ph2);
J46=0;

%Gamma5 with respect to xA yA Phi, ph1, ph2, ph3
J51=1;
J52=0;
J53=-h*sin(Phi+pi/3);
J54=0;
J55=0;
J56=b3*sin(th3+ph3);

%Gamma6 with respect to xA yA Phi, ph1, ph2, ph3
J61=0;
J62=1;
J63=h*cos(Phi+pi/3);
J64=0;
J65=0;
J66=-b3*cos(th3+ph3);

J=[J11 J12 J13 J14 J15 J16;
   J21 J22 J23 J24 J25 J26;
   J31 J32 J33 J34 J35 J36;
   J41 J42 J43 J44 J45 J46;
   J51 J52 J53 J54 J55 J56;
   J61 J62 J63 J64 J65 J66];

repo: ARTE


text:
function found = roiExistName(view,roiName,fullMatch)
%
%      found = roiExistName(view,roiName,[fullMatch])
%
% Author: Wandell
% Purpose:
%
%   Check whether an ROI with roiName exists in the view.
%
%   If a full match to the roiName is not found, then we return
%   found=0.
%
%   If a full match is found, view.ROIs(found) is the roi.
%
%   If a partial match (fullMatch = 0) is requested, then roiExistName() can
%   find several matches (i.e., whenever the first part of the roiName matches
%   the roiName string). In this case, found returns a list of indices
%   whose names match the dataTypeName.  If no matches are found, then
%   found = 0 is returned.
%
%   Examples:
%       roiExistName(FLAT{1},'Quad',0) returns a list of partial matches
%       FLAT{1}.ROIs(found).name all begin with the phrase 'Quad', such as
%       'Quad-1'.
%

if ~exist('fullMatch','var'), fullMatch = 1; end
ROIs = view.ROIs;

nROIs = length(ROIs);

if fullMatch
    found = 0;
    for ii=1:nROIs
        if strcmp(ROIs(ii).name,roiName);
            found = ii;
            return;
        end
    end
else
    found = zeros(1,nROIs);
    for ii=1:nROIs
        found(ii) =  strncmp(ROIs(ii).name,roiName,length(roiName));
    end
    found = find(found);
    if isempty(found), found = 0; end
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%MDL_BAXTER Kinematic model of Baxter dual-arm robot
%
% MDL_BAXTER is a script that creates the workspace variables left and
% right which describes the kinematic characteristics of the two 7-joint
% arms of a Rethink Robotics Baxter robot using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%   qd         lower arm horizontal as per data sheet
%
% Notes::
% - SI units of metres are used.
%
% References::
% "Kinematics Modeling and Experimental Verification of Baxter Robot"
% Z. Ju, C. Yang, H. Ma, Chinese Control Conf, 2015.
%
% See also mdl_nao, SerialLink.

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

% MODEL: Baxter, Rethink Robotics, 7DOF, standard_DH

%     th  d  a  alpha

links = [
        Revolute('d', 0.27,        'a', 0.069, 'alpha', -pi/2)
        Revolute('d', 0,           'a', 0, 'alpha', pi/2, 'offset', pi/2)
        Revolute('d', 0.102+0.262, 'a', 0.069, 'alpha', -pi/2)
        Revolute('d', 0,           'a', 0, 'alpha', pi/2)
        Revolute('d', 0.103+0.271, 'a', 0.010, 'alpha', -pi/2)
        Revolute('d', 0,           'a', 0, 'alpha', pi/2)
        Revolute('d', 0.28,        'a', 0, 'alpha', 0)
];

left =  SerialLink(links, 'name', 'Baxter LEFT', 'manufacturer', 'Rethink Robotics');
right = SerialLink(links, 'name', 'Baxter RIGHT', 'manufacturer', 'Rethink Robotics');

left.base = transl(0.064614, 0.25858, 0.119)*rpy2tr(0, 0, pi/4, 'xyz');
right.base = transl(0.063534, -0.25966, 0.119)*rpy2tr(0, 0, -pi/4, 'xyz');

% define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%   qstretch   arm is stretched out in the X direction
%   qn         arm is at a nominal non-singular configuration
%
qz = [0 0 0 0 0 0 0]; % zero angles, L shaped pose
qr = [0 -pi/2 -pi/2 0 0 0 0]; % ready pose, arm up
qs = [0 0 -pi/2 0 0 0 0];
qn = [0 pi/4 pi/2 0 pi/4  0 0];

author: petercorke

repo: robotics-toolbox-matlab


text:
function [state_red,inv_state_red,g0,g1,c,pi,psi] = clean_G0(g0,g1,c,pi,psi)
% Solves out static constraints of the linear model using SVD
%    Adapted from Chris Sims's gensys code
%
% Input/output/References: It will be faster to read the codes below
%
% by SeHyoun Ahn, June 2016
%
% SYNTAX:
% [state_red,inv_state_red,g0,g1,c,pi,psi] = clean_G0(g0,g1,c,pi,psi)
tmp=(max(abs([g0,psi]),[],2)==0);
redundant=find(tmp);
keep=find(1-tmp);
inv_state_red=sparse(null(full(g1(redundant,:))));

g0=inv_state_red'*g0*inv_state_red;
g1=inv_state_red'*g1*inv_state_red;
g1=g0\g1;
psi=g0\inv_state_red'*psi;
pi=g0\inv_state_red'*pi;
c=g0\inv_state_red'*c;

state_red = inv_state_red';

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function model=setParam(model, paramType, rxnList, params, var)
% setParam
%   Sets parameters for reactions
%
%   model       a model structure
%   paramType   the type of parameter to set:
%               'lb'    Lower bound
%               'ub'    Upper bound
%               'eq'    Both upper and lower bound (equality
%                       constraint)
%               'obj'   Objective coefficient
%               'rev'   Reversibility
%               'var'   Variance around measured bound
%   rxnList     a cell array of reaction IDs or a vector with their
%               corresponding indexes
%   params      a vector of the corresponding values
%   var         percentage of variance around measured value, if 'var' is
%               set as paramType. Defining 'var' as 5 results in lb and ub
%               at 97.5% and 102.5% of the provide params value (if params
%               value is negative, then lb and ub are 102.5% and 97.5%).
%
%   model       an updated model structure
%
%   Usage: model=setParam(model, paramType, rxnList, params)

paramType=convertCharArray(paramType);
if ~any(strcmpi(paramType,{'lb','ub','eq','obj','rev','var'}))
    EM=['Incorrect parameter type: "' paramType '"'];
    dispEM(EM);
end
if isnumeric(rxnList) || islogical(rxnList)
    rxnList=model.rxns(rxnList);
else
    rxnList=convertCharArray(rxnList);
end

%Allow to set several parameters to the same value
if numel(rxnList)~=numel(params) && numel(params)~=1
    EM='The number of parameter values and the number of reactions must be the same';
    dispEM(EM);
end

if length(rxnList)>1
    if length(paramType)==1
        paramType(1:length(rxnList))=paramType;
    end
    if length(params)==1
        params(1:length(rxnList))=params;
    end
end

%Find the indexes for the reactions in rxnList. Do not use getIndexes
%as we do not want to throw errors if matches fail
indexes=zeros(numel(rxnList),1);

for i=1:numel(rxnList)
    index=find(strcmp(rxnList{i},model.rxns),1);
    if ~isempty(index)
        indexes(i)=index;
    else
        indexes(i)=-1;
        EM=['Reaction ' rxnList{i} ' is not present in the reaction list'];
        dispEM(EM,false);
    end
end

%Remove the reactions that were not found
params(indexes==-1)=[];
indexes(indexes==-1)=[];
paramType(indexes==-1)=[];
%Change the parameters

if ~isempty(indexes)
    if contains(paramType,'obj')
        model.c=zeros(numel(model.c),1); % parameter is changed, not added
    end
    for j=1:length(paramType)
        if strcmpi(paramType{j},'eq')
            model.lb(indexes(j))=params(j);
            model.ub(indexes(j))=params(j);
        end
        if strcmpi(paramType{j},'lb')
            model.lb(indexes(j))=params(j);
        end
        if strcmpi(paramType{j},'ub')
            model.ub(indexes(j))=params(j);
        end
        if strcmpi(paramType{j},'obj')
            model.c(indexes(j))=params(j);
        end
        if strcmpi(paramType{j},'rev')
            %Non-zero values are interpreted as reversible
            model.rev(indexes(j))=params(j)~=0;
        end
        if strcmpi(paramType{j},'var')
            if params(j) < 0
                model.lb(indexes(j)) = params(j) * (1+var/200);
                model.ub(indexes(j)) = params(j) * (1-var/200);
            else
                model.lb(indexes(j)) = params(j) * (1-var/200);
                model.ub(indexes(j)) = params(j) * (1+var/200);
            end
        end
    end
end
end

author: SysBioChalmers

repo: RAVEN


text:
function [t,s]=zerocros(y,m,x)
%ZEROCROS finds the zeros crossings in a signal [T,S]=(Y,M,X)
% Inputs:  y = input waveform
%          m = mode string containing:
%              'p' - positive crossings only
%              'n' - negative crossings only
%              'b' - both (default)
%              'r' - round to sample values
%          x = x-axis values corresponding to y [default 1:length(y)]
%
% Outputs: t = x-axis positions of zero crossings
%          s = estimated slope of y at the zero crossing
%
% This routine uses linear interpolation to estimate the position of a zero crossing
% A zero crossing occurs between y(n) and y(n+1) iff (y(n)>=0) ~= (y(n+1)>=0)

% Example: y=sin(2*pi*(0:1000)/200); y(1:100:1001)=0; zerocros(y);
% Note that we get a zero crossing at the end but not at the start.

%	   Copyright (C) Mike Brookes 2003-2015
%      Version: $Id: zerocros.m 6563 2015-08-16 16:56:24Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 || ~numel(m)
    m='b';
end
s=y>=0;
k=s(2:end)-s(1:end-1);
if any(m=='p')
    f=find(k>0);
elseif any(m=='n')
    f=find(k<0);
else
    f=find(k~=0);
end
s=y(f+1)-y(f);
t=f-y(f)./s;
if any(m=='r')
    t=round(t);
end
if nargin>2
    tf=t-f; % fractional sample
    t=x(f).*(1-tf)+x(f+1).*tf;
    s=s./(x(f+1)-x(f));
end
if ~nargout
    n=length(y);
    plot(1:n,y,'-',t,zeros(length(t),1),'o');
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [ori, Tori,i] = project(obj,ori)
% project an embedding back onto the manifold of orientations.
%
% Syntax
%   [ori, Tori] = project(e)
%
% Input
%  e - @embedding
%
% Output
%  ori - @orientation
%  Tori - projected @embedding
%

% get the embedding of the identity
Tid = embedding.id(obj.CS);

% normalize obj correctly
obj = obj .* norm(Tid)./norm(obj);

%get weights beta
%[~,~,weights]= embedding.coefficients(obj.CS);

% ensure obj is symmetric
%obj = obj.sym;

% special case for triclinic symmetry
if 0 && obj.CS.Laue.id ==2

  %weighted sum in Horn
  for i = 1:length(obj.u)
    r(i,:) = obj.u{i}(:);
  end

  ori = orientation(fit(obj.l,r,obj.CS),obj.CS);

  if nargout == 2, Tori = ori * Tid; end
  return

end

% initial guess - we need it to be sufficently close to avoid local extrema
if nargin == 1
  ori = equispacedSO3Grid(obj.CS,'points',10);

  d = norm(reshape(obj,[],1) - (ori * Tid).');

  [~,id] = min(d,[],2);

  ori = reshape(ori(id),[],1);
  %ori = orientation.id(obj.CS);
end

% basis of the tangential space
t = spinTensor([xvector,yvector,zvector]);

% perform steepest descent iteration to maximize dot(ori * Tid, obj)
maxIter = 200;
for i = 1:maxIter

  % compute the gradient in ori
  Tori = rotate_outer(Tid,ori);
  for k = 1:length(obj.u)
    Tori.u{k} = obj.rank(k) * EinsteinSum(t,[1,-1],Tori.u{k},[-1 2:obj.rank(k)]);
  end

  g = vector3d(dot(Tori,obj).').';

  % eradicate normalizing of embedding: adapt length ofs gradient
  %g = g * obj.rho^2;

  % stop if gradient is sufficently small
  if all(norm(g)<1e-10), break; end
  %disp([xnum2str(max(norm(g))) ' ' char(ori(1)) ' ' char(g(1))]);

  % update ori
  ori = exp(ori,g,'left');

end
end

author: mtex-toolbox

repo: mtex

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


text:
function i4_is_even_test ( )

%*****************************************************************************80
%
%% I4_IS_EVEN_TEST tests I4_IS_EVEN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 December 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_IS_EVEN_TEST\n' );
  fprintf ( 1, '  I4_IS_EVEN reports whether an I4 is even.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I     I4_IS_EVEN(I)\n' );
  fprintf ( 1, '\n' );

  for i = -2 : 25
    fprintf ( 1, '  %6d  %1d\n', i, i4_is_even ( i ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [mGal] = mmps22mGal(mmps2)
% Convert acceleration from millimeters per square-second to milligalileos
% Chad A. Greene 2012
mGal = mmps2*1e+2;


text:
function dimacs = computedimacs(b,c,A,xin,y,s,K);
% COMPUTEDIMACS
%
% min <C,X> s.t     AX = b, X > 0
% max b'y   s.t S-C+A'y =0, S > 0

% If no primal exist, fake till later
if isempty(xin)
    x = c*0;
else
    x = xin;
end

if isempty(s)
    s = c-A'*y;
end

xres = inf;
sres = inf;

% Not officially defined in DIMACS
if K.f>0
    sres = -min(norm(s(1:K.f),inf));
end

% Errors in linear cone
if K.l>0
    xres = min(x(1+K.f:K.f+K.l));
    sres = min(s(1+K.f:K.f+K.l));
end

% Errors in quadratic cone
if K.q(1)>0
    top = K.f+K.l;
    for i = 1:length(K.q)
        X = x(1+top:top+K.q(i));
        S = s(1+top:top+K.q(i));
        xres = min(xres,X(1)-norm(X(2:end)));
        sres = min(sres,S(1)-norm(S(2:end)));
        top = top + K.q(i);
    end
end

% Errors in semidefinite cone
if K.s(1)>0
    top = K.f+K.l+K.q+K.r;
    for i = 1:length(K.s)
        X = reshape(x(1+top:top+K.s(i)^2),K.s(i),K.s(i));
        S = reshape(s(1+top:top+K.s(i)^2),K.s(i),K.s(i));
        xres = min(xres,min(eig(full(X))));
        sres = min(sres,min(eig(full(S))));
        top = top + K.s(i)^2;
    end
end

err1 = norm(b-A*x)/(1 + norm(b,inf));
err2 = max(0,-xres)/(1 + norm(b,inf));
err3 = conenorm(s-(c-A'*y),K)/(1+norm(c,inf));
%err3 = norm(s-(c-A'*y))/(1+norm(c,inf)); % Used by some solvers
err4 = max(0,-sres)/(1+max(abs(c)));
err5 = (c'*x-b'*y)/(1+abs(c'*x)+abs(b'*y));
err6 = x'*(c-A'*y)/(1+abs(c'*x)+abs(b'*y));

% No primal was computed
if isempty(xin)
    err1 = nan;
    err2 = nan;
    err5 = nan;
    err6 = nan;
end
dimacs = [err1 err2 err3 err4 err5 err6];

function t = conenorm(s,K)

% Implementation of the norm described on
% http://plato.asu.edu/dimacs/node3.html

t = 0;

if K.f + K.l>0
    t = t + norm(s(1:K.f+K.l));
end

top = 1+K.f+K.l;
if K.q(1)>0
    for i = 1:length(K.q)
        t = t + norm(s(top:top+K.q(i)-1));
        top  = top + K.q(i);
    end
end

if K.s(1)>0
    for i = 1:length(K.s)
        S = reshape(s(top:top+K.s(i)^2-1),K.s(i),K.s(i));
        t = t + norm(S,'fro');
        top  = top + K.s(i)^2;
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Kernel Affine Projection Subgradient Method
%
% K. Slavakis, S. Theodoridis, and I. Yamada, "Online kernel-based
% classification using adaptive projection algorithms," IEEE Transactions
% on Signal Processing, Vol. 56, No. 7, pp. 2781-2796, 2008.
% http://dx.doi.org/10.1109/TSP.2008.917376
%
% Remark: implemented with L2-ball forgetting. Code contributed by
% Pantelis Bouboulis.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef kapsm < kernel_adaptive_filter

    properties (GetAccess = 'public', SetAccess = 'private') % parameters
        M = 200; % dictionary size
        Delta = 10; % L2-ball forgetting radius
        mu = 1.8, % learning rate
        Q = 10; % number of subgradients
        loss = 'l2'; % loss function type
        loss_param = 2; % Huber loss parameter
        kerneltype = 'gauss'; % kernel type
        kernelpar = .5; % kernel parameter
    end

    properties (GetAccess = 'public', SetAccess = 'private') % variables
        dict = []; % dictionary
        alpha = []; % expansion coefficients
        norm_f = []; % function norm for L2-ball forgetting
        xmem = []; % input memory
        ymem = []; % output memory
    end

    methods
        function kaf = kapsm(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(kaf))
                        kaf.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(kaf,x) % evaluate the algorithm
            if size(kaf.dict,1)>0
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                y_est = k'*kaf.alpha;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(kaf,x,y) % train the algorithm
            if size(kaf.xmem,1)<kaf.Q
                % grow memory
                kaf.xmem = [kaf.xmem; x];
                kaf.ymem = [kaf.ymem; y];
            else
                % slide memory
                kaf.xmem = [kaf.xmem(2:kaf.Q,:); x];
                kaf.ymem = [kaf.ymem(2:kaf.Q); y];
            end
            q = size(kaf.xmem,1);

            d_hat = kaf.evaluate(kaf.xmem);
            e = kaf.ymem - d_hat;

            L_n = kaf.loss_fun(e,kaf.loss,kaf.loss_param);
            subgrad = kaf.compute_subgrad_coef(e,kaf.loss,kaf.loss_param);
            norm_subgrad = subgrad.^2 .* ...
                kernel(kaf.xmem(end-q+1:end,:),...
                kaf.xmem(end-q+1:end,:),...
                [kaf.kerneltype '-diag'],kaf.kernelpar);

            omega = 1/q*ones(q,1); % convex weights

            beta = omega.*L_n.*subgrad./(norm_subgrad+eps);

            % % extrapolation parameter
            % nominator = (omega.*L_n)'*(L_n./(norm_subgrad+eps));
            % K = kernel(kaf.xmem,kaf.xmem,kaf.kerneltype,kaf.kernelpar);
            % denominator = beta'*K*beta;
            % Mn = nominator/(denominator+eps);
            % mu = .5*Mn;

            kaf.alpha = [kaf.alpha; 0];
            kaf.alpha(end-q+1:end) = kaf.alpha(end-q+1:end) - kaf.mu*beta;
            kaf.dict = [kaf.dict; x];

            % L2-ball forgetting: update function norm
            K = kernel(kaf.xmem,kaf.xmem,kaf.kerneltype,kaf.kernelpar);
            R_sum = beta(end-q+1:end)'*K*beta(end-q+1:end);
            sum2 = beta(1:end-q)'*d_hat(1:end-q);
            kaf.norm_f = sqrt(kaf.norm_f*kaf.norm_f + 2*sum2 + R_sum);
            if kaf.norm_f > kaf.Delta
                kaf.alpha = kaf.Delta/kaf.norm_f * kaf.alpha;
                kaf.norm_f = kaf.Delta;
            end

            % sliding-window pruning
            if length(kaf.alpha)>kaf.M
                k = kernel(kaf.dict,kaf.dict(1,:),...
                    kaf.kerneltype,kaf.kernelpar);
                kaf.norm_f = sqrt(kaf.norm_f.^2 + kaf.alpha(1)^2*k(1) - ...
                    2*kaf.alpha(1)*k'*kaf.alpha);
                kaf.dict(1,:) = [];
                kaf.alpha(1) = [];
            end
        end

    end

    methods (Static = true)

        function loss = loss_fun(ksi,loss_type,loss_param)
            switch loss_type
                case 'l2'
                    loss = ksi.^2;
                case 'l1'
                    loss = abs(ksi);
                case 'huber'
                    sigma = loss_param;
                    loss = zeros(length(ksi),1);
                    for i=1:length(ksi)
                        if abs(ksi(i)) <= sigma
                            loss(i) = ksi(i)^2/(2*sigma);
                        else
                            loss(i) = abs(ksi(i)) - sigma/2;
                        end
                    end
            end
        end

        function subgrad_coef = ...
                compute_subgrad_coef(e,loss_type,loss_param)
            subgrad_coef = zeros(length(e),1);
            for i=1:length(e)
                switch loss_type
                    case 'l2'
                        subgrad_coef(i) = -2*e(i);
                    case 'l1'
                        subgrad_coef(i) = -sign(e(i));
                    case 'huber'
                        sigma = loss_param;
                        if abs(e(i)) <= sigma
                            subgrad_coef(i) = -e(i)/sigma;
                        else
                            subgrad_coef(i) = -sign(e(i));
                        end
                end
            end
        end
    end
end

repo: kafbox


text:
function clearex(varargin)
%% This function clear all workspace
% except for one or more selected variable
%
% File created on Feb 15th 2008
%
% Last modified: Feb 15th 2008
%
% Author: Arnaud Laurent
%
%   Inputs: name of variables to keep (e.g. 'a','b','c')
%   Note: It is possible to use wildcard (e.g. 'a*')

a = evalin('base','who');
var = cell(size(varargin));
for i=1:nargin
    var{i}=varargin{i};
end
assignin('base','ClEaRGsJioU',var);
var = evalin('base','who(ClEaRGsJioU{:})');
clearvar = a(~ismember(a,var));
assignin('base','ClEaRGsJioU',clearvar);
evalin('base','clear(ClEaRGsJioU{:},''ClEaRGsJioU'')')

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [fluxROOM, solutionROOM, totalFluxDiff] = linearROOM(model, fluxWT, rxnKO, varargin)
% Performs a LP version of the ROOM (Regulatory on/off minimization of
% metabolic flux changes) approach
%
% USAGE:
%
%    [fluxROOM, solutionROOM, totalFluxDiff] = linearROOM(model, WTflux, rxnKO, delta, epsilon, printLevel)
%
% INPUTS:
%    model:            Metabolic model
%    fluxWT:           Numeric array with flux distribution of wild type
%    rxnKO:            List of perturbations performed to the model
%                      (reactions that are eliminated)
%
% OPTIONAL INPUTS:
%    delta:            Multiplicative tol for flux change (Default = 0.03)
%    epsilon:          Additive tolerance for flux change (Default = 0.001)
%    printLevel:       Verbose output (Default = 1)
%
% OUTPUTS:
%    fluxROOM:         Flux distribution after ROOM calculation
%    solutionROOM:     Solution structure
%    totalFluxDiff:    Euclidean distance of ROOM objective, i.e.
%                      :math:`\sum (v_{wt}-v_{del})^2`
%
% Solve the following problem:
%
% .. math::
%      min ~&~ \sum y_{i} \\
%          ~&~ S_{del}v_{del} = 0 \\
%          ~&~ lb_{del} \leq v_{del} \leq ub_{del} \\
%          ~&~ for i=1:nRxns\\
%          ~&~ v_{i} - y_{i}(v_{max,i}-w_{wt,i}^u) \leq w_{wt,i}^u \\
%          ~&~ v_{i} - y_{i}(v_{min,i}-w_{wt,i}^l) \geq w_{wt,i}^l \\
%          ~&~ 0 \leq y_{i} \leq 1 \\
%          ~&~ w_{wt,i}^u = w_{wt,i} + \delta |w_{wt,i}| + \epsilon \\
%          ~&~ w_{wt,i}^l = w_{wt,i} - \delta |w_{wt,i}| - \epsilon \\
%
% NOTE::
%
%    The code here has been based on:
%    Shlomi, T., Berkman, O., & Ruppin, E. (2005). Regulatory on/off
%    minimization of metabolic flux changes after genetic perturbations.
%    Proceedings of the National Academy of Sciences, 102(21), 7695-7700
%
% .. Authors:
%       - Luis V. Valcarcel, 23/01/2019, University of Navarra, CIMA & TECNUN School of Engineering.

p = inputParser;
% check required arguments
addRequired(p, 'model');
addRequired(p, 'WTflux', @(x)isnumeric(x)&&isvector(x));
addRequired(p, 'rxnKO', @(x)iscell(x));
% Check optional arguments
addParameter(p, 'delta', 0.03, @(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'epsilon', 0.001, @(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'printLevel', 1, @(x)isnumeric(x)&&isscalar(x));
% extract variables from parser
parse(p, model, fluxWT, rxnKO, varargin{:});
delta = p.Results.delta;
epsilon = p.Results.epsilon;
printLevel = p.Results.printLevel;

% LP solution tolerance
global CBT_LP_PARAMS
if (exist('CBT_LP_PARAMS', 'var'))
    if isfield(CBT_LP_PARAMS, 'objTol')
        tol = CBT_LP_PARAMS.objTol;
    else
        tol = 1e-6;
    end
else
    tol = 1e-6;
end

% Check the inputs
fluxWT = reshape(fluxWT,[],1); % reshape flux vector
assert(numel(model.rxns)==numel(fluxWT), 'This flux distribution has different number of reactions than the model')
assert(norm(model.S * fluxWT)< 10*tol, 'This flux distribution cannot exist in this model')
assert(all(ismember(rxnKO, model.rxns)), 'Some reactions are not in the model')

% Eliminate almost-zero fluxes
fluxWT(abs(fluxWT)<tol) = 0;

% generate auxiliary variables
WT_upperTol = fluxWT + delta*abs(fluxWT) + epsilon;
WT_lowerTol = fluxWT - delta*abs(fluxWT) - epsilon;

%auxiliary variables
[nMets, nRxns] = size(model.S);

% Generate model
ROOMmodel = model;
new_variables = strcat('binary_',model.rxns);
% add novel variables
ROOMmodel = addCOBRAVariables(ROOMmodel,new_variables,'lb',0,'ub',1);
% add novel constrains
for i=1:nRxns
    ROOMmodel = addCOBRAConstraints(ROOMmodel, {model.rxns{i}, strcat('binary_',model.rxns{i})}, WT_upperTol(i),...
        'c', [+1, -(model.ub(i)-WT_upperTol(i))], 'dsense', 'L');
end
for i=1:nRxns
    ROOMmodel = addCOBRAConstraints(ROOMmodel, {model.rxns{i}, strcat('binary_',model.rxns{i})}, WT_lowerTol(i),...
        'c', [+1, -(model.lb(i)-WT_lowerTol(i))], 'dsense', 'G');
end
% Block reactions in model
ROOMmodel = changeRxnBounds(ROOMmodel, rxnKO, 0, 'b');
% change objective function
ROOMmodel.c(:) = 0;
ROOMmodel.evarc(cellfun(@length,regexp(ROOMmodel.evars,'^binary_'))==1) = 1;
ROOMmodel.osenseStr = 'min';

solutionROOM = optimizeCbModel(ROOMmodel);

if solutionROOM.stat == 1
    fluxROOM = solutionROOM.v(1:nRxns);
    totalFluxDiff = norm(fluxROOM-fluxWT);
else
    fluxROOM = nan(nRxns,1);
    totalFluxDiff = nan;
end

end

author: opencobra

repo: cobratoolbox


text:
function [ S, itype, stypes ] = cvx_create_structure( varargin )

%CVX_CREATE_STRUCTURE Construct a basis for a structured matrix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parse the arguments, if needed %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S = [];
itype = '';
stypes = {};
if isstruct( varargin{1} ),
    args = varargin{1};
    orig = { args.orig };
    name = { args.name };
    args = { args.args };
    nargs = numel( orig );
else
    nargs = nargin;
    orig = varargin;
    name = orig;
    args = cell( 1, nargs );
    if ischar( varargin{1} ),
        amin = 1;
    else
        args{1} = varargin{1};
        name{1} = '';
        varargin{1} = '';
        amin = 2;
    end
    toks = regexp( varargin, '^([a-zA-Z]\w*)(\(.*\))?$', 'tokens' );
    for k = amin : nargs,
        tok = toks{k};
        if isempty( tok ),
            if k == 1, type = 'variable'; else type = 'structure'; end
            cvx_throw( 'Invalid %s specification: %s', type, varargin{k} );
        end
        tok = tok{1};
        name{k} = tok{1};
        if length(tok) > 1 && ~isempty( tok{2} ),
            try
                args{k} = evalin( 'caller', [ '{', tok{2}(2:end-1), '};' ] );
            catch exc
                cvx_throw( 'Error evaluating structure type: %s\n    %s', orig{k}, exc.message );
            end
        else
            args{k} = {};
        end
    end
end
if nargs == 1 && nargout > 0,
    return
end
sz = args{1};
if iscell( sz ),
    sz = [ sz{:} ];
end
sz(end+1:2) = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scan the structure strings for symmetry and complex cases. We now handle %
% these cases specially for improved performance.                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uplo = '';
do_semi = 0;
do_skew = false;
do_conj = false;
do_comp = false;
do_symm = false;
do_nneg = false;
do_dble = false;
is_toep = false;
do_strx = false;
hflags  = false(1,nargs); n_hflags = 0;
iflags  = false(1,nargs); n_iflags = 0;
sflags  = false(1,nargs); n_sflags = 0;
bflags  = false(1,nargs); n_bflags = 0;
pflags  = false(1,nargs); n_pflags = 0;
switch length(sz),
    case 0, bands = [1,1];
    case 1, bands = sz * [1,1];
    otherwise, bands = sz(1:2);
end
for k = 2 : nargs,
    amin = 0; amax = 0;
    nm = lower( name{k} );
    if ~isempty(uplo),
        nm = [ uplo, '_', nm ]; %#ok
        name{k} = nm;
        orig{k} = [ uplo, '_', orig{k} ];
        uplo = '';
    end
    if nm(end) == '_',
        stypes{end+1} = nm; %#ok
        continue;
    end
    switch nm,
        case { 'binary', 'integer' },
            itype = nm;
            iflags(k) = true;
        case { 'upper', 'lower', 'skew', 'doubly' },
            uplo = nm;
            continue;
        case 'nonnegative',
            stypes{end+1} = nm; %#ok
            do_nneg = true;
        case 'complex',
            do_comp = true;
            if do_semi, do_conj = true; end
            do_strx = true;
        case 'complex_if',
            amin = 1; amax = 1;
            if length(args{k}) == 1,
                arg = args{k}{1};
                if isa(arg,'logical') && arg || ~isreal(arg),
                    do_comp = true;
                    if do_semi, do_conj = true; end
                    do_strx = true;
                end
            end
        case 'symmetric',
            sflags(k) = true;                 do_symm = true;
            do_strx = true;
        case 'symmetric_ut',
            sflags(k) = true;
            pflags(k) = true;
            do_strx = true;
        case 'hermitian',
            sflags(k) = true; do_comp = true; do_symm = true; do_conj = true;
            do_strx = true;
        case 'hermitian_if',
            amin = 1; amax = 1;
            sflags(k) = true; do_symm = true;
            if length(args{k}) == 1,
                arg = args{k}{1};
                if isa(arg,'logical') && arg || ~isreal(arg),
                    do_comp = true; do_conj = true;
                    do_strx = true;
                end
            end
        case { 'skew_symmetric', 'skew-symmetric' },
            sflags(k) = true;                 do_symm = true;                 do_skew = true;
            do_strx = true;
        case { 'skew_hermitian', 'skew-hermitian' }
            sflags(k) = true; do_comp = true; do_symm = true; do_conj = true; do_skew = true;
            do_strx = true;
        case {'hankel','upper_hankel'},
            pflags(k) = true;
            sflags(k) = true;
            hflags(k) = true;
            do_strx = true;
        case 'toeplitz',
            hflags(k) = true;
            is_toep = true;
            do_strx = true;
        case { 'semidefinite', 'doubly_nonnegative' },
            do_semi = true; do_symm = true;
            if do_comp, do_conj = true; end
            stypes{end+1} = 'semidefinite'; %#ok
            if nm(1) == 'd',
                stypes{end+1} = 'nonnegative'; %#ok
                do_nneg = true; do_dble = true;
            end
            do_strx = true;
        case 'upper_bidiagonal',
            bflags(k) = true; sflags(k) = true;
            bands = min(bands,[0,1]);
            do_strx = true;
        case 'lower_bidiagonal',
            bflags(k) = true; sflags(k) = true;
            bands = min(bands,[1,0]);
            do_strx = true;
        case 'upper_hessenberg',
            bflags(k) = true; sflags(k) = true;
            bands = min(bands,[1,Inf]);
            do_strx = true;
        case 'lower_hessenberg',
            bflags(k) = true; sflags(k) = true;
            bands = min(bands,[Inf,1]);
            do_strx = true;
        case 'upper_triangular',
            bflags(k) = true; sflags(k) = true;
            bands = min(bands,[0,Inf]);
            do_strx = true;
        case 'lower_triangular',
            bflags(k) = true; sflags(k) = true;
            bands = min(bands,[Inf,0]);
            do_strx = true;
        case 'diagonal',
            bflags(k) = true;
            bands = min(bands,[0,0]);
            do_strx = true;
        case 'tridiagonal',
            bflags(k) = true;
            bands = min(bands,[1,1]);
            do_strx = true;
        case 'scaled_identity',
            bflags(k) = true; is_toep = true;
            bands = min(bands,[0,0]);
            do_strx = true;
        case 'banded',
            amin = 1; amax = 2;
            sflags(k) = length(args{k}) == 2 && ~isequal(args{k}{1},args{k}{2});
            bflags(k) = true;
            lb = args{k}{1};
            switch length(args{k}),
                case 0, cvx_throw( 'Not enough arguments for "banded()" structure.' );
                case 1, ub = args{k}{1};
                case 2, ub = args{k}{2};
                otherwise, cvx_throw( 'Too many arguments for "banded()" structure.' );
            end
            if ~isnumeric( lb ) || length( lb ) ~= 1 || lb < 0 || lb ~= floor( lb ),
                cvx_throw( 'Bandwidth arguments must be nonnegative integers.' );
            elseif ~isnumeric( ub ) || length( ub ) ~= 1 || ub < 0 || ub ~= floor( ub ),
                cvx_throw( 'Bandwidth arguments must be nonnegative integers.' );
            end
            bands = min(bands,[lb,ub]);
            do_strx = true;
        otherwise,
            pflags(k) = true;
            if ~exist( [ 'cvx_s_', nm ], 'file' ),
                cvx_throw( 'Undefined matrix structure type: %s\nTrying to declare multiple variables? Use the VARIABLES keyword instead.', orig{k} );
            end
    end
    n_bflags = n_bflags + bflags(k);
    n_iflags = n_iflags + iflags(k);
    n_sflags = n_sflags + sflags(k);
    n_hflags = n_hflags + hflags(k);
    n_pflags = n_pflags + pflags(k);
    if length( args{k} ) < amin,
        cvx_throw( 'Not enough arguments: %s', orig{k} );
    elseif length( args{k} ) > amax,
        cvx_throw( 'Too many arguments: %s', orig{k} );
    end
end
if ~isempty(uplo),
    cvx_throw( 'Invalid structure type: %s', orig{end} );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verify integer consistency %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if n_iflags,
    if n_iflags > 1,
        cvx_throw( 'Multiple integer keywords used:\n   %s', sprintf(' %s', orig{iflags} ) );
    elseif do_comp,
        xflags = iflags | findstr( name, 'complex', 'hermitian', 'skew_hermitian' );
        cvx_throw( 'Integer variables must also be real:\n   %s', sprintf(' %s', orig{xflags}) );
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quick exit for no structure %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~do_strx && nargout > 0
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verify structure consistency %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xflags = 0;
if n_sflags > 1, xflags = xflags | sflags; end
if n_hflags > 1, xflags = xflags | hflags; end
if n_iflags > 1, xflags = xflags | iflags; end
if n_bflags > 1, xflags = xflags | bflags; end
if do_semi
    if do_skew
        xflags = xflags | sflags | findstr( name, 'semidefinite', 'doubly_nonnegative' );
    elseif do_nneg && ~do_dble,
        xflags = xflags | findstr( name, 'nonnegative', 'semidefinite' );
    end
end
if do_comp && ( n_iflags || do_nneg )
    xflags = xflags | iflags | findstr( name, 'complex', 'hermitian', 'nonnegative', 'doubly_nonnegative' );
end
if do_nneg && do_skew
    xflags = xflags | findstr( name, 'nonnegative', 'skew_symmetric' );
end
if any( xflags ),
    cvx_throw( 'These forms of structure may not be specified simultaneously:\n   %s', sprintf(' %s', orig{xflags} ) );
end
if do_symm && sz(1) ~= sz(2),
    xflags = sflags | findstr( name, 'semidefinite', 'doubly_nonnegative' );
    cvx_throw( 'This matrix structure requires square matrices:%s', sprintf(' %s', orig{xflags} ) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the basis matrices for the remaining structure elements        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

strs = {};
if is_toep || any( bands < sz(1:2) ),
    [ strs{end+1}, do_symm ] = cvx_s_banded( sz(1), sz(2), [ do_symm, is_toep ], bands(1), bands(2) );
    do_symm = do_symm(1);
end
if n_pflags,
    for k = find(pflags)',
        try
            [ strs{end+1}, do_symm ] = feval( [ 'cvx_s_', lower(name{k}) ], sz( 1 ), sz( 2 ), do_symm, args{k}{:} ); %#ok
        catch exc
            cvx_throw( 'Error constructing structure: %s\n   %s', orig{k}, exc.message );
        end
    end
end
if do_symm,
    [ strs{end+1}, do_symm ] = cvx_s_symmetric( sz(1), sz(2), do_symm ); %#ok
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If multiple structures have been requested (e.g., toeplitz and banded),  %
% combine them together by finding bases for their orthogonal complements, %
% concatenating, and taking the orthogonal complement of that. This should %
% be used much less frequently than before---if ever---now that we handle  %
% symmetry as a special case for improved performance.                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch length( strs ),
    case 0,
        if ~do_comp && ~do_skew && nargout ~= 0, return; end
        sz(end+1:2) = 1;
        nel = sz( 1 ) * sz( 2 );
        S = sparse( 1 : nel, 1 : nel, 1, nel, nel );
    case 1,
        S = strs{ 1 };
    otherwise,
        for k = 1 : length(strs),
            strs{k} = cvx_orthog_structure( strs{k} ); %#ok
        end
        S = cvx_orthog_structure( vertcat(strs{:}), true );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Handle complex, skew-symmetric, and Hermitian structures.                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if do_comp,
    S = [ S ; +1j * S ];
    S = S([1:end/2;end/2+1:end],:);
end
if do_skew || do_conj,
    r = (0:sz(1)-1)'; r = r(:,ones(1,sz(2)));
    c = 0:sz(2)-1; c = c(ones(1,sz(1)),:);
    ut = r < c; dg = r == c;
    if do_skew,
        S(:,ut) = - S(:,ut);
        S(:,dg) = 0;
    end
    if do_conj,
        S(:,ut) = conj(S(:,ut));
        S(:,dg) = real(S(:,dg));
    end
    S = S(any(S,2),:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Report an error of the structure is empty                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty( S ),
    cvx_throw( 'Incompatible structure modifiers:%s', sprintf( ' %s', args.orig ) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replicate structure for N-D arrays                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if length( sz ) > 2,
    S = cvx_replicate_structure( S, sz( 3 : end ) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the structure if called with no output arguments                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout == 0,
    [ii,jj,vv] = find( S );
    if isempty(sz), sz = [1,1]; end
    Z = reshape( full( sparse( jj, 1, ii .* vv, prod(sz), 1 ) ), sz );
    temp = sprintf( ',%d', sz );
    fprintf( 'Name: %s\nSize: %s\nStructure:', name{1}, temp(2:end) );
    if ~isempty( itype ),
        stypes = [ stypes, itype ];
    end
    for k = 2 : nargs,
        if ~any( strcmp(orig{k},stypes) )
            fprintf( ' %s', orig{k} );
        end
    end
    fprintf( '\n' );
    if ~isempty( stypes ),
        fprintf( 'Passthrough:' );
        for k = 1 : length(stypes),
            fprintf( ' %s', stypes{k} );
        end
    end
    fprintf( '\n\n' );
    fmt = get(0,'format');
    set(0,'format','rational');
    disp( Z );
    set(0,'format',fmt);
    clear S
end

function z = findstr( x, varargin )
q = varargin;
z = cellfun( @(z)any(strcmpi(z,q)), x );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for orig 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:
function [means,Nmeans] = simple_kmeans(X,K,maxerr)

% function [medias,Nmedias] = simple_kmedias(X,K,maxerr)
%   Finds K prototypes representing the samples in data matrix X,
%   where each row of X represents a sample.
%   Iterates until maximum norm difference between
%   prototypes found in successive iterations is < maxerr
%
%   This script uses square Euclidean distance,
%   but can be easily modified to use other metrics
%
% Output arguments
%   means: matrix with each row a cluster prototype
%   Nmeans: Number of samples in each cluster
%
% Example:
%   X = [randn(100,1) ;  2+randn(100,1)];
%   K = 2;
%   [means Nmeans] = simple_kmeans(X,K,0)
%
%   Mauricio Martinez-Garcia, 2003,2007

[Ndata, dims] = size(X);
dist = zeros(1,K);

% Initial prototype assignment (arbitrary)
for i=1:K-1
   means(i,:) = X(i,:);
end
means(K,:) = mean(X(K:Ndata,:));

cmp = 1 + maxerr;
while (cmp > maxerr)
   % Sums (class) and data counters (Nclass) initialization
   class = zeros(K,dims);
   Nclass = zeros(K,1);

   % Groups each elements to the nearest prototype
   for i=1:Ndata
      for j=1:K
         % Euclidean distance from data to each prototype
         dist(j) = norm(X(i,:)-means(j,:))^2;
      end
      % Find indices of minimum distance
      index_min = find(~(dist-min(dist)));
      % If there are multiple min distances, decide randomly
      index_min = index_min(ceil(length(index_min)*rand));
      class(index_min,:) = class(index_min,:) + X(i,:);
      Nclass(index_min) = Nclass(index_min) + 1;
   end
   for i=1:K
      class(i,:) = class(i,:) / Nclass(i);
   end

   % Compare results with previous iteration
   cmp = 0;
   for i=1:K
      cmp = norm(class(i,:)-means(i,:));
   end

   % Prototype update
   means = class;
end

Nmeans = Nclass;


text:
function prob = my_mvnpdf(x,meanV,varV,diagCov)
if nargin < 3
    diagCov = 0;
end
if diagCov == 0
    prob = mvnpdf(x,meanV,varV);
else
    % meanV, varV are row vectors
    % each row of x is a sample
    [N_vector,Dim] = size(x);
    norm_term = ( (2*pi)^(Dim/2) ) * prod( sqrt(varV) );
    if 1
        tmp = bsxfun(@minus, x, meanV);
        tmp = tmp.*tmp;
        tmp = bsxfun(@times, tmp, 1./varV);
    end
    if 0
        tmp = (x - repmat(meanV,N_vector,1)).^2;
        tmp = tmp ./ repmat(varV,N_vector,1);
    end
    tmp = sum(tmp,2);
    prob = exp(-tmp/2) / norm_term;
end

end

author: singaxiong

repo: SignalGraph


text:
function this = realign(this)
% Realigns all 3D images in 4D data to each other, then to the mean
% Uses SPM's realign: estimate+rewrite functionality
%
%   MrSeries = realign(MrSeries)
%
% This is a method of class MrSeries.
%
% IN
%   MrImage.data
%   MrImage.parameters.realign.quality
%
% OUT
%
% EXAMPLE
%   realign
%
%   See also MrSeries MrImage MrImage.realign

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-01
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

quality = this.parameters.realign.quality;

this.init_processing_step('realign');
this.data = this.data.realign('quality', quality);
this.finish_processing_step('realign', this.data);

author: translationalneuromodeling

repo: tapas


text:
function NNF = sr_pass(imgTrg, imgSrcPyr, NNF, modelPlane, numIterLvl, opt)
%
% SR_PASS
%
% Nearest neighbor field estimation at the current level using the
% generalized PatchMatch algorithm
%
% Input:
%   - imgTrg:     target image
%   - imgSrcPyr:  source image pyramid
%   - NNF:        the current nearest neighbor field
%   - modelPlane: planar structure model
%   - numIterLvl: number of iterations at the current level
%   - opt:        parameters
% Output:
%   - NNF:        updated nearest neighbor field
% =========================================================================

% =========================================================================
% Compute the initial patch cost at the current level
% =========================================================================
% Prepare target patch
trgPatch = sr_prep_target_patch(imgTrg, opt.pSize);

% Prepare source patch
[srcPatch, srcPatchScale] = sr_prep_source_patch(imgSrcPyr, NNF.uvTformH.data, opt);
% Compute patch matching cost: appearance cost
[NNF.uvCost.data, NNF.uvBias.data] = sr_patch_cost_app(trgPatch, srcPatch, opt);
% Compute patch matching cost: scale cost
if(opt.useScaleCost)
    costScale = opt.lambdaScale*max(0, opt.scaleThres - srcPatchScale);
    NNF.uvCost.data = NNF.uvCost.data + costScale;
end
% Compute patch matching cost: plane compatibility cost
if(opt.usePlaneGuide)
    costPlane = sr_patch_cost_plane(modelPlane.mLogLPlaneProb, NNF.uvPlaneID.data, NNF.uvPix.ind, NNF.uvTformH.data(:,7:8));
    NNF.uvCost.data = NNF.uvCost.data + opt.lambdaPlane*costPlane;
end

% Update cost map
NNF.uvCost.map = sr_update_uvMap(NNF.uvCost.map, NNF.uvCost.data, NNF.uvPix.ind);

% Initialize update index map (for early termination)
% NNF.update.data = true(NNF.uvPix.numUvPix, 1);
% NNF.update.map  = true(NNF.imgH, NNF.imgW);

% Initialize visualization
if(opt.visFlag)
    uvTformMap    = zeros(NNF.imgH, NNF.imgW, 3, numIterLvl);
    uvCostMap     = zeros(NNF.imgH, NNF.imgW, numIterLvl);
    uvCostMapInit = NNF.uvCost.map;
    imgRecSet     = zeros(NNF.imgH, NNF.imgW, 3, numIterLvl);
end

% =========================================================================
% Update the nearest neighbor field using PatchMatch
% =========================================================================

for iter = 1 : numIterLvl
    [NNF, nUpdate] = sr_update_NNF(trgPatch, imgSrcPyr, NNF, modelPlane, opt);
    avgPatchCost = mean(NNF.uvCost.data);

    fprintf('    %3d\t%12d\t%12d\t%14f\n', iter, nUpdate(1), nUpdate(2), avgPatchCost);

    if(0)
        NNFVis = sr_vis_nnf(NNF, imgSrcPyr, opt);

%         uvTformMap(:,:,:, iter) = NNFVis.uvTfomMapVis;
%         uvCostMap(:,:, iter) = NNF.uvCost.map;
%         imgRecSet(:,:,:, iter) = imgRec;

%         figure(1);
%         subplot(2,4,1), imshow(imresize(imgSrcPyr{opt.origResLvl}, [NNF.imgH, NNF.imgW], 'bicubic'));
%         title('Current Image', 'fontsize', 16);
%         subplot(2,4,2), imshow(imgRec);
%         title('Super-resolution Image', 'fontsize', 16);
%         subplot(2,4,3), imshow(NNFVis.uvTformScaleVis);
%         title('Scale map', 'fontsize', 16);
%         subplot(2,4,4), imshow(NNFVis.uvTformRotVis); colormap jet
%         title('Rotation map', 'fontsize', 16);
%         subplot(2,4,5), imshow(NNFVis.uvTfomMapVis);
%         title('Nearest neighbor field', 'fontsize', 16);
%         subplot(2,4,6), imshow(NNFVis.uvPosMap); colormap jet
%         title('Position map', 'fontsize', 16);
%         subplot(2,4,7), imshow(NNFVis.uvCostMapVis); colormap jet
%         title('Cost map', 'fontsize', 16);
%         subplot(2,4,8), imshow(NNFVis.uvPlaneIDMapVis); colormap jet
%         title('Plane ID', 'fontsize', 16);
    end
end

% Visualization
if(opt.visFlag)
    resPath = 'result/nIter';
    if(~exist(resPath, 'dir'))
        mkdir(resPath);
    end
    imgID = opt.imgID;
    maxCost = max(uvCostMapInit(:));
    uvCostMap = uvCostMap/maxCost;

    uvCostMapInit = uvCostMapInit/maxCost;
    h2 = figure(2);
    imshow(uvCostMapInit); colormap jet;
    uvCostMapC = getframe(h2);
    imwrite(uvCostMapC.cdata,  ...
        fullfile(resPath, ['img_', num2str(imgID, '%03d'), '_costMap_iter_00.png']));

    % Position map
    uvPosMap = NNFVis.uvPosMap;

    imwrite(uvPosMap, fullfile(resPath, ...
        ['img_', num2str(imgID, '%03d'), '_uvPosMap.png']));

    for iter = 1: numIterLvl
        % Save uvTformMap
        imwrite(uvTformMap(:,:,:,iter), ...
            fullfile(resPath, ['img_', num2str(imgID, '%03d'),'_tformMap_iter_', num2str(iter, '%02d'), '.png']));

        % Save reconstructed image
        imwrite(imgRecSet(:,:,:,iter),  ...
            fullfile(resPath, ['img_', num2str(imgID, '%03d'),'_rec_iter_', num2str(iter, '%02d'), '.png']));

        % Save matching cost map
        h2 = figure(2);
        imshow(uvCostMap(:,:,iter)); colormap jet;
        uvCostMapC = getframe(h2);
        imwrite(uvCostMapC.cdata,  ...
            fullfile(resPath, ['img_', num2str(imgID, '%03d'), '_costMap_iter_', num2str(iter, '%02d'), '.png']));
    end
end

end

repo: SelfExSR