text:
function [u,p,edgeC,A,eqn,info] = StokesisoP2P1(node,elem,pde,bdFlag,option)
%% STOKESISOP2P1 Stokes equation: isoP2-P1 modified Taylor-Hood elements.
%
%  [u,p] = STOKESisoP2P1(node,elem,pde,bdFlag) use constinous P1 element
%  on grid h and continous P1 element on grid H = 2*h to approximate
%  velocity u and pressure p, repectively.
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%   with
%       Dirichlet boundary condition        u = g_D  on \Gamma_D,
%       Neumann boundary condition du/dn - np = g_N  on \Gamma_N.
%
% Created by Ming Wang at Aug., 2012.
%
% See also StokesisoP2P0, StokesP2P1
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end

%% Refine grid for P1
nodeC = node; elemC = elem; bdFlagC = bdFlag;
[tempvar,edgeC] = dofP2(elemC);
[node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
NC = size(nodeC,1); NTC = size(elemC,1);
N = size(node,1);  NT = size(elem,1); Nu = N; Np = NC;

tic;
%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);
areaC = sum(reshape(area,NTC,4),2);

%% Assemble stiffness matrix for Laplace operator
A = sparse(Nu,Nu);
for i = 1:3
    for j = i:3
        Aij = (Dlambda(:,1,i).*Dlambda(:,1,j) + ...
            Dlambda(:,2,i).*Dlambda(:,2,j)).*area;
        if isfield(pde,'mu') && (pde.mu~=1)
            Aij = pde.mu*Aij;
        end
        if (j==i)
            A = A + sparse(elem(:,i),elem(:,j),Aij,N,N);
        else
            A = A + sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],...
                [Aij; Aij],N,N);
        end
    end
end
clear Aij
A = blkdiag(A,A);

%% Assemble the matrix for divergence operator
% idea: Basis on coarse grids can be expanded from fine grids, i.e.,
%       lambda_{i,c} = lambda_{i,f} + 1/2*sum_{j \ in V(i)} lambda_{j,f}
%        where V(i) is the index of points on edges surrounding point i.
Dx = sparse(N,N);
Dy = sparse(N,N);
for j = 1:3 % loop for u index
    Dx = Dx + sparse(elem(:),repmat(elem(:,j),3,1),...
                     repmat(1/3*Dlambda(:,1,j).*area,3,1),N,N);
    Dy = Dy + sparse(elem(:),repmat(elem(:,j),3,1),...
                     repmat(1/3*Dlambda(:,2,j).*area,3,1),N,N);
end
Bf = [-Dx -Dy];
Bv = Bf(1:NC,:);
Be = Bf(NC+1:end,:);
D = 1/2*icdmat(double(edgeC),[1,1]);
B = Bv + D'*Be;

%% Assemble right hand side by 4-points quadrature rule
f1 = zeros(Nu,1);
f2 = zeros(Nu,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,weight] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
    ft1 = zeros(NT,3);
    ft2 = zeros(NT,3);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        % function values at quadrature points
        fp = pde.f(pxy);
        % evaluate fp outside.
        for j = 1:3
            ft1(:,j) = ft1(:,j) + fp(:,1).*lambda(p,j)*weight(p);
            ft2(:,j) = ft2(:,j) + fp(:,2).*lambda(p,j)*weight(p);
        end
    end
    ft1 = ft1.*repmat(area,1,3);
    ft2 = ft2.*repmat(area,1,3);
    f1 = accumarray(elem(:),ft1(:),[Nu 1]);
    f2 = accumarray(elem(:),ft2(:),[Nu 1]);
end

[AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesisoP2P1;

%% Record assembeling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(ufreeDof), return; end
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if length(f)+length(g) <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else          % Multigrid-type  solver for large size systems
        option.solver = 'asmg';
    end
end
solver = option.solver;

%% Solver
switch solver
    case 'direct'
        tic;
        bigA = [AD, BD'; ...
                BD, sparse(Np,Np)];
        bigF = [f; g];
        bigu = [u; p];
        bigFreeDof = [ufreeDof; 2*Nu+pDof];
        bigu(bigFreeDof) = bigA(bigFreeDof,bigFreeDof)\bigF(bigFreeDof);
        u = bigu(1:2*Nu);
        p = bigu(2*Nu+1:end);
        residual = norm(bigF - bigA*bigu);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
        option.solver  = 'WCYCLE';
        [u(ufreeDof),p,info] = mgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                        u(ufreeDof),p,elemC,ufreeDof,option);
    case 'asmg'
        [u(ufreeDof),p,info] = asmgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                          u,p,nodeC,elemC,bdFlagC,ufreeDof,option);
end

%% Post-process
if length(pDof)~=Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(mean(p(elemC),2).*areaC)/sum(areaC);
    p = p - c;
end

%% Output information
eqn = struct('A',AD,'B',BD,'f',f,'g',g,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdStokesisoP2P1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesisoP2P1
        %% Initial set up
        g = zeros(Np,1);
        u = zeros(2*Nu,1);
        p = zeros(Np,1);
        ufreeDof = (1:Nu)';
        pDof = (1:Np)';
        if ~exist('bdFlag','var'), bdFlag = []; end
        if ~isfield(pde,'g_D'), pde.g_D = []; end
        if ~isfield(pde,'g_N'), pde.g_N = []; end
        if ~isfield(pde,'g_R'), pde.g_R = []; end

        %% Part 1: Find Dirichlet dof and modify the matrix
        % Find Dirichlet boundary dof: fixedDof and pDof
        isFixedDof = false(Nu,1);
        if ~isempty(bdFlag) % case: bdFlag is not empty
            allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
            Dirichlet = allEdge((bdFlag(:) == 1),:);
            isFixedDof(Dirichlet(:)) = true;
            fixedDof = find(isFixedDof);
            ufreeDof = find(~isFixedDof);
        end
        if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N)
            fixedDof = findboundary(elem);
            isFixedDof(fixedDof) = true;
            ufreeDof = find(~isFixedDof);
        end
        if isempty(fixedDof) % pure Neumann boundary condition
            % pde.g_N could be empty which is homogenous Neumann boundary condition
            fixedDof = 1;
            ufreeDof = 2:Nu;    % eliminate the kernel by enforcing u(1) = 0;
        end

        % Modify the matrix
        % Build Dirichlet boundary condition into the matrix AD by enforcing
        % AD(fixedDof,fixedDof)=I, AD(fixedDof,ufreeDof)=0, AD(ufreeDof,fixedDof)=0.
        % BD(:,fixedDof) = 0 and thus BD'(fixedDof,:) = 0.
        bdidx = zeros(2*Nu,1);
        bdidx(fixedDof) = 1;
        bdidx(Nu+fixedDof) = 1;
        Tbd = spdiags(bdidx,0,2*Nu,2*Nu);
        T = spdiags(1-bdidx,0,2*Nu,2*Nu);
        AD = T*A*T + Tbd;
        BD = B*T;

        %% Part 2: Find boundary edges and modify the right hand side f and g
        % Find boundary edges: Neumann and Robin
        Neumann = []; Robin = []; %#ok<*NASGU>
        if ~isempty(bdFlag)
            allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
            Neumann = allEdge((bdFlag(:)==2)|(bdFlag(:) == 3),:);
            Robin = allEdge((bdFlag(:) == 3),:);
        end
        if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
            % no bdFlag, only pde.g_N or pde.g_R is given in the input
            [tempvar,Neumann] = findboundary(elem);
            if ~isempty(pde.g_R)
                Robin = Neumann;
            end
        end

        % Neumann boundary condition
        if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
            [lambda,w] = quadpts1(3);
            nQuad = size(lambda,1);
            ve = node(Neumann(:,1),:) - node(Neumann(:,2),:); % length of edge
            edgeLength = sqrt(sum(ve.^2,2));
            % update RHS
            gex = zeros(size(Neumann,1),2);   % x-component
            gey = zeros(size(Neumann,1),2);   % y-component
            for pp = 1:nQuad
                pxy = lambda(pp,1)*node(Neumann(:,1),:)+lambda(pp,2)*node(Neumann(:,2),:);
                gp = pde.g_N(pxy);
                gex(:,1) = gex(:,1) + w(pp)*edgeLength.*gp(:,1)*lambda(pp,1);
                gex(:,2) = gex(:,2) + w(pp)*edgeLength.*gp(:,1)*lambda(pp,2);
                gey(:,1) = gey(:,1) + w(pp)*edgeLength.*gp(:,2)*lambda(pp,1);
                gey(:,2) = gey(:,2) + w(pp)*edgeLength.*gp(:,2)*lambda(pp,2);
            end
            f1(1:N) = f1(1:N) + accumarray(Neumann(:), gex(:),[N,1]);
            f2(1:N) = f2(1:N) + accumarray(Neumann(:), gey(:),[N,1]);
        end
        f = [f1; f2];
        % The case non-empty Neumann but g_N=[] corresponds to the zero flux
        % boundary condition on Neumann edges and no modification is needed.

        % Dirichlet boundary conditions
        if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
            u1 = zeros(Nu,1);
            u2 = zeros(Nu,1);
            uD = pde.g_D(node(fixedDof,:));
            u1(fixedDof) = uD(:,1);
            u2(fixedDof) = uD(:,2);
            u = [u1;u2];
            f = f - A*u;  % bring affect of nonhomgenous Dirichlet bd condition
            g = g - B*u;  % to the right hand side
            g = g - mean(g);
            f(fixedDof)    = u1(fixedDof);
            f(fixedDof+Nu) = u2(fixedDof);
        end
        % The case non-empty Dirichlet but g_D=[] corresponds to the zero Dirichlet
        % boundary condition and no modification is needed.

        % modfiy pressure dof for pure Dirichlet
        if isempty(Neumann)
            pDof = (1:Np-1)';
        end

        ufreeDof = [ufreeDof; Nu+ufreeDof];
    end % end of function getbdStokesisoP2P1
end

repo: Aero-matlab


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

function total_cost = computeConstraintCustomizedCost(X, theta_x, theta_y, theta_z, T, target_size, box_width)
    cost_x = 0;
    cost_y = 0;
    cost_z = 0;
    R = rotx(theta_x) * roty(theta_y) * rotz(theta_z);
    x_prime = R(1, :) * X(1:3, :) + T(1);
    y_prime = R(2, :) * X(1:3, :) + T(2);
    z_prime = R(3, :) * X(1:3, :) + T(3);
    for i = 1:size(X, 2)
        cost_z = cost_z + checkCost(z_prime(i), -target_size/2, target_size/2);
        cost_y = cost_y + checkCost(y_prime(i), -target_size/2, target_size/2);
        cost_x = cost_x + checkCost(x_prime(i), -box_width, box_width);
    end
%     total_cost = sqrt(cost_x + cost_y + cost_z);
    total_cost = cost_x + cost_y + cost_z;
end

author: UMich-BipedLab


text:
function [P_intersect,distances] = lineIntersect3D(PA,PB)
% Find intersection point of lines in 3D space, in the least squares sense.
% PA :          Nx3-matrix containing starting point of N lines
% PB :          Nx3-matrix containing end point of N lines
% P_Intersect : Best intersection point of the N lines, in least squares sense.
% distances   : Distances from intersection point to the input lines
% Anders Eikenes, 2012

Si = PB - PA; %N lines described as vectors
ni = Si ./ (sqrt(sum(Si.^2,2))*ones(1,3)); %Normalize vectors
nx = ni(:,1); ny = ni(:,2); nz = ni(:,3);
SXX = sum(nx.^2-1);
SYY = sum(ny.^2-1);
SZZ = sum(nz.^2-1);
SXY = sum(nx.*ny);
SXZ = sum(nx.*nz);
SYZ = sum(ny.*nz);
S = [SXX SXY SXZ;SXY SYY SYZ;SXZ SYZ SZZ];
CX  = sum(PA(:,1).*(nx.^2-1) + PA(:,2).*(nx.*ny)  + PA(:,3).*(nx.*nz));
CY  = sum(PA(:,1).*(nx.*ny)  + PA(:,2).*(ny.^2-1) + PA(:,3).*(ny.*nz));
CZ  = sum(PA(:,1).*(nx.*nz)  + PA(:,2).*(ny.*nz)  + PA(:,3).*(nz.^2-1));
C   = [CX;CY;CZ];
P_intersect = (S\C)';

if nargout>1
    N = size(PA,1);
    distances=zeros(N,1);
    for i=1:N %This is faster:
        ui=(P_intersect-PA(i,:))*Si(i,:)'/(Si(i,:)*Si(i,:)');
        distances(i)=norm(P_intersect-PA(i,:)-ui*Si(i,:));
    end
    %for i=1:N %http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html:
    %    distances(i) = norm(cross(P_intersect-PA(i,:),P_intersect-PB(i,:))) / norm(Si(i,:));
    %end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function result=BF0NLS(speechSignal,samplingFreq,plot_flag,prew_flag,Par_set)
% -----------------------------------------------------------------------
% format of the input and output arguments
%
% input:
% input arguments --> as the names suggested
% plot_flag: 0 when you do not want a plot  (0 in default)
% prew_flag: 0 when you do not want to prewhiten the signal (0 in default)
% Par_set:   the parameter set you can set, including the segment time,
% segment shift and F0 boundaries.
%
% Par_set.segmentTime is the segmentTime in seconds  (25 ms in default)
% Par_set.segmentShift is the segmentShift in seconds (10 ms in default)
% Par_set.f0Bounds is the F0 boundaries in Hz ([70 400] Hz in default)
%
% output:
% result.tt         -->   time vector
% result.ff         -->   Fundamental frequency estimates
%                         (when all the frames are considered as voiced)
% result.oo         -->   order estimates
%                         (when all the frames are considered as voiced)
% result.vv         -->   voicing probability
% result.best_ff    -->   Best fundamental frequency estimates (setting the F0=nan
%                         when voicing probability is less than .5)
% result.best_order -->   Best harmonic order estimates (setting the order=nan
%                         when voicing probability is less than .5)
%
% Example for for customize Par_set
%
% Par_set.segmentTime = 0.025; 25 ms for each segment (default value)
% Par_set.segmentShift = 0.01; 10 ms for segment shift (default value)
% Par_set.f0Bounds =[70, 400]; pitch is bounded between 70 to 400 Hz (default value)
%
%
% Written by Liming Shi, Aalborg 9000, Denmark
% Email: ls@create.aau.dk
% -----------------------------------------------------------------------

if nargin<3
%   do not use prewhitening in default, and do not plot
    plot_flag=0;
    prew_flag=0;
    Par_set=[];
end
if nargin<4
%   do not use prewhitening in default
    prew_flag=0;
    Par_set=[];
end
if nargin<5
%   use default parameter set
    Par_set=[];
end

%% resample to 16 KHz for tuned parameters. However, the sampling frequency can be changed.
fs_fine_tuned=16000;
speechSignal=resample(speechSignal,fs_fine_tuned,samplingFreq);
samplingFreq=fs_fine_tuned;
if isempty(Par_set)
    segmentTime = 0.025; %     seconds
    segmentShift = 0.010; % seconds
    f0Bounds = [70, 400]; % cycles/sample
    std_pitch=2;
else
    if isfield(Par_set,'segmentTime')
        segmentTime = Par_set.segmentTime; %     seconds
    else
        segmentTime = 0.025; %     seconds
    end
    if isfield(Par_set,'segmentShift')
        segmentShift = Par_set.segmentShift;
        % std_pitch --> standard deviation for the pitch for 10 ms frame shift.
        % This value shouold be modified if frame shift is shorter or longer, with a lower bound of .5 Hz.
        std_pitch=max(2/10*segmentShift*1000,0.5);
    else
        segmentShift = 0.010; % seconds
        % std_pitch --> standard deviation for the pitch for 10 ms frame shift.
        std_pitch=2;
    end
    if isfield(Par_set,'f0Bounds')
        f0Bounds = Par_set.f0Bounds; % cycles/sample
    else
        f0Bounds = [70, 400]; % cycles/sample
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normalization step
sumE = sqrt(speechSignal'*speechSignal/length(speechSignal));
scale = sqrt(3.1623e-5)/sumE; % scale to -45-dB loudness level
speechSignal=speechSignal*scale;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nData = length(speechSignal);
% set up
f0Bounds=f0Bounds/samplingFreq;
segmentLength = round(segmentTime*samplingFreq/2)*2; % samples
nShift = round(segmentShift*samplingFreq); % samples
nSegments = floor((nData+segmentLength/2-segmentLength)/nShift)+1;

if prew_flag==0
    maxNoHarmonics = 10;
else
    if prew_flag==1
        maxNoHarmonics = 30;
    end
end
f0Estimator = BayesianfastF0NLS(segmentLength, maxNoHarmonics, f0Bounds, std_pitch/samplingFreq,.7);
speechSignal_padded=[zeros(segmentLength/2,1);speechSignal];
% do the analysis
idx = 1:segmentLength;
f0Estimates = nan(nSegments,1); % cycles/sample
order=nan(nSegments,1);
voicing_prob=nan(nSegments,1);
if prew_flag==0
    for ii = 1:nSegments
        speechSegment = speechSignal_padded(idx);
        [f0Estimates(ii),order(ii),voicing_prob(ii)]=f0Estimator.estimate(speechSegment,0);
        idx = idx + nShift;
    end
else
    if prew_flag==1
        for ii = 1:nSegments
            speechSegment = speechSignal_padded(idx);
            [f0Estimates(ii),order(ii),voicing_prob(ii)]=f0Estimator.estimate(speechSegment,1,segmentShift);
            idx = idx + nShift;
        end
    end
end

f0Estimates_remove_unvoiced=f0Estimates;
unvoiced_indicator=voicing_prob<.5;;
f0Estimates_remove_unvoiced(unvoiced_indicator)=nan;
order_remove_unvoiced=order;
order_remove_unvoiced(unvoiced_indicator)=nan;
timeVector = [(0:nSegments-1)*segmentShift]';

result.tt=timeVector;
result.ff=f0Estimates*samplingFreq;
result.oo=order;
result.vv=voicing_prob;
result.best_ff=f0Estimates_remove_unvoiced*samplingFreq;
result.best_order=order_remove_unvoiced;
if plot_flag==1
    figure;
    subplot(4,1,4)
    plot([0:length(speechSignal_padded)-1]/samplingFreq,speechSignal_padded/max(abs(speechSignal_padded)));
    xlim([0,(length(speechSignal_padded)-1)/samplingFreq])
    xlabel('Time [s]');
    ylabel('Amplitude');
    subplot(4,1,3)
%   plot the spectrogram
    window = gausswin(segmentLength);
    nOverlap = segmentLength-nShift;
    nDft = 2048;
    [stft, stftFreqVector] = ...
        spectrogram(speechSignal_padded, window, nOverlap, nDft, samplingFreq);
    powerSpectrum = abs(stft).^2;
    imagesc(timeVector, stftFreqVector, ...
    10*log10(dynamicRangeLimiting(powerSpectrum, 60)));
    axis xy;
    ylim([0,f0Bounds(2)*samplingFreq+100])
    hold on;
    plot(timeVector, result.best_ff, 'r-', 'linewidth',2);
    ylabel('Frequency [Hz]');
    subplot(4,1,2)
    plot(timeVector,result.best_order,'r.', 'linewidth',2)
    xlim([timeVector(1),timeVector(end)])
    ylabel('Order estimate')
    subplot(4,1,1)
    plot(timeVector,result.vv,'r-', 'linewidth',2)
    xlim([timeVector(1),timeVector(end)])
    ylabel('Voicing probability')
end

end

author: LimingShi


text:
function [ responses ] = CCNF_ncc_response( patches, patch_experts, normalisation_options, window_size, patch_length)
%PATCHRESPONSESVM Computing a patch response from a CCNF patch expert
%   Using convolution, for testing purposes and not for actual speed

    SigmaInv = patch_experts.SigmaInv;

    patchSize = normalisation_options.patchSize;

    if(~iscell(patches))
        patches = {patches};
    end

    num_modalities = numel(patches);

    responses = zeros(size(patches{1},1), patch_length);

    % prepare the patches by normalising them to zscore (if used)
    if(normalisation_options.zscore)
        for i=1:num_modalities
            patches{i} = zscore(patches{i});
        end
    end

    for i = 1:size(patches{1},1)

        norm_cross_corr = normalisation_options.useNormalisedCrossCorr == 1;

        b = zeros(patch_length,1);

        hl_per_modality = size(patch_experts.thetas,1);

        for p=1:num_modalities
            smallRegionVec = patches{p}(i,:);
            smallRegion = reshape(smallRegionVec, window_size(1), window_size(2));

            for hls = 1:hl_per_modality

                % because the normalised cross correlation calculates the
                % responses from a normalised template and a normalised image,
                % normalise the thetas here and then apply the normalisation to
                % the response

                w = patch_experts.thetas(hls, 2:end, p);
                norm_w = norm(w);
                w = w/norm(w);
                w = reshape(w, patchSize);

                response = -norm_w * Cross_corr_resp(smallRegion, w, norm_cross_corr, patchSize) - patch_experts.thetas(hls,1,p);

                % here we include the bias term as well, as it wasn't added
                % during the response calculation
                h1 = 1./(1 + exp(response(:)));
                b = b + (2 * patch_experts.alphas((p-1)*hl_per_modality + hls) * h1);

            end
        end
        response = SigmaInv \ b;

        responses(i,:) = response(:);

    end
    responses = responses';
    responses = responses(:);
end

function response = Cross_corr_resp(region, patchExpert, normalise_x_corr,patchSize)

    if(normalise_x_corr)
        [response] = normxcorr2(patchExpert, region);
        response = response(patchSize(1):end-patchSize(1)+1,patchSize(2):end-patchSize(2)+1);
    else
        % this assumes that the patch is already normed, so just use
        % cross-correlation
        template = rot90(patchExpert,2);
        response = conv2(region, template, 'valid');

    end
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function deg=degree(p,y,e)
%DEGREE Polynomial degree
%
% DEG = DEGREE(p,x,e)
%
% p : SDPVAR object.
% x : Degree w.r.t linear SDPVAR objects, can be [].

% e : If e=1, returns degree of each element in p
%
% Examples
% x1 = sdpvar(1,1);x2 = sdpvar(1,1);
% p = [x1;x1*x2+x2^2];
%
% degree(p) returns 2
%
% degree(p,x1) returns 1
%
% degree(p,[x1 x2]) returns [1 2]
%
% degree(p,[x1 x2],1) returns [1 0;1 2]
%
% degree(p,[],1) returns [1;3]

if isa(p,'double')
    if nargin==1
        deg = 0;
    else
        deg = zeros(1,length(y));
    end
    return
end

if nargin<2
    y = recover(depends(p));
end

if nargin<3 | (nargin==3 & e==0)
    exponent_p = exponents(p,y);
    switch nargin
        case 1
            deg = full(max(sum(exponent_p,2)));
        case {2,3}
            deg = full(max(exponent_p,[],1));
        otherwise
            error('Too many arguments. Wadda ya mean?')
    end
else
    p = p(:);
    if isempty(y)
        yy = recover(depends(p));
    else
        yy = y;
    end

    for i = 1:length(p)
        z.type = '()';
        z.subs{1} = i;
        exponent_p = exponents(subsref(p,z),yy);
        switch nargin
            case 1
                deg(i,:) = full(max(sum(exponent_p,2)));
            case {2,3}
                deg(i,:) = full(max(exponent_p,[],1));
            otherwise
               error('Too many arguments. Wadda ya mean?')
        end
    end
    if isempty(y)
        deg = sum(deg,2);
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function definput=arg_firwin(definput)

truncgaussflags = {'truncgauss'};

definput.flags.wintype=[{...
    'hanning','hann','sine','cosine', 'sqrthann','hamming', 'square',...
    'rect', 'tria','bartlett', 'triangular','sqrttria','blackman',...
    'blackman2', 'nuttall','nuttall10','nuttall01','nuttall20',...
    'nuttall11','nuttall03', 'nuttall12','nuttall21', 'nuttall30',...
    'ogg', 'itersine', 'nuttall02'}, truncgaussflags ];

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%%******************************************************************
%% checkdepconst: compute AAt to determine if the
%%             constraint matrices Ak are linearly independent.
%%
%% [At,b,y,idxB,neardepconstr,feasible,AAt] =
%%         checkdepconstr(blk,At,b,y,rmdepconstr);
%%
%% rmdepconstr = 1, if want to remove dependent columns in At
%%             = 0, otherwise.
%%
%% idxB = indices of linearly independent columns of original At.
%% neardepconstr = 1 if there is nearly dependent columns in At
%%            = 0, otherwise.
%% Note: the definition of "nearly dependent" is dependent on the
%%       threshold used to determine the small diagonal elements in
%%       the LDLt factorization of A*At.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function  [At,b,y,idxB,neardepconstr,feasible,AAt] = ...
              checkdepconstr(blk,At,b,y,rmdepconstr);

   global existlowrank printlevel
   global nnzmatold
%%
%% compute AAt
%%
   m = length(b);
   AAt = sparse(m,m);  numdencol = 0; UU = [];
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s');
         m1 = size(At{p,1},2); m2 = m - m1;
         if (m2 > 0)
            if (m2 > 0.3*m); AAt = full(AAt); end
            dd = At{p,3};
            lenn = size(dd,1);
            idxD = [0; find(diff(dd(:,1))); lenn];
            ss = [0, cumsum(pblk{3})];
            if (existlowrank)
               AAt(1:m1,1:m1) = AAt(1:m1,1:m1) + At{p,1}'*At{p,1};
               for k = 1:m2
                  idx = [ss(k)+1 : ss(k+1)];
                  idx2 = [idxD(k)+1: idxD(k+1)];
                  ii = dd(idx2,2)-ss(k); %% undo cumulative indexing
                  jj = dd(idx2,3)-ss(k);
                  len = pblk{3}(k);
                  Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                  tmp = svec(pblk,At{p,2}(:,idx)*Dk*At{p,2}(:,idx)');
                  tmp2 = AAt(1:m1,m1+k) + (tmp'*At{p,1})';
                  AAt(1:m1,m1+k) = tmp2;
                  AAt(m1+k,1:m1) = tmp2';
               end
            end
            DD = spconvert([dd(:,2:4); sum(pblk{3}),sum(pblk{3}),0]);
            VtVD = (At{p,2}'*At{p,2})*DD;
            VtVD2 = VtVD'.*VtVD;
            for k = 1:m2
               idx0 = [ss(k)+1 : ss(k+1)];
               %%tmp = VtVD(idx0,:)' .* VtVD(:,idx0);
               tmp = VtVD2(:,idx0);
               tmp = tmp*ones(length(idx0),1);
               tmp3 = AAt(m1+[1:m2],m1+k) + mexqops(pblk{3},tmp,ones(length(tmp),1),1);
               AAt(m1+[1:m2],m1+k) = tmp3;
            end
         else
            AAt = AAt + abs(At{p,1})'*abs(At{p,1});
         end
      else
         decolidx = checkdense(At{p,1}');
         if ~isempty(decolidx);
            n2 = size(At{p,1},1);
            dd = ones(n2,1);
            len= length(decolidx);
            dd(decolidx) = zeros(len,1);
            AAt = AAt + (abs(At{p,1})' *spdiags(dd,0,n2,n2)) *abs(At{p,1});
            tmp = At{p,1}(decolidx,:)';
            UU = [UU, tmp];
            numdencol = numdencol + len;
         else
            AAt = AAt + abs(At{p,1})'*abs(At{p,1});
         end
      end
   end
   if (numdencol > 0) & (printlevel)
      fprintf('\n number of dense column in A = %d',numdencol);
   end
   numdencol = size(UU,2);
%%
%%
%%
   feasible = 1; neardepconstr = 0;
   if ~issparse(AAt); AAt = sparse(AAt); end
   nnzmatold = mexnnz(AAt);
   rho = 1e-15;
   diagAAt = diag(AAt);
   mexschurfun(AAt,rho*max(diagAAt,1));
   [L.R,indef,L.perm] = chol(AAt,'vector');
   L.d = full(diag(L.R)).^2;
   if (indef)
      msg = 'AAt is not pos. def.';
      idxB = [1:m];
      neardepconstr = 1;
      if (printlevel); fprintf('\n checkdepconstr: %s',msg); end
      return;
   end
%%
%% find independent rows of A
%%
   dd = zeros(m,1);
   idxB = [1:m]';
   dd(L.perm) = abs(L.d);
   idxN = find(dd < 1e-13*mean(L.d));
   ddB = dd(setdiff([1:m],idxN));
   ddN = dd(idxN);
   if ~isempty(ddN) & ~isempty(ddB) & (min(ddB)/max(ddN) < 10)
      %% no clear separation of elements in dd
      %% do not label constraints as dependent
      idxN = [];
   end
   if ~isempty(idxN)
      neardepconstr = 1;
      if (printlevel)
         fprintf('\n number of nearly dependent constraints = %1.0d',length(idxN));
      end
      if (numdencol==0)
         if (rmdepconstr)
            idxB = setdiff([1:m]',idxN);
            if (printlevel)
               fprintf('\n checkdepconstr: removing dependent constraints...');
            end
            [W,resnorm] = findcoeffsub(blk,At,idxB,idxN);
   	    tol = 1e-8;
            if (resnorm > sqrt(tol))
               idxB = [1:m]';
               neardepconstr = 0;
               if (printlevel)
                  fprintf('\n checkdepconstr: basis rows cannot be reliably identified,');
                  fprintf(' abort removing nearly dependent constraints');
               end
               return;
            end
            tmp = W'*b(idxB) - b(idxN);
            nnorm = norm(tmp)/max(1,norm(b));
            if (nnorm > tol)
               feasible = 0;
               if (printlevel)
                  fprintf('\n checkdepconstr: inconsistent constraints exist,');
                  fprintf(' problem is infeasible.');
               end
            else
               feasible = 1;
               for p = 1:size(blk,1)
                  At{p,1} = At{p,1}(:,idxB);
               end
               b = b(idxB);
               y = y(idxB);
               AAt = AAt(idxB,idxB);
            end
	 else
            if (printlevel)
               fprintf('\n To remove these constraints,');
               fprintf(' re-run sqlp.m with OPTIONS.rmdepconstr = 1.');
            end
         end
      else
         if (printlevel)
            fprintf('\n warning: the sparse part of AAt may be nearly singular.');
         end
      end
   end
%%*****************************************************************************
%% findcoeffsub:
%%
%% [W,resnorm] = findcoeffsub(blk,At,idXB,idXN);
%%
%% idXB = indices of independent columns of At.
%% idxN = indices of   dependent columns of At.
%%
%% AB = At(:,idxB); AN = At(:,idxN) = AB*W
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

   function [W,resnorm] = findcoeffsub(blk,At,idxB,idxN);

   AB = []; AN = [];
   for p = 1:size(blk,1)
      AB = [AB; At{p,1}(:,idxB)];
      AN = [AN; At{p,1}(:,idxN)];
   end
   [m,n] = size(AB);
%%
%%-----------------------------------------
%% find W so that AN = AB*W
%%-----------------------------------------
%%
   [L,U,P,Q] = lu(sparse(AB));
   rhs  = P*AN;
   Lhat = L(1:n,:);
   W = Q*( U \ (Lhat \ rhs(1:n,:)));
   resnorm = norm(AN-AB*W,'fro')/max(1,norm(AN,'fro'));
%%*****************************************************************************

author: yu-jiang

repo: radpbook

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


text:
function volume = cuboidVolume(bb)

dis = (bb([1 2 5 6],:)-bb([3 4 3 4],:)).^2;

volume = (bb(10,:)-bb(9,:)).*sqrt((dis(1,:)+dis(2,:)).*(dis(3,:)+dis(4,:)));

author: zouchuhang

repo: LayoutNet


text:
%features
% load 'facefeature_lw6lr1_512.mat'
% features1=features;

%mirror feature
% load 'Mfacefeature_lw6lr1_512.mat'
% features=[features,features1];

tgt_method='facefeature_512'
tgtDir='./LFWTest/';

load pair_lfw
features=double(feature)';
features= bsxfun(@rdivide, features, sqrt(sum(features.^2,2)));

image_path = list;

for i=1:10
	intra_pairs=pair_lfw.IntraPersonPair{i};
	extra_pairs=pair_lfw.ExtraPersonPair{i};

        train_features=features;
        train_features([intra_pairs(:);extra_pairs(:)],:)=[];
	[coeff,~,latent] = pca(double(train_features)');
        m_feature=mean(double(train_features));
        accuracies = zeros(10,1);

%test the result with the range of PCA components
%	for numDim=64:64:1024
    for numDim=256:256
    %	for numDim=79:104

        test_features1=double(features([intra_pairs(:,1);extra_pairs(:,1)],:));
        test_features2=double(features([intra_pairs(:,2);extra_pairs(:,2)],:));

        test_features1=(test_features1-ones([size(test_features1,1) 1])*m_feature)*coeff(:,1:numDim);
        test_features2=(test_features2-ones([size(test_features2,1) 1])*m_feature)*coeff(:,1:numDim);

        test_features1= bsxfun(@rdivide, test_features1, sqrt(sum(test_features1.^2,2)));
        test_features2= bsxfun(@rdivide, test_features2, sqrt(sum(test_features2.^2,2)));

        scores=diag(test_features1*test_features2');

        test_path1=image_path([intra_pairs(:,1);extra_pairs(:,1)]);
        test_path2=image_path([intra_pairs(:,2);extra_pairs(:,2)]);
        [tgtDir, '/',tgt_method, '_', num2str(numDim),  '/TR',num2str(i)]
        mkdir([tgtDir, '/',tgt_method, '_', num2str(numDim)]);
        mkdir([tgtDir, '/',tgt_method, '_', num2str(numDim),  '/TR',num2str(i)]);

        fid=fopen([tgtDir, '/' ,tgt_method, '_', num2str(numDim), '/TR',num2str(i),'/NCMNCM52evaC.dat'],'w');
        for j=1:length(scores)
            idx=0;
            if ~isempty(strmatch(test_path1{j}(1:end-9),test_path2{j}(1:end-9)))
                idx=1;
            end
            fprintf(fid,'%s\t%s\t1\t%d\t%f\n',test_path1{j}(1:end-4),test_path2{j}(1:end-4),idx,-scores(j));
        end
        fclose(fid);
        accuracy = Match_measure([tgtDir, '/' ,tgt_method, '_', num2str(numDim), '/TR',num2str(i),'/NCMNCM52evaC.dat'])
        accuracies(i,1) = accuracy;
    end
end
mean(accuracy)

author: happynear

repo: FaceVerification


text:
classdef ORBoost < Algorithm
    %ORBoost ORBoost Boosting ensemble for Ordinal Regression [1]. This class
    %uses orensemble implementation at http://www.work.caltech.edu/~htlin/program/orensemble/
    %
    %   ORBoost methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   References:
    %
    %     [1] H.-T. Lin and L. Li, Large-margin thresholded ensembles for
    %         ordinal regression: Theory and practice, in Proc. of the 17th
    %         Algorithmic Learning Theory International Conference, ser.
    %         Lecture Notes in Artificial Intelligence (LNAI), J. L. Balcazar,
    %         P. M. Long, and F. Stephan, Eds., vol. 4264. Springer-Verlag,
    %         October 2006, pp. 319-333.
    %     [1] P.A. Gutiérrez, M. Pérez-Ortiz, J. Sánchez-Monedero,
    %         F. Fernández-Navarro and C. Hervás-Martínez
    %         Ordinal regression methods: survey and experimental study
    %         IEEE Transactions on Knowledge and Data Engineering, Vol. 28. Issue 1
    %         2016
    %         http://dx.doi.org/10.1109/TKDE.2015.2457911
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html
    properties
        description = 'OR Ensemble with perceptrons';
        parameters = [];
        weights = true;
    end

    methods

        function obj = ORBoost(varargin)
            %ORBoost constructs an object of the class ORBoost and sets its default
            %   characteristics
            %   OBJ = ORBoost() builds ORBoost object
            obj.parseArgs(varargin);
        end

        function [projectedTrain, predictedTrain] = privfit(obj,train,parameters)
            %PRIVFIT trains the model for the ORBoost method with TRAIN data and
            %vector of parameters PARAMETERS.

            % Output model file
            modelFile = tempname();
            % Write train file
            trainFile = tempname();
            dlmwrite(trainFile,[train.patterns train.targets],'delimiter',' ','precision',10);

            % Prepare command line
            if ispc
                bin_train = fullfile('Algorithms','orensemble', 'boostrank-train.exe');
                execute_train = sprintf('%s %s %d %d %d1 204 %d 2000 %s',...
                    bin_train, trainFile,size(train.patterns,1),size(train.patterns,2),(3+obj.weights),max(unique(train.targets)),modelFile);
            else
                execute_train = sprintf('%s/orensemble/hack.sh %s/orensemble/boostrank-train %s %d %d %d1 204 %d 2000 %s',...
                    fileparts(which('Algorithm.m')), fileparts(which('Algorithm.m')),...
                    trainFile,size(train.patterns,1),size(train.patterns,2),(3+obj.weights),max(unique(train.targets)),modelFile);
            end

            % Execute train
            system(execute_train);

            % Extract contents of model file
            fid = fopen(modelFile);
            if ~(exist ('OCTAVE_VERSION', 'builtin') > 0) && verLessThan('matlab','8.4')
                s = textscan(fid,'%s','Delimiter','\n','bufsize', 2^18-1);
            else
                s = textscan(fid,'%s','Delimiter','\n');
            end
            s = s{1};
            fclose(fid);

            model.textInformation = s;
            model.weights = obj.weights;
            obj.model = model;
            [projectedTrain,predictedTrain] = obj.predict(train.patterns);
            % Delete temp files
            delete(trainFile);
            delete(modelFile);
        end

        function [projected, predicted]= privpredict( obj,test)
            %PREDICT predicts labels of TEST patterns labels. The object needs to be fitted to the data first.

            % Write test file
            testFile = tempname();
            dlmwrite(testFile,[test ones(size(test,1),1)],'delimiter',' ','precision',10);
            % Write model file
            modelFile = tempname();
            fid = fopen(modelFile,'w');
            nLines = size(obj.model.textInformation,1);
            for iLine = 1:nLines
                fprintf(fid,'%s\n',obj.model.textInformation{iLine});
            end
            fclose(fid);
            % Name of prediction file
            predictFile = tempname();

            % Prepare command line
            if ispc
                bin_predict = fullfile(fileparts(which('Algorithm.m')),'orensemble', 'boostrank-predict.exe');
                execute_test = sprintf('%s %s %d %d %s 2000 %s',...
                    bin_predict,testFile,size(test,1),...
                    size(test,2),modelFile,predictFile);
            else
                execute_test = ...
                    sprintf('%s/orensemble/hack.sh %s/orensemble/boostrank-predict %s %d %d %s 2000 %s',...
                    fileparts(which('Algorithm.m')),fileparts(which('Algorithm.m')),testFile,size(test,1),size(test,2),modelFile,predictFile);
            end

            % Execute test
            system(execute_test);
            % Extract predictions (targets and projections)
            all = load(predictFile);
            predicted = all(:,1);
            projected = all(:,2);
            % Delete temp files
            delete(predictFile);
            delete(testFile);
            delete(modelFile);
        end

    end
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function dfsetbinwidthrules(dataset, v0, v1, v2, v3, v4, v5, v6)
% DFSETBINWIDTHRULES Helper function for the dfittool set default bin width panel

%   Copyright 2003-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.3 $  $Date: 2004/01/24 09:35:47 $

binDlgInfo.rule = v0;
binDlgInfo.nbinsExpr = v1;
binDlgInfo.nbins = str2num(v1);
binDlgInfo.widthExpr = v2;
binDlgInfo.width = str2num(v2);
binDlgInfo.placementRule = v3;
binDlgInfo.anchorExpr = v4;
binDlgInfo.anchor = str2num(v4);
binDlgInfo.applyToAll = v5;
binDlgInfo.setDefault = v6;

if (v5 == true) % apply to all
    dsdb = getdsdb;
    ds = down(dsdb);
    while(~isempty(ds))
        ds.binDlgInfo = binDlgInfo;
        ds = right(ds);
    end
    dfupdateallplots(true, false, true);
elseif ~isempty(dataset)
    ds = handle(dataset);
    ds.binDlgInfo = binDlgInfo;
    clearplot(ds);
    updateplot(ds);
end

% Update the main fig axis limits to fit the new histograms
dfupdatexlim;
dfupdateylim;

if (v6 == true) % set default
    dfgetset('binDlgInfo', binDlgInfo);
end


text:
function [mri] = ft_read_mri(filename, varargin)

% FT_READ_MRI reads anatomical and functional MRI data from different file formats.
% The output data is structured in such a way that it is compatible with
% FT_DATATYPE_VOLUME.
%
% Use as
%   [mri] = ft_read_mri(filename)
%
% Additional options should be specified in key-value pairs and can be
%   'dataformat' = string specifying the file format, determining the low-
%                  level reading routine to be used. If no explicit format
%                  is given, it is determined automatically from the filename.
%
% The following values apply for the dataformat
%   'afni_head'/'afni_brik'      uses AFNI code
%   'analyze_img'/'analyze_hdr'  uses SPM code
%   'analyze_old'                uses Darren Webber's code
%   'asa_mri'
%   'ctf_mri'
%   'ctf_mri4'
%   'ctf_svl'
%   'dicom'                      uses FreeSurfer code
%   'dicom_old'                  uses own code
%   'freesurfer_mgh'             uses FreeSurfer code
%   'freesurfer_mgz'             uses FreeSurfer code
%   'matlab'                     assumes a MATLAB *.mat file containing a mri structure according to FT_DATATYPE_VOLUME
%   'minc'                       uses SPM (<= version SPM5)
%   'neuromag_fif'               uses MNE toolbox
%   'neuromag_fif_old'           uses meg-pd toolbox
%   'nifti'                      uses FreeSurfer code
%   'nifti_fsl'                  uses FreeSurfer code
%   'nifti_spm'                  uses SPM
%   'yokogawa_mri'
%   'mrtrix_mif'                 uses mrtrix code
%
% The following MRI file formats are supported
%   CTF (*.svl, *.mri version 4 and 5)
%   NIFTi (*.nii) and zipped NIFTi (*.nii.gz)
%   Analyze (*.img, *.hdr)
%   DICOM (*.dcm, *.ima)
%   AFNI (*.head, *.brik)
%   FreeSurfer (*.mgz, *.mgh)
%   MINC (*.mnc)
%   Neuromag/Elekta (*.fif)
%   ANT - Advanced Neuro Technology (*.mri)
%   Yokogawa (*.mrk, incomplete)
%   Mrtrix image format (*.mif)
%
% If you have a series of DICOM files, please provide the name of any of the files
% in the series (e.g. the first one). The other files will be found automatically.
%
% The output MRI may have a homogenous transformation matrix that converts
% the coordinates of each voxel (in xgrid/ygrid/zgrid) into head
% coordinates.
%
% See also FT_DATATYPE_VOLUME, FT_WRITE_MRI, FT_READ_DATA, FT_READ_HEADER, FT_READ_EVENT

% Copyright (C) 2008-2013, Robert Oostenveld & Jan-Mathijs Schoffelen
%
% 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$

% optionally get the data from the URL and make a temporary local copy
filename = fetch_url(filename);

% get the options
dataformat = ft_getopt(varargin, 'dataformat');

% the following is added for backward compatibility of using 'format' rather than 'dataformat'
format    = ft_getopt(varargin, 'format');
if ~isempty(format)
  ft_warning('the option ''format'' will be deprecated soon, please use ''dataformat'' instead');
  if isempty(dataformat)
    dataformat  = format;
  end
end

if isempty(dataformat)
  % only do the autodetection if the format was not specified
  dataformat = ft_filetype(filename);
end

if strcmp(dataformat, 'compressed')
  % the file is compressed, unzip on the fly
  inflated = true;
  filename = inflate_file(filename);
  dataformat = ft_filetype(filename);
else
  inflated = false;
end

% test whether the file exists
if ~exist(filename, 'file')
  ft_error('file ''%s'' does not exist', filename);
end

% test for the presence of some external functions from other toolboxes
hasspm2  = ft_hastoolbox('spm2');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm5  = ft_hastoolbox('spm5');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm8  = ft_hastoolbox('spm8');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm12 = ft_hastoolbox('spm12');   % see http://www.fil.ion.ucl.ac.uk/spm/

switch dataformat
  case 'ctf_mri'
    [img, hdr] = read_ctf_mri(filename);
    transform = hdr.transformMRI2Head;
    coordsys  = 'ctf';

  case 'ctf_mri4'
    [img, hdr] = read_ctf_mri4(filename);
    transform = hdr.transformMRI2Head;
    coordsys  = 'ctf';

  case 'ctf_svl'
    [img, hdr] = read_ctf_svl(filename);
    transform = hdr.transform;

  case 'asa_mri'
    [img, seg, hdr] = read_asa_mri(filename);
    transform = hdr.transformMRI2Head;
    if isempty(seg)
      % in case seg exists it will be added to the output
      clear seg
    end

  case 'minc'
    if ~(hasspm2 || hasspm5)
      fprintf('the SPM2 or SPM5 toolbox is required to read *.mnc files\n');
      ft_hastoolbox('spm2', 1);
    end
    % use the functions from SPM
    hdr = spm_vol_minc(filename);
    img = spm_read_vols(hdr);
    transform = hdr.mat;

  case 'nifti_spm'
    if ~(hasspm5 || hasspm8 || hasspm12)
      fprintf('the SPM5 or newer toolbox is required to read *.nii files\n');
      ft_hastoolbox('spm12', 1);
    end
    % use the functions from SPM
    hdr = spm_vol_nifti(filename);
    img = double(hdr.private.dat);
    %img = spm_read_vols(hdr);
    transform = hdr.mat;

  case {'analyze_img' 'analyze_hdr'}
    if ~(hasspm8 || hasspm12)
      fprintf('the SPM8 or newer toolbox is required to read analyze files\n');
      ft_hastoolbox('spm8up', 1);
    end

    % use the image file instead of the header
    filename((end-2):end) = 'img';
    % use the functions from SPM to read the Analyze MRI
    hdr = spm_vol(filename);
    img = spm_read_vols(hdr);
    transform = hdr.mat;

  case 'analyze_old'
    % use the functions from Darren Weber's mri_toolbox to read the Analyze MRI
    ft_hastoolbox('mri', 1);     % from Darren Weber, see http://eeg.sourceforge.net/

    avw = avw_img_read(filename, 0); % returned volume is LAS*
    img = avw.img;
    hdr = avw.hdr;
    % The default Analyze orientation is axial unflipped (LAS*), which means
    % that the resulting volume is according to the radiological convention.
    % Most other fMRI and EEG/MEG software (except Mayo/Analyze) uses
    % neurological conventions and a right-handed coordinate system, hence
    % the first axis of the 3D volume (right-left) should be flipped to make
    % the coordinate system comparable to SPM
    ft_warning('flipping 1st dimension (L-R) to obtain volume in neurological convention');
    img = flipdim(img, 1);

    transform      = diag(hdr.dime.pixdim(2:4));
    transform(4,4) = 1;

  case {'afni_brik' 'afni_head'}
    % needs afni
    ft_hastoolbox('afni', 1);    % see http://afni.nimh.nih.gov/

    [err, img, hdr, ErrMessage] = BrikLoad(filename);
    if err
      ft_error('could not read AFNI file');
    end

    % FIXME: this should be checked, but I only have a single BRIK file
    % construct the homogenous transformation matrix that defines the axes
    ft_warning('homogenous transformation might be incorrect for AFNI file');
    transform        = eye(4);
    transform(1:3,4) = hdr.ORIGIN(:);
    transform(1,1)   = hdr.DELTA(1);
    transform(2,2)   = hdr.DELTA(2);
    transform(3,3)   = hdr.DELTA(3);

    % FIXME: I am not sure about the "RAI" image orientation
    img = flipdim(img,1);
    img = flipdim(img,2);
    dim = size(img);
    transform(1,4) = -dim(1) - transform(1,4);
    transform(2,4) = -dim(2) - transform(2,4);

  case 'neuromag_fif'
    % needs mne toolbox
    ft_hastoolbox('mne', 1);

    % use the mne functions to read the Neuromag MRI
    hdr = fiff_read_mri(filename);
    img_t = cat(3, hdr.slices.data);
    img = permute(img_t,[2 1 3]);
    hdr.slices = rmfield(hdr.slices, 'data'); % remove the image data to save memory

    % information below is from MNE - fiff_define_constants.m
    % coordinate system 4 - is the MEG head coordinate system (fiducials)
    % coordinate system 5 - is the MRI coordinate system
    % coordinate system 2001 - MRI voxel coordinates
    % coordinate system 2002 - Surface RAS coordinates (is mainly vertical
    %                                     shift, no rotation to 2001)
    % MEG sensor positions come in system 4
    % MRI comes in system 2001

    transform = eye(4);
    if isfield(hdr, 'trans') && issubfield(hdr.trans, 'trans')
      if (hdr.trans.from == 4) && (hdr.trans.to == 5)
        transform = hdr.trans.trans;
      else
        ft_warning('W: trans does not transform from 4 to 5.');
        ft_warning('W: Please check the MRI fif-file');
      end
    else
      ft_warning('W: trans structure is not defined.');
      ft_warning('W: Maybe coregistration is missing?');
    end
    if isfield(hdr, 'voxel_trans') && issubfield(hdr.voxel_trans, 'trans')
      % centers the coordinate system
      % and switches from mm to m
      if (hdr.voxel_trans.from == 2001) && (hdr.voxel_trans.to == 5)
        % matlab_shift compensates for the different index conventions
        % between C and matlab

        % the lines below is old code (prior to Jan 3, 2013) and only works with
        % 1 mm resolution MRIs
        %matlab_shift = [ 0 0 0 0.001; 0 0 0 -0.001; 0 0 0 0.001; 0 0 0 0];
        % transform transforms from 2001 to 5 and further to 4
        %transform = transform\(hdr.voxel_trans.trans+matlab_shift);

        % the lines below should work with arbitrary resolution
        matlab_shift = eye(4);
        matlab_shift(1:3,4) = [-1,-1,-1];
        transform = transform\(hdr.voxel_trans.trans * matlab_shift);

        coordsys  = 'neuromag';
        mri.unit  = 'm';
      else
        ft_warning('W: voxel_trans does not transform from 2001 to 5.');
        ft_warning('W: Please check the MRI fif-file');
      end
    else
      ft_warning('W: voxel_trans structure is not defined.');
      ft_warning('W: Please check the MRI fif-file');
    end

  case 'neuromag_fif_old'
    % needs meg_pd toolbox
    ft_hastoolbox('meg-pd', 1);

    % use the meg_pd functions to read the Neuromag MRI
    [img,coords] = loadmri(filename);
    dev = loadtrans(filename,'MRI','HEAD');
    transform  = dev*coords;
    hdr.coords = coords;
    hdr.dev    = dev;

  case 'dicom'
    % this seems to return a right-handed volume with the transformation
    % matrix stored in the file headers.

    % needs the freesurfer toolbox
    ft_hastoolbox('freesurfer', 1);
    [dcmdir,junk1,junk2] = fileparts(filename);
    if isempty(dcmdir),
      dcmdir = '.';
    end
    [img,transform,hdr,mr_params] = load_dicom_series(dcmdir,dcmdir,filename);
    transform = vox2ras_0to1(transform);

  case 'dicom_old'
    % this does not necessarily return a right-handed volume and only a
    % transformation-matrix with the voxel size

    % this uses the Image processing toolbox
    % the DICOM file probably represents a stack of slices, possibly even multiple volumes
    orig = dicominfo(filename);
    dim(1) = orig.Rows;
    dim(2) = orig.Columns;

    [p, f] = fileparts(filename);

    % this works for the Siemens scanners at the FCDC
    tok = tokenize(f, '.');
    for i=5:length(tok)
      tok{i} = '*';
    end
    filename = sprintf('%s.', tok{:});  % reconstruct the filename with wildcards and '.' between the segments
    filename = filename(1:end-1);       % remove the last '.'
    dirlist  = dir(fullfile(p, filename));
    dirlist  = {dirlist.name};

    if isempty(dirlist)
      % this is for the Philips data acquired at KI
      ft_warning('could not determine list of dicom files, trying with *.dcm');
      dirlist  = dir(fullfile(p, '*.dcm'));
      dirlist  = {dirlist.name};
    end

    if isempty(dirlist)
      ft_warning('could not determine list of dicom files, trying with *.ima');
      dirlist  = dir(fullfile(p, '*.ima'));
      dirlist  = {dirlist.name};
    end

    if length(dirlist)==1
      % try something else to get a list of all the slices
      dirlist = dir(fullfile(p, '*'));
      dirlist = {dirlist(~[dirlist.isdir]).name};
    end

    keep = false(1, length(dirlist));
    for i=1:length(dirlist)
      filename = char(fullfile(p, dirlist{i}));
      if ~strcmp(dataformat, 'dicom_old')
        keep(i) = false;
        fprintf('skipping ''%s'' because of incorrect filetype\n', filename);
      end
      % read the header information
      info     = dicominfo(filename);
      if info.SeriesNumber~=orig.SeriesNumber
        keep(i) = false;
        fprintf('skipping ''%s'' because of different SeriesNumber\n', filename);
      else
        keep(i) = true;
        hdr(i)  = info;
      end
    end
    % remove the files that were skipped
    hdr     = hdr(keep);
    dirlist = dirlist(keep);

    % pre-allocate enough space for the subsequent slices
    dim(3) = length(dirlist);
    img    = zeros(dim(1), dim(2), dim(3));
    for i=1:length(dirlist)
      filename = char(fullfile(p, dirlist{i}));
      fprintf('reading image data from ''%s''\n', filename);
      img(:,:,i) = dicomread(hdr(i));
    end

    % reorder the slices
    [z, indx]   = sort(cell2mat({hdr.SliceLocation}));
    hdr = hdr(indx);
    img = img(:,:,indx);

    try
      % construct a homgeneous transformation matrix that performs the scaling from voxels to mm
      dx = hdr(1).PixelSpacing(1);
      dy = hdr(1).PixelSpacing(2);
      dz = hdr(2).SliceLocation - hdr(1).SliceLocation;
      transform = eye(4);
      transform(1,1) = dx;
      transform(2,2) = dy;
      transform(3,3) = dz;
    end

  case {'nifti', 'freesurfer_mgz', 'freesurfer_mgh', 'nifti_fsl'}
    if strcmp(dataformat, 'freesurfer_mgz') && ispc
      ft_error('Compressed .mgz files cannot be read on a PC');
    end

    ft_hastoolbox('freesurfer', 1);
    tmp = MRIread(filename);
    ndims = numel(size(tmp.vol));
    if ndims==3
      img = permute(tmp.vol, [2 1 3]); %FIXME although this is probably correct
      %see the help of MRIread, anecdotally columns and rows seem to need a swap
      %in order to match the transform matrix (alternatively a row switch of the
      %latter can be done)
    elseif ndims==4
      img = permute(tmp.vol, [2 1 3 4]);
    end
    hdr = rmfield(tmp, 'vol');
    transform = tmp.vox2ras1;

  case 'yokogawa_mri'
    ft_hastoolbox('yokogawa', 1);
    fid = fopen_or_error(filename, 'rb');
    mri_info = GetMeg160MriInfoM(fid);
    patient_info = GetMeg160PatientInfoFromMriFileM(fid);
    [data_style, model, marker, image_parameter, normalize, besa_fiducial_point] = GetMeg160MriFileHeaderInfoM(fid);
    fclose(fid);

    % gather all meta-information
    hdr.mri_info = mri_info;
    hdr.patient_info = patient_info;
    hdr.data_style = data_style;
    hdr.model = model;
    hdr.marker = marker;
    hdr.image_parameter = image_parameter;
    hdr.normalize = normalize;
    hdr.besa_fiducial_point = besa_fiducial_point;

    ft_error('FIXME yokogawa_mri implementation is incomplete');

  case 'matlab'
    mri = loadvar(filename, 'mri');

  case {'mif' 'mrtrix_mif'}
    ft_hastoolbox('mrtrix', 1);
    tmp = read_mrtrix(filename);

    mri.anatomy = tmp.data;
    mri.dim     = tmp.dim;
    mri.transform = tmp.transform;
    mri.transform(1:3,1:3) = diag(tmp.vox(1:3))*mri.transform(1:3,1:3);

  otherwise
    ft_error('unrecognized filetype ''%s'' for ''%s''', dataformat, filename);
end

if exist('img', 'var')
  % set up the axes of the volume in voxel coordinates
  %nx = size(img,1);
  %ny = size(img,2);
  %nz = size(img,3);
  mri.dim = size(img); %[nx ny nz];
  % store the anatomical data
  mri.anatomy = img;
end

if exist('seg', 'var')
  % store the segmented data
  mri.seg = seg;
end

if exist('hdr', 'var')
  % store the header with all file format specific details
  mri.hdr = hdr;
end

try
  % store the homogenous transformation matrix if present
  mri.transform = transform;
end

try
  % try to determine the units of the coordinate system
  mri = ft_determine_units(mri);
end

try
  % try to add a descriptive label for the coordinate system
  mri.coordsys = coordsys;
end

if inflated
  % compressed file has been unzipped on the fly, clean up
  delete(filename);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = loadvar(filename, varname)
var = whos('-file', filename);
if length(var)==1
  filecontent = load(filename); % read the one variable in the file, regardless of how it is called
  value       = filecontent.(var.name);
  clear filecontent
else
  filecontent = load(filename, varname);
  value       = filecontent.(varname);  % read the variable named according to the input specification
  clear filecontent
end

author: spm


text:
function n = gqu_order ( l )

%*****************************************************************************80
%
%% GQU_ORDER computes the order of a GQU rule from the level.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2012
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, integer L, the level of the rule.
%    1 <= L.
%
%    Output, integer N, the order of the rule.
%
  if ( l < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GQU_ORDER - Fatal error!\n' );
    fprintf ( 1, '  1 <= L required.\n' );
    fprintf ( 1, '  Input L = %d\n', l );
    error (  'GQU_ORDER - Fatal error!' );
  elseif ( l <= 25 )
    n = l;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GQU_ORDER - Fatal error!\n' );
    fprintf ( 1, '  L <= 25 required.\n' );
    fprintf ( 1, '  Input L = %d\n', l );
    error (  'GQU_ORDER - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function prob_test132 ( )

%*****************************************************************************80
%
%% TEST132 tests RAYLEIGH_MEAN, RAYLEIGH_SAMPLE, RAYLEIGH_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST132\n' );
  fprintf ( 1, '  For the Rayleigh PDF:\n' );
  fprintf ( 1, '  RAYLEIGH_MEAN computes the mean;\n' );
  fprintf ( 1, '  RAYLEIGH_SAMPLE samples;\n' );
  fprintf ( 1, '  RAYLEIGH_VARIANCE computes the variance.\n' );

  a = 2.0;

  check = rayleigh_check ( a );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST132 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = rayleigh_mean ( a );
  variance = rayleigh_variance ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = rayleigh_sample ( a, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function value = r8mat_max ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_MAX returns the maximum entry of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the M by N matrix.
%
%    Output, real VALUE, the maximum entry of A.
%
  value = max ( max ( a(1:m,1:n) ) );

  return
end


text:
function polygon_properties_test ( )

%*****************************************************************************80
%
%% POLYGON_PROPERTIES_TEST tests the POLYGON_PROPERTIES library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_PROPERTIES_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the POLYGON_PROPERTIES library.\n' );

  polygon_properties_test01 ( );
  polygon_properties_test02 ( );
  polygon_properties_test03 ( );
  polygon_properties_test04 ( );
  polygon_properties_test05 ( );
  polygon_properties_test06 ( );
  polygon_properties_test07 ( );
  polygon_properties_test08 ( );
  polygon_properties_test09 ( );

  polygon_properties_test10 ( );
  polygon_properties_test11 ( );
  polygon_properties_test12 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_PROPERTIES_PRB\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%% OPT01_RUN
%
%  Modified:
%
%    09 January 2008
%
   %---------------------------------------------------------------------
   %  Running the first testcase, from D+S, pp. 100, 202; x_* = [2;-1].
   %  This should finish in 10 iterations for Newton's method and in
   %  12 iterations for BFGS.
   %---------------------------------------------------------------------

   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_1:  exact solution (2,-1)\n')
   fprintf('---------------------------------------------------------\n')

   fname = 'opt01_fgh';

   options = [];
   options.verbose            = 0;
   options.step_tolerance     = 1.e-11;
   options.gradient_tolerance = 1.e-11;
   options.max_iterations     = 15;

   x0 = [ 1; 1];
%
%  Run with Newton method.
%
   fprintf('Newton method:\n')
   options.method             = 'newton';
   options.globalization      = 'none';

   x = entrust(fname, x0, options);

   fprintf('Newton method produced (%10.7e,%10.7e)\n\n',x(1),x(2))
   f = opt01_fgh ( x, 'f' );
   fprintf('Value of F(X) = %f\n', f );
%
%  Run with secant method.
%
   fprintf('Secant method:\n')
   options.method             = 'secant';
   options.initial_hessian    = [ 14., -4. ; -4.  4. ];
   x = entrust(fname, x0, options);
   fprintf('BFGS   method produced (%10.7e,%10.7e)\n\n',x(1),x(2))
   f = opt01_fgh ( x, 'f' );
   fprintf('Value of F(X) = %f\n', f );
   %---------------------------------------------------------------------
   %  Test Gauss-Newton strategies.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_1 as least squares problem: \n')
   fprintf('Exact solution (2,-1)\n')
   fprintf('---------------------------------------------------------\n')

   fname = 'opt01_rj';
   options = [];
   options.verbose            = 0;
   options.method             = 'gauss_newton';
   options.step_tolerance     = 1.e-15;
   options.globalization      = 'none';
   options.gradient_tolerance = 1.e-10;
   options.max_iterations     = 40;

   x0 = [1,1];

   x = entrust(fname, x0, options );

   fprintf('Gauss-Newton produced  (%10.7e, %10.7e)\n\n',x(1),x(2))
   [ res, jac ] = opt01_rj ( x, 'f' );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );


text:
classdef SettingsHomogenizedVarComputer < AbstractSettings

    properties (Access = public)
        type
        designVariable
        dim
    end

    methods (Access = public, Static)

        function obj = create(cParams)
            factory = SettingsHomogenizedVarComputerFactory();
            obj     = factory.create(cParams);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% GEODEMO_4 Compare horizontal slices from two different CF compliant structured
% grid models (CH3D and ROMS) at a particular time step and depth
clear

%% CH3D salinity
url{1}='http://testbedapps.sura.org/thredds/dodsC/estuarine_hypoxia/ch3d/agg.nc';
var{1}='salinity';
titl{1}='CH3D';

%% ROMS salt
url{2}='http://testbedapps.sura.org/thredds/dodsC/estuarine_hypoxia/chesroms_1tdo/agg.nc';
var{2}='salt';
titl{2}='CHESROMS';

%% Choose a date,depth, and spatial extent
dat=[2005 1 10 0 0 0];  % Jan 10, 2005 00:00 UTC
depth=-5;  % horizontal slice 5 m from surface
ax=[  -77.4572  -75.4157   36.7224   39.6242]; %lon/lat range
cax=[0 33];  %color range
lat_mid=38; % for scaling plots

%% Perform analysis without using dataset dependant code
% Access datasets, get 3d field at given time data, interpolate data to a constant z, plot results at z depth

figure;
for i=1:length(url);
  nc{i}=ncgeodataset(url{i});
  jd{i}=nj_time(nc{i},var{i});  % using nj_time
  itime=date_index(jd{i},dat);
  disp(['reading data from ' titl{i} '...'])

  % using nj_tslice here, which doesn't allow for subsetting or
  % striding:  you get the whole 3D field at a particular time step:
  [s{i},g{i}]=nj_tslice(nc{i},var{i},itime);
  sz{i}=zsliceg(s{i},g{i}.z,depth);
  a{i}=subplot(1,length(url),i);
  pcolorjw(g{i}.lon,g{i}.lat,double(sz{i}));colorbar
  axis(ax);
  caxis(cax);
  title(sprintf('%s, depth=%f: %s',titl{i},depth,datestr(g{i}.time)));
  set (a{i}, 'DataAspectRatio', [1 cos(lat_mid*pi/180) 1000] );

end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function y = nansum(x,dim)
%Replacement for Matlab NANSUM Sum, ignoring NaNs.
%
x(isnan(x)) = 0;
if nargin == 1 % let sum figure out which dimension to work along
    y = sum(x);
else           % work along the explicitly given dimension
    y = sum(x,dim);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [uV] = eeg_volt2uv(Volts);

% eeg_volt2uv - Convert volts to microvolts
%
%   [uV] = eeg_volts2uv(Volts)
%
%   Simply, uV = Volts .* 10^6
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:50 $

% Licence:  GNU GPL, no implied or express warranties
% History:  11/2001, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uV = Volts .* 10^6;

author: PatternRecognition

repo: OpenBMI


text:
function [metabolite_structure_rBioNet] = createrBioNetStructure(metabolite_structure_rBioNet,start,stop)
% for the moment rBioNet will serve as universal database, but I will
% update this
warning off;
mkdir('data/');
currentPath = pwd;

if ~exist('metabolite_structure_rBioNet','var')
    websave('data/MetaboliteDatabase.txt','https://raw.githubusercontent.com/opencobra/COBRA.papers/master/2021_demeter/input/MetaboliteDatabase.txt');
    websave('data/ReactionDatabase.txt','https://raw.githubusercontent.com/opencobra/COBRA.papers/master/2021_demeter/input/ReactionDatabase.txt');
   % websave('data/compartments.mat','https://raw.githubusercontent.com/opencobra/COBRA.papers/master/2021_demeter/input/compartments.mat');
   % cd data;
    createRBioNetDBFromVMHDB('rBioNetDBFolder','data');
    load('data/metab.mat');
    cd(currentPath);
end

molFileDirectory = 'C:\Users\0123322S\Documents\GitHub\chemTableFiles\ctf\mets\molFiles';

annotationType = 'automatic';

if ~exist('metabolite_structure_rBioNet','var')
    fprintf('Create metabolite structure \n');
    % create metabolite structure out of metab
    metab2={'VMHId' 'metNames' 'neutralFormula' 'chargedFormula' 'charge' 'keggId' 'pubChemID' 'cheBIId' 'inchiKey' 'smiles' 'hmdb','subsystem'};
    metab(:,end-1)=[];%remove date
    metab2 = [metab2;metab];
    [metabolite_structure_rBioNet] =createNewMetaboliteStructure(metab2,'rBioNet');

    %[metabolite_structure_rBioNet,hit] = searchMultipleUnknownMetOnline(metabolite_structure_rBioNet);

    % these files run at this point as they are all offline, loading files
    % deposited in /data/. It would take more time to open these files for each
    % metabolite.
    [IDsStart,IDcountStart,TableStart] = getStatsMetStruct(metabolite_structure_rBioNet);
end

F = fieldnames(metabolite_structure_rBioNet);
if ~exist('start','var')
    start = 1;
end
if ~exist('stop','var')
    stop = length(F);
end

if start == 1
% apparently rBioNet does not contain all recon metabolites, so I will add
% them here.
% However, I will retrieve VMH information on the metabolites at the next
% stage
% this function also assigns the presence/absence in AGORA/Recon to enable
% quick identification of currently human or microbe only metabolites
fprintf('Assign AGORA Recon Presence \n');
[metabolite_structure_rBioNet] = assignAGORAReconPresence(metabolite_structure_rBioNet);
fprintf('Metabolon mapping \n');
[metabolite_structure_rBioNet] = VMH2Metabolon(metabolite_structure_rBioNet);
fprintf('Add manually refined/corrected annotations \n');
[metabolite_structure_rBioNet] = replaceVMHIds(metabolite_structure_rBioNet);
fprintf('Assign Seed IDs \n');
[metabolite_structure_rBioNet] = VMH2Seed(metabolite_structure_rBioNet);
% offline file provided by seed/kbase
fprintf('Add Kegg IDs using Seed  \n');
[metabolite_structure_rBioNet] = getSeed2Kegg(metabolite_structure_rBioNet);
save met_strc_rBioNet_new_30_06_2022
end
inchiKey = 1;
smiles = 1;
formula = 1;

% do not query HMDB by name as it takes too much time right now.
retrievePotHMDB = 0;
retrievePotHMDB2 = 0;

for i = start:stop
    startSearch = i;
    endSearch = i;
    i
    progress = i/stop;
    fprintf([num2str(progress) '  ... Annotating metabolites from different resources ... \n']);

    % obtain info from the VMH api, this also fills the gaps from earlier
    % where recon only metabolites were added without additional
    % information
    fprintf('Connect to VMH \n');
    [metabolite_structure_rBioNet] = parseVMH4IDs(metabolite_structure_rBioNet,startSearch,endSearch);

    fprintf('Add information from mol files \n');
    [metabolite_structure_rBioNet] = addInfoFromMolFiles(metabolite_structure_rBioNet,molFileDirectory,startSearch,endSearch);

    fprintf('Find HMDB Ids \n');
    [metabolite_structure_rBioNet] = getMissingHMDBMolForm(metabolite_structure_rBioNet,molFileDirectory,retrievePotHMDB,startSearch,endSearch);

    fprintf('Collect information from different databases ... \n');
    [metabolite_structure_rBioNet] = parseDBCollection(metabolite_structure_rBioNet,startSearch,endSearch);

    fprintf('Collect more mol files \n');
    [metabolite_structure_rBioNet] = getMolFilesMultipleSources(metabolite_structure_rBioNet,molFileDirectory,startSearch,endSearch);

    fprintf('Generate database independent id from mol files \n');
    [metabolite_structure_rBioNet] = generateInchiFromMol(metabolite_structure_rBioNet,molFileDirectory, inchiKey, smiles,formula,startSearch,endSearch);

    fprintf('Collect more mol files - 2 \n');
    % now repeat the seearch for mol files again but skip the guessing of hmdb
    % ids from name

    [metabolite_structure_rBioNet] = getMissingHMDBMolForm(metabolite_structure_rBioNet,molFileDirectory,retrievePotHMDB2,startSearch,endSearch);
    [metabolite_structure_rBioNet] = getMissingDrugMolForm(metabolite_structure_rBioNet,molFileDirectory,startSearch,endSearch);
    [metabolite_structure_rBioNet] = getMolFilesMultipleSources(metabolite_structure_rBioNet,molFileDirectory,startSearch,endSearch);
    fprintf('Generate database independent id from mol files - 2 \n');

    [metabolite_structure_rBioNet] = generateInchiFromMol(metabolite_structure_rBioNet,molFileDirectory, inchiKey, smiles,formula,startSearch,endSearch);
    fprintf('Assign metabolite classification \n');

    [metabolite_structure_rBioNet] = assignClassyFire(metabolite_structure_rBioNet,startSearch,endSearch);
    %fprintf('Check that HMDB Ids are valid \n');
    % works only for HMDB right now
    % do it later, the last time only 16 out of 7900 were wrong - it takes
    % too much time
    % [metabolite_structure_rBioNet] = checkLinkValidity(metabolite_structure_rBioNet,startSearch,endSearch);
    if mod(i,10)==1
        save met_strc_rBioNet_new_30_06_2022
    end
end

% offline files
fprintf('Collecting information from Echa \n')
 [metabolite_structure_rBioNet] = getCas2Echa(metabolite_structure_rBioNet);
fprintf('Collecting information from CTD \n')
[metabolite_structure_rBioNet] = getCas2CTD(metabolite_structure_rBioNet);

% getIDsFromBIGG; --> I should use the online version --> neither one is
% accounted for in the last rBioNetX hence the final numbers are still
% lower

% sort the fields in the structure
F = fieldnames(metabolite_structure_rBioNet);
for i = 1 : length(F)
metabolite_structure_rBioNet.(F{i}) = orderfields(metabolite_structure_rBioNet.(F{i}));
end
[IDsEnd2,IDcountEnd2,TableEnd2] = getStatsMetStruct(metabolite_structure_rBioNet);
save met_strc_rBioNet_new_07_10_2021
writecell(TableEnd2,'rbionetMetStructure.xls');

author: opencobra

repo: cobratoolbox


text:
function plot_recall_ILSVRC2013()
  % Plot the recall of ILSVRC2013 validation set ground truth for all methods.
  %
  % This function requires the proposals to already be saved to disk. It will
  % compute a matching between ground truth and proposals (if the result is not
  % yet found on disk) and then plot recall curves. The plots are saved to
  % figures/.

  val = load('data/ILSVRC2013_val_annotations.mat');
  methods = get_method_configs();
  methods([14 16 19:24]) = [];

  valid_methods = compute_best_candidates(val, methods);
  methods = methods(valid_methods);

  plot_legend(methods);
  printpdf('figures/imagenet_recall_legend.pdf');

  fprintf('\n\n');
  for i = 1:numel(methods)
    fprintf('Valid data for plots: %s\n', methods(i).name);
  end

  fh = figure;
  plot_overlap_recall_curve({methods.best_imagenet_candidates_file}, methods, 100, fh, true, 'none');
  hei = 10;
  wid = 10;
  set(gcf, 'Units','centimeters', 'Position',[0 0 wid hei]);
  set(gcf, 'PaperPositionMode','auto');
  printpdf('figures/imagenet_recall_100.pdf')

  fh = figure;
  plot_overlap_recall_curve({methods.best_imagenet_candidates_file}, methods, 1000, fh, false, 'none');
  hei = 10;
  wid = 10;
  set(gcf, 'Units','centimeters', 'Position',[0 0 wid hei]);
  set(gcf, 'PaperPositionMode','auto');
  printpdf('figures/imagenet_recall_1000.pdf')

  fh = figure;
  plot_overlap_recall_curve({methods.best_imagenet_candidates_file}, methods, 10000, fh, false, 'none');
  hei = 10;
  wid = 10;
  set(gcf, 'Units','centimeters', 'Position',[0 0 wid hei]);
  set(gcf, 'PaperPositionMode','auto');
  printpdf('figures/imagenet_recall_10000.pdf')

  plot_num_candidates_auc({methods.best_imagenet_candidates_file}, methods, 'imagenet_');
end

function valid_methods = compute_best_candidates(testset, methods)
  num_annotations = numel(testset.pos);
  candidates_thresholds = round(10 .^ (0:0.5:4));
  num_candidates_thresholds = numel(candidates_thresholds);

  valid_methods = false(numel(methods), 1);
  for method_idx = 1:numel(methods)
    method = methods(method_idx);
    fprintf('%s\n', method.name);
    try
      load(method.best_imagenet_candidates_file, 'best_candidates');
      valid_methods(method_idx) = true;
      continue;
    catch
    end

    % preallocate
    best_candidates = [];
    best_candidates(num_candidates_thresholds).candidates_threshold = [];
    best_candidates(num_candidates_thresholds).best_candidates = [];
    for i = 1:num_candidates_thresholds
      best_candidates(i).candidates_threshold = candidates_thresholds(i);
      best_candidates(i).best_candidates.candidates = zeros(num_annotations, 4);
      best_candidates(i).best_candidates.iou = zeros(num_annotations, 1);
      best_candidates(i).image_statistics(numel(testset.impos)).num_candidates = 0;
    end

    files_missing = false;
    pos_range_start = 1;
    for j = 1:numel(testset.impos)
      tic_toc_print('evalutating %s: %d/%d\n', method.name, j, numel(testset.impos));
      pos_range_end = pos_range_start + size(testset.impos(j).boxes, 1) - 1;
      assert(pos_range_end <= num_annotations);

%       tic_toc_print('sampling candidates for image %d/%d\n', j, numel(testset.impos));
      [~,img_id,~] = fileparts(testset.impos(j).im);

      for i = 1:num_candidates_thresholds
        try
          [candidates, scores] = get_candidates(method, img_id, ...
            candidates_thresholds(i));
        catch
          fprintf('%s candidates for %s missing\n', method.name, img_id);
          files_missing = true;
          break;
        end
        if isempty(candidates)
          impos_best_ious = zeros(size(testset.impos(j).boxes,1),1);
          impos_best_boxes = zeros(size(testset.impos(j).boxes,1),4);
        else
          [impos_best_ious, impos_best_boxes] = closest_candidates(...
            testset.impos(j).boxes, candidates);
        end

        best_candidates(i).best_candidates.candidates(pos_range_start:pos_range_end,:) = impos_best_boxes;
        best_candidates(i).best_candidates.iou(pos_range_start:pos_range_end) = impos_best_ious;
        best_candidates(i).image_statistics(j).num_candidates = size(candidates, 1);
      end
      if files_missing
        break;
      end

      pos_range_start = pos_range_end + 1;
    end
    valid_methods(method_idx) = ~files_missing;
    if ~files_missing
      save(method.best_imagenet_candidates_file, 'best_candidates');
    end

  end
end

author: hosang

repo: detection-proposals


text:
function t = model_cmp(m1, m2)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

[v1, b1, map1] = model2blocks(m1);
[v2, b2, map2] = model2blocks(m2);

e = sum(abs(v1-v2));
if e > 0
  fprintf('error: %.14f\n', e);
  for i = 1:length(b1)
    I = find(b1{i} ~= b2{i});
    if ~isempty(I)
      for j = 1:length(I)
        fprintf('at %s : index %d [%.5f vs %.5f]\n', ...
                map1{i}, I(j), b1{i}(I(j)), b2{i}(I(j)));
      end
    end
  end
else
  fprintf('no error\n');
end

function [m, blocks, map] = model2blocks(model)

blocks = cell(model.numblocks, 1);

% filters
for i = 1:model.numfilters
  if model.filters(i).flip == 0
    bl = model.filters(i).blocklabel;
    w = my_get_block(model, model.filters(i));
    blocks{bl} = w(:);
    map{bl} = ['filter ' num2str(i)];
  end
end

% offsets
for i = 1:length(model.rules)
  for j = 1:length(model.rules{i})
    bl = model.rules{i}(j).offset.blocklabel;
    w = my_get_block(model, model.rules{i}(j).offset);
    blocks{bl} = w;
    map{bl} = ['offset rule ' num2str(i) ' ind ' num2str(j)];
  end
end

% deformation models
for i = 1:length(model.rules)
  for j = 1:length(model.rules{i})
    if model.rules{i}(j).type == 'D' && model.rules{i}(j).def.flip == 0
      bl = model.rules{i}(j).def.blocklabel;
      w = my_get_block(model, model.rules{i}(j).def);
      blocks{bl} = w(:);
      map{bl} = ['def rule ' num2str(i)];
    end
  end
end

% concatenate
m = [];
for i = 1:model.numblocks
  m = [m; blocks{i}];
end

function w = my_get_block(model, obj)

if isfield(model, 'blocks')
  w = model_get_block(model, obj);
else
  w = obj.w;
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function combo_test16 ( )

%*****************************************************************************80
%
%% COMBO_TEST16 tests MARRIAGE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 5;
%
%  PREFER(M,W) is the index of women W on man M's list.
%
  prefer = [ ...
    2, 1, 2, 1, 5; ...
    5, 2, 3, 3, 3; ...
    1, 3, 5, 2, 2; ...
    3, 4, 4, 4, 1; ...
    4, 5, 1, 5, 4 ]';
%
%  RANK(W,M) is the index of man M on woman W's list.
%
  rank = [ ...
    2, 4, 1, 4, 5; ...
    4, 3, 3, 2, 2; ...
    5, 5, 4, 1, 3; ...
    3, 1, 2, 3, 1; ...
    1, 2, 5, 5, 4 ]';

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST16\n' );
  fprintf ( 1, '  MARRIAGE arranges a set of stable marriages\n' );
  fprintf ( 1, '  given a set of preferences.\n' );

  [ fiancee, next ] = marriage ( n, prefer, rank );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Man, Wife''s rank, Wife\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d\n', i, next(i), prefer(i,next(i)) );
  end

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Woman, Husband''s rank, Husband\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d\n', i, rank(i,fiancee(i)), fiancee(i) );
  end

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Correct result:\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '  M:W 1  2  3  4  5\n' );
  fprintf ( 1, '   1  +  .  .  .  .\n' );
  fprintf ( 1, '   2  .  .  .  +  .\n' );
  fprintf ( 1, '   3  .  .  .  .  +\n' );
  fprintf ( 1, '   4  .  .  +  .  .\n' );
  fprintf ( 1, '   5  .  +  .  .  .\n' );

  return
end


text:
% This is a demo for segmentation using local gaussian distribution (LGD)
% fitting energy
%
% Reference: <Li Wang, Lei He, Arabinda Mishra, Chunming Li.
% Active Contours Driven by Local Gaussian Distribution Fitting Energy.
% Signal Processing, 89(12), 2009,p. 2435-2447>
%
% Please DO NOT distribute this code to anybody.
% Copyright (c) by Li Wang
%
% Author:       Li Wang
% E-mail:       li_wang@med.unc.edu
% URL:          http://www.unc.edu/~liwa/
%
% 2010-01-02 PM

clc;clear all;close all;

Img=imread('1.bmp');
Img = double(Img(:,:,1));

NumIter = 300; %iterations
timestep=0.1; %time step
mu=0.1/timestep;% level set regularization term, please refer to "Chunming Li and et al. Level Set Evolution Without Re-initialization: A New Variational Formulation, CVPR 2005"
sigma = 3;%size of kernel
epsilon = 1;
c0 = 2; % the constant value
lambda1=1.0;%outer weight, please refer to "Chunming Li and et al,  Minimization of Region-Scalable Fitting Energy for Image Segmentation, IEEE Trans. Image Processing, vol. 17 (10), pp. 1940-1949, 2008"
lambda2=1.0;%inner weight
%if lambda1>lambda2; tend to inflate
%if lambda1<lambda2; tend to deflate
nu = 0.001*255*255;%length term
alf = 30;%data term weight

figure,imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal
[Height Wide] = size(Img);
[xx yy] = meshgrid(1:Wide,1:Height);
phi = (sqrt(((xx - 65).^2 + (yy - 40).^2 )) - 20);
phi = sign(phi).*c0;

Ksigma=fspecial('gaussian',round(2*sigma)*2 + 1,sigma); %  kernel
ONE=ones(size(Img));
KONE = imfilter(ONE,Ksigma,'replicate');
KI = imfilter(Img,Ksigma,'replicate');
KI2 = imfilter(Img.^2,Ksigma,'replicate');

figure,imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal,
hold on,[c,h] = contour(phi,[0 0],'r','linewidth',1); hold off
pause(0.5)

tic
for iter = 1:NumIter
    phi =evolution_LGD(Img,phi,epsilon,Ksigma,KONE,KI,KI2,mu,nu,lambda1,lambda2,timestep,alf);

    if(mod(iter,25) == 0)
        figure(2),
        imagesc(uint8(Img),[0 255]),colormap(gray),axis off;axis equal,title(num2str(iter))
        hold on,[c,h] = contour(phi,[0 0],'r','linewidth',1); hold off
        pause(0.02);
    end

end
toc

author: Sable


text:
function varargout = svd(varargin)
%SVD    Singular value decomposition of a CHEBFUN2.
%   SVD(F) returns the singular values of F. The number of singular values
%   returned is equal to the rank of the CHEBFUN2.
%
%   S = SVD(F) returns the singular values of F. S is a vector of singular
%   values in decreasing order.
%
%   [U, S, V] = SVD(F) returns the SVD of F. U and V are quasi-matrices of
%   orthogonal CHEBFUN objects and S is a diagonal matrix with the singular
%   values on the diagonal.
%
%   The length and rank of a CHEBFUN2 are slightly different quantities.
%   LENGTH(F) is the number of pivots used by the constructor, and
%   RANK(F) is the number of significant singular values of F. The relation
%   RANK(F) <= LENGTH(F) should always hold.

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

[varargout{1:nargout}] = svd@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function train_full_net(varargin)
% load data & label
data = load('./data/train/SDR_youtube_80.mat') ;
label = load('./data/train/HDR_youtube_80.mat') ;
imdb.images.data = data.SDR_data;
imdb.images.label = label.HDR_data;
imdb.images.set = cat(2, ones(1, size(data.SDR_data, 4)-500), 2*ones(1, 500));

% set CNN model
net = net_full();
netstruct = load('./net/net_base/net-epoch-200.mat');
net_init = dagnn.DagNN.loadobj(netstruct.net);

for i = 1:148
    % learning rate and weight decay for biases
    if mod(i, 2) == 0
        net.params(i).learningRate = 0.1;
        net.params(i).weightDecay = 0;
    end
    % initialize net_full with pre-trained values of net_base
    if i>6 && i< 17
        net.params(i).value = net_init.params(i-6).value;
    elseif i>20 && i<29
        net.params(i).value = net_init.params(i-10).value;
    elseif i>32 && i<41
        net.params(i).value = net_init.params(i-14).value;
    elseif i>50 && i<63
        net.params(i).value = net_init.params(i-24).value;
    elseif i>66 && i<79
        net.params(i).value = net_init.params(i-28).value;
    elseif i>82 && i<95
        net.params(i).value = net_init.params(i-32).value;
    elseif i>98
        net.params(i).value = net_init.params(i-36).value;
    end
end
net.conserveMemory = true;

% options
opts.solver = @adam;
opts.train.batchSize = 16;
opts.train.continue = false;
opts.train.gpus = 1;
opts.train.prefetch = false ;
opts.train.expDir = './net/net_full' ;
opts.train.learningRate = [1e-7*ones(1, 250) 1e-8*ones(1, 10) 1e-9*ones(1, 10)];
opts.train.weightDecay = 0.0005;
opts.train.numEpochs = numel(opts.train.learningRate) ;
opts.train.derOutputs = {'objective', 1} ;
[opts, ~] = vl_argparse(opts.train, varargin) ;

% record
if(~isdir(opts.expDir))
    mkdir(opts.expDir);
end

% call training function
[net,info] = cnn_train_dag(net, imdb, @getBatch, opts) ;

function inputs = getBatch(imdb, batch, opts)
image = imdb.images.data(:, :, :, batch) ;
label = imdb.images.label(:, :, :, batch) ;

image = single(image)/255;
label = single(label)/1023;
inputs = {'input', gpuArray(image), 'label', gpuArray(label)};

author: sooyekim

repo: Deep-SR-ITM


text:
% Fig. 8.8   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami

%    (requires fig8_08c.mdl)

clear all;
%close all;

% response comparison of continuous and digital control.

clear
clf
[tout,yout]=sim('fig8_08c');
r=[1 1];   %reference input
t=[0 2];
figure(1)
subplot(2,1,1)
plot(t,r,'r')
hold on
plot(ycd(:,1),ycd(:,2))
plot(ycd(:,1),ycd(:,3),'m')
title('Figure 8.8 Step Responses of Digital and Continuous Controllers')
ylabel('Position, y')
xlabel('Time (sec)')
text(.33,.9, '\leftarrow continuous  controller')
text(.7,1.3, 'digital  controller')
nicegrid
hold off
subplot(2,1,2)
plot(ycd(:,1),ycd(:,4))
hold on
plot(ycd(:,1),ycd(:,5),'m')
ylabel('Control, u')
xlabel('Time (sec)')
text(.55,10, ' continuous  controller')
text(.08,29, '\leftarrow digital  controller')
nicegrid
hold off

repo: mcbench-benchmarks


text:
function varargout = xyz2lab(varargin)
% VL_XYZ2LAB  Convert XYZ color space to LAB
%   J = VL_XYZ2LAB(I) converts the image from XYZ format to LAB format.
%
%   VL_XYZ2LAB(I,IL) uses one of the illuminants A, B, C, E, D50, D55,
%   D65, D75, D93. The default illuminatn is E.
%
%   See also: VL_XYZ2LUV(), VL_HELP().
[varargout{1:nargout}] = vl_xyz2lab(varargin{:});

author: yihui-he

repo: panorama

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


text:
%% MKL Numerical Jacobian Install for OPTI Toolbox
% Supplied binaries are built from MKL 10.3 Release 11

%   Copyright (C) 2011 Jonathan Currie (I2C2)

% This file will help you compile the Intel Math Kernel Library (MKL)
% djacobi function for use with MATLAB. NOTE you must NOT link the threaded
% MKL libraries as the MATLAB callback function is not thread safe!

% My build platform:
% - Windows 7 SP1 x64
% - Visual Studio 2010
% - Intel Math Kernel Library 10.3

% To recompile you will need to get / do the following:

% 1) Get and Install Intel MKL
% http://software.intel.com/en-us/articles/intel-mkl/

% 2) MKL JAC MEX Interface
% The MKL JAC MEX Interface is a simple MEX interface I wrote to use this
% function and is supplied in the Utilities\NumDiff folder.

% 6) Compile the MEX File
% The code below will automatically include all required libraries and
% directories to build the MKL JAC MEX file. Once you have completed all
% the above steps, simply run this file to compile! You MUST BE in the
% base directory of OPTI!

clear mklJac

% Modify below function if it cannot find Intel MKL on your system.
mkl_link = opti_FindMKL('seq'); %NOTE sequential only build!

fprintf('\n------------------------------------------------\n');
fprintf('MKL JAC MEX FILE INSTALL\n\n');

%Get MKL Includes & MKL Libraries (for DJACOBI)
post = mkl_link;

%CD to Source Directory
cdir = cd;
cd 'Utilities/Differentiation/Numerical/Source';

%Compile & Move
pre = 'mex -v -largeArrayDims mklJac.c';
try
    eval([pre post])
    movefile(['mklJac.' mexext],'../','f')
    fprintf('Done!\n');
catch ME
    cd(cdir);
    error('opti:mkljac','Error Compiling MKL JAC!\n%s',ME.message);
end
cd(cdir);
fprintf('------------------------------------------------\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [model, removedMets, removedRxns] = removeDeadEnds(model)
% Removes all dead end metabolites and reactions from the
% model
%
% USAGE:
%
%    [model, removedMets, removedRxns] = removeDeadEnds(model)
%
% INPUT:
%     model:          COBRA model structure
%
% OUTPUTS:
%    model:          COBRA model structure w/o dead end metabolites and
%                    reactions
%    removedMets:    List of removed metabolites
%    removedRxns:    List of removed reactions
%
% .. Author: - Markus Herrgard 8/29/06

removedMets = {};
removedRxns = {};

while (1)

    deadEnd = detectDeadEnds(model);

    if (length(deadEnd) == 0)
        break;
    end

    removedMets = union(removedMets,model.mets(deadEnd));
    if (length(deadEnd) > 1)
        deadRxns = model.rxns(find(sum(model.S(deadEnd,:) ~= 0) > 0));
    else
        deadRxns = model.rxns(find(model.S(deadEnd,:) ~= 0));
    end

    model = removeRxns(model,deadRxns);

    removedRxns = union(removedRxns,deadRxns);

end


text:
function updateTSDF(camRtC2W)

global ViewFrustumC;
ViewFrustumW = transformRT(ViewFrustumC, camRtC2W);

range2test = [min(ViewFrustumW, [], 2) max(ViewFrustumW, [], 2)];

global XYZcam;
global voxel;
global K;
global tsdf_value;
global tsdf_weight;

% choose a bounding box to contain viewing frustum
rangeGrid = (range2test - voxel.range(:,[1 1])) / voxel.unit + 1;
rangeGrid(:,1) = max(1,floor(rangeGrid(:,1)));
rangeGrid(:,2) = min(ceil (rangeGrid(:,2)),voxel.size_grid);
rangeGrid = int32(rangeGrid);

% get the grid there
disp('meshgrid');
tic;
[Y,X,Z]=meshgrid(rangeGrid(1,1):rangeGrid(1,2),rangeGrid(2,1):rangeGrid(2,2),rangeGrid(3,1):rangeGrid(3,2)); % strange matlab syntax
toc;
disp('grid to world');
tic;
X = X(:)'; Y = Y(:)'; Z = Z(:)';
gridIndex = sub2ind(voxel.size_grid',X,Y,Z);
gridCoordinateW = [single(X)*voxel.unit + voxel.range(1,1); single(Y)*voxel.unit + voxel.range(2,1); single(Z)*voxel.unit + voxel.range(3,1)];
clear X Y Z;
toc;

% transform the grid
disp('transform');
tic;
gridCoordinateC = transformRT(gridCoordinateW, camRtC2W, true);
toc;

% select: in front of camera
disp('select: in front of camera');
tic;
isValid = find(gridCoordinateC(3,:)>0);
gridCoordinateC = gridCoordinateC(:,isValid);
gridIndex = gridIndex(isValid);
toc;

% select: project
disp('select: project');
tic;
px = round(K(1,1)*(gridCoordinateC(1,:)./gridCoordinateC(3,:)) + K(1,3));
py = round(K(2,2)*(gridCoordinateC(2,:)./gridCoordinateC(3,:)) + K(2,3));
isValid = (1<=px & px <= 640 & 1<=py & py<= 480);
gridCoordinateC = gridCoordinateC(:,isValid);
gridIndex = gridIndex(isValid);
py = py(isValid);
px = px(isValid);
toc;

% select: valid depth
disp('select: valid depth');
tic;
ind = sub2ind([480 640],py,px);
isValid = XYZcam(ind+640*480*3)~=0;
gridCoordinateC = gridCoordinateC(:,isValid);
gridIndex = gridIndex(isValid);
ind = ind(isValid);
toc;

% compare distance between measurement and the grid
disp('compare distance between measurement and the grid');
%{

        const float diff = (depth[px] - cameraX.z) * sqrt(1+sq(pos.x/pos.z) + sq(pos.y/pos.z));
        if(diff > -mu){
            const float sdf = fminf(1.f, diff/mu);
            float2 data = vol[pix];
            data.x = clamp((data.y*data.x + sdf)/(data.y + 1), -1.f, 1.f);
            data.y = fminf(data.y+1, maxweight);
            vol.set(pix, data);
        }

%}
tic;
eta = (XYZcam(ind+640*480*2)- gridCoordinateC(3,:)) .* ((1+ (gridCoordinateC(1,:)./gridCoordinateC(3,:)).^2 + (gridCoordinateC(2,:)./gridCoordinateC(3,:)).^2 ).^0.5);
toc;

%tic;
%XYZmeasure = [XYZcam(ind);XYZcam(ind+640*480);XYZcam(ind+640*480*2)];
%Dmeasure = sqrt(sum(XYZmeasure.^2,1));
%Dtran    = sqrt(sum(gridCoordinateC.^2,1));
%eta = Dmeasure-Dtran;
%toc;

% select: > - mu
disp('select: > - mu');
tic;
isValid = eta>-voxel.mu;
eta = eta(isValid);
gridIndex = gridIndex(isValid);
ind = ind(isValid);
new_value = min(1,eta/voxel.mu);
toc;

disp('read write tsdf');
tic;
old_weight = tsdf_weight(gridIndex);
new_weight = old_weight + 1;
tsdf_weight (gridIndex)= new_weight;
tsdf_value (gridIndex) = (tsdf_value(gridIndex).*old_weight +new_value)./new_weight;
toc;

disp('coloring');
tic;
global tsdf_color;
global IMGcam;
if ~isempty(tsdf_color)
    isValid = eta<=voxel.mu;
    eta = eta(isValid);
    gridIndex = gridIndex(isValid);
    ind = ind(isValid);

    IMGmeasure = [IMGcam(ind);IMGcam(ind+640*480);IMGcam(ind+640*480*2)];

    colorArray = reshape(typecast(tsdf_color(gridIndex),'uint8'),4,[]);

    old_weight = colorArray(4,:);
    old_weight(double(old_weight)>=255) = 254;
    new_weight = old_weight + 1;

    newArray = uint8(round((double(colorArray(1:3,:)) .* repmat(double(old_weight),3,1) + double(IMGmeasure)) ./ repmat(double(new_weight),3,1)));
    newArray(4,:) = new_weight;

    tsdf_color(gridIndex) = typecast(newArray(:),'uint32');
end
toc;

%{
__global__ void integrate( Volume vol, const Image<float> depth, const Matrix4 invTrack, const Matrix4 K, const float mu, const float maxweight){
    uint3 pix = make_uint3(thr2pos2());
    float3 pos = invTrack * vol.pos(pix);
    float3 cameraX = K * pos;
    const float3 delta = rotate(invTrack, make_float3(0,0, vol.dim.z / vol.size.z));
    const float3 cameraDelta = rotate(K, delta);

    for(pix.z = 0; pix.z < vol.size.z; ++pix.z, pos += delta, cameraX += cameraDelta){
       if(pos.z < 0.0001f) // some near plane constraint
            continue;
        const float2 pixel = make_float2(cameraX.x/cameraX.z + 0.5f, cameraX.y/cameraX.z + 0.5f);
        if(pixel.x < 0 || pixel.x > depth.size.x-1 || pixel.y < 0 || pixel.y > depth.size.y-1)
            continue;
        const uint2 px = make_uint2(pixel.x, pixel.y);
        if(depth[px] == 0)
            continue;
        const float diff = (depth[px] - cameraX.z) * sqrt(1+sq(pos.x/pos.z) + sq(pos.y/pos.z));
        if(diff > -mu){
            const float sdf = fminf(1.f, diff/mu);
            float2 data = vol[pix];
            data.x = clamp((data.y*data.x + sdf)/(data.y + 1), -1.f, 1.f);
            data.y = fminf(data.y+1, maxweight);
            vol.set(pix, data);
        }
    }
}
%}

author: jianxiongxiao

repo: ProfXkit