text:
function vw = initMontageSlider(vw)
%
% vw = initMontageSlider(vw)
%
% initializes the slider that sets the size of
% an inplane montage.
%
% ras 09/04, off initSliceSlider
if ~strcmp(vw.viewType,'Inplane')
   error('initMontageSlider: Only used for INPLANE window');
end

ui = viewGet(vw,'ui');
sliderHandle = ui.montageSize.sliderHandle;

% Set the step size for clicking on the arrows and the trough
nSlices = viewGet(vw, 'numSlices');

% if there's only one slice, this slider should be turned off
if nSlices < 2
    nSlices=2;
    set(sliderHandle,'Enable','off');
else
    set(sliderHandle,'Enable','on');
end

% set the range and step size of the slider
set(sliderHandle,'Min',1,'Max',nSlices);
sliderStep = [1/(nSlices-1) , 3/(nSlices-1)];
set(sliderHandle,'sliderStep',sliderStep);

% Update the slider callback to be an integer.
sliderCb = 'val=round(get(gcbo,''Value''));';
sliderCb = sprintf('%s \n setSlider(%s,%s.ui.montageSize,val,0);',...
               sliderCb,vw.name,vw.name);
sliderCb = sprintf('%s \n %s=refreshScreen(%s);',...
                sliderCb,vw.name,vw.name);
set(vw.ui.montageSize.sliderHandle,'CallBack',sliderCb);

% Update the edit callback to be an integer.
editCb = 'val=round(str2num(get(gcbo,''String'')));';
editCb = sprintf('%s \n setSlider(%s,%s.ui.montageSize,val,0);',...
               editCb,vw.name,vw.name);
editCb = sprintf('%s \n %s=refreshScreen(%s);',...
                editCb,vw.name,vw.name);
set(vw.ui.montageSize.labelHandle,'CallBack',editCb);

% If we are first opening the window, create the text box.
if ~isfield(vw.ui.montageSize,'textHandle')

   % Position the text box to the right of the slider
   pos = get(sliderHandle,'Position');
   l2 = pos(1) + pos(3) + .01;
   w2 = pos(3)/4; b = pos(2); h = pos(4);
   position = [l2 b w2 h];

   vw.ui.montageSize.textHandle = ...
      uicontrol('Style','text',...
      'Units','normalized',...
      'Position',position);
end

% This text tells the user the max number of slices being
% shown in the montage (starting with the selected slice)
str = sprintf('%.0f',nSlices);
set(vw.ui.montageSize.textHandle, 'String', str);

vw = setSlider(vw, vw.ui.montageSize, 1, 0);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [K, sk] = wienerKernCompute(kern, x, x2)

% WIENERKERNCOMPUTE Compute the WIENER kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the wiener
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the wiener
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : wienerKernParamInit, kernCompute, kernCreate, wienerKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2009

% KERN

  if any(x<0)
    error('WIENER kernel only valid for time greater than zero')
  end
  if nargin < 3
    K = repmat(x, 1, size(x, 1));
    sk = min(K, K');
  else
    if any(x2<0)
      error('WIENER kernel only valid for time greater than zero')
    end
    sk = min(repmat(x, 1, size(x2, 1)), repmat(x2', size(x, 1), 1));
  end
  K = kern.variance*sk;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% test for CS minoration

path(path, 'toolbox/');

% number of measurements
n = 280;
% number of atoms
m = 300;
% sparsity
s = 5;

% random measurements matrix
D = compute_compressed_sensing_matrix(n,m);

Delta = compute_cs_bounds(D,s)

author: gpeyre

repo: matlab-toolboxes


text:
clear all; close all; clc;
%% Basic Generative Adversarial Network
%% Load Data
load('mnistAll.mat')
trainX = preprocess(mnist.train_images);
trainY = mnist.train_labels;
testX = preprocess(mnist.test_images);
testY = mnist.test_labels;
%% Settings
settings.latent_dim = 100;
settings.batch_size = 32; settings.image_size = [28,28,1];
settings.lrD = 0.0002; settings.lrG = 0.0002; settings.beta1 = 0.5;
settings.beta2 = 0.999; settings.maxepochs = 50;

%% Initialization
%% Generator
paramsGen.FCW1 = dlarray(...
    initializeGaussian([256,settings.latent_dim],.02));
paramsGen.FCb1 = dlarray(zeros(256,1,'single'));
paramsGen.BNo1 = dlarray(zeros(256,1,'single'));
paramsGen.BNs1 = dlarray(ones(256,1,'single'));
paramsGen.FCW2 = dlarray(initializeGaussian([512,256]));
paramsGen.FCb2 = dlarray(zeros(512,1,'single'));
paramsGen.BNo2 = dlarray(zeros(512,1,'single'));
paramsGen.BNs2 = dlarray(ones(512,1,'single'));
paramsGen.FCW3 = dlarray(initializeGaussian([1024,512]));
paramsGen.FCb3 = dlarray(zeros(1024,1,'single'));
paramsGen.BNo3 = dlarray(zeros(1024,1,'single'));
paramsGen.BNs3 = dlarray(ones(1024,1,'single'));
paramsGen.FCW4 = dlarray(initializeGaussian(...
    [prod(settings.image_size),1024]));
paramsGen.FCb4 = dlarray(zeros(prod(settings.image_size)...
    ,1,'single'));

stGen.BN1 = []; stGen.BN2 = []; stGen.BN3 = [];

%% Discriminator
paramsDis.FCW1 = dlarray(initializeGaussian([1024,...
     prod(settings.image_size)],.02));
paramsDis.FCb1 = dlarray(zeros(1024,1,'single'));
paramsDis.BNo1 = dlarray(zeros(1024,1,'single'));
paramsDis.BNs1 = dlarray(ones(1024,1,'single'));
paramsDis.FCW2 = dlarray(initializeGaussian([512,1024]));
paramsDis.FCb2 = dlarray(zeros(512,1,'single'));
paramsDis.BNo2 = dlarray(zeros(512,1,'single'));
paramsDis.BNs2 = dlarray(ones(512,1,'single'));
paramsDis.FCW3 = dlarray(initializeGaussian([256,512]));
paramsDis.FCb3 = dlarray(zeros(256,1,'single'));
paramsDis.FCW4 = dlarray(initializeGaussian([1,256]));
paramsDis.FCb4 = dlarray(zeros(1,1,'single'));

stDis.BN1 = []; stDis.BN2 = [];

% average Gradient and average Gradient squared holders
avgG.Dis = []; avgGS.Dis = []; avgG.Gen = []; avgGS.Gen = [];
%% Train
numIterations = floor(size(trainX,2)/settings.batch_size);
out = false; epoch = 0; global_iter = 0;
while ~out
    tic;
    trainXshuffle = trainX(:,randperm(size(trainX,2)));
    fprintf('Epoch %d\n',epoch)
    for i=1:numIterations
        global_iter = global_iter+1;
        noise = gpdl(randn([settings.latent_dim,...
            settings.batch_size]),'CB');
        idx = (i-1)*settings.batch_size+1:i*settings.batch_size;
        XBatch=gpdl(single(trainXshuffle(:,idx)),'CB');

        [GradGen,GradDis,stGen,stDis] = ...
                dlfeval(@modelGradients,XBatch,noise,...
                paramsGen,paramsDis,stGen,stDis);

        % Update Discriminator network parameters
        [paramsDis,avgG.Dis,avgGS.Dis] = ...
            adamupdate(paramsDis, GradDis, ...
            avgG.Dis, avgGS.Dis, global_iter, ...
            settings.lrD, settings.beta1, settings.beta2);

        % Update Generator network parameters
        [paramsGen,avgG.Gen,avgGS.Gen] = ...
            adamupdate(paramsGen, GradGen, ...
            avgG.Gen, avgGS.Gen, global_iter, ...
            settings.lrG, settings.beta1, settings.beta2);

        if i==1 || rem(i,20)==0
            progressplot(paramsGen,stGen,settings);
%             if i==1 || (epoch>=0 && i==1)
%                 h = gcf;
%                 % Capture the plot as an image
%                 frame = getframe(h);
%                 im = frame2im(frame);
%                 [imind,cm] = rgb2ind(im,256);
%                 % Write to the GIF File
%                 if epoch == 0
%                   imwrite(imind,cm,'GANmnist.gif','gif', 'Loopcount',inf);
%                 else
%                   imwrite(imind,cm,'GANmnist.gif','gif','WriteMode','append');
%                 end
%             end
        end

    end

    elapsedTime = toc;
    disp("Epoch "+epoch+". Time taken for epoch = "+elapsedTime + "s")
    epoch = epoch+1;
    if epoch == settings.maxepochs
        out = true;
    end
end
%% Helper Functions
%% preprocess
function x = preprocess(x)
x = double(x)/255;
x = (x-.5)/.5;
x = reshape(x,28*28,[]);
end
%% extract data
function x = gatext(x)
x = gather(extractdata(x));
end
%% gpu dl array wrapper
function dlx = gpdl(x,labels)
dlx = gpuArray(dlarray(x,labels));
end
%% Weight initialization
function parameter = initializeGaussian(parameterSize,sigma)
if nargin < 2
    sigma = 0.05;
end
parameter = randn(parameterSize, 'single') .* sigma;
end
%% Generator
function [dly,st] = Generator(dlx,params,st)
% fully connected
%1
dly = fullyconnect(dlx,params.FCW1,params.FCb1);
dly = leakyrelu(dly,0.2);
% if isempty(st.BN1)
%     [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,params.BNs1);
% else
%     [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,...
%         params.BNs1,st.BN1.mu,st.BN1.sig);
% end
%2
dly = fullyconnect(dly,params.FCW2,params.FCb2);
dly = leakyrelu(dly,0.2);
% if isempty(st.BN2)
%     [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,params.BNs2);
% else
%     [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,...
%         params.BNs2,st.BN2.mu,st.BN2.sig);
% end
%3
dly = fullyconnect(dly,params.FCW3,params.FCb3);
dly = leakyrelu(dly,0.2);
% if isempty(st.BN3)
%     [dly,st.BN3.mu,st.BN3.sig] = batchnorm(dly,params.BNo3,params.BNs3);
% else
%     [dly,st.BN3.mu,st.BN3.sig] = batchnorm(dly,params.BNo3,...
%         params.BNs3,st.BN3.mu,st.BN3.sig);
% end
%4
dly = fullyconnect(dly,params.FCW4,params.FCb4);
% tanh
dly = tanh(dly);
end
%% Discriminator
function [dly,st] = Discriminator(dlx,params,st)
% fully connected
%1
dly = fullyconnect(dlx,params.FCW1,params.FCb1);
dly = leakyrelu(dly,0.2);
dly = dropout(dly);
% if isempty(st.BN1)
%     [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,params.BNs1);
% else
%     [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,...
%         params.BNs1,st.BN1.mu,st.BN1.sig);
% end
%2
dly = fullyconnect(dly,params.FCW2,params.FCb2);
dly = leakyrelu(dly,0.2);
dly = dropout(dly);
% if isempty(st.BN2)
%     [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,params.BNs2);
% else
%     [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,...
%         params.BNs2,st.BN2.mu,st.BN2.sig);
% end
%3
dly = fullyconnect(dly,params.FCW3,params.FCb3);
dly = leakyrelu(dly,0.2);
dly = dropout(dly);
%4
dly = fullyconnect(dly,params.FCW4,params.FCb4);
% sigmoid
dly = sigmoid(dly);
end
%% modelGradients
function [GradGen,GradDis,stGen,stDis]=modelGradients(x,z,paramsGen,...
    paramsDis,stGen,stDis)
[fake_images,stGen] = Generator(z,paramsGen,stGen);
d_output_real = Discriminator(x,paramsDis,stDis);
[d_output_fake,stDis] = Discriminator(fake_images,paramsDis,stDis);

% Loss due to true or not
d_loss = -mean(.9*log(d_output_real+eps)+log(1-d_output_fake+eps));
g_loss = -mean(log(d_output_fake+eps));

% For each network, calculate the gradients with respect to the loss.
GradGen = dlgradient(g_loss,paramsGen,'RetainData',true);
GradDis = dlgradient(d_loss,paramsDis);
end
%% progressplot
function progressplot(paramsGen,stGen,settings)
r = 5; c = 5;
noise = gpdl(randn([settings.latent_dim,r*c]),'CB');
gen_imgs = Generator(noise,paramsGen,stGen);
gen_imgs = reshape(gen_imgs,28,28,[]);

fig = gcf;
if ~isempty(fig.Children)
    delete(fig.Children)
end

I = imtile(gatext(gen_imgs));
I = rescale(I);
imagesc(I)
title("Generated Images")
colormap gray

drawnow;
end
%% dropout
function dly = dropout(dlx,p)
if nargin < 2
    p = .3;
end
n = p*10;
mask = randi([1,10],size(dlx));
mask(mask<=n)=0;
mask(mask>n)=1;
dly = dlx.*mask;

end

author: zcemycl

repo: Matlab-GAN

save_path: github-repos/MATLAB/zcemycl-Matlab-GAN


text:
function [accgvs, accpps] = iccvEvaluateGeometry(imsegs)

resdir = './data4/smallsegs/3/';

totalgvs = 0;
accgvs = 0;

totalpps = 0;
accpps = 0;

for f = 1:numel(imsegs)
    fn = imsegs(f).imname;
    load([resdir strtok(fn, '.') '_seg.mat']);

    gvsgt = imsegs(f).vert_labels(imsegs(f).segimage);
    ppsgt = imsegs(f).horz_labels(imsegs(f).segimage);
    ppsgt = (ppsgt<=3) + (ppsgt==4)*2 + (ppsgt==5)*3;

    pg = bndinfo.result.geomProb;
    gvspg = [pg(:, 1) sum(pg(:, 2:4),2) pg(:, 5)];
    ppspg = pg(:, 2:4);
    [tmp, gvsres] = max(gvspg, [], 2);
    [tmp, ppsres] = max(ppspg, [], 2);

    gvsres = gvsres(bndinfo.wseg);
    ppsres = ppsres(bndinfo.wseg);

    accgvs = accgvs + sum(gvsres(:)==gvsgt(:)) / sum(gvsgt(:)>0);
    totalgvs = totalgvs + mean(gvsgt(:)>0);

    accpps = accpps + sum(ppsres(:)==ppsgt(:)) / sum(ppsgt(:)>0);
    totalpps = totalpps + mean(ppsgt(:)>0);

    disp(num2str([accgvs/totalgvs accpps/totalpps]));

end

author: Cloud-CV

repo: object-proposals


text:
function y = cvx_isvalid( x, full )
persistent remap
if isempty( remap ),
    remap = cvx_remap( 'valid' );
end
y = remap( cvx_classify( x ) );
if nargin < 2 || ~full,
	y = all( y );
else
	y = reshape( y, x.size_ );
end

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

author: yu-jiang

repo: radpbook

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


text:
%% (Internal) Display a message framed to a string
%
%     This function format a string for fprintf-like functions in order to
%     present a message string2display entitled with '-' characters below.
%
%   str_out = disp_string_title(style, string2display )
%
% Arguments:
%
%      + style: if numeric, could be the stream identifier for fprintf
%      function. If char, can be the "style" argument for cprintf function.
%
%      + string2display: The string to frame.
%
% Output:
%
%      + str_out: resulting string
%
% Example:
%
%
% See also disp_string_framed, disp_option_enumeration
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function str_out = disp_string_title(style, string2display )

    str_out = [];

    if( nargin == 1 || isempty(string2display) )
        string2display = style;
        style = 'text';
        fid = 1;
    end

    if( isempty(style) )
        style = 'text';
        fid = 1;
    elseif( isnumeric(style) )
        fid = style;
        style = 'text';
    else
        if( ischar(style) )
            fid = 1;
        else
            style = 'text';
            fid = 1;
        end
    end

    str_aux = sprintf(' %s ', string2display);
    str_aux1 = repmat('-', 1, length(str_aux));

    if( fid == 0)
        string2display = strrep(string2display, '\', '\\');
        str_out = sprintf('\n%s\n# %s #\n%s\n\n', str_aux1, string2display, str_aux1);
    elseif( fid == 1)
        fprintf(fid, '\n');
        cprintf(style, '%s', str_aux);
        fprintf(fid, '\n%s\n\n', str_aux1);
    else
        fprintf(fid, '\n%s\n%s\n\n', str_aux, str_aux1);
    end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [maxSum,startIdx,endIdx]=maxSubsequenceSum(a)
%%MAXSUBSEQUENCESUM  Given a vector of numbers that might be positive or
%                    find the maximum value that can be obtained by summing
%                    a subsequence of the set a.
%
%INPUTS: a  A 1XN or NX1 vector of real numbers.
%
%OUTPUTS: maxSum The maximum value of sums possible by summing consecutive
%                values in a. If a only contains negative values, then this
%                will be zero as the minimum value is obtained by summing
%                nothing in a.
%startIdx,endIdx The beginning and ending indices of the maximum
%                subsequence of a. If a is all negative, these will both be
%                zero.
%
%The solution to the maximum subsequence sum is described in Chapter 2.4
%of [1]. Multiple approaches to the problem are described in [2]. The
%algorithm has O(N) complexity as noted in its presentation in [3].
%
%The algorithm has been slightly modified to handle arrays that contain all
%negative or zero elements.
%
%REFERENCES:
%[1] M.A.Weiss, Data Structures and Algorithm Analysis in C++, 2nd ed.
%    Reading, MA: Addison-Wesley, 1999.
%[2] J. Bentley, "Algorithm design techniques," Communications of the ACM,
%    vol. 27, no. 9, pp. 865-871, Sep. 1984.
%[3] D. Gries, "A note on the standard strategy for developing loop
%    invariants and loops," Department of Computer Science, Cornell
%    University, Ithaca, NY, Tech. Rep. TR 82-531, Oct. 1982.
%
%July 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numElA=length(a);

maxSum=0;
thisSum=0;
thisStartIdx=1;

startIdx=0;
endIdx=0;

maxVal=-Inf;
maxIdx=0;

for j=1:numElA
    %Finding the maximum value is needed so that something can be returned if
    %an array of only zero or negative numbers is passed.
    if(a(j)>maxVal)
        maxVal=a(j);
        maxIdx=j;
    end

    thisSum=thisSum+a(j);

    if(thisSum>maxSum)
        maxSum=thisSum;
        startIdx=thisStartIdx;
        endIdx=j;
    elseif(thisSum<0)
        thisSum=0;
        thisStartIdx=j+1;
    end
end

if(maxVal<=0)
    %If all of the entries are negative or zero, then return the index of
    %the maximum element.
    maxSum=maxVal;
    startIdx=maxIdx;
    endIdx=maxIdx;
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function addScaleBar(in,barLength,barWidth, barPosition,barColor)
% Adds scale bar into an image.
% out=addScaleBar(in,barLength,barWidth,barPosition,barColor)
% in: image to which the scale bar should be added.
% barLength: length of the scale bar in pixels (defaults barLength=10);
% barWidth: thickness of the scale bar in pixesl (default: barWidth = 3);
% barPosition: position of the scale bar in the image (default: barPosition = size(in)*9/10;
% barColor: color of the scale bar (default: 'white')

if and(exist('in','var'),~isempty(in));
    dipshow(in);
else
    % added to a current picture
    c=get(gca);
    in = zeros(floor([c.XLim(2),c.YLim(2)])); % just to get size right for barPosition
end;
if ~exist('barLength','var'); barLength = 10; end;
if ~exist('barWidth','var'); barWidth = 3; end;
if ~exist('barPosition','var'); barPosition = size(in)*9/10; end;
if ~exist('barColor','var'); barColor = [1 1 .99]; end;

line([barPosition(1) barPosition(1)-barLength], [barPosition(2) barPosition(2)],'color',barColor,'linewidth',barWidth)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function y=ill(t,x)
a=1;b=0.3;
y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]';

author: sfvsfv

repo: Mathematical-modeling


text:
function dz = doublePendulumDynamics(~,z,u,P)
%DZ = DOUBLEPENDULUMDYNAMICS(T,Z,U,P)
%
%FUNCTION:  This function computes the dynamics of a double
%    pendulum, and is designed to be called from ode45. The
%    model allows for arbitrary mass and inertia for each
%    link, but no friction or actuation
%
%INPUTS:
%    t = time. Dummy input for ode45. Not used.
%    z = [4xn] matrix of states.
%    u = [2xn] matrix of inputs
%    P = struct of parameters
%OUTPUTS:
%    dz = [4xn] matrix of state derivatives
%
%NOTES:
%    This file was automatically generated by writeDoublePendulumDynamics

m1 = P.m1; %link one mass
m2 = P.m2; %link two mass
g  = P.g ; %gravity
l1 = P.l1; %link one length
l2 = P.l2; %link two length
I1 = P.I1; %link one moment of inertia about its center of mass
I2 = P.I2; %link two moment of inertia about its center of mass
d1 = P.d1; %distance between link one center of mass and parent joint
d2 = P.d2; %distance between link two center of mass and parent joint

th1 = z(1,:); %link one absolute angle
dth1 = z(2,:); %link one angular rate
th2 = z(3,:); %link two absolute angle
dth2 = z(4,:); %link two angular rate

u1 = u(1,:); %torque acting on link 1 wrt ground
u2 = u(2,:); %torque acting on link 2 wrt ground

f1 = u1 - m2.*(d2.*dth2.^2.*cos(th2) + dth1.^2.*l1.*cos(th1)).*(d2.*sin(th2) + l1.*sin(th1)) + m2.*(d2.*dth2.^2.*sin(th2) + dth1.^2.*l1.*sin(th1)).*(d2.*cos(th2) + l1.*cos(th1)) - g.*m2.*(d2.*cos(th2) + l1.*cos(th1)) - d1.*g.*m1.*cos(th1);
f2 = u2 - d2.*g.*m2.*cos(th2) + d2.*dth1.^2.*l1.*m2.*sin(th1 - th2);

M11 = - I1 - d1.^2.*m1.*cos(th1).^2 - d1.^2.*m1.*sin(th1).^2 - l1.*m2.*cos(th1).*(d2.*cos(th2) + l1.*cos(th1)) - l1.*m2.*sin(th1).*(d2.*sin(th2) + l1.*sin(th1));
M12 = - I2 - d2.*m2.*cos(th2).*(d2.*cos(th2) + l1.*cos(th1)) - d2.*m2.*sin(th2).*(d2.*sin(th2) + l1.*sin(th1));
M21 = -d2.*l1.*m2.*cos(th1 - th2);
M22 = - I2 - d2.^2.*m2;

D = M11.*M22 - M12.*M21;

ddth1 = (f2.*M12 - f1.*M22)./D;
ddth2 = -(f2.*M11 - f1.*M21)./D;

dz = [...
    dth1; %derivative of link one absolute angle
    ddth1; %derivative of link one angular rate
    dth2; %derivative of link two absolute angle
    ddth2; %derivative of link two angular rate
];
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [L_BG]=uiThreshErode(varargin)

% function [L_BG]=uiThreshErode(M,thresholdInitial,blurKernelSize,groupCropOption,inConvexHull)

%% Parse input

thresholdInitial=[];
blurKernelSize=[];
groupCropOpt=[];
inConvexHullOpt=[];
switch nargin
    case 1
        M=varargin{1};
    case 2
        M=varargin{1};
        thresholdInitial=varargin{2};
    case 3
        M=varargin{1};
        thresholdInitial=varargin{2};
        blurKernelSize=varargin{3};
    case 4
        M=varargin{1};
        thresholdInitial=varargin{2};
        blurKernelSize=varargin{3};
        groupCropOpt=varargin{4};
    case 5
        M=varargin{1};
        thresholdInitial=varargin{2};
        blurKernelSize=varargin{3};
        groupCropOpt=varargin{4};
        inConvexHullOpt=varargin{5};
end

%Set defaults
if isempty(thresholdInitial)
     thresholdInitial=0.1;
end

%Quick fix for low thresholds
if thresholdInitial<0.01
    thresholdInitial=0.01;
end

if isempty(blurKernelSize)
     blurKernelSize=ndims(M);
end

if isempty(groupCropOpt)
     groupCropOpt=0;
end

if isempty(inConvexHullOpt)
     inConvexHullOpt=0;
end
%% PLOT SETTINGS

fontSize=10;
cMap=gray(250);
falpha=1;

logicVoxels=false(size(M));
logicVoxels(round(size(M,1)/2),:,:)=1;
logicVoxels(:,round(size(M,2)/2),:)=1;
logicVoxels(:,:,round(size(M,3)/2))=1;

%%
M_original=M;
M=double(M);

%% Normalizing the image data
M=dataNorm(M,1);

%% APPLY PRE-BLURRING
if blurKernelSize>0
    %Defining erosion/dilation kernel
    k=blurKernelSize;
    p=k-round(k./2);

    hb=gauss_kernel(k,ndims(M),1.5,'width'); %Gaussian kernel

    M_rep=zeros(size(M)+(2.*p));
    M_rep(p+(1:size(M,1)),p+(1:size(M,2)),p+(1:size(M,3)))=M;
    M = convn(M_rep,hb,'valid');
end

%% Defining erosion/dilation kernel
k=3;
hg=ones(k,k,k);
p=k-round(k./2);

%% Start loop
done=0;
while done==0
    %Get initial threshold level
    T_threshold=thresholdInitial; %default start value
    done=0;
    qc=1;
    runMode=1;
    while done==0

        switch runMode
            case 1 %Determining threshold
                L_BG=M>T_threshold;
            case 2 %Setting erosion/dilation
                %"Blurring current mask"
                L_BG_rep=zeros(size(M)+(2.*p));
                L_BG_rep(p+(1:size(M,1)),p+(1:size(M,2)),p+(1:size(M,3)))=L_BG;
                L_BG_blur = convn(double(L_BG_rep),hg,'valid');
        end

        if qc==1 %Plot figure on first iteration

            %Cropped image
            [Fs,Vs,Cs]=ind2patch(logicVoxels&L_BG,M_original,'v');

            hf1=cFigure;
            title(['Threshold is ',num2str(T_threshold),'*mean, press up to increase or down to decrease (by 10%), press space to keep and continue'],'FontSize',fontSize);
            hold on; xlabel('X-J','FontSize',fontSize);ylabel('Y-I','FontSize',fontSize);zlabel('Z-K','FontSize',fontSize);
            hs=patch('Faces',Fs,'Vertices',Vs,'EdgeColor','none', 'CData',Cs,'FaceColor','flat');
            set(hs,'FaceAlpha',falpha); hold on; grid on;  view(3); axis equal; axis tight; axis vis3d;  colormap(cMap); colorbar;
            set(gca,'FontSize',fontSize);
            drawnow;

        else %redefine patch data
            delete(hs); %remove patch data from figure

            %redefine patch data for the cropped image
            [Fs,Vs,Cs]=ind2patch(logicVoxels&L_BG,M_original,'v');

            %patch again
            if runMode==1
                title(['Threshold is ',num2str(T_threshold),'*mean, press up to increase or down to decrease (by 10%), press space to keep and continue'],'FontSize',fontSize);
            else
                title('Current cropping press up/down to dilate/erode, WARNING EROSION MAY REMOVE ENTIRE SLICES','FontSize',fontSize);
            end

            hs=patch('Faces',Fs,'Vertices',Vs,'EdgeColor','none', 'CData',Cs,'FaceColor','flat');
            set(hs,'FaceAlpha',falpha); hold on; grid on;  view(3); axis equal; axis tight; axis vis3d;  colormap(cMap); colorbar;
            set(gca,'FontSize',fontSize);
            drawnow;
        end

        [~,~,b]=qginput(1);
        switch b
            case 30 %  up arrow key => increase threshold by 10% OR dilate
                if runMode==1
                    T_threshold=T_threshold.*1.1;
                else
                    L_BG=(L_BG_blur>0); %blurred boundary elements that increased due to blur are now set to 1 as well so added to mask
                end
            case 31 %down arrow key => decrease threshold by 10% OR erode
                if runMode==1
                    T_threshold=T_threshold.*0.9;
                else
                    L_BG=(L_BG_blur==sum(hg(:))); %blurred boundary elements that decreased due to blur are now set to 0 as well so removed from mask
                end
            case 32 %spacebar key => confirms current threshold result, exits while loop
                if runMode==1
                    runMode=2;
                else
                    done=1;
                end
        end
        qc=qc+1;
    end

    %% GROUPING BASED CROPPING
    if groupCropOpt==1

        GROUP_STRUCT = bwconncomp(L_BG,6);
        IMAGE_OBJECTS=GROUP_STRUCT.PixelIdxList;

        %Filtering out groups that are too small
        group_sizes = cell2mat(cellfun(@(x) max(size(x)), IMAGE_OBJECTS,'UniformOutput',0)'); %get group sizes
        [~,ind_max]=max(group_sizes); %Index of largest group
        IND_OBJECT=IMAGE_OBJECTS{ind_max}; %indices of voxels for largest group

        %Redefining background logic using group
        L_BG=false(size(L_BG));
        L_BG(IND_OBJECT)=1;

        %Plotting the result
        delete(hs); %remove patch data from figure

        %redefine patch data for the cropped image
        [Fs,Vs,Cs]=ind2patch(logicVoxels&L_BG,M_original,'v');

        title('Grouping result','FontSize',fontSize);

        hs=patch('Faces',Fs,'Vertices',Vs,'EdgeColor','none', 'CData',Cs,'FaceColor','flat');
        set(hs,'FaceAlpha',falpha); hold on; grid on;  view(3); axis equal; axis tight; axis vis3d;  colormap(cMap);
        set(gca,'FontSize',fontSize);
        drawnow;
    end

    %% Convex hull based filling
    if inConvexHullOpt==1
        [I,J]=ndgrid(1:1:size(M,1),1:1:size(M,2));
        for q=1:1:size(M,3);
            l_bg=L_BG(:,:,q);
            if nnz(l_bg)>0
                if nnz(l_bg)>3
                    try
                        P=[I(l_bg) J(l_bg)];
                        K = convhull(P(:,1),P(:,2));
                        indNot=find(~l_bg);
                        [Lp,~] = inpolygon(I(indNot),J(indNot),P(K,1),P(K,2));
                        indAdd=indNot(Lp);
                        L=false(size(M,1),size(M,2));
                        L(indAdd)=1;
                        L_BG(:,:,q)=L_BG(:,:,q) | L;
                    end
                end
            else
               L_BG(:,:,q)=0;
            end
        end

        title('Convexhull fill result','FontSize',fontSize);
        delete(hs); %remove patch data from figure
        %redefine patch data for the cropped image
        [Fs,Vs,Cs]=ind2patch(logicVoxels&L_BG,M_original,'v');

        hs=patch('Faces',Fs,'Vertices',Vs,'EdgeColor','none', 'CData',Cs,'FaceColor','flat');
        set(hs,'FaceAlpha',falpha); hold on; grid on;  view(3); axis equal; axis tight; axis vis3d;  colormap(cMap);
        set(gca,'FontSize',fontSize);
        drawnow;
    end

    %%
    Q = questdlg('Would you like to keep this result?','Keep result?','YES','NO','YES');
    switch Q
        case 'YES'
            done=1;
        case 'NO'
            done=0;
    end
    close(hf1)

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function area = triangle_area_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_AREA_2D computes the area of a triangle in 2D.
%
%  Discussion:
%
%    If the triangle's vertices are given in counterclockwise order,
%    the area will be positive.  If the triangle's vertices are given
%    in clockwise order, the area will be negative!
%
%    An earlier version of this routine always returned the absolute
%    value of the computed area.  I am convinced now that that is
%    a less useful result!  For instance, by returning the signed
%    area of a triangle, it is possible to easily compute the area
%    of a nonconvex polygon as the sum of the (possibly negative)
%    areas of triangles formed by node 1 and successive pairs of vertices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real AREA, the area of the triangle.
%
  area = 0.5 * ( ...
      t(1,1) * ( t(2,2) - t(2,3) ) ...
    + t(1,2) * ( t(2,3) - t(2,1) ) ...
    + t(1,3) * ( t(2,1) - t(2,2) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% function to get rcnn detected part boxes
%% Written by Ning Zhang

function detect_boxes = get_rcnn_detections(config)
try
  load('caches/rcnn_detect_boxes.mat');
catch
  try
    load('caches/rcnn_part_models.mat');
  catch
    disp('Train part detectors in rcnn and save it in the caches');
    exit(-1);
  end
  detect_boxes = test_rcnn_parts(part_models, config);
  save('caches/rcnn_detect_boxes.mat', 'detect_boxes');
end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function model = gpReconstruct(kern, noise, gpInfo, X, y)

% GPRECONSTRUCT Reconstruct an GP form component parts.
% FORMAT
% DESC takes component parts of an GP model and reconstructs the
% GP model. The component parts are normally retrieved from a
% saved file.
% ARG kern : a kernel structure for the GP.
% ARG noise : a noise structure for the GP (currently ignored).
% ARG gpInfo : the active set and other information stored in a structure.
% ARG X : the input training data for the GP.
% ARG y : the output target training data for the GP.
% RETURN model : an GP model structure that combines the component
% parts.
%
% SEEALSO : gpDeconstruct, gpCreate
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2009

% GP

  model = gpInfo;
  model.X = X;
  model.y = y;
  model.kern = kern;
  if ~isempty(noise)
    model.noise = noise;
  end
  model.m = gpComputeM(model);

  if isfield(model, 'computeS') && model.computeS
    model.S = model.m*model.m';
  end
  params = gpExtractParam(model);
  model = gpExpandParam(model, params);

end


text:
function [clrmap,clrIX_x,clrIX_y] = MapXYto2Dcolormap(gIX_in,X,Y,Xrange,Yrange,cmap2D)
% gIX_in corresponds to the 2D values X and Y.
% The output clrmap is a linear colormap of the same length as numK
% (=max(gIX_in)), and same as clrIX_x and clrIX_y.
% clrIX_x and clrIX_y are the coordinates in the 2D colormap for the range of values in gIX_in.

if nargin < 4
    Xrange = [min(X),max(X)];
end
if nargin < 5
    Yrange = [min(Y),max(Y)];
end

if nargin < 6
    cmap2D = MakeDiagonal2Dcolormap(huex,satmin,pw);
end
res = size(cmap2D,1);

% check input dimensions
if size(X,1)~=size(Y,1)
    if size(X,1)>1
        Y = Y';
    else
        X = X';
    end
end
assert(isequal(size(X),size(Y)));

% set data range
X(X<Xrange(1)) = Xrange(1);
Y(Y<Yrange(1)) = Yrange(1);
X(X>Xrange(2)) = Xrange(2);
Y(Y>Yrange(2)) = Yrange(2);

clrIX_x = round((X-Xrange(1))/(Xrange(2)-Xrange(1))*(res-1))+1;
clrIX_y = round((Y-Yrange(1))/(Yrange(2)-Yrange(1))*(res-1))+1;

%%
clrmap = zeros(length(clrIX_x),3);
clrmap_2D_flat = reshape(cmap2D,res*res,3); % for efficient indexing

U = unique(gIX_in);

ix = (clrIX_x(U)-1)*res+clrIX_y(U);
clrmap(U,:) = clrmap_2D_flat(ix,:);
% for i = 1:length(U)
% %     ix = sub2ind([res,res],clrIX_y(U(i))',clrIX_x(U(i))');
%     ix = (clrIX_x(U(i))'-1)*res+clrIX_y(U(i))';
% %     center_ix = (M_xyz(j,2)-1)*dimv(1)+M_xyz(j,1); % linear pixel index, faster equivalent of:
%         %     center_ix = sub2ind([dimv(1),dimv(2)],M_xyz(j,1),M_xyz(j,2));
%
%
%     clrmap(U(i),:) = clrmap_2D_flat(ix,:);
% end
end

author: xiuyechen

repo: FishExplorer


text:
function dt = datatypes
% Dictionary of datatypes
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: datatypes.m 7147 2017-08-03 14:07:01Z spm $

persistent dtype
if isempty(dtype)
    t = true;
    f = false;
    table = {...
        0   ,'UNKNOWN'   ,'uint8'   ,@uint8  ,1,1  ,t,t,f
        1   ,'BINARY'    ,'uint1'   ,@logical,1,1/8,t,t,f
        256 ,'INT8'      ,'int8'    ,@int8   ,1,1  ,t,f,t
        2   ,'UINT8'     ,'uint8'   ,@uint8  ,1,1  ,t,t,t
        4   ,'INT16'     ,'int16'   ,@int16  ,1,2  ,t,f,t
        512 ,'UINT16'    ,'uint16'  ,@uint16 ,1,2  ,t,t,t
        8   ,'INT32'     ,'int32'   ,@int32  ,1,4  ,t,f,t
        768 ,'UINT32'    ,'uint32'  ,@uint32 ,1,4  ,t,t,t
        1024,'INT64'     ,'int64'   ,@int64  ,1,8  ,t,f,f
        1280,'UINT64'    ,'uint64'  ,@uint64 ,1,8  ,t,t,f
        16  ,'FLOAT32'   ,'float32' ,@single ,1,4  ,f,f,t
        64  ,'FLOAT64'   ,'double'  ,@double ,1,8  ,f,f,t
        1536,'FLOAT128'  ,'float128',@error  ,1,16 ,f,f,f
        32  ,'COMPLEX64' ,'float32' ,@single ,2,4  ,f,f,f
        1792,'COMPLEX128','double'  ,@double ,2,8  ,f,f,f
        2048,'COMPLEX256','float128',@error  ,2,16 ,f,f,f
        128 ,'RGB24'     ,'uint8'   ,@uint8  ,3,1  ,t,t,f};
    dtype = struct(...
        'code'     ,table(:,1),...
        'label'    ,table(:,2),...
        'prec'     ,table(:,3),...
        'conv'     ,table(:,4),...
        'nelem'    ,table(:,5),...
        'size'     ,table(:,6),...
        'isint'    ,table(:,7),...
        'unsigned' ,table(:,8),...
        'min',-Inf,'max',Inf',...
        'supported',table(:,9));
    for i=1:length(dtype)
        if dtype(i).isint
            if dtype(i).unsigned
                dtype(i).min =  0;
                dtype(i).max =  2^(8*dtype(i).size)-1;
            else
                dtype(i).min = -2^(8*dtype(i).size-1);
                dtype(i).max =  2^(8*dtype(i).size-1)-1;
            end
        end
    end
end

dt = dtype;

author: spm


text:
function GCO_BuildLib(Options)
% GCO_BuildLib    Attempt to compile and link the GCO_MATLAB library.
%    GCO_BuildLib is used internally by all other GCO_MATLAB commands
%    to recompile the wrapper library if it is not yet built.
%
%    YOU DO NOT NEED TO EXPLICITLY CALL THIS FUNCTION, unless you want to
%    customise the build settings via GCO_BuildLib(Options).
%    Default options:
%      Options.Debug=0        % optimised, detailed checking disabled
%      Options.EnergyType32=0 % 64-bit energy counter, 32-bit energy terms
%
%    Example:
%      % Enable detailed assertions (e.g. than energy does not go up
%      % during expansion) and use 32-bit energy counters (slightly faster)
%      GCO_BuildLib(struct('Debug',1,'EnergyType32',1));
%

if (nargin < 1)
    Options = struct();
end
if (~isfield(Options,'Debug')), Options.Debug = 0; end
if (~isfield(Options,'EnergyType32')), Options.EnergyType32 = 0; end
if (~isfield(Options,'Force')), Options.Force = 1; end

MEXFLAGS = '';
if (strcmp(computer(),'GLNXA64') || strcmp(computer(),'PCWIN64') || strcmp(computer(),'MACI64'))
    MEXFLAGS = [MEXFLAGS ' -largeArrayDims -DA64BITS'];
end
if (Options.Debug)
    MEXFLAGS = [MEXFLAGS ' -g'];
end
if (Options.EnergyType32)
    MEXFLAGS = [MEXFLAGS ' -DGCO_ENERGYTYPE32'];
end
if (strcmp(computer(),'PCWIN')) % link with libut for user interruptibility
    MEXFLAGS = [MEXFLAGS ' -D_WIN32 "' matlabroot() '\extern\lib\win32\microsoft\libut.lib"' ];
elseif (strcmp(computer(),'PCWIN64'))
    MEXFLAGS = [MEXFLAGS ' -D_WIN64 "' matlabroot() '\extern\lib\win64\microsoft\libut.lib"' ];
else
    MEXFLAGS = [MEXFLAGS ' -lut' ];
end

LIB_NAME = 'gco_matlab';
GCOMATDIR = fileparts(mfilename('fullpath'));
GCODIR = fileparts(GCOMATDIR);
OUTDIR = [ GCOMATDIR filesep 'bin' ];
[status msg msgid] = mkdir(GCOMATDIR, 'bin'); % Create bin directory
addpath(OUTDIR);                              % and add it to search path
if (~Options.Force && exist('gco_matlab')==3)
    return;
end
clear gco_matlab;

mexcmd = ['mex ' MEXFLAGS ' -outdir ''' OUTDIR ''' -output ' LIB_NAME ' ' ];

% Append all source file names to the MEX command string
SRCCPP = {
    [GCOMATDIR filesep 'gco_matlab.cpp'],
    [GCODIR filesep 'GCoptimization.cpp'],
    [GCODIR filesep 'graph.cpp'],
    [GCODIR filesep 'maxflow.cpp'],
    [GCODIR filesep 'LinkedBlockList.cpp']
    };
for f=1:length(SRCCPP)
    mexcmd = [mexcmd ' ''' SRCCPP{f} ''' '];
end

eval(mexcmd);  % compile and link in one step

end

author: thfylsty


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef cableClass<handle
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The ``cableClass`` creates a ``cable`` object saved to the MATLAB
    % workspace. The ``cableClass`` includes properties and methods used
    % to define cable connections relative to other bodies.
    % It is suggested that the ``cableClass`` be used for connections between
    % joints or ptos.
    %
    %.. autoattribute:: objects.cableClass.cableClass
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (SetAccess = 'public', GetAccess = 'public') %input file
        damping                 = 0                                         % (`float`) Cable damping coefficient (N/(m/s)). Default = `0`.
        inertia                 = [1 1 1];                                  % (`1 x 3 float vector`) body inertia kg-m^2, default [1 1 1]
        initial                 = struct(...                                % (`structure`) Defines the initial displacement of the body.
            'displacement',         [0 0 0],...                             %
            'axis',                 [0 1 0], ...                            %
            'angle',                0)                                      % (`structure`) Defines the initial displacement of the body. ``displacement`` (`3x1 float vector`) is defined as the initial displacement of the body center of gravity (COG) [m] in the following format [x y z], Default = [``0 0 0``]. ``axis`` (`3x1 float vector`) is defined as the axis of rotation in the following format [x y z], Default = [``0 1 0``]. ``angle`` (`float`) is defined as the initial angular displacement of the body COG [rad], Default = ``0``.
        length                  = 0                                         % (`float`) Cable equilibrium length (m), calculated from rotloc and preTension. Default =`0`.
        linearDamping           = [0 0 0 0 0 0];                            % (`1 x 6 float vector`)linear damping aplied to body motions
        mass                    = 1;                                        % (`float`) mass in kg, default 1
        name                    = 'NOT DEFINED'                             % (`string`) Defines the Cable name. Default = ``NOT DEFINED``.
        orientation             = struct(...                                % (`structure`) Defines the orientation axis of the pto.
            'z',                    [0, 0, 1], ...                          %
            'y',                    [0, 1, 0], ...                          %
            'x',                    [], ...                                 %
            'rotationMatrix',       [])                                     % (`structure`) Defines the orientation axis of the pto. ``z`` (`3x1 float vector`) defines the direciton of the Z-coordinate of the pto, Default = [``0 0 1``]. ``y`` (`3x1 float vector`) defines the direciton of the Y-coordinate of the pto, Default = [``0 1 0``]. ``x`` (`3x1 float vector`) internally calculated vector defining the direction of the X-coordinate for the pto, Default = ``[]``. ``rotationMatrix`` (`3x3 float matrix`) internally calculated rotation matrix to go from standard coordinate orientation to the pto coordinate orientation, Default = ``[]``.
        paraview                = 1;                                        % (`integer`) Flag for visualisation in Paraview either 0 (no) or 1 (yes). Default = ``1`` since only called in paraview.
        preTension              = 0                                         % (`float`) Cable pretension (N).
        quadDrag                = struct(...                                % (`structure`) Defines the viscous quadratic drag forces.
            'area',                 [0 0 0 0 0 0], ...                      %
            'drag',                 zeros(6), ...                           %
            'cd',               [0 0 0 0 0 0]);                             % (`structure`) Defines the viscous quadratic drag forces. First option define ``drag``, (`6x6 float matrix`), Default = ``zeros(6)``. Second option define ``cd``, (`6x1 float vector`), Default = ``zeros(6,1)``, and ``area``, (`6x1 float vector`), Default = ``zeros(6,1)``.
        stiffness               = 0                                         % (`float`) Cable stiffness (N/m). Default = `0`.
        viz                     = struct(...                                % (`structure`) Defines visualization properties in either SimScape or Paraview.
            'color',                [1 0 0.5], ...                          %
            'opacity',              1)                                      % (`structure`) Defines visualization properties in either SimScape or Paraview. ``color`` (`3x1 float vector`) is defined as the body visualization color, Default = [``1 1 0``]. ``opacity`` (`integer`) is defined as the body opacity, Default = ``1``.
    end

    properties (SetAccess = 'public', GetAccess = 'public')%internal
        base                    = struct(...                                % (`structure`) Defines the base connection. `centerBuoyancy` (`1 x 3 float vector`) center of buoyancy location of the base drag body, Default = `[0 0 0]`. `centerGravity` (`1 x 3 float vector`) center of gravity location of the base drag body, Default = `[0 0 0]`. `location` (`3x1 float vector`) base location [m], Defined in the following format [x y z], Default = ``[]``. `name` (`string`) name of the base constraint or PTO, Default = `'NOT DEFINED'`;
            'centerBuoyancy',   	[0 0 0], ...                            %
            'centerGravity',        [0 0 0], ...                            %
            'location',             [], ...                                 %
            'name',                 'NOT DEFINED')                          % (`structure`) Defines the base connection. `centerBuoyancy` (`1 x 3 float vector`) center of buoyancy location of the base drag body, Default = `[0 0 0]`. `centerGravity` (`1 x 3 float vector`) center of gravity location of the base drag body, Default = `[0 0 0]`. `location` (`3x1 float vector`) base location [m], Defined in the following format [x y z], Default = ``[]``. `name` (`string`) name of the base constraint or PTO, Default = `'NOT DEFINED'`;
        follower                  = struct(...                              % (`structure`) Defines the follower connection. `centerBuoyancy` (`1 x 3 float vector`) center of buoyancy location of the base drag body, Default = `[0 0 0]`. `centerGravity` (`1 x 3 float vector`) center of gravity location of the base drag body, Default = `[0 0 0]`. `location` (`3x1 float vector`) base location [m], Defined in the following format [x y z], Default = ``[]``. `name` (`string`) name of the base constraint or PTO, Default = `'NOT DEFINED'`;
            'centerBuoyancy',   	[0 0 0], ...                            %
            'centerGravity',        [0 0 0], ...                            %
            'location',             [], ...                                 %
            'name',                 'NOT DEFINED')                          % (`structure`) Defines the follower connection. `centerBuoyancy` (`1 x 3 float vector`) center of buoyancy location of the base drag body, Default = `[0 0 0]`. `centerGravity` (`1 x 3 float vector`) center of gravity location of the base drag body, Default = `[0 0 0]`. `location` (`3x1 float vector`) base location [m], Defined in the following format [x y z], Default = ``[]``. `name` (`string`) name of the base constraint or PTO, Default = `'NOT DEFINED'`;
        location                = [999 999 999]                             % (`3x1 float vector`) pto location [m]. Defined in the following format [x y z]. Default = ``[999 999 999]``.
        volume                  = [];                                       % (`float`) displacement volume, defaults to neutral buoyancy
    end

    properties (SetAccess = 'private', GetAccess = 'public') %internal
        number                  = []                                        % Cable number.
    end

    %%
    methods

        function obj = cableClass(name,baseName,followerName)
            % This method initilizes the ``cableClass`` and creates a
            % ``cable`` object.
            %
            % Parameters
            % ------------
            %     name : string
            %         String specifying the name of the cable
            %
            %     baseConnection : string
            %         Variable for the base constraint/pto as a string
            %
            %     followerConnection : string
            %         Variable for the follower constraint/pto as a string
            %
            % Returns
            % ------------
            %     cable : obj
            %         cableClass object
            %
            if exist('name','var')
                obj.name = name;
                obj.base.name = baseName;
                obj.follower.name= followerName;
                obj.base.location = evalin('caller',[baseName '.location']);
                obj.follower.location = evalin('caller',[followerName '.location']);
            else
                error('The cable class number(s) in the wecSimInputFile must be specified in ascending order starting from 1. The cableClass() function should be called first to initialize each cable with a name.')
            end
        end

        function setTransPTOLoc(obj)
            % This method specifies the translational PTO location as half-
            % way between the fixed ends of the cable if not previously
            % set.
            if any(obj.location == [999 999 999])
                rotDiff = obj.base.location - obj.follower.location;
                obj.location = obj.follower.location + rotDiff/2;
                fprintf('\n\t location undefined, set halfway between follower.location and base.location \n')
            end
        end

        function obj = checkLoc(obj,action)
            % This method checks WEC-Sim user inputs and generate an error
            % message if the constraint location is not defined in constraintClass.
            % Checks if location is set and outputs a warning or error. Used in mask Initialization.
            switch action
                case 'W'
                    if obj.location == 999 % Because "Allow library block to modify its content" is selected in block's mask initialization, this command runs twice, but warnings cannot be displayed during the first initialization.
                        obj.location = [888 888 888];
                    elseif obj.location == 888
                        obj.location = [0 0 0];
                    end
                case 'E'
                    try
                        if obj.location == 999
                            s1 = ['For ' obj.name ': pto(#).location needs to be specified in the WEC-Sim input file.'...
                                ' pto.location is the [x y z] location, in meters, for the rotational PTO.'];
                            error(s1)
                        end
                    catch exception
                        throwAsCaller(exception)
                    end
            end
        end

        function obj = setOrientation(obj)
            % This method calculates the constraint ``x`` vector and
            % ``rotationMatrix`` matrix in the ``orientation`` structure
            % based on user input.
            obj.orientation.z = obj.orientation.z / norm(obj.orientation.z);
            obj.orientation.y = obj.orientation.y / norm(obj.orientation.y);
            z = obj.orientation.z;
            y = obj.orientation.y;
            if abs(dot(y,z))>0.001
                error('The Y and Z vectors defining the constraint''s orientation must be orthogonal.')
            end
            x = cross(y,z)/norm(cross(y,z));
            x = x(:)';
            obj.orientation.x = x;
            obj.orientation.rotationMatrix  = [x',y',z'];

        end

        function setInitDisp(obj, relCoord, axisAngleList, addLinDisp)
            % Function to set a body's initial displacement
            %
            % This function assumes that all rotations are about the same relative coordinate.
            % If not, the user should input a relative coordinate of 0,0,0 and
            % use the additional linear displacement parameter to set the cg or location
            % correctly.
            %
            % Parameters
            % ------------
            %    relCoord : [1 3] float vector
            %        Distance from x_rot to the body center of gravity or the constraint
            %        or pto location as defined by: relCoord = cg - x_rot. [m]
            %
            %    axisAngleList : [nAngle 4] float vector
            %        List of axes and angles of the rotations with the
            %        format: [n_x n_y n_z angle] (angle in rad)
            %        Rotations applied consecutively in order of dimension 1
            %
            %    addLinDisp : [1 3] float vector
            %        Initial linear displacement (in addition to the
            %        displacement caused by rotation) [m]
            %

            % initialize quantities before for loop
            axisList = axisAngleList(:,1:3);
            angleList = axisAngleList(:,4);
            nAngle = size(axisList,1);
            rotMat = eye(3);

            % Loop through all axes and angles.
            for i=1:nAngle
                rotMat = axisAngle2RotMat(axisList(i,:),angleList(i))*rotMat;
            end
            % calculate net axis-angle rotation
            [netAxis, netAngle] = rotMat2AxisAngle(rotMat);
            % calculate net displacement due to rotation
            rotatedRelCoord = relCoord*(rotMat');
            linDisp = rotatedRelCoord - relCoord;
            % apply rotation and displacement to object
            obj.initial.displacement = linDisp + addLinDisp;
            obj.initial.axis = netAxis;
            obj.initial.angle = netAngle;
        end

        function setVolume(obj, rho)
            % This method mades the mass of the cable drag bodies neutrally bouyant
            obj.volume = obj.mass/rho;
        end

        function dragForcePre(obj,rho)
            % This method performs Drag Body pre-processing calculations,
            % similar to hydroForcePre, but only loads in the necessary
            % values to calculate linear damping and viscous drag. Note
            % that body DOF is inherited from the length of the drag
            % coefficients.
            if  any(any(obj.quadDrag.drag)) ~= 1  %check if obj.quadDrag.drag is not defined
                obj.quadDrag.drag = diag(0.5*rho.*obj.quadDrag.cd.*obj.quadDrag.area);
            end

        end

        function linearDampingMatrix(obj)
            % This method makes the linear damping vector (as specified)
            % into a 6 x 6 matrix with this damping along the diagonal (as
            % required for calculation). Operates on the drag bodies
            % representing the cable dynamics.
            obj.linearDamping = diag(obj.linearDamping);
        end

        function setCb(obj)
            % This method sets the buoyancy center to equal the center of
            % gravity, if the center of buoyancy is not defined
            obj.base.centerBuoyancy = obj.base.centerGravity;
            obj.follower.centerBuoyancy = obj.follower.centerGravity;
        end

        function setCg(obj)
            % This method specifies the Cg of the drag bodies as coincident
            % with the fixed ends of the cable, if not otherwise specied.
            obj.base.centerGravity = obj.base.location;
            obj.follower.centerGravity = obj.follower.location;
        end

        function setLength(obj)
            % This method specifies length as the distance between cable fixed
            % ends (i.e. pretension = 0), if not otherwise specified.
            if ~any(obj.length) && ~any(obj.preTension)
                obj.length = sqrt((obj.base.location(1)-obj.follower.location(1)).^2 ...
                + (obj.base.location(2)-obj.follower.location(2)).^2 ...
                + (obj.base.location(3)-obj.follower.location(3)).^2);
                fprintf('\n\t cable(i).length undefined and cable(i).preTension undefined. \n \r',...
                    'cable(i).length set equal to distance between follower.location and base.location \n and cable(i).preTension set equal to zero \n');
            elseif ~any(obj.length) && any(obj.preTension)
                obj.length = sqrt((obj.base.location(1)-obj.follower.location(1)).^2 ...
                + (obj.base.location(2)-obj.follower.location(2)).^2 ...
                + (obj.base.location(3)-obj.follower.location(3)).^2) + obj.preTension/obj.stiffness;
            elseif ~any(obj.preTension) && any(obj.length)
                obj.preTension = obj.stiffness * (sqrt((obj.base.location(1)-obj.follower.location(1)).^2 ...
                + (obj.base.location(2)-obj.follower.location(2)).^2 ...
                + (obj.base.location(3)-obj.follower.location(3)).^2) - obj.length);
            elseif any(obj.preTension) && any(obj.length)
                error('System overdefined. Please define cable(i).preTension OR cable(i).length, not both.')
            end
        end

        function listInfo(obj)
            % This method prints cable information to the MATLAB Command Window.
            fprintf('\n\t***** Cable Name: %s *****\n',obj.name)
            fprintf('\tCable Stiffness           (N/m;Nm/rad) = %G\n',obj.stiffness)
            fprintf('\tCable Damping           (Ns/m;Nsm/rad) = %G\n',obj.damping)
        end

        function setNumber(obj,number)
            % Method to set the private number property
            obj.number = number;
        end
    end

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function geometry_test0367 ( )

%*****************************************************************************80
%
%% TEST0367 tests SEGMENT_POINT_DIST_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  test_num = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0367\n' );
  fprintf ( 1, '  SEGMENT_POINT_NEAR_2D computes the nearest point\n' );
  fprintf ( 1, '    from a line segment to a point in 2D.\n' );

  for test = 1 : test_num

    [ p1, seed ] = r8vec_uniform_01 ( dim_num, seed );
    [ p2, seed ] = r8vec_uniform_01 ( dim_num, seed );
    [ p, seed ] = r8vec_uniform_01 ( dim_num, seed );

    [ pn, dist, t ] = segment_point_near_2d ( p1, p2, p );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  TEST = %d', test );
    fprintf ( 1, '  P1 =   %12f  %12f\n', p1(1:dim_num) );
    fprintf ( 1, '  P2 =   %12f  %12f\n', p2(1:dim_num) );
    fprintf ( 1, '  P =    %12f  %12f\n', p(1:dim_num) );
    fprintf ( 1, '  PN =   %12f  %12f\n', pn(1:dim_num) );
    fprintf ( 1, '  DIST = %12f\n', dist );
    fprintf ( 1, '  T =    %12f\n', t );

  end

  return
end


text:
function chol_example
%CHOL_EXAMPLE simple Cholesky factorization example
% Example
%   chol_example
% See also: testall

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

D = 10 ;
X = 1 ;
o = 0 ;

A = sparse ([
D o X o o o o X o o
o D o o X o o o o X
X o D o o o X o o o
o o o D o o o o X X
o X o o D o o o X X
o o o o o D X X o o
o o X o o X D o o o
X o o o o X o D X X
o o o X X o o X D o
o X o X X o o X o D ]) ;

disp ('A = ') ; disp (A) ;

L = chol(A)' ;

disp ('L = ') ; disp (L) ;

clf
subplot (1,2,1) ; spy (A) ;
subplot (1,2,2) ; spy (L+L') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [angle_rot, bCheck] = WardSimpleRot(img1, img2)
%
%
%       [angle_rot, bCheck] = WardSimpleRot(img1, img2)
%
%       This function computes the Ward's MTB.
%
%       Input:
%           -img1: the target image
%           -img2: the image that needs to be aligned to img1
%
%       Output:
%           -rot: rotation angle (degree) for aligning img2 into img1.
%
%     Copyright (C) 2013-15  Francesco Banterle. A big thank to Greg J. Ward
%     for help during the implementation.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[r, c, ~] = size(img1);

blocksY = 3;
blocksX = 4;
sizeY = round(r / blocksY);
sizeX = round(c / blocksX);

angle = [];

%Analyzing blocks
for i=1:blocksY
    rect  = [sizeY * (i - 1) + 1, sizeY * i, 1, sizeX];
    [tmpAngle, check] = WardSimpleRotAux(img1, img2, rect);
    if(check == 1)
        angle = [angle, tmpAngle];
    end
end

%Final Merging
if(isempty(angle))
    angle_rot = 0.0;
    bCheck = 0;
else
    rotThreshold = (0.07 * 180.0) / pi;
    npos = 0;
    nneg = 0;
    for i=1:length(angle)

        if(angle(i) >  rotThreshold)
            npos = npos + 1;
        end

        if(angle(i) < -rotThreshold)
            nneg = nneg + 1;
        end
    end

    if(bitand(nneg, npos))
        angle_rot = 0.0;
        bCheck = 0;
    else
        angle_rot = mean(angle);
        bCheck = 1;
    end
end

end

author: banterle


text:
function testResults = runTestsOnModel(model, microbeID, inputDataFolder)
% Part of the DEMETER pipeline. This function performs all quality
% control/quality assurance tests on a model.
%
% USAGE:
%
%   testResults = runTestsOnModel(model, microbeID, inputDataFolder)
%
% INPUTS
% modelFolder           Folder with COBRA models (draft or refined
%                       reconstructions) to analyze
% inputDataFolder       Folder with experimental data and database files
%                       to load
%
% OUTPUT
% testResults           Structure with results of the test run
%
% .. Author:
%   - Almut Heinken, 03/2021

fields = {
    'Mass_imbalanced'
    'Charge_imbalanced'
    'Mets_without_formulas'
    'Leaking_metabolites'
    'ATP_from_O2'
    'Blocked_reactions'
    'RefinedReactionsCarryingFlux'
    'BlockedRefinedReactions'
    'Incorrect_Gene_Rules'
    'Incorrect_Compartments'
    'Carbon_sources_TruePositives'
    'Carbon_sources_FalseNegatives'
    'Fermentation_products_TruePositives'
    'Fermentation_products_FalseNegatives'
    'growsOnDefinedMedium'
    'growthOnKnownCarbonSources'
    'Biomass_precursor_biosynthesis_TruePositives'
    'Biomass_precursor_biosynthesis_FalseNegatives'
    'Metabolite_uptake_TruePositives'
    'Metabolite_uptake_FalseNegatives'
    'Secretion_products_TruePositives'
    'Secretion_products_FalseNegatives'
    'Bile_acid_biosynthesis_TruePositives'
    'Bile_acid_biosynthesis_FalseNegatives'
    'Drug_metabolism_TruePositives'
    'Drug_metabolism_FalseNegatives'
    'PutrefactionPathways_TruePositives'
    'PutrefactionPathways_FalseNegatives'
    'AromaticAminoAcidDegradation_TruePositives'
    'AromaticAminoAcidDegradation_FalseNegatives'
    };

testResults=struct();
% iterate over the field
for k=1:length(fields)
    testResults.(fields{k}){1, 1} = microbeID;
end

% Load reaction and metabolite database
database=loadVMHDatabase;
try
biomassReaction = model.rxns{strncmp('bio', model.rxns, 3)};

%
%% relax enforced uptake of some vitamins-causes infeasibility problems
relaxConstraints=model.rxns(find(model.lb>0));
model=changeRxnBounds(model,relaxConstraints,0,'l');

%% Mass and charge balance
% Tests whether the reactions in the reconstruction are mass and charge balanced.
% The function will report any reactions that have imbalanced mass or charge,
% along with their reaction formulas and imbalanced elements. The function will
% also report the metabolite formulas of all metabolites that occur in any imbalanced
% reactions. Finally, any metabolites without metabolite formulas are reported.
%
% *Inputs:*
%
% * *model*: COBRA model structure
% * *excludeExchanges*: Exclude exchange reactions when reporting imblanced
% reactions (default: true).
% * *biomassReaction*: Give biomass reaction abbreviation to exclude biomass
% reaction from analysis
%
% *Outputs:*
%
% * *massImbalancedRxns*: Cell array listing (col 1) all mass imbalanced reactions
% in the model, (col 2) the reaction formulas, and (col 3) the imbalanced elements.
% * *chargeImbalancedRxns*: Cell array listing (col 1) all charge imbalanced
% reactions in the model, (col 2) the reaction formulas, and (col 3) the imblanced
% elements.
% * *imbalancedRxnMets*: Cell array listing (col 1) all metabolites involved
% in any mass or charge imblanced reaction, (col 2) the metabolite name, (col
% 3) the metabolite formula, and (col 4) the metabolite charge.
% * *metsMissingFormulas*: Cell array listing (col 1) all metabolites that do
% not have metabolite formulas and (col 2) the metabolite names.
try
[massImbalancedRxns, chargeImbalancedRxns, ~, metsMissingFormulas] = testModelMassChargeBalance(model, true, biomassReaction);
testResults.Mass_imbalanced(1, 2:size(massImbalancedRxns,1)) =massImbalancedRxns(2:end,1);
testResults.Charge_imbalanced(1, 2:size(chargeImbalancedRxns,1))= chargeImbalancedRxns(2:end,1);
testResults.Mets_without_formulas(1, 2:length(metsMissingFormulas)) = metsMissingFormulas(2:end,1);
catch
    warning('Mass and charge balance could not be tested for draft reconstructions!')
end
%%
% A metabolic reconstruction should contain no mass or charge imbalanced
% reactions. However, some reactions will always be imbalanced by definition (exchange,
% demand, sink, and biomass reactions), and some reactions will be imbalanced
% when formulas of participating metabolites are not known.
%% Leaking metabolites
% Tests whether the model can produce any metabolites when no mass enters the
% model (all exchange and sink reactions blocked).
%
% *Inputs:*
%
% * *model*: COBRA model structure
%
% *Outputs:*
%
% * *leakingMets*: Cell array listing any metabolites that can be be produced
% from nothing and any reactions that are active in the simulation (along with
% reaction formula and flux value).

leakingMets = testLeakingMetabolites(model);
testResults.Leaking_metabolites(1, 2:length(leakingMets)+1) = leakingMets;

% The model should not leak any metabolites. Leaking metabolites require
% manual refinement.

%% ATP from oxygen
% Tests if models can produce ATP and carbon just from water,
% phosphate, and oxygen without an energy source. This is
% thermodynamically infeasible and should not happen.
basicCompounds={
    'EX_h2o(e)' '-100'
    'EX_h(e)' '-100'
    'EX_pi(e)' '-10'
    'EX_o2(e)' '-10'
    };
modelTest = useDiet(model,basicCompounds);
modelTest.lb(find(strncmp(modelTest.rxns,'sink_',5)))=0;
try
    modelTest=changeObjective(modelTest,'DM_atp_c_');
catch
    modelTest = addReaction(modelTest,'DM_atp_c','atp[c] + h2o[c] -> adp[c] + h[c] + pi[c]');
end
FBA=optimizeCbModel(modelTest);
testResults.ATP_from_O2{1, 2} = FBA.f;

%% Blocked reactions
% Computes all reactions in the reconstruction that can never carry flux.
% set "unlimited" constraints
model = changeRxnBounds(model, model.rxns(strncmp('EX_', model.rxns, 3)), -1000, 'l');
model = changeRxnBounds(model, model.rxns(strncmp('EX_', model.rxns, 3)), 1000, 'u');
BlockedRxns = identifyFastBlockedRxns(model);
testResults.Blocked_reactions(1, 2:length(BlockedRxns)+1) = BlockedRxns;

%% Test whether reactions added through refinement of gene annotations carry flux
[RefinedReactionsCarryingFlux, BlockedRefinedReactions] = testRefinedReactions(microbeID, BlockedRxns);
testResults.RefinedReactionsCarryingFlux(1, 2:length(RefinedReactionsCarryingFlux)+1) = RefinedReactionsCarryingFlux;
testResults.BlockedRefinedReactions(1, 2:length(BlockedRefinedReactions)+1) = BlockedRefinedReactions;

%% Entries in gene rules that have incorrect nomenclature
incorrectGeneRules = testGeneRules(model);
testResults.Incorrect_Gene_Rules(1, 2:length(incorrectGeneRules)+1) = incorrectGeneRules;

%% Reactions and metabolites with non-microbial compartments
[incorrectRxns,incorrectMets] = validateCompartments(model);
testResults.Incorrect_Compartments(1, 2:length(incorrectRxns)+1) = incorrectRxns;
testResults.Incorrect_Compartments(1, 2:length(incorrectMets)+1) = incorrectMets;

%% Carbon sources
% Tests which carbon sources that are reported in _in vitro_ data can be taken
% up by the model. The test requires a low flux through the biomass objective
% function.
% *Note:* Because we do not have data on carbon sources that should _not_
% be taken up by the model, we can only compare those that are known to be taken
% up by the microbe. Therefore, it is not possible to evaluate true negatives
% and false positives.

[TruePositives, FalseNegatives] = testCarbonSources(model, microbeID, biomassReaction, database, inputDataFolder);
testResults.Carbon_sources_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.Carbon_sources_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;
%%
% The function |*carbonSourceGapfill|* aims to enable to uptake of all _in
% vitro_ reported carbon sources. In the case of any false negatives, the reconstruction
% needs to be manually refined to make sure that uptaken carbon source can enter
% an active metabolic pathway within the reconstruction.
%% Fermentation pathways
% Tests which fermentation pathways that are reported in _in vitro_ data are
% active in the model. The test requires a low flux through the biomass objective
% function.
% *Note:* Because we do not have data on carbon sources that should _not_
% be taken up by the model, we can only compare those that are known to be taken
% up by the microbe. Therefore, it is not possible to evaluate true negatives
% and false positives.

[TruePositives, FalseNegatives] = testFermentationProducts(model, microbeID, biomassReaction, database, inputDataFolder);
testResults.Fermentation_products_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.Fermentation_products_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;
%%
% The function |*fermentationPathwayGapfill|* aims to enable to secretion
% of all _in vitro_ reported fermentation productions. In the case of any false
% negatives, the reconstruction needs to be manually refined to make sure that
% fermentation pathways are active within the reconstruction.

%% test the defined media according to literature
[growsOnDefinedMedium,constrainedModel,growthOnKnownCarbonSources] = testGrowthOnDefinedMedia(model, microbeID, biomassReaction, inputDataFolder);
testResults.growsOnDefinedMedium{1, 2} = growsOnDefinedMedium;
testResults.growthOnKnownCarbonSources{1, 2} = growthOnKnownCarbonSources;

% In the case of any false negatives, the reconstruction needs to be manually
% refined to make sure that essential metabolites cannot be synthesized _de novo_
% by the model or that non-essential metabolites are not required in the growth
% media.

%% Metabolite consumption
% Performs an FVA and reports those metabolites (exchange reactions)
% that can be taken up by the model and should be taken up according to
% data (true positives) and those metabolites that cannot be taken up by
% the model but should be taken up according to in vitro data (false
% negatives).
% Based on literature search of reported vitamin-secreting probiotics
% performed 11/2017, and a compendium of uptake and secretion products
% (PMID:28585563)

[TruePositives, FalseNegatives] = testMetaboliteUptake(model, microbeID, biomassReaction, database, inputDataFolder);
testResults.Metabolite_uptake_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.Metabolite_uptake_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;

%% Metabolite secretion
% Performs an FVA and reports those secretion products (exchange reactions)
% that can be secreted by the model and should be secreted according to
% data (true positives) and those secretion products that cannot be secreted by
% the model but should be secreted according to in vitro data (false
% negatives).
% Based on literature search of reported vitamin-secreting probiotics
% performed 11/2017, and a compendium of uptake and secretion products
% (PMID:28585563)

[TruePositives, FalseNegatives] = testSecretionProducts(model, microbeID, biomassReaction, database, inputDataFolder);
testResults.Secretion_products_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.Secretion_products_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;

%% Bile acid biotransformation
% Performs an FVA and reports those bile acid metabolites (exchange reactions)
% that can be secreted by the model and should be secreted according to
% data (true positives) and those bile acid metabolites that cannot be secreted by
% the model but should be secreted according to in vitro data (false
% negatives).

[TruePositives, FalseNegatives] = testBileAcidBiosynthesis(model, microbeID, biomassReaction, database);
testResults.Bile_acid_biosynthesis_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.Bile_acid_biosynthesis_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;
%% Drug biotransformation
% Performs an FVA and reports those drug metabolites (exchange reactions)
% that can be secreted by the model and should be taken up and/or secreted according to
% data (true positives) and those drug metabolites that cannot be secreted by
% the model but should be secreted according to in vitro data (false
% negatives).

if exist('drugTable.txt','File')==2
    [TruePositives, FalseNegatives] = testDrugMetabolism(model, microbeID, biomassReaction, database);
    testResults.Drug_metabolism_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
    testResults.Drug_metabolism_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;
end
%% Putrefaction pathways
% Performs an FVA and reports those putrefaction pathway end reactions (exchange reactions)
% that can carry flux in the model and should carry flux according to
% data (true positives) and those putrefaction pathway end reactions that
% cannot carry flux in the model but should be secreted according to in
% vitro data (false negatives).

[TruePositives, FalseNegatives] = testPutrefactionPathways(model, microbeID, biomassReaction, database);
testResults.PutrefactionPathways_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.PutrefactionPathways_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;
%% Aromatic amino acid degradation
% Performs an FVA and reports those AromaticAA pathway end reactions (exchange reactions)
% that can carry flux in the model and should carry flux according to
% data (true positives) and those AromaticAA pathway end reactions that
% cannot carry flux in the model but should be secreted according to in
% vitro data (false negatives).
%
% INPUT
% model             COBRA model structure
% microbeID         Microbe ID in carbon source data file
% biomassReaction   Biomass objective functions (low flux through BOF
%                   required in analysis)
%
% OUTPUT
% TruePositives     Cell array of strings listing all aromatic amino acid
% degradation products
% that can be secreted by the model and in comnparative genomic data.
% FalseNegatives    Cell array of strings listing all aromatic amino acid
% degradation products
% that cannot be secreted by the model but should be secreted according to comparative genomic data.

[TruePositives, FalseNegatives] = testAromaticAADegradation(model, microbeID, biomassReaction, database);
testResults.AromaticAminoAcidDegradation_TruePositives(1, 2:length(TruePositives)+1) = TruePositives;
testResults.AromaticAminoAcidDegradation_FalseNegatives(1, 2:length(FalseNegatives)+1) = FalseNegatives;
end
end

author: opencobra

repo: cobratoolbox


text:
function structure = addLevelVersion(structure, level, version)
%  structure = addLevelVersion(structure, level, version)
%
% this function is used internally by TranslateSBML
%
% Takes
%
% 1. structure - any MATLAB_SBML structure
% 2. level - an integer for the value of the level field to be added
% 3. version - an integer for the value of the version field to be added
%
% Returns
%
% 1. structure - the original structure with additional fields 'level' and
%                           'version' added to this and every child
%                           structure
%
% Note:
%             The structure must contain a 'typecode' field

%<!---------------------------------------------------------------------------
% This file is part of libSBML.  Please visit http://sbml.org for more
% information about SBML, and the latest version of libSBML.
%
% Copyright (C) 2013-2018 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%     3. University of Heidelberg, Heidelberg, Germany
%
% Copyright (C) 2009-2013 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%
% Copyright (C) 2006-2008 by the California Institute of Technology,
%     Pasadena, CA, USA
%
% Copyright (C) 2002-2005 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. Japan Science and Technology Agency, Japan
%
% This library is free software; you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation.  A copy of the license
% agreement is provided in the file named "LICENSE.txt" included with
% this software distribution and also available online as
% http://sbml.org/software/libsbml/license.html
%----------------------------------------------------------------------- -->

if ~isfield(structure, 'typecode')
    return;
end;
if length(structure) ~= 1 ||strcmp(structure.typecode, 'SBML_MODEL') ~= 1
    structure = addData(structure, level, version);
end;
for i = 1:length(structure)
    f = fieldnames(structure(i));
    for j = 1:length(f)
        if isstruct(structure(i).(f{j}))
            substr = structure(i).(f{j});
            structure(i).(f{j}) = addLevelVersion(substr, level, version);
        end;
    end;
end;

function str = addData(str, l, v)
f = fieldnames(str);
if sum(ismember(f, 'level')) > 0
    return;
end;
for i=1:length(str)
    str(i).level = l;
    str(i).version = v;
end;


text:
function wex(directory)
% WEX opens the Windows Explorer in selected Directory
%
%   WEX does call the Windows Explorer in the wished Directory, which the
%   user does enter as input.
%   One Option is only to enter a single letter from the Root Drive or to
%   enter a complete Path Directory information.
%
%   Options: >> wex                   - Opens in the Current Directory
%            >> wex d                 - Opens Explorer in path 'd:\'
%            >> wex d:\applications   - Opens Path in 'd:\applications'
%            >> wex; wex('d'); wex('d:\applications') - Opens all
%

%% AUTHOR    : Frank Gonzalez-Morphy
%% $DATE     : 13-Mar-2001 18:10:27 $
%% $Revision : 1.12 $
%% DEVELOPED : R12.0
%% FILENAME  : wex.m

if ~ispc
    error('  :: Sorry these is only available under Windows !')
end

if nargin == 0
    choice = 'curdir';   % Current Directory
elseif nargin >= 1
    choice = 'adr';
end

switch choice
    case 'curdir'
        % Opens the Windows Explorer in current Directory (= pwd)
        !explorer /e, .\ &
    case 'adr'
        % If input is only one character, like C D M then the path
        % will be completed like c:\, d:\ or M:\.
        % If Input has more then just one letter, then the addrss
        % will be keep and Windows Explorer opened there.
        if length(directory) == 1
            ewrd = strcat('explorer /e,', directory);     % explorer /e,d
            exw = ':\ &';                                   % :\
            comd = [ewrd exw];                            % explorer /e,d:\
            dos(comd);
        else
            comd = strcat('explorer /e,', directory);     % explorer /e, D:\applications
            dos(comd);
        end
    otherwise
        disp('  :: Wrong Input, Type ... help wex ... for correct Options !')
end

% Created with NEWFCN.m by Frank González-Morphy
% Contact...: frank.gonzalez-morphy@mathworks.de
% ===== EOF ====== [wex.m] ======

author: Sable

repo: mcbench-benchmarks


text:
function candidates = compute_mcg_cands(params, rf, n_cands, mcg_cache_obj, D, RD)
% function candidates = compute_mcg_cands(params, rf, n_cands, mcg_cache_obj)

% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------

  n_hiers = length(params.hier_dirs);

  % Load trained random forest
  if(isstr(rf)), rf = loadvar(rf,'rf'); end

  % Load Pareto parameters
  if(isstr(n_cands)), n_cands = loadvar(n_cands,'n_cands'); end

  % Number of regions per candidate
  assert(n_hiers==size(n_cands,2));

  f_lp = mcg_cache_obj.f_lp; f_ms = full(mcg_cache_obj.f_ms); feats = mcg_cache_obj.feats;
  bboxes = mcg_cache_obj.bboxes; red_cands = mcg_cache_obj.red_cands;
  b_feats_intersections = mcg_cache_obj.b_feats_intersection;

  if(params.depth_features)
    tt = tic();
    sp = f_lp;
    sp2regC = cands2labels(red_cands, f_ms);
    sp2reg = false(max(sp(:)), length(sp2regC));
    for j = 1:length(sp2regC), sp2reg(sp2regC{j},j) = true; end
    D = double(D)./1000;
    C = params.camera_matrix;
    missingMask = RD == 0;
    fdepth = depthFeatures(sp, sp2reg, D, missingMask, C);
    feats = cat(2, feats, fdepth');
    fprintf('Time for depth features: %0.3f\n', toc(tt));
  end

  % Rank candidates
  class_scores = regRF_predict(feats,rf);
  [scores, ids] = sort(class_scores,'descend');
  red_cands = red_cands(ids,:);
  bboxes = bboxes(ids,:);
  if isrow(scores), scores = scores'; end

  % Max margin
  candidates=[];
  [new_ids, mm_scores] = mex_max_margin(red_cands-1,scores, full(b_feats_intersections), params.theta); %#ok<NASGU>
  cand_labels = red_cands(new_ids,:);
  candidates.scores = scores(new_ids);
  bboxes = bboxes(new_ids,:);

  % Change the coordinates of bboxes to be coherent with
  % other results from other sources (sel_search, etc.)
  candidates.bboxes = [bboxes(:,2) bboxes(:,1) bboxes(:,4) bboxes(:,3)];

  % Get the labels of leave regions that form each candidates
  candidates.superpixels = f_lp;
  candidates.labels = cands2labels(cand_labels,f_ms);
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function anal = recomputeChopTseriesMulti(tSeries);
% anal = recomputeChopTseriesMulti(tSeries);
%
% For time course UI, put up a dialog to get
% parameters such as baseline period, peak period,
% onsetDelta, etc., for chopping up a tSeries
% according to condition onsets, then re-chop up
% the time course and add it back to the UI.
%
% If the optional useDefaults flag is set to 1 (default
% 0), then it will forego the dialog and just recompute with
% the default values.
%
% 07/04 ras.
if nargin < 2
    useDefaults = 0;
end
defaults{1} = num2str(tc.timeWindow);
defaults{2} = num2str(tc.bslPeriod);
defaults{3} = num2str(tc.peakPeriod);
defaults{4} = num2str(tc.settings.alpha);
defaults{5} = num2str(tc.settings.onsetDelta);
defaults{6} = num2str(tc.settings.snrConds);
defaults{7} = num2str(tc.settings.glmHRF);
defaults{8} = num2str(tc.settings.normBsl);

if useDefaults==0
	prompt{1} = 'Time Window, in seconds (incl. pre-stimulus onset period):';
	prompt{2} = 'Baseline Period, in seconds:';
	prompt{3} = 'Peak Period, in seconds:';
	prompt{4} = 'Alpha for significant activation (peak vs baseline):';
	prompt{5} = 'Shift onsets relative to time course, in seconds:';
    prompt{6} = 'Use which conditions for computing SNR?';
    prompt{7} = 'HRF to use for GLM? (0 = Boynton gamma, 1 = avg TC for ROI):';
	prompt{8} = 'Normalize all trials during the baseline period? (0 for no, 1 for yes)';

	AddOpts.Resize = 'on';
	AddOpts.Interpreter = 'tex';
	AddOpts.WindowStyle = 'Normal';
	answers = inputdlg(prompt,'Re-chop tSeries Settings...',1,defaults,AddOpts);
else
    answers = defaults;
end

% exit if cancel is selected
if isempty(answers)
    return;
end

% parse the user responses / defaults
twin = str2num(answers{1});
bsl = str2num(answers{2});
peak = str2num(answers{3});
alpha = str2num(answers{4});
onsetDelta = str2num(answers{5});
snrConds = str2num(answers{6});
glmHRF = str2num(answers{7});
normBsl = str2num(answers{8});

anal = er_chopTSeries2(tc.wholeTc,tc.trialInfo,'timeWindow',twin,...
                       'bslPeriod',bsl,'peakPeriod',peak,...
                       'alpha',alpha,'onsetDelta',onsetDelta,...
                       'snrConds',snrConds,'normBsl',normBsl);

% assign over to tc struct (ugly, but there are other fields
% in tc I don't want to erase -- is there a way to eat the whole
% struct in one go?)
tc.wholeTc = anal.wholeTc;
tc.allTcs = anal.allTcs;
tc.meanTcs = anal.meanTcs;
tc.sems = anal.sems;
tc.timeWindow = anal.timeWindow;
tc.peakPeriod = anal.peakPeriod;
tc.bslPeriod = anal.bslPeriod;
tc.Hs = anal.Hs;
tc.ps = anal.ps;
tc.amps = anal.amps;
tc.relamps = anal.relamps;
tc.SNR = anal.SNR;
tc.SNRdb = anal.SNRdb;
tc.settings.alpha = alpha;
tc.settings.onsetDelta = onsetDelta;
tc.settings.snrConds = snrConds;
tc.settings.glmHRF = glmHRF;
tc.settings.normBsl = normBsl;

% update dataTYPES parameters for these scans
if isfield(tc.params,'viewName')
	mrGlobals;
    tmp = findstr('{',tc.params.viewName);
    varName = tc.params.viewName(1:tmp(1)-1);
    if exist(varName,'var')
        view = eval(tc.params.viewName);

		settings.timeWindow = anal.timeWindow;
		settings.peakPeriod = anal.peakPeriod;
		settings.bslPeriod = anal.bslPeriod;
		settings.normBsl = normBsl;
		settings.alpha = alpha;
		settings.onsetDelta = onsetDelta;
		settings.glmHRF = glmHRF;

		for s = tc.params.scans
            er_setParams(view,settings,s);
        end
    end
end

% assign back to the UI and refresh
set(gcf,'UserData',tc);
timeCourseUI;

return

author: vistalab

repo: vistasoft


text:
function qmf = wavelet_qmf(w)
	qmf = w.qmf;

author: thomaskuestner


text:
function plotlowerbound(lb,peval,savethis, xlow)
% plotlowerbound(lb,peval,savethis, xlow)
% xlow : (optional) lower bound on x axes (not to show first few pointss in
% lb...)
if ~exist('xlow', 'var') xlow = 1; end

figure
plot(lb(xlow:end))
xlabel('iteration');
ylabel('lower bound');
grid on
setfontsizefigure(12)
if savethis
    SaveImageFULL([peval.res_dir '/lowerbound'])
end

author: aludnam