text:
function varargout = BrowserGUI(varargin)
%BrowserGUI GUI code for OlyDat Browser
%
%   See also OlyDat.Browser.

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
if ismac
    gui_Name = 'BrowserGUI_Mac';
else
    gui_Name = 'BrowserGUI_PC'; % use this for both PC and Unix
end
gui_State = struct('gui_Name',       gui_Name, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @browserGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @browserGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% browserGUI = BrowserGUI(browserObj)
function browserGUI_OpeningFcn(hObject, eventdata, handles, varargin) %#ok<*INUSL>
handles.output = hObject;
handles.browser = varargin{1};
guidata(hObject, handles);
centerGUI(hObject);

function varargout = browserGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

function pumAuxVar_Callback(hObject, eventdata, handles) %#ok<*DEFNU>
handles.browser.auxVarChanged();

function lbBStats_Callback(hObject, eventdata, handles)
handles.browser.statChanged();

function pumBStats_Callback(hObject, eventdata, handles)
handles.browser.statChanged();

function lbGroupVars_Callback(hObject, eventdata, handles)
handles.browser.groupOrGroup2SelectionChanged();

function pumGroup_Callback(hObject, eventdata, handles)
handles.browser.groupChanged();

function lbGroupVars2_Callback(hObject, eventdata, handles)
handles.browser.groupOrGroup2SelectionChanged();

function pumGroup2_Callback(hObject, eventdata, handles)
handles.browser.group2Changed();

function pbUseAllGroupVals_Callback(hObject, eventdata, handles)
strs = get(handles.lbGroupVars,'String');
set(handles.lbGroupVars,'Value',1:numel(strs));
handles.browser.groupOrGroup2SelectionChanged();

function pbUseAllGroupVals2_Callback(hObject, eventdata, handles)
strs = get(handles.lbGroupVars2,'String');
set(handles.lbGroupVars2,'Value',1:numel(strs));
handles.browser.groupOrGroup2SelectionChanged();

function pbPlot_Callback(hObject, eventdata, handles)
handles.browser.plot();

function pumPlotType_Callback(hObject, eventdata, handles)
handles.browser.newPlot();

function cbRemoveExp_Callback(hObject, eventdata, handles)
handles.browser.setRemoveExp(get(hObject,'Value'));

function pbCompare_Callback(hObject, eventdata, handles)
handles.browser.addComparePlot();

function figure1_CloseRequestFcn(hObject, eventdata, handles)
handles.browser.browserGUIWantsToClose();

function pbEnterCurationInfo_Callback(hObject, eventdata, handles)
handles.browser.curationOpenCurationInfoEntry();

function pbNextExpDetail_Callback(hObject, eventdata, handles)
handles.browser.nextExpDetail;

function pbPrevExpDetail_Callback(hObject, eventdata, handles)
handles.browser.prevExpDetail;

function pbCreateCustomGroup_Callback(hObject, eventdata, handles)
grpname = inputdlg('Specify a name for the group');
if ~isempty(grpname)
    grpname = grpname{1};
    handles.browser.createCustomGroup(grpname);
end

function pbLineReport_Callback(hObject, eventdata, handles)
handles.browser.openLineReport();

function pbExpDetail_Callback(hObject, eventdata, handles)
handles.browser.openExperimentDetailPlots();

function pbFileSystem_Callback(hObject, eventdata, handles)
handles.browser.openExperimentFileSystem();

% MENU
function mnuCuration_Load_Callback(hObject, eventdata, handles)
handles.browser.curationLoadCurationFile();

function mnuCuration_Append_Callback(hObject, eventdata, handles)
handles.browser.curationAppendCurationDataToFile();

function mnuCuration_Write_Callback(hObject, eventdata, handles)
handles.browser.curationWriteCurationDataToNewFile();

function mnuCuration_OpenCurationManager_Callback(hObject, eventdata, handles)
handles.browser.curationOpenCurationManager();

function pbExport_Callback(hObject, eventdata, handles)
handles.browser.exportPlot();

author: kristinbranson

repo: JAABA


text:
function update_spatial_batch(obj, use_parallel)
%% update the the spatial components for all batches and then summarize the results
% input:
%   use_parallel: boolean, do initialization in patch mode or not.
%       default(true); we recommend you to set it false only when you want to debug the code.

%% Author: Pengcheng Zhou, Columbia University, 2017
%% email: zhoupc1988@gmail.com

%% update spatial components for all batches
nbatches = length(obj.batches);

[d, K] = size(obj.A);
if K==0
    fprintf('You have to initialize neurons first\n');
    return;
end
A = zeros(d, K);
cc = zeros(nbatches, K);

for mbatch=1:nbatches
    batch_k = obj.batches{mbatch};
    neuron_k = batch_k.neuron;

    fprintf('\nprocessing batch %d/%d\n', mbatch, nbatches);

    % update temporal components
    neuron_k.update_spatial_parallel(use_parallel);
    cc(mbatch, :) = sum(neuron_k.C.^2, 2);
    A = A + bsxfun(@times, neuron_k.A, cc(mbatch, :));
end

obj.A = bsxfun(@times, A, 1./sum(cc, 1));
obj.post_process_spatial();

%% spread the same A to all batches and update temporal components
for mbatch=1:(length(obj.batches)-1)
    batch_k = obj.batches{mbatch};
    neuron_k = batch_k.neuron;

    neuron_k.A = obj.A;

    batch_k.neuron = neuron_k;
    obj.batches{mbatch} = batch_k;
end

author: zhoupc


text:
function table = hen_power_product ( p, e )

%*****************************************************************************80
%
%% HEN_POWER_PRODUCT: power products, x^e*Hen(i,x)*Hen(j,x).
%
%  Discussion:
%
%    Hen(i,x) is the normalized probabilist's Hermite polynomial of degree I.
%
%    For polynomial chaos applications, it is of interest to know the
%    value of the integrals of products of X with every possible pair
%    of basis functions.  That is, we'd like to form
%
%      Tij = Integral ( -oo < X < +oo ) X^E * Hen(I,X) * Hen(J,X) exp(-0.5*X*X) dx
%
%    We will estimate these integrals using Gauss-Hermite quadrature.
%
%    When E is 0, the computed table should be the identity matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polyonomial factors.
%    0 <= P.
%
%    Input, integer E, the exponent of X in the integrand.
%    0 <= E.
%
%    Output, real TABLE(P+1,P+1), the table of integrals.  TABLE(I,J)
%    represents the weighted integral of X^E * Hen(I+1,X) * Hen(J+1,X).
%
  table(1:p+1,1:p+1) = 0.0;

  order = p + 1 + floor ( ( e + 1 ) / 2 );

  [ x_table, w_table ] = he_quadrature_rule ( order );

  for k = 1 : order

    x = x_table(k);
    h_table = hen_polynomial_value ( 1, p, x );
%
%  The following formula is an outer product in H_TABLE.
%
    if ( e == 0 )
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) * h_table(1:p+1)' * h_table(1:p+1);
    else
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) * x^e * h_table(1:p+1)' * h_table(1:p+1);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function tests = test_imEuler1dEstimate(varargin)
%TEST_IMEULER1DESTIMATE  One-line description here, please.
%
%   output = test_imEuler1dEstimate(input)
%
%   Example
%   test_imEuler1dEstimate
%
%   See also
%

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

tests = functiontests(localfunctions);

function testHoriz(testCase)

img = logical([0 1 1 1 0 0]);
exp = 1;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([0 1 1 0 1 0]);
exp = 2;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([0 0 1 1 1 1]);
exp = .5;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([1 1 1 0]);
exp = .5;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([1 1 1 1 1 1]);
exp = 0;
assertEqual(testCase, exp, imEuler1dEstimate(img));

function testVertical(testCase)

img = logical([0 1 1 1 0 0])';
exp = 1;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([0 1 1 0 1 0]');
exp = 2;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([0 0 1 1 1 1]');
exp = .5;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([1 1 1 0]');
exp = .5;
assertEqual(testCase, exp, imEuler1dEstimate(img));

img = logical([1 1 1 1 1 1]');
exp = 0;
assertEqual(testCase, exp, imEuler1dEstimate(img));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% column sum
% function Z=csum(X)

function Z=csum(X)

N=length(X(:,1));
if (N>1)
  Z=sum(X);
else
  Z=X;
end;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function z = power( x, y )

%   Disciplined convex programming information for POWER (.^):
%      When used in CVX expressions, either X or Y must be constant. Only
%      certain convex or concave branches are accepted as valid:
%         --- if both X and Y are constant, then Z = X.^Y is interpreted
%             precisely as the MATLAB built-in version.
%         --- if Y is constant and 0 < Y < 1, then Z = X.^Y is concave and
%             nondecreasing in X. Therefore, X must be concave, and is
%             implicitly constrained to be nonnegative.
%         --- if Y is constant and Y == 1, then Z = X.
%         --- if Y is constant and a positive even integer, then Z = X.^Y
%             is convex and nonmonotonic in X. Therefore, X must be affine.
%         --- if Y is constant and Y > 1, but *not* an integer, then
%             Z = X.^Y is convex and nonmonotonic in X. Therefore, X must
%             be affine (and real), and is implicitly constrained to be
%             nonnegative.
%      In expert mode, additional cases are handled:
%         --- if X is constant and 0 < X < 1, then Z = X.^Y is convex and
%             nonincreasing in X. Therefore, Y must be concave.
%         --- if X is constant and X == 1, then Z = 1.
%         --- if X is constant and X > 1, then Z = X.^Y is convex and
%             nondecreasing in X. Therefore, Y must be convex.
%      All other combinations are rejected as invalid. For instance, if Y
%      is an odd integer, then X .^ Y is neither convex nor concave, so it
%      is rejected. In such cases, consider using POW_P, POW_POS, or
%      POW_ABS instead.
%
%   Disciplined geometric programming information for POWER (.^):
%      In disciplined geometric programs, the power operation Z=X.^Y is
%      valid only if Y is a real constant. There are no restrictions on
%      X. Note that a negative exponent Y reverses curvature; that is, Z
%      is log-convex if X is log-concave, and vice versa.

%
% Check sizes
%

sx = size( x ); xs = all( sx == 1 );
sy = size( y ); ys = all( sy == 1 );
if xs,
    sz = sy;
elseif ys || isequal( sx, sy ),
    sz = sx;
else
    error( 'Matrix dimensions must agree.' );
end

%
% Determine the expression types
%

if cvx_isconstant( y ),

    z = pow_cvx( x, y, 'power' );
    return

elseif ~cvx_isconstant( x ),

    error( 'Disciplined convex programming error:\n   In an expression X .^ Y, either X or Y must be constant.', 1 ); %#ok

end

%
% Now handle constant .^ non-constant
%

persistent remap
if isempty( remap ),
	remap_y1 = cvx_remap( 'real-affine' );
	remap_y2 = cvx_remap( 'convex' )    & ~remap_y1;
	remap_y3 = cvx_remap( 'concave' )   & ~remap_y1;
	remap_y4 = cvx_remap( 'log-valid' ) & ~( remap_y1 | remap_y2 | remap_y3 );
	remap    = [0;0;2;2;2] * remap_y1 + ...
	           [0;0;0;2;2] * remap_y2 + ...
	           [0;0;2;2;0] * remap_y3 + ...
	           [0;1;0;2;0] * remap_y4;
end
x  = cvx_constant( x );
vy = cvx_classify( y );
vx = 1 + isreal( x ) .* ( ( x >= 0 ) + ( x > 0 ) + ( x >= 1 ) + ( x > 1 ) );
vr = remap( vx + size( remap, 1 ) * ( vy - 1 ) );
vu = sort( vr(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );

%
% Perform the individual computations and combine
%

x = cvx( x ); xt = x;
y = cvx( y ); yt = y;
if nv ~= 1,
    z = cvx( sz, [] );
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    if nv ~= 1,
        t = vr == vu( k );
        if ~xs, xt = cvx_subsref( x, t ); sz = size( xt ); end
        if ~ys, yt = cvx_subsref( y, t ); sz = size( yt ); end
    end

    %
    % The computational kernels
    %

    switch vu( k ),
        case 0,
            % Invalid
            error( 'Disciplined convex programming error:\n    Cannot perform the operation {%s}.^{%s}', cvx_class( xt, true, true, true ), cvx_class( yt, true, true, true ) );
        case 1,
            % zero .^ convex
            cvx_optval = cvx( zeros( sz ) );
        case 2,
            % (0<x<1) .^ concave, (x>1) .^ convex
            cvx_optval = exp( log( cvx_constant( xt ) ) .* yt );
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    %
    % Store the results
    %

    if nv == 1,
        z = cvx_optval;
    else
        z = cvx_subsasgn( z, t, cvx_optval );
    end

end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [mask, weight]=mesh2mask(node,face,xi,yi,hf)
%
% [mask weight]=mesh2mask(node,face,Nxy)
%   or
% [mask weight]=mesh2mask(node,face,[Nx,Ny])
%   or
% [mask weight]=mesh2mask(node,face,xi,yi,hf)
%
% fast rasterization of a 2D mesh to an image with triangle index labels
%
% author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
% date for initial version: July 18,2013
%
% input:
%      node: node coordinates, dimension N by 2 or N by 3 array
%      face: a triangle surface, N by 3 or N by 4 array
%      Nx,Ny,Nxy: output image in x/y dimensions, or both
%      xi,yi: linear vectors for the output pixel center positions in x/y
%      hf: (optional) the handle of a pre-created figure window, for faster
%          rendering
%
% output:
%      mask: a 2D image, the value of each pixel is the index of the
%            enclosing triangle, if the pixel is outside of the mesh, NaN
%      weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions for
%            the mask
%
% note: This function only works in MATLAB when the DISPLAY is not
%       disabled. The maximum size of the mask output is limited by the
%       screen size.
%
% example:
%
%   [no,fc]=meshgrid6(0:5,0:5);
%   [mask weight]=mesh2mask(no,fc,-1:0.1:5,0:0.1:5);
%   imagesc(mask);
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin==3 && length(xi)==1 && xi>0)
    mn=min(node);
    mx=max(node);
    df=(mx(1:2)-mn(1:2))/xi;
elseif(nargin==3 && length(xi)==2 && all(xi>0))
    mn=min(node);
    mx=max(node);
    df=(mx(1:2)-mn(1:2))./xi;
elseif(nargin==4 || nargin==5)
    mx=[max(xi) max(yi)];
    mn=[min(xi) min(yi)];
    df=[min(diff(xi(:))) min(diff(yi(:)))];
else
    error('you must give at least xi input');
end
if(size(node,2)<=1 || size(face,2)<=2)
    error('node must have 2 or 3 columns; face can not have less than 2 columns');
end

if(nargin<5)
    hf=figure('visible','on');
else
    clf(hf);
end
patch('Vertices',node,'Faces',face,'linestyle','none','FaceColor','flat',...
 'FaceVertexCData',(1:size(face,1))','CDataMapping', 'scaled');
set(gca, 'Position', [0 0 1 1]);
cm=jet(size(face,1));
colormap(cm);
axis off
set(gca,'xlim',[mn(1) mx(1)]);
set(gca,'ylim',[mn(2) mx(2)]);
set(gca,'clim',[1 size(face,1)]);

output_size = round((mx(1:2)-mn(1:2))./df); %Size in pixels

if(isoctavemesh || isempty(getenv('DISPLAY')))
    resolution = 300; %Resolution in DPI
    set(gcf,'PaperPositionMode','manual')
    set(gcf,'paperunits','inches','paperposition',[0 0 output_size/resolution]);
    deletemeshfile(mwpath('post_mesh2mask.png'));
    print(mwpath('post_mesh2mask.png'),'-dpng',['-r' num2str(resolution)]);
    mask=imread(mwpath('post_mesh2mask.png'));
else
    pos=get(hf,'position');
    pos(3:4)=max(pos(3:4),output_size+20);
    set(hf,'position',pos);
    set(gca, 'Units','pixels','position',[1, 1, output_size(1), output_size(2)]);
    mask=getframe(gca);
    if(any(size(mask.cdata)<[output_size([2 1]) 3]))
        error('the requested rasterization grid is larger than the screen resolution');
    end
    mask=mask.cdata(1:output_size(2),1:output_size(1),:);
end
if(nargin<5)
    close(hf);
end
mask=int32(reshape(mask,[size(mask,1)*size(mask,2) size(mask,3)]));
[isfound,locb]=ismember(mask,floor(cm*255),'rows');
locb(isfound==0)=nan;

mask=rot90(reshape(locb,output_size([2 1]))');

if(nargout>=2)
    xi=mn(1)+df(1)/2:df(1):mx(1);
    yi=mn(2)+df(2)/2:df(2):mx(2);
    weight=barycentricgrid(node,face,xi,yi,mask);
    if(size(face,2)>=4)
        badidx=find(weight(1,:,:)<0 | weight(2,:,:)<0 | weight(3,:,:)<0);
        badidx=badidx(face(mask(badidx),3)~=face(mask(badidx),4));
        weight2=barycentricgrid(node,face(:,[1 3 4]),xi,yi,mask);
        weight(:,badidx)=0;
        weight([1 3 4],badidx)=weight2(:,badidx);
    end
end

function weight=barycentricgrid(node,face,xi,yi,mask)
[xx,yy]=meshgrid(xi,yi);
idx=find(~isnan(mask));
eid=mask(idx);
t1=node(face(:,1),:);
t2=node(face(:,2),:);
t3=node(face(:,3),:);
tt=(t2(:,2)-t3(:,2)).*(t1(:,1)-t3(:,1))+(t3(:,1)-t2(:,1)).*(t1(:,2)-t3(:,2));
w(:,1)=(t2(eid,2)-t3(eid,2)).*(xx(idx)-t3(eid,1))+(t3(eid,1)-t2(eid,1)).*(yy(idx)-t3(eid,2));
w(:,2)=(t3(eid,2)-t1(eid,2)).*(xx(idx)-t3(eid,1))+(t1(eid,1)-t3(eid,1)).*(yy(idx)-t3(eid,2));
w(:,1)=w(:,1)./tt(eid);
w(:,2)=w(:,2)./tt(eid);
w(:,3)=1-w(:,1)-w(:,2);
weight=zeros(3,size(mask,1),size(mask,2));
ww=zeros(size(mask));
ww(idx)=w(:,1);
weight(1,:,:)=ww;
ww(idx)=w(:,2);
weight(2,:,:)=ww;
ww(idx)=w(:,3);
weight(3,:,:)=ww;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i = 1:length(idx)
    distance = zeros(K, 1);
    for j = 1:K
        % TODO(SaveTheRbtz@): Can be vectorized as diff * diff'
        distance(j) = sum(sum((X(i, :) - centroids(j, :)) .^ 2 ));
    endfor
    [value, idx(i)] = min(distance);
endfor

% =============================================================

end

repo: machine-learning


text:
function test_example_effects_of_tapering

% MEM 4gb
% WALLTIME 00:10:00

%
%% Effects of tapering for power estimates
%
% A simple way of looking at how (multi-)tapering affects the estimate of power in your signal is by using a very simple simulated signal. The **[ft_freqsimulation](https://github.com/fieldtrip/fieldtrip/blob/release/ft_freqsimulation.m)** function allows you to quickly create a simulated signal with a well-defined frequency component in it. Subsequently you can use **[ft_freqanalysis ](https://github.com/fieldtrip/fieldtrip/blob/release/ft_freqanalysis.m)** with different taper settings to see the effect of tapering on your power estimate.
%
% Create and plot a simulated signal. The simulated data contains only one trial, with a length of one second and a 50Hz sine wave.
%
close all

cfg = [];
cfg.method  = 'superimposed';
cfg.fsample = 1000;
cfg.numtrl  = 1;
cfg.trllen  = 1;
cfg.s1.freq = 50;
cfg.s1.ampl = 1;
cfg.s1.phase = 0;
cfg.noise.ampl = 0;
data = ft_freqsimulation(cfg);
figure
plot(data.time{1}, data.trial{1}(1,:))

%
% Compare the power estimate using a single Hanning taper with a single dpss taper (cfg.taper).
%
cfg        = [];
cfg.method = 'mtmfft';
cfg.output = 'pow';
cfg.pad    = 'maxperlen';
cfg.foilim = [0 100];
cfg.taper  = 'hanning';
freq       = ft_freqanalysis(cfg, data);
figure
semilogy(freq.freq, freq.powspctrm(1,:), 'b-');

cfg.taper     = 'dpss';
cfg.tapsmofrq = 1;   % i.e. no real smoothing
freq          = ft_freqanalysis(cfg, data);
hold on
semilogy(freq.freq, freq.powspctrm(1,:), 'g-');
legend({'hanning', 'dpss'});

%
% Look at the effect of multitapering with various amounts of smoothing (cfg.tapsmofrq).
%
cfg           = [];
cfg.method    = 'mtmfft';
cfg.output    = 'pow';
cfg.pad       = 'maxperlen';
cfg.foilim    = [0 100];
cfg.taper     = 'dpss';
cfg.tapsmofrq = 1;
freq          = ft_freqanalysis(cfg, data);
figure
semilogy(freq.freq, freq.powspctrm(1,:), 'b-');

cfg.tapsmofrq = 5;
freq          = ft_freqanalysis(cfg, data);
hold on
semilogy(freq.freq, freq.powspctrm(1,:), 'g-');

cfg.tapsmofrq = 10;
freq          = ft_freqanalysis(cfg, data);
hold on
semilogy(freq.freq, freq.powspctrm(1,:), 'r-');

legend({'1 Hz', '5 Hz', '10 Hz'});

%
% Look at the effect of spectral leakage of frequencies that are in between the natural frequencies of your time segment (cfg.pad). Note that the power spectral density per sqrt(Hz) decreases due to the zero-padding (compare the green and the blue). Multiplying the time series with a scaling factor fixes this, and makes the power spectral density estimates easier to compare (red and blue).
%
cfg           = [];
cfg.method    = 'mtmfft';
cfg.output    = 'pow';
cfg.foilim    = [0 100];
cfg.taper     = 'dpss';
cfg.tapsmofrq = 2;
cfg.pad       = 1;  % this is the same as the actual length of the data segment
freq          = ft_freqanalysis(cfg, data);
figure
semilogy(freq.freq, freq.powspctrm(1,:), 'b.-');

cfg.pad       = 100;
freq          = ft_freqanalysis(cfg, data);
hold on
semilogy(freq.freq, freq.powspctrm(1,:), 'g-');

data.trial{1} = data.trial{1}*10;           % equate the power spectral density to the first one
freq          = ft_freqanalysis(cfg, data);
hold on
semilogy(freq.freq, freq.powspctrm(1,:), 'r-');

legend({'1s', '5s', '10s'});


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
%  deconvolution

%vectors y and x and the time step are defined in the m-file c422.m

hh=deconv(y,x)*(1/step);

plot(t,hh)
ylim([-.1 1.1]);
legend('impulse response h(t)')

author: Sable

repo: mcbench-benchmarks


text:
function value = r4_normal_01_cdf_inverse ( p )

%*****************************************************************************80
%
%% R4_NORMAL_01_CDF_INVERSE inverts the standard normal CDF.
%
%  Discussion:
%
%    The result is accurate to about 1 part in 10**7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2004
%
%  Author:
%
%    Original FORTRAN77 version by Michael Wichura.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Michael Wichura,
%    The Percentage Points of the Normal Distribution,
%    Algorithm AS 241,
%    Applied Statistics,
%    Volume 37, Number 3, pages 477-484, 1988.
%
%  Parameters:
%
%    Input, real P, the value of the cumulative probability densitity function.
%    0 < P < 1.  If P is not in this range, an "infinite" value is returned.
%
%    Output, real VALUE, the normal deviate value with the property that
%    the probability of a standard normal deviate being less than or
%    equal to the value is P.
%
  a = [ 3.3871327179, 50.434271938, 159.29113202, 59.109374720 ];
  b = [ 1.0, 17.895169469, 78.757757664, 67.187563600 ];
  c = [ 1.4234372777, 2.7568153900, 1.3067284816, 0.17023821103 ];
  const1 = 0.180625;
  const2 = 1.6;
  d = [ 1.0, 0.73700164250, 0.12021132975 ];
  e = [ 6.6579051150, 3.0812263860, 0.42868294337, 0.017337203997 ];
  f = [ 1.0, 0.24197894225, 0.012258202635 ];
  split1 = 0.425;
  split2 = 5.0;

  if ( p <= 0.0 )
    value = -r4_huge ( );
    return
  end

  if ( 1.0 <= p )
    value = r4_huge ( );
    return
  end

  q = p - 0.5;

  if ( abs ( q ) <= split1 )

    r = const1 - q * q;
    value = q * r4poly_value ( 4, a, r ) / r4poly_value ( 4, b, r );

  else

    if ( q < 0.0 )
      r = p;
    else
      r = 1.0 - p;
    end

    if ( r <= 0.0 )
      value = -1.0;
      error ( 'R4_NORMAL_CDF_INVERSE - Fatal error!' );
    end

    r = sqrt ( -log ( r ) );

    if ( r <= split2 )

      r = r - const2;

      value = r4poly_value ( 4, c, r ) / r4poly_value ( 3, d, r );

    else

      r = r - split2;

      value = r4poly_value ( 4, e, r ) / r4poly_value ( 3, f, r );

    end

    if ( q < 0.0 )
      value = -value;
    end

  end

  return
end


text:
function [metabolite_structure] = queryExposomeExplorer(metabolite_structure)
%
% the function will search for metabolite names
% http://exposome-explorer.iarc.fr/search?utf8=%E2%9C%93&query=2-aminophenol+sulfate&button=
%
%
% INPUT
% metabolite_structure  metabolite structure
% startSearch           specify where the search should start in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% endSearch             specify where the search should end in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
%
% OUTPUT
% metabolite_structure  updated metabolite structure
%
%
% Ines Thiele, 09/2021

Mets = fieldnames(metabolite_structure);

if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(Mets);
end

annotationSource = 'Exposome Explorer website by name search';
annotationType = 'automatic';

for i = startSearch : endSearch
    if isempty(metabolite_structure.(Mets{i}).lipidmaps) || length(find(isnan(metabolite_structure.(Mets{i}).lipidmaps)))>0
        metN = metabolite_structure.(Mets{i}).metNames;
        % make name fit for internet
        metN = regexprep(metN,'"','');
        metInt = regexprep(metN,' ','+') ;
       % http://exposome-explorer.iarc.fr/search?utf8=%E2%9C%93&query=2-aminophenol+sulfate&button=
        url=['http://exposome-explorer.iarc.fr/search?utf8=%E2%9C%93&query=' metInt];
        syst = urlread(url);
        exposomeExplorer = '';
        if exist('syst','var') && ~isempty(strfind(syst,'returned 1 result'))
            tok = split(syst,'a href="/compounds/');
            tok2 = split(tok{2},'>');
            tok3 = regexprep(tok2{1},'"','');
            exposomeExplorer = tok3;
            % no real information on there - could be useful for getting mol
            % files
            if ~isempty(exposomeExplorer)
                metabolite_structure.(Mets{i}).exposomeExplorer = exposomeExplorer;
                metabolite_structure.(Mets{i}).exposomeExplorer_source = [annotationSource,':',annotationType,':',':',datestr(now)];
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
close all;
clear all;
clc;
%rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

% Signal space
N = 1000;
% Number of measurements
M = 200;
K = 50;

omp_success = 0;
mc_omp_success = 0;
for nt=1:100
    % Sensing matrix
    Phi = spx.dict.simple.gaussian_dict(M, N);
    % Construct the signal generator.
    gen  = spx.data.synthetic.SparseSignalGenerator(N, K);
    % Generate bi-uniform signals
    x = gen.gaussian;
    % Measurement vectors
    y = Phi.apply(x);

    % OMP solver instance
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(Phi, K);
    % Solve the sparse recovery problem
    omp_result = solver.solve(y);
    % Solution vector
    z = omp_result.z;
    omp_stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
    %fprintf('OMP\n');
    %spx.commons.sparse.print_recovery_performance(omp_stats);

    % MC-OMP solver instance
    solver = spx.pursuit.single.MC_OMP(Phi, K);
    %solver.BranchingFactor = 4;
    %solver.MaxCandidatesToRetain = 8;
    % Solve the sparse recovery problem
    mcomp_result = solver.solve(y);
    % Solution vector
    z = mcomp_result.z;
    mcomp_stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
    %fprintf('MC-OMP\n');
    %spx.commons.sparse.print_recovery_performance(mcomp_stats);

    fprintf('K=%d, Trial: %d, OMP: %s, MC-OMP: %s\n', ...
        K, nt, spx.io.true_false_short(omp_stats.success), ...
        spx.io.true_false_short(mcomp_stats.success));
    omp_success = omp_success  + omp_stats.success;
    mc_omp_success = mc_omp_success  + mcomp_stats.success;
end
fprintf('TOTAL: OMP : %d, MC-OMP: %d\n', omp_success, mc_omp_success);

author: indigits

repo: sparse-plex


text:
function cube_plots ( )

%*****************************************************************************80
%
%% CUBE_PLOTS plots the surface and the R(Theta) function for a cube.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CUBE_PLOTS:\n' );
  fprintf ( 1, '  For a 3D cube defined by a 0/1 characteristic function,\n' );
  fprintf ( 1, '  plot the surface, and R(Theta),\n' );
  fprintf ( 1, '  using a centered point, and then an offcentered point.\n' );
%
%  1): Plot the surface, using a centered base point.
%
  n1 = 31;
  n2 = 61;
  theta1 = linspace ( 0.0,       pi, n1 );
  theta2 = linspace ( 0.0, 2.0 * pi, n2 );

  [ t1, t2 ] = meshgrid ( theta2, theta1 );
  r = zeros ( n1, n2 );

  x0 = [ 0.0, 0.0, 0.0 ];
  for i = 1 : n1
    for j = 1 : n2
      r(i,j) = bisect_characteristic ( x0, [ t1(i,j); t2(i,j) ], @cube_characteristic );
    end
  end

  x1 = r .* cos ( t1 );
  x2 = r .* sin ( t1 ) .* cos ( t2 );
  x3 = r .* sin ( t1 ) .* sin ( t2 );

  figure ( 1 )
  mesh ( x3, x2, x1, 'FaceColor', 'interp' );
  axis equal
  grid on
  xlabel ( '<---X--->', 'FontSize', 16 );
  ylabel ( '<---Y--->', 'FontSize', 16 );
  zlabel ( '<---Z--->', 'FontSize', 16 );
  title ( 'Cube transition surface', 'FontSize', 24 )
  hold off
  filename =  'cube_centered_surface.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );
%
%  2): Plot R as a function of Theta for a centered base point.
%
  figure ( 2 )
  mesh ( t1, t2, r, 'FaceColor', 'Interp' );
  axis ( [0, 2.0 * pi, 0.0, pi, 0.0, 1.5 ] );
  xlabel ( '<---Theta2--->', 'FontSize', 16 )
  ylabel ( '<---Theta1--->', 'FontSize', 16 );
  zlabel ( '<---R(Theta1,Theta2)--->', 'FontSize', 16 )
  title ( 'R(*,*) from base point to surface', 'FontSize', 24 )
  grid on
  filename = 'cube_centered_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );
%
%  3): Plot the surface, using an offcentered base point.
%
  n1 = 31;
  n2 = 61;
  theta1 = linspace ( 0.0,       pi, n1 );
  theta2 = linspace ( 0.0, 2.0 * pi, n2 );

  [ t1, t2 ] = meshgrid ( theta2, theta1 );
  r = zeros ( n1, n2 );

  x0 = [ +0.7, -0.2, -0.5 ];
  for i = 1 : n1
    for j = 1 : n2
      r(i,j) = bisect_characteristic ( x0, [ t1(i,j); t2(i,j) ], @cube_characteristic );
    end
  end

  x1 = r .* cos ( t1 );
  x2 = r .* sin ( t1 ) .* cos ( t2 );
  x3 = r .* sin ( t1 ) .* sin ( t2 );

  figure ( 3 )
  c = zeros ( n1, n2 );
  colormap ( 'gray' )
  mesh ( x3, x2, x1, c );
  axis equal
  grid on
  xlabel ( '<---X--->', 'FontSize', 16 );
  ylabel ( '<---Y--->', 'FontSize', 16 );
  zlabel ( '<---Z--->', 'FontSize', 16 );
  title ( 'Cube transition surface', 'FontSize', 24 )
  hold off
  filename =  'cube_offcentered_surface.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );
%
%  4): Plot R as a function of Theta for an offcentered base point.
%
  figure ( 4 )
  mesh ( t1, t2, r, 'FaceColor', 'Interp' );
  xlabel ( '<---Theta2--->', 'FontSize', 16 )
  ylabel ( '<---Theta1--->', 'FontSize', 16 );
  zlabel ( '<---R(Theta1,Theta2)--->', 'FontSize', 16 )
  title ( 'R(*,*) from base point to surface', 'FontSize', 24 )
  grid on
  filename = 'cube_offcentered_plot.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created plotfile "%s".\n', filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Note variation in R for the offcenter case:\n' );
  fprintf ( 1, '  Min R = %g\n', min ( min ( r ) ) );
  fprintf ( 1, '  Max R = %g\n', max ( max ( r ) ) );
  return
end


text:
function nameList = makeNamesUnique(nameList, varargin)
% Identifies duplicate names in a list of names (cell array)
% and prompts the user to input new names.
% Called by `cleanTmodel`, `verifyModel`, `addSEEDInfo`, calls `countUnique`.
%
% USAGE:
%
%    nameList = makeNamesUnique(nameList, [nameInfo])
%
% INPUTS:
%    nameList:    Cell array of entires that contains duplicates.
%
% OPTIONAL INPUTS:
%    nameInfo:    Cell array the same length as nameList that contains
%                 information pertaining to the matching entry in `nameList`. This
%                 array can help the renaming process.
%
% OUTPUTS:
%    nameList:    Version of nameList with all unique entires.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

if (nargin > 1) % Declare variables.
    nameInfo = varargin{1} ;
end

% List of unique names, sorted.
uniqList = unique(nameList) ;

% Automatically rename entities?
[names, cnt] = countUnique(nameList);
nameInd = find(cnt > 1);
if ~isempty(nameInd)

    fprintf('%d non-unique names. ', length(nameInd))
    renameFlag = input(['Rename automatically?\n' ...
                     'Otherwise manual renaming will proceed. (y/n): '],...
                      's') ;
    if strcmpi(renameFlag, 'y') || strcmpi(renameFlag, 'yes')
        for i = 1:length(nameInd)
           thisName = names{nameInd(i)} ;
           fprintf('Renaming non-unique name %s (%d instances): ',...
                thisName,cnt(nameInd(i)));
           searchString = regexprep(thisName, '\[', '\\[') ;
           searchString = regexprep(searchString,'\]', '\\]') ;
           searchString = ['^' searchString '$'] ;
           IDs = find(~cellfun(@isempty, regexp(nameList, searchString))) ;
           for j = 1:length(IDs)
               if regexp(nameList{IDs(j)}, '\[\w\]$')
                   nameList{IDs(j)} = [thisName(1:end - 3) '_' num2str(j) ...
                                       thisName(end - 2:end)] ;
               else
                   nameList{IDs(j)} = [thisName '_' num2str(j)] ;
               end
               fprintf('%s\t', nameList{IDs(j)}) ;
           end
           fprintf('\n')
        end
    end
end

%% Find duplicate names.
for iName = 1:length(uniqList)
    dupPos = find(strcmp(uniqList{iName}, nameList)) ;
    % If there are duplicates.
    if length(dupPos) > 1
        dupNo = length(dupPos) ;
        % Duplicate name
        dupName = uniqList{iName} ;
        fprintf('%s has duplicates:\n', dupName)
        % Go through each duplicate and print info.
        for iDup = 1:dupNo
            if exist('nameInfo','var')
                string = [num2str(dupPos(iDup)), '\t', ...
                          nameInfo{dupPos(iDup)}, '\n'] ;
                fprintf(string)
            else
                fprintf('%d\n', num2str(dupPos(iDup)))
            end
        end
        % Ask for how names should be renamed.
        for iDup = 1:dupNo
            while 1
                prompt = ['Rename ' num2str(dupPos(iDup)) ' as: '] ;
                newName = input(prompt,'s') ;
                alreadytherepos = strcmp(newName, nameList) ;
                alreadytherepos(dupPos) = false ;
                % Check to see if name is already taken
                if sum(alreadytherepos) == 0
                    nameList{dupPos(iDup)} = newName ;
                    break
                else
                    errorPrompt = 'Name already taken. Proceed (y/n)?: ' ;
                    proceed = input(errorPrompt, 's') ;
                    if strcmpi(proceed, 'y') || strcmpi(proceed, 'yes')
                        nameList{dupPos(iDup)} = newName ;
                        break
                    else
                        continue
                    end
                end
            end
        end
    end
end

%% Check to make sure it worked.
if length(unique(nameList)) ~= length(nameList)
    fprintf('ERROR: Duplicate names still exist!\n')
end


text:
function err = getL2error3NE1(node,elem,exactE,Eh,markedElem)
%% GETL2ERROR3NE1 L2 norm of the linear Nedelec edge element.
%
%  err = getL2error3NE(node,elem,exactE,Eh,markedElem);
%
% Example
%
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     maxIt = 4;
%     L2Err = zeros(maxIt,1);
%     N = zeros(maxIt,1);
%     for k = 1:maxIt
%         [node,elem] = uniformbisect3(node,elem);
%         [elem2dof,edge] = dof3edge(elem);
%         pde = Maxwelldata2;
%         uI = edgeinterpolate1(pde.exactu,node,edge);
%         L2Err(k) = getL2error3NE1(node,elem,pde.exactu,uI);
%         N(k) = length(uI);
%     end
%     r = showrate(N,L2Err,1,'b-+');
%     legend('||u-u_I||',['N^{' num2str(r) '}'],'LOCATION','Best');
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Sort elem to ascend ordering
elem = sort(elem,2);

%% Construct Data Structure
elem2dof = dof3edge(elem);
NT = size(elem,1); Ndof = max(elem2dof(:)); %N = size(node,1);
[Dlambda,volume] = gradbasis3(node,elem);
locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];

%% Compute square of the L2 error element-wise
[lambda,w] = quadpts3(3); % quadrature order is 3
nQuad = size(lambda,1);
err = zeros(NT,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);
    Ep = exactE(pxy);
    Ehp = zeros(NT,3);
    for k = 1:6 % for each basis
        i = locEdge(k,1); j = locEdge(k,2);
        % phi_k = lambda_iDlambda_j - lambda_jDlambda_i;
        Ehp = Ehp + repmat(Eh(elem2dof(:,k)),1,3).*...
                   (lambda(p,i)*Dlambda(:,:,j)-lambda(p,j)*Dlambda(:,:,i));
        % psi_k = lambda_iDlambda_j + lambda_jDlambda_i;
        Ehp = Ehp + repmat(Eh(elem2dof(:,k)+Ndof),1,3).*...
                   (lambda(p,i)*Dlambda(:,:,j)+lambda(p,j)*Dlambda(:,:,i));
    end
    err = err + w(p)*sum((Ep - Ehp).^2,2);
end
err = err.*volume;
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % L2 err on some marked region
end
err = sqrt(sum(err));

repo: Aero-matlab


text:
function matrix_exponential_test ( )

%*****************************************************************************80
%
%% MATRIX_EXPONENTIAL_TEST tests the MATRIX_EXPONENTIAL library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2013
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../c8lib' )
  addpath ( '../r8lib' )
  addpath ( '../test_matrix_exponential' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MATRIX_EXPONENTIAL_TEST:\n'  );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the MATRIX_EXPONENTIAL library.\n' );
  fprintf ( 1, '  The R8LIB library is needed.\n' );
  fprintf ( 1, '  The test needs the TEST_MATRIX_EXPONENTIAL library.\n' );

  matrix_exponential_test01 ( );
  matrix_exponential_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MATRIX_EXPONENTIAL_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../c8lib' )
  rmpath ( '../r8lib' )
  rmpath ( '../test_matrix_exponential' )

  return
end


text:
function [ x, seed ] = erlang_sample ( a, b, c, seed )

%*****************************************************************************80
%
%% ERLANG_SAMPLE samples the Erlang PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, integer C, the parameters of the PDF.
%    0.0D+00 < B.
%    0 < C.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  a2 = 0.0;
  b2 = b;
  x = a;
  for i = 1 : c
    [ x2, seed ] = exponential_sample ( a2, b2, seed );
    x = x + x2;
  end

  return
end


text:
function H = circlePlane3D( center, normal, radious, theintv, normalon, color, style )
%CIRCLEPLANE3D Summary of this function goes here
%--------------------------------------------------------------------------
%Generate a circle plane in 3D with the given center and radious
%The plane is defined by the normal vector
%theintv is the interval theta which allow you to control your polygon
%shape
% Example:,
%
%   circlePlane3D([0 0 0], [1 -1 2], 5, 0.2, 1, [0 0 1], '-');
%   circlePlane3D([3 3 -3],[0 1 1], 3, 0.1, 1, 'y', '-');
%
%   Cheng-Yuan Wu <ieda_wind@hotmail.com>
%   Version 1.00
%   Aug, 2012
%--------------------------------------------------------------------------
%generate circle polygon
t = 0:theintv:2*pi;
x = radious*cos(t);
y = radious*sin(t);
z = zeros(size(x));
%compute rotate theta and axis
zaxis = [0 0 1];
normal = normal/norm(normal);
ang = acos(dot(zaxis,normal));
axis = cross(zaxis, normal)/norm(cross(zaxis, normal));
% A skew symmetric representation of the normalized axis
axis_skewed = [ 0 -axis(3) axis(2) ; axis(3) 0 -axis(1) ; -axis(2) axis(1) 0];
% Rodrigues formula for the rotation matrix
R = eye(3) + sin(ang)*axis_skewed + (1-cos(ang))*axis_skewed*axis_skewed;
fx = R(1,1)*x + R(1,2)*y + R(1,3)*z;
fy = R(2,1)*x + R(2,2)*y + R(2,3)*z;
fz = R(3,1)*x + R(3,2)*y + R(3,3)*z;
%translate center
fx = fx+center(1);
fy = fy+center(2);
fz = fz+center(3);
H = fill3(fx, fy, fz, color);
if normalon == 1
    hold on;
    H = plot3([center(1) center(1)+normal(1)],[center(2) center(2)+normal(2)],[center(3) center(3)+normal(3)],style);
end

end


text:
function varargout = SIGN(varargin)
%SIGN (overloaded)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/SIGN CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'
        switch varargin{1}
            case {'graph','exact'}

                t = varargin{2};
                X = varargin{3};

                d1 = binvar(1);
                d2 = binvar(1);
                d3 = binvar(1);
                [M,m] = derivebounds(X);
                if isequal(getbase(X),[0 1]) && ismember(getvariables(X),yalmip('tempintvariables'))
                    % Numerically unstable case with integer X
                    % d1: Negative <= -1
                    % d2: 0
                    % d3: Postive >= 1
                    F = [X >= m*d1, X <= M*d3, X <= (1-d1)*(M+1)-1,X >=(1-d3)*(m-1)+1, m*(1-d2) <= X <= M*(1-d2), t == -d1 + d3,d1+d2+d3==1];
                else
                    % Numerically unstable case with continuous X
                    F = [X >= m*d1, X <= M*d3, X <= (1-d1)*M,X >=(1-d3)*m, m*(1-d2) <= X <= M*(1-d2), t == -d1 + d3,d1+d2+d3==1];
                end

                varargout{1} = F;
                varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none');
                varargout{3} = X;

            otherwise
                error('SDPVAR/SIGN called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/SIGN');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function showmatrix(A,node,theta)
%% SHOWMATRIX displays a matrix.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('theta','var'), theta = 0; end
hold on
[i,j,aij] = find(A);
% line([node(i,1)'; node(j,1)'],[node(i,2)'; node(j,2)'],'Color','k');
nnz = length(i);
u = zeros(nnz,1);
v = u;
x = 0.5*node(i,1)+0.5*node(j,1);
y = 0.5*node(i,2)+0.5*node(j,2);
if (size(node,2)==3) % 3-D
    z = 0.5*node(i,3)+0.5*node(j,3);
else   % 2-D
    z = zeros(nnz,1);
end
idx = ((aij<0) & (abs(aij)>theta*mean(abs(aij))));
plot(x(idx), y(idx), 'g*', 'MarkerSize', 6);
idx = ((aij>0) & (abs(aij)>theta*mean(abs(aij))));
plot(x(idx), y(idx), 'c*', 'MarkerSize', 6);
quiver3(x,y,z,u,v,aij,'lineWidth',2);
%text(x+0.002,y+0.002,z+0.01,num2str(aij,3),'FontSize',14);
hold off


text:
function [FRresultsTable, FRresults] = makeFRresultsTable(FRresults, resultsDirectory, resultsFileName, modelMetaData, tableFilename)
% Makes a table of FR results
%
% USAGE:
%
%    [FRresultsTable, FRresults] = makeFRresultsTable(FRresults, resultsDirectory, resultsFileName, modelMetaData, tableFilename)
%
% INPUT:
%    FRresults:           output of `checkRankFRdriver`
%
% OPTIONAL INPUTS:
%    resultsDirectory:    directory where output of `checkRankFRdriver` has been saved
%    resultsFileName:     filename where output of `checkRankFRdriver` has been saved
%    modelMetaData:       Cell array, where each row is metadata for one model
%                         with five columns: `species`, `version`, `fileName`, `PMID`, `doi`.
%                         See function `modelMetaData = modelCitations()`
%    tableFilename:       If provided, a the table of results is written out
%                         to a csv file, with specified filename
%
% OUTPUTS:
%    FRresultsTable:      table displaying the results of `checkRankFRdriver`
%    FRresults:           output of `checkRankFRdriver`
%
% .. Author: - Ronan M.T. Fleming 'Conditions for duality between fluxes and concentrations in biochemical networks

if isempty(FRresults)
    if ~exist('resultsFileName','var')
        %resultsFileName='FRresults_20150128T225813';
        resultsFileName='FRresults_20150130T011200';
    end
    %results directory
    if ~exist('resultsDirectory','var')
        resultsDirectory='/home/rfleming/Dropbox/modelling/results/FRresults/';
    end
    cd(resultsDirectory)
    load([resultsDirectory resultsFileName])

    nResults=length(FRresults);
    %filename order of results structure
    for k=1:nResults
        tmp=FRresults(k).modelFilename;
        FRresults(k).modelID=tmp(1:end-4);%take off .mat
    end
else
    nResults=length(FRresults);
end

%extra column and extra row for headings
FRresultsTable=cell(29,nResults+1); %todo, come back and set correct number

if exist('modelMetaData','var')
    if isempty(modelMetaData)
        clear modelMetaData;
    end
end

firstColumn=1;
k=1;
while k<=nResults
    if ~firstColumn %&& 0
        disp(FRresults(k).modelID)
    end
    i=1;
    %the results may not be in alphabetical order but the table should be
    %so find the alphabetical position in the table for result(k)
    %modelMetaData is sorted alphabetically by species
    %Each row of modelMetaData: species, modelID, PMID, doi;
    if ~firstColumn
        %search against the filename specific to each model
        if exist('modelMetaData','var')
            bool=strcmp(FRresults(k).modelID,modelMetaData(:,3));
            if any(bool)
                n=find(bool);
            end
        else
            n=k;
            warning('no metadata found for for model')
        end
    end

    if firstColumn
        FRresultsTable{i,1}='Species';
    else
        if exist('modelMetaData','var')
            FRresultsTable{i,n+1}=modelMetaData{bool,1};
        else
            FRresultsTable{i,n+1}='';
        end
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='Version';
    else
        if exist('modelMetaData','var')
            FRresultsTable{i,n+1}=modelMetaData{bool,2};
        else
            FRresultsTable{i,n+1}='';
        end
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='modelID';
    else
        if exist('modelMetaData','var')
            FRresultsTable{i,n+1}=modelMetaData{bool,3};
        else
            FRresultsTable{i,n+1}=FRresults(k).modelID;
        end
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='PMID or DOI';
    else
        if exist('modelMetaData','var')
            if ~isempty(modelMetaData{bool,4})
                %pubmed id
                FRresultsTable{i,n+1}=modelMetaData{bool,4};
            else
                %doi
                FRresultsTable{i,n+1}=modelMetaData{bool,5};
            end
        else
            FRresultsTable{i,n+1}='';
        end
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Reactants = # Rows of [S S_e]';
    else
        disp(k)
        FRresultsTable{i,n+1}=size(FRresults(k).model.S,1);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rank of reconstruction [F,R]';
    else
        FRresultsTable{i,n+1}=FRresults(k).rankFRvanilla;
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rank [S S_e]';
    else
        FRresultsTable{i,n+1}=FRresults(k).rankS;
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='min coefficient magnitude [S S_e]';
    else
        bool=FRresults(k).model.S~=0;
        A=abs(FRresults(k).model.S);
        FRresultsTable{i,n+1}=min(min(A(A~=0)));
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='max coefficient magnitude [S S_e]';
    else
        bool=FRresults(k).model.S~=0;
        A=abs(FRresults(k).model.S);
        FRresultsTable{i,n+1}=max(max(A(A~=0)));
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='min/max coefficient magnitude [S S_e]';
    else
        bool=FRresults(k).model.S~=0;
        A=abs(FRresults(k).model.S);
        FRresultsTable{i,n+1}=min(min(A(A~=0)))/max(max(A(A~=0)));
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Elementally balanced rows (given formulae)';
    else
        if isfield(FRresults(k).model,'balancedMetBool')
            FRresultsTable{i,n+1}=nnz(FRresults(k).model.balancedMetBool);
        else
            FRresultsTable{i,n+1}=NaN;
        end
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Stoich. consistent rows';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentMetBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Nontrivial stoich. consistent rows';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentMetBool & FRresults(k).model.FRuniqueRowBool & FRresults(k).model.FRnonZeroRowBool1);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Nontrivial stoich. and flux consistent rows of [F,R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentMetBool & FRresults(k).model.FRuniqueRowBool & FRresults(k).model.FRnonZeroRowBool1 & FRresults(k).model.fluxConsistentMetBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Unique stoich. and nontrivial flux consistent rows of [F,R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentMetBool & FRresults(k).model.FRuniqueRowBool & FRresults(k).model.FRnonZeroRowBool1 & FRresults(k).model.fluxConsistentMetBool & FRresults(k).model.FRnonZeroRowBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rows of proper [F,R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.FRrows);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='  Rank of nontrivial stoich. and flux consistent [F,R]';
    else
        FRresultsTable{i,n+1}=FRresults(k).rankFR;
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rows of bilinear [F,R]';
    else
        FRresultsTable{i,n+1}=size(FRresults(k).model.Frb,1);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rank of bilinear [F,R]';
    else
        FRresultsTable{i,n+1}=FRresults(k).model.rankBilinearFrRr;
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Largest connected rows of [F,R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.largestConnectedRowsFRBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Exchange rows = # Rows exclusive to S_e';
    else
        FRresultsTable{i,n+1}=nnz(~FRresults(k).model.SIntMetBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Reactions = # Cols of [S S_e]';
    else
        FRresultsTable{i,n+1}=size(FRresults(k).model.S,2);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Exchange cols = # Cols of S_e';
    else
        FRresultsTable{i,n+1}=nnz(~FRresults(k).model.SIntRxnBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Elementally balanced cols';
    else
        if isfield(FRresults(k).model,'balancedRxnBool')
            FRresultsTable{i,n+1}=nnz(FRresults(k).model.balancedRxnBool);
        else
            FRresultsTable{i,n+1}=NaN;
        end
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Stoichiometrially consistent cols';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentRxnBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Unique and stoichiometrially consistent cols';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentRxnBool & FRresults(k).model.FRuniqueColBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Unique stoich. and flux consistent cols of [F;R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentRxnBool & FRresults(k).model.FRuniqueColBool & FRresults(k).model.fluxConsistentRxnBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Nontrivial stoich. and flux consistent cols of [F;R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.SConsistentRxnBool & FRresults(k).model.FRuniqueColBool & FRresults(k).model.fluxConsistentRxnBool & FRresults(k).model.FRnonZeroColBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Largest connected cols of [F;R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.largestConnectedColsFRVBool);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Cols of proper [F;R]';
    else
        FRresultsTable{i,n+1}=nnz(FRresults(k).model.FRVcols);
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rank of proper [F;R]';
    else
        FRresultsTable{i,n+1}=FRresults(k).rankFRV;
    end
    i=i+1;
    if firstColumn
        FRresultsTable{i,1}='# Rank of vanilla [F;R]';
    else
        FRresultsTable{i,n+1}=FRresults(k).rankFRVvanilla;
    end
    i=i+1;
%     if firstColumn
%         FRresultsTable{i,1}='Coherence of S';
%     else
%         FRresultsTable{i,n+1}=FRresults(k).coherenceS;
%     end

    %now move to columns for results
    if firstColumn
        firstColumn=0;
    else
        k=k+1;
    end
end

FRresultsTable=cell2table(FRresultsTable);

if exist('tableFilename','var')
    writetable(FRresultsTable,tableFilename,'Delimiter','\t')
end


text:
function [x fval] =  runps(a,b,c,x0)
[x, fval] = patternsearch(@nestedfun,x0);
% Nested function that computes the objective function
    function y = nestedfun(x)
        y = (a - b*x(1)^2 + x(1)^4/3)*x(1)^2 + x(1)*x(2) + ...
         (-c + c*x(2)^2)*x(2)^2;
    end
end

author: hliangzhao

repo: Mathematical-Model-Implementation


text:
function test_subplot_2D
%TEST_SUBPLOT_2D Test subplot 2D with DRAGZOOM

figure('Units', 'pixels');
x = -pi*2:0.1:pi*2;
y1 = sin(x);
y2 = cos(x);
y3 = sin(x).^3 - cos(x).^2;

hax1 = subplot(3, 1, 1);
plot(hax1, x, y1, '.-r')

hax2 = subplot(3, 1, 2);
plot(hax2, x, y2, 'o-b')

hax3 = subplot(3, 1, 3);
plot(hax3, x, y3, '*-g')
legend('test')

% dragzoom([hax1; hax3]); % manage only axes 1 and axes 3
dragzoom()


text:
function test_bug1833

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_mesh ft_read_headshape

% at this moment (20 November 2012) this test script is known not to work
% and the bugzilla report is still open
return

% this test script borrows some data from another bug
cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1820'));

mesh = ft_read_headshape('tet_4layer_127_127_127.1.node');
mesh.tet = mesh.tet(1:10000,:); % prune the number of elements
ft_plot_mesh(mesh, 'vertexcolor', 'none', 'edgecolor', 'none');

mesh = ft_read_headshape('cube2mm3layervorwerk_ns_127_127_127.v');
mesh.hex = mesh.hex(1:1000,:); % prune the number of elements
ft_plot_mesh(mesh, 'vertexcolor', 'none', 'edgecolor', 'none');


text:
function [V,D,bnd,j,work] = laneig(A,nin,k,sigma,options)

%LANEIG  Compute a few eigenvalues and eigenvectors.
%   LANEIG solves the eigenvalue problem A*v=lambda*v, when A is
%   real and symmetric using the Lanczos algorithm with partial
%   reorthogonalization (PRO).
%
%   [V,D] = LANEIG(A)
%   [V,D] = LANEIG('Afun',N)
%
%   The first input argument is either a real symmetric matrix, or a
%   string containing the name of an M-file which applies a linear
%   operator to the columns of a given matrix.  In the latter case,
%   the second input argument must be N, the order of the problem.
%
%   The full calling sequence is
%
%   [V,D,ERR] = LANEIG(A,K,SIGMA,OPTIONS)
%   [V,D,ERR] = LANEIG('Afun',N,K,SIGMA,OPTIONS)
%
%   On exit ERR contains the computed error bounds.  K is the number of
%   eigenvalues desired and SIGMA is numerical shift or a two letter string
%   which specifies which part of the spectrum should be computed:
%
%   SIGMA            Specified eigenvalues
%
%   'AL'            Algebraically Largest
%   'AS'            Algebraically Smallest
%   'LM'            Largest Magnitude   (default)
%   'SM'            Smallest Magnitude  (does not work when A is an m-file)
%   'BE'            Both Ends.  Computes k/2 eigenvalues
%                   from each end of the spectrum (one more
%                   from the high end if k is odd.)
%
%   The OPTIONS structure specifies certain parameters in the algorithm.
%
%    Field name      Parameter                              Default
%
%    OPTIONS.tol     Convergence tolerance                  16*eps
%    OPTIONS.lanmax  Dimension of the Lanczos basis.
%    OPTIONS.v0      Starting vector for the Lanczos        rand(n,1)-0.5
%                    iteration.
%    OPTIONS.delta   Level of orthogonality among the       sqrt(eps/K)
%                    Lanczos vectors.
%    OPTIONS.eta     Level of orthogonality after           10*eps^(3/4)
%                    reorthogonalization.
%    OPTIONS.cgs     reorthogonalization method used        0
%                    '0' : iterated modified Gram-Schmidt
%                    '1' : iterated classical Gram-Schmidt
%    OPTIONS.elr     If equal to 1 then extended local      1
%                    reorthogonalization is enforced.
%
%   See also LANPRO, EIGS, EIG.

% References:
% R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998.
%
% B. N. Parlett, ``The Symmetric Eigenvalue Problem'',
% Prentice-Hall, Englewood Cliffs, NJ, 1980.
%
% H. D. Simon, ``The Lanczos algorithm with partial reorthogonalization'',
% Math. Comp. 42 (1984), no. 165, 115--142.

% Rasmus Munk Larsen, DAIMI, 1998

%%%%%%%%%%%%%%%%%%%%% Parse and check input arguments. %%%%%%%%%%%%%%%%%%%%%%

if ~isstr(A)
  if nargin<1
    error('Not enough input arguments.');
  end
  [m n] = size(A);
  Aisfunc = 0;
  if m~=n | ~isequal(A,A') | ~isreal(A)
    error('A must be real symmetric')
  end
  if nargin < 4 | isempty(sigma)
    options = [];
  else
    options = sigma;
  end
  if nargin < 3 | isempty(k), sigma = 'LM'; else, sigma = k; end
  if nargin < 2 | isempty(nin), k = min(n,5); else, k = nin; end
else
  if nargin<2
    error('Not enough input arguments.');
  end
  Aisfunc = 1;
  n = nin;
  if nargin < 5 | isempty(options)
    options.tol = 16*eps;
    options.lanmax = n;
    options.v0 = rand(n,1)-0.5;
  end
  if nargin < 4 | isempty(sigma), sigma = 'LM'; end
  if nargin < 3 | isempty(k), k = min(n,5);  end
end

if ~isnumeric(k) | real(abs(fix(k)))~=k | ~isnumeric(n) | real(abs(fix(n)))~=n
  error('Input arguments N and K must be positive integers.')
end

% Quick return for n<2  or k<1
if n < 1 | k<1
  if nargout < 2
    V = zeros(k,1);
  else
    V = eye(n,k);
    D = zeros(k,k);
    bnd =zeros(k,1);
  end
  return
end
if n == 1
  if ~Aisfunc
    D = A;
    V = 1;
    bnd = 0;
  else
    D = feval(A,1);
    V = 1;
    dnb = 0;
  end
  if nargout<2
    V=D;
  end
  return
end

% A is the matrix of all zeros (not detectable if A is a string)
if ~Aisfunc
  if nnz(A)==0
    if nargout < 2
      V = zeros(k,1);
    else
      V = eye(n,k);
      D = zeros(k,k);
      bnd =zeros(k,1);
    end
    return
  end
end

lanmax = n;
tol = 16*eps;
r = rand(n,1)-0.5;
part = sigma;
% Parse options struct
if ~isempty(options) & isstruct(options)
  c = fieldnames(options);
  for i=1:length(c)
    if strmatch(c(i),'v0'), r = getfield(options,'v0'); r=r(:); end
    if strmatch(c(i),'tol'), tol = getfield(options,'tol'); end
    if strmatch(c(i),'lanmax'), lanmax = getfield(options,'lanmax'); end
  end
end

% Protect against absurd arguments.
tol = max(tol,eps);
lanmax = min(lanmax,n);
if size(r,1)~=n
  error('v0 must be a vector of length n')
end

lanmax = min(lanmax,n);
if k>lanmax
  error('K must satisfy  K <= LANMAX <= N.');
end
ksave = k;

if strcmp(sigma,'SM') & ~isstr(A)
  sigma = 0;
end

% Prepare for shift-and-invert if sigma is numeric.
if  isnumeric(sigma)
  part = 'LM';
  if isstr(A)
    error('Shift-and-invert works only when the matrix A is given explicitly.');
  else
    pmmd = symmmd(A);
    A = A(pmmd,pmmd);
    [S.L,S.U] = lu(A - sigma*speye(n));
    condU = condest(S.U);
    dsigma = n * full(max(max(abs(A)))) * eps;
    if sigma < 0
      sgnsig = -1;
    else
      sgnsig = 1;
    end
    sigitr = 1;
    while condU > 1/eps & ((dsigma <= 1 & sigitr <= 10) | ~isfinite(condU))
      disps1 = sprintf(['sigma = %10e is near an exact eigenvalue of A,\n' ...
			'so we cannot use the LU factorization of (A-sigma*I): ' ...
			' condest(U) = %10e.\n'],sigma,condU);
      if abs(sigma) < 1
	sigma = sigma + sgnsig * dsigma;
	disps2 = sprintf('We are trying sigma + %10e = %10e instead.\n', ...
			 sgnsig*dsigma,sigma);
      else
	sigma = sigma * (1 + dsigma);
	disps2 = sprintf('We are trying sigma * (1 + %10e) = %10e instead.\n', ...
			 dsigma,sigma);
      end
      %     if nargout < 3 & dispn ~= 0
      disp([disps1 disps2])
      %     end
      [S.L,S.U] = lu(A - sigma*speye(n));
      condU = condest(S.U);
      dsigma = 10 * dsigma;
      sigitr = sigitr + 1;
    end
  end
  A = S;
end

neig = 0; nrestart=-1;
if ~strcmp(part,'BE')
  j = min(2*k+2,lanmax);
else
  j = min(k+1,lanmax);
end

%%%%%%%%%%%%%%%%%%%%% Here begins the computation  %%%%%%%%%%%%%%%%%%%%%%

V = []; T = []; anorm = []; work = zeros(1,2); rnorm=-1;

while neig < k
  %%%%%%%%%%%%%%%%%%%%% Compute Lanczos tridiagonalization %%%%%%%%%%%%%%%%%
  j = min(lanmax,j+1-mod(j,2));
  % "Trick" to avoid unwanted zero eigenvalues when laneig is used for
  % SVD calculations. (Nothing to if lanmax is odd, though.)

  if  ~isstr(A)
    [V,T,r,anorm,ierr,w] = lanpro(A,j,r,options,V,T,anorm);
  else
    [V,T,r,anorm,ierr,w] = lanpro(A,n,j,r,options,V,T,anorm);
  end
  work= work + w;

  if ierr<0 % Invariant subspace of dimension -ierr found.
    j = -ierr;
  end

  %%%%%%%%%%%%%%%%%% Compute eigenvalues and error bounds %%%%%%%%%%%%%%%%%%
  % Analyze T
  [D,top,bot,err] = tqlb([full(diag(T))],full([0;diag(T,1)]));
  %  if err>0
  %    printf(['TQLB failed. Eigenvalue no. %i did not converge in 30', ...
  %	  ' iterations'],err);
  %  end
  %  full(T)
  %  [P,D] = eig(full(T));
  %  D = diag(D);
  %  bot = P(end,:)';
  %  [P(1,:)' P(end,:)']
  [D,I] = sort(D);
  bot = bot(I);

  % Set simple error bounds
  rnorm = norm(r);
  bnd = rnorm*abs(bot);

  % Use Largest Ritz value to estimate ||A||_2. This might save some
  % reorth. in case of restart.
  anorm = max(abs(D));

  % Estimate gap structure and refine error bounds
  bnd = refinebounds(D,bnd,n*eps*anorm);

  %%%%%%%%%%%%%%%%%%% Check convergence criterion %%%%%%%%%%%%%%%%%%%%
  % Reorder eigenvalues according to SIGMA
  switch part
   case 'AS'
    IPART = 1:j;
   case 'AL'
    IPART = j:-1:1;
   case 'LM'
    [dummy,IPART] = sort(-abs(D));
   case 'BE'
    if j<k
      IPART=1:j;
    else
      mid = floor(k/2);
      par = rem(k,1);
      IPART = [1:mid,(j-mid-par):j]';
    end
   otherwise
    error(['Illegal value for SIGMA: ',part]);
  end
  D = D(IPART);  bnd = bnd(IPART);
  if isnumeric(sigma)
    D = sigma + 1./D;
  end

  % Check if enough have converged.
  neig = 0;
  for i=1:min(j,k)
    if bnd(i) <= tol*abs(D(i))
      neig = neig + 1;
    end
  end

  %%%%%%%%%%% Check whether to stop or to extend the Krylov basis? %%%%%%%%%%
  if ierr<0 % Invariant subspace found
    if j<k
      warning(['Invariant subspace of dimension ',num2str(j-1),' found.'])
    end
    break;
  end
  if j>=lanmax % Maximal dimension of Krylov subspace reached => Bail out!
    if neig<ksave
      warning(['Maximum dimension of Krylov subspace exceeded prior',...
	       ' to convergence.']);
    end
    break;
  end

  % Increase dimension of Krylov subspace and try again.
  if neig>0
    %    j = j + ceil(min(20,max(2,((j-1)*(k-neig+1))/(2*(neig+1)))));
    j = j + min(100,max(2,0.5*(k-neig)*j/(neig+1)));
  elseif neig<k
    %    j = j + ceil(min(20,max(8,(k-neig)/2)));
    j = max(1.5*j,j+10);
  end
  j = min(j+1,lanmax);
  nrestart = nrestart + 1;
end

%%%%%%%%%%%%%%%% Lanczos converged (or failed). Prepare output %%%%%%%%%%%%%%%
k = min(ksave,j);

if nargout>1
  j = size(T,1);
  [Q,D] = eig(full(T)); D = diag(D);
  [D,I] = sort(D);
  % Compute and normalize Ritz vectors (overwrite V to save memory).
  V = V*Q(:,I(IPART(1:k)));
  for i=1:k
    nq = norm(V(:,i));
    if isfinite(nq) & nq~=0 & nq~=1
      V(:,i) = V(:,i)/nq;
    end
  end
  [D,I] = sort(D);
  D = D(IPART(1:k));
  if isnumeric(sigma)
    D = sigma + 1./D;
    V(pmmd,:) = V;
  end
end

% Pick out desired part of the spectrum
if length(D)~=k
  D = D(1:k);
  bnd = bnd(1:k);
end

if nargout<2
  V = D;
else
  D = diag(D);
end

author: andrewssobral

repo: lrslibrary


text:
function [cl,cd] = bpolar(BladeState,AC)
% To make your own bpolar:
% BladeState is the state of the blade element with matrix components:
% BladeState =
%            M: Mach number
%        alpha: angle of attack (radians)
%            r: radius location (nondimensional)
%          psi: azimuth location (radians)
%     alphadot: rate of change of angle of attack (rad/s)
%           Re: Reynolds number
% all components on BladeState inputs are m x n, where n is the number of
% blade elements, and m is the uniformly distributed number of azimuth
% elements.
%
% NB: this function should run VERY fast, so avoid loops and extensive
% interpolation.
% see also lininterp1f
%
% AC is a structure with AC data (see documentation)
%
% cl and cd are the m x n arrays containing each blade element's lift and
% drag coefficients.
% note that output coefficients will be non-dimensionalized by a constant
% blade chord, typically the average blade chord.
%
% The blade polar function is responsible for capturing all important rotor
% phenomena, including reverse flow, retreating blade stall, and advancing
% tip critical mach number.
%
% to capture reverse flow effects, the drag polar should be able to handle
% alphas +/- 360 deg
%
% critical mach drag rise and stall should also be modeled.

% This file is a simple but reasonable drag polar
%alphadot is not used for this simple polar, but can be used for more
%complex historesis stall behavior, etc.

RootCutout = .15;
Clmax = 1.35;

%%
alpha = BladeState.alpha;
r = BladeState.r;
M = BladeState.M;

tau = 2*pi;

astallmax = Clmax/tau;

% put all values between 0 and 360
ind = alpha<0;
while any(ind(:))
    alpha(ind) = alpha(ind)+tau;
    ind = alpha<0;
end
ind = alpha>tau;
while any(ind(:))
    alpha(ind) = alpha(ind)-tau;
    ind = alpha>tau;
end

% flat plate:
% fpClmax = 1.05;
% fpCl = @(a) fpClmax*sin(2*a);
% fpCd = @(x) -.9*cos(2*x) +.92;
cl = 1.05*sin(2*alpha);
cd= -.9*cos(2*alpha)+.92;

ind = alpha<=astallmax;
cl(ind) = tau*alpha(ind);
cd(ind)=.006+.1737*alpha(ind).^2;

ind = alpha>=tau-astallmax;
cl(ind) = tau*(alpha(ind)-tau);
cd(ind)=.006+.1737*(alpha(ind)-tau).^2;

% Critical Mach
absalpha = abs(alpha);
absalpha(alpha>pi/2) = abs(alpha(alpha>pi/2)-pi);
absalpha(alpha>3*pi/2) = abs(alpha(alpha>3*pi/2)-tau);
Mcrit = max(.55,-.02*absalpha+.765);
% from ESDU AERO W.00.03.01

dM = abs(M)-Mcrit;
drise = 801.6*exp(-((dM-1.441)/.4146).^2);
drise(dM>.3147) = .5;
cd = cd + drise;

%tip loss
ind = r>.97;
cl(ind)=cl(ind).*(1-(r(ind)-.97)/.03);

%% root cutout
cl(r<RootCutout) = 0;
cd(r<RootCutout) = .1;


text:
function [GrayCode,GrayCodeLength] = ImprovedGenerateGrayCode( Range )
% Input: Range is vector storing the maximum state in each digit.
% Output; GrayCode is the generated Graycode,GrayCodeLength is the number
% of generated Gray code.
Len = size( Range,2 );
if isempty( find( Range>1 ) ) == 1 %#ok<EFIND>
    GrayCode = zeros(1,Len);  GrayCodeLength = 1;
    return
end
[ LocalGrayCode,GrayCodeLength ] = GenerateGrayCode( Range( Range > 1 ) );

 t = 1;GrayCode  = zeros( GrayCodeLength,Len );
for p = 1:Len
   if Range( p ) == 1
      GrayCode(:,p) = zeros( GrayCodeLength,1 );
   else
      GrayCode(:,p) = LocalGrayCode(:,t);
      t = t + 1;
   end
end

end


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

%%
clear; close all; clc;

%% Syntax
% |[M]=affineTransformationMatrixDirect(V1,V2);|

%% Description
% This function computes the affine transformation matrix (translation,
% rotation, and scaling). For the two point matches input sets V1 and V2.

%% Examples

%%
% Plot settings
fontSize=15;
faceAlpha=1;
edgeColor='k';

%% Example:

%% Determine the affine transformation to overlay two surfaces
% Below is a demonstration of the features of the |affineTransformationMatrixDirect| function

% Load example patch data
[F,Vd]=parasaurolophus;

%Translation
OR=[-1 2 -3]; %Translations
T  = [1 0 0 OR(1);...
      0 1 0 OR(2);...
      0 0 1 OR(3);...
      0 0 0 1];
%Rotation
a=[-0.25*pi 0.75*pi 0.1*pi]; %Euler angles
R  = eye(4,4); R(1:3,1:3)=euler2DCM(a);

%Scaling
v=[2 3 1]; %Scaling factors
S  = [v(1) 0    0    0;...
      0    v(2) 0    0;...
      0    0    v(3) 0;...
      0    0    0    1];

M_true  = T * R * S; %The true transformation matrix

%Point set 1
V1=Vd+0.5;

%Point set 2
V2=tform(M_true,V1);

%%
% Plotting input data

hf=cFigure;
title('The untransformed surfaces','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);

hp=patch('Faces',F,'Vertices',V1,'FaceColor','g','FaceAlpha',faceAlpha);

hp=patch('Faces',F,'Vertices',V2,'FaceColor','r','FaceAlpha',faceAlpha,'edgeColor',edgeColor);

camlight headlight;
set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal; box on;
drawnow;

%%
% Get the transformation matrix for the point matched data using |affineTransformationMatrixDirect|

[M_fit]=affineTransformationMatrixDirect(V1,V2);

V1f=tform(M_fit,V1);

%%
% Plotting results

hf=cFigure;
title('The green surfaces transformed towards the red','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);

hp=patch('Faces',F,'Vertices',V2,'FaceColor','r','FaceAlpha',0.5,'edgeColor',edgeColor);

hp=patch('Faces',F,'Vertices',V1f,'FaceColor','none','FaceAlpha',0.5,'edgeColor','g');
% hp=patch('Faces',F,'Vertices',V2ff,'FaceColor','none','FaceAlpha',0.5*2,'edgeColor','b');

set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal; box on;
drawnow;

%%
%
% <<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 varargout=rotmc2qc(varargin)
%
% For calling details please see v_rotmc2qc.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_rotmc2qc.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: rotmc2qc.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_rotmc2qc(varargin{:});
else
    v_rotmc2qc(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function varargout = myComplement(MENUS,menu_defined)
tf = logical(zeros(1,length(MENUS)));
for menu=1:length(menu_defined)
    tf = tf + strcmp(menu_defined(menu),MENUS);
end
varargout{1} = MENUS(logical(1-tf));


text:
function [ t_stop, y_stop ] = p30_stop ( neqn )

%*****************************************************************************80
%
%% P30_STOP returns the stopping point for problem p30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = 1.0;

  return
end


text:
function h = m_annotation( varargin )
% M_ANNOTATION Creates an annotation object
%   M_ANNOTATION Generally passes its arguments through to ANNOTATION,
%   bur accepts lon/lat arguments.
%
%   ANNOTATION(ANNOTATIONTYPE) creates a default annotation of type
%   ANNOTATIONTYPE in the current figure.  ANNOTATIONTYPE may be one of the
%   following:
%       'rectangle'
%       'ellipse'
%       'textbox'
%       'line'
%       'arrow'
%       'doublearrow' = two headed arrow
%       'textarrow' = arrow with text at tail end
%
%   ANNOTATION('rectangle',POSITION) creates a rectangle annotation at the
%   position specified in normalized figure units by the vector POSITION
%   ANNOTATION('ellipse',POSITION) creates an ellise annotation at the
%   position specified in normalized figure units by the vector POSITION
%   ANNOTATION('textbox',POSITION) creates a textbox annotation at the
%   position specified in normalized figure units by the vector POSITION
%
%   POSTITION = [LON_LEFT LAT_BOTTOM WIDTH HEIGHT] where WIDTH and HEIGHT
%   are normalized to the size of the axis.
%
%   ANNOTATION('line',LON,LAT) creates a line annotation with endpoints
%   specified in normalized figure coordinates by the vectors LON and LAT
%   ANNOTATION('arrow',LON,LAT) creates an arrow annotation with endpoints
%   specified in normalized figure coordinates by the vectors LON and LAT.
%   LON(1) and LAT(1) specify the position of the tail end of the arrow
%   and LON(2) and LAT(2) specify the position at the tip of the arrow head.
%   ANNOTATION('doublearrow',LON,LAT) creates a doublearrow annotation with
%   endpoints specified in normalized figure coordinates by the vectors LON
%   and LAT.
%   ANNOTATION('textarrow',LON,LAT) creates a textarrow annotation with
%   endpoints specified in normalized figure coordinates by the vectors LON
%   and LAT. LON(1) and LAT(1) specify the position of the tail end of the
%   arrow and LON(2) and LAT(2) specify the position at the tip of the arrow
%   head.
%
%   You should (optionally) call ORIENT first, and then WYSIWYG to see how
%   the printed version will look, because otherwise the location of the
%   annotation may not look right.
%
%   M_ANNOTATION(HANDLE,...) creates the annotation in the  AXES specified
%   by HANDLE (note - ANNOTATION requires a FIGURE handle)
%
%   H=M_ANNOTATION(...) returns a handle to the annotation object
%
%   The arguments to ANNOTATION can be followed by parameter/value pairs to
%   specify additional properties of the annotation object. The X and Y or
%   POSITION arguments to ANNOTATION can be omitted entirely, and all
%   properties specified using parameter/value pairs.
%
%   Examples: rh=annotation('rectangle',[-123 48 .3 .3]);
%             ah=annotation('arrow',[10 30],[-5 10],'Color','r');
%             th=annotation('textarrow',[-50 -30],[-89 -70],'String','ABC');
%

% R. Pawlowicz Nov/2017
%
%

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

% Is first argument an axis handle?
if nargin>0 && length(varargin{1})==1 && ishandle(varargin{1})
   if  strcmp(get(varargin{1},'type'),'axes')
      useax=varargin{1};
      varargin(1)=[];
   else
      error(['map: ' mfilename ':invalidAxesHandle'],...
            ' First argument must be an axes handle ');
   end
else
    useax=gca;
end
figH=get(useax,'parent');

annotype=[varargin{1} '   '];

switch lower(annotype(1:5))
   case {'recta','ellip','textb'}
       POS=varargin{2};
       LON=POS(1);
       LAT=POS(2);
       [X,Y]=m_ll2xy(LON,LAT,'clip','off');
       [fX,fY]=axistofig(useax,X,Y,MAP_VAR_LIST.xlims,MAP_VAR_LIST.ylims);
       varargin{2}=[fX fY POS(3) POS(4)];

   case {'line ','arrow','doubl','texta'}
       LON=varargin{2};
       LAT=varargin{3};
       [X,Y]=m_ll2xy(LON,LAT,'clip','off');
       [fX,fY]=axistofig(useax,X,Y,MAP_VAR_LIST.xlims,MAP_VAR_LIST.ylims);
       varargin{2}=fX;
       varargin{3}=fY;
   otherwise  % Just a pass through!
end

h=annotation(figH,varargin{:});

set(h,'tag','m_annotation');

if nargout==0
    clear h
end

end

% Transform into normalized FIGURE coordinates
function [fX,fY]=axistofig(useax,X,Y,xlims,ylims)

axun = get(useax,'Units');
set(useax,'Units','normalized');
%%axpos = get(useax,'Position');  % This doesn't work if dataspectratio set
axpos =plotboxpos(useax);         %...but this does the right thing.

fX= (X-xlims(1))/diff(xlims)*axpos(3) + axpos(1);
fY= (Y-ylims(1))/diff(ylims)*axpos(4) + axpos(2);

%% Restore axes units
set(useax,'Units',axun);

end

% The following function came from MATLAB File exchange
% https://www.mathworks.com/matlabcentral/fileexchange/9615-plotboxpos
%
function pos = plotboxpos(h)
%PLOTBOXPOS Returns the position of the plotted axis region
%
% pos = plotboxpos(h)
%
% This function returns the position of the plotted region of an axis,
% which may differ from the actual axis position, depending on the axis
% limits, data aspect ratio, and plot box aspect ratio.  The position is
% returned in the same units as the those used to define the axis itself.
% This function can only be used for a 2D plot.
%
% Input variables:
%
%   h:      axis handle of a 2D axis (if ommitted, current axis is used).
%
% Output variables:
%
%   pos:    four-element position vector, in same units as h

% Copyright 2010 Kelly Kearney

% Check input

if nargin < 1
    h = gca;
end

if ~ishandle(h) || ~strcmp(get(h,'type'), 'axes')
    error('Input must be an axis handle');
end

% Get position of axis in pixels

currunit = get(h, 'units');
set(h, 'units', 'pixels');
axisPos = get(h, 'Position');
set(h, 'Units', currunit);

% Calculate box position based axis limits and aspect ratios

darismanual  = strcmpi(get(h, 'DataAspectRatioMode'),    'manual');
pbarismanual = strcmpi(get(h, 'PlotBoxAspectRatioMode'), 'manual');

if ~darismanual && ~pbarismanual

    pos = axisPos;

else

    xlim = get(h, 'XLim');
    ylim = get(h, 'YLim');

    % Deal with axis limits auto-set via Inf/-Inf use

    if any(isinf([xlim ylim]))
        hc = get(h, 'Children');
        hc(~arrayfun( @(h) isprop(h, 'XData' ) & isprop(h, 'YData' ), hc)) = [];
        xdata = get(hc, 'XData');
        if iscell(xdata)
            xdata = cellfun(@(x) x(:), xdata, 'uni', 0);
            xdata = cat(1, xdata{:});
        end
        ydata = get(hc, 'YData');
        if iscell(ydata)
            ydata = cellfun(@(x) x(:), ydata, 'uni', 0);
            ydata = cat(1, ydata{:});
        end
        isplotted = ~isinf(xdata) & ~isnan(xdata) & ...
                    ~isinf(ydata) & ~isnan(ydata);
        xdata = xdata(isplotted);
        ydata = ydata(isplotted);
        if isempty(xdata)
            xdata = [0 1];
        end
        if isempty(ydata)
            ydata = [0 1];
        end
        if isinf(xlim(1))
            xlim(1) = min(xdata);
        end
        if isinf(xlim(2))
            xlim(2) = max(xdata);
        end
        if isinf(ylim(1))
            ylim(1) = min(ydata);
        end
        if isinf(ylim(2))
            ylim(2) = max(ydata);
        end
    end

    dx = diff(xlim);
    dy = diff(ylim);
    dar = get(h, 'DataAspectRatio');
    pbar = get(h, 'PlotBoxAspectRatio');

    limDarRatio = (dx/dar(1))/(dy/dar(2));
    pbarRatio = pbar(1)/pbar(2);
    axisRatio = axisPos(3)/axisPos(4);

    if darismanual
        if limDarRatio > axisRatio
            pos(1) = axisPos(1);
            pos(3) = axisPos(3);
            pos(4) = axisPos(3)/limDarRatio;
            pos(2) = (axisPos(4) - pos(4))/2 + axisPos(2);
        else
            pos(2) = axisPos(2);
            pos(4) = axisPos(4);
            pos(3) = axisPos(4) * limDarRatio;
            pos(1) = (axisPos(3) - pos(3))/2 + axisPos(1);
        end
    elseif pbarismanual
        if pbarRatio > axisRatio
            pos(1) = axisPos(1);
            pos(3) = axisPos(3);
            pos(4) = axisPos(3)/pbarRatio;
            pos(2) = (axisPos(4) - pos(4))/2 + axisPos(2);
        else
            pos(2) = axisPos(2);
            pos(4) = axisPos(4);
            pos(3) = axisPos(4) * pbarRatio;
            pos(1) = (axisPos(3) - pos(3))/2 + axisPos(1);
        end
    end
end

% Convert plot box position to the units used by the axis

hparent = get(h, 'parent');
hfig = ancestor(hparent, 'figure'); % in case in panel or similar
currax = get(hfig, 'currentaxes');

temp = axes('Units', 'Pixels', 'Position', pos, 'Visible', 'off', 'parent', hparent);
set(temp, 'Units', currunit);
pos = get(temp, 'position');
delete(temp);

set(hfig, 'currentaxes', currax);

end

author: goGPS-Project


text:
%LP_MAKER  Makes mixed integer linear programming problems.
%
%  SYNOPSIS: lp_handle = lp_maker(f,a,b,e,vlb,vub,xint,scalemode,setminim)
%     make the MILP problem
%       max v = f'*x
%         a*x <> b
%           vlb <= x <= vub
%           x(int) are integer
%
%  ARGUMENTS: The first four arguments are required:
%           f: n vector of coefficients for a linear objective function.
%           a: m by n matrix representing linear constraints.
%           b: m vector of right sides for the inequality constraints.
%           e: m vector that determines the sense of the inequalities:
%                     e(i) < 0  ==> Less Than
%                     e(i) = 0  ==> Equals
%                     e(i) > 0  ==> Greater Than
%         vlb: n vector of non-negative lower bounds. If empty or omitted,
%              then the lower bounds are set to zero.
%         vub: n vector of upper bounds. May be omitted or empty.
%        xint: vector of integer variables. May be omitted or empty.
%   scalemode: Autoscale flag. Off when 0 or omitted.
%    setminim: Set maximum lp when this flag equals 0 or omitted.
%
%  OUTPUT: lp_handle is an integer handle to the lp created.

function lp_handle = lp_maker(f,a,b,e,vlb,vub,xint,scalemode,setminim)

if nargin == 0
        help lp_maker;
        return;
end

[m,n] = size(a);
lp_handle = mxlpsolve('make_lp', m, n);
mxlpsolve('set_verbose', lp_handle, 3);
mxlpsolve('set_mat', lp_handle, a);
mxlpsolve('set_rh_vec', lp_handle, b);
mxlpsolve('set_obj_fn', lp_handle, f);
mxlpsolve('set_maxim', lp_handle); % default is solving minimum lp.

for i = 1:length(e)
  if e(i) < 0
        con_type = 1;
  elseif e(i) == 0
        con_type = 3;
  else
        con_type = 2;
  end
  mxlpsolve('set_constr_type', lp_handle, i, con_type);
end

if nargin > 4
  for i = 1:length(vlb)
    mxlpsolve('set_lowbo', lp_handle, i, vlb(i));
  end
end

if nargin > 5
  for i = 1:length(vub)
    mxlpsolve('set_upbo', lp_handle, i, vub(i));
  end
end

if nargin > 6
  for i = 1:length(xint)
    mxlpsolve('set_int', lp_handle, xint(i), 1);
  end
end

if nargin > 7
  if scalemode ~= 0
    mxlpsolve('set_scaling', lp_handle, scalemode);
  end
end

if nargin > 8
  if setminim ~= 0
    mxlpsolve('set_minim', lp_handle);
  else
    mxlpsolve('set_maxim', lp_handle);
  end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function cost = L2_norm(out, labels, batch_size)
    cost = sum((labels(:) - out(:)).^2) / 2;
    cost = cost / batch_size;
end

author: jimmy-ren


text:
function [R,eff] = randmio_und_connected(R, ITER)
%RANDMIO_UND_CONNECTED     Random graph with preserved degree distribution
%
%   R = randmio_und_connected(W,ITER);
%   [R eff] = randmio_und_connected(W, ITER);
%
%   This function randomizes an undirected network, while preserving the
%   degree distribution. The function does not preserve the strength
%   distribution in weighted networks. The function also ensures that the
%   randomized network maintains connectedness, the ability for every node
%   to reach every other node in the network. The input network for this
%   function must be connected.
%
%   Input:      W,      undirected (binary/weighted) connection matrix
%               ITER,   rewiring parameter
%                       (each edge is rewired approximately ITER times)
%
%   Output:     R,      randomized network
%               eff,    number of actual rewirings carried out
%
%   References: Maslov and Sneppen (2002) Science 296:910
%
%
%   2007-2012
%   Mika Rubinov, UNSW
%   Jonathan Power, WUSTL
%   Olaf Sporns, IU

%   Modification History:
%   Jun 2007: Original (Mika Rubinov)
%   Apr 2008: Edge c-d is flipped with 50% probability, allowing to explore
%             all potential rewirings (Jonathan Power)
%   Mar 2012: Limit number of rewiring attempts, count number of successful
%             rewirings (Olaf Sporns)

n=size(R,1);
[i,j]=find(tril(R));
K=length(i);
ITER=K*ITER;

% maximal number of rewiring attempts per 'iter'
maxAttempts= round(n*K/(n*(n-1)));
% actual number of successful rewirings
eff = 0;

for iter=1:ITER
    att=0;
    while (att<=maxAttempts)                                     %while not rewired
        rewire=1;
        while 1
            e1=ceil(K*rand);
            e2=ceil(K*rand);
            while (e2==e1)
                e2=ceil(K*rand);
            end
            a=i(e1); b=j(e1);
            c=i(e2); d=j(e2);

            if all(a~=[c d]) && all(b~=[c d])
                break           %all four vertices must be different
            end
        end

        if rand>0.5
            i(e2)=d; j(e2)=c; 	%flip edge c-d with 50% probability
            c=i(e2); d=j(e2); 	%to explore all potential rewirings
        end

        %rewiring condition
        if ~(R(a,d) || R(c,b))
            %connectedness condition
            if ~(R(a,c) || R(b,d))
                P=R([a d],:);
                P(1,b)=0; P(2,c)=0;
                PN=P;
                PN(:,d)=1; PN(:,a)=1;

                while 1
                    P(1,:)=any(R(P(1,:)~=0,:),1);
                    P(2,:)=any(R(P(2,:)~=0,:),1);
                    P=P.*(~PN);
                    if ~all(any(P,2))
                        rewire=0;
                        break
                    elseif any(any(P(:,[b c])))
                        break
                    end
                    PN=PN+P;
                end
            end %connectedness testing

            if rewire               %reassign edges
                R(a,d)=R(a,b); R(a,b)=0;
                R(d,a)=R(b,a); R(b,a)=0;
                R(c,b)=R(c,d); R(c,d)=0;
                R(b,c)=R(d,c); R(d,c)=0;

                j(e1) = d;          %reassign edge indices
                j(e2) = b;
                eff = eff+1;
                break;
            end %edge reassignment
        end %rewiring condition
        att=att+1;
    end %while not rewired
end %iterations

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function modelPlane = sc_extract_plane(imgName, img, mask, optA)

% SC_EXTRACT_PLANE:
%
% Extract plane model from an image
%
% Output: modelPlane
% The model plane has the following data structure
%
% modelPlane
% modelPlane.vpData:                        Vanishing point data
% modelPlane.numPlane                       Number of planes
% modelPlane.plane{indPlane}.vLine          The vanishing line of the plane
% modelPlane.plane{indPlane}.imgPlaneProb;  The planar location density
% modelPlane.plane{indPlane}.sourceVP       Which two VPs form the plane
% modelPlane.plane{indPlane}.rotPar(vpInd)  Rotation parameters for aligning
%                                           two sets of lines with the x-axis                                           x-axis
% modelPlane.plane{indPlane}.postProb       Posteior probability of the plane

modelPlane = [];

%% VP detection
imgPath = 'data';
vpFilePath = 'cache\vpdetection';
vpFileName  = [imgName(1:end-4), '-vanishingpoints.txt'];

if(~exist(fullfile(vpFilePath, 'text', vpFileName), 'file'))
    vpExeFile = fullfile('source', 'vpdetection.exe');
    vpDetectCMD = [vpExeFile, ' -indir ', imgPath, ' -infile ', imgName, ' -outdir ', vpFilePath];
    system(vpDetectCMD);
end
% Read vanishing point data
vpData = sc_read_vpdata(fullfile(vpFilePath, 'text', vpFileName));

modelPlane = sc_detect_plane_from_vp(vpData, img, mask, optA);

end

function vpData = sc_read_vpdata(fileName)

% SC_READ_VPDATA: read the data from vanishing point detection algorithm
% Input:
%   - fileName: the txt file containing the pre-computed vanishing point
%   detection code
% Output:
%   - vpData
%   The data structure of vpData
%       - vpData.numVP: number of detected vanishing points
%       - vpData.vp{i}.pos: the vanishing point position in the homogenous coordiante
%       - vpData.vp{i}.score: the score of the vanishing point
%       - vpData.vp{i}.numLines: number of lines supporting the vanishing point
%       - vpData.vp{i}.lines{j}.p1: (x1, y1): starting position
%       - vpData.vp{i}.lines{j}.p2: (x2, y2): ending position
%       - vpData.vp{i}.lines{j}.length: length of the line segment

vpData = [];

% Read data
fid = fopen(fileName);

%% Parse VP positions
temp = fscanf(fid, '%s ', [1 5]);
numVP = 0;
readVPFlag = 1;
VP = [];
while(readVPFlag)
    numVP = numVP + 1;
    vpCurr = fscanf(fid, '%g %g %g %g %g', [5 1]);
    if(~isempty(vpCurr))
        VP(:,numVP) = vpCurr;
    else
        temp = fscanf(fid, '%s ', [1 6]);
        readVPFlag = 0;
    end
end
VP = VP';

vpData.numVP = size(VP, 1);

% Save VP position data
for i = 1: vpData.numVP
    vpData.vp{i}.pos = VP(i, 1:3);
    vpData.vp{i}.score = VP(i, 4);
    vpData.vp{i}.numLines = VP(i, 5);
end

%% Parse each set of line segments for the corresponding VP
for i = 1: vpData.numVP
    numLine = fscanf(fid, '%d ', [1 1]);
    lines = fscanf(fid, '%g %g %g %g %g', [5 numLine]);
    vpData.vp{i}.lines = lines';
end

fclose(fid);
end

function modelPlane = sc_detect_plane_from_vp(vpData, img, mask, optA)

% SC_DETECT_PLANE_FROM_VP: simple plane detection algorithm
% Input:
%     - vpData: vanishing point data
%     - img: input image
%     - mask: hole mask
% Output:
%     - modelPlane

%%

modelPlane = [];

% === Setting up ===
[imgH, imgW, ch] = size(img);
HfilterX = fspecial('gaussian', [1, optA.filterSize], optA.filterSigma);
HfilterY = HfilterX';
% fspecial('gaussian', optA.filterSize, optA.filterSigma);

img = im2double(img);

% === Supporting lines spatial support estimation ===
shapeInserter = vision.ShapeInserter('Shape', 'Lines','BorderColor', 'White');

for i = 1: vpData.numVP
    % The support lines
    imgLines = zeros(imgH, imgW);
    imgLines = step(shapeInserter, imgLines, int16(round(vpData.vp{i}.lines(:,1:4))));
    % Spatial density estimation via blurring
    imgLinesPosMap = imgLines;
    for k = 1:optA.numFilterIter
        imgLinesPosMap = imfilter(imgLinesPosMap, HfilterX, 'conv', 'replicate');
    end
    for k = 1:optA.numFilterIter
        imgLinesPosMap = imfilter(imgLinesPosMap, HfilterY, 'conv', 'replicate');
    end

    % Save results
    modelPlane.vp{i}.imgLines = imgLines;
    modelPlane.vp{i}.imgLinesPosMap = imgLinesPosMap;
end

% === Estimate plane support and plane parameters ===
numPlane = (vpData.numVP)*(vpData.numVP-1)/2;
% Initialize plane data
modelPlane.plane = cell(numPlane, 1);

indPlane = 1;
% A pair of vanishing points forms a plane hypothesis
for i = 1: vpData.numVP - 1
    for j = i+1: vpData.numVP
        % Compute the vanishing line
        modelPlane.plane{indPlane}.vLine = vLineFromTwoVP(vpData.vp{i}.pos, vpData.vp{j}.pos);
        % Element-wise product of two support line density
        modelPlane.plane{indPlane}.imgPlaneProb = modelPlane.vp{i}.imgLinesPosMap.*modelPlane.vp{j}.imgLinesPosMap; % Product of two probability maps
        modelPlane.plane{indPlane}.score = sum(modelPlane.plane{indPlane}.imgPlaneProb(:));
        modelPlane.plane{indPlane}.sourceVP = [i, j];

        indPlane = indPlane + 1;
    end
end

% === Compute rectified rotation parameters ===

for i = 1: numPlane
    for vpInd = 1: 2

        linesCurr = vpData.vp{modelPlane.plane{i}.sourceVP(vpInd)}.lines;
        invalidLineInd = linesCurr(:,5) == 0;
        linesCurr = linesCurr(~invalidLineInd,:);
        numLines = size(linesCurr, 1);

        vLineCurr = modelPlane.plane{i}.vLine;

        % Rectified homography
        H = eye(3);
        H(3,:) = vLineCurr;

        linesStart = cat(2, linesCurr(:,1), linesCurr(:,2),ones(numLines, 1))';
        linesEnd = cat(2, linesCurr(:,3), linesCurr(:,4),ones(numLines, 1))';

        linesStartRect = H*linesStart;
        linesStartRect = linesStartRect./repmat(linesStartRect(3,:), 3, 1);

        linesEndRect = H*linesEnd;
        linesEndRect = linesEndRect./repmat(linesEndRect(3,:), 3, 1);

        linesVec = linesStartRect(1:2, :) - linesEndRect(1:2, :);
        linesSign = linesEndRect(2,:) > linesStartRect(2, :);
        linesSign = 2*linesSign - 1;
        linesLength = sqrt(sum(linesVec.^2, 1));
        linesCos = linesSign.*linesVec(1,:)./linesLength; % repmat(linesLength, 2, 1);

        theta = acos(linesCos);

        % Estimate average theta so that all the supporting lines aligned
        % with the x-axis
        thetaAvg = mean(theta, 2);
        for iter = 1: 5
            thetaDiff = theta - thetaAvg;
            indLargeTheta = thetaDiff > pi/2;
            theta(indLargeTheta) = pi - theta(indLargeTheta);

            indSmallTheta = thetaDiff < -pi/2;
            theta(indSmallTheta) = pi + theta(indSmallTheta);
            thetaAvg = mean(theta, 2);
        end

        thetaEst = thetaAvg;

        modelPlane.plane{i}.rotPar(vpInd) = thetaEst;
    end
end

% === Add a fronto-parallel plane ===

modelPlane.plane{indPlane}.vLine = [0 0 1];
modelPlane.plane{indPlane}.imgPlaneProb = optA.fpPlaneProb*ones(imgH, imgW);
modelPlane.plane{indPlane}.score = sum(modelPlane.plane{indPlane}.imgPlaneProb(:));
modelPlane.plane{indPlane}.rotPar(1) = 0;
modelPlane.plane{indPlane}.rotPar(2) = 0;

numPlane = numPlane + 1;

modelPlane.numPlane = numPlane;

% === Compute posterior probability ===

planeProb = zeros(imgH, imgW, numPlane);
for i = 1 : numPlane
    planeProb(:,:,i) = modelPlane.plane{i}.imgPlaneProb;
end
planeProbSum = sum(planeProb, 3);
planeProb = (planeProb)./repmat(planeProbSum, [1, 1, numPlane]);

modelPlane.postProbHole = planeProb;

% === Propagate posterior probability into the hole region ===

% Get border pixels
%     borderImg = im2double(mask);
%     borderImg = imfilter(borderImg, fspecial('gaussian'));
%     borderImg = (borderImg~=0) & (borderImg~=1);
%     borderImg = borderImg & ~mask;
[distMap, idMap] = bwdist(~mask, 'euclidean');

maskInt = mask;
maskInt(1,:) = 0;   maskInt(end,:) = 0;
maskInt(:,1) = 0;   maskInt(:,end) = 0;

for i = 1:numPlane
    planeProbCh = planeProb(:,:,i);
    planeProb(:,:,i) = planeProbCh(idMap);
    planeProb(:,:,i) = roifill(planeProb(:,:,i), maskInt);
end

planeProbSum = sum(planeProb, 3);
planeProb = (planeProb)./repmat(planeProbSum, [1, 1, numPlane]);

planeProbSum = 1 + numPlane*optA.probConst;
planeProb = (planeProb + optA.probConst)/planeProbSum;

if(0)

    % Get border pixels
    borderImg = im2double(mask);
    borderImg = imfilter(borderImg, fspecial('gaussian'));
    borderImg = (borderImg~=0) & (borderImg~=1);
    borderImg = borderImg & ~mask;

    holeImg = zeros(imgH, imgW) + inf;
    holeImg(borderImg) = 0;
    [~, neighbors] = vl_imdisttf(single(holeImg)) ;
    [v_, u_] = ind2sub([imgH, imgW], neighbors);

    [u, v] = meshgrid(1:imgW,1:imgH);

    holeProb = zeros(imgH, imgW, numPlane);
    for i = 1:numPlane
        planeProbCh = planeProb(:,:,i);
        holeProb(:,:,i) = planeProbCh(neighbors);
    end
    holeProbSum = sum(holeProb, 3);
    holeProb = holeProb./repmat(holeProbSum, [1, 1, numPlane]);

    maskPlane = repmat(mask, [1,1,numPlane]);
    planeProb = (~maskPlane).*planeProb + maskPlane.*holeProb;

end

modelPlane.postProb = planeProb;

end

function vLine = vLineFromTwoVP(vp1, vp2)

A = cat(1, vp1, vp2);

[U S V] = svd(A, 0);
vLine = V(:,end);
vLine = vLine/vLine(3); % [h7, h8, 1]

end

repo: StructCompletion