text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
function [sys,x0,str,ts] = quadrotor_plot(t,x,u,flag,s,plot,enable,vehicle)
    % Flyer plot, lovingly coded by Paul Pounds, first coded 17/4/02
    % version 2 2004 added scaling and ground display
    % version 3 2010 improved rotor rendering and fixed mirroring bug
    %
    % Displays X-4 flyer position and attitude in a 3D plot.
    % GREEN ROTOR POINTS NORTH
    % BLUE ROTOR POINTS EAST

    % PARAMETERS
    % s defines the plot size in meters
    % swi controls flyer attitude plot; 1 = on, otherwise off.

    % INPUTS
    % 1 Center X position
    % 2 Center Y position
    % 3 Center Z position
    % 4 Yaw angle in rad
    % 5 Pitch angle in rad
    % 6 Roll angle in rad

    % OUTPUTS
    %   None
    ts = [-1 0];

    if ~isfield(vehicle, 'nrotors')
        vehicle.nrotors = 4;    % sensible default for quadrotor function
    end

    switch flag
        case 0
            [sys,x0,str,ts] = mdlInitializeSizes(ts,plot,enable); % Initialization
        case 3
            sys = mdlOutputs(t,u,s,plot,enable, vehicle); % Calculate outputs
        case {1,2, 4, 9} % Unused flags
            sys = [];
        otherwise
            error(['unhandled flag = ',num2str(flag)]); % Error handling
    end

    % Initialize
function [sys,x0,str,ts] = mdlInitializeSizes(ts,plot,enable)
    % Call simsizes for a sizes structure, fill it in, and convert it
    % to a sizes array.
    sizes = simsizes;
    sizes.NumContStates  = 0;
    sizes.NumDiscStates  = 0;
    sizes.NumOutputs     = 0;
    sizes.NumInputs      = 6;
    sizes.DirFeedthrough = 1;
    sizes.NumSampleTimes = 1;
    sys = simsizes(sizes);
    x0  = [];
    str = [];          % Set str to an empty matrix.
    ts = [0.05 0];

    if enable == 1
        figure(plot);
        clf;
        %colordef(1,'none');
        %global anim
        %anim = Animate('quad-movie');

        %TODO enable animation saving from the block
    end
    % End of mdlInitializeSizes.

function sys = mdlOutputs(t,u,s, plot, enable, quad)
    global a1s b1s

    % not quite sure what this is about -- PIC
    if numel(a1s) == [0];
        a1s = zeros(1, quad.nrotors);
        b1s = zeros(1, quad.nrotors);
    end

    % vehicle dimensons
    d = quad.d; %Hub displacement from COG
    r = quad.r; %Rotor radius

    for i = 1:quad.nrotors
        theta = (i-1)/quad.nrotors*2*pi;
        %   Di      Rotor hub displacements (1x3)
        % first rotor is on the x-axis, clockwise order looking down from above
        D(:,i) = [ d*cos(theta); d*sin(theta); 0];
        scal = s(1)/4;
        %Attitude center displacements
        C(:,i) = [ scal*cos(theta); scal*sin(theta); 0];
    end

    if enable == 1
        %draw ground
        figure(plot);
        clf;
        if length(s) == 1
            axis([-s s -s s 0 s]);
        else
            axis([-s(1) s(1) -s(1) s(1) 0 s(2)])
            s = s(1);
        end
        hold on;

        % plot the ground boundaries and the big cross
        plot3([-s -s],[s -s],[0 0],'-b')
        plot3([-s s],[s s],[0 0],'-b')
        plot3([s -s],[-s -s],[0 0],'-b')
        plot3([s s],[s -s],[0 0],'-b')
        plot3([s -s],[-s s],[0 0],'-b')
        plot3([-s s],[-s s],[0 0],'-b')
        grid on

        %READ STATE
        z = [u(1);u(2);u(3)];
        n = [u(4);u(5);u(6)];

        %PREPROCESS ROTATION MATRIX
        phi = n(1);    %Euler angles
        the = n(2);
        psi = n(3);

        R = [cos(the)*cos(phi) sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi) cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi);   %BBF > Inertial rotation matrix
            cos(the)*sin(phi) sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi) cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi);
            -sin(the)         sin(psi)*cos(the)                            cos(psi)*cos(the)];

        %Manual Construction
        %Q3 = [cos(psi) -sin(psi) 0;sin(psi) cos(psi) 0;0 0 1];   %Rotation mappings
        %Q2 = [cos(the) 0 sin(the);0 1 0;-sin(the) 0 cos(the)];
        %Q1 = [1 0 0;0 cos(phi) -sin(phi);0 sin(phi) cos(phi)];
        %R = Q3*Q2*Q1;    %Rotation matrix

        %CALCULATE FLYER TIP POSITONS USING COORDINATE FRAME ROTATION
        F = [1 0 0;0 -1 0;0 0 -1];

        %Draw flyer rotors
        t = [0:pi/8:2*pi];
        for j = 1:length(t)
            circle(:,j) = [r*sin(t(j));r*cos(t(j));0];
        end

        for i = 1:quad.nrotors
            hub(:,i) = F*(z + R*D(:,i)); %points in the inertial frame

            q = 1; %Flapping angle scaling for output display - makes it easier to see what flapping is occurring
            Rr = [cos(q*a1s(i))  sin(q*b1s(i))*sin(q*a1s(i)) cos(q*b1s(i))*sin(q*a1s(i));   %Rotor > Plot frame
                0              cos(q*b1s(i))               -sin(q*b1s(i));
                -sin(q*a1s(i)) sin(q*b1s(i))*cos(q*a1s(i)) cos(q*b1s(i))*cos(q*a1s(i))];

            tippath(:,:,i) = F*R*Rr*circle;
            plot3([hub(1,i)+tippath(1,:,i)],[hub(2,i)+tippath(2,:,i)],[hub(3,i)+tippath(3,:,i)],'b-')
        end

        %Draw flyer
        hub0 = F*z;  % centre of vehicle
        for i = 1:quad.nrotors
            % line from hub to centre plot3([hub(1,N) hub(1,S)],[hub(2,N) hub(2,S)],[hub(3,N) hub(3,S)],'-b')
            plot3([hub(1,i) hub0(1)],[hub(2,i) hub0(2)],[hub(3,i) hub0(3)],'-b')

            % plot a circle at the hub itself
            plot3([hub(1,i)],[hub(2,i)],[hub(3,i)],'o')
        end

        % plot the vehicle's centroid on the ground plane
        plot3([z(1) 0],[-z(2) 0],[0 0],'--k')
        plot3([z(1)],[-z(2)],[0],'xk')

        % label the axes
        xlabel('x');
        ylabel('y');
        zlabel('z (height above ground)');
        view(3)
    end

    sys = [];
    % End of mdlOutputs.

author: petercorke

repo: robotics-toolbox-matlab


text:
function output = sqrt(input)

output = cellfun(@sqrt, input, 'uniformoutput', false);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function outTable = convertToTableCore(featureLexicon,inParams)

% deal with optional arguments
if nargin<2 || isempty(inParams)
  inParams=struct();
end

%featureLexicon = ReadXMLParams(featureconfigfile);
allpf = fieldnames(featureLexicon.perframe);
pftypeList = {};
for pfndx = 1:numel(allpf)
  curpf = allpf{pfndx};
  curtypes  = featureLexicon.perframe.(curpf).type;
  if ischar(curtypes)
    curT = curtypes;
    if ~any(strcmp(pftypeList,curT))
      pftypeList{end+1} = curT;
    end
  else
    for tndx = 1:numel(curtypes)
      curT = curtypes{tndx};
      if ~any(strcmp(pftypeList,curT))
        pftypeList{end+1} = curT;
      end
    end
  end
end

types = fieldnames(inParams);
outTable = cell(numel(pftypeList),3);
for ndx = 1:numel(pftypeList)
  matches = strcmp(types,pftypeList{ndx});
  if any(matches)
    mndx = find(matches,1);
    outTable{ndx,1} = types{mndx};
    outTable{ndx,2} = inParams.(types{mndx}).mode{1};
    outTable{ndx,3} = inParams.(types{mndx}).selection{1};
  else
    outTable{ndx,1} = pftypeList{ndx};
    outTable{ndx,2} = 'None';
    outTable{ndx,3} = 'normal';
  end
end

author: kristinbranson

repo: JAABA


text:
function [Y, problem, S] = elliptope_SDP_complex(A, p, Y0)
% Solver for complex semidefinite programs (SDP's) with unit diagonal.
%
% function [Y, problem, S] = elliptope_SDP_complex(A)
% function [Y, problem, S] = elliptope_SDP_complex(A, p)
% function [Y, problem, S] = elliptope_SDP_complex(A, p, Y0)
%
% A is a Hermitian matrix of size n.
%
% This function uses a local optimization method in Manopt to solve the SDP
%
%   min_X trace(A*X) s.t. diag(X) = 1, X is complex, positive semidefinite.
%
% In practice, the Hermitian matrix X of size n is parameterized as
% X = Y*Y', where Y has size n x p. By default, p is taken large enough
% (that is, sqrt(n)) to ensure that there exists an optimal X whose rank is
% smaller than p. This ensures that the SDP is equivalent to the new
% problem in Y:
%
%   min_Y  trace(Y'*A*Y)  s.t.  diag(Y*Y') = 1, Y complex
%
% The constraints on Y require each row of Y to have unit norm, which is
% why Manopt is appropriate software to solve this problem. An optional
% initial guess can be specified via the input Y0.
%
% See the paper below for theory, specifically, for a proof that, for
% almost all A, second-order critical points of the problem in Y are
% globally optimal. In other words: there are no local traps in Y, despite
% non-convexity.
%
% Outputs:
%
%       Y: is the best point found (an nxp matrix with unit norm rows.)
%          To find X, form Y*Y' (or, more efficiently, study X through Y.)
%
%       problem: is the Manopt problem structure used to produce Y.
%
%       S: is a dual optimality certificate (a Hermitian matrix of size n,
%          sparse if A is sparse). The optimality gap (in the cost
%          function) is at most n*min(eig(S)), for both Y and X = Y*Y'.
%          Hence, if min(eig(S)) is close to zero, Y is close to globally
%          optimal. This can be computed via eigs(S, 1, 'SR').
%
% Paper: https://arxiv.org/abs/1606.04970
%
% @inproceedings{boumal2016bmapproach,
%   author  = {Boumal, N. and Voroninski, V. and Bandeira, A.S.},
%   title   = {The non-convex {B}urer-{M}onteiro approach works on smooth semidefinite programs},
%   booktitle={Neural Information Processing Systems (NIPS 2016)},
%   year    = {2016}
% }
%
% See also: maxcut elliptope_SDP

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Oct. 21, 2016
% Contributors:
% Change log:
%
%    Xiaowen Jiang Aug. 20, 2021
%       Added AD to compute the egrad and the ehess

    % If no inputs are provided, since this is an example file, generate
    % a random complex matrix. This is for illustration purposes only.
    if ~exist('A', 'var') || isempty(A)
        n = 100;
        A = randn(n) + 1i*randn(n);
        A = (A+A')/sqrt(2*n);
    end

    n = size(A, 1);
    assert(n >= 2, 'A must be at least 2x2.');
    assert(size(A, 2) == n, 'A must be square.');

    % Force A to be Hermitian
    A = (A+A')/2;

    % By default, pick a sufficiently large p (number of columns of Y).
    if ~exist('p', 'var') || isempty(p)
        p = floor(sqrt(n)+1);
    end

    assert(p >= 1 && p == round(p), 'p must be an integer >= 1.');

    % Pick the manifold of complex n-by-p matrices with unit norm rows.
    manifold = obliquecomplexfactory(p, n, true);

    problem.M = manifold;

    % These three, quick commented lines of code are sufficient to define
    % the cost function and its derivatives. This is good code to write
    % when prototyping. Below, a more advanced use of Manopt is shown,
    % where the redundant computation A*Y is avoided between the gradient
    % and the cost evaluation.
    % % problem.cost  = @(Y) .5*sum(sum(real((A*Y).*conj(Y))));
    % % problem.egrad = @(Y) A*Y;
    % % problem.ehess = @(Y, Ydot) A*Ydot;

    % Products with A dominate the cost, hence we store the result.
    % This allows to share the results among cost, grad and hess.
    % This is completely optional.
    function store = prepare(Y, store)
        if ~isfield(store, 'AY')
            AY = A*Y;
            store.AY = AY;
            store.diagAYYt = sum(real(AY .* conj(Y)), 2);
        end
    end

    % Define the cost function to be /minimized/.
    problem.cost = @cost;
    function [f, store] = cost(Y, store)
        store = prepare(Y, store);
        f = .5*sum(store.diagAYYt);
    end

    % Define the Riemannian gradient.
    problem.grad = @grad;
    function [G, store] = grad(Y, store)
        store = prepare(Y, store);
        G = store.AY - bsxfun(@times, Y, store.diagAYYt);
    end

    % If you want to, you can specify the Riemannian Hessian as well.
    problem.hess = @hess;
    function [H, store] = hess(Y, Ydot, store)
        store = prepare(Y, store);
        SYdot = A*Ydot - bsxfun(@times, Ydot, store.diagAYYt);
        H = manifold.proj(Y, SYdot);
    end

    % An alternative way to compute the egrad and the ehess is to use
    % automatic differentiation provided in the deep learning toolbox
    % (slower). AD does not support complex numbers if the Matlab version
    % is R2021a or earlier. The cost function should be defined differently
    % In this case. See complex_example_AD.m and manoptADhelp.m for more
    % information.
    % problem.cost = @cost_AD;
    %    function f = cost_AD(Y)
    %        AY = cprod(A, Y);
    %        diagAYYt = csum(creal(cdottimes(AY, cconj(Y))), 2);
    %        f = .5*csum(diagAYYt);
    %    end
    % Call manoptAD to automatically obtain egrad and ehess:
    % problem = manoptAD(problem);

    % If the version of Matlab installed is R2021b or later, specify the
    % cost function in the normal way and call manoptAD.
    % problem.cost = @cost_AD;
    %    function f = cost_AD(Y)
    %        AY = A*Y;
    %        diagAYYt = sum(real(AY .* conj(Y)), 2);
    %        f = .5*sum(diagAYYt);
    %    end
    % problem = manoptAD(problem);

    % If no initial guess is available, tell Manopt to use a random one.
    if ~exist('Y0', 'var') || isempty(Y0)
        Y0 = [];
    end

    % Call your favorite solver.
    opts = struct();
    opts.verbosity = 0;      % Set to 0 for no output, 2 for normal output
    opts.maxinner = 500;     % maximum Hessian calls per iteration
    opts.tolgradnorm = 1e-6; % tolerance on gradient norm
    Y = trustregions(problem, Y0, opts);

    % If required, produce an optimality certificate.
    if nargout >= 3
        S = A - spdiags(sum(real((A*Y).*conj(Y)), 2), 0, n, n);
    end

end

author: NicolasBoumal

repo: manopt


text:
function [t, stats] = upsampledttest(x,ufac,w);
% upsampledttest - t test for upsampled data
%
% [t, stat_struct] = upsampledttest(x,upsamplefactor,weight)
%
% T test for upsampled data (e.g. V1 in 1x1x1mm3 while the actual data was
% collected at some larger voxel size). The correction is based upon the
% assumption that upsampling duplicates the original dataset
% (upsamplefactor). This is true for nearest-neighbor interpolation and
% roughly true for linear interpolation schemes. If this correction is not
% taken into account the t-values will be overestimated.
% Option to specify weight for each value.
%
% example:
%   a = rand(100,1);
%   b = a*ones(1,100);b=b(:);
%    upsampledttest(a,1)
%   equals,
%    upsampledttest(b,100)
%
% 2006/06/01 SOD: wrote it.

if nargin < 2 || isempty(ufac),
  ufac = 1;
end;

if size(x,1) == 1, x = x(:); end;

% mean and n
xmean = mean(x);
n     = size(x,1);

% define corrected df
stats.df    = n./ufac-1;
% this could go below 0
if stats.df < 0,
    disp(sprintf(['[%s]:WARNING:corrected df < 0, setting to something ' ...
                  'small. Use data with caution.'],mfilename));
    stats.df  = 0.00001;
end;

% variance corrected for upsampling
xvar  = ( sum((x-xmean).^2) ./ufac) ./ stats.df;

% t (df corrected for upsampling)
t     = xmean ./ sqrt(xvar ./ (n./ufac));

% compute more stats if requested
if nargout > 1,
  stats.t     = t;
  stats.mean  = xmean;
  stats.stdev = sqrt(xvar);
  stats.sterr = sqrt(xvar)./sqrt(n./ufac);
  stats.p     = 0.5*betainc(stats.df./(stats.df+t.^2),stats.df/2,0.5);
end;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
filename='Bridge_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% X = NPAIRSK(N, K)
%
% Number of pair combinations for N objects to K pairs.
%
% For instance,
%
%   npairsk(2,1) = 1
%
% because two elements can be divided into a pair in one way.
%
%   npairsk(4,2) = 3
%
% because four elements (A,B,C,D) can be divided into two pairs in three
% ways: [(A,B);(C,D)], [(A,C);(B,D)], [(A,D);(B,C)].
%
%   npairsk(3,1) = 3
%
% because [(A,B)], [(A,C)], [(B,C)].

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

function x = npairsk(n, k)

x = factorial(n) ./ (factorial(2*k) .* factorial(n-2*k)) .* ngroupsk(k,2);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function varargout = TakeVideoAndPlayMovie(varargin)
% This program uses Image Acquisition toolbox to take video at any time
% length as user defined.  The function imaqtool has the limitation in taking
% long duration videos.
% By Jax Cao
% February 2008

% TAKEVIDEOANDPLAYMOVIE M-file for TakeVideoAndPlayMovie.fig
%      TAKEVIDEOANDPLAYMOVIE, by itself, creates a new TAKEVIDEOANDPLAYMOVIE or raises the existing
%      singleton*.
%
%      H = TAKEVIDEOANDPLAYMOVIE returns the handle to a new TAKEVIDEOANDPLAYMOVIE or the handle to
%      the existing singleton*.
%
%      TAKEVIDEOANDPLAYMOVIE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in TAKEVIDEOANDPLAYMOVIE.M with the given input arguments.
%
%      TAKEVIDEOANDPLAYMOVIE('Property','Value',...) creates a new TAKEVIDEOANDPLAYMOVIE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before TakeVideoAndPlayMovie_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to TakeVideoAndPlayMovie_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help TakeVideoAndPlayMovie

% Last Modified by GUIDE v2.5 24-Feb-2008 21:31:51

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @TakeVideoAndPlayMovie_OpeningFcn, ...
                   'gui_OutputFcn',  @TakeVideoAndPlayMovie_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

% --- Executes just before TakeVideoAndPlayMovie is made visible.
function TakeVideoAndPlayMovie_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to TakeVideoAndPlayMovie (see VARARGIN)

% Choose default command line output for TakeVideoAndPlayMovie
delete('*asv');
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes TakeVideoAndPlayMovie wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = TakeVideoAndPlayMovie_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in MovieSetupPushButton.
function MovieSetupPushButton_Callback(hObject, eventdata, handles)
% hObject    handle to MovieSetupPushButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Choose default command line output for TakeVideoAndPlayMovie
% This function determine which device is available and return appropriate
% paramters

% determine the camera divice ID
imaqreset;
infom = imaqhwinfo;
for k = 1:length(infom.InstalledAdaptors)
    info = imaqhwinfo(infom.InstalledAdaptors{k});
    if ~isempty(info.DeviceIDs)
        Installedadaptor = infom.InstalledAdaptors{k};
        break;
    end
end

if ~exist('Installedadaptor', 'var')
    h = warndlg('Video camera not connected !','!! Warning !!', 'replace');
    pause (5);
    if ishandle(h), close (h); end; beep;
    return;
end

% user input camera and video directory
prompt={'Enter the video directory name:',...
        'Enter the video filename prefix:', ...
        'Enter the video total length (seconds): ', ...
        'Enter the adaptorname: ', ...
        'Enter the frames per trigger: ', ...
        'Enter the time length (seconds) for each section', ...
        'Enter the device ID (the maximum should be the total number of cameras connected): '};
name='Input for video setup';
numlines=1;
defaultanswer={'c:/tempfiles/videotest/','myvedio', '400', Installedadaptor, '1500', ...
    '10', '1'};
options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';
answer=inputdlg(prompt,name,numlines,defaultanswer, options);

handles.loopnumber = ceil(str2double(answer{3})/str2double(answer{6}));
handles.directory = answer{1};
handles.filenameprefix = answer{2};
handles.adaptorname = answer{4};
handles.framepertrigger = str2double(answer{5});
handles.sectiontime = str2double(answer{6});
handles.deviceID = str2double(answer{7});
if ~exist(handles.directory , 'dir')
    mkdir(handles.directory);
else
    rmdir(handles.directory, 's');
    mkdir(handles.directory);
end;

info = imaqhwinfo(answer{4}, handles.deviceID);
[s,v] = listdlg('PromptString','Select a format:',...
    'SelectionMode','single',...
    'ListString',info.SupportedFormats);

if v
    handles.format = info.SupportedFormats{s};
else
    h = warndlg('Video format not selected, redo the work','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
end

beep;

handles.output = hObject;

% Update handles structure
guidata(hObject, handles);
% --- Executes on button press in TakeVideoPushButton.
function TakeVideoPushButton_Callback(hObject, eventdata, handles)
% hObject    handle to TakeVideoPushButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Choose default command line output for TakeVideoAndPlayMovie

if ~isfield(handles, 'format')
    h = warndlg('Video setting is not ready,  redo it','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
    beep;  return;
end

if isfield(handles, 'vid')
    h = warndlg('Video in preview mode stop preview first!','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
    beep; return;
end

% the following sections takes video according to the time specified by the
% user
k = 1;
while 1
    my_log = [handles.directory  handles.filenameprefix num2str(k) '.avi'];
    aviobj = avifile(my_log, 'compression', 'None');

    vid = videoinput(handles.adaptorname,handles.deviceID, handles.format);
    vid.LoggingMode = 'disk&memory';
    vid.DiskLogger = aviobj;
    vid.TriggerRepeat = Inf;
    set(vid,'FramesPerTrigger', handles.framepertrigger)

    start(vid);
    pause (handles.sectiontime)
    aviobj = close(vid.DiskLogger);
    flushdata(vid)
    delete(vid); clear vid;
    clear my_log;

    k = k +1;
    if k > handles.loopnumber, break; end;
end

beep;
% update handles
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% --- Executes on button press in PlayMoviePushButton.
function PlayMoviePushButton_Callback(hObject, eventdata, handles)
% hObject    handle to PlayMoviePushButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% This function plays the recorded videos.  However, it is recommended to
% use the windows movie maker to play the number of recorded videos
% exceeding 5.

if ~isfield(handles, 'format')
    h = warndlg('Video setting is not ready, redo it','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
    beep; return;
end

if isempty(ls([handles.directory, '*avi']))
    h = warndlg('Video not take yet, take video first!!!','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
    beep; return;
end

h = warndlg(['May be very slow, and it is suggested using Windows Move Maker' ...
    'to view the movie over 30 seconds!!'],'!! Warning !!');
pause (5);
if ishandle(h), close (h); end;

axes(handles.movieshow)
% mov = cell(1, loopnumber);
for k=1:handles.loopnumber
    mov = aviread([handles.directory  handles.filenameprefix num2str(k) '.avi']);
    movie(mov);
    clear mov
end

beep;

% --- Executes on button press in Previewpushbutton.
function Previewpushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Previewpushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isfield(handles, 'format')
    h = warndlg('Video setting is not ready,  redo it','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
    beep;     return;
end

handles.vid = videoinput(handles.adaptorname,handles.deviceID, handles.format);

% update handles
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

vidRes = get(handles.vid, 'VideoResolution');
nBands = get(handles.vid, 'NumberOfBands');
hImage = image( zeros(vidRes(2), vidRes(1), nBands) );
preview(handles.vid, hImage);

beep;

% --- Executes on button press in StopPreviewpushbutton.
function StopPreviewpushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to StopPreviewpushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% This function stops the video preview
if ~isfield(handles, 'format') || ~isfield(handles, 'vid')
    h = warndlg('Video setting is not ready or preview not started!','!! Warning !!');
    pause (5);
    if ishandle(h), close (h); end;
    beep; return;
end

closepreview(handles.vid)

handles = rmfield(handles, 'vid'); clear handles.vid;
beep;

% update handles
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

author: Sable

repo: mcbench-benchmarks


text:
% PROMAX - perform Promax oblique rotation after orthogonal Varimax
%            rotation of the rows of the input data. A method for
%            linear decomposition by "rotating to simple structure."
% Usage:
%        >> [R]   = promax(data,ncomps);
%        >> [R,V] = promax(data,ncomps,maxit);
%
% Inputs:
%     data    - Promax operates on rows of the input data matrix
%     ncomps  - operate on the N largest PCA components (default|0 -> all)
%     maxit   - maximum number of iterations {default|0 -> 5}
%
% Outputs:
%     R       - is the non-orthogonal Promax rotation matrix
%                 i.e.,  >> promax_rotated_data = R*data;
%     V       - is the orthogonal Varimax rotation matrix
%                 i.e.,  >> varimax_rotated_data = V*data;
%
% Author: Colin Humphries, CNL / Salk Institute, 1998
%
% See also: RUNICA

% Copyright (C) Colin Humphries, CNL / Salk Institute, June 1998
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% formatted and modified to return V by Scott Makeig, 6/23/98
% reset maxit default to 5, added ncomps -sm 7/8/98
% 01-25-02 reformated help & license, added links -ad
%
% Reference:
%
% Hendrickson AE and White PO (1964) Promax: A quick method for rotation
% to oblique simple structure, Br J of Stat Psych, X:xxx-xxx.

function [R,V] = promax(data,ncomps,maxit)

v = version;
indp = find(v == '.');
v = str2num(v(1:indp(2)-1));
if v >= 8.1
    disp('Note: for some unknown reason, this function does not return')
    disp('      NaN for a singular matrix under Matlab 2013a and later versions.')
    disp('      Promax on other matrices seems to as in previous revisions though.');
end

DEFAULT_POWER     = 4;
DEFAULT_TOLERANCE = 1e-5;
MAX_ITERATIONS    = 5;
NEAR_ZERO         = 1e-8;

powr = DEFAULT_POWER;
tol = DEFAULT_TOLERANCE;

if nargin < 1
  help promax
  return
end
if isempty(data)
  help promax
  return
end

if nargin < 2
   ncomps = 0;
end
chans = size(data,1)
if ncomps == 0
   ncomps = chans
end
if ncomps > chans
    error(sprintf('promax(): components must be <= number of data rows (%d).\n',chans));
end
if nargin < 3
  maxit = 0;
end
if maxit == 0
  maxit = MAX_ITERATIONS;
end

if ncomps < chans
  [eigenvectors,eigenvalues,compressed,datamean] = pcsquash(data,ncomps);
  data = compressed;
  clear compressed;
  eigenvectors = eigenvectors(:,1:ncomps); % make non-square
  eigenwts = pinv(eigenvectors); % find forward (non-square) weight matrix
end

R = varimax(data); % run Varimax on the (projected) data
B = R*data;        % compute rotated data
V = R;             % save Varimax matrix as V
if ncomps < chans
 V = V*eigenwts;   % include PCA reduction matrix
end
B = B';            % transpose
R = R';
cont = 1;
fprintf(...
  'Finding oblique Promax rotation using exponent %g and tolerance %g\n',...
                                                  powr,tol)
it = 1;
Pz = zeros(size(B));
while cont & it <= maxit
  P = Pz;
  ii = find(abs(B) > NEAR_ZERO); % avoid division by 0
  P(ii) = (abs(B(ii).^(powr+1)))./B(ii);
  tmp = inv(B'*B)*B'*P;
  tmp = normalcol(tmp);
  Rn = R*tmp;
  B = B*tmp;
  distnew = dot(Rn(:),R(:));
  if it > 1
    delta = abs(distnew-distold);
    if delta < tol
      cont = 0;
    end
    fprintf('#%d delta %f\n',it,delta)
    if isnan(delta)
      cont = 0;
    end
  end
  R = Rn;
  distold = distnew;
  it = it+1;
end
B = B';
R = R';
if ncomps < chans
   R = R*eigenwts; % include the pcsquash() compression
end

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

function n=normalcol(m)

if isempty(m)
  fprintf('normalcol() has empty input!\n');
  return
end
[mr,mc] = size(m);
n = sqrt(ones./sum(m.*m));
n = ones(mr,1)*n;
n = n.*m;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [P] = spm_mci_update_cov (P)
% Update covariance matrix of proposal density using Robbins-Monro
% FORMAT [P] = spm_mci_update_cov (P)
%
% See e.g.
% H. Haario, E. Saksman, and J. Tamminen. An adaptive Metropolis algorithm.
% Bernoulli, 7(2):223-242, 2001.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_update_cov.m 7679 2019-10-24 15:54:07Z spm $

Np=size(P.theta,1);
gamma=1/P.adapt_its;

x=P.theta(:,end);
dx=x-P.mu;

P.mu=P.mu+gamma*dx;
P.Ct=P.Ct+gamma*(dx*dx'-P.Ct);

author: spm


text:
function img = discretePolygon(varargin)
%DISCRETEPOLYGON Discretize a planar polygon
%
%   IMG = discretePolygon(DIM, POINTS)
%   DIM is the size of image, with the format [x0 dx x1;y0 dy y1]
%   POINTS is a N-by-2 array containing coordinate of polygon vertices.
%   Returns an image containing a discrete approximation of the polygon.
%
%   IMG = discretePolygon(LX, LY, ...);
%   Specifes the pixels coordinates with the two row vectors LX and LY.
%
%  See Also
%  imShapes, discretePolyline

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-05-16
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   04/03/2009: use meshgrid
%   29/05/2009: use more possibilities for specifying grid

% compute coordinate of image voxels
[lx, ly, varargin] = parseGridArgs(varargin{:});
[x, y]   = meshgrid(lx, ly);

% process input parameters
if length(varargin)==1
    var = varargin{1};
    px = var(:,1);
    py = var(:,2);

elseif length(varargin)==2
    px = varargin{1};
    py = varargin{2};
end

% compute discrete version of the polygon
img = inpolygon(x, y, px, py);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function mbasis = basis_matrix_overhauser_nonuni ( alpha, beta )

%*****************************************************************************80
%
%% BASIS_MATRIX_OVERHAUSER_NONUNI sets up the nonuniform Overhauser spline basis matrix.
%
%  Discussion:
%
%    This basis matrix assumes that the data points P1, P2, P3 and
%    P4 are not uniformly spaced in T, and that P2 corresponds to T = 0,
%    and P3 to T = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, BETA.
%    ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || )
%    BETA  = || P3 - P2 || / ( || P4 - P3 || + || P3 - P2 || ).
%
%    Output, real MBASIS(4,4), the basis matrix.
%
  mbasis(1,1) = - ( 1.0 - alpha ) * ( 1.0 - alpha ) / alpha;
  mbasis(1,2) =   beta + ( 1.0 - alpha ) / alpha;
  mbasis(1,3) =   alpha - 1.0 / ( 1.0 - beta );
  mbasis(1,4) =   beta * beta / ( 1.0 - beta );

  mbasis(2,1) =   2.0 * ( 1.0 - alpha ) * ( 1.0 - alpha ) / alpha;
  mbasis(2,2) = ( - 2.0 * ( 1.0 - alpha ) - alpha * beta ) / alpha;
  mbasis(2,3) = ( 2.0 * ( 1.0 - alpha ) ...
    - beta * ( 1.0 - 2.0 * alpha ) ) / ( 1.0 - beta );
  mbasis(2,4) = - beta * beta / ( 1.0 - beta );

  mbasis(3,1) = - ( 1.0 - alpha ) * ( 1.0 - alpha ) / alpha;
  mbasis(3,2) =   ( 1.0 - 2.0 * alpha ) / alpha;
  mbasis(3,3) =   alpha;
  mbasis(3,4) =   0.0;

  mbasis(4,1) =   0.0;
  mbasis(4,2) =   1.0;
  mbasis(4,3) =   0.0;
  mbasis(4,4) =   0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function results = testModelsAlgorithmically(solvers,matFolder,modelNames,printLevel)
%test a set of models for the feasibility of a nonzero objective(>1e-4)
%by default, do not run this script with testAll
%
%OPTIONAL INPUT
% solvers       z x 1 cell array of solvers to test e.g. {'gurobi5','quadMinos'}
% matFolder     absolute path to the folder where k model .mat files are
% modelNames    k x 2 cell array of {Species, modelFilename} setting order
%               of results structure when returned
%               e.g. modelNames={...
%                    'Lactococcus lactis','LLACTIS';
%                    'Mus musculus','iMM302'};
%
%printLevel    {(0),1} choose the amount to be printed to screen
%
%OUTPUT
% results      k x (2 + z) cell array of results

% 39 Jan 2016 Ronan M.T. Fleming

global CBT_LP_SOLVER
if 1
    %%parameters

    %select solver to use
    if ~exist('solvers','var')
        if ~isunix || ismac
            solvers={CBT_LP_SOLVER};
        else
            [status,cmdout]=system('which minos');
            if isempty(cmdout)
                [status,cmdout]=system('echo $PATH');
                disp(cmdout);
                warning('Minos not installed or not on system path.');
            else
                solvers={CBT_LP_SOLVER,'dqqMinos'};
            end
        end

        %gurobi (http://www.gurobi.com/) gurboi versions 5 and 6 are supported
        %solvers={'gurobi'};
        %optionally also test with quadMinos (a quadruple precision solver by Michael Saunders, Stanford Systems Optimization Laboratory)
        %solvers={'gurobi','dqqMinos'};
        solvers={'dqqMinos'};
    end

    pth=which('initCobraToolbox.m');
    global CBTDIR
    CBTDIR = pth(1:end-(length('initCobraToolbox.m')+1));

    if ~exist('matFolder','var')
        matFolder=pwd;
        % set the folder within the folder 'testModels' where the .xml files are located
        %matFolder=[CBTDIR '/testing/testModels/m_model_collection_mat']; %selection of mat files already parsed from https://github.com/opencobra/m_model_collection.git'
        %matFolder=['/usr/local/bin/cobratoolbox_master' '/testing/testModels/m_model_collection_mat'];
        %matFolder=[CBTDIR '/testing/testModels/testedModels'];
        %matFolder=[CBTDIR '/testing/testModels/testBiGG'];
        %matFolder='~/Dropbox/modelling/natureComm/source/BiGG_Models/mat';
        matFolder='~/Dropbox/modelling/natureComm/source/MSP_ME';
        matFolder='~/Dropbox/modelling/natureComm/source/MSP_ReconX/';
    end

    %modelNameStructure='MSP_ReconX';
    if ~exist('modelNames','var')
        modelNames=[];
    else
        if ~isempty(modelNames)
        if ~isstruct(modelNames)

        %choose the ordering of the table of results, the second column contains a
        %unique abbreviation for each model
        switch modelNames
            case ''
                %use the filenames
                modelNames=[];
            case 'MSP_ME'
                modelNames={...
                    'E. coli ME (Thiele et al)','ME_matrix_GlcAer_WT'};
            case 'MSP_ReconX'
                modelNames={...
                    'Homo sapiens','121114_Recon2betaModel';
                    'Homo sapiens','Recon2.v04'};
            case 'mongoose'
                modelNames={...
                    'Lactococcus lactis','LLACTIS';
                    'Halobacterium salinarum','';
                    'Mus musculus','';
                    'Streptomyces coelicolor','S_coilicolor';
                    'Mus musculus ','';
                    'Thermotoga maritima','T_Maritima';
                    'Clostridium acetobutylicum ','';
                    'Clostridium acetobutylicum ','';
                    'Corynebacterium glutamicum ','';
                    'Corynebacterium glutamicum ','';
                    'Geobacter metallireducens ','';
                    'Geobacter sulfurreducens ','';
                    'Mus musculus','mus_musculus';
                    'Lactobacillus plantarum','';
                    'Mannheimia succiniciproducens','';
                    'Mannheimia succiniciproducens','';
                    'Neisseria meningitidis','';
                    'Rhodoferax ferrireducens','';
                    'Staphylococcus aureus','';
                    'Streptococcus thermophilus','';
                    'Streptomyces coelicolor','';
                    'Natronomonas pharaonis','Natronomonas_pharaonis';
                    'Synechocystis sp. PCC6803 ','';
                    'Plasmodium falciparum','PlasmoNet';
                    'Acinetobacter baumannii','AbyMBEL891';
                    'Aspergillus oryzae','AORYZAE_COBRA';
                    'Arabidopsis thaliana','AraGEM';
                    'Mycobacterium tuberculosis','GSMN-TB';
                    'Acinetobacter baylyi','iAbaylyiV4';
                    'Acinetobacter baylyi','iAbaylyiV4';
                    'Acinetobacter baylyi','iAbaylyiV4';
                    'Acinetobacter baylyi ','iAbaylyiV4';
                    'Leishmania major','iAC560';
                    'Escherichia coli ','iAF1260';
                    'Methanosarcina barkeri','iAF692';
                    'Dehalococcoides ethenogenes','iAI549';
                    'Yersinia pestis','iAN818m';
                    'Bacillus subtilis','iBsu1103';
                    'Escherichia coli','iCA1273';
                    'Clostridium beijerinckii','iCB925';
                    'Helicobacter pylori','iCS291';
                    'Haemophilus influenzae','iCS400';
                    'Saccharomyces cerevisiae','iFF708';
                    'Buchnera aphidicola','iGT196';
                    'Aspergillus nidulans','iHD666 ';
                    'Saccharomyces cerevisiae','iIN800';
                    'Helicobacter pylori','iIT341';
                    'Synechocystis sp. PCC6803','iJN678';
                    'Escherichia coli','iJO1366';
                    'Pseudomonas putida','iJP815';
                    'Escherichia coli','iJR904';
                    'Burkholderia cenocepacia','iKF1028';
                    'Pichia pastoris','iLC915';
                    'Saccharomyces cerevisiae','iLL672';
                    'Aspergillus niger','iMA871';
                    'Salmonella typhimurium','iMA945';
                    'Methanosarcina acetivorans','iMB745';
                    'Staphylococcus aureus','iMH551';
                    'Saccharomyces cerevisiae','iMH805/775';
                    'Mus musculus','iMM1415';
                    'Saccharomyces cerevisiae','iMM904';
                    'Pseudomonas aeruginosa','iMO1056';
                    'Saccharomyces cerevisiae','iND750';
                    'Mycobacterium tuberculosis','iNJ661';
                    'Mycobacterium tuberculosis','iNJ661m';
                    'Pseudomonas putida','iNJ746';
                    'Cryptosporidium hominis','iNV213';
                    'Chromohalobacter salexigens','iOA584';
                    'Rhizobium etli','iOR363';
                    'Pichia pastoris','iPP668';
                    'Mycoplasma genitalium','iPS189';
                    'Chlamydomonas reinhardtii','iRC1080';
                    'Salmonella typhimurium','iRR1083';
                    'Zea mays','iRS1563';
                    'Arabidopsis thaliana','iRS1597';
                    'Francisella tularensis','iRS605';
                    'Rhodobacter sphaeroides','iRsp1095';
                    'Staphylococcus aureus','iSB619';
                    'Shewanella oneidensis','iSO783';
                    'Clostridium thermocellum','iSR432';
                    'Pichia stipitis','iSS884';
                    'Synechocystis sp. PCC6803','iSyn669';
                    'Plasmodium falciparum','iTH366';
                    'Scheffersomyces stipitis','iTL885';
                    'Porphyromonas gingivalis','iVM679';
                    'Methanosarcina acetivorans','iVS941';
                    'Ketogulonicigenium vulgare','iWZ663';
                    'Klebsiella pneumoniae','iYL1228';
                    'Zymomonas mobilis','iZM363';
                    'Bacillus subtilis ','model_v3';
                    'Pichia pastoris','PpaMBEL1254 ';
                    'Pseudomonas putida','PpuMBEL1071';
                    'Homo sapiens','121114_Recon2betaModel';
                    'Saccharomyces pombe','SpoMBEL1693';
                    'Salmonella typhimurium','STM_v1.0';
                    'Escherichia coli','textbook';
                    'Vibrio vulnificus','VvuMBEL943';
                    'Zymomonas mobilis','ZmobMBEL601'};
        end
        end
        end
    end
    %optionally convert the batch of .xml files into .mat files
    if 0
        %the folder where the SBML .xml files are located
        xmlFolder='m_model_collection';

        %folder where he .mat files are to be located
        matFolder=[pathContainingOpencobra 'cobratoolbox_master/testing/testModels/m_model_collection_mat'];

        %parsing via cobra toolbox and sbml toolbox and libsbml with matlab bindings
        sbmlTestModelToMat(xmlFolder,matFolder);
    end

    if ~exist('printLevel','var')
        printLevel=0;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%

    %The remainder of this code should run uninterupted,
    if any(strcmp(solvers,'quadMinos')) && isunix
        %quadMinos
        %addpath(pathContainingQuadMinos)
        %try to find the minos solver
        [status,cmdout]=system('which minos');
        if isempty(cmdout)
            [status,cmdout]=system('echo $PATH');
            disp(cmdout);
            warning('Minos not installed or not on system path.');
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % batch of mat files models the matFolder directory
    cd(matFolder)

    %batch of models in .mat format in a directory
    %assumes that each .mat file is a model
    matFiles=dir(matFolder);
    nModels=length(matFiles)-2;

    modelResults=cell(nModels,2+length(solvers));
    j=1;
    for k=3:length(matFiles) %loop through the mat files in the directory
        if strcmp(matFiles(k).name(end-2:end),'mat')

            if printLevel>0
                %disp(j)
                disp(matFiles(k).name)
            end

            %get the id of the model from the filename
            whosFile=whos('-file',matFiles(k).name);
            modelID=matFiles(k).name(1:end-4);
            modelResults{j,1}=modelID;
            if printLevel>-1
                fprintf('%20s\n',modelID);
            end
            load(matFiles(k).name);
            model=eval(whosFile.name);
            %stamp the model with the ID of the file
            model.modelID=modelID;
            model.description=modelID;

            %find the exchange reactions
            [m,n]=size(model.S);
            try
                if length(find(model.c~=0))~=1 % if more than one objective vector are defined in the model, only the one of the them is chosen for the testing.
                    ind=find(model.c~=0);
                    if isempty(ind); % if no objective vector is defined, set the first reaction as the objective for the testing.
                        warning('all zero entries in model.c')
                        model.c(1)=1;
                    else
                        warning('multiple nonzero entries in model.c')
                        model.c(ind)=0;
                        model.c(ind(1))=1; % by default, the first objective vector is used for the testing
                    end
                end
                model=findSExRxnInd(model,m,printLevel-1);
            catch
                disp('good');
            end

            if 0
                %check if stoichiometrically consistent
                [inform,m,model]=checkStoichiometricConsistency(model,printLevel-1);
            end

            %record the reaction to be optimized
            if any(model.c~=0) && isfield(model,'rxns')
                modelResults{j,2}=model.rxns{model.c~=0};
            else
                modelResults{j,2}='?';
            end

            %test with different solvers
            [out,solutions{j}]=testDifferentLPSolvers(model,solvers,printLevel);
            %save results
            for z=1:length(solvers)
                modelResults{j,2+z}=solutions{j}{z}.obj;
            end
            j=j+1;%used below
        end
    end

    %depending how the data is to be ordered, results structure is
    %different
    if isempty(modelNames)
        %create results structure
        results=cell(j,2+length(solvers));
        %add results in the order given by modelNames
        for k=1:j-1
            %Abbreviation of model
            results{k+1,1}=modelResults{k,1};
            results{k+1,2}=modelResults{k,2};
            %Abbreviation of reaction optimised'
            %solver objectives
            for z=1:length(solvers)
                results{k+1,2+z}=modelResults{k,2+z};
            end
        end
        %add headings to the results structure
        results{1,1}='Model';
        results{1,2}='Rxn maximised';
        for z=1:length(solvers)
            results{1,2+z}=solvers{z};
        end
    else
        %create results structure
        results=cell(size(modelNames,1)+1,3+length(solvers));
        %add results in the order given by modelNames
        for k=1:size(modelNames,1)
            %Abbreviation of model
            results{k+1,1}=modelNames{k,1};
            results{k+1,2}=modelNames{k,2};
            bool=strcmp(modelNames{k,2},modelResults(:,1));
            if any(bool)
                %Abbreviation of reaction optimised'
                results{k+1,3}=modelResults{bool,2};
                %solver objectives
                for z=1:length(solvers)
                    results{k+1,3+z}=modelResults{bool,2+z};
                end
            end
        end
        %add headings to the results structure
        results{1,1}='Species';
        results{1,2}='Model';
        results{1,3}='Rxn maximised';
        for z=1:length(solvers)
            results{1,3+z}=solvers{z};
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function dat= prep_powerspectrum (data, band, varargin)
%PROC_SPECTRUM -  calculate the power spectrum
%
%dat= proc_spectrum(dat, band, <win/N, step>)
%dat= proc_spectrum(dat, band, <opts>)
% IN   dat  - data structure of continuous or epoched data
%      band - frequency band
%      win  - window for FFT
%      N    - window width for FFT -> square window, default dat.fs
%      step - step for window (= # of overlapping samples), default N/2
%      opt  - struct of options:
%       .win       - window for FFT, default ones(dat.fs, 1)
%       .step      - step for window, default N/2
%       .db_scaled - boolean, if true values are db scaled (10*log10),
%                    default true

%%
% data
dat=data;
ban=band;
opt=opt_cellToStruct(varargin{:});
epo=struct('win',[],'N',[],'step',[],'scale',[]);

if isempty(dat)
    warning('Warning! data is empty');
end
if isempty(ban)
    warning('Band is not exist.');
end

if ~isfield(opt,'win')
   epo.win=dat.fs;
else
    epo.win=opt.win;
end

if ~isfield(opt,'N')
   epo.N=dat.fs;
else
    epo.N=opt.N;
end

if ~isfield(opt,'step')
   epo.step=dat.fs/2;
else
   epo.N=opt.N;
end

if ~isfield(opt,'scale')
   epo.scale='db';
else
    epo.scale=opt.scale;
end

[T, nEvents , nChans]= size(dat.x);

if length(epo.win)==1
    if epo.win>T
        warning('window legth is higher than signal')
    end
    epo.win=ones(epo.win,1);
end
N=length(epo.win);
normwin=norm(epo.win);
Freq=(0:N)/2*dat.fs/N;

%%
XX= zeros(N, nChans*nEvents);
nWindows= 1 + max(0, floor((T-N)/epo.step));
iv= 1:min(N, T);
Win= repmat(epo.win(:), [1 nChans*nEvents]);
bInd= band(1): ban(2);

%%calculate file

switch(lower(epo.scale)),
    case 'db',
        for iw= 1:nWindows,
            XX= XX + abs(fft(dat.x(iv,:).*Win, N)).^2;
            iv= iv + epo.step;
        end
        XX = XX/(nWindows*normwin^2);
        dat.x= reshape( 10*log10( XX(bInd,:)+eps ), [length(bInd), nChans, nEvents]);
        dat.yUnit= 'dB';
    case 'power',
        for iw= 1:nWindows,
            XX= XX + abs(fft(dat.x(iv,:).*Win, N).^2);
            iv= iv + epo.step;
        end
        dat.x= reshape(XX(bInd,:)/(nWindows*normwin^2), [length(bInd), nChans, nEvents]);
        dat.yUnit= 'power';
end

end

author: PatternRecognition

repo: OpenBMI


text:
%This function calculates BETA probabilities at each stage for all states,
%using GAMMA probabilities obtained previously. Uses recursion formula for
%BETA to calculate it for the previous stage. Each column is for states 00,10,01
%and 11 respectively. As we move backward in the block BETA will become very
%small, as the 1st term corresponding to gamma can be very less(of the order
%of 10^(-15)). Hence BETA will keep on decreasing and  will become very small.
%After some stages it will become exactly 0. So to avoid that we can
%multiply each BETA by 10^(-20) at a stage where they all become less than
%10^(-20). As we need BETA in calculation of LAPPR. So scaling wont affect the ratio

function [BETA]=beta_1(GAMMA,N)

    BETA=zeros(4,N);
    %Initialization assuming the final stage to be 00
    BETA(1,N)=1;BETA(2,N)=0;BETA(3,N)=0;BETA(4,N)=0;

    j=2*N-1;
    for i=N-1:-1:1
       BETA(1,i)=(GAMMA(1,j)*BETA(1,i+1))+(GAMMA(1,j+1)*BETA(2,i+1));
       BETA(2,i)=(GAMMA(2,j)*BETA(3,i+1))+(GAMMA(2,j+1)*BETA(4,i+1));
       BETA(3,i)=(GAMMA(3,j)*BETA(2,i+1))+(GAMMA(3,j+1)*BETA(1,i+1));
       BETA(4,i)=(GAMMA(4,j)*BETA(4,i+1))+(GAMMA(4,j+1)*BETA(3,i+1));
       j=j-2;

       if (BETA(1,i)<10^(-20) && BETA(2,i)<10^(-20) &&...
               BETA(3,i)<10^(-20) && BETA(4,i)<10^(-20) )
           BETA(:,i)=10^(20)*BETA(:,i);         %Scaling beta if became very less
       end
    end

end


text:
classdef MimSetWLTool < MimGuiPlugin
    % MimSetWLTool. Gui Plugin for using a preset bone window/level
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimSetWLTool is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Window / Level'
        SelectedText = 'Window / Level'
        ToolTip = 'Adjust window and level by dragging mouse over image'
        Category = 'Tools'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'windowlevel.png'
        Location = 12
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            if strcmp(gui_app.ImagePanel.Mode, MimModes.View3DMode)
                gui_app.ChangeMode([]);
            end
            gui_app.ImagePanel.SetControl('W/L');
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = strcmp(gui_app.ImagePanel.SelectedControl, 'W/L');
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function pathS = MC_VGGOU(S0,r,d,T,C,G,M,lambda,a,b,NTime,NSim,NBatches)
    intNt = 20;                 % steps in between orginial grid points
    allsteps = intNt * NTime;   % grid containing all necsseary steps
    dT = T / NTime;             % time step
    time = 0 : dT : T;          % tiem for martingale correction

    pathS = zeros(NSim,NTime+1,NBatches); % output
    lnS = zeros(NSim,NTime+1);      % used in batch
    lnS(:,1) = log(S0*exp(-d*T));   % set S(0) dividend adjusted

    %precompute the constants used for martingale correction
    y0 = 1;
    psiVG = (-1i)*C*log(G*M/(G*M+(M-G)-1));     % char exp                                     %characteristic exponent
    phiGOU = 1i*psiVG*y0/lambda*(1-exp(-lambda*time)) + ...
         lambda*a./(1i*psiVG-lambda*b).* ...
         (b*log(b./(b-1i*psiVG/lambda* ...
         (1-exp(-lambda*time))))-1i*psiVG*time); % char func
    omegaT = -phiGOU;           % martingale correction                                                         %martingale correction value
    omegaT(1) = 0;              % martingale correction in 0

for l = 1:NBatches              % batch loop
    yy = ones(NSim,allsteps+1); % init stochastic clock
    Np = poissrnd(a*lambda/allsteps*T,[NSim,allsteps]);

    for k = 1 : NSim
        for j = 1 : allsteps     % generating OU process
            if Np(k,j) > 0
                Ex = -log(rand(Np(k,j),1))/b; % exponential law
                U = exp(-lambda * T / allsteps * rand(Np(k,j),1));          % Uniforms
                yy(k,j+1) = (1-lambda*T/allsteps)*yy(k,j) ...
                    + sum(Ex .* U);
         else
                yy(k,j+1) = (1-lambda*T/allsteps)*yy(k,j);
            end
        end
    end

    ZZ = T*cumsum(yy,2)/allsteps;       %Integrated Time
    Y = zeros(NSim,NTime+1);
    for m = 2:NTime+1
        Y(:,m) = ZZ(:,(m-1)*intNt);
    end

    Intensity = C * (Y(:,2:end)-Y(:,1:end-1));      % intensity
    DGam = gamrnd(Intensity,1/M) - gamrnd(Intensity,1/G);
    diffomegaT = omegaT(2:end) - omegaT(1:end-1);   % maringale corr

    for m=2:NTime+1             % time loop
        lnS(:,m) = lnS(:,m-1) + (r-d)*dT ...
            + diffomegaT(m-1) + DGam(:,m-1);
    end
    pathS(:,:,l) = exp(lnS);
end


text:
function [CartPoint,latLong]=randEllipsoidLoc(N,a,f)
%%RANDELLIPSOIDLOC Generate a random point uniformly distributed on the
%              surface of a reference ellipsoid, such as the WGS-84
%              reference ellipsoid. The point is provided in Cartesian
%              coordinates as well as in terms of latitude and longitude.
%
%INPUTS: N The number of random samples on the reference ellipsoid to draw.
%        a The semi-major axis of the reference ellipsoid. If this argument
%          is omitted, the value in Constants.WGS84SemiMajorAxis is used.
%        f The flattening factor of the reference ellipsoid. If this
%          argument is omitted, the value in Constants.WGS84Flattening is
%          used.
%
%OUTPUTS: CartPoint A 3XN matrix of random Cartesian [x;y;z] points on the
%                   ellipsoid with the given semi-major axis and flattening
%                   factor.
%           latLong A 2XN matrix of ellipsoidal latitude and longitude
%                   corresponding to points in CartPoint.
%
%The problem of generating a uniformly distributed sample on a reference
%ellipsoid is number 28 in Chapter 3.4.1 of [1], where the algorithm to do
%so is given in the back of the book. That algorithm is implemented here,
%where the general ellipsoid formulation has been modified to support the
%typical parameterization in terms of a semimajor axis and a flattening
%factor. The algorithm is a rejection sampling method.
%
%Large flattening factor values can make the algorithm slow.
%
%REFERENCES:
%[1] D. Knuth, The Art of Computer Programming: Seminumerical Algorithms,
%    3rd ed. Reading, MA: Addison-Wesley, 1998, vol. 2.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3)
    f=Constants.WGS84Flattening;
end

if(nargin<2)
    a=Constants.WGS84SemiMajorAxis;
end

%The equation for a reference ellipsoid is
%(x^2+y^2)/a^2+z^2/(a^2*(1-f)^2)=1.
%The algorithm describe in Knuth's book wants things paramterized in terms
%of c(1)*x^2+c(2)*y^2+c(3)*z^2=1. Thus, the coefficients are
c(1,1)=1/a^2;
c(2,1)=1/a^2;
c(3,1)=1/(a^2*(1-f)^2);

%The c's have to be sorted. This means that the ordering of the components
%in x, y, and z will be wrong, so indices to undo the sorting when saving
%the final result need to be found.
[c,sortIdx]=sort(c,'descend');
unsortIdx=1:3;
unsortIdx=unsortIdx(sortIdx);

n=3;%The number of dimensions.

%Allocate space for the return variables.
CartPoint=zeros(3,N);

for curSamp=1:N
    while(1)
        %First, generate a random point on the unit sphere
        y=randDirVec(3);
        rho=sqrt(sum(c.*y.^2));

        %Find K
        if(n*c(n)>=c(1))
            K=sqrt(c(n)^(n-1));
        else
            K=sqrt(((n+1)/(c(1)+c(n)))^(n+1)*(c(1)*c(n)/n)^n);
        end

        %Generate a uniform random variable U
        U=rand(1);

        if(rho^(n+1)*U<K*sqrt(sum(c.^2.*y.^2)))
            CartPoint(unsortIdx,curSamp)=y/rho;
            break;
        end
    end
end

if(nargout>1)
    point=Cart2Ellipse(CartPoint,[],a,f);
    latLong=point(1:2,:);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef Compliance3DVoigtTensor < SymmetricFourthOrder3DVoigtTensor ...
                                 & ComplianceDescriptor

end


text:
function kpu_sparse_test ( )

%*****************************************************************************80
%
%% KPU_SPARSE_TEST uses the KPU function to build a sparse grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2012
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%    This MATLAB version by John Burkardt.
%
%  Local parameters:
%
%    Local, integer D, the spatial dimension.
%
%    Local, integer MAXK, the maximum level to check.
%
  d = 10;
  maxk = 4;
  func = 'prod( exp(-(x/2).^2/2)/2/sqrt(2*pi), 2)';
  trueval = fu_integral ( d );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'KPU_SPARSE_TEST:\n' );
  fprintf ( 1, '  KPU sparse grid:\n' );
  fprintf ( 1, '  Sparse nested, unweighted quadrature over [0,1].\n' );
  fprintf ( 1, '  Exact integral is %g\n', trueval );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   D  Level   Nodes    SG error    MC error\n' );
  fprintf ( 1, '\n' );

  for k = 2 : maxk
%
%  Compute sparse grid estimate.
%
    [ x, w ] = nwspgr ( 'kpu', d, k );
    fx = eval ( func );
    SGappr = w' * fx;
    SGerror = sqrt((SGappr - trueval).^2)/trueval;
%
%  Average 1000 Monte Carlo estimates.
%
    numnodes = length ( w );
    sim = zeros(1000,1);
    for r = 1 : 1000
      x = rand ( numnodes, d );
      fx = eval ( func );
      sim(r) = mean ( fx );
    end
    simerror = sqrt ( mean( ( sim - trueval).^2) ) / trueval;

    fprintf( '  %2d     %2d  %6d  %10.5g  %10.5g\n', ...
      d, k, numnodes, SGerror, simerror )

  end

  return
end


text:
function cm = gold(n)

if nargin == 0,
  n = 64;
end

golds = [255,215,0
  255,215,0
  238,201,0
  205,173,0
  139,117,0]/255;

cm = zeros(n,3);
for i = 1:3,
  cm(:,i) = interp1(linspace(1,n,size(golds,1)),golds(:,i),1:n,'spline');
end
cm = max(0,min(cm,1));


text:
%%
% Test for OT color transfer.

addpath('../toolbox/');
addpath('../colors_functions/');
addpath('../image_blur/');
addpath('../blur_functions//');
addpath('../convolutional_wasserstein/');
addpath('../../data/meshes/');
addpath('../mesh_functions/');

% #bins
N = 40;
N = 60;

%%
% helpers

mmin = @(x)min(x(:));
mmax = @(x)max(x(:));
normalize = @(h)h/sum(h(:));
setfigname = @(name)set(gcf, 'Name', name, 'NumberTitle','off');
% plot histograms
delta = .1/N^2;
dispHist = @(x)-log(x+delta);
dispCell = @(H)cellfun(dispHist, H, 'UniformOutput', false);
Entropy = @(x)-sum(x(x>0).*log(x(x>0)));

%%
% Blur kernel

mu = N/50;
mu = N/40;
mu = N/25;
blur = load_filtering('imgaussian', N);
K = @(x)blur(x,mu);
Kv = @(x)apply_3d_func(K,x);

%%
% Load densities

if 1 % not(exist('names'))
    names = {'spiky', 'spheres'};
    names = {'spiky', 'sphere'};
    names = {'sphere', 'boxes'};
    names = {'duck', 'horse'};
    names = {'duck' 'horse' 'shears' 'moomoo_s0'};
    names = {'spot_flipped','duck','torus'};
    names = {'spiky','sphere','boxes','cone_rotated'};
    names = {'duck','spiky','moomoo_s0','double-torus'};
    names = {'shears' 'chair' 'wolf0' 'duck' };
    names = {'mushroom' 'torus' 'hand1' 'trim-star'}; % 'star_subdivided'
    names = {'torus'  'duck' 'mushroom' 'spiky'}; % 'star_subdivided'
end
p = length(names);

f = {};
for i=1:p
    name = names{i};
    f{i} = normalize( load_volume(names{i}, N) );
    if isempty(f{i})
        %% try to load a mesh %%
        [V,F] = read_off([name '.off']);
        V = rescale(V, .03, .97);
        % check if rotation is needed
        A = eye(3);
        switch name
            case 'duck'
                A = [-0.0193 -0.2900 0.9568; 0.9677 0.2352 0.0908; -0.2514  0.9277 0.2761];
            case 'moomoo_s0'
                A = [-0.0971 0.7483 -0.6562; -0.9928 -0.1191 0.0111; -0.0699 0.6525 0.7545];
            case 'homer'
                A = [-0.4287  0.5566 -0.7116; -0.7577 -0.6505 -0.0524; -0.4921 0.5167 0.7006];
            case 'dinosaur'
                A = [-0.8362    0.0082   -0.5483   -0.2390   -0.9054    0.3509   -0.4935    0.4245    0.7591];
            case 'chair'
                A = [-0.8165    0.2728   -0.5089    0.3447   -0.4768   -0.8086   -0.4632   -0.8356    0.2952];
            case 'trim-star'
                A = [-0.7297    0.3243    0.6020   -0.6814   -0.4187   -0.6003    0.0574   -0.8483    0.5265];
            case 'hand1'
                A = [-0.4030   -0.7111   -0.5761    0.8945   -0.4391   -0.0837   -0.1934   -0.5491    0.8131];
            case 'mushroom'
                A = [-0.1212   -0.9430    0.3099    0.9227    0.0081    0.3855   -0.3660    0.3326    0.8691];
        end
        V = reshape(A,[3 3])*(V-.5)+.5;
        %  V0 = V; [A,~] = qr(randn(3)); V = A*(V0-.5)+.5;
        t = linspace(0,1,N);
        f{i} = inpolyhedron(F',V',t,t,t);
        opts.isolevel = .5;
        clf; plot_isosurface(f{i},opts);
        a=1;
    end
end

opts.alpha = 1; % transparency
opts.color = [0 1 0];
opts.isolevel = .5;
clf;
for i=1:p
    subplot(1,length(names),i);
    plot_isosurface(f{i},opts);
end

rep = ['../results/volumetric/' names{1}];
for i=2:p
    rep = [rep '-' names{i}];
end
rep = [rep '/'];

if not(exist(rep))
    mkdir(rep);
end

%%
% Compute transport coupling
%   pi = diag(w1)*K*diag(w0)
% and
%   pi*1 = w1.*K(w0) = p1
% and
%   pi'*1 = w1.*K(w0) = p0

options.tol = 1e-9;
options.tol = 0;
options.niter = 100;
options.verb = 2;
slicing = @(x)x(:,:,end/2);
options.disp = @(w0,w1)imageplot( slicing(w0.*K(w1)) );

% clf;
% [distances,w0,w1] = convolutionalDistance(f{1}, f{2}, [], K,[], options);

%%
% Compute displacement interpolation

q = 5;
switch p
    case 2
        % displacement interpolation
        t = linspace(0,1,q);
        W = [t;1-t];
    case 3
        % triangle interpolation
        W = [ ...
            [0, 0, 1]; ...
            [1, 0, 3]; [0, 1, 3]; ...
            [1,0,1]; [1,1,2]; [0,1,1]; ...
            [3,0,1]; [2,1,1]; [1,2,1]; [0,3,1]; ...
            [1,0,0]; [3,1,0]; [1,1,0]; [1,3,0]; [0,1,0] ...
            ]';
    case 4
        % bilinear interpolation
        t = linspace(0,1,q);
        [T,S] = meshgrid(t,t); S = S(:); T = T(:);
        W = [(1-S).*(1-T) S.*(1-T) (1-S).*T S.*T]';
end
Q = size(W,2);

%%
% Do the computations.

for i=1:Q
    w = W(:,i)'; w = w/sum(w);
    % select entropy bound
	entropyLimit = [];
    % store as 2D matrix
    Hv = [];
    for k = 1:p
        Hv = [Hv f{k}(:)/sum(f{k}(:))];
    end
    % do the computation
    options.disp = @(x)plot_isosurface(reshape(x,[N N N]));
    [B,u] = convolutionalBarycenter( Hv, w, [], Kv, [],entropyLimit, options);
    B = reshape(B, [N N N]);
    B = B/max(B(:));
    Bsvg{i} = B;
    if isnan(max(B(:)))
        warning('NaN problem.');
    end
end

%%
% Do the rendering.

col = [ [1 0 0]; [0 1 0]; [0 0 1]; [1 1 0]; [1/2 1 1/2]; [1/2 1/2 1] ];

for i=1:Q
    w = W(:,i)'; w = w/sum(w);
    B = Bsvg{i};
    % display
    clf;
    opts.alpha = 1; % transparency
	opts.color = sum(col(1:p,:) .* repmat(w(:),[1 3]) );
    opts.isolevel = median(B(:));
    opts.isolevel = (max(B(:))-min(B(:)))/2;
    F = plot_isosurface(B,opts);
    % save as image
    if p==4
        str = ['barycenter-' num2str(S(i)*(q-1)) '-' num2str(T(i)*(q-1))];
    else
        str = ['barycenter-' num2str(i)];
    end
    saveas(gcf, [rep str '.png'], 'png');
    % save as object
    if 0
    mesh.vertices = F.vertices;
    mesh.triangles = F.faces;
    mesh.numVertices = size(F.vertices, 1);
    vertexTexture = zeros(mesh.numVertices,1);
    writeTexturedObj([rep 'barycenter-' str '.obj'], mesh, vertexTexture);
    end
end

author: gpeyre


text:
function []=ps_calc_tca
%PS_CALC_TCA calculate topo-correlated atmosphere
%   execute in according INSAR (for PS) or SMALL_BASELINES (for SB) ifg
%   correction after running stamps(1,7)
%
%   The script will load already existing correction
%   Keep in mind that it is better to apply no
%   correction if the correlation is low !
%
%   Hannes Bathke, November 2012
%
%   ================================================================
%   11/2012 HB added more precise pixel selection, bugfixes
%   11/2012 HB drop ifg index included
%   01/2013 AH converted to function and name changed
%   08/2014 DB Append the correction in case the tca variable esits
%   08/2014 DB Change such can be ran without having ran step 7
%   ================================================================

clear all
close all

plot_flag = 1;    % 0 --> no plot/ 1 -> plot uph, correction, corrected ifg after each correction
%% load standard ps indata check if sbas or ps folder
% load interferograms

load psver
psname=['./ps',num2str(psver)];
sclaname=['./scla',num2str(psver) '.mat'];
sclasbname=['./scla_sb',num2str(psver) '.mat'];
phuwname=['./phuw',num2str(psver) '.mat'];
phuwsbname=['./phuw_sb',num2str(psver) '.mat'];

ps=load(psname);

num_ifg = ps.n_ifg;
num_pix = ps.n_ps;

small_baseline_flag = getparm('small_baseline_flag');
if ~strcmpi(small_baseline_flag,'y')
    disp('working on the PS interferograms');
    if exist([phuwname],'file')
        stratname=['./tca',num2str(psver)];
        u_ph = load(phuwname');

    end
    if exist(sclasbname,'file')==2
        scla=load(sclaname);
        % subtract dem error
        ph_all=u_ph.ph_uw - scla.ph_scla;

        fprintf('DEM error removed \n')
    else
        ph_all=u_ph.ph_uw  ;
    end
else
    disp('working on the SB interferograms');
    if exist([phuwsbname],'file')
        stratname=['./tca_sb',num2str(psver)];
        u_ph = load(phuwsbname);
     end
    if exist(sclasbname,'file')==2
        scla=load(sclasbname);
        % subtract dem error
        ph_all=u_ph.ph_uw - scla.ph_scla;
        fprintf('DEM error removed \n')

    else
        ph_all=u_ph.ph_uw;
    end
end

% if existent load strat correction
strat_corr = zeros(num_pix,num_ifg);

% load dem data
if exist('hgt2.mat','file')
    dem = load('hgt2.mat');
else
    disp('dem height not existing yet, please run StaMPS(1,7) once')
    exit
end

% select interferograms for correction
drop_ifg_index = getparm('drop_ifg_index');
disp(['dropped ifgs: ',num2str(drop_ifg_index)]);
disp(['Loop through ',num2str(num_ifg),' ifgs:']);
ifg = input('Start with which ifg? ');
end_ifg = input('End with which ifg? ');

for k = ifg : end_ifg
    disp(['interferogram number ',num2str(k)])
    % check if ifg is on drop index
    if k == drop_ifg_index;
        disp('ifg dropped, nothing done, next');
    else
        figure;
        plot(dem.hgt,ph_all(:,k),'.k')
        hold on
        xlabel('DEM Elevation [m]');
        ylabel('unwrapped phase')
        do_corr = input('correlation existent? (yes - 1; no(next ifg) - 0 ) ');
        if do_corr == 1
            disp('Click 2 edges of a rectangle within the linear function\n should be fitted to the pixels ');
            bounds = ginput(2);
            low_hgt = min(bounds(:,1));
            high_hgt = max(bounds(:,1));
            low_uph = min(bounds(:,2));
            high_uph = max(bounds(:,2));
            %find pixels within bounds
            did = find(low_hgt < dem.hgt & dem.hgt < high_hgt);
            phid = find(low_uph < ph_all(:,k) & ph_all(:,k) < high_uph);
            pix_int = intersect(did,phid);
            plot(dem.hgt(pix_int),ph_all(pix_int,k),'.b')

            % build design matrix
            dem_range = low_hgt:1:high_hgt;
            A = [dem.hgt(pix_int) ones(length(pix_int),1)];
            Afull = [dem_range' ones(length(dem_range),1)];

            % do the inversion
            coef = (A'*A)\A'*ph_all(pix_int,k);
            strat_fun = Afull * coef;

            % display fitted function
            plot(dem_range,strat_fun,'-r')
            hold off
            strat_corr(:,k) = coef(1)*dem.hgt;

            if (plot_flag)
                % plot results
                figure; scatter(ps.lonlat(:,1),ps.lonlat(:,2),5,strat_corr(:,k),'filled');
                title('estimated stratified atm comp'); colorbar
                caxis([min(ph_all(:,k)) max(ph_all(:,k))])

                figure; scatter(ps.lonlat(:,1),ps.lonlat(:,2),5,ph_all(:,k),'filled');
                title(['unwrapped phase of ifg ',num2str(k)]); colorbar
                caxis([min(ph_all(:,k)) max(ph_all(:,k))])

                figure; scatter(ps.lonlat(:,1),ps.lonlat(:,2),5,...
                    (ph_all(:,k)-strat_corr(:,k)),'filled');
                title(['corrected ifg ',num2str(k)]); colorbar
                caxis([min(ph_all(:,k)) max(ph_all(:,k))])

                next = input('next ifg? (1 --> yes) (0 --> EXIT) ');
                if next == 1
                    close all
                else
                    break
                end
            end
            close all
        else
            disp(['no correlation in interferogram ',num2str(k),' no correction applied--> 0']);
            close all
        end
    end
end

% save matrix
if exist([stratname,'.mat'],'file')~=2
    save([stratname,'.mat'],'strat_corr')
else
    save([stratname,'.mat'],'-append','strat_corr')
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function r1 = rOneCalculation(fracFreq)
%method for caculating r1 used in identifying noise type for a giving tau
%value. input argument is array of fractrional frequency values for a given
%tau.
%r1 = <(Zi - Zavg)*(Zi+1 - Zavg)> / <(Zi - Zavg)^2>
%Z is fractional freq values. Numerator sum is from 1 to N-1 and demonator
%is from 1 to N

avgFF = mean(fracFreq); %get average value
N = length(fracFreq); %get number of frac freq value in array
n=0; %store numerator of r1
j=0; %store denomenator of r1

%This loop builds the numerator for r1
for i = 1:(N-1)
    if i == 1
        n = (fracFreq(i) - avgFF)*(fracFreq(i+1) - avgFF);
    else
        n = n + (fracFreq(i) - avgFF)*(fracFreq(i+1) - avgFF);
    end
end
%this loop builds the denomentor for r1
for i = 1:N
    if i == 1
        j = (fracFreq(i) - avgFF)^2;
    else
        j = j + (fracFreq(i) - avgFF)^2;
    end
end

%calculate and return r1
r1 = n/j;


text:
function [A,uE2F,uE] = facet_adjacency_matrix(F,varargin)
  % FACET_ADJACENCY_MATRIX  Adjacency matrix between facets determined by
  % whether two facets share an edge.
  %
  % A = facet_adjacency_matrix(F)
  % A = facet_adjacency_matrix(F,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   F  #F by 3 list of triangles
  %   Optional:
  %     'ManifoldOnly' followed by whether to only consider adjacency across
  %     manifold edges (valence <=2) {false}
  % Outputs:
  %   A  #F by #F adjacency matrix
  %   uE2F  #E by #F matrix so that (e,f) = 1 means face f is adjacent to
  %     unique edge e
  %   uE  #E by 2 list of unique edges
  %
  % See also: adjacency_matrix

  manifold_only = false;
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'ManifoldOnly'},{'manifold_only'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  ss = size(F,2);
  switch ss
  case 3
    % List of all "half"-edges: 3*#F by 2
    allE = [F(:,[2 3]); F(:,[3 1]); F(:,[1 2])];
    % Sort each row
    sortallE = sort(allE,2);
    % IC(i) tells us where to find sortallE(i,:) in uE:
    % so that sortallE(i,:) = uE(IC(i),:)
    [uE,~,IC] = unique(sortallE,'rows');
    % uE2F(e,f) = 1 means face f is adjacent to unique edge e
    uE2F = sparse(IC(:),repmat(1:size(F,1),1,ss)',1);
  case 2
    % We're really dealing with edges, so E-->vertices, F-->edges
    uE = 1:max(F(:));
    uE2F = sparse(F,repmat(1:size(F,1),2,1)',1);
  end

  % kill non-manifold edges
  if manifold_only
    uE2F(sum(uE2F,2)>2,:) = 0;
  end
  % Face-face Adjacency matrix
  A = uE2F'*uE2F;
  % All ones
  A = A>0;
end

author: alecjacobson

repo: gptoolbox


text:
function [G1, impact] = sparse_solver(g1,psi,n_v,n_g)
%% NOTE: THIS SOLVER WILL BE UPDATED TO MATCH SCHUR_SOLVER EVENTUALLY
% Solves the rational expectation model with sparse iteration
%    method. You should use this method of the number of stable and
%    unstable solution differ a lot. Otherwise, full matrix
%    decomposition method should be used.
%
% by SeHyoun Ahn, June 2016
%
% PARAMETERS/OUTPUTS:
%     Read attached documentation
%
% SYNTAX:
% [G1, impact] = sparse_solver(g1,psi,n_v,n_g)

%-100 is set to find negative eigenvalues. This value can be changed
[x,v,flag] = eigs(g1,n_g,-100);
impact = real(x*(x(n_v+1:end,:)\psi(n_v+1:end,:)));
G1 = real(x*v*(x(n_v+1:end,:)\[sparse(n_g,n_v),speye(n_g,n_g)]));

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
% 3D rotatin about Y
function rot_y = rotMatY_3D(ty)

	rot_y = [cos(ty)	 0.0		sin(ty);
	         0.0		 	 1.0		0.0;
		    -sin(ty)	 0.0		cos(ty) ];

end

repo: MOTBeyondPixels


text:
function fmincon_test02 ( )

%*****************************************************************************80
%
%% FMINCON_TEST02 minimizes the Niederreiter-McCurley function 2.
%
%  Discussion:
%
%    The value determined by FMINCON seems highly sensitive to the
%    starting condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FMINCON_TEST02\n' );
  fprintf ( 1, '  Minimize the negative Niederreiter-McCurley #2 function\n' );

  n = 4;
  lb = [ 0, 0, 0, 0 ];
  ub = [ 1, 1, 1, 1 ];
  x0 = ( lb + ub ) / 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial data:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  ---X(%d)---', j );
  end
  fprintf ( 1, '  ----F(X)----\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %10f', x0(j) );
  end
  fprintf ( 1, '  %12e\n', - niederreiter_mccurley2 ( x0 ) );

  [ x, fval, exitflag, output ] = fmincon ( @niederreiter_mccurley2, x0, ...
    [], [], [], [], lb, ub );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed optimum after %d function evaluations:\n', ...
    output.funcCount );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  ---X(%d)---', j );
  end
  fprintf ( 1, '  ----F(X)----\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %10f', x(j) );
  end
  fprintf ( 1, '  %12e\n', - niederreiter_mccurley2 ( x ) );

  return
end


text:
% objdall1n - objective function of DAL with non-negative L1 regularization
%
% Copyright(c) 2009-2011 Ryota Tomioka
%                   2011 Shigeyuki Oba
% This software is distributed under the MIT license. See license.txt

function varargout=objdall1n(aa, info, prob, ww, uu, A, B, lambda, eta)

m = length(aa);
n = length(ww);

if isempty(info.ATaa)
  info.ATaa=A.Ttimes(aa);
end

vv = ww+eta(1)*info.ATaa;

if nargout<=3
  [floss, gloss, hmin]=prob.floss.d(aa, prob.floss.args{:});
else
  [floss, gloss, hloss, hmin]=prob.floss.d(aa, prob.floss.args{:});
end

[vsth,ss] = l1n_softth(vv,eta(1)*lambda);
info.wnew = vsth;
info.spec = ss;

fval = floss+0.5*sum(vsth.^2)/eta(1);
if ~isempty(uu)
  u1   = uu+eta(2)*(B'*aa);
  fval = fval + 0.5*sum(u1.^2)/eta(2);
end

varargout{1}=fval;

if nargout<=2
  varargout{2} = info;
else
  gg  = gloss+A.times(vsth);
  soc = sum(((vsth-ww)/eta(1)).^2);
  if ~isempty(uu)
    gg  = gg+B*u1;
    soc = soc+sum((B'*aa).^2);
  end

  if soc>0
    info.ginfo = norm(gg)/(sqrt(min(eta)*hmin*soc));
  else
    info.ginfo = inf;
  end
  varargout{2} = gg;

  if nargout==3
    varargout{3} = info;
  else
    I = find(vv>eta(1)*lambda);
    len = length(I);
    AF = A.slice(I);

    switch(info.solver)
     case 'cg'
      prec=hloss+spdiag(eta(1)*sum(AF.^2,2));
      if ~isempty(uu)
        prec =prec+spdiag(eta(2)*sum(B.^2,2));
      end
      varargout{3} = struct('hloss',hloss,'AF',AF,'I',I,'n',n,'prec',prec,'B',B);
     otherwise
      if length(I)>0
        varargout{3} = hloss+eta(1)*AF*AF';
      else
        varargout{3} = hloss;
      end
      if ~isempty(uu)
        varargout{3} = varargout{3}+eta(2)*B*B';
      end
    end
    varargout{4} = info;
  end
end

repo: REST


text:
function Y = cosh(X)
    % Symbolic hyperbolic cosine.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % construct the operation string
    sstr = ['Cosh[' X.s ']'];

    % create a new object with the evaluated string
    Y = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

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


text:
function fixed_num = p11_fixed_num ( )

%*****************************************************************************80
%
%% P11_FIXED_NUM returns the number of fixed points in problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer FIXED_NUM, the number of fixed points.
%
  fixed_num = 6;

  return
end


text:
function[varargout]=cellgrid(varargin)
%CELLGRID  Interpolate a cell array of numeric arrays onto a regular grid.
%
%   [TO,Y]=CELLGRID(T,X,DT) where T is a cell array of time arrays, and
%   X is of the same size as T, linearly interpolates the elements of X
%   within T at times TO, which are regularly spaced with interval DT.
%
%   DT may either be a scalar, or an array of the same length as T and X.
%
%   CELLGRID does not modify bad data points, such as those marked with
%   NaNs.  See CELLFILL to interpolate over bad data points.
%
%   As an example, with T and X given by
%
%       T{1} = [1 2 3  5  6]';   T{2}=[3 7  9 10]';
%       X{1} = [2 4 6 10 12]';   X{2}=[5 9 11 12]';
%
%   [TO,Y]=CELLGRID(T,X,1) will return
%
%       TO{1} = [1 2 3 4  5  6]';   TO{2} = [3 4 5 6 7 8  9  10]';
%       Y{1}  = [2 4 6 8 10 12]';   Y{2}  = [5 6 7 8 9 10 11 12]';
%
%   By default, CELLGRID uses INTERP1 with the 'pchip' method of
%   interpolation.  CELLGRID(...,STR) instead uses the method specified by
%   STR, e.g. STR='linear'.  See INTERP1 for details.
%
%   [TO,Y1,Y2,...YN]=CELLGRID(T,X1,X2,...XN,DT) with multiple input
%   arguments also works provided the XN are all the same size.
%
%   CELLGRID(X1,X2,...XN); with no output arguments overwrites the
%   original input variables.
%   __________________________________________________________________
%
%   Specifying interpolated times
%
%   [TO,Y]=CELLGRID(T,X,DT,A,B) will set TO to TO=A:DT:B, where A and B
%   are either scalars or an arrays of the same size as X.
%
%   The default behavior is equivalent to choosing A and B as the first and
%   last elements of TO, that is, to setting A=CELLMIN(T) and B=CELLMAX(T).
%   __________________________________________________________________
%
%   Parallelization
%
%   CELLGRID(...,'parallel') parallelizes the computation using a PARFOR
%   loop over the various input variables.  This requires that Matlab's
%   Parallel Computing Toolbox be installed.
%   __________________________________________________________________
%
%   'cellgrid --t' runs a test.
%
%   Usage: [to,y]=cellgrid(t,x,dt);
%          [to,y1,y2,y3]=cellgrid(t,x1,x2,x3,dt);
%          [to,y1,y2,y3]=cellgrid(t,x1,x2,x3,dt,a,b);
%          cellgrid(t,x1,x2,x3,dt);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015--2019 J.M. Lilly --- type 'help jlab_license' for details

if strcmp(varargin{1}, '--t')
    cellgrid_test,return
end

str='pchip';
cores='serial';

for i=1:2
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'ser')||strcmpi(varargin{end}(1:3),'par')
            cores=varargin{end};
        else
            str=varargin{end};
        end
        varargin=varargin(1:end-1);
    end
end

if strcmpi(cores(1:3),'par')
    if exist('parpool')~=2
        disp('Sorry, parallel algorithm requires the Parallel Computing Toolbox.')
        disp('Defaulting to the serial algorithm.')
        str='serial';
    end
end

t=varargin{1};

if ~iscell(varargin{end-2})&&~ischar(varargin{end-2})
    a=varargin{end-1};
    b=varargin{end};
    if length(a)==1
        a=a+zeros(size(t));
    end
    if length(b)==1
        b=b+zeros(size(t));
    end
    varargin=varargin(1:end-2);
else
    a=cellmin(t);
    b=cellmax(t);
end

dt=varargin{end};

varargin=varargin(1:end-1);

%Size check
for i=2:length(varargin)
    if size(varargin{i})~=size(varargin{1});
        error('All input arguments must be the same size.')
    end
end

if length(dt)==1
    dt=dt+0*a;
end

for i=1:length(t)
   to{i,1}=[a(i):dt(i):b(i)]';
end

if strcmpi(cores(1:3),'par')
    disp('CELLGRID employing parallel algorithm.')
    parfor j=2:length(varargin)
        for i=1:length(t)
            if length(t{i})>1
                bool=~isnan(t{i})&~isnan(varargin{j}{i});
                if length(find(bool))~=0
                    args{j}{i,1}=interp1(t{i}(bool),varargin{j}{i}(bool),to{i},str);
                else
                    args{j}{i,1}=nan*to{i};
                end
            else
                args{j}{i,1}=varargin{j}{i};
            end
        end
    end
else
    for j=2:length(varargin)
        for i=1:length(t)
            if length(t{i})>1
                bool=~isnan(t{i})&~isnan(varargin{j}{i});
                if length(find(bool))>=2
                    %i,j,length(find(bool))
                    args{j}{i,1}=interp1(t{i}(bool),varargin{j}{i}(bool),to{i},str);
                else
                    args{j}{i,1}=nan*to{i};
                end
            else
                args{j}{i,1}=varargin{j}{i};
            end
        end
    end
end

varargout=args;
if ~isempty(varargout{2})
    for i=1:length(to)
        if isempty(varargout{2}{i})
            to{i}=[];
        end
    end
else
    to=[];
end
varargout{1}=to;

eval(to_overwrite(length(varargin)))

function[]=cellgrid_test

t{1}=[1 2 3  5  6]';
t{2}=[3 7  9 10]';
x{1}=[2 4 6 10 12]';
x{2}=[5 9 11 12]';

to{1}=[1 2 3 4  5  6]';
to{2}=[3 4 5 6 7 8  9  10]';
y{1}=[2 4 6 8 10 12]';
y{2}=[5 6 7 8 9 10 11 12]';

cellgrid(t,x,1);

reporttest('CELLGRID',aresame(x,y)&&aresame(to,t))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
classdef TestAdd
    %TestAdd

    methods (Static)
        function test_rgb_image
            img = cv.imread(fullfile(mexopencv.root(),'test','img001.jpg'), 'ReduceScale',2);
            [h,w,~] = size(img);
            mask = false(h,w);
            mask(100:h-100,100:w-100) = true;

            out = cv.add(img, img);
            validateattributes(out, {class(img)}, {'size',size(img)});
            if mexopencv.require('images')
                expected = my_add(img, img);
                assert(isequal(out, expected));
            end

            out = cv.add(img, img, 'DType','double');
            validateattributes(out, {'double'}, {'size',size(img)});
            if mexopencv.require('images')
                expected = my_add(img, img, 'double');
                assert(isequal(out, expected));
            end

            out = cv.add(img, img, 'Mask',mask);
            validateattributes(out, {class(img)}, {'size',size(img)});
            if mexopencv.require('images')
                expected = my_add(img, img, class(img), mask);
                assert(isequal(out, expected));
            end

            out = cv.add(img, img, ...
                'Dest',img, 'Mask',mask, 'DType',class(img));
            validateattributes(out, {class(img)}, {'size',size(img)});
            if mexopencv.require('images')
                expected = my_add(img, img, class(img), mask, img);
                assert(isequal(out, expected));
            end
        end

        function test_2d_matrix
            M = uint8([100 150; 200 250]);
            mask = true(size(M));
            mask(4) = false;

            out = ones(size(M), 'uint8');
            out = cv.add(M, M, 'Mask',mask, 'Dest',out);
            validateattributes(out, {'uint8'}, {'size',size(M)});
            expected = uint8([200 255; 255 1]);
            assert(isequal(out, expected));

            out = cv.add(M, M, 'Mask',mask);
            validateattributes(out, {'uint8'}, {'size',size(M)});
            expected = uint8([200 255; 255 0]);
            assert(isequal(out, expected));

            out = cv.add(M, M);
            validateattributes(out, {'uint8'}, {'size',size(M)});
            expected = M + M;
            assert(isequal(out, expected));

            out = cv.add(M, 1);
            validateattributes(out, {'uint8'}, {'size',size(M)});
            expected = M + 1;
            assert(isequal(out, expected));

            out = cv.add(1, M);
            validateattributes(out, {'uint8'}, {'size',size(M)});
            expected = 1 + M;
            assert(isequal(out, expected));

            out = cv.add(1, 1);
            validateattributes(out, {'double'}, {'scalar'});
            expected = 1 + 1;
            assert(isequal(out, expected));
        end

        function test_output_depth
            out = cv.add(uint16(1), int8(1), 'DType','single');
            validateattributes(out, {'single'}, {'scalar'});
            assert(isequal(out, 2));
        end

        function test_saturation
            out = cv.add(uint8(200), uint8(200), 'DType','uint8');
            validateattributes(out, {'uint8'}, {'scalar'});
            assert(isequal(out, 255));

            out = cv.add(uint8(200), uint8(200), 'DType','uint16');
            validateattributes(out, {'uint16'}, {'scalar'});
            assert(isequal(out, 400));

            % saturation not applied for int32 (result overflows)
            out = cv.add(int32(1073741824), int32(1073741824), 'DType','int32');
            %validateattributes(out, {'int32'}, {'scalar', '<',0});

            % unlike MATLAB where result saturates
            %out = int32(1073741824) + int32(1073741824);
            %assert(isequal(out, intmax('int32')));
        end

        function test_error_argnum
            try
                cv.add();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function out = my_add(src1, src2, dtype, mask, dst)
    %MY_ADD  Similar to cv.add using imadd from IPT

    if nargin < 3, dtype = class(src1); end
    if nargin < 4, mask = true(size(src1,1),size(src1,2)); end
    if nargin < 5, dst = zeros(size(src1), dtype); end

    % add two images with specified output class
    out = imadd(src1, src2, dtype);

    % apply masking
    if ~isempty(mask)
        for k=1:size(out,3)
            out_slice = out(:,:,k);
            dst_slice = dst(:,:,k);
            out_slice(~mask) = dst_slice(~mask);
            out(:,:,k) = out_slice;
        end
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [fBValue, fStdDev, fMc, fAValue, nNumberQuakes] = calc_BandMc(mCatalog, nMinimumNumber, mcCalcMethod, fBinning, bConstrainMc, fMcMin, fMcMax, fMcAdd)
%  b-values (with stddev) and the magnitudes of completeness of a given catalog
%
  % [fBValue, fStdDev, fMc, fAValue, nNumberQuakes] = calc_BandMc(mCatalog, nMinimumNumber, mcCalcMethod, fBinning, bConstrainMc, fMcMin, fMcMax, fMcAdd)
% --------------------------------------------------------------------------------------------------------------------------------------------------------------
% Calculation of the b-values, its standard deviation and the magnitudes of completeness of a given catalog
%
% Input parameters:
%   mCatalog            Earthquake catalog
%   nMinimumNumber      Minimum number of earthquakes in the catalog for calculating the output values
%   mcCalcMethod        Method to determine the magnitude of completeness (see also: McMethods)
%   fBinning            Magnitude binning of the catalog (default 0.1)
%   bConstrainMc        Constrain Mc to [fMcMin, fMcMax] if set to 1 (default 0)
%   fMcMin              see bConstrainMc
%   fMcMax              see bConstrainMc
%   fMcAdd              Value to be added to Mc after potential constraining and prior to
%                       b-value computation for conservative computations
%
% Output parameters:
%   fBValue             b-value of the catalog with respect to magnitude of completeness
%   fStdDev             Standard deviation of the b-value (Shi & Bolt)
%   fMc                 Magnitude of completeness (Mc)
%   fAValue            a-value of the catalog
%   nNumberQuakes       Number of quakes used to compute the b-value (M > Mc)
%
% Copyright (C) 2003 by Danijel Schorlemmer
%
% 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 2 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, write to the
% Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

report_this_filefun();

% Magnitude binning
if ~exist('fBinning', 'var')
  fBinning = 0.1;
end

% Constrain magnitude of completeness
if ~exist('bConstrainMc', 'var')
  bConstrainMc = false;
end
if ~exist('fMcMin', 'var')
  fMcMin = 0;
end
if ~exist('fMcMax', 'var')
  fMcMax = 1;
end

% Conservative adding of a constant
if ~exist('fMcAdd', 'var')
  fMcAdd = 0;
end

% Init output variables
fBValue = nan;
fStdDev = nan;
fMc = nan;
fAValue = nan;

try
  % Determine magnitude of completeness
  fMc = calc_Mc(mCatalog, mcCalcMethod, fBinning);
  % Constrain magnitude of completeness
  if bConstrainMc
    if fMc < fMcMin
      fMc = fMcMin;
    elseif fMc > fMcMax
      fMc = fMcMax;
    end
  end
  % Conservative Mc
  fMc = fMc + fMcAdd;
  % Calculate the b-value of the learning period
  vSel_ = mCatalog.Magnitude >= fMc;
  mCatalog = mCatalog.subset(vSel_);
  nNumberQuakes = mCatalog.Count;
  if nNumberQuakes >= nMinimumNumber
    [fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog.Magnitude, fBinning);
  end
  if isempty(fMc)
    fMc = nan;
  end
catch
end

author: CelsoReyes


text:
function str=disp(tt,name)
%Command window display of a TT-tensor.
%
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if ~exist('name','var')
    name = 'ans';
end
r=tt.r; n=tt.n;
d=tt.d;
str=[];
str=[str,sprintf('%s is a %d-dimensional TT-tensor, ranks and mode sizes: \n',name,d)];
%fprintf('r(1)=%d \n', r(1));
for i=1:d
    str=[str,sprintf('r(%d)=%d \t n(%d)=%d \n',i,r(i),i,n(i))];
%   fprintf(' \t n(%d)=%d \n',i,n(i));
%   fprintf('r(%d)=%d \n',i+1,r(i+1));
end
 str=[str,sprintf('r(%d)=%d \n',d+1,r(d+1))];
 if ( nargout == 0 )
    fprintf(str);
 end
end

author: oseledets

repo: TT-Toolbox