text:
% Select target active features.

% Input:
%   feature_groups    -  a cell array with each cell is one kind of
%                        features of the search image window
%   filter_sz         -  size of the filter [height wieth depth], also the
%                        target size in the feature map.
%   sel_feat_num      -  selected feature numbers

% Output:
%   feat_weight -   weights of each features.
%   IAF_model   -   initialized IAF model

% by Xin Li, April 2 2018
function [feat_weights,w_balance] = TAF_model( feature_groups, filter_sz)

assert(iscell(feature_groups),'the input feature should be a cell array!\n');
feat_num     = numel(feature_groups);
feat_weights = cell(1,feat_num);
channel_num  = [300,80]; % we select 300 and 80 filters from the conv43 and conv41 layer features.
nz_num       = 0;
nz_num_min   = 250;       % the minimum filter number.
w_balance    = gpuArray(single([1 1 1])); % weights leveraging these two features.

%peform feature selection on the conv41 and conv43 layer features.
for layer_i=feat_num:-1:1
    %% initialization
        filter_sz(3) = size(feature_groups{layer_i},3);
        [model.net] = TAF_net_init(filter_sz);

        %generate the Gaussian label map
        output_sigma = filter_sz([1:2])*0.1; % output_sigma_factor=0.1
        [h,w,~] = size(feature_groups{layer_i}); sw_sz = [h w]; % search window size [height width]
        model.label_gaussian = single(gaussian_shaped_labels(output_sigma, sw_sz));
        model.label_binary = single(binary_labels(filter_sz([1:2]), sw_sz, [0.3 0.7]));
    %% First train the network with ridge loss
        trainOpts.numEpochs=100;
        trainOpts.learningRate=5e-7;
        trainOpts.weightDecay= 1;
        trainOpts.batchSize = 1 ;
        trainOpts.gpus = 1 ;
        trainOpts.train = 1;
        input = {feature_groups{layer_i}, gpuArray(model.label_gaussian)};

        [net_trained, ~, ders_iter] = cnn_train_dag_ridge(model.net, [], input, [], trainOpts);
        [v,i] = sort(sum(sum(net_trained.params(1).value)),'descend'); % GAP.
        % The value of the parameters equals to the sum of the gradients in
        % all BP processes. And we found that using the converged parameters is more
        % stable.
    %% mapping gradients values to weights
        feat_weight=zeros(1,1,filter_sz(3), 'like', feature_groups{layer_i});
        feat_weight(i(v>0))=1;     % here we select the filters with positive gradient values, which have a positive relationship with the target loss.
        if layer_i>1 % we perform scale sensitive feature selection on the conv41 feature, as it retains more spatial information.
            filter_sz(3) = size(feature_groups{layer_i},3);
            % initialize the ranking loss network
            [IAF_net_rank, filter_sz_new] = TAF_net_init_rank(filter_sz);
            % compute the scale sensitive features
            [temp_weight, temp_i] = rank_selection(feature_groups{layer_i}, filter_sz_new, IAF_net_rank);
            % select the features not only are active to the target but also are scale-sensitieve.
            feat_weight = feat_weight.*temp_weight;
        end
        feat_weight(i(channel_num(layer_i)+1:end)) = 0; % we only select the first N (channel_num) filters.
        nz_num=nz_num+sum(feat_weight);

        % In case, there are too less features, we set a minimum feature
        % number. If the total number is less than the minimum number, then select more from conv43
        if layer_i==1 && nz_num<nz_num_min
            added_inds=(i(sum(feat_weight)+1:sum(feat_weight)+nz_num_min-nz_num));
            feat_weight(added_inds) = 1;
        end

        % To combine the features from conv41 and conv43 effective, we define
        % balance weights based on the maximum value in their features.
        feat_weights{layer_i}=feat_weight;
        w_balance(layer_i) = max(max(sum(feature_groups{layer_i}(:,:,i(1:50)),3)));
        clear net_trained;

end
        w_balance=w_balance(1)./w_balance;
end

function [feat_weight, i_index] = rank_selection(feat, filter_sz, net)
% this function selects the scale sensitve features based on the ranking
% loss.

    channel_num= size(feat,3);
    [scale_samples, pair_label, pairs, locs] = gen_ranked_samples(feat, filter_sz(1:2));

% gradients of samples with different scales
    inputs = {'input', gpuArray(scale_samples), 'scale_label', gpuArray(pair_label)};
    net.eval(inputs,{'objective_scale',1});
    feat_der = sum(net.vars(1).der,4);

    [v_scale,i_index] = sort(sum(sum(sum(feat_der)),4), 'descend');
    feat_weight=ones(1,1,channel_num, 'like', feat)*0;
    feat_weight(i_index(v_scale>0))=1;  % we select the features, which have positive relationships with the ranking loss.
end

function [IAF_net_rank,filter_sz] = TAF_net_init_rank(filter_sz)
rng('default');

channel=filter_sz(3);

rw=ceil(filter_sz(2)/2);
rh=ceil(filter_sz(1)/2);

fw=2*rw+1;
fh=2*rh+1;

filter_sz = [fh,fw,channel];

IAF_net_rank=dagnn.DagNN();

%% rank loss net

IAF_net_rank.addLayer('conv12', dagnn.Conv('size', [fw,fh,channel,1],...
    'hasBias', true, 'pad',...
    [0,0,0,0], 'stride', [1,1]), 'input', 'conv12', {'conv21_f', 'conv21_b'});

f = IAF_net_rank.getParamIndex('conv21_f') ;
IAF_net_rank.params(f).value=single(randn(fh,fw,channel,1) /...
    sqrt(rh*rw*channel))/1e8;
IAF_net_rank.params(f).learningRate=1;
IAF_net_rank.params(f).weightDecay=1e3;

f = IAF_net_rank.getParamIndex('conv21_b') ;
IAF_net_rank.params(f).value=single(zeros(1,1));
IAF_net_rank.params(f).learningRate=2;
IAF_net_rank.params(f).weightDecay=1e3;
IAF_net_rank.addLayer('rank_loss',...
    RankingLoss(),{'conv12','scale_label'},'objective_scale');

IAF_net_rank.move('gpu');
end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function [newMap] = changeRxnColorAndWidth(map, rxnList, color, areaWidth)
% Change color and areaWidth of reactions from a list of names
%
% USAGE:
%
%   [newMap] = changeRxnColorAndareaWidth(map, rxnList, color, areaWidth);
%
% INPUTS:
%   map:            File from CellDesigner parsed to MATLAB format
%   rxnList:        List of reactions
%
% OPTIONAL INPUTS:
%   color:          New color of reactions from list (default: 'RED')
%   areaWidth:      New areaWidth of reactions from list (default: 8)
%
% OUTPUT:
%   newMap:         Matlab structure of map with reaction modifications
%
% .. Authors:
%       - A.Danielsdottir 17/07/2017 LCSB. Belval. Luxembourg
%       - N.Sompairac - Institut Curie, Paris, 17/07/2017.

    if nargin < 4
        areaWidth = 8;
    end
    if nargin < 3
        color = 'RED';
    end

    newMap = map;
    colors = createColorsMap;

    index = find(ismember(newMap.rxnName, rxnList));
    for j = index'
        newMap.rxnColor{j, 1} = colors(color);
        newMap.rxnWidth{j, 1} = areaWidth;

        % Use the existence of reactant lines to check if the newMap has the
        % complete structure.
        if any(strcmp('rxnReactantLineColor', fieldnames(newMap))) == 1
            if ~isempty(newMap.rxnReactantLineColor{j})
                for k = 1:length(newMap.rxnReactantLineColor{j})
                    newMap.rxnReactantLineColor{j, 1}{k, 1} = colors(color);
                    newMap.rxnReactantLineWidth{j, 1}{k, 1} = areaWidth;
                end
            end
            if ~isempty(newMap.rxnProductLineColor{j})
                for m = 1:1:length(newMap.rxnProductLineColor{j})
                    newMap.rxnProductLineColor{j, 1}{m, 1} = colors(color);
                    newMap.rxnProductLineWidth{j, 1}{m, 1} = areaWidth;
                end
            end
        end
    end

end

author: opencobra

repo: cobratoolbox


text:
function [vertices, faces] = readMesh_off(fileName)
%READMESH_OFF Read mesh data stord in OFF format
%
%   [VERTICES FACES] = readMesh_off(FILNAME)
%
%   Example
%     [v, f] = readMesh_off('mushroom.off');
%     figure; drawMesh(v, f, 'faceColor', [0 1 0], 'edgeColor', 'none')
%     view([5 80]); light; lighting gouraud
%
%   See also
%     meshes3d, drawMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-12-20,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% open file
f = fopen(fileName, 'r');
if f == -1
    error('matGeom:readMesh_off:FileNotFound', ...
        ['Could not find file: ' fileName]);
end

% check format
line = fgets(f);   % -1 if eof
if ~strcmp(line(1:3), 'OFF')
    error('matGeom:readMesh_off:FileFormatError', ...
        'Not a valid OFF file');
end

% number of faces and vertices
line = fgets(f);
vals = sscanf(line, '%d %d');
nv = vals(1);
nf = vals(2);

% read vertex data
[vertices, count] = fscanf(f, '%f ', [3 nv]);
if count ~= nv*3
    error('matGeom:readMesh_off:FileFormatError', ...
        ['Could not read all the ' num2str(nv) ' vertices']);
end
vertices = vertices';

% read face data (face start by index)
[faces, count] = fscanf(f, '%d %d %d %d\n', [4 nf]);
if count ~= nf * 4
    error('matGeom:readMesh_off:FileFormatError', ...
        ['Could not read all the ' num2str(nf) ' faces']);
end

% clean up: remove index, and use 1-indexing
faces = faces(2:4, :)' + 1;

% close the file
fclose(f);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
imageRootPath = '/n/fs/vision/datasets/SUN397/';

%% get image list
disp('generating image list...');
folders = genpath(imageRootPath);
folders = regexp(folders,':','split');
folders = folders(1:end-1);
imageList = {};
for f=1:length(folders)
    imagefiles = dir(fullfile(folders{f},'*.jpg'));
    if ~isempty(imagefiles)
        for i=1:length(imagefiles)
            imageList{end+1} = fullfile(folders{f},imagefiles(i).name);
        end
    end
end

%% read image, crop, and resize
disp('reading and resizing images...');

tWidth =32;
tHeight=32;
imArray = uint8(zeros(tHeight,tWidth,3,length(imageList)));

for i=1:length(imageList)
    im = imread(imageList{i});
    w = size(im,2);
    h = size(im,1);
    ch = w/tWidth*tHeight;
    if ch<h
        srow = max(1,floor(h/2-ch/2));
        erow = srow + round(ch)-1;
        im = im(srow:erow,:,:);
    elseif ch>h
        cw = h/tHeight*tWidth;
        if cw<w
            scol = max(1,floor(w/2-cw/2));
            ecol = scol + round(cw)-1;
            im = im(:,scol:ecol,:);
        end
    end

    if size(im,3)==1
        im = im(:,:,[1 1 1]);
    end

    imArray(:,:,:,i) = imresize(im,[tHeight tWidth]);

    if mod(i,100)==1
        fprintf('%d/%d\n',i,length(imageList));
    end
end

save('sun397.mat','-v7.3');

%%
%
load sun397.mat
imArray = double(imArray)/255;

baseImage = 'base9.png';
baseImage = im2double(imread(baseImage));
baseImage = baseImage(:,:,1:3);
nH = floor(size(baseImage,1)/tHeight);
nW = floor(size(baseImage,2)/tWidth);
baseImage = baseImage(1:nH*tHeight,1:nW*tWidth,:);

nI = size(imArray,4);
usedVector = false(1,nI);

mosaicImage = zeros(nH*tHeight,nW*tWidth,3);

for y=randperm(nH)
    ey = y*tHeight;
    sy = ey-tHeight+1;
    for x=randperm(nW)
        ex = x*tWidth;
        sx = ex-tWidth+1;

        sourcePath = baseImage(sy:ey,sx:ex,:);

        diff = sum(sum(sum(abs(imArray - repmat(sourcePath,[1,1,1,nI])),1),2),3);
        diff = reshape(diff,1,[]) + usedVector .* 100000;

        [minVal,minIdx]= min(diff);

        usedVector(minIdx) = true;

        mosaicImage(sy:ey,sx:ex,:) = imArray(:,:,:,minIdx);
    end
end

imwrite(mosaicImage,'mosaic9.png');
imshow(mosaicImage);

author: jianxiongxiao

repo: ProfXkit


text:
function files_out = read_folder_structure(root,files_in)

[curr_files,isdir] = my_dir(root);

files_out = files_in;
for i=1:length(curr_files)
    if isdir(i) == 1
        files_out = read_folder_structure([root '\' curr_files{i}],files_out);
    else
        files_out{length(files_out)+1} = [root '\' curr_files{i}];
    end
end

author: singaxiong

repo: SignalGraph


text:
function [ outp, sdp_mode ] = cvx_pushcnstr( x, y, op, sdp_mode )

global cvx___
try
    pstr = cvx___.problems(end);
catch
    cvx_throw( 'No CVX model is present.' );
end

%
% Check problem
%

y_orig = y;
if isa( x, 'cvxcnst' ),
    x = rhs( x );
end

%
% Check for a dual reference
%

dx = cvx_getdual( x );
dy = cvx_getdual( y );
if isempty( dx ),
    dx = dy;
elseif ~isempty( dy ),
    cvx_throw( 'Two dual variable references found: %s %s', dx, dy );
end
if ~isempty( dx ),
    duals = pstr.duals;
    try
        dual = builtin( 'subsref', duals, dx );
    catch
        nm = cvx_subs2str( dx );
        cvx_throw( 'Dual variable "%s" has not been declared.', nm(2:end) );
    end
    if ~isempty( dual ),
        nm = cvx_subs2str( dx );
        cvx_throw( 'Dual variable "%s" is already in use.', nm(2:end) );
    end
end

%
% Check arguments
%

cx = isnumeric( x ) | isa( x, 'cvx' );
cy = isnumeric( y ) | isa( y, 'cvx' );
if ~cx,
    x = cvx_collapse( x, false, true );
    cx = isnumeric( x ) | isa( x, 'cvx' );
end
if ~cy,
    y = cvx_collapse( y, false, true );
    cy = isnumeric( y ) | isa( y, 'cvx' );
end
if ~cx || ~cy,
    if cx || cy || op(1) ~= '=',
        cvx_throw( 'Invalid CVX constraint: {%s} %s {%s}', class( x ), op, class( y ) );
    end
    sx = size( x );
    if ~isequal( sx, size( y ) ),
        cvx_throw( 'The left- and right-hand sides have incompatible sizes.' );
    else
        if ~isempty( dx ),
            pstr.duals = builtin( 'subsasgn', pstr.duals, dx, cell(sx) );
            cvx___.problems( end ) = pstr;
        end
        for k = 1 : prod( sx ),
            cvx_pushcnstr( x{k}, y{k}, op );
        end
        if nargout,
            outp = cvxcnst( y_orig );
        end
        return
    end
end

persistent param_eq param_ge param_le
if isempty( param_eq ),
    param_eq.map = cvx_remap( ...
            { { 'affine',  } }, ...
            { { 'l_affine' } } );
    param_eq.funcs = { @minus_nc, @eq_2 };
    param_eq.name = '=';
    param_eq.constant = [];
    param_ge.map = cvx_remap( ...
            { { 'constant' } }, ...
            { { 'l_concave' }, { 'l_convex' } }, ...
            { { 'concave' }, { 'convex' } }, ...
            [ 1, 2, 1 ] );
    param_ge.funcs = { @minus_nc, @ge_2 };
    param_ge.name = '>=';
    param_ge.constant = false;
    param_le.map = param_ge.map';
    param_le.funcs = { @le_1, @le_2 };
    param_le.name = '<=';
    param_le.constant = false;
end

if nargin < 5,
    sdp_mode = pstr.sdp;
end
switch op(1),
case '=',
    P = param_eq;
    sdp_mode = false;
case '<',
    P = param_le;
case '>',
    P = param_ge;
end
P.name = op;
P.sdp = sdp_mode;
[ z, sdp_mode ] = cvx_binary_op( P, x, y );

%
% Handle LMIs
%

zR = [];
sz = size( z );
zr = isreal( z );
if sdp_mode,
    z  = cvx_basis( z );
    dg = false;
    if zr,
        n  = sz(1);
        zq = any( z, 1 );
        nq = nnz( zq );
        if nq <= n,
            qn = bsxfun( @plus, (1:n+1:n*n)', 0:n*n:prod(sz)-1 );
            dg = nq == nnz( zq( qn ) );
            if dg,
                [ tx, dummy ] = find( cvx_basis( nonnegative( nq ) ) ); %#ok
                zz = sparse( tx, qn(zq(qn)), 1, tx(end), prod(sz) );
            end
        end
    end
    if ~dg,
        zt = z(:,permute(reshape(1:prod(sz),sz),[2,1,3:length(sz)]));
        if ~zr, zt = conj( zt ); end
        err = z - zt;
        if nnz(err),
            z   = 0.5 * ( z + zt );
            err = max( sum(abs(err),1) ./ max(realmin,sum(abs(z),1)) );
            if err > 8 * eps,
                if isreal(z), str = 'symmetric'; else str = 'Hermitian'; end
                cvx_throw( [ ...
                    'SDP constraints are expected to be %s (deviation: %g).\n', ...
                    '--- If this number is small (<1e-6), it may simply be due to roundoff error.\n', ...
                    '    This can be corrected by applying the SYM(X) function.\n', ...
                    '--- Otherwise, this is likely due to a modeling error. Did you declare the\n', ...
                    '    relevant matrix variables to be "symmetric" or "hermitian"?' ], str, full(err) );
            end
        end
        zz = cvx_basis( semidefinite( sz, ~zr ) );
    end
    z(size(zz,1),:) = 0;
    z = z - zz;
    iseq = true;
    if ~dg,
        [ zR, z ] = cvx_bcompress( z, 'full' );
    end
else
    iseq = op(1) == '=';
    z = cvx_basis( z );
end

%
% Add the (in)equalities, create the dual variable if request
%

need_dual = ~isempty( dx );
zD = cvx_newcnstr( z, iseq, need_dual );
if need_dual
    if isempty( cvx_getdual( z ) ),
        if ~isempty( zR ), zD = zD * cvx_invert_structure( zR ); end
        pstr.duals = builtin( 'subsasgn', pstr.duals, dx, cvx( sz, zD ) );
        cvx___.problems( end ) = pstr;
    else
        warning( 'CVX:NoDualVariables', 'Dual variables are not availble for log-convex constraints.' )
    end
end

%
% Create the output object
%

if nargout,
    outp = cvxcnst( y_orig );
end

function z = eq_2( x, y )
z = cvx_setdual( minus_nc( log( x ), log( y ) ), 1 );

function z = ge_2( x, y )
z = cvx_setdual( minus_nc( log( x ), log( y ) ), 1 );

function z = le_1( x, y )
z = minus_nc( y, x );

function z = le_2( x, y )
z = cvx_setdual( minus_nc( log( y ), log( x ) ), 1 );

% 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:
function H = RS(sequence,isplot)
%
% 'RS' estimate the hurst parameter of a given sequence with R/S method.
%
% Inputs:
%     sequence: the input sequence for estimate
%     isplot: whether display the plot. without a plot if isplot equal to 0
% Outputs:
%     H: the estimated hurst coeffeient of the input sequence

%  Author: Chu Chen
%  Version 1.0,  03/10/2008
%  chen-chu@163.com
%

if nargin == 1
    isplot = 0;
end

N = length(sequence);
dlarge = floor(N/5);
dsmall = max(10,log10(N)^2);
D = floor(logspace(log10(dsmall),log10(dlarge),50));
D = unique(D);
n = length(D);
x = zeros(1,n);
y = zeros(1,n);

R = cell(1,n);
S = cell(1,n);
for i = 1:n
    d = D(i);
    m = floor(N/d);
    R{i} = zeros(1,m);
    S{i} = zeros(1,m);
    matrix_sequence = reshape(sequence(1:d*m),d,m);

    Z1 = cumsum(matrix_sequence);
    Z2 = cumsum(repmat(mean(matrix_sequence),d,1));
    R{i} = (max(Z1-Z2)-min(Z1-Z2));
    S{i} = std(matrix_sequence);

    if min(R{i})==0 || min(S{i}) ==0
        continue;
    end

    x(i) = log10(d);
    y(i) = mean(log10(R{i}./S{i}));
end

% fit a line with middle part of sequence
index = x~=0;
x = x(index);
y = y(index);
n2 = length(x);
cut_min = ceil(3*n2/10);
cut_max = floor(9*n2/10);

X = x(cut_min:cut_max);
Y = y(cut_min:cut_max);
p1 = polyfit(X,Y,1);
Yfit = polyval(p1,X);
H = (Yfit(end)-Yfit(1))/(X(end)-X(1));

if isplot ~= 0
    figure,hold on;
    bound = ceil(log10(N));
    axis([0 bound 0 0.75*bound]);

    temp = (1:n).*index;
    index = temp(index);
    for i = 1:n2
        plot(x(i),log10(R{index(i)}./S{index(i)}),'b.');
    end

    x = linspace(0,bound,10);
    y1 = 0.5*x;
    y2 = x;
    h1 = plot(x,y1,'b--','LineWidth',2);
    h2 = plot(x,y2,'b-.','LineWidth',2);
    plot(X,Yfit,'r-','LineWidth',3);
    legend([h1,h2],'slope 1/2','slope 1',4)
    xlabel('log10(blocks of size m)'),ylabel('log10(R/S)'),title('R/S Method');
end

author: Sable

repo: mcbench-benchmarks


text:
%rf_goia.m
%Jamie Near, McGill University 2020.  Based on old code from 2007.
%
% USAGE:
% [RF,FM,mv,sc]=rf_goia(N,Tp,dx,tbw,gmax,xshift)
%
% DESCRIPTION:
% this funciton creates a GOIA gradient modulated adiabatic pulse, using
% the method of gradient modulated offset-independent adiabaticity as first
% described by Tannus and Garwood in NMR Biomed 1997; 10:423-434.
%
% INPUTS:
% N              = Number of points in RF waveform (must be an even number).
% Tp             = Duration of the RF pulse in [ms].
% dx             = Desired slice thickness in [cm]
% tbw            = Time bandwidth product.
% gmax           = Maximum allowed gradient strength [G/cm] (Optional.
%                  Default = 4 G/cm == 40 mT/m)
% xshift         = Desired shift of the selected slice from isocentre in [cm].
%                  (Optional. Default=0);
%

% OUTPUTS:
% RF             = Output rf waveform for GOIA pulse, in FID-A rf pulse
%                  structure format.
% FM             = Frequency modulation waveform (in Hz).
% mv             = Simulated magnetization vector in three columns (x,y,z)
%                  as a function of frequency.
% sc             = Frequency scale (in kHz) corresponding to the simulated
%                  mv vectors.

function [RF,FM,mv,sc]=rf_goia(N,Tp,dx,tbw,gmax,xshift)

if nargin<7
    xshift=0;
    if nargin<6
        gmax=4;
    end
end

%convert Tp from ms to s;
Tp=Tp/1000;

%make sure N is even
if mod(N,2)~=0
    N=N+1;
end

%initialize the time vectors
%ta has N steps from 0 to Tp.
t=[0:Tp/(N-1):Tp];

%tau has N steps from -1 to 1. (useful for defining our AM and GM
%functions.)
tau=t*2/Tp-1;

%create time vector from 0 to 1 that is N/2 in length (useful for creation
%of FM).
tau2=tau(end/2+1:end);

%create truncation factor.
B=asech(0.01);

%Find Bandwith Factor A
bw=tbw/Tp;
A=bw/2;

%find time step size
dt=t(2)-t(1);

%define Gyromagnetic Ratio (Hz/G)
gyro=4257.7;

%First define the AM function:
F1=sech(B*(tau.^4));

%now calculate the FM function based on the assumption of a constant
%gradient by integrating the AM function:
F2a=zeros(1,length(F1));
F2a(length(F1)/2+1:length(F1))=cumsum(F1(length(F1)/2+1:length(F1)));
F2a(1:length(F1)/2)=-F2a(end:-1:length(F1)/2+1);
F2a=F2a/max(F2a);
%plot(tau,F2a);

%calculate a static gradient value based on desired slice thickness and
%bandwidth.  Ans in Gauss/cm.
G=bw/(gyro*dx);

%calculate an x of t function based on a constant gradient using the
%desired slice thickness.
xoft=A*F2a/(gyro*G);

%now offset the x of t function by the amount given by xoffset:
xoft=xoft+xshift;
%figure
%plot(tau,xoft);

%now define the GM function:
F3=(1-0.90*sech(B*(tau.^2)));

%now define the FM function that was derived from the solution to the
%differential equation in Tannus et al, NMR Biomed, 10:423-434 (1997).  The
%differential equation was solved by Jamie Near using Maple.
G=cosh(B*(tau2.^2))./((5*cosh((B*(tau2.^4))+(B*(tau2.^2))))+(5*cosh((B*(tau2.^4))-(B*(tau2.^2))))-(9*cosh((B*(tau2.^4)))));
h=(10*cosh(B*(tau2.^2))-9)./cosh(B*(tau2.^2));

%the above funciton G must be integrated from 0 to t and divided by the
%expression h to obtain the value of the FM function at time t for t>0.
%For t<0, simply mirror the function about the origin.

%integrate first half of G and multiply by h.
halfg=cumsum(G);
halfg=halfg*(tau2(2)-tau2(1)).*h;

fullg=[-halfg(end:-1:1) halfg];

F2=fullg/max(fullg);

%now calculate max gradient strength based on slice thickness desired: (assume
%slice located at x=0.
Ag=A.*(1/gyro).*(1/(dx/2));

%GM is the maximum gradient value, and it cannot exceed 4G/cm.  Therefore,
%check if AG is greated than 4, and if so, ask permission to reduce the A
%value of the pulse.
if Ag >= gmax
    chg = input('Gradient too high for such a narrow slice.  Allow reduction of tbw? ([y] or n):  ','s');
    if isempty(chg)
        chg='y';
    end
    if chg == 'y' || chg == 'Y'
        A = gmax*gyro*(dx/2);
        disp(['Reducing max gradient to' num2str(gmax) ' G/cm.']);
        tbw=2*A*Tp/1000;
    end
end

GM=Ag*F3;
AM=F1;
FM=(A*F2);

%create new FM function based on x of t function, and gradient function:
FM2=xoft.*gyro.*GM;

%If xshift is not zero, then we need to adjust the FM function.
FM=FM +(GM*gyro*xshift);

%now create phase modulation function using FM
ph=cumsum(FM)*dt*360;

rf(:,1)=ph;
rf(:,2)=AM;
rf(:,3)=1;
rf(:,4)=GM;

%Now find the b1max required to get full inversion:
%Since the pulse is phase modulated, so we will need to run some test to find
%out the w1max;  To do this, we can plot Mz as a function of w1 and
%find the value of w1 that results in the desired flip angle.
Tp_sim=0.005;
[mv,sc]=bes(rf,Tp_sim*1000,'b',0,0,5,40000);
plot(sc,mv(3,:));
xlabel('w1 (kHz)');
ylabel('mz');
w1max=input('Input desired w1max in kHz:  ');
w1max=w1max*1000; %convert w1max to [Hz]
tw1=Tp_sim*w1max;

RF.waveform=rf;
RF.type='inv';
RF.f0=xshift;
RF.tw1=tw1;
RF.tbw='N/A - gradient modulated pulse';
RF.isGM=true;
RF.tthk=Tp*dx;

[sc,mv]=bes(RF.waveform,Tp*1000,'f',tw1/Tp/1000,xshift-(2*dx),xshift+(2*dx),40000);

author: CIC-methods

repo: FID-A

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


text:
function gf=comp_filterbankresponse(g,a,L,do_real)

M=numel(g);

if size(a,2)>1
    % G1 is done this way just so that we can determine the data type.
    G1=comp_transferfunction(g{1},L);
    gf=abs(G1).^2*(L/a(1,1)*a(1,2));

    for m=2:M
        gf=gf+abs(comp_transferfunction(g{m},L)).^2*(L/a(m,1)*a(m,2));
    end;

else
    % G1 is done this way just so that we can determine the data type.
    G1=comp_transferfunction(g{1},L);
    gf=abs(G1).^2*(L/a(1));

    for m=2:M
        gf=gf+abs(comp_transferfunction(g{m},L)).^2*(L/a(m));
    end;

end;

if do_real
    gf=gf+involute(gf);
end;

gf=gf/L;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% A demo for evaluation on SYSU-MM01 dataset using features learned by deep
% zero padding in "RGB-Infrared Cross-Modality Person Re-identification"

% Features of each cameras are saved in seperated mat files named "name_cam#.mat"
% In the mat files, feature{id}(i,:) is a row feature vector of the i-th image of id
feature_info.name = 'feat_deep_zero_padding';
feature_info.dir = './feature';
result_dir = './result'; % directory for saving result

setting.mode = 'all_search'; %'all_search' 'indoor_search'
setting.number_shot = 1; % 1 for single shot, 10 for multi-shot

model.test_fun = @euclidean_dist; % Similarity measurement function
% (You could define your own function in the same way as euclidean_dist function)
model.name = 'euclidean'; % model name
model.para = []; % No parameter is needed for euclidean distance here
% (If mahalanobis distance is used, the parameter can be the metric M learned from training data)

% load data split
content = load('./data_split/test_id.mat'); % fixed testing person IDs
data_split.test_id = content.id;

content = load('./data_split/rand_perm_cam.mat'); % fixed permutation of samples in each camera
data_split.rand_perm_cam = content.rand_perm_cam;

% evaluation
performance = evaluation_SYSU_MM01(feature_info, data_split, model, setting, result_dir);

author: wuancong


text:
function make_roi_paths(pial, outersmoothed, radius, stepsize, outdir, flagf)

% Original Author: Marie Schaer
% Date: 2007/11/15
%
% This function creates multiple circular regions of interest on the outer
% surface, measures their area (saved in matlab format) and save path file
% containing a set of
% corresponding points on the pial surface. The path files are required as
% an input to the mri_path2label function
%
% Parameters:
% radius: radius of the circular region of interest on the outer surface.
% Most appropriate values between 20 and 25 mm, further documentation
% available in the validation publication: "A surface-based approach to
% quantify local cortical gyrification" by Schaer M, Bach Cuadra M,
% Tamarit L, Eliez E, Thiran JP, IEEE Transactions on Medical Imaging, 2007
%
% Utilities to transfer points on the pial surface calls "mesh_vertex_nearest",
% a function written by Darren Weber. Utilities to compute matrix adjacency
% for the pial surface calls "mesh_adjacency", a program first written by
% Darren Weber, then modified by G. Peyre (source: toolbox_graph on
% mathworks)
%
% Example: make_roi_paths('lh.pial','lh.outer-pial-smoothed',25,100,/tmp)

t0 = cputime;

disp('loading datas ...')
[mesh_total.vertices, mesh_total.faces]=freesurfer_read_surf(pial);
[mesh_outer.vertices, mesh_outer.faces]=freesurfer_read_surf(outersmoothed);

disp('preparing outer mesh structure ...')
mesh_outer = createMeshFacesOfVertex (mesh_outer.vertices, mesh_outer.faces)

disp('preparing pial mesh structure ...')
A = mesh_adjacency(mesh_total);

% Keep track of the area for each individual region of interest on the
% outer surface:
nbVertices = size(mesh_outer.vertices,1);
areasOuterROI = zeros(nbVertices,1);
areasOuterROI = sparse (areasOuterROI);
clear nbVertices

% Circular regions of interest are defined each 100 vertex on the outer
% surface. Due to the high resolution of the outer mesh (average face area
% of ~0.3mm2), calculations are computed each 1 on 100 vertex, to avoid
% high redundancies and optimize calculation time.

for iV = 1:stepsize: length(mesh_outer.vertices)

    disp(['... creating path file for vertex ',num2str(iV),' / ',num2str(length(mesh_outer.vertices))])

    % ------ Part 1: find the ROI on the enveloppe --------------

    [verticesInROI, facesInROI] = getVerticesAndFacesInSphere(mesh_outer, iV, radius);

    % In rare case, the intersection of the enveloppe with the sphere result
    % in two regions of interest: the one that we expect, circular and
    % centered in iV; and an aberrant one at the bottom of the sphere
    % (e.g. near the superior sagittal vault). The next step is used to
    % control for that and keep only the radial region centered at iV.
    facesListOuterGeo = MakeGeodesicOuterROI (mesh_outer, iV, facesInROI);

    % Find the perimeter of the outer ROI
    verticesOuterGeo=unique(facesListOuterGeo(:));
    perim=setdiff(verticesOuterGeo,verticesInROI);

    % Measure its area
    areaOuterROI = getFacesArea(mesh_outer,facesListOuterGeo);

    areasOuterROI(iV) = areaOuterROI;

    if (mod(iV,50*stepsize) == 1)
        fname = [outdir '/' pial '.area_intermed.mat'];
        save(fname, 'areasOuterROI') ;
        disp(['area file for outer ROIs saved at ',num2str(iV)])
    end

    clear verticesInROI facesInROI facesListOuterGeo verticesOuterGeo areaOuterROI

    % ----- Part2: Define the corresponding set of points on the pial surface:
    % Transfer a few points of the perimeter of the hulls ROI from
    % the envelope to the pial surface, and save them as path file.
    step = 7;
    [verticeslist] = SearchProjectionOnPial(mesh_total, mesh_outer, perim, step);

    % reorganize the set of points in the right order to input them to
    % mri_path2label
    reorglist = reorganize_verticeslist (mesh_total, A, mesh_outer, perim, verticeslist, step);
    lindex = reorglist';
    lxyz = mesh_total.vertices(lindex,:); lindex= lindex-1;

    p = sprintf ('%06d', iV);
    write_path(lxyz, lindex,[outdir '/' pial '.' p '.path']);

    clear verticeslist lindex lxyz

end

disp(['saving area file in matlab format... '])
fname = [outdir '/' pial '.outer_ROIs_area.mat'];
save(fname, 'areasOuterROI')

disp('DONE.')
deltaT = cputime - t0

% indicate successful completion by deleting flagfile
delete(flagf);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

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

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=35;
lineWidth=3;
cMap=viridis(20); %colormap

%% Control parameters

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

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

%Specifying dimensions and number of elements
cubeSize=10;
sampleWidth=cubeSize; %Width
sampleThickness=cubeSize; %Thickness
sampleHeight=cubeSize; %Height
pointSpacings=2*ones(1,3); %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

%Define applied displacement
appliedStrain=0.4; %Linear strain (Only used to compute applied stretch)
loadingOption='compression'; % or 'tension'
switch loadingOption
    case 'compression'
        stretchLoad=1-appliedStrain; %The applied stretch for uniaxial loading
    case 'tension'
        stretchLoad=1+appliedStrain; %The applied stretch for uniaxial loading
end
displacementMagnitude=(stretchLoad*sampleHeight)-sampleHeight; %The displacement magnitude

%Material parameter set
E_youngs1=0.1; %Material Young's modulus
nu1=0.4; %Material Poisson's ratio

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

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
cubeDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
cubeElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(cubeDimensions,cubeElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E=meshStruct.elements; %The elements
V=meshStruct.nodes; %The nodes (vertices)
Fb=meshStruct.facesBoundary; %The boundary faces
Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E,1),1); %Element material indices

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node sets
logicFace=Cb==1; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_X=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==3; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_Y=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==5; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_Z=unique(Fr(:)); %Node set part of selected face

%Prescribed displacement nodes
logicPrescribe=Cb==6; %Logic for current face set
Fr=Fb(logicPrescribe,:); %The current face set
bcPrescribeList=unique(Fr(:)); %Node set part of selected face

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList_X,:),'r.','MarkerSize',markerSize);
hl(2)=plotV(V(bcSupportList_Y,:),'g.','MarkerSize',markerSize);
hl(3)=plotV(V(bcSupportList_Z,:),'b.','MarkerSize',markerSize);
hl(4)=plotV(V(bcPrescribeList,:),'k.','MarkerSize',markerSize);

legend(hl,{'BC x support','BC y support','BC z support','BC z prescribe'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

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

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

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

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

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

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='neo-Hookean';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.E=E_youngs1;
febio_spec.Material.material{1}.v=nu1;

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

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

% -> NodeSets
nodeSetName1='bcSupportList_X';
nodeSetName2='bcSupportList_Y';
nodeSetName3='bcSupportList_Z';
nodeSetName4='bcPrescribeList';

febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList_X(:);

febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcSupportList_Y(:);

febio_spec.Mesh.NodeSet{3}.ATTR.name=nodeSetName3;
febio_spec.Mesh.NodeSet{3}.node.ATTR.id=bcSupportList_Z(:);

febio_spec.Mesh.NodeSet{4}.ATTR.name=nodeSetName4;
febio_spec.Mesh.NodeSet{4}.node.ATTR.id=bcPrescribeList(:);

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

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x';

febio_spec.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{2}.dofs='y';

febio_spec.Boundary.bc{3}.ATTR.type='fix';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName3;
febio_spec.Boundary.bc{3}.dofs='z';

febio_spec.Boundary.bc{4}.ATTR.type='prescribe';
febio_spec.Boundary.bc{4}.ATTR.node_set=nodeSetName4;
febio_spec.Boundary.bc{4}.dof='z';
febio_spec.Boundary.bc{4}.scale.ATTR.lc=1;
febio_spec.Boundary.bc{4}.scale.VAL=displacementMagnitude;
febio_spec.Boundary.bc{4}.relative=0;

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

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

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

febio_spec.Output.logfile.element_data{2}.ATTR.file=febioLogFileName_stress_prin;
febio_spec.Output.logfile.element_data{2}.ATTR.data='s1;s2;s3';
febio_spec.Output.logfile.element_data{2}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{3}.ATTR.file=febioLogFileName_stiffness;
febio_spec.Output.logfile.element_data{3}.ATTR.data=['cxxxx;',...
                                                     'cxxyy;','cyyyy;',...
                                                     'cxxzz;','cyyzz;','czzzz;',...
                                                     'cxxxy;','cyyxy;','czzxy;','cxyxy;',...
                                                     'cxxyz;','cyyyz;','czzyz;','cxyyz;','cyzyz;',...
                                                     'cxxxz;','cyyxz;','czzxz;','cxyxz;','cyzxz;','cxzxz'];
febio_spec.Output.logfile.element_data{3}.ATTR.delim=',';

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

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

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

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

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

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

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

%% Import FEBio results

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

    %%

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

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

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

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

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

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1,2); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(cMap); colorbar;
    caxis([0 max(DN_magnitude)]); caxis manual;
    axis(axisLim(V_DEF)); %Set axis limits statically
    view(140,30);
    camlight headlight;

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

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

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress),1,1);

    %Access data
    E_stress_mat=dataStruct.data;

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

    [CV]=faceToVertexMeasure(E,V,E_stress_mat(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure  /usr/local/MATLAB/R2020a/bin/glnxa64/jcef_helper: symbol lookup error: /lib/x86_64-linux-gnu/libpango-1.0.so.0: undefined symbol: g_ptr_array_copy

    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$\sigma_{zz}$ [MPa]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',1,2); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(cMap); colorbar;
    caxis([min(E_stress_mat(:)) max(E_stress_mat(:))]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    view(140,30);
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments

        [CV]=faceToVertexMeasure(E,V,E_stress_mat(:,:,qt));

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

    %%
    % Calculate metrics to visualize stretch-stress curve

    DZ_set=N_disp_mat(bcPrescribeList,end,:); %Z displacements of the prescribed set
    DZ_set=mean(DZ_set,1); %Calculate mean Z displacements across nodes
    stretch_sim=(DZ_set(:)+sampleHeight)./sampleHeight; %Derive stretch
    stress_cauchy_sim=mean(squeeze(E_stress_mat(:,end,:)),1)';

    %%
    % Visualize stress-stretch curve

    cFigure; hold on;
    title('Uniaxial stress-stretch curve','FontSize',fontSize);
    xlabel('$\lambda$ [.]','FontSize',fontSize,'Interpreter','Latex');
    ylabel('$\sigma_{zz}$ [MPa]','FontSize',fontSize,'Interpreter','Latex');

    plot(stretch_sim(:),stress_cauchy_sim(:),'r-','lineWidth',lineWidth);

    view(2); axis tight;  grid on; axis square; box on;
    set(gca,'FontSize',fontSize);
    drawnow;

    %%
    % Importing element principal stresses from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress_prin),1,1);

    %Access data
    E_stress_prin_mat=dataStruct.data;

    %Compute pressure
    P = squeeze(-1/3*mean(sum(E_stress_prin_mat,2),1));

    %%
    % Visualize pressure-stretch curve

    cFigure; hold on;
    title('Pressure-stretch curve','FontSize',fontSize);
    xlabel('$\lambda$ [.]','FontSize',fontSize,'Interpreter','Latex');
    ylabel('$p$ [MPa]','FontSize',fontSize,'Interpreter','Latex');

    plot(stretch_sim(:),P(:),'r-','lineWidth',lineWidth);

    view(2); axis tight;  grid on; axis square; box on;
    set(gca,'FontSize',fontSize);
    drawnow;

    %% IMPORTING ELEMENT STIFFNESS MATRICES
    % Importing element stiffness tensors from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stiffness),0,1); %Nodal forces

    stiffness_mat=dataStruct.data(:,1:end,end); %Final stiffness state

    stiffness_mat_voigt=stiffness_mat(:,[1  2  4  11 16 7;...
                                         2  3  5  12 17 8;...
                                         4  5  6  13 18 9;...
                                         11 12 13 15 20 14;...
                                         16 17 18 20 21 19;...
                                         7  8  9  14 19 10]);
    stiffness_mat_voigt=reshape(stiffness_mat_voigt',6,6,size(stiffness_mat_voigt,1));
    stiffness_mat_voigt=reshape(mat2cell(stiffness_mat_voigt,6,6,...
        ones(size(stiffness_mat_voigt,3),1)),[size(stiffness_mat,1),1]);

    stiffness_mat_kelvin=stiffness_mat_voigt;
    for q=1:1:numel(stiffness_mat_voigt)
        cVoigt=stiffness_mat_voigt{q};
        c=voigtUnMap(cVoigt);
        cKelvin=kelvinMap(c);
        stiffness_mat_kelvin{q}=cKelvin;
    end

    %%
    % Visualize an element stiffness tensor

    viewFourthOrderTensor(c); %Visualize tensor C

end

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [C, S] = fcs(x)
%FCS Computes Fresnel integrals.
%[C, S] = fcs(x) returns Fresnel integrals C ans S for argument x,
% x must be double and real.
% F = fcs(x) returns complex F = C+i*S
%
% Algorithm:
% This function uses an improved method for computing Fresnel integrals
% with an error of less then 1x10-9, described in:
%
% Klaus D. Mielenz, Computation of Fresnel Integrals. II
% J. Res. Natl. Inst. Stand. Technol. 105, 589 (2000), pp 589-590
%
% Copyright (c) 2002 by Peter L. Volegov (volegov@unm.edu)
% All Rights Reserved.

error('You need to compile fcs.c to a mex file for your platform.');


text:
% Vm = vsh_modified_in_fast(theta,phi,l,m,p0)
%
% Vector spherical harmonic funtion modified from Hill's function
% V_lm for quantum numbers l and m. Angles theta and phi should be
% given in radians. p0 contains the values of the Legendre polynomial.
%
function Vm = vsh_modified_in_fast(theta,phi,l,m,p0)

% Copyright (c) 2016, Elekta Oy
% ---------------------------------------
%
% Redistribution and use of the Software in source and binary forms, with or without
% modification, are permitted for non-commercial use.
%
% The Software is provided "as is" without warranties of any kind, either express or
% implied including, without limitation, warranties that the Software is free of defects,
% merchantable, fit for a particular purpose. Developer/user agrees to bear the entire risk
% in connection with its use and distribution of any and all parts of the Software under this license.
%

scale = 1;

%scale_sph = (-1)^m*sqrt((2*l+1)*prod(1:(l-m))/(4*pi*prod(1:(l+m))));
%scale_minus = 1/((-1)^(m-1)*sqrt((2*l+1)*prod(1:(l-(m-1)))/(4*pi*prod(1:(l+(m-1))))));
%scale_plus = 1/((-1)^(m+1)*sqrt((2*l+1)*prod(1:(l-(m+1)))/(4*pi*prod(1:(l+(m+1))))));

Y = spharm_fast(theta,phi,l,m,p0);
if m > -l
  Yminus = spharm_fast(theta,phi,l,m-1,p0);
else
  Yminus = 0;
end
if m < l
  Yplus = spharm_fast(theta,phi,l,m+1,p0);
else
  Yplus = 0;
end

cosp = cos(phi); sinp = sin(phi);
%dY = 0.5*scale_sph*((l+m)*(l-m+1)*scale_minus*Yminus*complex(cos(phi),sin(phi))-scale_plus*Yplus*complex(cos(-phi),sin(-phi)));  % dY/dtheta
%dY = 0.5*(-sqrt((l+m)*(l-m+1))*Yminus*complex(cos(phi),sin(phi))+sqrt((l-m)*(l+m+1))*Yplus*complex(cos(-phi),sin(-phi)));
dY = 0.5*(-sqrt((l+m)*(l-m+1))*Yminus*complex(cosp,sinp)+sqrt((l-m)*(l+m+1))*Yplus*complex(cosp,-sinp));
if theta == 0 % To prevent division by zero
  theta = eps;
end
Vm(1) = scale*(-(l+1))*Y;
Vm(2) = scale*dY;
Vm(3) = scale*(i*m/sin(theta))*Y;

author: spm


text:
%% Simple rotation model images stitcher
% A basic example on image stitching.
%
% In this demo, we show how to use the high-level stitching API provided by
% |cv.Stitcher|, and we learn how to use preconfigured stitcher configurations
% to stitch images using different camera models.
%
% Sources:
%
% * <https://docs.opencv.org/3.3.1/d8/d19/tutorial_stitcher.html>
% * <https://github.com/opencv/opencv/blob/3.3.1/samples/cpp/stitching.cpp>
%

%% Input images (two or more)
imgs = {
    imread(fullfile(mexopencv.root(),'test','b1.jpg')), ...
    imread(fullfile(mexopencv.root(),'test','b2.jpg'))
};
for i=1:numel(imgs)
    subplot(1,numel(imgs),i), imshow(imgs{i})
end

%% Options

% Try to use GPU. The default value is false.
% All default values are for CPU mode.
try_use_gpu = false;

% Determines configuration of stitcher. The default is 'Panorama' mode
% suitable for creating photo panoramas. Option 'Scans' is suitable for
% stitching materials under affine transformation, such as scans.
smode = 'Panorama';

% Internally create three chunks of each image to increase stitching success
divide_images = false;

%%
if divide_images
    for i=1:numel(imgs)
        sz = size(imgs{i});
        imgs{i} = {
            cv.Rect.crop(imgs{i}, [0       0 sz(2)/2 sz(1)]), ...
            cv.Rect.crop(imgs{i}, [sz(2)/3 0 sz(2)/2 sz(1)]), ...
            cv.Rect.crop(imgs{i}, [sz(2)/2 0 sz(2)/2 sz(1)])
        };
    end
    imgs = [imgs{:}];
end

%% Stitch
stitcher = cv.Stitcher('Mode',smode, 'TryUseGPU',try_use_gpu);
tic
pano = stitcher.stitch(imgs);
toc

%% Panorama result
figure, imshow(pano)

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function info = optiSolverInfo(solver,probtype,nlprob,opts)
%OPTISOLVERINFO  Return Configuration Information about a particular solver
%
%   info = optiSolverInfo(solver,probtype)

%Default input args and checks
if(nargin < 4), opts = []; end
if(nargin < 3), nlprob = []; end
if(nargin < 2 || isempty(probtype)), probtype = ''; end
if(nargin < 1 || ~ischar(solver)), error('You must supply the name of the solver to this function'); end

%Create default information structure
info.der1 = 0;
info.der2 = 0;
info.global = 0;
info.sparse = 0;
info.multialg = 0;
info.parallel = 0;
info.whitebox = 0;
info.con = struct('con',0,'bnd',0,'lineq',0,'leq',0,'qineq',0,'qeq',0,'sdp',0,'nlineq',0,...
                  'nleq',0,'int',0);
info.opt = struct('miter',0,'meval',0,'mnode',0,'mtime',0,'tolr',0,...
                  'tola',0,'tolint',0,'disp',0,'opts',0,'ctrlc',0,...
                  'iterf',0);
info.unq = struct('x0',0,'qpform',0,'hform',0,'jstr',0,'hstr',0);

%Constants
YES = 1;
SOME = -1;

%Fill in information based on selected solver
switch(lower(solver))
    case 'baron'
        info.global = YES;
        info.sparse = YES;
        info.whitebox = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,yes,no,yes,yes,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.qineq = YES;
        info.con.nlineq = YES;
        info.con.nleq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(some,no,yes,yes,yes,no,no,yes,no,yes,no);
        info.opt.miter = SOME;
        info.opt.mnode = YES;
        info.opt.mtime = YES;
        info.opt.tola = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = SOME;

    case 'bonmin'
        info.der1 = YES;
        info.der2 = SOME;
        info.sparse = YES;
        info.multialg = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,some,some,no,no,yes,yes,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.nlineq = YES;
        info.con.nleq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(some,no,yes,yes,some,no,yes,yes,yes,yes,no);
        info.opt.miter = SOME;
        info.opt.mnode = YES;
        info.opt.mtime = YES;
        info.opt.tolr = SOME;
        info.opt.tolint = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        %unique
        info.unq.x0 = YES;
        info.unq.hform = 'tril';
        info.unq.jstr = YES;
        info.unq.hstr = YES;

    case 'clp'
        info.sparse = YES;
        info.multialg = YES;
        info.parallel = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,no,no,yes,no,yes,no);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp =  YES;
        info.opt.ctrlc = YES;
        info.opt.opts = YES;
        %unique
        info.unq.qpform = 'tril';

    case 'cbc'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(no,no,yes,yes,no,no,yes,yes,no,yes,no);
        info.opt.mnode = YES;
        info.opt.mtime = YES;
        info.opt.tolint = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;

    case 'cplex'
        info.sparse = YES;
        info.parallel = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,yes,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.qineq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,yes,yes,yes,no,yes,yes,yes,yes,no);
        info.opt.miter = YES;
        info.opt.mnode = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.tolint = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        %unique
        info.unq.x0 = SOME;
        info.unq.qpform = 'sym';

    case 'csdp'
        info.sparse = YES;
        info.parallel = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,some,no,yes,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = SOME;
        info.con.sdp = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,no,no,yes,yes,yes,no);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        %unique
        info.unq.x0 = SOME;

    case 'dsdp'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,some,no,yes,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = SOME;
        info.con.sdp = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,no,yes,no,yes,yes,yes,no);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tola = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        %unique
        info.unq.x0 = SOME;

    case 'filtersd'
        info.der1 = YES;
        info.sparse = YES;
        info.multialg = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,some,some,no,no,yes,yes,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = SOME;
        info.con.leq = SOME;
        info.con.nlineq = YES;
        info.con.nleq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,no,yes,no,some,no,yes,no,yes,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tola = SOME;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;
        info.unq.jstr = YES;

    case 'glpk'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,no,yes,yes,no,yes,no);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.tolint = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;

    case 'gmatlab'
        info.global = YES;
        info.sparse = YES;
        info.multialg = YES;
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.nlineq = YES;
        info.con.int = YES;
        switch(lower(probtype))
            case {'uno','nlp',''}
                %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,yes,yes,yes);
                info.con.leq = YES;
                info.con.nleq = YES;
        end
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,some,yes,yes,no,some,yes,yes,yes,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.mnode = SOME;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.tolint = SOME;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = SOME;

    case 'hybrj'
        info.der1 = SOME;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(no,yes,no,yes,no,no,no,yes,no,yes,yes);
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'ipopt'
        info.der1 = YES;
        info.der2 = SOME;
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,some,some,no,no,yes,yes,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.nlineq = YES;
        info.con.nleq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,no,no,yes,yes,yes,yes);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;
        info.unq.hform = 'tril';
        info.unq.jstr = YES;
        info.unq.hstr = YES;

    case 'lbfgsb'
        info.der1 = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,no,no,no,no,no,no,no);
        info.con.con = YES;
        info.con.bnd = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,no,no,yes,no,yes,yes);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'levmar'
        info.der1 = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,no);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,no,no,no,no,yes,no,some,yes);
        info.opt.miter = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = SOME;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'lipsol'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = SOME;
        info.con.leq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,no,no,yes,no,yes,no);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp =  YES;
        info.opt.ctrlc = YES;

    case 'lmder'
        info.der1 = SOME;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(no,yes,no,yes,yes,no,no,yes,no,yes,yes);
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'lp_solve'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(no,no,some,yes,no,no,yes,yes,no,no,no);
        info.opt.mnode = SOME;
        info.opt.mtime = YES;
        info.opt.tolint = YES;
        info.opt.disp = YES;

    case 'm1qn3'
        info.der1 = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,no,yes,no,some,no,yes,no,yes,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tola = SOME;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'matlab'
        info.der1 = SOME;
        info.der2 = SOME;
        info.sparse = YES;
        info.multialg = YES;
        info.con.con = YES;
        info.con.bnd = YES;
        switch(lower(probtype))
            case {'lp','qp','bilp'}
                %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
                info.con.lineq = YES;
                info.con.leq = YES;
                %unique
                info.unq.x0 = SOME;
            case {'nlp',''}
                info.der1 = YES;
                info.der2 = SOME;
                %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,yes,yes,yes);
                info.con.lineq = YES;
                info.con.leq = YES;
                info.con.nlineq = YES;
                info.con.nleq = YES;
                %unique
                info.unq.x0 = YES;
                info.unq.hform = 'sym';
                info.unq.jstr = YES;
                info.unq.hstr = YES;
        end
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,some,yes,yes,no,some,yes,yes,yes,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;

    case 'mkltrnls'
        %If we don't have options, assume not called by OPTI and don't allow internal mklJac
        if(isempty(opts))
           info.der1 = SOME;
        end
        %bnd = yes;
        info.con.con = YES;
        info.con.bnd = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,yes,yes,yes,no,yes,no,yes,yes);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.tola = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'mosek'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,yes,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.qineq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,yes,yes,yes,no,yes,yes,yes,yes,no);
        info.opt.miter = YES;
        info.opt.mnode = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.tolint = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        %unique
        info.unq.x0 = SOME;
        info.unq.qpform = 'tril';

    case 'mumps'
        %Nothing!

    case 'nl2sol'
        info.der1 = SOME;
        %bnd = yes;
        info.con.con = YES;
        info.con.bnd = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,no,no,yes,yes,no,yes,no,some,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.tolr = YES;
        info.opt.tola = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = SOME;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'nlopt'
        %Some NLOPT algorithms don't require derivatives, see if we have problem data
        if(~isempty(opts))
            if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts) && isfield(opts.solverOpts,'algorithm'))
                [~,~,~,~,~,info.der1] = nloptSolver(opts.solverOpts.algorithm);
            elseif(~isempty(nlprob) && isfield(nlprob,'algorithm'))
                [~,~,~,~,~,info.der1] = nloptSolver(nlprob.algorithm);
            end %else assume not required (default based on solver)
        else
            info.der1 = SOME; %generally some do
        end
        info.global = SOME;
        info.multialg = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,some,some,no,no,yes,yes,no);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = SOME;
        info.con.leq = SOME;
        info.con.nlineq = YES;
        info.con.nleq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(no,yes,no,yes,yes,yes,no,yes,yes,some,yes);
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.tola = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = SOME;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'nomad'
        info.global = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,some,no,no,no,yes,no,no);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = SOME;
        info.con.nlineq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,no,yes,yes,no,no,yes,yes,yes,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'ooqp'
        info.sparse = YES;
        info.multialg = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(no,no,no,no,no,no,no,yes,no,no,no);
        info.opt.miter = YES;
        info.opt.mtime = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        %unique
        info.unq.qpform = 'tril';

    case 'pswarm'
        info.global = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,no,no,no,no,no,no);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,yes,no,yes,yes,no,no,yes,yes,yes,yes);
        info.opt.miter = YES;
        info.opt.meval = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.opts = YES;
        info.opt.ctrlc = YES;
        info.opt.iterf = YES;
        %unique
        info.unq.x0 = YES;

    case 'qsopt'
        info.sparse = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,no,no,no,no,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(yes,no,no,some,no,no,no,some,no,no,no);
        info.opt.miter = YES;
        info.opt.mtime = SOME;
        info.opt.disp = SOME;

    case 'scip'
        info.global = YES;
        info.sparse = YES;
        info.whitebox = YES;
        %[bnd,lineq,leq,qineq,sdp,nlineq,nleq,sp] = fillFields(yes,yes,yes,yes,no,yes,yes,yes);
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = YES;
        info.con.qineq = YES;
        info.con.qeq = YES;
        info.con.nlineq = YES;
        info.con.nleq = YES;
        info.con.int = YES;
        %[miter,meval,mnode,mtime,tolr,tola,tolint,disp,opts,ctrlc,iterf] = fillCFields(some,no,yes,yes,yes,no,no,yes,no,yes,no);
        info.opt.miter = SOME;
        info.opt.mnode = YES;
        info.opt.mtime = YES;
        info.opt.tolr = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;

    case 'sedumi'
        info.sparse = YES;
        info.multialg = YES;
        %Constraints
        info.con.con = YES;
        info.con.bnd = YES;
        info.con.lineq = YES;
        info.con.leq = SOME;
        info.con.sdp = YES;
        %Options
        info.opt.miter = YES;
        info.opt.disp = YES;
        info.opt.ctrlc = YES;

    case 'auto'
        %Assume we require derivatives
        info.der1 = YES;

    otherwise
        error('Unknown solver: %s',solver);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%FtsZPolymerization process test case
%
% Author: Jayodita Sanghvi, jayodita@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/1/2010
classdef FtsZPolymerization_Test < edu.stanford.covert.cell.sim.ProcessTestCase
    methods
        function this = FtsZPolymerization_Test(name)
            this = this@edu.stanford.covert.cell.sim.ProcessTestCase(name);
        end

        function setUp(this)
            this.setUp@edu.stanford.covert.cell.sim.ProcessTestCase();
            warning('off', 'WholeCell:warning');
        end

        function loadTestFixture(this)
            this.loadTestFixture@edu.stanford.covert.cell.sim.ProcessTestCase();
            m = this.process;
            m.geometry.volume = 1e-18;
        end
    end

    %tests
    methods
        function testInitializedToSteadyState(this)
            m = this.process;

            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ) = 1000;
            m.boundEnzymes(:);
            this.setGTP(1e5);
            this.setGDP(1e4);

            m.initializeState();
            initialEnzymes = m.enzymes;

            for i = 1:5
                m.evolveState();
            end
            assertNonnegativeCounts(m.enzymes);
            assertEqual(countFtsZ(initialEnzymes), countFtsZ(m.enzymes));
            assertElementsAlmostEqual(initialEnzymes', m.enzymes', 'relative', 0.2, 3);
        end

        function testFtsZAndGtpBalance(this)
            m = this.process;
            this.setGTP(400);
            this.setGDP(0);
            m.enzymes(m.enzymeIndexs_FtsZ) = 300;
            m.enzymes(m.enzymeIndexs_FtsZ_GDP) = 100;
            m.enzymes(m.enzymeIndexs_FtsZ_activated) = [200 101 0 0 0 0 0 0 22]';

            for i = 1:20
                m.evolveState();
                assertNonnegativeCounts(m.enzymes);
                assertEqual(1000, countFtsZ(m.enzymes));
                assertEqual(...
                    m.substrates(m.substrateIndexs_gtp),...
                    sum(m.enzymes([m.enzymeIndexs_FtsZ m.enzymeIndexs_FtsZ_GDP])));
            end
        end

        function testNoFtsZ(this)
            m = this.process;
            m.enzymes(:) = 0;
            this.setGTP(1e5);

            m.evolveState();
            assertEqual(zeros(size(m.enzymes)), m.enzymes);
            assertEqual(1e5, m.substrates(m.substrateIndexs_gtp));
        end

        function testNoGtp(this)
            m = this.process;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ) = 1000;
            m.substrates(:) = 1e6;
            this.setGTP(0);
            this.setGDP(1e5);

            m.evolveState();

            assertTrue(all(m.substrates >= 0));
            assertTrue(all(m.enzymes >= 0));
            assertEqual(1000, [1 1 1:9] * m.enzymes);
            assertAllEqual(0, m.enzymes(m.enzymeIndexs_FtsZ_activated));
            assertEqual(0, m.substrates(m.substrateIndexs_gtp));
        end

        function testNoGdp(this)
            m = this.process;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ) = 1000;
            m.substrates(:) = 1e6;
            this.setGTP(1e5);
            this.setGDP(0);

            m.evolveState();

            assertTrue(all(m.substrates >= 0));
            assertTrue(all(m.enzymes >= 0));
            assertEqual(1000, [1 1 1:9] * m.enzymes);
            assertAllEqual(0, m.enzymes(m.enzymeIndexs_FtsZ_GDP));
            assertEqual(0, m.substrates(m.substrateIndexs_gdp));
        end

        function testNoWater(this)
            m = this.process;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ) = 1000;
            m.substrates(:) = 1e6;
            m.substrates(m.substrateIndexs_water) = 0;
            this.setGDP(0);

            m.evolveState();

            assertTrue(all(m.substrates >= 0));
            assertTrue(all(m.enzymes >= 0));
            assertEqual(1000, [1 1 1:9] * m.enzymes);
            assertAllEqual(0, m.enzymes(m.enzymeIndexs_FtsZ_GDP));
        end

        function testForwardActivationOnly(this)
            m = this.process;
            m.enzymes(:) = [1000 0 0 0 0 0 0 0 0 0 0]';  %inactive monomer only
            this.setGTP(1000);

            this.setAllReactionRatesToZero();
            m.activationFwd = 10;

            m.evolveState();
            assertEqual(0, m.substrates(m.substrateIndexs_gtp));
            assertEqual(0, m.enzymes(m.enzymeIndexs_FtsZ));
            assertEqual(0, m.enzymes(m.enzymeIndexs_FtsZ_GDP));
            assertEqual([1000 0 0 0 0 0 0 0 0], m.enzymes(m.enzymeIndexs_FtsZ_activated)');
        end

        function testReverseActivationOnly(this)
            m = this.process;
            m.enzymes(:) = [0 0 1000 0 0 0 0 0 0 0 0]';  %activated monomer only
            this.setGTP(0);

            this.setAllReactionRatesToZero();
            m.activationRev = 10;

            m.evolveState();
            assertEqual(1000, m.substrates(m.substrateIndexs_gtp));
            assertEqual(1000, m.enzymes(m.enzymeIndexs_FtsZ));
            assertEqual(0, m.enzymes(m.enzymeIndexs_FtsZ_GDP));
            assertEqual([0 0 0 0 0 0 0 0 0], m.enzymes(m.enzymeIndexs_FtsZ_activated)');
        end

        function testForwardExchangeOnly(this)
            m = this.process;
            m.enzymes(:) = [0 1000 0 0 0 0 0 0 0 0 0]';  %deactivated monomer only
            this.setGTP(1000);
            this.setGDP(0);

            this.setAllReactionRatesToZero();
            m.exchangeFwd = 1e6;

            m.evolveState();
            assertEqual(0, m.substrates(m.substrateIndexs_gtp));
            assertEqual(1000, m.substrates(m.substrateIndexs_gdp));
            assertEqual(0, m.enzymes(m.enzymeIndexs_FtsZ));
            assertEqual(0, m.enzymes(m.enzymeIndexs_FtsZ_GDP));
            assertEqual([1000 0 0 0 0 0 0 0 0], m.enzymes(m.enzymeIndexs_FtsZ_activated)');
        end

        function testReverseExchangeOnly(this)
            m = this.process;
            m.enzymes(:) = [0 0 1000 0 0 0 0 0 0 0 0]';  %activated monomer only
            this.setGTP(0);
            this.setGDP(1000);

            this.setAllReactionRatesToZero();
            m.exchangeRev = 1e6;

            m.evolveState();
            assertEqual(1000, m.substrates(m.substrateIndexs_gtp));
            assertEqual(0, m.substrates(m.substrateIndexs_gdp));
            assertEqual(0, m.enzymes(m.enzymeIndexs_FtsZ));
            assertEqual(1000, m.enzymes(m.enzymeIndexs_FtsZ_GDP));
            assertEqual([0 0 0 0 0 0 0 0 0], m.enzymes(m.enzymeIndexs_FtsZ_activated)');
        end

        function testForwardNucleationOnly(this)
            m = this.process;
            m.enzymes = [0 0 1000 0 0 0 0 0 0 0 0]';  %activated monomer only

            this.setAllReactionRatesToZero();
            m.nucleationFwd = 1e6;

            m.evolveState();
            assertElementsAlmostEqual(...
                [0 0 0 500 0 0 0 0 0 0 0], m.enzymes', 'absolute', 2);
            assertEqual(1000, countFtsZ(m.enzymes));
        end

        function testReverseNucleationOnly(this)
            m = this.process;
            m.enzymes = [0 0 0 500 0 0 0 0 0 0 0]';  %activated dimer only

            this.setAllReactionRatesToZero();
            m.nucleationRev = 10;

            m.evolveState();
            assertEqual([0 0 1000 0 0 0 0 0 0 0 0], m.enzymes');
        end

        function testForwardElongationOnly(this)
            m = this.process;
            m.enzymes = [0 0 700 100 0 0 0 0 0 0 0]';  %activated monomers and dimers

            this.setAllReactionRatesToZero();
            m.elongationFwd = 1e7;

            m.evolveState();
            assertElementsAlmostEqual(...
                [0 0 0 0 0 0 0 0 0 0 100], m.enzymes', 'absolute', 9);
            assertEqual(900, countFtsZ(m.enzymes));
        end

        function testReverseElongationOnly(this)
            m = this.process;
            m.substrates(:) = 1e6;
            m.enzymes = [0 0 0 0 0 0 0 0 0 0 100]';  %9mers only

            this.setAllReactionRatesToZero();
            m.elongationRev = 100;

            m.evolveState();
            assertElementsAlmostEqual(...
                [0 0 700 100 0 0 0 0 0 0 0], m.enzymes', 'absolute', 2);
            assertEqual(900, countFtsZ(m.enzymes));
        end

        function testApplySubstrateLimits(this)
            m = this.process;

            %% Case 1: ODEs trying to make FtsZ-GTP polymers from FtsZ-GDP, but unable to because no GTP
            %state before running ODEs
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ_GDP) = sum(1:9);
            m.substrates(:) = 1e6;
            m.substrates(m.substrateIndexs_gtp) = 0;

            %state after running ODEs
            substrates = m.substrates;
            enzymes = zeros(size(m.enzymes));
            enzymes(m.enzymeIndexs_FtsZ_activated) = 1;

            %final state
            [finEnzymes, finSubstrates] = m.applySubstrateLimits(enzymes, substrates);

            %assertions
            assertTrue(all(finSubstrates >= 0)); %counts are non-negative
            assertTrue(all(finEnzymes >= 0)); %counts are non-negative
            assertEqual([1 1 1:9] * m.enzymes, [1 1 1:9] * finEnzymes); %number of FtsZ monomer is constant
            assertEqual([0 1 1:9] * m.enzymes + [1 1] * m.substrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp]), ... %number of guanosine (free + bound) is constant
                [0 1 1:9] * finEnzymes + [1 1] * finSubstrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp]));
            assertIn([0 0 1:9] * finEnzymes + [0 1] * finSubstrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp]), ... %number of GTP molecules (free + bound) doesn't increase / energy is not created
                [0 [0 0 1:9] * m.enzymes + [0 1] * m.substrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp])]);
            assertEqual(m.substrates' * m.substrateMolecularWeights + m.enzymes' * m.enzymeMolecularWeights, ... %mass is conserved
                finSubstrates' * m.substrateMolecularWeights + finEnzymes' * m.enzymeMolecularWeights);

            %% Case 2: Making FtsZ-GDP from free FtsZ and GTP
            %state before running ODEs
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ) = 100;
            m.substrates(:) = 1e6;
            m.substrates(m.substrateIndexs_gdp) = 0;

            %state after running ODEs
            substrates = m.substrates;
            enzymes = zeros(size(m.enzymes));
            enzymes(m.enzymeIndexs_FtsZ_GDP) = 100;

            %final state
            [finEnzymes, finSubstrates] = m.applySubstrateLimits(enzymes, substrates);

            %assertions
            assertTrue(all(finSubstrates >= 0)); %counts are non-negative
            assertTrue(all(finEnzymes >= 0)); %counts are non-negative
            assertEqual([1 1 1:9] * m.enzymes, [1 1 1:9] * finEnzymes); %number of FtsZ monomer is constant
            assertEqual([0 1 1:9] * m.enzymes + [1 1] * m.substrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp]), ... %number of guanosine (free + bound) is constant
                [0 1 1:9] * finEnzymes + [1 1] * finSubstrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp]));
            assertIn([0 0 1:9] * finEnzymes + [0 1] * finSubstrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp]), ... %number of GTP molecules (free + bound) doesn't increase / energy is not created
                [0 [0 0 1:9] * m.enzymes + [0 1] * m.substrates([m.substrateIndexs_gdp; m.substrateIndexs_gtp])]);
            assertElementsAlmostEqual(m.substrates' * m.substrateMolecularWeights + m.enzymes' * m.enzymeMolecularWeights, ... %mass is conserved
                finSubstrates' * m.substrateMolecularWeights + finEnzymes' * m.enzymeMolecularWeights, ...
                'relative', 1e-8, 0);
        end

        function testGeneEssentiality(this)
            m = this.process;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_FtsZ) = 100;
            m.boundEnzymes(:) = 0;
            this.setGTP(1e5);
            this.setGDP(1e4);
            this.helpTestGeneEssentiality({
                'MG_224'},...  %ftsZ
                @(m,i) m.enzymes(m.enzymeIndexs_FtsZ_9mer) > ...
                       i.enzymes(m.enzymeIndexs_FtsZ_9mer), ...
                       struct('lengthSec', 20));
        end
    end

    %helper methods
    methods
        function setAllReactionRatesToZero(this)
            m = this.process;
            m.activationFwd = 0;
            m.activationRev = 0;
            m.exchangeFwd = 0;
            m.exchangeRev = 0;
            m.nucleationFwd = 0;
            m.nucleationRev = 0;
            m.elongationFwd = 0;
            m.elongationRev = 0;
        end

        function setGTP(this, numMolecules)
            m = this.process;
            m.substrates(m.substrateIndexs_gtp) = numMolecules;
        end

        function setGDP(this, numMolecules)
            m = this.process;
            m.substrates(m.substrateIndexs_gdp) = numMolecules;
        end
    end
end

function assertNonnegativeCounts(values)
    assertTrue(all(values >= 0));
    assertEqual(round(values), values);
end

function result = countFtsZ(counts)
    result = [1 1 1:length(counts)-2] * counts;
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [values, points] = sample(f, varargin)
%SAMPLE   Sample a DELTAFUN on an "appropriate" grid.
%   VALUES = SAMPLE(F, N) returns a vector VALUES of the values of F at N
%   "appropriately-chosen" points.  What "appropriately-chosen" means depends
%   on the type of SMOOTHFUN object on which F is ultimately based.
%
%   [VALUES, POINTS] = SAMPLE(F, N) returns also the vector POINTS at which the
%   values were computed.
%
%   [...] = SAMPLE(F) uses N = LENGTH(F).

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

    [varargout{1:nargout}] = sample(f.funPart, varargin{:});
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function stats = plot_obs_spa_for_ranks( fig_ind, saveResults )
% plot obs % vs corruption % for a few ranks

load ..\..\data\syn-lr-tnames   % load tnames

marker = { 'r-*', 'b-*', 'k-*', 'g-*' };
Ntens = length(tnames);
Itens = [ 1,2,3 ];
Ntens = length(Itens);
Nobs = 20;
Nnoise = 9;
Nrep = 1;
Imag = 1;
alg = 2;
rRatio = 1;
lambdaS = 1;
IsTC = true;
verbose = false;
err_bar = 0.01;

stats.obs = [];
stats.spa = zeros( Ntens, Nobs );
stats.exp = 'obs-spa-for-ranks';

figure( fig_ind );

for t = 1:Ntens
    tname = tnames{Itens(t)};
    load( ['..\..\data\',tname] );
    maxNoise = Nnoise;

    for i = Nobs:-1:1

        for j = maxNoise:-1:1
            err = zeros( 1, Nrep );
            for k = 1:Nrep
                dtemp = gen_syn_data( data, i, j, Imag, k );
                results = test_trpca( dtemp, alg, rRatio, lambdaS, IsTC, verbose );
                err(k) = results.rel_err1;
            end
            avg_err = mean(err);
            fprintf( '%s,  obs: %1.2f,  spa: %1.2f,  err: %1.2e\n', tname, data.obsPct(i), data.noisePct(j), avg_err );
            if avg_err <= err_bar
                % record spa
                stats.spa( t, i ) = data.noisePct(j);
                break;  % found the largest corruption % tolerated
            end
        end

        maxNoise = j;
        if j == 1 && stats.spa(t,i) == 0    % no noise tolerated from now on
            break;
        end
    end

    stats.obs = data.obsPct(1:Nobs);
    plot( stats.obs, stats.spa(t,:), marker{t}, 'LineWidth', 2 );
    hold on;
end

title( 'obs % v.s. max corruption % for various ranks' );
xlabel( 'obs %' );
ylabel( 'max corruption %' );
legend( tnames(1:Ntens) );
hold off;

if exist( 'saveResults', 'var' ) && ~isempty(saveResults) && saveResults
    save ..\..\data\syn-obs-spa stats
end

end

author: andrewssobral

repo: lrslibrary


text:
function [crfw, priors] = crfPairwiseTrain(imsegs, pvSP, phSP, adjlist, pE, params, edgelen)

global itercount;
itercount = 0;

vlab = {imsegs(:).vert_labels};
hlab = {imsegs(:).horz_labels};
spweight = {imsegs(:).npixels};
for f = 1:numel(imsegs)
    spweight{f} = spweight{f} / sum(spweight{f});
end

% get f1(y1, y2) = log(P(y1|x)) + log(P(y2|x))
% and f2(y1, y2) = I(y1=y2)*log(P(y1=y2|x)) + I(y1~=y2)*log(P(y1~=y2|x))

npair = 0;
for f = 1:numel(adjlist)
    npair = npair + size(adjlist{f}, 1);
end

f1 = zeros(npair, 49);
f2 = zeros(npair, 49);
f3 = zeros(npair, 49);
lab = zeros(npair, 1);
w = zeros(npair, 1);
%nbpair = zeros(npair, 2);

c = 0;
for f = 1:numel(adjlist)
    npairf = size(adjlist{f}, 1);
    pg = [pvSP{f}(:, 1)  repmat(pvSP{f}(:, 2), 1, 5).*phSP{f}  pvSP{f}(:, 3)];
    labf = (vlab{f}==1)*1 + (vlab{f}==2).*(hlab{f}>0).*(hlab{f}+1) + (vlab{f}==3)*7;

    s1 = adjlist{f}(:, 1);
    s2 = adjlist{f}(:, 2);

    edgeperc = zeros(numel(s1), 2);
    nb = zeros(size(pg, 1), 1);
    for k = 1:npairf
        nb(s1(k)) = nb(s1(k))+1;
        nb(s2(k)) = nb(s2(k))+1;
        edgeperc(k, 1) = edgelen{f}(s1(k), s2(k)) / (sum(edgelen{f}(s1(k), :))+sum(edgelen{f}(:, s1(k))));
        edgeperc(k, 2) = edgelen{f}(s1(k), s2(k)) / (sum(edgelen{f}(s2(k), :))+sum(edgelen{f}(:, s2(k))));
        w(c+k) = sum(spweight{f}([s1(k) s2(k)]));
    end
    %nbpair(c+1:c+npairf, :) = [nb(s1) nb(s2)];
%    min(1./edgeperc(:, 1) + 1./edgeperc(:, 2))
%    max(1./edgeperc(:, 1) + 1./edgeperc(:, 2))
    % unaries normalize for number of edges, and pairwise normalize for
    % edge length
    for k1 = 1:7
        for k2 = 1:7
            f1(c+1:c+npairf, (k1-1)*7+k2) = log(pg(s1, k1))./nb(s1) + log(pg(s2, k2))./nb(s2);
            kp = params(k1 + (k2-1)*7, :);
            if k1==k2
                f2(c+1:c+npairf, (k1-1)*7+k2) = ...
                    log(kp(3)./(1+exp(-kp(1)-kp(2)*(log(pE{f})-log(1-pE{f}))))).*...
                    (edgeperc(:, 1) + edgeperc(:, 2));
            else
                f3(c+1:c+npairf, (k1-1)*7+k2) = ...
                    log(kp(3)./(1+exp(-kp(1)-kp(2)*(log(pE{f})-log(1-pE{f}))))).*...
                    (edgeperc(:, 1) + edgeperc(:, 2));
            end
        end
    end

    lab(c+1:c+npairf) = (labf(s1)>0 & labf(s2)>0).*((labf(s1)-1)*7 + labf(s2));

    c = c + npairf;
end

ind = find(lab==0);
lab(ind) = [];
f1(ind, :) = [];
f2(ind, :) = [];
f3(ind, :) = [];
w(ind, :) = [];

w = w/sum(w);
%nbpair(ind, :) = [];

f2 = f2 - repmat(log(sum(exp(f2), 2)), [1 size(f2, 2)]);
f3 = f3 - repmat(log(sum(exp(f3), 2)), [1 size(f3, 2)]);

fnc{1} = f1;
fnc{2} = f2;
fnc{3} = f3;
%disp(numel(lab))

disp('uniform w')
w = ones(size(w))/numel(w);

[p, err] = avep([1 0 0 zeros(1,27)], fnc, lab, w);
disp(num2str([p err]))

initval = [3 0 1 3 -1 -1 -1 0 0 -1 2 -1 -1 -1 -1 2 2 -1 -1 -1 2 2 -2 -1 2 1 -1 2 2 1];
crfw = fminunc(@(x) objective(x, fnc, lab, w), initval, optimset('TolFun', 0.01));

[p, err] = avep(crfw, fnc, lab, w);
disp(num2str([p err]))

priors = priormat(7, crfw(numel(fnc)+1:end));
crfw = crfw(1:numel(fnc));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function lpw = objective(w, f, lab, weight)
global itercount;
correctind = [1:numel(lab)]' + (lab-1)*numel(lab);

%w(1:3) = max(w(1:3), 0);

pmat = priormat(7, w(numel(f)+1:end));
lpw = pmat(lab);
tmpz = repmat(pmat', [numel(lab) 1]);
for k = 1:numel(f)
    lpw = lpw + w(k)*f{k}(correctind);
    tmpz = tmpz + w(k)*f{k};
end
z = log(sum(exp(tmpz), 2));

lpw = -sum((lpw - z).*weight);

itercount = itercount + 1;
if mod(itercount, 25)==0
    disp(num2str([lpw w pmat(end)]))
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [p err] = avep(w, f, lab, weight)

correctind = [1:numel(lab)]' + (lab-1)*numel(lab);

%w(1:3) = max(w(1:3), 0);

pmat = priormat(7, w(numel(f)+1:end));
lpw = pmat(lab);
tmpz = repmat(pmat', [numel(lab) 1]);
for k = 1:numel(f)
    lpw = lpw + w(k)*f{k}(correctind);
    tmpz = tmpz + w(k)*f{k};
end
z = log(sum(exp(tmpz), 2));

p = sum(exp(lpw-z).*weight);

[tmp, guess] = max(tmpz, [], 2);
err = sum((guess~=lab).*weight);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pmat = priormat(nlab, param);

param(end+1) = -sum(param); % sum of free params equals zero
pmat = zeros(nlab, nlab);
c=0;
for k1 = 1:nlab
    for k2=k1:nlab
        c = c + 1;
        pmat(k1, k2) = param(c);
        pmat(k2, k1) = param(c);
    end
end
pmat = pmat(:);
pmat(end) = -sum(pmat(1:end-1));

author: Cloud-CV

repo: object-proposals