text:
% Run file to produce the CME experiment from the paper
%   S. Dolgov, D. Savostyanov,
%   "Alternating minimal energy methods for linear systems in higher
%   dimensions"

try
    maxNumCompThreads(1);
catch
    % Just skip. Sometimes if you specify -singleCompThread in the command
    % line, MATLAB will fail at maxNumCompThreads with scary, so tell him
    % it's okay.
end;

n = 64*ones(20,1);
T = 10;
Lt = 12;
tol = 1e-6;
kickrank = 4;
nswp = 20;
trunc_norm = 'fro';
verb = 3;

tol_gmres = 1e-3; % I would not set anything lower

d = numel(n);

Ap = cell(d,1);
% Cascadic MPO -- production reactions
Ap{1} = zeros(1,n(1),n(1),3);
Ap{1}(1,:,:,1) = eye(n(1));
Ap{1}(1,:,:,2) = diag((0:n(1)-1)./(5+(0:n(1)-1)));
Ap{1}(1,:,:,3) = (diag(ones(n(1)-1,1),-1)-eye(n(1)))*0.7*diag([ones(n(1)-1,1);0]);
for i=2:d-1
    Ap{i} = zeros(3,n(i),n(i),3);
    Ap{i}(1,:,:,1)=eye(n(i));
    Ap{i}(3,:,:,3)=eye(n(i));
    Ap{i}(1,:,:,2)=diag((0:n(i)-1)./(5+(0:n(i)-1)));
    Ap{i}(2,:,:,3)=(diag(ones(n(i)-1,1),-1)-eye(n(i)))*diag([ones(n(i)-1,1);0]);
end;
Ap{d} = zeros(3,n(d),n(d));
Ap{d}(2,:,:) = (diag(ones(n(d)-1,1),-1)-eye(n(d)))*diag([ones(n(d)-1,1);0]);
Ap{d}(3,:,:) = eye(n(d));
% Laplace MPO -- destruction reactions
Ad = cell(d,1);
Ad{1} = zeros(1,n(1),n(1),2);
Ad{1}(1,:,:,1) = (diag(ones(n(1)-1,1),1)-eye(n(1)))*diag((0:n(1)-1))*0.07;
Ad{1}(1,:,:,2) = eye(n(1));
for i=2:d-1
    Ad{i} = zeros(2,n(i),n(i),2);
    Ad{i}(1,:,:,1)=eye(n(i));
    Ad{i}(2,:,:,2)=eye(n(i));
    Ad{i}(2,:,:,1)=(diag(ones(n(i)-1,1),1)-eye(n(i)))*diag((0:n(i)-1))*0.07;
end;
Ad{d} = zeros(2,n(d),n(d),1);
Ad{d}(1,:,:) = eye(n(d));
Ad{d}(2,:,:) = (diag(ones(n(d)-1,1),1)-eye(n(d)))*diag((0:n(d)-1))*0.07;

% QTT-tize each dimension separately, it is more stable due to smaller size
Apq = [];
Adq = [];
for i=1:d
    r1 = size(Ap{i},1);
    r2 = size(Ap{i},4);
    Aloc = cell2core(tt_matrix, Ap(i));
    Aloc = tt_reshape(Aloc, factor(n(i))'*[1,1], 1e-13, r1, r2);
    Apq = tkron(Apq, Aloc);

    r1 = size(Ad{i},1);
    r2 = size(Ad{i},4);
    Aloc = cell2core(tt_matrix, Ad(i));
    Aloc = tt_reshape(Aloc, factor(n(i))'*[1,1], 1e-13, r1, r2);
    Adq = tkron(Adq, Aloc);
end;

% Now the final CME operator
A = Apq+Adq;
I = tt_eye(A.n);

% Global time scheme
tau = T/2^Lt;
Gt = IpaS(Lt,-1)/tau;
iGt = tt_qtoepl(tkron(tt_ones(2,Lt), tt_tensor([0;1])), Lt)*tau;
Mt = tt_eye(2,Lt);
e1t = tt_unit(2,Lt,1);
et = tt_ones(2,Lt);
eNt = tt_unit(2,Lt,2^Lt);

% Global matrix
B = tkron(I,tt_eye(2,Lt)) - tkron(A,iGt*Mt);

% Initial state
u = [];
for i=1:d
    u = tkron(u, tt_unit(factor(n(i))', numel(factor(n(i))), 1)); %QTT
end;
U0 = tkron(u,et);

% RHS
f = tkron(u, et);

% Symmetrized matrix and RHS
B2 = B'*B;
B2 = round(B2, 1e-13);
f2 = B'*f;
f2 = round(f2, 1e-13);

% Use MEX library. It is safe, since in the last TT-Toolbox it switches to
% pure Matlab automatically, if you did not compile the MEXes.
ismex = true;

% Solvers
% Reference solution
tic;
U_ex = amen_solve2(B, f, tol*1e-3, 'x0', U0, 'nswp', nswp*3, 'kickrank', kickrank, 'trunc_norm', trunc_norm, 'ismex', true, 'verb', 1, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2);
toc;

% ALS(SD)
tic;
[U_alstpz,td_alstpz] = alstpz_solve(B,f,tol, 'x0',U0, 'max_full_size', 50, 'kickrank', kickrank, 'nswp', nswp, 'kicktype', 'svd', 'symm', false, 'ismex', ismex, 'verb', verb);
toc;
% ALS(SD)+SYMM
tic;
[U_alstpz_s,td_alstpz_s] = alstpz_solve(B2,f2,tol, 'x0',U0, 'max_full_size', 50, 'kickrank', kickrank, 'nswp', nswp, 'kicktype', 'svd', 'symm', false, 'ismex', ismex, 'verb', verb);
toc;
% AMEN+SVD
tic;
[U_amen_svd,td_amen_svd] = amen_solve2(B, f, tol, 'x0', U0, 'nswp', nswp, 'kicktype', 'svd', 'kickrank', kickrank, 'trunc_norm', trunc_norm, 'ismex', ismex, 'verb', verb, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2);
toc;
% AMEN+SVD+SYMM
tic;
[U_amen_svd_s,td_amen_svd_s] = amen_solve2(B2, f2, tol, 'symm', false, 'x0', U0, 'nswp', nswp, 'kicktype', 'svd', 'kickrank', kickrank, 'trunc_norm', trunc_norm, 'ismex', ismex, 'verb', verb, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2);
toc;
% AMEN+ALS
tic;
[U_amen_als,td_amen_als] = amen_solve2(B, f, tol, 'x0', U0, 'nswp', nswp, 'kickrank', kickrank, 'trunc_norm', trunc_norm, 'ismex', ismex, 'verb', verb, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2);
toc;
% AMEN+ALS+SYMM
tic;
[U_amen_als_s,td_amen_als_s] = amen_solve2(B2, f2, tol, 'symm', false, 'x0', U0, 'nswp', nswp, 'kickrank', kickrank, 'trunc_norm', trunc_norm, 'ismex', ismex, 'verb', verb, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2);
toc;
% DMRG
tic;
[U_dmrg,td_dmrg] = dmrg_solve3(B, f, tol, 'x0', U0, 'nswp', nswp, 'kickrank', 0, 'trunc_norm', trunc_norm, 'ismex', ismex, 'verb', verb, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2, 'step_drank', 0, 'step_dpow', 0, 'min_dpow', 0.5, 'dirfilter', 1);
toc;
% DMRG+SYMM
tic;
[U_dmrg_s,td_dmrg_s] = dmrg_solve3(B2, f2, tol, 'symm', false, 'x0', U0, 'nswp', nswp, 'kickrank', 0, 'trunc_norm', trunc_norm, 'ismex', ismex, 'verb', verb, 'max_full_size', 50, 'local_restart', 50, 'local_iters', 2, 'resid_damp', 2, 'step_drank', 0, 'step_dpow', 0, 'min_dpow', 0.5, 'dirfilter', 1);
toc;

% TT-GMRES
[U_gmres,td_gmres] = tt_gmres(core(B), core(f), tol_gmres, 10, 15, tol_gmres, tol_gmres, [], [], [], [], verb);

% This is a time-dep problem => check the last snapshots
if (~isempty(whos('U_ex')))
    u_ex = chunk(U_ex, 1, f.d-Lt)*dot(eNt, chunk(U_ex, f.d-Lt+1, f.d)).';
    u_ex = tt_reshape(u_ex, n);
end;
fprintf('Snapshot errors: \n');
if (~isempty(whos('U_alstpz')))
    u_alstpz = chunk(U_alstpz, 1, f.d-Lt)*dot(eNt, chunk(U_alstpz, f.d-Lt+1, f.d)).';
    u_alstpz = tt_reshape(u_alstpz, n);
    fprintf('alstpz:\t\t%3.5e\n', norm(u_alstpz-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_alstpz_s')))
    u_alstpz_s = chunk(U_alstpz_s, 1, f.d-Lt)*dot(eNt, chunk(U_alstpz_s, f.d-Lt+1, f.d)).';
    u_alstpz_s = tt_reshape(u_alstpz_s, n);
    fprintf('alstpz_s:\t\t%3.5e\n', norm(u_alstpz_s-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_amen_svd')))
    u_amen_svd = chunk(U_amen_svd, 1, f.d-Lt)*dot(eNt, chunk(U_amen_svd, f.d-Lt+1, f.d)).';
    u_amen_svd = tt_reshape(u_amen_svd, n);
    fprintf('amen_svd:\t\t%3.5e\n', norm(u_amen_svd-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_amen_svd_s')))
    u_amen_svd_s = chunk(U_amen_svd_s, 1, f.d-Lt)*dot(eNt, chunk(U_amen_svd_s, f.d-Lt+1, f.d)).';
    u_amen_svd_s = tt_reshape(u_amen_svd_s, n);
    fprintf('amen_svd_s:\t\t%3.5e\n', norm(u_amen_svd_s-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_amen_als')))
    u_amen_als = chunk(U_amen_als, 1, f.d-Lt)*dot(eNt, chunk(U_amen_als, f.d-Lt+1, f.d)).';
    u_amen_als = tt_reshape(u_amen_als, n);
    fprintf('amen_als:\t\t%3.5e\n', norm(u_amen_als-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_amen_als_s')))
    u_amen_als_s = chunk(U_amen_als_s, 1, f.d-Lt)*dot(eNt, chunk(U_amen_als_s, f.d-Lt+1, f.d)).';
    u_amen_als_s = tt_reshape(u_amen_als_s, n);
    fprintf('amen_als_s:\t\t%3.5e\n', norm(u_amen_als_s-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_dmrg')))
    u_dmrg = chunk(U_dmrg, 1, f.d-Lt)*dot(eNt, chunk(U_dmrg, f.d-Lt+1, f.d)).';
    u_dmrg = tt_reshape(u_dmrg, n);
    fprintf('dmrg:\t\t%3.5e\n', norm(u_dmrg-u_ex)/norm(u_ex));
end;
if (~isempty(whos('U_dmrg_s')))
    u_dmrg_s = chunk(U_dmrg_s, 1, f.d-Lt)*dot(eNt, chunk(U_dmrg_s, f.d-Lt+1, f.d)).';
    u_dmrg_s = tt_reshape(u_dmrg_s, n);
    fprintf('dmrg_s:\t\t%3.5e\n', norm(u_dmrg_s-u_ex)/norm(u_ex));
end;

% Compare with KSL
Nksl = 50;
% Set the ranks to the proper valus returned by AMEN
u_ksl = chunk(U_amen_als, 1, f.d-Lt)*dot(eNt, chunk(U_amen_als, f.d-Lt+1, f.d)).';
u_ksl = round(u_ksl, tol*0.1);
r_ksl = u_ksl.r;
tic;
u_ksl = u+0*tt_rand(u.n,u.d,r_ksl,-1);
for j=1:Nksl
    u_ksl = tt_ksl_ml(u_ksl, A, tt_zeros(u_ksl.n), T/Nksl);
end;
ttimes_ksl=toc;
u_ksl = tt_reshape(u_ksl, n);
fprintf('ksl:\t\t%3.5e\n', norm(u_ksl-u_ex)/norm(u_ex));
fprintf('CPU Time of the KSL: %g\n', ttimes_ksl);

% % Result history processing
err_alstpz = zeros(nswp,1);
err_alstpz_s = zeros(nswp,1);
err_amen_svd = zeros(nswp,1);
err_amen_svd_s = zeros(nswp,1);
err_amen_als = zeros(nswp,1);
err_amen_als_s = zeros(nswp,1);
err_dmrg = zeros(nswp,1);
err_dmrg_s = zeros(nswp,1);
% Measure the errors
for i=1:nswp
    if (~isempty(whos('td_alstpz')))
        err_alstpz(i) = norm(td_alstpz{2}{end,i}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_alstpz_s')))
        err_alstpz_s(i) = norm(td_alstpz_s{2}{end,i}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_amen_svd')))
        err_amen_svd(i) = norm(td_amen_svd{2}{end,i}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_amen_svd_s')))
        err_amen_svd_s(i) = norm(td_amen_svd_s{2}{end,i}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_amen_als')))
        err_amen_als(i) = norm(td_amen_als{2}{end,i}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_amen_als_s')))
        err_amen_als_s(i) = norm(td_amen_als_s{2}{end,i}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_dmrg')))
        err_dmrg(i) = norm(td_dmrg{2}{end-1,i*2-1}-U_ex)/norm(U_ex);
    end;
    if (~isempty(whos('td_dmrg_s')))
        err_dmrg_s(i) = norm(td_dmrg_s{2}{end-1,i*2-1}-U_ex)/norm(U_ex);
    end;
end;

% Prepare the data in the TikZ-readable form
dats = [(1:nswp)', td_alstpz{1}(end,:)', err_alstpz, td_alstpz_s{1}(end,:)', err_alstpz_s, ...
    td_amen_svd{1}(end,:)', err_amen_svd, td_amen_svd_s{1}(end,:)', err_amen_svd_s, ...
    td_amen_als{1}(end,:)', err_amen_als, td_amen_als_s{1}(end,:)', err_amen_als_s, ...
    td_dmrg{1}(end-1,1:2:end)', err_dmrg, td_dmrg_s{1}(end-1,1:2:end)', err_dmrg_s];
% dats layout:
%   iter(1)     t_alstpz(2)     e_alstpz(3)     t_alstpz_s(4)       e_alstpz_s(5)
%   t_as(6)     e_as(7)         t_as_s(8)       e_as_s(9)
%   t_aa(10)    e_aa(11)        t_aa_s(12)      e_aa_s(13)
%   t_d(14)     e_d(15)         t_d_s(16)       e_d_s(17)

% Process the GMRES output. It is different...
if (~isempty(whos('td_gmres')))
    iter_gmres = min([find(td_gmres{1}==0, 1), numel(td_gmres{1})+1])-1;
    err_gmres = zeros(iter_gmres, 1);
    for i=1:iter_gmres
        err_gmres(i) = norm(tt_tensor(td_gmres{2}{i})-U_ex)/norm(U_ex);
    end;
    dat_gmres = [(1:iter_gmres)', td_gmres{1}(1:iter_gmres)', err_gmres];
end;

% Draw 'em for humans
% iter
figure(1);
semilogy(dats(:,1), dats(:,[3,5,7,9,11,13,15,17]), dat_gmres(1:min(iter_gmres,nswp),1), dat_gmres(1:min(iter_gmres,nswp),3));
legend('alstz', 'alstz-s', 'amen-svd', 'amen-svd-s', 'amen-als', 'amen-als-s', 'dmrg', 'dmrg-s', 'gmres');
% time
figure(2);
loglog(dats(:,2),dats(:,3), dats(:,4),dats(:,5), ...
    dats(:,6),dats(:,7), dats(:,8),dats(:,9), ...
    dats(:,10),dats(:,11), dats(:,12),dats(:,13), ...
    dats(:,14),dats(:,15), dats(:,16),dats(:,17), ...
    dat_gmres(:,2), dat_gmres(:,3));

% % Uncomment this if you want to draw the data elsewhere
% save('cme20_err.dat', '-ascii', 'dats');
% save('cme20_err_gmres.dat', '-ascii', 'dat_gmres');

author: oseledets

repo: TT-Toolbox


text:
%% DEMO_febio_0054_lattice_hydrostatic_01
% Below is a demonstration for:
%
% * Building geometry for a cube with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement and stress results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * uniaxial loading
% * compression, tension, compressive, tensile
% * displacement control, displacement boundary condition
% * hexahedral elements, hex8
% * cube, box, rectangular
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=1;
edgeColor=0.25*ones(1,3);
edgeWidth=1.5;
markerSize=25;
cMap=gjet(4);

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress
febioLogFileName_stiffness=[febioFebFileNamePart,'_stiffness_out.txt']; %Log file name for exporting stiffness

%Specifying dimensions and number of elements
sampleSize=10;
latticeType=1;

%Define applied displacement
J_final=0.7; %Final Jacobian or volume ration
lambdaFinal=J_final^(1/3); %Stretch values in all directions
displacementMagnitude=((lambdaFinal*sampleSize)-sampleSize)/2; %The displacement magnitude

%Material parameter set
c1=1; %Shear-modulus-like parameter
m1=2;
k=50*c1;

% FEA control settings
numTimeSteps=50; %Number of time steps desired
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=25; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=(1/numTimeSteps); %Maximum time step size
min_residual=1e-20;
symmetric_stiffness=0;
runMode='external';

%%

%Specifying dimensions and number of elements
r=0.5; %Radii, results in a width of 1
n=4;
nCopies=n*ones(1,3); %Number of offset copies
d=2*r; %Diameter
w=(n-1)*d; %sampleSize

%% Create lattice

switch latticeType
    case 1 %Octet truss
        [Er,Vr,Cr,Fr,CFr]=rhombicDodecahedronMesh(r,nCopies);
        Vr=Vr./(n-1);
        Vr=Vr*sampleSize;

        [indBoundary]=tesBoundary(Fr,Vr);
        cPar.shrinkFactor=0.15; %Strut sides are formed by shrinking the input mesh faces by this factor
        cPar.meshType='hex'; %desired output mesh type
        cPar.indBoundary=indBoundary; %indices of the boundary faces
        cPar.hexSplit=2;
        cPar.latticeSide=2; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
        [E,V,C]=element2lattice(Er,Vr,cPar); %Get lattice structure

        logicKeep1=~(V(:,1)<=-1e-3);
        logicKeep2=~(V(:,2)<=-1e-3);
        logicKeep3=~(V(:,3)<=-1e-3);
        logicKeep4=~(V(:,1)>=sampleSize+1e-3);
        logicKeep5=~(V(:,2)>=sampleSize+1e-3);
        logicKeep6=~(V(:,3)>=sampleSize+1e-3);

        logicKeepEs=sum(logicKeep1(E),2)>=4 &...
            sum(logicKeep2(E),2)>=4 &...
            sum(logicKeep3(E),2)>=4 &...
            sum(logicKeep4(E),2)>=4 &...
            sum(logicKeep5(E),2)>=4 &...
            sum(logicKeep6(E),2)>=4;

        E=E(logicKeepEs,:);
        C=C(logicKeepEs,:);
        [E,V,indFix]=patchCleanUnused(E,V);

        % [Es,Vs,~,~]=subHex(Es,Vs,1,1);
        % Cs=repmat(Cs,8,1);

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(E,C); %Patch data for plotting

        %Get new boundary set
        indB=tesBoundary(Fs,V);
        Fb=Fs(indB,:);
    case 2 %Rhombic dodecahedron mesh ("dual" of octet truss lattice)
        [Er,Vr,Cr,Fr,CFr]=rhombicDodecahedronMesh(r,nCopies);
        Vr=Vr./(n-1);
        Vr=Vr*sampleSize;

        [indBoundary]=tesBoundary(Fr,Vr);
        cPar.shrinkFactor=0.15; %Strut sides are formed by shrinking the input mesh faces by this factor
        cPar.meshType='hex'; %desired output mesh type
        cPar.indBoundary=indBoundary; %indices of the boundary faces
        cPar.hexSplit=2;
        cPar.latticeSide=1; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
        [E,V,C]=element2lattice(Er,Vr,cPar); %Get lattice structure

        logicKeep1=~(V(:,1)<=-1e-3);
        logicKeep2=~(V(:,2)<=-1e-3);
        logicKeep3=~(V(:,3)<=-1e-3);
        logicKeep4=~(V(:,1)>=sampleSize+1e-3);
        logicKeep5=~(V(:,2)>=sampleSize+1e-3);
        logicKeep6=~(V(:,3)>=sampleSize+1e-3);

        logicKeepEs=sum(logicKeep1(E),2)>=4 &...
            sum(logicKeep2(E),2)>=4 &...
            sum(logicKeep3(E),2)>=4 &...
            sum(logicKeep4(E),2)>=4 &...
            sum(logicKeep5(E),2)>=4 &...
            sum(logicKeep6(E),2)>=4;

        E=E(logicKeepEs,:);
        C=C(logicKeepEs,:);
        [E,V,indFix]=patchCleanUnused(E,V);

        % [Es,Vs,~,~]=subHex(Es,Vs,1,1);
        % Cs=repmat(Cs,8,1);

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(E,C); %Patch data for plotting

        %Get new boundary set
        indB=tesBoundary(Fs,V);
        Fb=Fs(indB,:);
    case 3
        boxDim=sampleSize*[1 1 1];
        boxEl=[2 2 2];
        [meshStruct]=hexMeshBox(boxDim,boxEl);
        Er=meshStruct.E;
        Vr=meshStruct.V;
        minV=min(Vr,[],1); %Get lower left front corner
        Vr=Vr-minV(ones(size(Vr,1),1),:); %Set corner as origin
        [Er,Vr,~]=hex2tet(Er,Vr,[],1); %Convert to tetrahedral elements
        [Fr,Cr]=element2patch(Er,[]); %Patch data for plotting
        [indBoundary]=tesBoundary(Fr,Vr);

        % Create lattice structure
        controlParameter.latticeSide=1;
        controlParameter.numDigitKeep=5; %used for merging nodes
        controlParameter.indBoundary=indBoundary; %indices of the boundary faces
        controlParameter.shrinkFactor=0.15;
        controlParameter.meshType='hex';
        controlParameter.hexSplit=2;

        [E,V,C]=element2lattice(Er,Vr,controlParameter); %Get lattice structure

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(E,C); %Patch data for plotting

        indB=tesBoundary(Fs,V);
        Fb=Fs(indB,:);
end
%%
% Visualizing input mesh and lattic structures

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(Fr,Vr,0.5*ones(1,3),'k',0.5);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fb,V,'bw','k',1);
% patchNormPlot(Fs,Vs);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%% DEFINE BC's

% Define node set logics
indAll=(1:1:size(V,1))';
logicBoundary=ismember(indAll,Fb);

Z=V(:,3);
logicTop=Z>=(sampleSize-eps(sampleSize))& logicBoundary;
logicBottom=Z<=eps(sampleSize) & logicBoundary;

X=V(:,1);
logicSide1=X>=(sampleSize-eps(sampleSize))& logicBoundary;
logicSide2=X<=eps(sampleSize)& logicBoundary;

Y=V(:,2);
logicSide3=Y>=(sampleSize-eps(sampleSize))& logicBoundary;
logicSide4=Y<=eps(sampleSize)& logicBoundary;

%Prescribed force nodes
bcPrescribeListCell{1}=find(logicSide1)';
bcPrescribeListCell{2}=find(logicSide2)';
bcPrescribeListCell{3}=find(logicSide3)';
bcPrescribeListCell{4}=find(logicSide4)';
bcPrescribeListCell{5}=find(logicTop)';
bcPrescribeListCell{6}=find(logicBottom)';

%% Smoothing lattice

% indKeep=unique([bcPrescribeListCell{:}]);
% [Fb_clean,Vb_clean,indFix]=patchCleanUnused(Fb,Vs);
%
% cPar.Method='HC';
% cPar.n=6;
%
% cPar.RigidConstraints=indFix(indKeep);
% % cPar.RigidConstraints=cPar.RigidConstraints(cPar.RigidConstraints>0);
%
% [Vb_clean]=tesSmooth(Fb_clean,Vb_clean,[],cPar);
% ind=Fb(:);
% ind=unique(ind(:));
% Vs(ind,:)=Vb_clean;

% cFigure; hold on;
% gpatch(Fb,Vs,'bw','k',1);
% % patchNormPlot(Fs,Vs);
% % plotV(Vs(indKeep,:),'k.','MarkerSize',25)
% axisGeom(gca,fontSize);
% camlight headlight; lighting flat;
% drawnow;

%%
% Visualizing input mesh and lattic structures

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(Fr,Vr,0.5*ones(1,3),'k',0.5);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fb,V,'bw');
% patchNormPlot(Fs,Vs);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%%
% Visualize BC's

cFigure; hold on;
title('Boundary conditions','FontSize',fontSize);
gpatch(Fb,V,'kw','none',0.4);
hl=gobjects(1,6);
plotColors=gjet(6);
for q=1:1:numel(bcPrescribeListCell)
    hl(q)=plotV(V(bcPrescribeListCell{q},:),'k.','MarkerSize',markerSize);
    hl(q).Color=plotColors(q,:);
end

legend(hl,{'BC 1','BC 2','BC 3','BC 4','BC 5','BC 6'});
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.symmetric_stiffness=symmetric_stiffness;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
% febio_spec.Material.material{1}.c2=c1;
% febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

% Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='Object1_lattice'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1_lattice';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
for q=1:1:numel(bcPrescribeListCell)
    febio_spec.Mesh.NodeSet{q}.ATTR.name=['bcPrescribeList_',num2str(q)];
    febio_spec.Mesh.NodeSet{q}.node.ATTR.id=bcPrescribeListCell{q}';
end

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section

% -> Prescribe boundary conditions
directionStringSet={'x','x','y','y','z','z'};
displacementMagnitudeDir=[1 -1 1 -1 1 -1];
for q=1:1:numel(bcPrescribeListCell)
    nodeSetName=febio_spec.Mesh.NodeSet{q}.ATTR.name;
    febio_spec.Boundary.bc{q}.ATTR.type='prescribe';
    febio_spec.Boundary.bc{q}.ATTR.node_set=nodeSetName;
    febio_spec.Boundary.bc{q}.dof=directionStringSet{q};
    febio_spec.Boundary.bc{q}.scale.ATTR.lc=1;
    febio_spec.Boundary.bc{q}.scale.VAL=displacementMagnitudeDir(q).*displacementMagnitude;
    febio_spec.Boundary.bc{q}.relative=0;
end

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;

febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sz';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode
%system(['gedit ',febioFebFileName,' &']);

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode=runMode;
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%

    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%

    % Importing nodal forces
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_force),1,1);
    N_force_mat=dataStruct.data;

    indicesSide=bcPrescribeListCell{1};
    areaSide=sampleSize.^2;

    stressVal=mean(squeeze(N_force_mat(indicesSide,1,:))./areaSide,1);
    J_Val=1-((1-J_final).*timeVec(:));

    %%

    cFigure; hold on;
    xlabel('$J$','Interpreter','Latex'); ylabel('$\sigma^*$','Interpreter','Latex');
    plot(J_Val(:),stressVal(:),'r.-','LineWidth',3,'MarkerSize',15);
    axis square; axis tight; grid on; box on;
    set(gca,'FontSize',fontSize);
    drawnow;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    %     gpatch(Fb,Vs,'kw','none',0.25); %A static graphics object
    hp.FaceColor='interp';

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    %     view(130,25); %Set view direction
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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:
%% Subfunction: get a parameter in CSA series ASC header: MrPhoenixProtocol
function val = asc_header(s, key, dft)
if nargin>2, val = dft; else, val = []; end
csa = 'CSASeriesHeaderInfo';
if ~isfield(s, csa) % in case of multiframe
    try s.(csa) = s.SharedFunctionalGroupsSequence.Item_1.(csa).Item_1; end
end
if isfield(s, 'Private_0029_1020') && isa(s.Private_0029_1020, 'uint8')
    str = char(s.Private_0029_1020(:)');
    str = regexp(str, 'ASCCONV BEGIN(.*)ASCCONV END', 'tokens', 'once');
    if isempty(str), return; end
    str = str{1};
elseif isfield(s, 'MrPhoenixProtocol') % X20A
    str = s.MrPhoenixProtocol;
elseif ~isfield(s, csa), return; % non-siemens
elseif isfield(s.(csa), 'MrPhoenixProtocol') % most Siemens dicom
    str = s.(csa).MrPhoenixProtocol;
elseif isfield(s.(csa), 'MrProtocol') % older version dicom
    str = s.(csa).MrProtocol;
else, return;
end

% tSequenceFileName  = ""%SiemensSeq%\gre_field_mapping""
expr = ['\n' regexptranslate('escape', key) '\s*=\s*(.*?)\n'];
str = regexp(str, expr, 'tokens', 'once');
if isempty(str), return; end
str = strtrim(str{1});

if strncmp(str, '""', 2) % str parameter
    val = str(3:end-2);
elseif strncmp(str, '"', 1) % str parameter for version like 2004A
    val = str(2:end-1);
elseif strncmp(str, '0x', 2) % hex parameter, convert to decimal
    val = sscanf(str(3:end), '%x', 1);
else % decimal
    val = sscanf(str, '%g', 1);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function y = cvx_readlevel( x )

global cvx___
s = size( x.basis_ );
[ r, c ] = find( x.basis_ );
y = max( sparse( r, c, cvx___.readonly( r ), s(1), s(2) ), [], 1 );
y = cvx_reshape( y, x.size_ );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [ y2, m2, d2, f2 ] = ymdf_inc_common ( y1, m1, d1, f1, days )

%*****************************************************************************80
%
%% YMDF_INC_COMMON increments a Common YMDF date by DAYS days.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1, the
%    YMDF date.
%
%    Input, real DAYS, the number of days to advance the date.
%
%    Output, integer Y2, M2, D2, real F2, the
%    incremented YMDF date.
%

%
%  Copy the parameters.
%
  y2 = y1;
  m2 = m1;
  d2 = d1;
  f2 = f1 + days;
%
%  Check the parameters.
%
  [ y2, m2, d2, f2, ierror ] = ymdf_check_common ( y2, m2, d2, f2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function File = GetFullPath(File, Style)
% GetFullPath - Get absolute canonical path of a file or folder
% Absolute path names are safer than relative paths, when e.g. a GUI or TIMER
% callback changes the current directory. Only canonical paths without "." and
% ".." can be recognized uniquely.
% Long path names (>259 characters) require a magic initial key "\\?\" to be
% handled by Windows API functions, e.g. for Matlab's FOPEN, DIR and EXIST.
%
% FullName = GetFullPath(Name, Style)
% INPUT:
%   Name:  String or cell string, absolute or relative name of a file or
%          folder. The path need not exist. Unicode strings, UNC paths and long
%          names are supported.
%   Style: Style of the output as string, optional, default: 'auto'.
%          'auto': Add '\\?\' or '\\?\UNC\' for long names on demand.
%          'lean': Magic string is not added.
%          'fat':  Magic string is added for short names also.
%          The Style is ignored when not running under Windows.
%
% OUTPUT:
%   FullName: Absolute canonical path name as string or cell string.
%          For empty strings the current directory is replied.
%          '\\?\' or '\\?\UNC' is added on demand.
%
% NOTE: The M- and the MEX-version create the same results, the faster MEX
%   function works under Windows only.
%   Some functions of the Windows-API still do not support long file names.
%   E.g. the Recycler and the Windows Explorer fail even with the magic '\\?\'
%   prefix. Some functions of Matlab accept 260 characters (value of MAX_PATH),
%   some at 259 already. Don't blame me.
%   The 'fat' style is useful e.g. when Matlab's DIR command is called for a
%   folder with les than 260 characters, but together with the file name this
%   limit is exceeded. Then "dir(GetFullPath([folder, '\*.*], 'fat'))" helps.
%
% EXAMPLES:
%   cd(tempdir);                    % Assumed as 'C:\Temp' here
%   GetFullPath('File.Ext')         % 'C:\Temp\File.Ext'
%   GetFullPath('..\File.Ext')      % 'C:\File.Ext'
%   GetFullPath('..\..\File.Ext')   % 'C:\File.Ext'
%   GetFullPath('.\File.Ext')       % 'C:\Temp\File.Ext'
%   GetFullPath('*.txt')            % 'C:\Temp\*.txt'
%   GetFullPath('..')               % 'C:\'
%   GetFullPath('..\..\..')         % 'C:\'
%   GetFullPath('Folder\')          % 'C:\Temp\Folder\'
%   GetFullPath('D:\A\..\B')        % 'D:\B'
%   GetFullPath('\\Server\Folder\Sub\..\File.ext')
%                                   % '\\Server\Folder\File.ext'
%   GetFullPath({'..', 'new'})      % {'C:\', 'C:\Temp\new'}
%   GetFullPath('.', 'fat')         % '\\?\C:\Temp\File.Ext'
%
% COMPILE:
%   Automatic: InstallMex GetFullPath.c uTest_GetFullPath
%   Manual:    mex -O GetFullPath.c
%   Download:  http://www.n-simon.de/mex
% Run the unit-test uTest_GetFullPath after compiling.
%
% Tested: Matlab 2009a, 2015b(32/64), 2016b, 2018b, Win7/10
%         Compiler: LCC2.4/3.8, BCC5.5, OWC1.8, MSVC2008/2010
% Assumed Compatibility: higher Matlab versions
% Author: Jan Simon, Heidelberg, (C) 2009-2019 matlab.2010(a)n(MINUS)simon.de
%
% See also: CD, FULLFILE, FILEPARTS.

% $JRev: R-M V:038 Sum:C/6JMzUYsYsc Date:19-May-2019 17:25:55 $
% $License: BSD (use/copy/change/redistribute on own risk, mention the author) $
% $UnitTest: uTest_GetFullPath $
% $File: Tools\GLFile\GetFullPath.m $
% History:
% 001: 20-Apr-2010 22:28, Successor of Rel2AbsPath.
% 010: 27-Jul-2008 21:59, Consider leading separator in M-version also.
% 011: 24-Jan-2011 12:11, Cell strings, '~File' under linux.
%      Check of input types in the M-version.
% 015: 31-Mar-2011 10:48, BUGFIX: Accept [] as input as in the Mex version.
%      Thanks to Jiro Doke, who found this bug by running the test function for
%      the M-version.
% 020: 18-Oct-2011 00:57, BUGFIX: Linux version created bad results.
%      Thanks to Daniel.
% 024: 10-Dec-2011 14:00, Care for long names under Windows in M-version.
%      Improved the unittest function for Linux. Thanks to Paul Sexton.
% 025: 09-Aug-2012 14:00, In MEX: Paths starting with "\\" can be non-UNC.
%      The former version treated "\\?\C:\<longpath>\file" as UNC path and
%      replied "\\?\UNC\?\C:\<longpath>\file".
% 032: 12-Jan-2013 21:16, 'auto', 'lean' and 'fat' style.
% 038: 19-May-2019 17:25, BUGFIX, Thanks HHang Li, "File(7:..." -> "File(8:..."

% Initialize: ==================================================================
% Do the work: =================================================================

% #############################################
% ### USE THE MUCH FASTER MEX ON WINDOWS!!! ###
% #############################################

% Difference between M- and Mex-version:
% - Mex does not work under MacOS/Unix.
% - Mex calls Windows API function GetFullPath.
% - Mex is much faster.

% Magic prefix for long Windows names:
if nargin < 2
   Style = 'auto';
end

% Handle cell strings:
% NOTE: It is faster to create a function @cell\GetFullPath.m under Linux, but
% under Windows this would shadow the fast C-Mex.
if isa(File, 'cell')
   for iC = 1:numel(File)
      File{iC} = GetFullPath(File{iC}, Style);
   end
   return;
end

% Check this once only:
isWIN    = strncmpi(computer, 'PC', 2);
MAX_PATH = 260;

% Warn once per session (disable this under Linux/MacOS):
persistent hasDataRead
if isempty(hasDataRead)
   % Test this once only - there is no relation to the existence of DATAREAD!
   %if isWIN
   %   Show a warning, if the slower Matlab version is used - commented, because
   %   this is not a problem and it might be even useful when the MEX-folder is
   %   not inlcuded in the path yet.
   %   warning('JSimon:GetFullPath:NoMex', ...
   %      ['GetFullPath: Using slow Matlab-version instead of fast Mex.', ...
   %       char(10), 'Compile: InstallMex GetFullPath.c']);
   %end

   % DATAREAD is deprecated in 2011b, but still available. In Matlab 6.5, REGEXP
   % does not know the 'split' command, therefore DATAREAD is preferred:
   hasDataRead = ~isempty(which('dataread'));
end

if isempty(File)  % Accept empty matrix as input:
   if ischar(File) || isnumeric(File)
      File = cd;
      return;
   else
      error(['JSimon:', mfilename, ':BadTypeInput1'], ...
         '*** %s: Input must be a string or cell string', mfilename);
   end
end

if ischar(File) == 0  % Non-empty inputs must be strings
   error(['JSimon:', mfilename, ':BadTypeInput1'], ...
      '*** %s: Input must be a string or cell string', mfilename);
end

if isWIN  % Windows: --------------------------------------------------------
   FSep = '\';
   File = strrep(File, '/', FSep);

   % Remove the magic key on demand, it is appended finally again:
   if strncmp(File, '\\?\', 4)
      if strncmpi(File, '\\?\UNC\', 8)
         % [BUGFIX] 19-May-2019, Thanks HHang Li, "File(7:..." -> "File(8:..."
         File = ['\', File(8:length(File))];  % Two leading backslashes!
      else
         File = File(5:length(File));
      end
   end

   isUNC   = strncmp(File, '\\', 2);
   FileLen = length(File);
   if isUNC == 0                        % File is not a UNC path
      % Leading file separator means relative to current drive or base folder:
      ThePath = cd;
      if File(1) == FSep
         if strncmp(ThePath, '\\', 2)   % Current directory is a UNC path
            sepInd  = strfind(ThePath, '\');
            ThePath = ThePath(1:sepInd(4));
         else
            ThePath = ThePath(1:3);     % Drive letter only
         end
      end

      if FileLen < 2 || File(2) ~= ':'  % Does not start with drive letter
         if ThePath(length(ThePath)) ~= FSep
            if File(1) ~= FSep
               File = [ThePath, FSep, File];
            else                        % File starts with separator:
               File = [ThePath, File];
            end
         else                           % Current path ends with separator:
            if File(1) ~= FSep
               File = [ThePath, File];
            else                        % File starts with separator:
               ThePath(length(ThePath)) = [];
               File = [ThePath, File];
            end
         end

      elseif FileLen == 2 && File(2) == ':'   % "C:" current directory on C!
         % "C:" is the current directory on the C-disk, even if the current
         % directory is on another disk! This was ignored in Matlab 6.5, but
         % modern versions considers this strange behaviour.
         if strncmpi(ThePath, File, 2)
            File = ThePath;
         else
            try
               backCD = cd;
               File   = cd(cd(File));
               cd(backCD);
            catch ME
               if exist(File, 'dir')  % No idea what could cause an error then!
                  rethrow(ME);
               else  % Reply "K:\" for not existing disk:
                  File = [File, FSep];
               end
            end
         end
      end
   end

else         % Linux, MacOS: ---------------------------------------------------
   FSep = '/';
   File = strrep(File, '\', FSep);

   if strcmp(File, '~') || strncmp(File, '~/', 2)  % Home directory:
      HomeDir = getenv('HOME');
      if ~isempty(HomeDir)
         File(1) = [];
         File    = [HomeDir, File];
      end

   elseif strncmpi(File, FSep, 1) == 0
      % Append relative path to current folder:
      ThePath = cd;
      if ThePath(length(ThePath)) == FSep
         File = [ThePath, File];
      else
         File = [ThePath, FSep, File];
      end
   end
end

% Care for "\." and "\.." - no efficient algorithm, but the fast Mex is
% recommended at all!
if ~isempty(strfind(File, [FSep, '.']))
   if isWIN
      if strncmp(File, '\\', 2)  % UNC path
         index = strfind(File, '\');
         if length(index) < 4    % UNC path without separator after the folder:
            return;
         end
         Drive            = File(1:index(4));
         File(1:index(4)) = [];
      else
         Drive     = File(1:3);
         File(1:3) = [];
      end
   else  % Unix, MacOS:
      isUNC   = false;
      Drive   = FSep;
      File(1) = [];
   end

   hasTrailFSep = (File(length(File)) == FSep);
   if hasTrailFSep
      File(length(File)) = [];
   end

   if hasDataRead
      if isWIN  % Need "\\" as separator:
         C = dataread('string', File, '%s', 'delimiter', '\\');  %#ok<REMFF1>
      else
         C = dataread('string', File, '%s', 'delimiter', FSep);  %#ok<REMFF1>
      end
   else  % Use the slower REGEXP, when DATAREAD is not available anymore:
      C = regexp(File, FSep, 'split');
   end

   % Remove '\.\' directly without side effects:
   C(strcmp(C, '.')) = [];

   % Remove '\..' with the parent recursively:
   R = 1:length(C);
   for dd = reshape(find(strcmp(C, '..')), 1, [])
      index    = find(R == dd);
      R(index) = [];
      if index > 1
         R(index - 1) = [];
      end
   end

   if isempty(R)
      File = Drive;
      if isUNC && ~hasTrailFSep
         File(length(File)) = [];
      end

   elseif isWIN
      % If you have CStr2String, use the faster:
      %   File = CStr2String(C(R), FSep, hasTrailFSep);
      File = sprintf('%s\\', C{R});
      if hasTrailFSep
         File = [Drive, File];
      else
         File = [Drive, File(1:length(File) - 1)];
      end

   else  % Unix:
      File = [Drive, sprintf('%s/', C{R})];
      if ~hasTrailFSep
         File(length(File)) = [];
      end
   end
end

% "Very" long names under Windows:
if isWIN
   if ~ischar(Style)
      error(['JSimon:', mfilename, ':BadTypeInput2'], ...
         '*** %s: Input must be a string or cell string', mfilename);
   end

   if (strncmpi(Style, 'a', 1) && length(File) >= MAX_PATH) || ...
         strncmpi(Style, 'f', 1)
      % Do not use [isUNC] here, because this concerns the input, which can
      % '.\File', while the current directory is an UNC path.
      if strncmp(File, '\\', 2)  % UNC path
         File = ['\\?\UNC', File(2:end)];
      else
         File = ['\\?\', File];
      end
   end
end

% return;

author: alecjacobson

repo: gptoolbox


text: [T,Y]=ode45('F',[0 1],[0;1;-1])

author: Eurus-Holmes


text:
function [varargout]=quad4_quad8(varargin)

% function [QUAD8,V8,VX8C]=quad4_quad8(QUAD4,V4,VXC)
%
% This function converts 4 node (e.g. linear) quadrilaterial elements into
% 8 node (e.g. quadratic) quadrilateral elements compatible with FEBio.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2018/10/19 Created based on tri3_tri6
%----------------------------------------------------------------------
%%

switch nargin
    case 2
        QUAD4=varargin{1};
        V4=varargin{2};
        VXC={};
    case 3
        QUAD4=varargin{1};
        V4=varargin{2};
        VXC=varargin{3};
end

%%

E=[QUAD4(:,[1 2]); QUAD4(:,[2 3]); QUAD4(:,[3 4]); QUAD4(:,[4 1])]; %Edges matrix
Es=sort(E,2); %Sorted edges matrix
[~,ind1,~]=unique(Es,'rows'); %Indices for unique edges
E=E(ind1,:); %The unique esges

numPoints = size(V4,1);
numEdges = size(E,1);

% Get indices of the three edges associated with each face
A = sparse(E(:,1),E(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
A = max(A,A'); %Copy symmetric

%Indices for A matrix
indA_12=QUAD4(:,1)+(QUAD4(:,2)-1)*numPoints;
indA_23=QUAD4(:,2)+(QUAD4(:,3)-1)*numPoints;
indA_34=QUAD4(:,3)+(QUAD4(:,4)-1)*numPoints;
indA_41=QUAD4(:,4)+(QUAD4(:,1)-1)*numPoints;

%Get indices for vertex array
indV_12=full(A(indA_12));
indV_23=full(A(indA_23));
indV_34=full(A(indA_34));
indV_41=full(A(indA_41));

%Create faces array
QUAD8=[QUAD4(:,1) QUAD4(:,2) QUAD4(:,3) QUAD4(:,4) indV_12 indV_23 indV_34 indV_41];

%Create vertex array
Vn=0.5*(V4(E(:,1),:)+V4(E(:,2),:)); %new mid-edge points
V8 = [V4; Vn]; %Join point sets

%%
varargout{1}=QUAD8;
varargout{2}=V8;

if nargout==3
    %Derive VX8C
    if ~isempty(VXC)
        VX8C=VXC;
        for q=1:1:numel(VXC)
            VX=VXC{q};
            VX_1_4=VX;
            VX_5_8=0.5*(VX(E(:,1),:)+VX(E(:,2),:)); %new mid-edge data
            VX8=[VX_1_4; VX_5_8];
            VX8C{q}=VX8;
        end
    else
        VX8C={};
    end
    varargout{3}=VX8C;
end

%%
% _*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/>.


text:
classdef OPBE < Algorithm
    %OPBE Ordinal Projection Based Ensemble (OPBE)[1]. This class derives
    %from the Algorithm Class and implements the OPBE method with the best
    %configuration found (product combiner, SVM base methodology, logit
    %function and equal distribution of probabilities). By default, this class uses
    %SVORIM implementation, but potentially any ORCA model can be used.
    %
    %   OPBE methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   References:
    %     [1] María Pérez-Ortiz, Pedro Antonio Gutiérrez and César
    %     Hervás-Martínez, Projection based ensemble learning for ordinal
    %     regression, IEEE Transactions on Cybernetics. Vol. 44 (5), 2014
    %     https://doi.org/10.1109/TCYB.2013.2266336

    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original author: María Pérez Ortiz
    %   Contributors: Pedro A. Gutiérrez, Javier Sánchez-Monedero
    %   Citation: If you use this code, please cite the associated papers
    %      - http://www.uco.es/grupos/ayrna/elor2013
    %      - http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html
    properties
        description = 'Ordinal Projection Based Ensemble';
        parameters = [];
        baseMethod = SVORIM;
        % Method name as str to allow loading from configuration file. The
        % constructor will instanciate it via feval()
        base_algorithm = 'SVORIM';
    end

    methods

        % TODO: Update and test parameter description
        function obj = OPBE(varargin)
            %OPBE constructs an object of the class OPBE and sets its
            %default properties.
            %   obj = OPBE('baseMethod', METHOD) sets METHOD as base
            %   algorithm.
            obj.parseArgs(varargin);
            % TODO: Pass varargin parameters to base algorithm?
            obj.baseMethod = feval(obj.base_algorithm);
            %obj.name_parameters = obj.baseMethod.getParameterNames();
            obj.parameters = obj.baseMethod.parameters;
        end

        function [projected, trainTargets] = privfit(obj, train, param)
            %PRIVFIT trains the model for the OPBE method with TRAIN data and
            %vector of parameters PARAMETERS.

            classes = unique(train.targets);
            nOfClasses = numel(classes);
            n = zeros(1,nOfClasses);
            for i=1:nOfClasses
                n(i) = sum(train.targets == i);
            end

            % Sort patterns
            orderedPatterns = train.patterns(train.targets==1,:);
            orderedTargets = train.targets(train.targets == 1);
            for i = 2:nOfClasses
                orderedPatterns = [orderedPatterns ; train.patterns(train.targets==i,:)];
                orderedTargets = [orderedTargets; train.targets(train.targets == i)];
            end

            train.patterns = orderedPatterns;
            train.targets = orderedTargets;
            classBelongingProbTrain = ones(nOfClasses, size(train.patterns,1));
            weights = zeros(nOfClasses,1);

            for i = 1:nOfClasses
                nLowerRankingClasses = sum(classes<i);
                nHigherRankingClasses = sum(classes>i);

                nPreviousClasses = numel(train.targets(train.targets < i));
                nFollowingClasses = numel(train.targets(train.targets > i));

                % Assign labels depending on the decomposition
                if nPreviousClasses == 0
                    currentTargets = [train.targets(train.targets==1); ones(size(train.targets(train.targets>1)))*2];
                elseif nFollowingClasses ==0
                    currentTargets = [ones(size(train.targets(train.targets<i)));  ones(size(train.targets(train.targets==i)))*2];
                else
                    currentTargets = [ones(size(train.targets(train.targets<i))); ones(size(train.targets(train.targets==i)))*2; ones(size(train.targets(train.targets>i)))*3];
                end

                auxtrain.patterns = train.patterns;
                auxtrain.targets = currentTargets;

                % Train each label decomposition
                [projectedTrain] = obj.baseMethod.fit(auxtrain, param);
                models(i) = obj.baseMethod.getModel();

                % Estimate probabilities
                probTrain = obj.calculateProbabilities(projectedTrain, models(i).thresholds');

                % Compute weights and fused probabilities
                for j = 1: nOfClasses
                    if nHigherRankingClasses~= 0 && nLowerRankingClasses ~=0
                        if(j<i)
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .* (probTrain(1,:)/nLowerRankingClasses);
                            weights(j) = weights(j) + 1/nLowerRankingClasses;
                        elseif (j>i)
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .* (probTrain(3,:)/nHigherRankingClasses);
                            weights(j) = weights(j) + 1/nHigherRankingClasses;
                        else
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .* probTrain(2,:);
                            weights(j) = weights(j) + 1;
                        end
                    elseif i==j
                        if nLowerRankingClasses == 0
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .*  probTrain(1,:);
                            weights(j) = weights(j) + 1;
                        else
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .* probTrain(2,:);
                            weights(j) = weights(j) + 1;
                        end
                    else
                        if nLowerRankingClasses == 0
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .* (probTrain(2,:)/nHigherRankingClasses);
                            weights(j) = weights(j) + 1/nHigherRankingClasses;
                        else
                            classBelongingProbTrain(j,:) = classBelongingProbTrain(j,:) .* (probTrain(1,:)/nLowerRankingClasses);
                            weights(j) = weights(j) + 1/nLowerRankingClasses;
                        end
                    end

                end
            end
            model.ensembleModels = models;
            model.parameters = param;
            obj.model = model;

            % Join weights and probabilities into a final
            % decision label
            classBelongingProbTrain = classBelongingProbTrain ./ (weights*ones(1,size(classBelongingProbTrain,2)));
            % There is not a single projection, so projected vector
            % should not be used (it is however needed in the framework)

            % Compute final prediction
            [projected, trainTargets] = max(classBelongingProbTrain);
            projected = projected';
            trainTargets = trainTargets';

        end

        function [projected, testTargets] = privpredict(obj, test)
            %PREDICT predicts labels of TEST patterns labels using model in MODEL.
            models = obj.model.ensembleModels;
            nOfClasses = size(models,2);
            classes = 1:nOfClasses;
            weights = zeros(nOfClasses,1);
            classBelongingProbTest = ones(nOfClasses, size(test,1));

            for i = 1:nOfClasses
                nLowerRankingClasses = sum(classes<i);
                nHigherRankingClasses = sum(classes>i);

                % Estimate probabilities
                obj.baseMethod.setModel(models(i));
                [projectedTest] = obj.baseMethod.predict(test);
                probTest = obj.calculateProbabilities(projectedTest, models(i).thresholds');

                % Compute weights and fused probabilities

                for j = 1: nOfClasses
                    if nHigherRankingClasses~= 0 && nLowerRankingClasses ~=0
                        if(j<i)
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .* (probTest(1,:)/nLowerRankingClasses);
                            weights(j) = weights(j) + 1/nLowerRankingClasses;
                        elseif (j>i)
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .* (probTest(3,:)/nHigherRankingClasses);
                            weights(j) = weights(j) + 1/nHigherRankingClasses;
                        else
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .* probTest(2,:);
                            weights(j) = weights(j) + 1;
                        end
                    elseif i==j
                        if nLowerRankingClasses == 0
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .*  probTest(1,:);
                            weights(j) = weights(j) + 1;
                        else
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .* probTest(2,:);
                            weights(j) = weights(j) + 1;
                        end
                    else
                        if nLowerRankingClasses == 0
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .* (probTest(2,:)/nHigherRankingClasses);
                            weights(j) = weights(j) + 1/nHigherRankingClasses;
                        else
                            classBelongingProbTest(j,:) = classBelongingProbTest(j,:) .* (probTest(1,:)/nLowerRankingClasses);
                            weights(j) = weights(j) + 1/nLowerRankingClasses;
                        end
                    end

                end
            end

            % Join weights and probabilities into a final
            % decision label
            classBelongingProbTest = classBelongingProbTest ./ (weights*ones(1,size(classBelongingProbTest,2)));
            % There is not a single projection, so projected vector
            % should not be used (it is however needed in the framework)

            % Compute final prediction
            [projected, testTargets] = max(classBelongingProbTest);
            projected = projected';
            testTargets = testTargets';

        end

        function [y] = cummulativeProb(obj, x, beta)
            %CUMMULATIVEPROB computes the cummulative probabilities for a
            %set of projected patterns and thresholds
            %
            %   [Y] = CUMMULATIVEPROB(X, BETA) compute cummulative
            %   probabilities Y of projected patterns X with thresholds
            %   BETA

            y =  1 ./ (1+exp((x-beta))); %Logit
            %    y =  1-exp(-exp(beta-x)); %complementary log log
            %    y = exp(-exp(x-beta)); %negative log log
            %    y = normcdf(beta-x); %probit
            %    y = atan(beta-x)/pi + 0.5; % cauchit
        end

        function [g] = calculateProbabilities(obj, z, theta)
            %CALCULATEPROBABILITIES computes the probabilities for a set of
            %projected patterns and thresholds
            %
            %   [G] = CALCULATEPROBABILITIES(OBJ, Z, THETA)
            %   compute probabilities G of projected patterns PROJECTED
            %   with thresholds THRESHOLDS

            % Numerical fix
            nOfClasses = numel(theta)+1;
            if (numel(theta)==2)
                desired=4.0;
                actual=abs(theta(2) - theta(1));
                if actual<4
                    z = z*(desired/actual);
                    theta = theta*(desired/actual);
                end
            end

            f = zeros(nOfClasses, numel(z));
            g = zeros(nOfClasses, numel(z));

            for i=1:(nOfClasses-1)
                f(i,:) = obj.cummulativeProb(z',theta(i));
            end
            f(nOfClasses,:) = ones(1, size(z',2));

            g(1,:) = f(1,:);
            for i=2:nOfClasses
                g(i,:)=f(i,:)-f(i-1,:);
            end

        end

    end

end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function H = gmm_entropy(Priors,Mu,Sigma)
%GMM_ENTROPY: Computes the differentinal entropy of a Gaussian Mixture
%Model.

D = size(Mu,1);
K = size(Mu,2);
H = 0;

if D == 1
    for k=1:K
        H = H  + Priors(k) * (-log(Priors(k)) + 0.5 * log(((2*pi*exp(1))^(D)) * det(Sigma(k)) ));
    end
else
    for k=1:K
        H = H  + Priors(k) * (-log(Priors(k)) + 0.5 * log(((2*pi*exp(1))^(D)) * det(Sigma(:,:,k)) ));
    end
end

end

author: epfl-lasa


text:
function [value, varname] = in_matlab_var(varname, vartype)
% IN_MATLAB_VAR: Reads the contents of a Matlab variable in base workspace
%
% USAGE:  value = in_matlab_var(varname, vartype)
%         value = in_matlab_var(varname, 'numeric')
%         value = in_matlab_var(varname)
%         value = in_matlab_var()        : Ask the variable name to the user

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2015

value = [];
% Variable type not specified
if (nargin < 2) || isempty(vartype)
    vartype = [];
end
% Ask user the name of the destination variable
if (nargin < 1) || isempty(varname)
    % Get the list of all the base workspace variables (structures only)
    allVars = evalin('base', 'whos');
    if isempty(allVars)
        bst_error('No variables in the base workspace.', 'Import from Matlab', 0);
        return;
    end
    % If variable type not specified, show all variables
    if isempty(vartype)
        listVar = {allVars.name};
    elseif isequal(vartype, 'numeric')
        iVar = find(ismember({allVars.class}, {'uint8','int8','uint16','int16','uint32','int32','uint64','int64','single','double'}));
        listVar = {allVars(iVar).name};
    else
        iVar = find(strcmpi({allVars.class}, vartype));
        listVar = {allVars(iVar).name};
    end
    % Check again that there are variables availables
    if isempty(listVar)
        bst_error('No valid variables in the base workspace.', 'Import from Matlab', 0);
        return;
    end
    % Show question
    varname = java_dialog('combo', 'Select a workspace variable:

', 'Import from Matlab', [], listVar);
    if isempty(varname)
        return
    end
end
% Get variable value
try
    value = evalin('base', varname);
catch
    bst_error(['Variable "' varname '" does not exist in Matlab workspace.'], 'Import from workspace', 0);
    value = [];
end

author: brainstorm-tools


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

%*****************************************************************************80
%
%% NCC_COMPUTE computes a Newton-Cotes Closed quadrature rule.
%
%  Discussion:
%
%    For the interval [-1,+1], the Newton-Cotes Closed quadrature rule
%    estimates
%
%      Integral ( -1 <= X <= +1 ) F(X) dX
%
%    using N abscissas X and weights W:
%
%      sum ( 1 <= I <= N ) W(I) * F ( X(I) ).
%
%    For the CLOSED rule, the abscissas are equally spaced, and include
%    the end points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = ncc_compute_points ( n );

  w = ncc_compute_weights ( n );

  return
end


text:
function [Xc_opt,R_opt,T_opt,err_opt,iter]=optimize_betas_gauss_newton(Km,Cw,Beta0,Alph,Xw,U,A)

% COMPUTE_BETAS_GAUSS_NEWTON
%
%       Km: vector of the kernel
%       Cw: position of the control point in world coordinates
%       Beta0: initial guess of the betas
%
%
% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% 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/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, October 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

n=size(Beta0,2);

[Beta_opt,err,iter]=gauss_newton(Km,Cw,Beta0);

%Extract control point camera coordinates from Betas and Kernel
X=zeros(12,1);
for i=1:n
   X=X+Beta_opt(i)*Km(:,i);
end

Cc=zeros(4,3);
for i=1:4
    Cc(i,:)=X(3*i-2:3*i);
end

%check sign of the determinant (keep orienation of the control points)
s_Cw=sign_determinant(Cw);
s_Cc=sign_determinant(Cc);
Cc=Cc*(s_Cw/s_Cc);

%Reconstruct and compute error=
Xc_opt=Alph*Cc; %reconstruction: points in camera coordinate system
[R_opt,T_opt]=getrotT(Xw,Xc_opt);  %solve exterior orientation

[err_opt,Urep_opt]=reprojection_error_usingRT(Xw,U,R_opt,T_opt,A);

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
%  INTERNAL FUNCTION: different kernel densities
%
%

author: jmaih


text:
function [ud,UD_up,UD_dist] = distort(up,dist)

% DISTORT  Distort projected point with radial distortion.
%   DISTORT(UP,DIST) computes the position in the image plane of the
%   distorted pixel corresponding to the projected pixel UP. DIST is the
%   vector of distortion parameters. Distortion is calculated with the
%   radial distortion model:
%
%     UD = UP * (1 + k2*r^2 + k4*r^4 + k6*r^6 + ...)
%
%   where r^2  = UP(1)^2 + UP(2)^2
%     and DIST = [k4 k4 k6 ...]
%
%   [ud,UD_up,UD_dist] = DISTORT(...) returns Jacobians wrt UP and DIST.
%
%   See also PINHOLE.
%

%   FIXME: RMAX works not so nice
%   DISTORT(...,RMAX) accepts the maximum image radius to limit the effect
%   of distortion to the image plane. The image radius is one half of the
%   image diagonal. In case norm(UP) > RMAX then UD = RMAX/norm(UP)*UP. The
%   Jacobians are trivially computed just for argument consistency: UDup =
%   eye(2) and UDdist = zeros(2,length(dist)).

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

r2 = sum(up.^2);

n = length(dist);

if n == 0  % No distortion
    ud = up;

    if size(up,2) == 1
        UD_up   = eye(2);
        UD_dist = zeros(2,0);
    else
        error('??? Jacobians not available for multiple points.')
    end

else       % Distortion

    if size(up,2) == 1   % One only pixel

        nn    = 1:n;           % orders vector
        r2n   = r2.^nn;        % powers vector
        ratio = 1 + r2n*dist;  % distortion ratio

        ud    = ratio*up;

        if nargout > 1 % jacobians

            dnn    = 2:2:2*n;
            dr2n   = r2.^(0:n-1);
            dratio = (dnn.*dr2n)*dist;

            UD_up = ratio*eye(2) + dratio*(up*up');

            UD_dist = kron(up,r2n);
        end

    else  % vectorized operation for multiple pixels.

        ratio = ones(size(r2));
        for i=1:n
            ratio = ratio + dist(i)*r2.^i;
        end

        ud = up.*[ratio;ratio];

        if nargout > 1 % jacobians

            error('??? Jacobians not available for multiple points.')

        end
    end
end

return

%% Build jacobians up to order 3
% For higher orders:
%   1/ add k8,k10,k(2n)... to syms line below for an order n.
%   2/ add k8,k10,k(2n)... to vector dist below for an order n.
%   3/ Execute this cell (Apple+return keys)
%   4/ Copy the results in UDup and UDdist above, caring for the matrices'
%   opening and closing brackets (do not forget them).

syms u1 u2 k2 k4 k6 k8 real
up   = [u1;u2];
dist = [k2;k4;k6;k8];
ud   = distort(up,dist);

UDup   = simple(jacobian(ud,up))
UDdist = simple(jacobian(ud,dist))

%% test

[ud, UD_up, UD_dist] = distort(up,dist);

simplify(UD_up - jacobian(ud,up))
simplify(UD_dist - jacobian(ud,dist))

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
%%
%  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File:      lownerjohn_inner.m
%
%  Purpose:
%  Computes the Lowner-John inner ellipsoidal
%  approximation of a polytope.
%
%
%  The inner ellipsoidal approximation to a polytope
%
%     S = { x \in R^n | Ax < b }.
%
%  maximizes the volume of the inscribed ellipsoid,
%
%     { x | x = C*u + d, || u ||_2 <= 1 }.
%
%  The volume is proportional to det(C)^(1/n), so the
%  problem can be solved as
%
%    maximize         t
%    subject to       t       <= det(C)^(1/n)
%                || C*ai ||_2 <= bi - ai^T * d,  i=1,...,m
%                  C is PSD
%
%  which is equivalent to a mixed conic quadratic and semidefinite
%  programming problem.
%
%
%  References:
%  [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
%
function [C, d] = lownerjohn_inner(A, b)

import mosek.fusion.*;
import mosek_lownerjohn.det_rootn;

M = Model('lownerjohn_inner');

[m, n] = size(A);

% Setup variables
t = M.variable('t', 1, Domain.greaterThan(0.0));
C = M.variable('C', NDSet(n,n), Domain.unbounded());
d = M.variable('d', n, Domain.unbounded());

% (bi - ai^T*d, C*ai) \in Q
for i=1:m,
    M.constraint( sprintf('qc%d', i),  ...
                  Expr.vstack(Expr.sub(b(i) ,Expr.dot(A(i,:),d) ) , Expr.mul(C,A(i,:) )), ...
                  Domain.inQCone() );
end

% t <= det(C)^{1/n}
det_rootn(M, C, t)

% Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t)

M.solve()

C = reshape(C.level(), n, n);
d = reshape(d.level(), n, 1);

M.dispose();

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function x = p12_start ( option, nvar )

%*****************************************************************************80
%
%% P12_START returns a starting point for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Output, real X(NVAR), the starting point.
%
  x = zeros ( nvar, 1 );

  return
end


text:
function spm_progress_bar(action,varargin)
% Display a 'Progress Bar' in the 'Interactive' window
% FORMAT spm_progress_bar('Init',height,xlabel,ylabel,flgs)
% Initialises the bar in the 'Interactive' window.
% If flgs contains a 't', then use tex interpreter for labels.
%
% FORMAT spm_progress_bar('Set',value)
% Sets the height of the bar itself.
%
% FORMAT spm_progress_bar('Clear')
% Clears the 'Interactive' window.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

if ~nargin, action = 'Init'; end

% Find the interactive window and exit if not
%-----------------------------------------------------------------------
Finter = spm_figure('FindWin','Interactive');
if isempty(Finter), return; end

switch lower(action)
    % Initialise
    %-------------------------------------------------------------------
    case 'init'
        error(nargchk(0,5,nargin));
        if nargin > 1, arg1 = varargin{1}; else arg1 = 1;           end
        if nargin > 2, arg2 = varargin{2}; else arg2 = 'Computing'; end
        if nargin > 3, arg3 = varargin{3}; else arg3 = '';          end
        if nargin > 4, arg4 = varargin{4}; else arg4 = ' ';         end
        if any(arg4 == 't'), interp = 'tex'; else interp = 'none';  end
        pb = struct('pointer',get(Finter,'Pointer'),...
                    'name'   ,get(Finter,'Name'),...
                    'buffer', get(Finter,'DoubleBuffer'));
        spm_progress_bar('Clear');
        set(Finter,'Pointer','watch');
        set(Finter,'Name',pb.name);
        set(Finter,'DoubleBuffer','on');
        pb.ax = axes('Position', [0.45 0.2 0.05 0.6],...
                     'XTick',    [],...
                     'Xlim',     [0 1],...
                     'Ylim',     [0 max([arg1 eps])],...
                     'Box',      'on',...
                     'Parent',   Finter);
        lab = get(pb.ax,'Xlabel');
        set(lab,'string',arg2,'FontSize',10,'Interpreter',interp);
        lab = get(pb.ax,'Ylabel');
        set(lab,'string',arg3,'FontSize',10,'Interpreter',interp);
        lab = get(pb.ax,'Title');
        set(lab,'string','0% Complete','Interpreter',interp);
        t = clock;
        str = sprintf('Began %2.0f:%02.0f:%02.0f',t(4),t(5),t(6));
        text(2,arg1/2,0,str,'FontSize',10,'Parent',pb.ax);
        l = line('Xdata',     [0.5 0.5],...
                 'Ydata',     [0 0],...
                 'LineWidth', 8,...
                 'Color',     [1 0 0],...
                 'Tag',       'ProgressBar',...
                 'Parent',    pb.ax);
        set(l,'UserData',pb);
        drawnow;

    % Set
    %-------------------------------------------------------------------
    case 'set'
        error(nargchk(1,2,nargin));
        if nargin == 1, value = 0; else  value = varargin{1}; end
        br = findobj(Finter,'Tag','ProgressBar');
        if ~isempty(br)
            pb = get(br,'UserData');
            set(br,'Ydata',[0 value]);
            lim = get(get(br,'Parent'),'Ylim');lim=lim(2);
            lab = get(pb.ax,'Title');
            set(lab,'string',sprintf('%.0f%% Complete',100*value/lim));
            drawnow;
        end

    % Clear
    %-------------------------------------------------------------------
    case 'clear'
        error(nargchk(1,1,nargin));
        pb = get(findobj(Finter,'Tag','ProgressBar'),'UserData');
        spm_figure('Clear',Finter);
        if isstruct(pb)
            set(Finter,'Pointer',     pb.pointer);
            set(Finter,'Name',        pb.name);
            set(Finter,'DoubleBuffer',pb.buffer);
        end
        drawnow;

    % Error
    %-------------------------------------------------------------------
    otherwise
        error('Unknown action string');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function M = spm_mesh_polyhedron(name)
% Return one of the Platonic solids with triangle faces
% FORMAT M = spm_mesh_polyhedron(name)
% name     - polyhedron name
%            (one of {'tetrahedron','octahedron','icosahedron'})
%
% M        - patch structure
%__________________________________________________________________________
%
% See https://www.wikipedia.org/wiki/Platonic_solid
%__________________________________________________________________________
% Copyright (C) 2012-2018 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_polyhedron.m 7383 2018-07-31 10:53:37Z guillaume $

switch lower(name)

    case 'tetrahedron'
        M.vertices = [  1   0  -1/sqrt(2) ;
                       -1   0  -1/sqrt(2) ;
                        0   1   1/sqrt(2) ;
                        0  -1   1/sqrt(2) ];
        M.faces    = [  1   2   3 ;
                        1   2   4 ;
                        2   3   4 ;
                        3   1   4];

    case 'octahedron'
        M.vertices = [  0   0   1 ;
                        1   0   0 ;
                        0   1   0 ;
                       -1   0   0 ;
                        0  -1   0 ;
                        0   0  -1 ];
        M.faces    = [  1   2   3 ;
                        1   3   4 ;
                        1   4   5 ;
                        1   5   2 ;
                        6   3   2 ;
                        6   4   3 ;
                        6   5   4 ;
                        6   2   5 ];

    case 'icosahedron'
        r = (1 + sqrt(5)) / 2;
        M.vertices = [  1   0  -r ;
                        0   r  -1 ;
                        r   1   0 ;
                        r  -1   0 ;
                        0  -r  -1 ;
                       -1   0  -r ;
                       -r   1   0 ;
                        0   r   1 ;
                        1   0   r ;
                        0  -r   1 ;
                       -r  -1   0 ;
                       -1   0   r ];
        M.faces    = [  1   2   3
                        3   4   1
                        1   4   5
                        5   6   1
                        1   6   2
                        2   6   7
                        7   8   2
                        2   8   3
                        3   8   9
                        9   4   3
                        4   9  10
                        4  10   5
                       10  11   5
                        5  11   6
                        6  11   7
                        7  12   8
                        8  12   9
                        9  12  10
                       10  12  11
                       11  12   7 ];

    otherwise
        error('''%s'' is not a Platonic solid with triangle faces.',name);
end

author: spm


text:
function dist = getDistPCKh(pred,gt,refDist)

assert(size(pred,1) == size(gt,1) && size(pred,2) == size(gt,2) && size(pred,3) == size(gt,3));
assert(size(refDist,1) == size(gt,3));

dist = nan(1,size(pred,2),size(pred,3));

for imgidx = 1:size(pred,3)

    % distance to gt joints
    dist(1,:,imgidx) = sqrt(sum((pred(:,:,imgidx) - gt(:,:,imgidx)).^2,1))./refDist(imgidx);

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function [  ] = tapas_BLT_setup(  )
% setup for Breathing Learning Task (BLT)
%
% author: Yu Yao, Translational Neuromodeling Unit (TNU), University Zurich
%         and ETH Zurich

% get path to tapas
dTool = mfilename('fullpath');
dTool = fileparts(dTool);
dTapas = fileparts(dTool);
dTapas = fileparts(dTapas);
% add downloader to path
addpath(fullfile(dTapas,'misc'));
addpath(fullfile(dTapas,'external'));
try
    % download example data
    [exdir] = tapas_download_example_data();
catch mx
    if strcmp(mx.identifier,'tapas:example_data:target_exists')
        exdir = mx.message;
        exdir(1:10) = [];
        exdir(end-15:end) = [];
    else
        rethrow(mx);
    end
end

% copy to task folder
copyfile(fullfile(exdir,'BLT','images'),fullfile(dTool,'images'));
disp('done')

end

author: translationalneuromodeling

repo: tapas


text:
function h = image_histogram(P,varargin)
% :Usage:
% ::
%
%    h = image_histogram(P,[method(string)],[range])
%
% eliminates 0, NaN voxels from either image
% red then blue
%
% :Inputs:
%
%   **P:**
%        is one or two images in string array
%
%   **Methods:**
%        'def'
%
%   **range:**
%        optional 3rd input, range of values to include in histogram
%
% :Example:
% ::
%
%    h = image_histogram('p_Omnibus.img','def',[0 1-eps]);
%
% ..
%    tor wager
% ..

V = spm_vol(P); %v = spm_read_vols(V);

v = iimg_read_vols(V);

if length(varargin) > 0, meth = varargin{1};  else, meth = 'def'; end
if length(varargin) > 1, range = varargin{2};  else, range = [-Inf Inf]; end

for f = 1:size(P,1), [dum,ff{f}] = fileparts(P(f,:)); end

%p2 = which('/Users/tor/Documents/tor_scripts/3DheadUtility/canonical_brains/scalped_avg152T1_graymatter.img')
%[p2,p2new] = reslice_imgs(P(1,:),p2,1);
%maskV = spm_read_vols(spm_vol(p2new)); maskV = maskV(:);

switch meth
    case 'def'

        disp('Plotting histograms')

        v1 = squeeze(v(:,:,:,1));
        if size(v,4) > 1
            v2 = squeeze(v(:,:,:,2));
        end

        v1 = v1(:);
        if size(v,4) > 1, v2 = v2(:); ,end

        if size(v,4) > 1
            wh = isnan(v1) | isnan(v2) | v1 == 0 | v2 == 0 ...
                | v1 < range(1) | v1 > range(2) | v2 < range(1) | v2 > range(2); %| ~maskV;
        else
            wh = isnan(v1) | v1 == 0 | v1 < range(1) | v1 > range(2);
        end

        v1(wh) = [];
        if size(v,4) > 1
            v2(wh) = [];
        else
            v2 = [];
        end

        [h,xx] = hist([v1;v2],max(10,length(v1)./300));
        h1 = hist(v1,xx);

        if size(v,4) > 1, h2 = hist(v2,xx); end

    otherwise
        error('unknown action string')
end

figure('Color','w');

if size(v,4) > 1, plot(xx,h1,'r',xx,h2,'b');
else
        plot(xx,h1,'r');
end

set(gca,'FontSize',16)

legend(ff)

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
%
% 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/>.

% -*- texinfo -*-
% @deftypefn {Function File} {@var{BW} = } octave_poly2mask (@var{x},@var{y},@var{m},@var{n})
% Convert a polygon to a region mask.
%
% BW=octave_poly2mask(x,y,m,n) converts a polygon, specified by a list of
% vertices in @var{x} and @var{y} and returns in a @var{m}-by-@var{n}
% logical mask @var{BW} the filled polygon. Region inside the polygon
% is set to 1, values outside the shape are set to 0.
%
% @var{x} and @var{y} should always represent a closed polygon, first
% and last points should be coincident. If they are not octave_poly2mask will
% close it for you. If @var{x} or @var{y} are fractional they are
% nearest integer.
%
% If all the polygon or part of it falls outside the masking area
% (1:m,1:n), it is discarded or clipped.
%
% This function uses scan-line polygon filling algorithm as described
% in http://www.cs.rit.edu/~icss571/filling/ with some minor
% modifications: capability of clipping and scan order, which can
% affect the results of the algorithm (algorithm is described not to
% reach ymax, xmax border when filling to avoid enlarging shapes). In
% this function we scan the image backwards (we begin at ymax and end
% at ymin), and we don't reach ymin, xmin, which we believe should be
% compatible with MATLAB.
% @end deftypefn

% TODO: check how to create a logical BW without any conversion

function BW = octave_poly2mask (x, y, m, n)
  if (nargin ~= 4)
    print_usage ();
  end

  % check x and y
  x = round (x (:).');
  y = round (y (:).');
  if (length (x) < 3)
    error ('octave_poly2mask: polygon must have at least 3 vertices.');
  end
  if (length (x) ~= length (y))
    error ('octave_poly2mask: length of x doesn''t match length of y.');
  end

  % create output matrix
  BW = false (m, n);

  % close polygon if needed
  if ((x (1) ~= x (length (x))) || (y (1) ~= y (length (y))))
    x = horzcat (x, x (1));
    y = horzcat (y, y (1));
  end

  % build global edge table
  ex = [x(1:length (x) - 1); x(1, 2:length (x))]; % x values for each edge
  ey = [y(1:length (y) - 1); y(1, 2:length (y))]; % y values for each edge
  idx = (ey(1, :) ~= ey(2, :));                 % eliminate horizontal edges
  ex = ex (:, idx);
  ey = ey (:, idx);
  eminy = min (ey);                               % minimum y for each edge
  emaxy = max (ey);                               % maximum y for each edge
  t = (ey == [eminy; eminy]);                     % values associated to miny
  exvec = ex(:);
  exminy = exvec(t);                            % x values associated to min y
  exmaxy = exvec(~t);                           % x values associated to max y
  emaxy = emaxy.';                                % we want them vertical now...
  eminy = eminy.';
  m_inv = (exmaxy - exminy)./(emaxy - eminy);     % calculate inverse slope
  ge = [emaxy, eminy, exmaxy, m_inv];             % build global edge table
  ge = sortrows (ge, [1, 3]);                     % sort on eminy and exminy

  % we add an extra dummy edge at the end just to avoid checking
  % while indexing it
  ge = [-Inf, -Inf, -Inf, -Inf; ge];

  % initial parity is even (0)
  parity = 0;

  % init scan line set to bottom line
  sl = ge (size (ge, 1), 1);

  % init active edge table
  % we use a loop because the table is sorted and edge list could be
  % huge
  ae = [];
  gei = size (ge, 1);
  while (sl == ge (gei, 1))
    ae = [ge(gei, 2:4); ae];
    gei = gei - 1;
  end

  % calc minimum y to draw
  miny = min (y);
  if (miny < 1)
    miny = 1;
  end

  while (sl >= miny)
    % check vert clipping
    if (sl <= m)
      % draw current scan line
      % we have to round because 1/m is fractional
      ie = round (reshape (ae (:, 2), 2, size (ae, 1)/2));

      % this discards left border of image (this differs from version at
      % http://www.cs.rit.edu/~icss571/filling/ which discards right
      % border) but keeps an exception when the point is a vertex.
      ie (1, :) = ie (1, :) + (ie (1, :) ~= ie (2, :));

      % we'll clip too, just in case m,n is not big enough
      ie (1, (ie (1, :) < 1)) = 1;
      ie (2, (ie (2, :) > n)) = n;

      % we eliminate segments outside window
      ie = ie (:, (ie (1, :) <= n));
      ie = ie (:, (ie (2, :) >= 1));
      for i = 1:size(ie,2)
        BW (sl, ie (1, i):ie (2, i)) = true;
      end
    end

    % decrement scan line
    sl = sl - 1;

    % eliminate edges that eymax==sl
    % this discards ymin border of image (this differs from version at
    % http://www.cs.rit.edu/~icss571/filling/ which discards ymax).
    ae = ae ((ae (:, 1) ~= sl), :);

    % update x (x1=x0-1/m)
    ae(:, 2) = ae(:, 2) - ae(:, 3);

    % update ae with new values
    while (sl == ge (gei, 1))
      ae = vertcat (ae, ge (gei, 2:4));
      gei = gei - 1;
    end

    % order the edges in ae by x value
    if (size(ae,1) > 0)
      ae = sortrows (ae, 2);
    end
  end
end

% This should create a filled octagon
%!demo
%! s = [0:pi/4:2*pi];
%! x = cos (s) * 90 + 101;
%! y = sin (s) * 90 + 101;
%! bw = octave_poly2mask(x, y, 200, 200);
%! imshow (bw);

% This should create a 5-vertex star
%!demo
%! s = [0:2*pi/5:pi*4];
%! s = s ([1, 3, 5, 2, 4, 6]);
%! x = cos (s) * 90 + 101;
%! y = sin (s) * 90 + 101;
%! bw = octave_poly2mask (x, y, 200, 200);
%! imshow (bw);

%!# Convex polygons

%!shared xs, ys, Rs, xt, yt, Rt
%! xs=[3,3,10,10];
%! ys=[4,12,12,4];
%! Rs=zeros(16,14);
%! Rs(5:12,4:10)=1;
%! Rs=logical(Rs);
%! xt=[1,4,7];
%! yt=[1,4,1];
%! Rt=[0,0,0,0,0,0,0;
%!     0,0,1,1,1,1,0;
%!     0,0,0,1,1,0,0;
%!     0,0,0,1,0,0,0;
%!     0,0,0,0,0,0,0];
%! Rt=logical(Rt);

%!assert(octave_poly2mask(xs,ys,16,14),Rs);          # rectangle
%!assert(octave_poly2mask(xs,ys,8,7),Rs(1:8,1:7));   # clipped
%!assert(octave_poly2mask(xs-7,ys-8,8,7),Rs(9:16,8:14)); # more clipping

%!assert(octave_poly2mask(xt,yt,5,7),Rt);            # triangle
%!assert(octave_poly2mask(xt,yt,3,3),Rt(1:3,1:3));   # clipped

%!# Concave polygons

%!test
%! x=[3,3,5,5,8,8,10,10];
%! y=[4,12,12,8,8,11,11,4];
%! R=zeros(16,14);
%! R(5:12,4:5)=1;
%! R(5:8,6:8)=1;
%! R(5:11,9:10)=1;
%! R=logical(R);
%! assert(octave_poly2mask(x,y,16,14), R);

%!# Complex polygons
%!test
%! x=[1,5,1,5];
%! y=[1,1,4,4];
%! R=[0,0,0,0,0,0;
%!    0,0,1,1,0,0;
%!    0,0,1,1,0,0;
%!    0,1,1,1,1,0;
%!    0,0,0,0,0,0];
%! R=logical(R);
%! assert(octave_poly2mask(x,y,5,6), R);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function r = parse_formula(varargin)

% PARSE_FORMULA Parses a chemical formula to form an atomic representation.
%
% SYNTAX
%
% r = parse_formula(str)
% r = parse_formula({str1,str2,str3, ...})
%
%   Parses chemical formulas and returns a structure array holding the an
%   atomic representation of the chemical forulas. The input is a string or
%   a cell array of strings.
%
%
% EXAMPLES
%
%   1. Chemical formulas of varying complexity
%
%       parse_formula('H2O');            % Water
%       parse_formula('NaHCO3');         % Sodium Bicarbonate
%       parse_formula('(CH4)8(H2O)46');  % Methane Clathrate
%       parse_formula('CH3COOCH2CH3');   % Ethyl Acetate
%       parse_formula('MnO4-');          % Negative Charge Ion
%
%       parse_formula('dCH4');           % Returns error message
%
%   2. Create an structure array of atomic representations for a set of
%      compounds
%
%       r = parse_formula({'CH4','O2','CO2','H2O'});
%
%
% USAGE NOTES
%
%   1. Formulas are made of up of sequences of elements followed by
%      integers  indicating the number of included atoms. Omitted integers
%      are assumed to be one.
%
%   2. Elements are the conventional one or two character abbreviations.
%      The character is captialized. If present, the second character is
%      lower case. In addition to the standard elements, the parser allows
%      for
%
%       Symbol  Entity                 Interpretation
%          e    electron               like an element with MW = 0
%          D    deuterium              an element
%          T    tritium                an element
%          M    any metal              like an element, mw = NaN
%          X    any halogen            like an element, mw = NaN
%          Me   methyl group (CH3)     CH3 substituted for Me
%          Et   ethyl group (C2H5)     C2H5 substituted for Et
%          Bu   butyl group (C4H9)     C4H9 substituted for Bu
%          Ph   phenol group (C6H5)    C6H5 substituted for Ph
%
%   3. Subgroups may be included between parenthesis or brackets followed
%      by an integer indicating number of repetitions. Two levels of
%      subgrouping are allowed.
%
%   4. A terminal lower case suffix denoting phases will be correctly
%      parsed. The phase must be one of (aq), (l), (g), or (s).
%
%   5. The charge on an ionic species is appended as a + or - followed by
%      an optional integer.  Examples are H+, OH-, or Ca+2.
%
%   6. The bare electron e- is used in balancing chemical half reactions.
%
%   7. Error messages are generated for invalid fomulas
%      Change for use with the cobra toolbox
%      A is accepted as the symbol for an element that could
%      represent an R group in an atom mapping scenario for an incompletely
%      defined metabolite.
%
%   8. str can be a cell array of chemical formula. The results is a
%      structure array. The elements of the output structure array are in
%      one-to-one correspondence with elements of the cell array. For
%      example
%
%          r = parse_formula({'CH4','CH3OH','CHOOH'})
%
%      r(1) holds the atomic formula for CH4, r(2) for CH3OH, and r(3) for
%      CHOOH.

% AUTHOR
%
%   Jeff Kantor
%   December 18, 2010

    assert(nargin > 0, 'parse_formula:input', ['No input. Expects a  ', ...
                        'string or cell array of chemical formulas.']);
    assert(nargin < 2, 'stoich:input', 'Unexpected extra inputs.');

    switch class(varargin{1})
        case 'char'                      % Single formula
            str = varargin;

        case 'cell'                      % Cell array of formulas
            str = varargin{1};

        otherwise
            error('parse_formula:input',['requires cell array of  ',...
              'chemical formulas.']);
    end

    assert(iscellstr(str), 'parse_formula:input', ...
        'Formulas must be strings.');

    % Trim any whitespace at front or back

    str = strtrim(str);

    % Remove phase information. This information is currently neglected. In
    % a later version we may wish to incorporate phase into a more complete
    % data structure for representing chemical formula.

    prex = '|\((aq|g|l|s)\)$';
    str = regexprep(str,prex,'');

    % Substitute for some common chemical abbreviations

    str = regexprep(str,'Bu','C4H9');    % Butyl
    str = regexprep(str,'Et','C2H5');    % Ethyl
    str = regexprep(str,'Me','CH3');     % Methyl
    str = regexprep(str,'Ph','C6H5');    % Phenol

    % Apply the main parser to every element of str

    q = cellfun(@(s)parse_formula_(s,3),str,'Uniform',false);

    % Union of all atomic species

    atoms = {};
    for i = 1:length(q(:))
        atoms = union(atoms, fields(q{i}));
    end

    % Add all atomic species to all structures.

    for i = 1:length(q(:))
        for j = 1:length(atoms)
            if ~ismember(atoms{j},fields(q{i}))
                q{i}.(atoms{j}) = 0;
            end
        end
    end

    % Form the structure array to have the same shape as str

    r = reshape([q{:}],size(str));

end % parse_formula

function r = parse_formula_(str,kdepth)

    assert(kdepth > 0, 'parse_formula_:Recursion', ...
        'Reached maximum recursion depth');

    r = struct([]);

    % Regular expression returning tokens for element and number
    % sexpr matches single elements followed by a digit, or a +/-
    % followed by a digit to denote charge

    persistent srex;  % Regexp pattern to match elements and charges
    persistent grex;  % Regexp pattern to match groups

    if isempty(srex) || isempty(grex)
        %this was the original line
        %srex = ['(A[lrsgutcm]|B[eraik]?|C[laroudsemf]?|D[y]?|E[urs]|',...
        srex = ['(A[lrsgutcm]?|B[eraik]?|C[laroudsemf]?|D[y]?|E[urs]|', ...
                'F[erm]?|G[aed]|H[eofgas]?|I[nr]?|Kr?|L[iaur]|', ...
                'M[gnodt]?|N[eaibdpos]?|Os?|P[drmtboau]?|R[buhenaf]|', ...
                'S[icernbmg]?|T[icebmalh]?|U|V|W|X[e]?|Yb?|Z[nr])', ...
                '(\d*\.\d+|\d*)', ...
                '|(e|+|-)(\d*)'];
        grex = '|\(([^\)]*)\)(\d*\.\d+|\d*)|\[([^\]]*)\](\d*\.\d+|\d*)';
    end

    % Parse formula for chemical groups. This picks out anything that looks
    % an element followed by a number, or a subgroup within parentheses.
    % The tokens are returned in the cell array u. Each u{k} has two
    % elements, the first is a string denoting the group, and the second is
    % number string of repetitions.

    [u,s,e] = regexp(str,[srex,grex],'tokens','start','end');

    % Report any parsing errors. A parse error occurs if there are any
    % characters not matched as tokens. We scan the start and end positions
    % of the tokens to determine if there are any gaps.

    g(1:length(str)) = '^';
    for i = 1:length(s);
        g(s(i):e(i)) = ' ';
    end

    assert(all(g ~= '^'), 'parse_formula:ParseError', ...
        'Could not parse formula:\n    %s\n    %s\n', str, char(g));

    % Extract atom tokens from the first part of each token

    tok = cellfun(@(v)v{1},u,'Uni',false);

    % Extract counts from the second part of each token, convert to
    % doubles, empty counts set to 1

    cnt = cellfun(@(v)v{2},u,'Uni',false);
    cnt = str2double(cnt);
    cnt(isnan(cnt)) = 1;

    % Loop over tokens

    for j = 1:length(u)

        % See if token matches an element

        if strcmp(tok{j},regexp(tok{j},srex,'match'))

            % The token exactly matches an element.
            % Change + or - tokens to 'Q'.

            tok{j} = regexprep(tok{j},'+','Q');

            if strcmp(tok{j}, '-')
                tok{j} = 'Q';
                cnt(j) = -cnt(j);
            end

            % Update atomic representation, adding a field if needed.

            if isfield(r,tok{j})
                r.(tok{j}) = r.(tok{j}) + cnt(j);
            else
                r(1).(tok{j}) = cnt(j);
            end

        else

            % The token must be a group, so do a recursion to find
            % an atomic represenation of the group.

            q = parse_formula_(tok{j},kdepth-1);

            % Updatethe  atomic representation to include the group.
            % Add fields if needed. Multiply by number of groups in the
            % formula we're parsing.

            f = fields(q);

            for k = 1:length(f)

                if isfield(r,f{k})
                    r.(f{k}) = r.(f{k}) + cnt(j)*q.(f{k});
                else
                    r(1).(f{k}) = cnt(j)*q.(f{k});
                end

            end
        end
    end

end % parse_formula_

author: opencobra

repo: cobratoolbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [y,dy] = affine2Dsparse(w,x,varargin)
%
% computes y = kron(I2,Q)*w = [w(1),w(2);w(4),w(5)] * x + w([3,6])
% and the derivative wrt. w.
% Q = {[x(:,1),x(:,2),1]}, dy = Q;
% if no argumanets are given, the parameters for the identity map are returned.
%
% see also transformations/contents.m, trafo.m
%==============================================================================

function [y,dy] = affine2Dsparse(w,x,varargin)

% the persitent variable stores the matrix
% Q(x) = kron( I_2 , [x(:,1),x(:,2),1] );
persistent Q

if nargin == 0,
  runMinimalExample;
  return;
else
  y = mfilename('fullfile');
  dy = [1;0;0;0;1;0];         % parameterization of identity
  if ischar(w),  Q  = []; w = []; end; % reset Q
  if isempty(w), return;          end;
end;

% test for need of updating Q
OK = ~any(isempty(Q));
if OK,
  m1 = size(Q{1});
  OK = (2*m1(1) == size(x,1)) && (2*m1(2) == numel(w));
end;

if ~OK,
  n = length(x)/2;
  Q = {[reshape(x,[],2),ones(n,1)]};
  if nargout == 0, return; end;
end;

w  = reshape(w,3,2);
% mimicing Qfull*w as [Q*w(:,1),Q*w(:,2)]
y  = [Q{1}*w(:,1);Q{1}*w(:,2)];
dy = Q;

%------------------------------------------------------------------------------
function runMinimalExample
help(mfilename);
fprintf('%s: minimal example\n',mfilename)

omega = [0,10,0,2]; m = [8,9];
w = 22/pi;c = (omega(2:2:end)-omega(1:2:end))'/2;
R = [ cos(w),-sin(w);sin(w),cos(w)];
g = (eye(2)-R)*reshape(c,2,1);
w = [R(1,1);R(1,2);g(1);R(2,1);R(2,2);g(2)]
x = getNodalGrid(omega,m);
z = feval(mfilename,w,x);
t = affine2D(w,x);
n = norm(z-t);
fprintf('diff between affine2D and affine2Dsparse is %s\n',num2str(n));

assert(n<1e-13,'difference to non-sparse version too big')
FAIRfigure(1); clf;
plotGrid(x,omega,m,'color','r'); axis image; hold on;
plotGrid(z,omega,m,'color','b'); axis image; hold off;
%==============================================================================