text:
function [Im_DN,D_] = Denoise(Im_N,Sigma,K,n,Algo)
% INPUT ARGUMENTS : Im_N - the noisy image (gray-level scale)
%                   Sigma - the s.d. of the noise (assume to be white Gaussian).
%                   K - the number of atoms in the representing dictionary.
%                   n - Block Size n x n
%% Data Stuff
Reduce_DC = 1;
[N1,N2] = size(Im_N);
C = 1.15;                       % Taken from the Paper Elad 06
E_T = C*Sigma;                  % Required Average target Error in OMP
if Sigma > 5
    noIt = 10;
else
    noIt = 5;
end
MaxTPatches = 62001;

%% DCT Dictionary Creation
% D_DCT = Dict_DCT(n,K);
if Reduce_DC
    D_DCT = odctdict(n^2,K+1);
    D_DCT = D_DCT(:,2:end);
else
    D_DCT = odctdict(n^2,K);
end
%% Block Patches to Columns
fprintf('Learning Dictionary using %s\n',Algo);
if(prod([N1,N2]-n+1)> MaxTPatches)
    randPermutation =  randperm(prod([N1,N2]-n+1));
    selectedBlocks = randPermutation(1:MaxTPatches);

    Y = zeros(n^2,MaxTPatches);
    for i = 1:MaxTPatches
        [row,col] = ind2sub(size(Im_N)-n+1,selectedBlocks(i));
        currBlock = Im_N(row:row+n-1,col:col+n-1);
        Y(:,i) = currBlock(:);
    end
else
    Y = im2col(Im_N,[n,n],'sliding');       % Signal Y with patches as Columns
end

%% Reducing DC component
if (Reduce_DC)
    vecOfMeans = mean(Y);
    Y = Y-ones(size(Y,1),1)*vecOfMeans;
end

%% Data Whitening
% Y = Data_Whiten(Y,1);

%% Going into Dictionary Learning Algo for Training
D_ = normc(D_DCT);
for it = 1:noIt
    W = omp2(D_,Y,D_'*D_,(n*E_T));    %
    switch lower(Algo)
        case 'ksvd'
            [D_,W] = Optimize_K_SVD(Y,D_,W);
        case 's1'
            alpha = .37;
            [D_,W] = Optimize_S1(Y,D_,W,alpha);
        case 's1svd'
            alpha = 100;    gamma = 1;
            [D_,W] = Optimize_S1SVD(Y,D_,W,alpha,gamma);
        otherwise
            error('Invalid Learning Method Specified');
            break;
    end
    D_ = I_clearDictionary(D_,W,Y);
    disp(['Iteration # ',num2str(it),' With average number of Coefficients = ',num2str(nnz(W)/size(W,2))])
end

%% DEnoising with the Trained Dictionary
disp('Denoising with the Trained Dictionary')
Y = im2col(Im_N,[n,n],'sliding');
if (Reduce_DC)
    vecOfMeans = mean(Y);
    Y = Y-ones(size(Y,1),1)*vecOfMeans;
end
Coefs = omp2(D_,Y,D_'*D_,(n*E_T));
if (Reduce_DC)
    Y = D_*Coefs + ones(size(Y,1),1) * vecOfMeans;
else
    Y = D_*Coefs;
end

%% Generation and Averaging of The signal from Y (columns)
count = 1;
Weight= zeros(N1,N2);
IMout = zeros(N1,N2);
idx = 1:size(Y,2);
[rows,cols] = ind2sub(size(Im_N)-n+1,idx);
for i  = 1:length(cols)
    col = cols(i); row = rows(i);
    Y_ = reshape(Y(:,count),[n,n]);
    IMout(row:row+n-1,col:col+n-1) = IMout(row:row+n-1,col:col+n-1)+Y_;
    Weight(row:row+n-1,col:col+n-1) = Weight(row:row+n-1,col:col+n-1)+ones(n);
    count = count+1;
end;
Im_DN = (Im_N+0.034*Sigma*IMout)./(1+0.034*Sigma*Weight);

end

%% KSVD implementation for Elad 2006
function [D_,W] = Optimize_K_SVD(Y1,D_,W)
    R = Y1 - D_*W;
    for k=1:size(D_,2)
        I = find(W(k,:));
        Ri = R(:,I) + D_(:,k)*W(k,I);
        [U,S,V] = svds(Ri,1,'L');
        D_(:,k) = U;
        W(k,I) = S*V';
        R(:,I) = Ri - D_(:,k)*W(k,I);
    end
end

%% S1 implementation
function [D,W] = Optimize_S1(Y,D,W,alpha)
    Ek = Y - D * W;
    for k = 1:size(D,2)
        Eki = Ek + D(:,k)*W(k,:);
        for j = 1:2
            G = D(:,k)'*Eki;   g = std(G);  G = G./g; % std = 1
            %alpha2 = std(G)*alpha;
            W(k,:) = g.*sign(G).*max(0,abs(G)-alpha);
            D(:,k) = (Eki * W(k,:)')/norm(Eki * W(k,:)');
        end
%         nnz(W(k,:))
        Ek = Eki - D(:,k)*W(k,:);
    end
end

%% S1 WIth SVD Dictionary Update Stage
function [D,W] = Optimize_S1SVD(Y,D,W,alpha,gamma)
    Ek = Y - D * W;
    for k = 1:size(D,2)
        Eki = Ek + D(:,k)*W(k,:);
        % SVD
%         [D(:,k),s,v] = svds(Eki,1,'L');
%         W(k,:) = s*v';
        % Power Iteration
        for i = 1:5
            W(k,:) = D(:,k)'*Eki;
            D(:,k) = Eki*W(k,:)';    D(:,k) = D(:,k)/norm(D(:,k));
        end
        for j = 1:2
            G = D(:,k)'*Eki;  g=1; %g = std(G);  G = G./g;
            alpha2 = std(G)*alpha;
            W(k,:) = g.* sign(G).*max(0,(abs(G) - g.* alpha2./(abs(W(k,:)))));  %.^gamma
            D(:,k) = Eki*W(k,:)'/norm(Eki*W(k,:)');
        end
%             nnz(W(k,:))
        Ek = Eki - D(:,k)*W(k,:);
    end
end


text:
function [results] = hyperHud(M, B, S)
% HYPERHUD Performs the hybrid unstructured detector (HUD) algorithm
%   Performs the hybrid unstructured detector algorithm for target
% detection.
%
% Usage
%   [results] = hyperHud(M, B, S)
% Inputs
%   M - 2d matrix of HSI data (p x N)
%   B - 2d matrix of background endmembers (p x q)
%   S - 2d matrix of target endmembers (p x #target_sigs)
% Outputs
%   results - vector of detector output (N x 1)
%
% References
%   J Broadwater & R Chellappa.  "Hybrid Detectors for Subpixel Targets."
% IEEE PAMI. Vol 29. No 11. November 2007.

[p, N] = size(M);
% Remove mean from data
u = mean(M.').';
M = M - repmat(u, 1, N);
S = S - repmat(u, 1, size(S,2));

numTargets = size(S,2);
%sigma = 1e-5;
E = [S B];
%E = [sigma*E; ones(1,size(E,2))];
q = size(E, 2);

R_hat = (M*M.')/N;
G = inv(R_hat);

results = zeros(1, N);

R = ones(q,1);
P = R - 1;
% TODO - put in the whitened version of fcls
a_hat_tmp = hyperNnls(M, E);
%a_hat_tmp = hyperFcls(M, E);
for k=1:N
    x = M(:,k);
    a_hat = a_hat_tmp(:,k);
    % Take the top r values from a_hat where r is number of targets.  We
    % are only interested in the abundances for the targets.  From J
    % Broadwater email 11/17/09.
    a_hat = a_hat(1:numTargets);
% %    x = [sigma*x; 1];
%     % FCLS optimzation
%     lambda = zeros(q,1);
%     aPrev = lambda;
%     for kk=1:100
%         a_hat = inv(E.'*G*E)*E.'*G*x - inv(E.'*G*E)*lambda;
%         norm(a_hat-aPrev)
%         lambda = E.'*G*(x-E*a_hat);
%         idx = find(a_hat>0);
%         P(idx) = 1;
%         R(idx) = 0;
%         aPrev = a_hat;
%     end
%     a_hat = a_hat(1:numTargets);
    results(k) = (x.'*G*S*a_hat) / (x.'*G*x);
end

author: davidkun

repo: HyperSpectralToolbox


text:
function T = argmin_g(w0, zeta, X, T)
     lhd= 1 ./  (w0 .^2 + zeta); % left hand

     % compute T for each channel
     for i = 1:size(X,3)
         T(:,:,i) = lhd .* X(:,:,i);
     end
end

repo: STRCF


text:
function [g1, g2, g3] = sdlfmaXsdlfmvKernGradientBlock(lfmKern1, lfmKern2, ...
    t1, t2, kyy, kyv, kvy, kvv, i, j, generalConst, generalConstGrad, ...
    covGrad)

% SDLFMAXSDLFMVKERNGRADIENTBLOCK Gradients of the parameters in block i,j
% FORMAT
% DESC computes the kernel parameters gradients for the SDLFM kernel
% function in the block specified at indeces i,j. It assumes the
% computation for functions that system 1 describes an acceleration and
% system 2 describes a velocity.
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG kyy : covariance for the initial conditions between position 1 and
% position 2 at block i,j
% ARG kyv : covariance for the initial conditions between position 1 and
% velocity 2 at block i,j
% ARG kvy : covariance for the initial conditions between velocity 1 and
% position 2 at block i,j
% ARG kvv : covariance for the initial conditions between velocity 1 and
% velocity 2 at block i,j
% ARG i : interval to be evaluated for system 1
% ARG j : interval to be evaluated for system 2
% ARG generalConstant : constants evaluated with sdlfmKernComputeConstant.m
% ARG generalConstGrad : derivatives of the constants computed with
% sdlfmKernGradientConstant.m
% ARG covGrad : partial derivatives of the objective function wrt portion
% of the kernel matrix in block i,j
% RETURN g1 : gradients of parameters for the system 1
% RETURN g2 : gradients of parameters for the system 2
% RETURN g3 : gradients of switching points
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin<11
    j = i;
    generalConst = [];
end

g3 = [];

% Compute derivatives of the mean terms with respect to the parameters

[g1Mean, g2Mean, gsp1Mean, gsp2Mean] = sdlfmKernGradientMean(lfmKern1(1), ...
    lfmKern2(1), t1, t2, kyy, kyv, kvy, kvv, covGrad, {'sdlfma','sdlfmv'}, ...
    {'sdlfmj', 'sdlfma'});

if i==j
    [g1, g2, g3t] = sdlfmaXsdlfmKernGradientBlockIEJ(lfmKern1, lfmKern2, t1, ...
        t2, covGrad, g1Mean, g2Mean, gsp1Mean, gsp2Mean, {'lfma', 'lfmv'}, ...
        {'lfmj', 'lfmv'}, {'lfma', 'lfma'});
    g3(i) = g3t;
else
    if i>j
        [g1, g2, g3] = sdlfmXsdlfmvKernGradientBlockIGJ(lfmKern1, lfmKern2, ...
            t1, t2, i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
            g2Mean, gsp1Mean, gsp2Mean, 'sdlfma', 'sdlfmj', 'lfmvXlfm', ...
            'lfmvXlfmv', {'lfmvXlfmv', 'lfmaXlfm'}, {'lfmaXlfmv', 'lfmaXlfmv'});
    else
        [g1, g2, g3] = sdlfmaXsdlfmKernGradientBlockILJ(lfmKern1, lfmKern2, ...
            t1, t2, i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
            g2Mean, gsp1Mean, gsp2Mean, 'sdlfmv', 'sdlfma', 'lfmaXlfm', ...
            'lfmaXlfmv',{'lfmjXlfm', 'lfmaXlfmv'}, {'lfmjXlfmv', 'lfmaXlfma'});
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%Copyright 2013 The MathWorks, Inc
if ii==1,
    %% Global Coordinate System
    figure('WindowStyle','docked');
    hold on; grid on
    antT = nan(3,N);      % Array trajectory
    tgtT = nan(3,N);      % Target trajectory
    antT = [antT(:,2:end) antpos/1e3];
    tgtT = [tgtT(:,2:end) tgtpos/1e3];
    hant = plot3(antT(1,:),antT(2,:),antT(3,:),'LineWidth',2);
    htgt = plot3(tgtT(1,:),tgtT(2,:),tgtT(3,:),'LineWidth',2,'Color','red');
    axis([-100 100 -100 50 0 50])
    xlabel('X'), ylabel('Y'), zlabel('Z')
    title('Global Coordinate System','FontWeight','bold')
    view(40,48)
    legend('Antenna','Target')
    %% Target Trajectory in Local Coordinate System (from antenna perspective)
    h2 = figure('WindowStyle','docked');
    hold on; grid on
    axis([-130 0 6 13 0 50])
    xlabel('Azimuth'), ylabel('Range'), zlabel('Elevation')
    title('Local Coordinate System','FontWeight','bold')
    view(-40,48)
    %% Array visualization
    figure('WindowStyle','docked');
    viewArray(sAnt,'ShowSubarray','None')
    bx = get(gca,'Children');
    hdots = bx(8);
    ndots = numel(get(hdots,'XData'));
    green = ones(ndots,1)*[0 1 0];
    view(40,48)
else
    %% Highlight active arrays
    mn=face(:)*ones(1,ndots);
    hlgt=repmat([mn(1,1:16) mn(2,17:32) mn(3,33:48) mn(4,49:64)],3,1);
    set(hdots,'CData',green.*hlgt');
    %% Update global trajectories
    antT = [antT(:,2:end) antpos/1e3];
    tgtT = [tgtT(:,2:end) tgtpos/1e3];
    set(hant,'XData',antT(1,:),'YData',antT(2,:),'ZData',antT(3,:))
    set(htgt,'XData',tgtT(1,:),'YData',tgtT(2,:),'ZData',tgtT(3,:))
end
%% Update local trajectory
figure(h2)
plot3(AzEl(1),range/1e4,AzEl(2),'ro','MarkerSize',20-round(range/1e4))

author: Sable

repo: mcbench-benchmarks


text:
classdef LvdPluginOptimVarWrapper < matlab.mixin.SetGet
    %LvdPluginOptimVarWrapper Summary of this class goes here
    %   Detailed explanation goes here

    properties
        name(1,:) char = 'Untitled Plugin Variable';
        value(1,1) double = 0;

        optVar PluginOptimizationVariable

        id(1,1) double = 0
    end

    methods
        function obj = LvdPluginOptimVarWrapper()
            obj.id = rand();
        end

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

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

        function [lb, ub] = getBounds(obj)
            lb = obj.optVar.lb;
            ub = obj.optVar.ub;
        end

        function setBounds(obj, lb, ub)
            obj.optVar.lb = lb;
            obj.optVar.ub = ub;
        end

        function useTf = isVariableActive(obj)
            useTf = obj.optVar.useTf;
        end

        function setIfVariableIsActive(obj, useTf)
            obj.optVar.useTf = useTf;
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            if(isempty(obj.id) || obj.id == 0)
                obj.id = rand();
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% GP_PREDICT_PSEUDO - Computes the predictive distribution of a Gaussian
%                     process using pseudo inputs.
%
% Usage:
%
%   F_MEAN = GP_PREDICT_PSEUDO(Y, V, K_FP, K_PP, K_PH, K_H)
%   [F_MEAN, F_VAR] = GP_PREDICT_PSEUDO(...)
%
% The marginal of a conditional model
%
%   Y|F ~ N(F, V)
%
% is approximated by
%
%   Y ~ N(0, K_FP*INV(K_PP)*K_FP' + V)
%
% Y : Nx1 vector of observations
% V : NxN diagonal matrix of noise covariance, or a function V(X) which
%     solves V\X
% K_FP : NxM covariance matrix of function values at observation inputs
%        and at pseudo inputs
% K_PP : MxM covariance matrix of function values at pseudo inputs
% K_PH : MxK covariance matrix of function values at pseudo inputs and at
%        inputs to be predicted
% K_H : Kx1 vector of variances of function values at inputs to be
%       predicted
%
% See also GP_LEARN_PSEUDO, GP_PREDICT.

% Last modified 2011-01-27
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function [f_mean, f_var] = gp_predict_pseudo(y, V, K_pf, K_pp, K_ph, k_h)

if isnumeric(V)
  % Get a solver for V\z
  inv_V = get_linsolve_cov(V);
else
  % A function handle which solves V\z
  inv_V = V;
end

% TODO:
%
% - Take K_pf instead of K_fp
%
% - Take V as a vector?

[L_p,p] = chol(K_pp, 'lower');
if p~=0
  figure
  imagesc(K_pp);
  error('Matrix not positive definite');
end

Z_f = linsolve_tril(L_p, K_pf);

Lambda = speye(size(K_pp)) + Z_f*inv_V(Z_f');
inv_Lambda = get_linsolve_cov(Lambda);

f_mean = K_ph' * linsolve_triu(L_p, inv_Lambda(Z_f*inv_V(y)), true);
if nargout >= 2
  Z_h = linsolve_tril(L_p, K_ph);
  f_var = k_h(:) - dot(Z_h, Z_h - inv_Lambda(Z_h),1)';
end

% $$$ Lambda = K_pp + K_fp' * inv_V(K_fp);
% $$$ inv_Lambda = get_linsolve_cov(Lambda);
% $$$ inv_Kpp = get_linsolve_cov(K_pp);
% $$$ f_mean = K_ph' * inv_Lambda(K_fp' * inv_V(y));
% $$$ f_var = k_h(:) - dot(K_ph, inv_Kpp(K_ph) - inv_Lambda(K_ph),1)';

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%% febioStruct2xml
% Below is a demonstration of the features of the |febioStruct2xml| function

%%
clear; close all; clc;

%% Syntax
% |[domNode]=febioStruct2xml(febio_spec,fileName,optionStruct);|

%% Description
% This function provides the basis for coding FEBio XML input file content
% and helps to generate the XML input file.

%%
% See also: |febioStructTemplate|

%% About FEBio XML input files
% FEBio input files (.feb) are XML files. The XML syntax bellow is for
% febio_spec 2.5. More on the input file formatting can be found in the
% FEBio user manual.

%% Coding XML files in MATLAB
% Coding XML files (or in this case .feb FEBio files) can be challenging
% and may look confusing. An alternative method is implemented here which
% is a trade-off between clarify, ease of programming, and fidelity with
% the generated XML file. In general XML files follow the following form:

%%

% <Globals>
%   <Constants>
%       <R>8.314e-6</R>
%       <T>298</T>
%       <Fc>96485e-9</Fc>
%   </Constants>
%   <Solutes>
%       <solute id="1" name="Na">
%           <charge_number>1</charge_number>
%           <molar_mass>22.99</molar_mass>
%           <density>0.97<density>
%       </solute>
% </Globals>

%%
% Here Globals is an element with the entries Constants and Solutes being
% child elements (and Constants and Solutes have R,T,Fc and solute as their
% respective child elements,... etc). The symbols <..> start an element and
% </..> closes that element. Components within <..> sometimes contain
% attributes. An example of this is id="1", here the attribute "id" has the
% value "1" associated with it, another example of an attribute is
% name="Na". Elements may contain child elements and may have a value. For
% instance the element Constants has children but the element R has the
% value 8.314e-6 associated with it.
% To code these features in MATLAB a particular syntax was developed.
% Structures are suitable as they are tree-like similar to XML. However
% structures lack the attribute feature, i.e. fields in structures have
% values but do not have attributes associated with them. To solve this
% keywords are used here to denote attributes and to distinguish them from
% values or elements. For instance node.ATTR.id=1 will become <node
% id="1">, where the ATTR keyword is here simply used as a target/trigger
% keyword for determining wheter the subsequent feature is an element,
% value or attribute. Futhermore the VAL keyword exists to denote a value,
% e.g. node.ATTR.id=1; node.VAL=2;  will become: <node id="1"> 2 </node>.

%%
% One may use arrays for values as well i.e. node.ATTR.id=nodeIdArray;
% node.VAL=nodeCoordinateArray; can be used to define sets of attributes
% and sets of values simultaneously.
%%
% In principle all FEBio functionality (based on XML content) can be coded
% for in a similar fashion. Study the FEBio user manual for more
% information.

%%
% The remainder of this documentation will describe how various XML
% components of the FEBio input file can be coded. A structure is created
% with the name "febio_spec". For more information on coding the febio
% structure see the various FEBio demos provided.

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

%% Module section
% XML syntax example:
% <Module type="Solid"/>

%%
% Implementation example
febio_spec.Module.ATTR.type='solid';

%% Control section

%% Example 1: Globals section

% <Globals>
%   <Constants>
%       <R>8.314e-6</R>
%       <T>298</T>
%       <Fc>96485e-9</Fc>
%   </Constants>
%   <Solutes>
%       <solute id="1" name="Na">
%           <charge_number>1</charge_number>
%           <molar_mass>22.99</molar_mass>
%           <density>0.97<density>
%       </solute>
%       <solute id="2" name="Cl">
%           <charge_number>-1</charge_number>
%           <molar_mass>35.45</molar_mass>
%           <density>3.21<density>
%       </solute>
%       <solute id="3" name="Glc">
%           <charge_number>0</charge_number>
%           <molar_mass>180.16</molar_mass>
%           <density>1.54<density>
%       </solute>
%   </Solutes>
% </Globals>

febio_spec.Globals.Constants.T=0;
febio_spec.Globals.Constants.R=0;
febio_spec.Globals.Constants.Fc=0;

febio_spec.Globals.Solutes{1}.ATTR.id=1;
febio_spec.Globals.Solutes{1}.ATTR.name='Na';
febio_spec.Globals.Solutes{1}.charge_number=1;
febio_spec.Globals.Solutes{1}.molar_mass=22.99;
febio_spec.Globals.Solutes{1}.density=0.97;

febio_spec.Globals.Solutes{2}.ATTR.id=2;
febio_spec.Globals.Solutes{2}.ATTR.name='Cl';
febio_spec.Globals.Solutes{2}.charge_number=-1;
febio_spec.Globals.Solutes{2}.molar_mass=35.45;
febio_spec.Globals.Solutes{2}.density=3.21;

%% Example 2: Material section

%%
% An uncouple material
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=1;
febio_spec.Material.material{1}.m1=6;
febio_spec.Material.material{1}.k=100;

%%
% A visco-elastic material

%Viscoelastic part
febio_spec.Material.material{2}.ATTR.type='uncoupled viscoelastic';
febio_spec.Material.material{2}.ATTR.Name='Block_material';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.g1=0.1;
febio_spec.Material.material{2}.t1=0.01;
febio_spec.Material.material{2}.density=1e-9;

%Elastic part
febio_spec.Material.material{2}.elastic{1}.ATTR.type='Ogden';
febio_spec.Material.material{2}.elastic{1}.c1=1;
febio_spec.Material.material{2}.elastic{1}.m1=2;
febio_spec.Material.material{2}.elastic{1}.k=100;
febio_spec.Material.material{2}.elastic{1}.density=1e-9;

%%
% A poro-elastic material

%Viscous part
febio_spec.Material.material{3}.ATTR.type='biphasic';
febio_spec.Material.material{3}.ATTR.name='Block_material';
febio_spec.Material.material{3}.ATTR.id=3;
febio_spec.Material.material{3}.phi0=0.5;
febio_spec.Material.material{3}.permeability.ATTR.type='perm-const-iso';
febio_spec.Material.material{3}.permeability.ATTR.name='permeability';
febio_spec.Material.material{3}.permeability.perm=7.41e-11; %hydraulic permeability
febio_spec.Material.material{3}.fluid_density=1e-9;

%Solid part
febio_spec.Material.material{3}.solid{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{3}.solid{1}.c1=1;
febio_spec.Material.material{3}.solid{1}.m1=2;
febio_spec.Material.material{3}.solid{1}.cp=100;
febio_spec.Material.material{3}.solid{1}.density=1e-9;

%%
% A poro-elastic and visco-elastic material

% 	<Material>
% 		<material id="1" name="Biphasic-Viscoelastic" type="biphasic">
% 			<phi0>0.2</phi0>
% 			<solid type="viscoelastic">
% 				<g1>0.5</g1>
% 				<t1>3000</t1>
% 				<elastic type="solid mixture">
% 					<solid type="neo-Hookean">
% 						<density>1</density>
% 						<E>0.5</E>
% 						<v>0</v>
% 					</solid>
% 					<solid type="ellipsoidal fiber distribution">
% 						<beta>2,2,2</beta>
% 						<ksi>5,5,5</ksi>
% 					</solid>
% 				</elastic>
% 			</solid>
% 			<permeability type="perm-Holmes-Mow">
% 				<perm>0.001</perm>
% 				<M>3</M>
% 				<alpha>2</alpha>
% 			</permeability>
% 		</material>

%Viscous part
febio_spec.Material.material{4}.ATTR.type='biphasic';
febio_spec.Material.material{4}.ATTR.name='Biphasic-Viscoelastic';
febio_spec.Material.material{4}.ATTR.id=4;
febio_spec.Material.material{4}.phi0=0.2;
febio_spec.Material.material{4}.permeability.ATTR.type='perm-Holmes-Mow';
febio_spec.Material.material{4}.permeability.perm=0.001;
febio_spec.Material.material{4}.permeability.M=3;
febio_spec.Material.material{4}.permeability.alpha=2;

%Solid part
%-> Viscoelastic part
febio_spec.Material.material{4}.solid{1}.ATTR.type='viscoelastic';
febio_spec.Material.material{4}.solid{1}.g1=0.5;
febio_spec.Material.material{4}.solid{1}.t1=3000;
febio_spec.Material.material{4}.solid{1}.density=1e-9;

%-> Elastic part
febio_spec.Material.material{4}.solid{1}.elastic{1}.ATTR.type='solid mixture';
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{1}.ATTR.type='neo-Hookean';
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{1}.density=1e-9;
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{1}.E=0.5;
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{1}.v=0;
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{2}.ATTR.type='ellipsoidal fiber distribution';
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{2}.beta=[2 2 2];
febio_spec.Material.material{4}.solid{1}.elastic{1}.solid{2}.ksi=[5 5 5];

%% Example 3: Mesh section

%Example geometry
[meshStruct]=hexMeshBox([2 2 2],[2 2 2]);
E=meshStruct.E; V=meshStruct.V;
bcSupportList=(1:size(V,1))';

% -> Nodes

% <Nodes name="Object1">
%   <node id="1">0,0,0,</node>
%   <node id="2">1,0,0,</node>
%   ...
% </Nodes>

febio_spec.Mesh.Nodes{1}.ATTR.name='Object1'; %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

% <Elements type="hex8" name="Part1">
%   <elem id="1">1,2,3,4,5,6,7,8</elem>
%   <elem id="2">1,2,3,4,9,10,11,12</elem>
%   <...>
% </Elements>

partName1='Part1';
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
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

%% Creating the FEBio input file
% You can use |febioStruct2xml| to write the xml data to a file AND/OR a
% domNode object (optional output). Leave the fileName variable empty to
% supress file export.

%Create file name for XML file
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');
fileName=fullfile(savePath,'tempModel.feb');

optionStruct.arrayParseMethod=1;
[domNode]=febioStruct2xml(febio_spec,fileName,optionStruct);

%% Viewing the FEBio input file
% The |febView| command can be used to render an XML file in a figure
% window. Alternatively the |textView| command can be used:
h=febView(fileName,2);

%%
%
% <<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:
function wavelet_test ( )

%*****************************************************************************80
%
%% WAVELET_TEST tests the WAVELET library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WAVELET_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the WAVELET library.\n' );

  wavelet_test01 ( );
  wavelet_test02 ( );
  wavelet_test03 ( );
  wavelet_test04 ( );
  wavelet_test05 ( );
  wavelet_test06 ( );
  wavelet_test07 ( );
  wavelet_test08 ( );
  wavelet_test09 ( );

  wavelet_test10 ( );
  wavelet_test11 ( );
  wavelet_test12 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WAVELET_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function tests = test_imSurfaceArea(varargin)
% Test case for the file imSurfaceArea.
%
%   output = test_imSurfaceArea(input)
%
%   Example
%   test_imSurfaceArea
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-10-08,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRAE - Cepia Software Platform.

tests = functiontests(localfunctions);

function testAddBorderD3(testCase)

img = ones([5 5 5]) > 0;
imgb = padarray(img, [1 1 1]) > 0;

nDir = 3;
b = imSurfaceArea(img, nDir);
bb = imSurfaceArea(imgb, nDir);
assertEqual(testCase, b, bb);

function testAddBorderD13(testCase)

img = ones([5 5 5]) > 0;
imgb = padarray(img, [1 1 1]) > 0;

nDir = 13;
b = imSurfaceArea(img, nDir);
bb = imSurfaceArea(imgb, nDir);
assertEqual(testCase, b, bb);

function test_Anisotropic(testCase)

img = ones([5 5 5]) > 0;
imgb = padarray(img, [1 1 1]) > 0;

nDir = 3;
b = imSurfaceArea(img, nDir, [1 2 3]);
bb = imSurfaceArea(imgb, nDir, [1 2 3]);
assertEqual(testCase, b, bb);

function test_Labels_D3(testCase)

img = zeros([5 5 5]);
img(1:3, 1:3, 1:3) = 1;
img(4:5, 1:3, 1:3) = 2;
img(1:3, 4:5, 1:3) = 3;
img(4:5, 4:5, 1:3) = 4;
img(1:3, 1:3, 4:5) = 5;
img(4:5, 1:3, 4:5) = 6;
img(1:3, 4:5, 4:5) = 7;
img(4:5, 4:5, 4:5) = 8;

b3 = imSurfaceArea(img, 3);
assertEqual(testCase, 8, length(b3));

imgb = padarray(img, [1 1 1]);
bb3 = imSurfaceArea(imgb, 3);

assertEqual(testCase, b3, bb3);

function test_Labels_D13(testCase)

img = zeros([5 5 5]);
img(1:3, 1:3, 1:3) = 1;
img(4:5, 1:3, 1:3) = 2;
img(1:3, 4:5, 1:3) = 3;
img(4:5, 4:5, 1:3) = 4;
img(1:3, 1:3, 4:5) = 5;
img(4:5, 1:3, 4:5) = 6;
img(1:3, 4:5, 4:5) = 7;
img(4:5, 4:5, 4:5) = 8;

b13 = imSurfaceArea(img, 13);
assertEqual(testCase, 8, length(b13));

imgb = padarray(img, [1 1 1]);
bb13 = imSurfaceArea(imgb, 13);
assertEqual(testCase, b13, bb13);

function test_Labels2_D3(testCase)

img = zeros([5 5 5]);
img(1:3, 1:3, 1:3) = 2;
img(4:5, 1:3, 1:3) = 3;
img(1:3, 4:5, 1:3) = 5;
img(4:5, 4:5, 1:3) = 7;
img(1:3, 1:3, 4:5) = 11;
img(4:5, 1:3, 4:5) = 13;
img(1:3, 4:5, 4:5) = 17;
img(4:5, 4:5, 4:5) = 19;

[b3, labels3] = imSurfaceArea(img, 3);

assertEqual(testCase, 8, length(b3));
assertEqual(testCase, 8, length(labels3));
assertEqual(testCase, labels3, [2 3 5 7 11 13 17 19]');

function test_Labels2_D13(testCase)

img = zeros([5 5 5]);
img(1:3, 1:3, 1:3) = 2;
img(4:5, 1:3, 1:3) = 3;
img(1:3, 4:5, 1:3) = 5;
img(4:5, 4:5, 1:3) = 7;
img(1:3, 1:3, 4:5) = 11;
img(4:5, 1:3, 4:5) = 13;
img(1:3, 4:5, 4:5) = 17;
img(4:5, 4:5, 4:5) = 19;

[b13, labels13] = imSurfaceArea(img, 13);

assertEqual(testCase, 8, length(b13));
assertEqual(testCase, 8, length(labels13));
assertEqual(testCase, labels13, [2 3 5 7 11 13 17 19]');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [W B] = box_decomp(U)
%BOX_DECOMP

% TODO: doc
% centering is recommended

% 2^n vertex!

M = [min(U); max(U)];
L = M(2,:) - M(1,:);

[n1 n2] = size(U);

%%% TODO: this is a hack
% if n2==2
% 	A = [linspace(0,1,n1)' linspace(1,0,n1)'];
% 	if norm(U*(U'*A)-A) < 1e-10
% 		% linear
% 		W = A;
% 		B = inv(U'*A);
% 		return
% 	end
% end

m = 2^n2;
B = zeros(m, n2);

% TODO: make it less ugly + efficient
r = 0:2:2*n2-2;
z = ones(1,n2);
z(1) = 0;
for i = 1:m
	z = nexti(z);
	B(i,:) = M(z+r);
end

W = zeros(n1,m);
for j = 1:n1
	q = [(M(2,:)-U(j,:))./L; (U(j,:)-M(1,:))./L];
	z = ones(1,n2);
	z(1) = 0;
	for i = 1:m
		z = nexti(z);
		W(j,i) = prod(q(z+r));
	end
end
%W = [U ones(size(U,1),1)]/[B ones(m,1)];

% next index
function z = nexti(z)
i = 1;
z(i) = z(i) + 1;
while z(i) > 2
	z(i) = 1;
	i = i + 1;
	z(i) = z(i) + 1;
end


text:
fclose all;
clear;

%I/O Files
DIRPATH='\Tight\sequential\data\';
FIMU=fopen([DIRPATH 'IMU.txt'],'r'); %%Assumed to be sampled at stricly 400Hz
FGPS=fopen([DIRPATH 'GPS.txt'],'r'); %%In order not to deal with program control stutucture it is assumed that GPS data is available at strictly 1Hz.
FCOV=fopen([DIRPATH 'Cov.txt'],'wt');
FOUT=fopen([DIRPATH 'Output.txt'],'wt');

%Load System parameters
[SenErrDef, ClkErrDef, IniErrDef, ObsErrDef, ini_pva, dt]=sys_def_v000();

%initial values
pos_e=ecef2geo_v000(ini_pva(:,1),1);
Cne=pos2Cne_v000(ini_pva(1,1),ini_pva(2,1));
Cbn=euler2dcm_v000(ini_pva(:,3));
Cbe=Cne*Cbn;
vel_e=Cne*ini_pva(:,2);

%Sensor error models
[Aimu, Bimu, Cimu, Dimu, sP0imu]=imu_modTI_v000(SenErrDef);
Qimu=Bimu*Bimu';
Rimu=Dimu*Dimu';
%Clock Error model
[Aclk, Bclk, Cclk, Dclk, sP0clk]=imu_modTI_v000(ClkErrDef);

%%In fact this part should be performed inside the propagate() routine.
%However, as I fixed the the propagation period, I simply discretize it with
%the known dt;
mx_a = 40*dt*[-Aclk,Bclk*Bclk';zeros(size(Aclk)),Aclk'];
mx_b = expm(mx_a);
STMclk = mx_b(3:4,3:4)';
Qclk = STMclk*mx_b(1:2,3:4);

%Initial covariance
nst=17;
P=zeros(nst);
P(1:3,1:3)=IniErrDef.pos_sP*IniErrDef.pos_sP';
P(4:6,4:6)=IniErrDef.vel_sP*IniErrDef.vel_sP';
P(7:9,7:9)=IniErrDef.att_sP*IniErrDef.att_sP';
P(10:15,10:15)=sP0imu*sP0imu';  %IMU errors
P(16:17,16:17)=sP0clk*sP0clk';  %Clock Errors

%Augmented states
Ximu=zeros(size(Aimu,1),1);
Xclk=zeros(size(Aclk,1),1);

pr_ctr=0;
imu_inc=zeros(6,1);
%Read the first IMU data
imu_data=fscanf(FIMU, '%f', 6);	%order=[Gyro;Acc]
while(~feof(FIMU))
    acc=imu_data(4:6)-Ximu(1:3);    %Correct them with error estimates
    gyro=imu_data(1:3)-Ximu(4:6);

    %Total IMU increments (used for low rate propagation)
    imu_inc=imu_inc+[acc;gyro]*dt;

    %INS calculations
    [Cbe, vel_e, pos_e]=strapdown_ecef_dcm_v000(Cbe, vel_e, pos_e, acc, gyro, dt);
    pr_ctr=pr_ctr+1;

    %Propagate the covariance and states at 10Hz
    if (~mod(pr_ctr,40))
        %Navigation model
        [STMNav QNav]=mdl_ecef_dcm_v000(Cbe, imu_inc(1:3)/0.1, Aimu, Qimu, Cimu, Rimu, dt*40);

        %Augment the clock model to nav model
        STM=diagmat_v000(STMclk,STMNav,1);
        Q=diagmat_v000(Qclk,QNav,1);

        %Propagate Covariance
        P=STM*P*STM'+Q;

        %Propagate IMU and Clock States
        Ximu=STMNav(10:15,10:15)*Ximu;
        Xclk=STMclk*Xclk;

        imu_inc=zeros(6,1);
    end

    %%Record the outputs
    if (~mod(pr_ctr,40))
        Llh=ecef2geo_v000(pos_e,0);
        Cne=pos2Cne_v000(Llh(1),Llh(2));
        vel_n=Cne'*vel_e;
        Cbn=Cne'*Cbe;
        eul=dcm2euler_v000(Cbn);
        fprintf(FOUT, '%12.12f ', [pr_ctr;Llh;vel_n;eul;Ximu;Xclk]);
        fprintf(FOUT, '\n');
        fprintf(FCOV, '%12.12f', [diag(P).^0.5]);
        fprintf(FCOV, '\n');
    end

    %Apply the GPS Pseudorange at 1Hz
    if (~mod(pr_ctr,400))
        gps_data=get_gps_data(FGPS);    %gps_data(i,:)=[sat_prn,pseudo distance,sat_x,sat_y,sat_z];

        %Start Sequential Processing of PseudoRanges
        dx=zeros(nst,1);
        for (in=1:size(gps_data,1))
            %Compute the innovation
            nom_range=norm(gps_data(in,3:5)'-pos_e);
            inno=gps_data(in,2)-(nom_range+Cclk*Xclk);

            %Observation Matrix
            H=zeros(1,nst);
            H(1,1:3)=(gps_data(in,3:5)-pos_e')/nom_range;
            H(1,16:17)=Cclk;

            Robs=ObsErrDef.sR*ObsErrDef.sR';

            %Update
            K=(P*H')/(H*P*H'+Robs);
            P=(eye(nst)-K*H)*P;
            dx=dx+K*(inno-H*dx);
        end

        %Correct states
        pos_e=pos_e-dx(1:3);
        vel_e=vel_e-dx(4:6);
        mx_a=eye(3)+skew(dx(7:9));
        Cbe=mx_a*Cbe;
        Ximu=Ximu+dx(10:15);
        Xclk=Xclk+dx(16:17);
    end

    %Read the next imu data from file
    imu_data=fscanf(FIMU, '%f', 6);
end

fclose(FIMU);
fclose(FGPS);
fclose(FCOV);
fclose(FOUT);

return;

%%plot the results
load([DIRPATH  'output.txt']);
plot(output(:,2),output(:,3));

repo: ekfmonoslam


text:
function a = getBestAction(Q, state, actions)

    % Select the best action a in state s

    n_actions = length(actions);
    status = find(Q(state,:), 1);

    if isempty(status)

        a = randi(n_actions);
    else

        [Q_value, a] = max(Q(state,:));
    end

end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
function labels = ind2label(obj, indices)
    % Returns the label of SymVariable specified by the position index
    %
    %
    % Parameters:
    %  indices: position indices of symbolic variable
    %
    % Return values:
    %  labels: a cellstr array of the labels
    %

    labels = obj.label(indices);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
sphere(30);
title('a sphere: x^2+y^2+z^2');
xlabel('x');
ylabel('y');
zlabel('z');
axis equal

author: thomas-koehler

repo: SupER


text:
function [vid_h_org, vid_l_bic] = load_video1(opt,part,param)

dirFile = opt.test_video.imgPath;
nFrames = opt.test_video.nFrame;
upscale = opt.res;

vid_h_org = [];
vid_l_bic = [];

for i = 1:nFrames
    filename = sprintf('original/Frame %03d.png', i);
    fullname = strcat(dirFile,filename);
    tmp = im2double(imread(fullname,'png'));
    [M, N] = size(tmp(:,:,1));
    switch part
        case 1 % whole frame
            vid_h_org{i} = tmp;
        case 2 % no use
            vid_h_org{i} = imresize(tmp, 1/2);
        case 3 % half x half (288x352)
            vid_h_org{i} = tmp(1:M/2,1:N/2,:);
        case 4 % quarter x quarter (144x176)
            vid_h_org{i} = tmp(1:M/4,1:N/4,:);
        case 5 % 64x64
            vid_h_org{i} = tmp(65:128,1:64,:);
        case 6 % 32x64
            vid_h_org{i} = tmp(65:96,1:64,:);
    end
    vid_l_bic{i} = imresize(vid_h_org{i}, 1/upscale);
end

if param.SHOW_IMAGE
    figure(1), imshow(vid_h_org{1}); title('original (frame1)');
end


text:
function inds = graphPeripheralVertices(v, e, l)
%GRAPHPERIPHERALVERTICES Peripheral vertices of a graph
%
%   INDS = graphPeripheralVertices(V, E)
%   Return indices of peripheral vertices, that is, vertices whose
%   eccentricity is maximal and equal to the diameter of the graph.
%
%   INDS = graphPeripheralVertices(V, E, L)
%   Specify length of each edge. Default is 1 for each edge.
%
%
%   Example
%     nodes = [20 20;20 50;20 80;50 50;80 20;80 50;80 80];
%     edges = [1 2;2 3;2 4;4 6;5 6;6 7];
%     figure; drawGraph(nodes, edges);
%     axis([0 100 0 100]); axis equal; hold on
%     INDP = graphPeripheralVertices(nodes, edges);
%     INDP =
%         1
%         3
%         5
%         7
%
%
%   See Also
%   grPropagateDistance, grVertexEccentricity
%   graphCenter, graphDiameter, graphRadius
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-09-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% ensure there is a valid length array
if nargin<3
    l = ones(size(e,1), 1);
end

g = grVertexEccentricity(v, e, l);

inds = find(g==max(g));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [trl] = read_fcdc_trl(fn)

% READ_FCDC_TRL reads trial definitions from a file
%
% Given a file which defines N trials, this function returns a Nx3
% matrix with the begin latency, end latency, and the latency offset
% of the first sample of each trial. The latencies are in seconds.
%
% [trl] = read_fcdc_trl(filename)
%
% An FCD trial definition file is formatted like
%   begin   end     offset
%   0.0000  1.0000  0.0000
%   3.0000  4.0000  0.0000
%   5.0000  5.5000  0.0000
%   ...
%
% The trial begin and end are given in seconds relative to the start
% of the recorded datafile. The offset is given in seconds and indicates
% the latency of the first sample, relative to the trial marker or
% trigger. E.g., given a trigger at 7000ms (relative to the recording
% begin), a trial of 1000ms with a pretrigger interval of 300ms would
% correspond to "6.700 7.700 -0.300".

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(fn, 'rt');

trl = [];
while ~feof(fid)
  tmp = fscanf(fid, '%f %f %f', 3);
  if ~isempty(tmp)
    trl = [trl; tmp'];
  end
end

fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Initialization
clear all; close all; clc;
gen = spx.data.synthetic.MultiToneSignalGenerator();
gen.TotalDuration = 0.35;
% Sampling frequency
gen.SamplingFrequency = 1000;
% Frequencies (Hz)
gen.Frequencies = [113 247 327 413];
% Each frequency will be present in a segment of the
% overall signal duration.
% The segment is described by its origin and its duration.
gen.Origins = [0 0 0.030 0.150];
gen.Durations = [0.350 0.050 0.200 0.200];
% Amplitudes for each frequency
gen.Amplitudes = [1 1.7 1.9 1.8];

[x , t] = gen.run();
plot(t, x);

author: indigits

repo: sparse-plex


text:
% EX_LIN_ELAST_CUBE_MP: solve linear elasticity problem in a multipatch cube.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_2cubesa.txt';
% You can see how the patches should match in the files
%  geo_2cubesa{b,c,d,e,f,g,h}.txt
% In all the eight cases the result must be the same

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2 3 4 5 6];

% Physical parameters
E  =  1; nu = 0.3;
problem_data.lambda_lame = @(x, y, z) ((nu*E)/((1+nu)*(1-2*nu)) * ones (size (x)));
problem_data.mu_lame = @(x, y, z) (E/(2*(1+nu)) * ones (size (x)));

% Source and boundary terms
fx = @(x, y, z) problem_data.mu_lame(x, y, z) .* (2*cos(x) - 4*x.*y.^2.*z) + ...
                problem_data.lambda_lame(x, y, z) .* (cos(x) - 4*x.*y.^2.*z);
fy = @(x, y, z) problem_data.mu_lame(x, y, z) .* (2*sin(y).*z - 4*x.^2.*y.*z) + ...
                problem_data.lambda_lame(x, y, z) .* (z.*sin(y) - 4*x.^2.*y.*z);
fz = @(x, y, z) problem_data.mu_lame(x, y, z) .* (-2*(y.*z).^2 - 2*(x.*z).^2 - cos(y) - 4*(x.*y).^2) - ...
                problem_data.lambda_lame(x, y, z) .* (2*(x.*y).^2 + cos(y));
problem_data.f = @(x, y, z) cat(1, ...
                    reshape (fx (x,y,z), [1, size(x)]), ...
                    reshape (fy (x,y,z), [1, size(x)]), ...
                    reshape (fz (x,y,z), [1, size(x)]));
problem_data.h = @(x, y, z, ind) cat(1, ...
                    reshape (cos(x), [1, size(x)]), ...
                    reshape (sin(y).*z, [1, size(x)]), ...
                    reshape ((x.*y.*z).^2, [1, size(x)]));

% Exact solution (optional)
uxex = @(x, y, z) cos(x);
uyex = @(x, y, z) sin(y).*z;
uzex = @(x, y, z) (x.*y.*z).^2;
problem_data.uex  = @(x, y, z) cat(1, ...
		      reshape (uxex (x,y,z), [1, size(x)]), ...
		      reshape (uyex (x,y,z), [1, size(x)]), ...
		      reshape (uzex (x,y,z), [1, size(x)]));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [2 2 2];     % Degree of the bsplines
method_data.regularity = [1 1 1];     % Regularity of the splines
method_data.nsub       = [2 2 2];     % Number of subdivisions
method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);

% 4) POST-PROCESSING.
% EXPORT TO PARAVIEW
output_file = 'lin_elast_cube_mp_Deg2_Reg1_Sub2';

vtk_pts = {linspace(0, 1, 10), linspace(0, 1, 10), linspace(0, 1, 10)};
fprintf ('results being saved in: %s.pvd\n \n', output_file)
sp_to_vtk (u, space, geometry, vtk_pts, output_file, {'displacement', 'stress'}, {'value', 'stress'}, ...
    problem_data.lambda_lame, problem_data.mu_lame)

% COMPARISON WITH THE EXACT SOLUTION
error_l2 = sp_l2_error (space, msh, u, problem_data.uex)

%!demo
%! ex_lin_elast_cube_mp

%!test
%! problem_data.geo_name = 'geo_2cubesa.txt';
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2 3 4 5 6];
%! E  =  1; nu = 0.3;
%! problem_data.lambda_lame = @(x, y, z) ((nu*E)/((1+nu)*(1-2*nu)) * ones (size (x)));
%! problem_data.mu_lame = @(x, y, z) (E/(2*(1+nu)) * ones (size (x)));
%! fx = @(x, y, z) problem_data.mu_lame(x, y, z) .* (2*cos(x) - 4*x.*y.^2.*z) + ...
%!                 problem_data.lambda_lame(x, y, z) .* (cos(x) - 4*x.*y.^2.*z);
%! fy = @(x, y, z) problem_data.mu_lame(x, y, z) .* (2*sin(y).*z - 4*x.^2.*y.*z) + ...
%!                 problem_data.lambda_lame(x, y, z) .* (z.*sin(y) - 4*x.^2.*y.*z);
%! fz = @(x, y, z) problem_data.mu_lame(x, y, z) .* (-2*(y.*z).^2 - 2*(x.*z).^2 - cos(y) - 4*(x.*y).^2) - ...
%!                 problem_data.lambda_lame(x, y, z) .* (2*(x.*y).^2 + cos(y));
%! problem_data.f = @(x, y, z) cat(1, ...
%!                     reshape (fx (x,y,z), [1, size(x)]), ...
%!                     reshape (fy (x,y,z), [1, size(x)]), ...
%!                     reshape (fz (x,y,z), [1, size(x)]));
%! problem_data.h = @(x, y, z, ind) cat(1, ...
%!                     reshape (cos(x), [1, size(x)]), ...
%!                     reshape (sin(y).*z, [1, size(x)]), ...
%!                     reshape ((x.*y.*z).^2, [1, size(x)]));
%! uxex = @(x, y, z) cos(x);
%! uyex = @(x, y, z) sin(y).*z;
%! uzex = @(x, y, z) (x.*y.*z).^2;
%! problem_data.uex  = @(x, y, z) cat(1, ...
%! 		      reshape (uxex (x,y,z), [1, size(x)]), ...
%! 		      reshape (uyex (x,y,z), [1, size(x)]), ...
%! 		      reshape (uzex (x,y,z), [1, size(x)]));
%! method_data.degree     = [2 2 2];     % Degree of the bsplines
%! method_data.regularity = [1 1 1];     % Regularity of the splines
%! method_data.nsub       = [2 2 2];     % Number of subdivisions
%! method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubesb.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubesc.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubesd.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubese.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubesf.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubesg.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)
%!
%! problem_data.geo_name = 'geo_2cubesh.txt';
%! [geometry, msh, space, u] = mp_solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (space.ndof, 504)
%! assert (error_l2, 2.42330537409875e-04, 1e-16)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ algs ] = gd_solver_list(category)
% Return list of solvers.
%
% Inputs:
%       category    category to be returned.
% Output:
%       algs        list of solvers in the category
%
% This file is part of GDLibrary.
%
% Created by H.Kasai on Nov. 02, 2016

    % supported algorithms by GDLibrary
    gd_algs = {'SD-STD','SD-BKT','SD-EXACT','SD-WOLFE','SD-SCALE-EXACT'};
    newton_algs = {'Newton-STD','Newton-DAMP','Newton-CHOLESKY'};
    cg_algs = {'CG-PRELIM','CG-BKT','CG-EXACT','CG-PRECON-EXACT'};
    ncg_algs = {'NCG-FR-BTK','NCG-FR-WOLFE','NCG-PR-BTK','NCG-PR-WOLFE'};
    bfgs_algs = {'BFGS-B-EXACT', 'BFGS-B-BKT','BFGS-H-EXACT','BFGS-H-BKT', ...
                    'DAMPED-BFGS-BKT','DAMPED-BFGS-EXACT', ...
                    'L-BFGS-BKT','L-BFGS-EXACT','L-BFGS-WOLFE'};
    other_algs = {'BB'};

    linesearch_algs = {'SD-BKT','SD-WOLFE','NCG-BTK','NCG-WOLFE','BFGS-B-BKT','BFGS-H-BKT','DAMPED-BFGS-BKT','L-BFGS-BKT','L-BFGS-WOLFE'};
    exact_algs = {'SD-EXACT','SD-SCALE-EXACT','CG-EXACT','CG-PRECON-EXACT','BFGS-B-EXACT','BFGS-H-EXACT','DAMPED-BFGS-EXACT','L-BFGS-EXACT'};

    switch category
        case 'SD'
            algs = gd_algs;
        case 'CG'
            algs = cg_algs;
        case 'NCG'
            algs = ncg_algs;
        case 'Newton'
            algs = newton_algs;
        case 'BFGS'
            algs = bfgs_algs;
        case 'LS'
            algs = linesearch_algs;
        case 'EXACT'
            algs = exact_algs;
        case 'ALL'
            algs = [gd_algs, newton_algs, cg_algs, bfgs_algs, other_algs];
        otherwise
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [ table ] = integerparttable(nfinal )
%Calculates a table of integer partions using numbers l.t. or e.t. k
%(or equivilently using k nonnegative integers which can be zero)
%Just calls with the largest value which has to recurrsively call all other
%possible values.
%
%this table states with first row "ways to partition 0" and first column
%"using only numbers <=0", this gives the one in the top left, these two
%rows are essentially a definition limit to finish the reccurance relation.
%Row k+1 column j+1 contains the ways to partition k
% using integers <= j (or equivilently j integers >= 0)
%
%to obtain the partitions of nfinal using numbers g.e. 2 and l.e. k do
%  [1;table(2:nfinal,k+1)-table(1:nfinal-1,k+1)]
%
%use of 32 bit integers means this is only good up to 2^32-1= 4294967295,
%for larger values change to 64 bit or dp
% Scales as ~ exp(pi*sqrt(2*n/3))/(4*n*sqrt(3)) for n>>1 and k>n

%Author: David Holdaway
%Last update: 27/04/2012

intpart(1,nfinal,1,1); %clears table

for n = 2:nfinal
    intpart(uint32(n),uint32(n));
end
table = intpart(1,nfinal,2); %gets table

for n = 1:nfinal
 table(n+1,n+1:nfinal+1) =  table(n+1,n+1);
end
end

function p = intpart(k,n,gettable,cleartable) %#ok<INUSD>
persistent table
if nargin == 4
    table = uint32(zeros(n+1)); %clears
%     table(1,:) = uint32(ones(1,1+n)); %n=0
%     table(2,2:n+1) = uint32(ones(1,n)); %n=1 k>1
%     table(:,2) = uint32(ones(n+1,1)); %k=0
    table(1,:) = uint32(1); %n=0
    table(2,2:n+1) = uint32(1); %n=1 k>1
    table(:,2) = uint32(1); %k=0
    p = 0; return
end
if nargin == 3
    p = table; return
end
k = uint32(k);
n = uint32(n);
if n < 0
    p = uint32(0); return
end
if k < 0
    p = uint32(0); return
end
if table(n+1,k+1) ~= 0
    p = table(n+1,k+1);  return
end

if n==0
    p = uint32(0); return
end

    table(n+1,k+1) = intpart(k,n-k) + intpart(k-1,n);
    p = table(n+1,k+1);
    for lp = n+2:length(table)
    table(n+1,lp) = table(n+1,n+1);
    end

end


text:
function newPrevPermute = removeDimPrevPermute(prevPermute, dimNumber)
    newPrevPermute = prevPermute;

    for iDim = 1:size(prevPermute, 2)
        dimensionNumber = prevPermute{1, iDim};
        if(dimensionNumber > dimNumber)
            newPrevPermute{1, iDim} = dimensionNumber - 1;
        elseif(dimensionNumber == dimNumber)
            removeIndex = iDim;
        end
    end
    newPrevPermute(:, removeIndex) = [];
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function aux_cs2(params, hParentFigure)
% function aux_FMD(params, hParentFigure);
%-------------------------------------------
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% Thomas van Stiphout, thomas@sed.ethz.ch
% last update: 7.9.2005

report_this_filefun(mfilename('fullpath'));

disp('here we go for a cross section');
view(0,90);
params2=params;
params.nDx=5;
params.nDy=5;
% ask for coordinates
% vSecLim1=[-153.60 61.10]
% vSecLim2=[-149.20 60.17]
% vSecLim1=[-153.00 62.00]
% vSecLim2=[-151.00 59.00]
vSecLim1=[-154.01 60.88]
vSecLim2=[-150.50 60.37]
% vSecLim1 = ginput(1);
% vSecLim2 = ginput(1);
dist_=deg2km(distance(vSecLim1(2),vSecLim1(1),vSecLim2(2),vSecLim2(1)));
azimuth_=azimuth(vSecLim1(2),vSecLim1(1),vSecLim2(2),vSecLim2(1));

% round length of cross section to multiple of nDx
[vSecLim2(2),vSecLim2(1)]=reckon('gc',vSecLim1(2),vSecLim1(1),km2deg(params.nDx*(round(dist_/params.nDx))),azimuth_);

vDist_=km2deg(0:params.nDx:params.nDx*(round(dist_/params.nDx)))';
vOnes_=ones(length(vDist_),1);
mPolygon(:,1)=vSecLim1(1)*ones(length(vDist_),1);
mPolygon(:,2)=vSecLim1(2)*ones(length(vDist_),1);
vAzimuth_=azimuth_*vOnes_;
[params.vY, params.vX]=reckon('gc',mPolygon(:,2),mPolygon(:,1),vDist_,vAzimuth_(:));
vDepth_=[40:5:120]';
params.mPolygon=reshape(ones(length(vDepth_),1)*params.vX',length(vDepth_)*length(mPolygon),1);
% vOrigin_=(max(params.mPolygon)+min(params.mPolygon))./2;
% mOrigin_=ones(size(params.mPolygon(:,1),1),1)*vOrigin_;
% distance(mOrigin_(:,2),params.mPolygon(:,1),mOrigin_(:,2),mOrigin_(:,1));
params.mPolygon(:,2)=reshape(ones(length(vDepth_),1)*params.vY',length(vDepth_)*length(mPolygon),1);
params.mPolygon(:,3)=reshape(vDepth_*ones(1,size(mPolygon,1)),size(mPolygon,1)*size(vDepth_,1),1);

parmas.mValueGrid=params.mPolygon(:,3);
params.mValueGrid(isnan(params.mPolygon(:,3)),:)=nan;
params.mX=reshape(params.mPolygon(:,1),length(vDepth_),size(params.mPolygon,1)/length(vDepth_));
params.mY=reshape(params.mPolygon(:,2),length(vDepth_),size(params.mPolygon,1)/length(vDepth_));
params.mZ=reshape(params.mPolygon(:,3),length(vDepth_),size(params.mPolygon,1)/length(vDepth_));
[Nx,Ny,Nz]=surfnorm(params.mX,params.mY,-params.mZ);
% vAzimuth2_=repmat(vAzimuth_,1,size(params.mX,1))';
% vAzimuth2_=repmat(vAzimuth_,size(params.mPolygon(:,1),1)/size(vAzimuth_,1),1);
% [Y_, X_]=reckon('gc',params.mY,params.mX,ones(size(vAzimuth2_))*km2deg(30),vAzimuth2_+90);
% Nx=params.mX-X_;
% Ny=params.mY-Y_;
params.mT=reshape(Nx,size(params.mPolygon,1),1);
params.mT(:,2)=reshape(Ny,size(params.mPolygon,1),1);
params.mT(:,3)=reshape(Nz,size(params.mPolygon,1),1);

% mT_1=reshape(params.mT(:,1),17,23)
% mT_2=reshape(params.mT(:,2),17,23)
% mT_3=reshape(params.mT(:,3),17,23)
% hold on;quiver3(params.mX(:),params.mY(:),-params.mZ(:),mT_1,mT_2,mT_3,2,'r')

params.mValueGrid=params.mPolygon(:,3);

params.vUsedNodes=ones(size(params.mPolygon,1),1);
Name=cellstr('Cylindrical Volume Sampling');
for ii=1:length(Name)
    if strcmp(Name(ii),'Cylindrical Volume Sampling')
        disp(Name(ii));
        params.sComment=strcat(Name(ii),' Cross-Section Vertical');
        [params.caNodeIndices2 params.vResolution2]=ex_CreateIndexCatalogCylinder(params.mCatalog,params.mPolygon,params.mT,...
            params.vUsedNodes,params.bCylSmpModeN,params.fCylSmpValue,params.fCylSmpBnd);
        if ((params.bCylSmpModeN==1) && (params.bCylSmpModeR==0))
            params.vcsGridNames(6) = cellstr(char('Resolution [R]'));
        else
            params.vcsGridNames(6) = cellstr(char('Resolution [N]'));
        end
    elseif strcmp(Name(ii),'3D Spherical Volume Sampling')
        disp(Name(ii));
        params.sComment=strcat(Name(ii),' Depth=',num2str(params.vPercentiles(nLayer)));
        % nGriddingMode = 3 means 3D sphere around polygonpoint
        params.nGriddingMode=3;
        % select
        % Create Indices to catalog and select quakes in time period
        [params.caNodeIndices2 params.vResolution2] = ex_CreateIndexCatalog3D(params.mCatalog, params.mPolygon, '1', params.n3DGriddingMode, ...
            params.f3DSmpValue, params.f3DSmpBnd, params.fSizeRectHorizontal, params.fSizeRectDepth);
        if ((params.b3DSmpModeN==1) && (params.b3DSmpModeR==0))
            params.vcsGridNames(6) = cellstr(char('Resolution [R]'));
        else
            params.vcsGridNames(6) = cellstr(char('Resolution [N]'));
        end
        % clear indices where not enough quakes for depth estimation
        for i=1:length(params.mPolygon)
            if isnan(params.mPolygon(i,3))
                params.caNodeIndices2{i}=[];
            end
        end
    else
        disp('Something is going wrong');
    end

% for testing sampling volumes
% hold on;i=300;plot3(params.mCatalog(params.caNodeIndices2{i},1),params.mCatalog(params.caNodeIndices2{i},2),-params.mCatalog(params.caNodeIndices2{i},7),'xr')

    % Calculation of bValue
    for i=1:length(params.mPolygon)
        if ~isempty(params.mCatalog(params.caNodeIndices2{i},:))
            % function [fBValue, fStdDev, fAValue, fMeanMag] =  calc_bvalue(mCatalog, fBinning)
            [params.mValueGrid(i,2),params.mValueGrid(i,3),params.mValueGrid(i,4),params.mValueGrid(i,5),] =  calc_bvalue(params.mCatalog(params.caNodeIndices2{i},:));
            params.mValueGrid(i,6)=max(params.vResolution2{i});
            % set values to NaN, when radius of sampling volume increase limit
            if (params.mValueGrid(i,6) > 20) %params.fCylSmp)
                params.mValueGrid(i,2:5)=nan;
            end
        else
            % create NaN's in mValueGrid, where strike is not defined
            params.mValueGrid(i,:)=nan;
        end

    end

    params.vcsGridNames(1:6) = cellstr(char(...
        'Depth Level',...               %   1
        'b-Value',...                   %   2
        'Standard Deviation',...        %   3
        'A-Value',...                   %   4
        'Mean Magnitude',...            %   5
        'Resolution [km]'));            %   6

    params.bMap=2;
    params.mValueGrid(isnan(params.mValueGrid(:,1)),:)=nan;
    params.mPlotValues=reshape(params.mValueGrid(:,2),size(params.mX,1),size(params.mX,2));
    hold on; surf(params.mX,params.mY,-params.mZ,params.mPlotValues);
    shading interp;
%     hold on;i=300;plot3(params.mCatalog(params.caNodeIndices2{i},1),params.mCatalog(params.caNodeIndices2{i},2),-params.mCatalog(params.caNodeIndices2{i},7),'k');
%     vResults((nLayer-1)*length(Name)+ii)=params;
end
% mPolygon1=ones(size(mPolygon(:,1),1),1)*mPolygon(:,1)';
% mPolygon2=(ones(size(mPolygon(:,1),1),1)*mPolygon(:,2)')';
% vPolygon1=reshape(mPolygon1,length(mPolygon1)*length(mPolygon2),1);
% vPolygon2=reshape(mPolygon2,length(mPolygon1)*length(mPolygon2),1);
% params.mZ_=(ones(size(vPolygon1,1),1))*vDepth_';
% vDepth_=reshape(params.mZ_,size(params.mZ_,1)*size(params.mZ_,2),1)
% mPoly1=ones(size(vDepth_,1)/size(vPolygon1,1),1)*vPolygon1';
% mPoly1=reshape(mPoly1',length(vDepth_),1);
% mPoly2=ones(size(vDepth_,1)/size(vPolygon2,1),1)*vPolygon2';
% mPoly2=reshape(mPoly2',length(vDepth_),1);

view(3);

author: CelsoReyes


text:
function jed = transition_to_jed_mayan_long ( )

%*****************************************************************************80
%
%% TRANSITION_TO_JED_MAYAN_LONG: Mayan long count calendar transition as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the date.
%
  jed = 2456282.5;

  return
end


text: any(cellfun(@exist,{'arraydeal', 'chopdeal'})==3)

repo: REST


text:
function value = p07_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P07_F evaluates the integrand for problem 07.
%
%  Dimension:
%
%    N arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    product ( pi / 2 ) * sin ( pi * x(1:dim_num) )
%
%  Exact Integral:
%
%    1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 0.0;

  for point = 1 : point_num
    value(point) = prod ( 0.5 * pi * sin ( pi * x(1:dim_num,point) ) );
  end

  p07_i4 ( 'I', '#', point_num );

  return
end


text:
classdef PostProcessDataBaseCreator < handle

    properties (Access = private)
        mesh
        outFileName
    end

    methods (Access = public)

        function obj = PostProcessDataBaseCreator(dI)
            obj.init(dI);
        end

        function data = create(obj)
            d.outFileName = obj.outFileName;
            d.coordinates = obj.mesh.coord;
            d.connectivities = obj.mesh.connec;
            d.nnode = size(obj.mesh.connec,2);
            d.npnod = size(obj.mesh.coord,1);
            d.gtype = obj.mesh.type;
            d.nelem = size(obj.mesh.connec,1);
            d.etype = obj.computeGiDElementType(d.gtype);
            data = d;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            if isfield(cParams,'outName')
                cParams.outFileName = cParams.outName;
            end
            obj.mesh        = cParams.mesh;
            obj.outFileName = cParams.outFileName;
        end

    end

    methods (Access = private, Static)

        function et = computeGiDElementType(gt)
            switch  gt
                case 'TRIANGLE'
                    et = 'Triangle';
                case 'QUAD'
                    et = 'Quadrilateral';
                case 'TETRAHEDRA'
                    et = 'Tetrahedra';
                case 'HEXAHEDRA'
                    et = 'Hexahedra';
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [] = selectStateSeqMinExpectedHamDist( jobID, taskIDs, MIN_ITER )
% Take all burned-in samples from multiple chains of MCMC
%   divide into Reference (3/4) and Test (1/4) sets
% Return the stateSeq from the Test set
%   that minimizes the expected hamming distance to all Reference samples

% jobID = 1121566;
% taskIDs = 1:5;
% MIN_ITER = 15000;
REF_FRAC = 9/10;

% ==================================================== Load Stored Samples
aa = 0;
for tt = 1:length( taskIDs )
    taskID = taskIDs(tt);
    X = loadSamplerOutput( jobID, taskID );
    for ss = 1:length( X.S )
       if X.iters.S(ss) <= MIN_ITER
           continue;
       end
       aa = aa+1;
       AllSamps(aa).stateSeq = X.S(ss).stateSeq;
       AllSamps(aa).taskID = taskID;
       AllSamps(aa).jobID  = jobID;
       AllSamps(aa).sampID = ss;
       AllSamps(aa).sampIter = X.iters.S(ss);
    end
end

nSamps = length( AllSamps );
permIDs = randperm( nSamps );
R = round( REF_FRAC * nSamps );
RefSamps = AllSamps( permIDs(1:R) );
TestSamps = AllSamps( permIDs(R+1:end)  );

fprintf( 'Searching for one test sample that best aligns to all reference samples\n');
fprintf( '  |Ref Set|=%d.  |Test Set|=%d\n', R, length(TestSamps)   );

MeanHamDist = zeros( 1, length( TestSamps ) );
tic;
for tt = 1:length( TestSamps )
    HDistToRefs = zeros( 1, R );
    for rr = 1:length( RefSamps )
        HDistToRefs(rr) = calcHammingDistance( RefSamps(rr).stateSeq,  TestSamps(tt).stateSeq );
    end
    MeanHamDist(tt) = mean( HDistToRefs );
    if tt==1 || mod( tt, 25 ) == 0 || tt == length(TestSamps)
        fprintf( ' %5.0f sec. %4d/%d test sequences examined\n', toc, tt, length(TestSamps) );
    end
end

[~, bestID] = min( MeanHamDist );

BestSample = TestSamps( bestID );
X = loadSamplerOutput( BestSample.jobID, BestSample.taskID );
BestSample.A = X.A( BestSample.sampID );

fid = fopen( '~/git/liv-video/SimulationResults.path' );
RESULTS_DIR = textscan( fid, '%s' );
RESULTS_DIR =RESULTS_DIR{1}{1};
fclose(fid);
outpath = fullfile( RESULTS_DIR, num2str(jobID), 'MinExpectedHamDistSample.mat' );
save( outpath, '-struct', 'BestSample' );
fprintf( '... wrote best sample to file %s\n', outpath );

author: michaelchughes

repo: NPBayesHMM


text:
function [Population,Range] = EnvironmentalSelection(Population,RefPoint,Range,N)
% The environmental selection of AR-MOEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    CV = sum(max(0,Population.cons),2);
    if sum(CV==0) > N
        %% Selection among feasible solutions
        Population = Population(CV==0);
        % Non-dominated sorting
        [FrontNo,MaxFNo] = NDSort(Population.objs,N);
        Next = FrontNo < MaxFNo;
        % Select the solutions in the last front
        Last   = find(FrontNo==MaxFNo);
        Choose = LastSelection(Population(Last).objs,RefPoint,Range,N-sum(Next));
        Next(Last(Choose)) = true;
        Population = Population(Next);
        % Update the range for normalization
        Range(2,:) = max(Population.objs,[],1);
        Range(2,Range(2,:)-Range(1,:)<1e-6) = 1;
    else
        %% Selection including infeasible solutions
        [~,rank]   = sort(CV);
        Population = Population(rank(1:N));
    end
end

function Remain = LastSelection(PopObj,RefPoint,Range,K)
% Select part of the solutions in the last front

    N  = size(PopObj,1);
    NR = size(RefPoint,1);

    %% Calculate the distance between each solution and point
    Distance    = CalDistance(PopObj-repmat(Range(1,:),N,1),RefPoint);
    Convergence = min(Distance,[],2);

    %% Delete the solution which has the smallest metric contribution one by one
    [dis,rank] = sort(Distance,1);
    Remain     = true(1,N);
    while sum(Remain) > K
        % Calculate the fitness of noncontributing solutions
        Noncontributing = Remain;
        Noncontributing(rank(1,:)) = false;
        METRIC = sum(dis(1,:)) + sum(Convergence(Noncontributing));
        Metric = inf(1,N);
        Metric(Noncontributing) = METRIC - Convergence(Noncontributing);
        % Calculate the fitness of contributing solutions
        for p = find(Remain & ~Noncontributing)
            temp = rank(1,:) == p;
            noncontributing = false(1,N);
            noncontributing(rank(2,temp)) = true;
            noncontributing = noncontributing & Noncontributing;
            Metric(p) = METRIC - sum(dis(1,temp)) + sum(dis(2,temp)) - sum(Convergence(noncontributing));
        end
        % Delete the worst solution and update the variables
        [~,del] = min(Metric);
        temp    = rank ~= del;
        dis     = reshape(dis(temp),sum(Remain)-1,NR);
        rank    = reshape(rank(temp),sum(Remain)-1,NR);
        Remain(del) = false;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function contrib_prev = reverse_conv(l_prev, contrib_curr, weights)

NON_NEG = true;
UNIT_LAYER = true;

new_weights = deconv_weights(weights);

    if NON_NEG && UNIT_LAYER
        new_weights(new_weights<0) = 0;
        l_prev(l_prev>0)=1;
        contrib_prev = my_vl_conv(contrib_curr, new_weights, []);
        contrib_prev = gather(gpuArray(contrib_prev) .* gpuArray(l_prev));
        contrib_prev = contrib_prev / max(abs(contrib_prev(:)));
    end

    if ~NON_NEG
        % For all positive elements in contrib_curr
        % positive weights means positive contribution
        % negative weights means negative contribution
        % so just do a normal conv pass with both positve and negative weights
        contrib_curr_pos = contrib_curr;
        contrib_curr_pos(contrib_curr_pos < 0) = 0;
        contrib_prev_pos = my_vl_conv(contrib_curr_pos, new_weights, []);
        ontrib_prev_pos = contrib_prev_pos .* l_prev;

        % For all non positive elements in contrib_curr
        % positive weights means negative contribution
        % negative weights means nothing (something that doesn't contribute to a negative pattern, can be positive or negative)
        % so we get the positive weights only to do the conv pass
        contrib_curr_neg = contrib_curr;
        contrib_curr_neg(contrib_curr_neg > 0) = 0;
        new_weights_pos = new_weights;
        new_weights_pos(new_weights_pos < 0) = 0;
        contrib_prev_neg = my_vl_conv(contrib_curr_neg, new_weights_pos, []);
        contrib_prev_neg = contrib_prev_neg .* l_prev;

        % Sum up the pos and neg results
        contrib_prev = contrib_prev_pos - sqrt((-1)*contrib_prev_neg);
        contrib_prev = contrib_prev / max(contrib_prev(:));
    end

end

    %{
    %{
    we can be sure that all l_prev is positive, and we don't need negative
    tensor flows to normalize the positive ones
    %}
    weights(weights<0) = 0;

    %normalize contrib_curr with its original input to get percentage
    %contribution

    %forward_pass = forward_conv(l_prev,weights);

    %to handle 0/0 case, we set all non-positive forward_pass element as -1
    %rationale being that we had removed all negative tensor flows by doing
    %weights(weights<0)=0, if the result is 0, then it must be originally
    %also 0 (because originally we didn't remove the negative flows so the
    %result will be negative, which is 0 after relu), so 0/0 will happen in
    %the element-wise division step below. To prevent this, set
    %non-positive elements in forward_pass to -1 so that it becomes 0/-1.
    %The answer is 0 which is the desired behavior (because 0% contribution
    %is what it deserves)
    % that way

    %forward_pass(forward_pass<=0) = -1;
    %contrib_curr = contrib_curr ./ forward_pass;

    %}
    %contrib_curr = contrib_curr .* sqrt(contrib_curr);
    %contrib_curr = contrib_curr / sum(contrib_curr(:));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Jiao Xianjun (putaoshu@msn.com; putaoshu@gmail.com)
% parameter setting via tcp to rtl_tcp
% A script of project: https://github.com/JiaoXianjun/multi-rtl-sdr-calibration

function tcp_obj = set_gain_tcp(tcp_obj, gain)

if gain
    fwrite(tcp_obj, 3, 'uint8');
    fwrite(tcp_obj, uint32(1), 'uint32');
    fwrite(tcp_obj, 4, 'uint8');
    fwrite(tcp_obj, uint32(gain), 'uint32');
else
    fwrite(tcp_obj, 3, 'uint8');
    fwrite(tcp_obj, uint32(0), 'uint32');
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function z = lambda_sum_largest( Y, k )

% LAMBDA_SUM_SMALLEST    Sum of the k smallest eigenvalues of a symmetric matrix.
%     For square matrix X, LAMBDA_SUM_LARGEST(X,K) is SUM_LARGEST(EIG(X),k)
%     if X is Hermitian or symmetric and real; and +Inf otherwise.
%
%     An error results if X is not a square matrix.
%
%     Disciplined convex programming information:
%         LAMBDA_SUM_LARGEST is convex and nonmonotonic (at least with
%         respect to elementwise comparison), so its argument must be affine.

error( nargchk( 2, 2, nargin ) );
if ndims( Y ) > 2 || size( Y, 1 ) ~= size( Y, 2 ),
    error( 'First input must be a square matrix.' );
elseif ~isnumeric( k ) || numel( k ) ~= 1 || ~isreal( k ),
    error( 'Second input must be a real scalar.' );
end
err = Y - Y';
Y   = 0.5 * ( Y + Y' );
if norm( err, 'fro' )  > 8 * eps * norm( Y, 'fro' ),
    z = Inf;
else
    z = sum_largest( eig( full( Y ) ), k );
end

% 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.


text:
 Dataflow {
     // This is an Eyeriss-like row-stationary dataflow”
    SpatialMap(1,1) Y';
    TemporalMap(1,1) X';
    TemporalMap(1,1) C;
    TemporalMap(16,16) K;
    TemporalMap(Sz(R),Sz(R)) R;
    TemporalMap(Sz(S),Sz(S)) S;
    Cluster(Sz(R),P);
    SpatialMap(1,1) Y;
    SpatialMap(1,1) R;
    TemporalMap(Sz(S),Sz(S)) S;
}

author: maestro-project

repo: maestro


text:
function [percent progressLength] = parfor_progress(N)
%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor.
%   PARFOR_PROGRESS works by creating a file called parfor_progress.txt in
%   your working directory, and then keeping track of the parfor loop's
%   progress within that file. This workaround is necessary because parfor
%   workers cannot communicate with one another so there is no simple way
%   to know which iterations have finished and which haven't.
%
%   PARFOR_PROGRESS(N) initializes the progress monitor for a set of N
%   upcoming calculations.
%
%   PARFOR_PROGRESS updates the progress inside your parfor loop and
%   displays an updated progress bar.
%
%   PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress
%   bar.
%
%   To suppress output from any of these functions, just ask for a return
%   variable from the function calls, like PERCENT = PARFOR_PROGRESS which
%   returns the percentage of completion.
%
%   Example:
%
%      N = 100;
%      parfor_progress(N);
%      parfor i=1:N
%         pause(rand); % Replace with real code
%         parfor_progress;
%      end
%      parfor_progress(0);
%
%   See also PARFOR.

% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/

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

if nargin < 1
    N = -1;
end

percent = 0;
progressLength = 0;

if usejava('desktop')==0
    return;
end
return;
w = 50; % Width of progress bar
saveDir = 'private';
filePath = [saveDir filesep 'parfor_progress.txt'];
if ~exist(saveDir,'dir');mkdir(saveDir);end

if N > 0
    f = fopen(filePath, 'w');
    if f<0
        error('Do you have write permissions for %s?', pwd);
    end
    fprintf(f, '%d\n', N); % Save N at the top of progress.txt
    fclose(f);

    if nargout == 0
        % disp(['  0%[>', repmat(' ', 1, w), ']']);
        dispstat(['  0%[>', repmat(' ', 1, w), ']']);
    end
elseif N == 0
    delete(filePath);
    percent = 100;

    if nargout == 0
        % disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']);
        dispstat([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']);
    end
else
    % if ~exist(filePath, 'file')
    %     error([filePath ' not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize ' filePath '.']);
    % end
    return;
    try
        f = fopen(filePath, 'A');
    catch
        error([filePath ' not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize ' filePath '.']);
    end
    fprintf(f, '1\n');
    fclose(f);

    f = fopen(filePath, 'r');
    progress = fscanf(f, '%d');
    progressLength = length(progress(2:end));
    fclose(f);
    try
        percent = (length(progress)-1)/progress(1)*100;
    catch
    end

    if nargout == 0
        perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage
        % reverseStr = repmat(sprintf('\b'), 1, length(perc)-1);
        % fprintf([reverseStr, perc]);drawnow;
        outputStr = [repmat(char(8), 1, (w+9)), char(10), perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']'];
        % disp(outputStr);
        dispstat(outputStr);
        % fprintf([repmat('\b', 1, (w+9)), '\n', perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']);
    end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [theta,alpha,minError,vertices3d,vertices2d,inliers]= ...
    select_from_range_fig(hObject, eventdata, handles)
% -------------------------------------------------------------------------
% FUNCTION: select_from_range_fig
% INPUT:
%  hObject : handle to current figure
%  eventdata : from calling function
%  handles : data common to GUI
% OUTPUT:
%  theta,alpha: parameters of selected plane as per equation
%               theta'*x=alpha , alpha>0 , norm(theta)=1
%  minError   : min fitting error using tls_robust routine
%  vertices3d : 3d coords of points corresponding to vertices of selected
%               polgon
%  inliers    : 3d coords of points considered inliers to the selected
%               plane
% -------------------------------------------------------------------------

%* Author: Ranjith Unnikrishnan                                          *
%* Carnegie Mellon University, Vision and Mobile Robotics Laboratory       *
%* THE MATERIAL EMBODIED IN THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"     *
%* AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,      *
%* INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR      *
%* FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL CARNEGIE MELLON  *
%* UNIVERSITY BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,            *
%* SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY         *
%* KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,        *
%* LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF     *
%* THIRD PARTIES, WHETHER OR NOT CARNEGIE MELLON UNIVERSITY HAS BEEN     *
%* ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON        *
%* ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE     *
%* POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.                      *
%

% Restrict to data with far > y > near
scan=handles.current_scan;

% Prune based on distance
dist_sq = scan(:,1).*scan(:,1)+scan(:,2).*scan(:,2);
scan(dist_sq < handles.nearVal*handles.nearVal ...
    | dist_sq > handles.farVal*handles.farVal,:) = [];

if (isempty(scan))
    return;
end

% Project scan to viewing vectors (up,right,forward)
P = [handles.right_vector handles.forward_vector handles.up_vector];
scanp = scan * P;
% Compute angle swept by data in xy plane
theta = atan2(scanp(:,2),scanp(:,1));

% Assume the angle swept by the data is less than 180 degrees
% Find angle corresponding to vector pointing at the middle of the cloud
midtheta = 0.5*(min(theta) + max(theta));
if (max(theta)-min(theta) > pi)
    if (midtheta > 0)
        midtheta = midtheta - pi;
    else
        midtheta = midtheta + pi;
    end
end

% Rotate the data by pi/2 - midtheta to put everything in view
R = [ cos(pi/2 - midtheta) sin(pi/2 - midtheta) 0;...
    -sin(pi/2 - midtheta) cos(pi/2 - midtheta) 0;...
    0 0 1];
scanp = scanp * R;

% Perspective projection (x/y, z/y)
scanp_y = scanp(:,2);
scanp=[scanp(:,1)./scanp(:,2) scanp(:,3)./scanp(:,2)];
%scanp(:,1) = scanp(:,1)/scanp(:,2);
%scanp(:,3) = scanp(:,3)/scanp(:,2);

% X,Y bounds of projection
xBounds=[ min(scanp(:,1)) max(scanp(:,1))];
yBounds=[ min(scanp(:,2)) max(scanp(:,2))];

% Size of a square pixel
pixel_size=max([range(xBounds) range(yBounds)])/150;

% Choose size of grid that will serve as the range image (n_cols n_rows)
nx=ceil(range(xBounds)/pixel_size) + 1;
ny=ceil(range(yBounds)/pixel_size) + 1;
grid=zeros(ny,nx);

% Add points to grid
ix=floor( (scanp(:,1)-xBounds(1))/pixel_size ) + 1;
iy=floor( (yBounds(2)-scanp(:,2))/pixel_size ) + 1;
ii=sub2ind(size(grid),iy,ix);
grid(ii)=max(grid(ii),scanp_y);

grid=-grid;
% Set values of zero to minimum value
grid(grid==0) = min(grid(:));
colormap(jet);
axes(handles.range_image_axis);
imagesc(grid);
axis('equal'); axis off;
title(handles.current_scanfilename);
%colorbar

% Draw a polygon and get its vertices
[x_polygon,y_polygon]=selectPolygon;
vertices2d=[x_polygon y_polygon];

% Choose points whose projections in range image
% are inside the polygon
region=scan(inpolygon(ix,iy,x_polygon,y_polygon),:);

fprintf(1,'Selected %d points\n',size(region,1));
% Fit plane to points
[theta,alpha,minError]=tls_robust(region);

% Compute the 3d points corresponding to the vertices of the polygon
vertices3d=[];
for i=1:4
    temp=[ix iy]-repmat([x_polygon(i) y_polygon(i)],size(ix,1),1);
    [m,j]=min(sum(temp.^2,2));
    %[m,j]=min(sum((scanp-repmat([x_polygon(i) y_polygon(i)],size(scanp,1),1)).^2,2));
    vertices3d=[vertices3d; scan(j,:)];
end

% Compute inliers = points within median distance from plane
dist=abs(theta'*region'-alpha);
inliers=region(dist<median(dist),:);

return;

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
function [xn, xmean, xstd] = normdata(x,xmean,xstd)
%NORMDATA  Normalize input to zero mean and unit variance
%
%  Description
%    [XN, XMEAN, XSTD] = NORMDATA(X) normalizes X to zero mean and
%    unit variance. Returns normalized XN, mean of X XMEAN and
%    standard deviance of X XSTD.
%
%    XN = NORMDATA(X,XMEAN,XSTD) normalizes X using precomputed
%    XMEAN and XSTD.
%
%  See also DENORMDATA
%
% Copyright (c) 2010 Aki Vehtari

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

  if nargin<=1
    xmean=nanmean(x);
    xstd=nanstd(x);
  end
  xn=bsxfun(@rdivide,bsxfun(@minus,x,xmean),xstd);

author: gpstuff-dev

repo: gpstuff

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


text:
function check_point = terminate(xi_1,xi,pi_1,pi)

temp1 = max(abs(xi_1(1,:)-xi(1,:))) + max(abs(pi_1(1,:)-pi(1,:)));
temp2 = max(abs(xi_1(2,:)-xi(2,:))) + max(abs(pi_1(2,:)-pi(2,:)));

check_point  = temp1+temp2;

author: Sable

repo: mcbench-benchmarks


text:
% EX_STOKES_ANNULUS_TH: solve the Stokes problem in one quarter of an annulus with generalized Taylor-Hood elements.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'annulus.mat';

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:4;
problem_data.nmnn_sides = [];

% Physical parameters
problem_data.viscosity = @(x, y) ones (size (x));

% Force term
fx = @(x, y) (4*(20*(98-57*x.^2).*y.^7+42*x.*(34*x.^2-75).*y.^6 ...
                 -3*(430*x.^4-1480*x.^2+1023).*y.^5 ...
                 +15*x.*(70*x.^4-310*x.^2+297).*y.^4 ...
                 +3*x.*(x.^4-5*x.^2+4).^2 ...
                 +2*(-290*x.^6+1500*x.^4-2079*x.^2+680).*y.^3 ...
                 +6*x.*(38*x.^6-255*x.^4+495*x.^2-260).*y.^2 ...
                 +(-75*x.^8+520*x.^6-1089*x.^4+720*x.^2-112).*y ...
                 +603*x.*y.^8-355*y.^9)-y./(x.^2+y.^2));

fy = @(x, y) (x./(x.^2+y.^2)+4*(5*x.^9-27*x.^8.*y+20*x.^7.*(15*y.^2-2) ...
                                +14*x.^6.*y.*(15-38*y.^2)+3*x.^5.*(290*y.^4-520*y.^2+33) ...
                                -15*x.^4.*y.*(70*y.^4-170*y.^2+33) ...
                                +x.^3.*(860*y.^6-3000*y.^4+2178*y.^2-80) ...
                                +18*x.^2.*y.*(-34*y.^6+155*y.^4-165*y.^2+20) ...
                                +x.*(285*y.^8-1480*y.^6+2079*y.^4-720*y.^2+16) ...
                                -67*y.^9+450*y.^7-891*y.^5+520*y.^3-48*y));
problem_data.f  = @(x, y) cat(1, ...
                reshape (fx (x,y), [1, size(x)]), ...
                reshape (fy (x,y), [1, size(x)]));

% Boundary terms
problem_data.h  = @(x, y, iside) zeros ([2, size(x)]); %Dirichlet

% Exact solution, to compute the errors
uxex = @(x, y) (2*(x-y).*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
                (x.^5-8*y-2*x.^4.*y+20*y.^3-6*y.^5+2*x.^2.*y.*(5-4*y.^2)+x.^3.*(-5 + 6*y.^2)+ ...
                 x.*(4+5*y.^2.*(-3+y.^2))));

uyex = @(x, y) (-2*(x - y).*y.^2.*(-4 + x.^2 + y.^2).*(-1 + x.^2 + y.^2).* ...
                (4 + 5*x.^2.*(-3 + x.^2) + 2*x.*(5 - 2*x.^2).*y + ...
                 (-5 + 6*x.^2).*y.^2 - 4*x.*y.^3 + y.^4));
problem_data.velex = @(x, y) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

problem_data.gradvelex = @test_stokes_annulus_graduex;

problem_data.pressex = @(x, y) (-(pi/8) + atan (y./x));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.element_name = 'th';     % Element type for discretization
method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
method_data.nsub         = [10 10];  % Number of subdivisions
method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = ...
                       solve_stokes (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) COMPARISON WITH EXACT SOLUTION
error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex)
[error_h1_v, error_l2_v] = ...
   sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex)

% 4.2) EXPORT TO PARAVIEW
output_file = 'ANNULUS_TH_Deg3_Reg2_Sub10';

fprintf ('The result is saved in the files %s \n and %s \n \n', ...
           [output_file '_vel'], [output_file '_press']);
vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
sp_to_vtk (press, space_p, geometry, vtk_pts, [output_file '_press'], 'press')
sp_to_vtk (vel,   space_v, geometry, vtk_pts, [output_file '_vel'  ], {'velocity', 'divergence'}, {'value', 'divergence'})

% 4.3) PLOT IN MATLAB
[eu, F] = sp_eval (vel, space_v, geometry, vtk_pts);
[X,  Y] = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));

figure()
subplot (1,2,2)
eu2 = problem_data.velex (X, Y);
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
axis equal
title('Exact solution')
subplot (1,2,1)
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
axis equal
title('Computed solution')

[div, F] = sp_eval (vel, space_v, geometry, vtk_pts, 'divergence');
figure()
surf (X, Y, div)
view(2)
axis equal
title('Computed divergence')
colorbar

%!demo
%! ex_stokes_annulus_th

%!test
%! problem_data.geo_name = 'annulus.mat';
%! problem_data.drchlt_sides = 1:4;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y) ones (size (x));
%! fx = @(x, y) (4*(20*(98-57*x.^2).*y.^7+42*x.*(34*x.^2-75).*y.^6 ...
%!                  -3*(430*x.^4-1480*x.^2+1023).*y.^5 ...
%!                  +15*x.*(70*x.^4-310*x.^2+297).*y.^4 ...
%!                  +3*x.*(x.^4-5*x.^2+4).^2 ...
%!                  +2*(-290*x.^6+1500*x.^4-2079*x.^2+680).*y.^3 ...
%!                  +6*x.*(38*x.^6-255*x.^4+495*x.^2-260).*y.^2 ...
%!                  +(-75*x.^8+520*x.^6-1089*x.^4+720*x.^2-112).*y ...
%!                  +603*x.*y.^8-355*y.^9)-y./(x.^2+y.^2));
%! fy = @(x, y) (x./(x.^2+y.^2)+4*(5*x.^9-27*x.^8.*y+20*x.^7.*(15*y.^2-2) ...
%!                                 +14*x.^6.*y.*(15-38*y.^2)+3*x.^5.*(290*y.^4-520*y.^2+33) ...
%!                                 -15*x.^4.*y.*(70*y.^4-170*y.^2+33) ...
%!                                 +x.^3.*(860*y.^6-3000*y.^4+2178*y.^2-80) ...
%!                                 +18*x.^2.*y.*(-34*y.^6+155*y.^4-165*y.^2+20) ...
%!                                 +x.*(285*y.^8-1480*y.^6+2079*y.^4-720*y.^2+16) ...
%!                                 -67*y.^9+450*y.^7-891*y.^5+520*y.^3-48*y));
%! problem_data.f  = @(x, y) cat(1, ...
%!                 reshape (fx (x,y), [1, size(x)]), ...
%!                 reshape (fy (x,y), [1, size(x)]));
%! problem_data.h  = @(x, y, iside) zeros ([2, size(x)]); %Dirichlet
%! uxex = @(x, y) (2*(x-y).*y.*(-4+x.^2+y.^2).*(-1+x.^2+y.^2).* ...
%!                 (x.^5-8*y-2*x.^4.*y+20*y.^3-6*y.^5+2*x.^2.*y.*(5-4*y.^2)+x.^3.*(-5 + 6*y.^2)+ ...
%!                  x.*(4+5*y.^2.*(-3+y.^2))));
%! uyex = @(x, y) (-2*(x - y).*y.^2.*(-4 + x.^2 + y.^2).*(-1 + x.^2 + y.^2).* ...
%!                 (4 + 5*x.^2.*(-3 + x.^2) + 2*x.*(5 - 2*x.^2).*y + ...
%!                  (-5 + 6*x.^2).*y.^2 - 4*x.*y.^3 + y.^4));
%! problem_data.velex = @(x, y) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.gradvelex = @test_stokes_annulus_graduex;
%! problem_data.pressex = @(x, y) (-(pi/8) + atan (y./x));
%! method_data.element_name = 'th';     % Element type for discretization
%! method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
%! method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
%! method_data.nsub         = [10 10];  % Number of subdivisions
%! method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule
%! [geometry, msh, space_v, vel, space_p, press] = ...
%!                        solve_stokes (problem_data, method_data);
%! error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex);
%! [error_h1_v, error_l2_v] = ...
%!    sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex);
%! assert (msh.nel, 100)
%! assert (space_v.ndof, 1058)
%! assert (space_p.ndof, 169)
%! assert (error_l2_p, 2.79382942386329e-05, 2e-14)
%! assert (error_h1_v, 0.0145880437598918, 1e-13)
%! assert (error_l2_v, 2.05350241308034e-04, 1e-14)


text:
function complete_symmetric_poly_test ( )

%*****************************************************************************80
%
%% COMPLETE_SYMMETRIC_POLY_TEST tests COMPLETE_SYMMETRIC_POLY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMPLETE_SYMMETRIC_POLY_TEST\n' );
  fprintf ( 1, '  COMPLETE_SYMMETRIC_POLY evaluates a complete symmetric.\n' );
  fprintf ( 1, '  polynomial in a given set of variables X.\n' );

  n = 5;
  x = [ 1.0, 2.0, 3.0, 4.0, 5.0 ];
  r8vec_print ( n, x, '  Variable vector X:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   N\\R     0       1       2       3       4       5\n' );
  fprintf ( 1, '\n' );

  for nn = 0 : n
    fprintf ( 1, '  %2d', nn );
    for rr = 0 : 5
      value = complete_symmetric_poly ( nn, rr, x );
      fprintf ( 1, '  %6d', value );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [out] = spm_eeg_displayECD(Pos,Orient,Var,Names,options)
% Plot dipole positions onto the SPM canonical mesh
% FORMAT [out] = spm_eeg_displayECD(Pos,Orient,Var,Names,options)
%
% IN (admissible choices):
%   - Pos: a 3xndip matrix containing the positions of the dipoles in
%   the canonical frame of reference
%   - Orient: the same with dipole orientations
%   - Var: the same with position variance
%   - Names: the same with dipole names
%   - options: an optional structure containing
%       .hfig: the handle of the display figure
%       .tag: the tag to be associated with the created UI objects
%       .add: binary variable ({0}, 1: just add dipole in the figure .hfig)
%
% OUT:
%   - out: a structure containing the handles of the object in the figure
%   (including the mesh, the dipoles, the transparency slider, etc...)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jean Daunizeau
% $Id: spm_eeg_displayECD.m 5737 2013-11-10 20:23:49Z karl $

% checks and defaults
%--------------------------------------------------------------------------
hfig       = [];
ParentAxes = [];
query      = [];
handles    = [];
tag        = '';
try, options; catch, options = [];      end
try, hfig        = options.hfig;        end
try, tag         = options.tag;         end
try, ParentAxes  = options.ParentAxes;  end
try, query       = options.query;       end
try, handles     = options.handles;     end
try
    figure(hfig);
catch
    hfig  = spm_figure('GetWin','Graphics');
    spm_figure('Clear',hfig);
    ParentAxes = axes('parent',hfig);
end
try
    markersize = options.markersize;
catch
    markersize = 20;
end
try
    meshsurf = options.meshsurf;
catch
    meshsurf = fullfile(spm('Dir'),'canonical','cortex_5124.surf.gii');
end

if isscalar(Var), Var = Pos*0 + Var^2;   end
try, Pos{1};    catch, Pos = {Pos};      end
try, Orient{1}; catch, Orient = {Orient};end
try, Var{1};    catch, Var = {Var};      end

ndip = size(Pos{1},2);
if ~exist('Names','var') || isempty(Names)
    for i=1:ndip
        Names{i} = num2str(i);
    end
end

col = ['b','g','r','c','m','y','k','w'];
tmp = ceil(ndip./numel(col));
col = repmat(col,1,tmp);
pa  = get(ParentAxes,'position');

if ndip > 0

    if isempty(query)
        opt.hfig = hfig;
        opt.ParentAxes = ParentAxes;
        opt.visible = 'off';
        pos2 = [pa(1),pa(2)+0.25*pa(4),0.03,0.5*pa(4)];
        out  = spm_eeg_render(meshsurf,opt);
        handles.mesh = out.handles.p;
        handles.BUTTONS.transp = out.handles.transp;
        handles.hfig = out.handles.fi;
        handles.ParentAxes  = out.handles.ParentAxes;
        set(handles.mesh,...
            'facealpha',0.1,...
            'visible','on')
        set(handles.BUTTONS.transp,...
            'value',0.1,...
            'position',pos2,...
            'visible','on')
    end

    set(ParentAxes,'nextplot','add')
    for j=1:length(Pos)
        for i =1:ndip
            try
                set(handles.hp(j,i),...
                    'xdata',Pos{j}(1,i),...
                    'ydata',Pos{j}(2,i),...
                    'zdata',Pos{j}(3,i));
            catch
                handles.hp(j,i) = plot3(handles.ParentAxes,...
                    Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                    [col(i),'.'],...
                    'markerSize',markersize,...
                    'visible','off');
            end
            try
                no = sqrt(sum(Orient{j}(:,i).^2));
                if no > 0
                    Oi = 1e1.*Orient{j}(:,i)./no;
                else
                    Oi = 1e-5*ones(3,1);
                end
                try
                    set(handles.hq(j,i),...
                        'xdata',Pos{j}(1,i),...
                        'ydata',Pos{j}(2,i),...
                        'zdata',Pos{j}(3,i),...
                        'udata',Oi(1),...
                        'vdata',Oi(2),...
                        'wdata',Oi(3))
                catch
                    handles.hq(j,i) = quiver3(handles.ParentAxes,...
                        Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                        Oi(1),Oi(2),Oi(3),col(i),...
                        'lineWidth',2,'visible','off');
                end
                if isequal(query,'add')
                    set(handles.hq(j,i),...
                        'LineStyle','--',...
                        'lineWidth',1)
                end
            end
            [x,y,z]= ellipsoid(Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                1.*sqrt(Var{j}(1,i)),1.*sqrt(Var{j}(2,i)),1.*sqrt(Var{j}(1,i)),20);
            try
                set(handles.hs(j,i),...
                    'xdata',x,...
                    'ydata',y,...
                    'zdata',z);
            catch
                handles.hs(j,i) = surf(handles.ParentAxes,...
                    x,y,z,...
                    'edgecolor','none',...
                    'facecolor',col(i),...
                    'facealpha',0.2,...
                    'visible','off');
            end
            try
                set(handles.ht(j,i),...
                    'position',Pos{j}(:,i));
            catch
                handles.ht(j,i) = text(...
                    Pos{j}(1,i),Pos{j}(2,i),Pos{j}(3,i),...
                    Names{i},...
                    'Parent',handles.ParentAxes,...
                    'visible','off');
            end
        end
    end

    if length(Pos) > 1

        try, set(handles.hp(end,:),'visible','on'); end
        try, set(handles.hq(end,:),'visible','on'); end
        try, set(handles.hs(end,:),'visible','on'); end
        try, set(handles.ht(end,:),'visible','on'); end

        handles.uic(1) = uicontrol(handles.fi,...
            'units','normalized',...
            'position',[0.45,0.5,0.2,0.03],...
            'style','radio','string','Show priors',...
            'callback',@doChange1,...
            'BackgroundColor',[1 1 1],...
            'tooltipstring','Display prior locations',...
            'userdata',handles,'value',0,...
            'BusyAction','cancel',...
            'Interruptible','off',...
            'tag','plotEEG');

        handles.uic(2) = uicontrol(handles.fi,...
            'units','normalized',...
            'position',[0.45,0.53,0.2,0.03],...
            'style','radio','string','Show posteriors',...
            'callback',@doChange2,...
            'BackgroundColor',[1 1 1],...
            'tooltipstring','Display posterior locations',...
            'userdata',handles,'value',1,...
            'BusyAction','cancel',...
            'Interruptible','off',...
            'tag','plotEEG');

    else

        try, set(handles.hp(1,:),'visible','on'); end
        try, set(handles.hq(1,:),'visible','on'); end
        try, set(handles.hs(1,:),'visible','on'); end
        try, set(handles.ht(1,:),'visible','on'); end

    end

end

try
    clear out
    out.handles = handles;
catch
    out = [];
end

%==========================================================================
function doChange1(i1,i2)
val = get(i1,'value');
handles = get(i1,'userdata');
if ~val
    try, set(handles.hp(1,:),'visible','off'); end
    try, set(handles.hq(1,:),'visible','off'); end
    try, set(handles.hs(1,:),'visible','off'); end
    try, set(handles.ht(1,:),'visible','off'); end
else
    try, set(handles.hp(1,:),'visible','on');  end
    try, set(handles.hq(1,:),'visible','on');  end
    try, set(handles.hs(1,:),'visible','on');  end
    try, set(handles.ht(1,:),'visible','on');  end
end

%==========================================================================
function doChange2(i1,i2)
val = get(i1,'value');
handles = get(i1,'userdata');
if ~val
    try, set(handles.hp(2,:),'visible','off'); end
    try, set(handles.hq(2,:),'visible','off'); end
    try, set(handles.hs(2,:),'visible','off'); end
    try, set(handles.ht(2,:),'visible','off'); end
else
    try, set(handles.hp(2,:),'visible','on');  end
    try, set(handles.hq(2,:),'visible','on');  end
    try, set(handles.hs(2,:),'visible','on');  end
    try, set(handles.ht(2,:),'visible','on');  end
end

author: spm


text:
function rtk=udtrop_rtkins(rtk,tt,bl) %#ok

global glc
INIT_ZWD=0.15; VAR_GRA=0.001^2;

for i=1:2
    if i==1,j=rtk.itr+1;else,j=rtk.itb+1;end

    if rtk.x(j)==0
        rtk=initx(rtk,INIT_ZWD,rtk.opt.std(3)^2,j);
        if rtk.opt.tropopt==glc.TROPOPT_ESTG
            for k=1:2
                rtk=initx(rtk,1e-6,VAR_GRA,j+k);
            end
        end
    else
        rtk.P(j,j)=rtk.P(j,j)+rtk.opt.prn(3)^2*abs(tt);
        if rtk.opt.tropopt==glc.TROPOPT_ESTG
            for k=1:2
                rtk.P(j+k,j+k)=rtk.P(j+k,j+k)+(rtk.opt.prn(3)*0.3)^2*abs(tt);
            end
        end
    end

end

return

repo: GINav