text:
function nVars = numVars(N)
%NUMVARS   Detect how many unknown variables a CHEBOP operates on.
%   NVARS = NUMVARS(N), where N is a CHEBOP, returns the integer NVARS which is
%   equal to the number of variables that N operates on.

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

% NUMVARS indicate how many unknown function we seek. This can be tricky of the
% operator is specified on CHEBMATRIX syntax, e.g. via
%   N.op = @(x, u) [diff(u{1}) + u{2}; u{1} + diff(u{2})];
% is in this case, nargin(N.op) does not match the number of variables we need.
% If nargin(N.op) is greater than two, we can however safely assume that N.op is
% specified on the form
%   N.op = @(x, u, v) [diff(u) + v; u + diff(v)];

% We begin by looking at whether we have the easy case!
narginN = nargin(N);
if ( narginN > 2 )
    % Need to subtract 1 since x is the first argument:
    nVars = narginN - 1;

elseif ( ~isempty(N.numVars) )
    % Now we know we're dealing with
    %   N.op = @(x, u) ...
    % But we don't know yet whether this is a system or not. In an ideal world,
    % the user has passed this information through the numVars property of the
    % CHEBOP.

    nVars = N.numVars; % Lucky us!

else
    % If that field is empty, we try to inspect the string representation of
    % N.op, and look for the highest index appearing inside {} -- this will only
    % work if N.op is specified as an anonymous function, not if it is a handle
    % to another function.

    % Obtain the function string:
    NopString = func2str(N.op);

    % Try to find out what the unknown function is. But first, we must check
    % whether we actually have the argument list available to us...
    firstRightPar = min(strfind(NopString, ')'));
    if ( isempty(firstRightPar) )
        % Don't have a list of arguments available. Take numVars == 1, and hope
        % for the best...
        warning('CHEBFUN:CHEBOP:numVars:numberOfArguments', ...
            ['Unable to determine the number of variable that the ', ...
            'chebop\noperates on. Assuming problem is a scalar problem,',...
            ' results might be\nunreliable. Please specify the number', ...
            ' of variables that the operator\n', ...
            'operates on via CHEBOP.numVars.'])
        nVars = 1;

    else
        % If nargin(N) == 2, the first variable appearing will be the
        % independent space variable, but if nargin(N) == 1, we only have the
        % unknown function in the argument list.
        if ( narginN == 1 )
            % The unknown variable appears between @( and the first ).
            firstLeftPar = min(strfind(NopString, '('));
            variableName = NopString(firstLeftPar+1 : firstRightPar - 1);

        else
            % The unknown variable appears between the first , and the first ),
            % e.g. @(x, ___ ).
            firstLeftComma = min(strfind(NopString, ','));
            variableName = NopString(firstLeftComma+1 : firstRightPar - 1);
        end

        % The regular expression we seek must include the variable name:
        expression = [variableName, '{[1-9]+}'];

        % Obtain all matches:
        match = regexp(NopString, expression, 'match');

        % If we don't have any matches, problem can't have been specified using
        % CHEBMATRIX syntax. But nargin(N) == 2, and we did indeed have the
        % argument list to the anonymous function passed, so we must be dealing
        % with a scalar problem!
        if ( isempty(match) )
            nVars = 1;
        else
            % Throw away the variable name and the { }, e.g. convert 'u{1}' to
            % '1':
            match = strrep(match, [variableName, '{'], '');
            match = strrep(match, '}', '');

            % We are now left with cell-array of strings that only contain
            % numbers. So convert to doubles!
            indx = str2double(match);

            % The number of variables that the CHEBOP operates on is the
            % greatest index that appears:
            nVars = max(indx);
        end

    end

end


text:
function sF = power(sF1,sF2)
%
% Syntax
%   sF = sF1.^a
%

if isnumeric(sF1)
  f = @(v) sF1 .^ eval(sF2, v);

  bw = sF1 * sF2.bandwidth;
  sF = S2FunHarmonic.quadrature(f,'bandwidth',min(bw,256));
elseif isnumeric(sF2)
  f = @(v) eval(sF1, v) .^ sF2;

  bw = sF1.bandwidth * sF2;
  sF = S2FunHarmonic.quadrature(f,'bandwidth',min(bw,256));
else
  f = @(v) eval(sF1, v) .^ eval(sF2, v);

  bw = max(sF1.bandwidth, sF2.bandwidth);
  sF = S2FunHarmonic.quadrature(f,'bandwidth',min(bw,256));
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverse dynamics for the 2dof planar robot
%
%   tau = inversedynamics_2dofplanar(robot, q, qd, qdd, gc, fext)
%
%   Where robot stores the kinematic and dynamic parameters for this robot.
%   q: joint positions.
%   qd: joint velocities.
%   qdd: joint accelerations.
%   fext: vector of external forces. Defined in the BASE reference system.
%   let gc=9.81 m/s^2 if the force is acting on the -Y0 direction, or gc =
%   -9.81 m/s^2 if acting on the opposite direction.
%
%   This function just executes the inverse dynamic model for this robot.
%   The equations to compute this dynamic model can be found in:
%   "ROBOT ANALYSIS. The mechanics of Serial and Parallel
%        manipulators". Lung Weng Tsai. John Wiley and Sons, inc. ISBN:
%        0-471-32593-7. page 405.
%
%
%   Author: Arturo Gil Aparicio arturo.gil@umh.es
%   Date: 08/03/2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2016, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.
function qdd = forwarddynamics_2dofplanar(robot, q, qd, tau, gc, fext)
a = eval(robot.DH.a);
a1=a(1);
a2=a(2);

m1=robot.dynamics.masses(1);
m2=robot.dynamics.masses(2);

%Express it as a general dynamic function:
%M*qdd + V + G = Q, where Q is a vector of generalized forces or moments
%M is a 2x2 manipulator inertia matrix
M = [(1/3)*m1*a1^2 + m2*(a1^2 + a1*a2*cos(q(2))+(1/3)*a2^2)  m2*((1/2)*a1*a2*cos(q(2))+(1/3)*a2^2);
     m2*((1/2)*a1*a2*cos(q(2))+(1/3)*a2^2)                      (1/3)*m2*a2^2];

V = [-m2*a1*a2*sin(q(2))*(qd(1)*qd(2)+(1/2)*qd(2)^2);
            (1/2)*m2*a1*a2*sin(q(2))*qd(1)^2        ];

G = [(1/2)*m1*gc*a1*cos(q(1))+ m2*gc*a1*cos(q(1)) + (1/2)*m2*gc*a2*cos(q(1)+q(2));
         (1/2)*m2*gc*a2*cos(q(1)+q(2))              ];

%assure fext is a column vector
fext=fext(:);

%External forces are propagated to every joint by using the manipulators
%Jacobian
J = manipulator_jacobian(robot, q);

% %Account for friction by summing in tau.
if robot.dynamics.friction
    for j=1:robot.DOF,
        tau(j) = tau(j) - friction(robot, qd, j);
    end
end

inv(M)
tau
V
G
J'*fext
%Finally use the general equation to compute the torques, considering the
%external forces.
qdd =  inv(M)*(tau - V - G - J'*fext);

repo: ARTE


text:
function [sos] = sosmatrixineq(sos,fM,option)
% SOSMATRIXINEQ --- Creates a SOS constraint from a matrix inequality
% constraint
%
% [SOSP] = sosmatrixineq(SOSP,fM)
%
% SOSP is the sum of squares program.
% fM is a polynomial matrix used to generate the polynomial inequality y'fMy>0

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% 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, see <http://www.gnu.org/licenses/>.

% AP - 16/4/2013
% JA - 6/6/2013

if nargin<3
    option='quadraticMineq';%sets the default: asks for sos expression v'M(x)v
end

[n,m] = size(fM);

if n~=m
    disp('ERROR: Matrix fM in inequality fM>0 must be square.');
    return
end

if isfield(sos,'symvartable')
    % Original Code

    if strcmp(option,'quadraticMineq')

        %creates the vector of variables Mvar to generate the quadratic expression
        %M_var'*fM*Mvar
        if n>sos.varmat.count
            s1 = 'syms ';
            for i=sos.varmat.count+1:n
                s1 = strcat(s1,sprintf(' Mvar_%d',i));
            end
            %eval(expression) evaluates the MATLAB code in the string expression.
            eval(s1);

            Mvar = sym(zeros(n-sos.varmat.count,1));
            for i=sos.varmat.count+1:n
                Mvar(i) = eval(['Mvar_',int2str(i)]);
            end

            %updates the vartable in the sos program
            sos.varmat.symvartable = [sos.varmat.symvartable; Mvar];
            Mvarctable = sym2chartable(Mvar);
            if sos.varmat.count==0
                sos.varmat.vartable = [Mvarctable(1:end)];
            else
                sos.varmat.vartable = [sos.varmat.vartable(1:end-1),',',Mvarctable(2:end)];
            end
            sos.varmat.count = n;
        end

        varMconst = sos.varmat.symvartable(1:n);

        %create the sosconstraint using the sparse multipartite option since it is
        %homogeneous in varMconst
        sos = sosineq(sos,varMconst.'*fM*varMconst,'sparsemultipartite',{sos.symvartable,varMconst}); %GV&JA 6/12/2013

    elseif strcmp(option,'Mineq')
        sos = sosineq(sos,fM); %GV&JA 10/01/2013
    end
else
    if strcmp(option,'quadraticMineq')
    % Multipoly Code: PJS 9/9/2013
    if n>sos.varmat.count
        Mvar = polynomial(zeros(n-sos.varmat.count,1));
        for i=sos.varmat.count+1:n
            Mvar(i) = pvar(['Mvar_' int2str(i)]);
        end

        %updates the vartable in the sos program
        sos.varmat.vartable = [sos.varmat.vartable; Mvar];
        sos.varmat.count = n;
    end

    %create the sosconstraint using the sparse multipartite option since it is
    %homogeneous in varMconst
    sos = sosineq(sos,Mvar.'*fM*Mvar,'sparsemultipartite',{sos.vartable,Mvar}); % PJS 9/9/2013

    elseif strcmp(option,'Mineq')
        sos = sosineq(sos,fM); %GV&JA 10/01/2013
    end

end

author: yu-jiang

repo: radpbook

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


text:
function [Mate1,Mate2,Mate3] = Neighbor_Pairing_Strategy(MatingPop,Zmin)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Jiawei Yuan

    Objs = MatingPop.objs;
    [Num,M] = size(Objs);
    Objs = (Objs - repmat(Zmin,Num,1));
    Objs = Objs./repmat(sqrt(sum(Objs.^2,2)),1,M);

    CosV = Objs * Objs';
    CosV = CosV - 3*eye(Num,Num);

    [~,SInd] = sort(-CosV,2);

    Nr = 10;
    Neighbor = SInd(:,1:Nr);

    Mate1 = MatingPop;

    P = ones(Num,2);
    for i = 1:Num
        P(i,1:2) = Neighbor(i,randsample(Nr,2));
        if rand>0.7
            P(i,2) = randsample(Num,1);
        end
    end

    Mate2 = MatingPop(P(:,1));
    Mate3 = MatingPop(P(:,2));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function pass = test_carrier_C2(pref)
% Test Carrier equation.
%
% Toby Driscoll / Asgeir Birkisson, Jume 2014.
if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;
pref.bvpTol = tol;
dom = [-1 1];

%%

pref.discretization = @chebcolloc2;
N = chebop(@(x,u) 0.01.*diff(u,2)+2.*(1-x.^2).*u+u.^2-1, dom);
rhs = 0;
N.bc = @(x,u) [u(-1); u(1)];
x = chebfun(@(x) x, dom);
N.init =  2.*(x.^2-1).*(1-2./(1+20.*x.^2));

u = solvebvp(N, rhs, pref);

xx = (-1:.25:1)';
hiquality_ans = [
    0
    -1.487429807540814
    -1.785617248281071
    1.572366197526305
    -1.539652044363185
    1.572366197526230
    -1.785617248281089
    -1.487429807540795
    0
    ];

pass(1) = norm(u(xx) - hiquality_ans) < tol;
pass(2) = norm(u([-1 1])) < tol ;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% SegmentalCylinder
% Concrete subclass of <GenericCylinder.html |GenericCylinder|>
% representing a cylinder with a cross section of a circular segment.

%%% Description
% |SegmentalCylinder| represents the shape of a segmental cylinder.  Its
% cross section is a circular segment (region surrounded by an arc and the
% chord connecting its end points).  The axis of the cylinder should be
% aligned with one of the axes of the Cartesian coordinate system.

%%% Construction
%  shape = SegmentalCylinder(normal_axis, height, center, radius, theta, d_theta)
%  shape = SegmentalCylinder(normal_axis, height, center, radius, theta, d_theta, dl_max)
%
% *Input Arguments*
%
% * |normal_axis|: axis of the cylinder.  It should be one of |Axis.x|,
% |Axis.y|, |Axis.z|.
% * |height|: size of the cylinder along its axis.
% * |center|: center of the cylinder in the format of |[x y z]|.  For
% |normal_axis = Axis.z|, |(x, y)| is the coordinate of the center of the
% circle.
% * |radius|: radius of the circle
% * |theta|: beginning angle of the segment in radian
% * |d_theta|: angular width of the segment in radian between -|pi| and |pi|
% (the angle of a segment should not be reflex)
% * |dl_max|: maximum grid size allowed in the cylinder.  It can be either |[dx
% dy dz]| or a single real number |dl| for |dx = dy = dz|.  If unassigned,
% |dl_max = Inf| is used.

%%% Example
%   % Create an instance of SegmentalCylinder.
%   shape = SegmentalCylinder(Axis.z, 100, [0 0 50], 50, pi/6, pi/3);
%
%   % Use the constructed shape in maxwell_run().
%   [E, H] = maxwell_run({INITIAL ARGUMENTS}, 'OBJ', {'vacuum', 'none', 1.0}, shape, {REMAINING ARGUMENTS});

%%% See Also
% <SectoralCylinder.html |CircularCylinder|>, <CircularCylinder.html
% |CircularCylinder|>, <CircularShellCylinder.html
% |CircularShellCylinder|>, <EllipticCylinder.html |EllipticCylinder|>,
% <PolyognalCylinder.html |PolygonalCylinder|>, <Shape.html |Shape|>,
% <maxwell_run.html |maxwell_run|>

classdef SegmentalCylinder < GenericCylinder

	methods
        function this = SegmentalCylinder(normal_axis, height, center, radius, theta, d_theta, dl_max)
			chkarg(istypesizeof(normal_axis, 'Axis'), '"normal_axis" should be instance of Axis.');
			chkarg(istypesizeof(height, 'real') && height > 0, '"height" should be positive.');
			chkarg(istypesizeof(center, 'real', [1, Axis.count]), ...
				'"center" should be length-%d row vector with real elements.', Axis.count);
			chkarg(istypesizeof(radius, 'real') && radius > 0, '"radius" should be positive.');
			chkarg(istypesizeof(theta, 'real'), '"theta" should be real.');
			chkarg(istypesizeof(d_theta, 'real') && (d_theta < pi || d_theta > -pi), '"d_theta" should be real between -pi and pi.');

			sectoral = SectoralCylinder(normal_axis, height, center, radius, theta, d_theta, dl_max);

			thetas = [theta, theta + d_theta];

			mc_dir = mean([cos(thetas.') sin(thetas.')]);  % direction bisecting the segment
			mc_dir = mc_dir / norm(mc_dir);

			[h, v, n] = cycle(normal_axis);
			midpt = center([h, v]) + mean(radius .* [cos(thetas.') sin(thetas.')]);  % midpoint of chord

			% lsf2d() can handle rho = [p q] with column vectors p and q.  The
			% level set function is the one for a rectangle defined in the
			% (theta, radius) domain.
			function level = lsf2d(p, q)
				chkarg(istypeof(p, 'real'), '"p" should be array with real elements.');
				chkarg(istypeof(q, 'real'), '"q" should be array with real elements.');
				chkarg(isequal(size(p), size(q)), '"p" and "q" should have same size.');

				level1 = sectoral.lsf2d(p, q);  % level set function of sectoral cylinder

				lm = {p - midpt(Dir.h), q - midpt(Dir.v)};
				level2 = zeros(size(p));
				for d = Dir.elems
					level2 = level2 + lm{d} .* mc_dir(d);
				end

				level = min(level1, level2);  % intersection of regions defined by level1 and level2
			end

			lprim = cell(1, Axis.count);
			lprim{h} = radius * cos(thetas) + center(h);
			lprim{v} = radius * sin(thetas) + center(v);
			lprim{n} = [-height height]/2 + center(n);

			if sectoral.lsf_th(0) > 0
				lprim{h} = [lprim{h}, center(h) + radius];
			end
			if sectoral.lsf_th(pi/2) > 0
				lprim{v} = [lprim{v}, center(v) + radius];
			end
			if sectoral.lsf_th(pi) > 0
				lprim{h} = [lprim{h}, center(h) - radius];
			end
			if sectoral.lsf_th(3*pi/2) > 0
				lprim{v} = [lprim{v}, center(v) - radius];
			end

			if nargin < 7  % no dl_max
				super_args = {normal_axis, @lsf2d, lprim};
			else
				super_args = {normal_axis, @lsf2d, lprim, dl_max};
			end

			this = this@GenericCylinder(super_args{:});
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rfAmp,rfPhase,rfFreq,rfCoil,GzAmp,GyAmp,GxAmp,ADC,Ext,uts,ts,flags]=PSD_AFI
global VCtl
global VVar
CV1=1.5e-3;
CV10=0;
CV11=0;
CV12=0;
CV13=0;
CV14=0;
CV2=2e-3;
CV3=20e-3;
CV4=0;
CV5=0;
CV6=0;
CV7=0;
CV8=0;
CV9=0;
rfAmpAll=[];
rfPhaseAll=[];
rfFreqAll=[];
rfCoilAll=[];
GzAmpAll=[];
GyAmpAll=[];
GxAmpAll=[];
ADCAll=[];
ExtAll=[];
rfTimeAll=[];
GzTimeAll=[];
GyTimeAll=[];
GxTimeAll=[];
ADCTimeAll=[];
ExtTimeAll=[];
SEtAll=[];
uts=[];
ts=[];
flags=[];
if VCtl.PlotSeq == 1
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
ts = [ts tS tE];
end
ts = [0 max(ts)-min(ts)];
return;
end
%==============Pulses 1==============
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
rfTime=[];
GzTime=[];
GyTime=[];
GxTime=[];
ADCTime=[];
ExtTime=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if isempty(tS) | isempty(tE) | (tS>=tE)
error('SE setting is incorrect for Pulses 1!');
end
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{1};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=CV3;
p.Duplicates=2;
p.FA=VCtl.FlipAng;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=0.5e-3;
p.tStart=0.2e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp1,rfPhase1,rfFreq1,rfCoil1,rfTime1]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil1(1)
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
else
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
end
p=[];
%--------------------
p.DupSpacing=CV3;
p.Duplicates=2;
p.Gz1Sign=1;
p.Gz2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV2;
p.t1Start=CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GzAmp1,GzTime1]=GzCartesian(p);
GzAmp=[GzAmp GzAmp1];
GzTime=[GzTime GzTime1];
end
p=[];
%--------------------
p.DupSpacing=CV3;
p.Duplicates=2;
p.Gy1Sign=1;
p.Gy2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV2;
p.t1Start=CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GyAmp1,GyTime1]=GyCartesian(p);
GyAmp=[GyAmp GyAmp1];
GyTime=[GyTime GyTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gx1Sign=-1;
p.Gx2Sign=1;
p.Gx3Sign=0;
p.Notes='cartesian frequency';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1Start=CV1;
p.t2Middle=VCtl.ME_TEs(1);
p.t3Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GxAmp1,GxTime1]=GxCartesian(p);
GxAmp=[GxAmp GxAmp1];
GxTime=[GxTime GxTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gx1Sign=-1;
p.Gx2Sign=1;
p.Gx3Sign=0;
p.Notes='cartesian frequency';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1Start=CV3+CV1;
p.t2Middle=CV3+VCtl.ME_TEs(2);
p.t3Start=CV3+VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GxAmp2,GxTime2]=GxCartesian(p);
GxAmp=[GxAmp GxAmp2];
GxTime=[GxTime GxTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Notes='cartesian readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tMiddle=VCtl.ME_TEs(1);
if strcmp(p.Switch,'on')
[ADC1,ADCTime1]=ADCCartesian(p);
ADC=[ADC ADC1];
ADCTime=[ADCTime ADCTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Notes='cartesian readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tMiddle=CV3+VCtl.ME_TEs(2);
if strcmp(p.Switch,'on')
[ADC2,ADCTime2]=ADCCartesian(p);
ADC=[ADC ADC2];
ADCTime=[ADCTime ADCTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=5;
p.Notes='calculate remaining scan time';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=0;
if strcmp(p.Switch,'on')
[Ext1,ExtTime1]=ExtBit(p);
Ext=[Ext Ext1];
ExtTime=[ExtTime ExtTime1];
end
p=[];
%--------------------
p.DupSpacing=CV3;
p.Duplicates=2;
p.Ext=1;
p.Notes='reset K space location';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV1*9/10;
if strcmp(p.Switch,'on')
[Ext2,ExtTime2]=ExtBit(p);
Ext=[Ext Ext2];
ExtTime=[ExtTime ExtTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV3*(9/10);
if strcmp(p.Switch,'on')
[Ext3,ExtTime3]=ExtBit(p);
Ext=[Ext Ext3];
ExtTime=[ExtTime ExtTime3];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=VCtl.TR*(9/10);
if strcmp(p.Switch,'on')
[Ext4,ExtTime4]=ExtBit(p);
Ext=[Ext Ext4];
ExtTime=[ExtTime ExtTime4];
end
p=[];
%--------------------
SEt=[tS tE];
rfAmp(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfPhase(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfFreq(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfCoil(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzAmp(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyAmp(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxAmp(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADC(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
Ext(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfTime(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzTime(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyTime(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxTime(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADCTime(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
ExtTime(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfAmp(abs(rfAmp)<eps) = 0;
rfTime = rfTime + SEt(1);
GzTime = GzTime + SEt(1);
GyTime = GyTime + SEt(1);
GxTime = GxTime + SEt(1);
ADCTime = ADCTime + SEt(1);
ExtTime = ExtTime + SEt(1);
rfAmpAll=[rfAmpAll rfAmp];
rfPhaseAll=[rfPhaseAll rfPhase];
rfFreqAll=[rfFreqAll rfFreq];
rfCoilAll=[rfCoilAll rfCoil];
GzAmpAll=[GzAmpAll GzAmp];
GyAmpAll=[GyAmpAll GyAmp];
GxAmpAll=[GxAmpAll GxAmp];
ADCAll=[ADCAll ADC];
ExtAll=[ExtAll Ext];
rfTimeAll=[rfTimeAll rfTime];
GzTimeAll=[GzTimeAll GzTime];
GyTimeAll=[GyTimeAll GyTime];
GxTimeAll=[GxTimeAll GxTime];
ADCTimeAll=[ADCTimeAll ADCTime];
ExtTimeAll=[ExtTimeAll ExtTime];
SEtAll=[SEtAll SEt];
end
%====================================
if isempty(rfTimeAll)
error('rf sequence line can not be empty! Master Tx coil element must be used.');
end
if isempty(GzTimeAll)
error('GzSS sequence line can not be empty!');
end
if isempty(GyTimeAll)
error('GyPE sequence line can not be empty!');
end
if isempty(GxTimeAll)
error('GxR sequence line can not be empty!');
end
if isempty(ADCTimeAll)
error('ADC sequence line can not be empty!');
end
if isempty(ExtTimeAll)
error('Ext sequence line can not be empty!');
end
SEflag=repmat([0 0 0 0 0 0]',[1 2]);
rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTimeAll))]);
Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTimeAll))]);
Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTimeAll))]);
Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTimeAll))]);
ADCflag=repmat([0 0 0 0 1 0]',[1 max(size(ADCTimeAll))]);
Extflag=repmat([0 0 0 0 0 1]',[1 max(size(ExtTimeAll))]);
ts=[[min(SEtAll) max(SEtAll)] rfTimeAll GzTimeAll GyTimeAll GxTimeAll ADCTimeAll ExtTimeAll]-min(SEtAll);
flags=[SEflag rfflag Gzflag Gyflag Gxflag ADCflag Extflag];
[ts,ind]=sort(ts);
uts=unique(ts);
flags=flags(:,ind);
[rfTime,ind]=sort(rfTimeAll-min(SEtAll));
rfAmp=rfAmpAll(:,ind);
rfPhase=rfPhaseAll(:,ind);
rfFreq=rfFreqAll(:,ind);
rfCoil=rfCoilAll(:,ind);
[GzTime,ind]=sort(GzTimeAll-min(SEtAll));
GzAmp=GzAmpAll(:,ind);
[GyTime,ind]=sort(GyTimeAll-min(SEtAll));
GyAmp=GyAmpAll(:,ind);
[GxTime,ind]=sort(GxTimeAll-min(SEtAll));
GxAmp=GxAmpAll(:,ind);
[ADCTime,ind]=sort(ADCTimeAll-min(SEtAll));
ADC=ADCAll(:,ind);
[ExtTime,ind]=sort(ExtTimeAll-min(SEtAll));
Ext=ExtAll(:,ind);
rfAmp(1) = 0;
rfPhase(1) = 0;
rfFreq(1) = 0;
GzAmp(1) = 0;
GyAmp(1) = 0;
GxAmp(1) = 0;
ADC(1) = 0;
Ext(1) = 0;
rfAmp(end) = 0;
rfPhase(end) = 0;
rfFreq(end) = 0;
GzAmp(end) = 0;
GyAmp(end) = 0;
GxAmp(end) = 0;
ADC(end) = 0;
Ext(end) = 0;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [ D ] = dctn( A, varargin )
%
% N-D Discrete Cosine Transform
%
% Computes the N-dimensional dct of an array, by applying the
% one-dimnesional dct along all dimensions.
%
% Input:
%
%   A        - full N-dimensional array
%   varargin - (optional flag in which dimensions to apply the dct)
%
% Output:
%
%   D        - discrete cosine transform of A
%
%==============================================================================

function [ D ] = dctn( A, varargin )

if nargin==0
    help(mfilename)
    return;
end

% size if the input array
m = size(A);
dim = length(m);

% default parameters
dimFlag = ones(1,dim);  % apply dct in all dimensions by default

% overwrites default parameter
for k=1:2:length(varargin),
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

m = size(A);
D = A;
md = m;
P = circshift(1:dim,[0 -1]);

for d=1:dim
    if dimFlag(d)
        D = reshape(D,md(1),prod(md(2:end)));
        D = dct(D);
        D = reshape(D,md);
    end
    md = circshift(md,[0 -1]);
    D = permute(D,P);
end


text:
function varargout = lvd_EditGeometricVectorConstraintGUI(varargin)
% LVD_EDITGEOMETRICVECTORCONSTRAINTGUI MATLAB code for lvd_EditGeometricVectorConstraintGUI.fig
%      LVD_EDITGEOMETRICVECTORCONSTRAINTGUI, by itself, creates a new LVD_EDITGEOMETRICVECTORCONSTRAINTGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITGEOMETRICVECTORCONSTRAINTGUI returns the handle to a new LVD_EDITGEOMETRICVECTORCONSTRAINTGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITGEOMETRICVECTORCONSTRAINTGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITGEOMETRICVECTORCONSTRAINTGUI.M with the given input arguments.
%
%      LVD_EDITGEOMETRICVECTORCONSTRAINTGUI('Property','Value',...) creates a new LVD_EDITGEOMETRICVECTORCONSTRAINTGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditGeometricVectorConstraintGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditGeometricVectorConstraintGUI_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 lvd_EditGeometricVectorConstraintGUI

% Last Modified by GUIDE v2.5 13-May-2021 10:35:04

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditGeometricVectorConstraintGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditGeometricVectorConstraintGUI_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 lvd_EditGeometricVectorConstraintGUI is made visible.
function lvd_EditGeometricVectorConstraintGUI_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 lvd_EditGeometricVectorConstraintGUI (see VARARGIN)

    % Choose default command line output for lvd_EditGeometricVectorConstraintGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    constraint = varargin{1};
    setappdata(hObject, 'constraint', constraint);

    lvdData = varargin{2};
    setappdata(hObject,'lvdData',lvdData);

	populateGUI(handles, constraint, lvdData);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, constraint, lvdData)
    handles.constraintTypeLabel.String = constraint.getConstraintType();
    handles.constraintTypeLabel.TooltipString = constraint.getConstraintType();
    [unit, ~, ~, usesLbUb, ~, ~] = constraint.getConstraintStaticDetails();

    if(usesLbUb)
        [lb, ub] = constraint.getBounds();
        handles.ubText.String = fullAccNum2Str(ub);
        handles.lbText.String = fullAccNum2Str(lb);
    else
        handles.ubText.String = fullAccNum2Str(0);
        handles.lbText.String = fullAccNum2Str(0);
        handles.lbText.Enable = 'off';
        handles.ubText.Enable = 'off';
    end

    handles.scaleFactorText.String = fullAccNum2Str(constraint.getScaleFactor());

    handles.ubUnitLabel.String = unit;
    handles.lbUnitLabel.String = unit;

    evtNum = lvdData.script.getNumOfEvent(constraint.event);
    handles.eventCombo.String = lvdData.script.getListboxStr();
    if(not(isempty(evtNum)))
        handles.eventCombo.Value = evtNum;
    else
        handles.eventCombo.Value = 1;
    end

    handles.vectorCombo.String = lvdData.geometry.vectors.getListboxStr();

    if(isempty(constraint.vector))
        vectorInd = 1;
    else
        vectorInd = lvdData.geometry.vectors.getIndsForVectors(constraint.vector);
    end
    handles.vectorCombo.Value = vectorInd;

    frame = constraint.frame;
    if(isempty(frame))
        frame = LvdData.getDefaultInitialBodyInfo(lvdData.celBodyData).getBodyCenteredInertialFrame();
        constraint.frame = frame;
    end
    setappdata(handles.lvd_EditGeometricVectorConstraintGUI,'frame',frame);

    handles.refFrameTypeCombo.String = ReferenceFrameEnum.getListBoxStr();
    [ind, ~] = ReferenceFrameEnum.getIndForName(frame.typeEnum.name);
    handles.refFrameTypeCombo.Value = ind;
    handles.setFrameOptionsButton.TooltipString = sprintf('Current Frame: %s', frame.getNameStr());
    handles.refFrameNameLabel.String = sprintf('%s', frame.getNameStr());
    handles.refFrameNameLabel.TooltipString = sprintf('%s', frame.getNameStr());

    handles.constraintEvalTypeCombo.String = ConstraintEvalTypeEnum.getListBoxStr();
    handles.constraintEvalTypeCombo.Value = ConstraintEvalTypeEnum.getIndForName(constraint.evalType.name);
    constraintEvalTypeCombo_Callback(handles.constraintEvalTypeCombo, [], handles);

    evtNum = lvdData.script.getNumOfEvent(constraint.stateCompEvent);
    handles.compEventCombo.String = lvdData.script.getListboxStr();
    if(not(isempty(evtNum)))
        handles.compEventCombo.Value = evtNum;
    else
        handles.compEventCombo.Value = 1;
    end

    handles.eventNodeCombo.String = ConstraintStateComparisonNodeEnum.getListBoxStr();
    handles.eventNodeCombo.Value = ConstraintStateComparisonNodeEnum.getIndForName(constraint.eventNode.name);

    handles.stateComparisonNodeCombo.String = ConstraintStateComparisonNodeEnum.getListBoxStr();
    handles.stateComparisonNodeCombo.Value = ConstraintStateComparisonNodeEnum.getIndForName(constraint.stateCompNode.name);

    handles.comparisonTypeCombo.String = ConstraintStateComparisonTypeEnum.getListBoxStr();
    handles.comparisonTypeCombo.Value = ConstraintStateComparisonTypeEnum.getIndForName(constraint.stateCompType.name);

    handles.constActiveCheckbox.Value = constraint.active;

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditGeometricVectorConstraintGUI_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
    if(isempty(handles))
        varargout{1} = false;
    else
        constraint = getappdata(hObject, 'constraint');
        lvdData = getappdata(hObject,'lvdData');

        constraint.lb = str2double(handles.lbText.String);
        constraint.ub = str2double(handles.ubText.String);

        constraint.setScaleFactor(str2double(handles.scaleFactorText.String));

        constraint.vector = lvdData.geometry.vectors.getVectorAtInd(handles.vectorCombo.Value);

        frame = getappdata(hObject,'frame');
        constraint.frame = frame;

        [~, enums] = ConstraintEvalTypeEnum.getListBoxStr();
        constraint.evalType = enums(handles.constraintEvalTypeCombo.Value);

        [~, enums] = ConstraintStateComparisonTypeEnum.getListBoxStr();
        constraint.stateCompType = enums(handles.comparisonTypeCombo.Value);

        [~, enums] = ConstraintStateComparisonNodeEnum.getListBoxStr();
        constraint.eventNode = enums(handles.eventNodeCombo.Value);

        [~, enums] = ConstraintStateComparisonNodeEnum.getListBoxStr();
        constraint.stateCompNode = enums(handles.stateComparisonNodeCombo.Value);

        numEvents = lvdData.script.getTotalNumOfEvents();
        if(numEvents > 0)
            eventNum = handles.eventCombo.Value;
            constraint.event = lvdData.script.getEventForInd(eventNum);

            if(constraint.evalType == ConstraintEvalTypeEnum.StateComparison)
                eventNum = handles.compEventCombo.Value;
                constraint.stateCompEvent = lvdData.script.getEventForInd(eventNum);
            else
                constraint.stateCompEvent = LaunchVehicleEvent.empty(1,0);
            end
        else
            constraint.event = LaunchVehicleEvent.empty(1,0);
            constraint.stateCompEvent = LaunchVehicleEvent.empty(1,0);
        end

        constraint.active = logical(handles.constActiveCheckbox.Value);

        varargout{1} = true;
        close(handles.lvd_EditGeometricVectorConstraintGUI);
    end

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

    if(isempty(errMsg))
        uiresume(handles.lvd_EditGeometricVectorConstraintGUI);
    else
        msgbox(errMsg,'Invalid Constraint Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    constraint = getappdata(handles.lvd_EditGeometricVectorConstraintGUI, 'constraint');
    [~, lbLim, ubLim, usesLbUb, ~, ~] = constraint.getConstraintStaticDetails();

    if(usesLbUb)
        lwrBnd = str2double(get(handles.lbText,'String'));
        enteredStr = get(handles.lbText,'String');
        numberName = 'Lower Bound';
        lb = lbLim;
        ub = Inf;
        isInt = false;
        errMsg = validateNumber(lwrBnd, numberName, lb, ub, isInt, errMsg, enteredStr);

        uprBnd = str2double(get(handles.ubText,'String'));
        enteredStr = get(handles.ubText,'String');
        numberName = 'Upper Bound';
        lb = -Inf;
        ub = ubLim;
        isInt = false;
        errMsg = validateNumber(uprBnd, numberName, lb, ub, isInt, errMsg, enteredStr);

        if(isempty(errMsg))
            uprBnd = str2double(get(handles.ubText,'String'));
            enteredStr = get(handles.ubText,'String');
            numberName = 'Upper Bound';
            lb = lwrBnd;
            ub = ubLim;
            isInt = false;
            errMsg = validateNumber(uprBnd, numberName, lb, ub, isInt, errMsg, enteredStr);
        end
    end

    sf = str2double(get(handles.scaleFactorText,'String'));
    enteredStr = get(handles.scaleFactorText,'String');
    numberName = 'Scale Factor';
    lb = 1E-12;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(sf, numberName, lb, ub, isInt, errMsg, enteredStr);

    [~, enums] = ConstraintEvalTypeEnum.getListBoxStr();
    evalTypeEnum = enums(handles.constraintEvalTypeCombo.Value);
    if(evalTypeEnum == ConstraintEvalTypeEnum.StateComparison && ...
       handles.eventCombo.Value == handles.compEventCombo.Value && ...
       handles.eventNodeCombo.Value == handles.stateComparisonNodeCombo.Value)
        errMsg{end+1} = 'When using the state comparison evaluation type, the Applicable Event and Node and Comparison Event and Node must be different.';
    end

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

function lbText_Callback(hObject, eventdata, handles)
% hObject    handle to lbText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of lbText as text
%        str2double(get(hObject,'String')) returns contents of lbText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function lbText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lbText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes during object creation, after setting all properties.
function constraintTypeLabel_CreateFcn(hObject, eventdata, handles)
% hObject    handle to constraintTypeLabel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

function ubText_Callback(hObject, eventdata, handles)
% hObject    handle to ubText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ubText as text
%        str2double(get(hObject,'String')) returns contents of ubText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function ubText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ubText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in vectorCombo.
function vectorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to vectorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns vectorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from vectorCombo

% --- Executes during object creation, after setting all properties.
function vectorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to vectorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in eventCombo.
function eventCombo_Callback(hObject, eventdata, handles)
% hObject    handle to eventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns eventCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from eventCombo

% --- Executes during object creation, after setting all properties.
function eventCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function scaleFactorText_Callback(hObject, eventdata, handles)
% hObject    handle to scaleFactorText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of scaleFactorText as text
%        str2double(get(hObject,'String')) returns contents of scaleFactorText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function scaleFactorText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to scaleFactorText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditGeometricVectorConstraintGUI or any of its controls.
function lvd_EditGeometricVectorConstraintGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditGeometricVectorConstraintGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditGeometricVectorConstraintGUI);
    end

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

% Hint: get(hObject,'Value') returns toggle state of constActiveCheckbox

% --- Executes on selection change in constraintEvalTypeCombo.
function constraintEvalTypeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to constraintEvalTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns constraintEvalTypeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from constraintEvalTypeCombo
    contents = cellstr(get(hObject,'String'));
    name = contents{get(hObject,'Value')};

    enum = ConstraintEvalTypeEnum.getEnumForListboxStr(name);
    switch enum
        case ConstraintEvalTypeEnum.FixedBounds
            handles.ubText.Enable = 'on';
            handles.lbText.Enable = 'on';
            handles.compEventCombo.Enable = 'off';
            handles.comparisonTypeCombo.Enable = 'off';
            handles.stateComparisonNodeCombo.Enable = 'off';

        case ConstraintEvalTypeEnum.StateComparison
            handles.ubText.Enable = 'off';
            handles.lbText.Enable = 'off';
            handles.compEventCombo.Enable = 'on';
            handles.comparisonTypeCombo.Enable = 'on';
            handles.stateComparisonNodeCombo.Enable = 'on';

        otherwise
            error('Unknown constraint evaluation type.');
    end

% --- Executes during object creation, after setting all properties.
function constraintEvalTypeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to constraintEvalTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in compEventCombo.
function compEventCombo_Callback(hObject, eventdata, handles)
% hObject    handle to compEventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns compEventCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from compEventCombo

% --- Executes during object creation, after setting all properties.
function compEventCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to compEventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in comparisonTypeCombo.
function comparisonTypeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to comparisonTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns comparisonTypeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from comparisonTypeCombo

% --- Executes during object creation, after setting all properties.
function comparisonTypeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to comparisonTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in refFrameTypeCombo.
function refFrameTypeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to refFrameTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns refFrameTypeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from refFrameTypeCombo
    updateFrameChange(handles);

% --- Executes during object creation, after setting all properties.
function refFrameTypeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to refFrameTypeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in setFrameOptionsButton.
function setFrameOptionsButton_Callback(hObject, eventdata, handles)
% hObject    handle to setFrameOptionsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    curFrame = getappdata(handles.lvd_EditGeometricVectorConstraintGUI,'frame');
    curFrame = curFrame.editFrameDialogUI(EditReferenceFrameContextEnum.ForState);

    setappdata(handles.lvd_EditGeometricVectorConstraintGUI,'frame',curFrame);

    handles.setFrameOptionsButton.Tooltip = sprintf('Current Frame: %s', curFrame.getNameStr());
    handles.refFrameNameLabel.String = sprintf('%s', curFrame.getNameStr());
    handles.refFrameNameLabel.TooltipString = sprintf('%s', curFrame.getNameStr());

function updateFrameChange(handles)
    lvdData = getappdata(handles.lvd_EditGeometricVectorConstraintGUI,'lvdData');
    celBodyData = lvdData.celBodyData;

    refEnumListBoxStr = ReferenceFrameEnum.getListBoxStr();
    selFrameTypeInd = handles.refFrameTypeCombo.Value;
    refFrameEnum = ReferenceFrameEnum.getEnumForListboxStr(refEnumListBoxStr{selFrameTypeInd});

%             if(not(isempty(newFrame)) && newFrame.typeEnum ~= refFrameEnum)
%                 refFrameEnum = newFrame.typeEnum;
%                 app.refFrameTypeCombo.Value = ReferenceFrameEnum.getIndForName(refFrameEnum.name);
%             end

    switch refFrameEnum
        case ReferenceFrameEnum.BodyCenteredInertial
            bodyInfo = getSelectedBodyInfo(handles);

            newFrame = bodyInfo.getBodyCenteredInertialFrame();

        case ReferenceFrameEnum.BodyFixedRotating
            bodyInfo = getSelectedBodyInfo(handles);

            newFrame = bodyInfo.getBodyFixedFrame();

        case ReferenceFrameEnum.TwoBodyRotating
            bodyInfo = getSelectedBodyInfo(handles);
            if(not(isempty(bodyInfo.childrenBodyInfo)))
                primaryBody = bodyInfo;
                secondaryBody = bodyInfo.childrenBodyInfo(1);
            else
                primaryBody = bodyInfo.getParBodyInfo();
                secondaryBody = bodyInfo;
            end

            originPt = TwoBodyRotatingFrameOriginEnum.Primary;

            newFrame = TwoBodyRotatingFrame(primaryBody, secondaryBody, originPt, celBodyData);

        case ReferenceFrameEnum.UserDefined
            numFrames = lvdData.geometry.refFrames.getNumRefFrames();
            if(numFrames >= 1)
                geometricFrame = AbstractGeometricRefFrame.empty(1,0);
                for(i=1:numFrames)
                    frame = lvdData.geometry.refFrames.getRefFrameAtInd(i);
                    if(frame.isVehDependent() == false)
                        geometricFrame = frame;
                        break;
                    end
                end

                if(not(isempty(geometricFrame)))
                    newFrame = UserDefinedGeometricFrame(geometricFrame, lvdData);
                else
                    bodyInfo = getSelectedBodyInfo(handles);
                    newFrame = bodyInfo.getBodyCenteredInertialFrame();

                    warndlg('There are no geometric frames available which are not dependent on vehicle properties.  A body-centered inertial frame will be selected instead.');
                end
            else
                bodyInfo = getSelectedBodyInfo(handles);
                newFrame = bodyInfo.getBodyCenteredInertialFrame();

                warndlg('There are no geometric frames available which are not dependent on vehicle properties.  A body-centered inertial frame will be selected instead.');
            end

        otherwise
            error('Unknown reference frame type: %s', string(refFrameEnum));
    end

    if(not(isempty(newFrame)) && newFrame.typeEnum ~= refFrameEnum)
        refFrameEnum = newFrame.typeEnum;
        handles.refFrameTypeCombo.Value = ReferenceFrameEnum.getIndForName(refFrameEnum.name);
    end

    setappdata(handles.lvd_EditGeometricVectorConstraintGUI,'frame', newFrame);
    handles.setFrameOptionsButton.TooltipString = sprintf('Current Frame: %s', newFrame.getNameStr());
    handles.refFrameNameLabel.String = sprintf('%s', newFrame.getNameStr());
    handles.refFrameNameLabel.TooltipString = sprintf('%s', newFrame.getNameStr());

function bodyInfo = getSelectedBodyInfo(handles)
    curFrame = getappdata(handles.lvd_EditGeometricVectorConstraintGUI,'frame');
    bodyInfo = curFrame.getOriginBody();

% --- Executes on selection change in eventNodeCombo.
function eventNodeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to eventNodeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns eventNodeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from eventNodeCombo

% --- Executes during object creation, after setting all properties.
function eventNodeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eventNodeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in stateComparisonNodeCombo.
function stateComparisonNodeCombo_Callback(hObject, eventdata, handles)
% hObject    handle to stateComparisonNodeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns stateComparisonNodeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from stateComparisonNodeCombo

% --- Executes during object creation, after setting all properties.
function stateComparisonNodeCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to stateComparisonNodeCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function u=rotrack(u0,x,y,nu,T)
%
%  Solves the linear equation u_t + y u_x - x u_y = 0 by "front tracking" +
%  dimensional splitting, using a time step dictated by dt/dx=nu for each
%  substep.
%  Neumann boundary condition are implicit.
%
dx=x(2)-x(1);
dt=nu*dx;
nstep=ceil(T/dt);
dt=T/nstep;
nu=dt/dx;
u=u0;
for i=1:nstep,
	u=transp(u,nu*y,1);
	u=transp(u,-nu*x,2);
end;

repo: Aero-matlab


text:
function deformed_image = PTKInterpolateToMatch(original_image, deformation_field, interpolation_type, reporting)
    % PTKInterpolateToMatch. Interpolates an image to match the origin and voxel
    %     size specified in a deformation field
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    % This is similar to PTKDeformImage but there is no actual image deformation;
    % we are just reinterpolating to match the template provided by the
    % deformation field
    deformed_image = deformation_field.BlankCopy;
    deformed_image.ImageType = PTKImageType.Grayscale;
    [i_o, j_o, k_o] = original_image.GetGlobalCoordinatesMm;
    [i_o, j_o, k_o] = ndgrid(i_o, j_o, k_o);

    [i_r, j_r, k_r] = deformation_field.GetGlobalCoordinatesMm;
    [i_r, j_r, k_r] = ndgrid(i_r, j_r, k_r);

   deformed_image_raw = interpn(i_o, j_o, k_o, single(original_image.RawImage), ...
       i_r , ...
       j_r , ...
       k_r , ...
       interpolation_type, 0);
   deformed_image.ChangeRawImage(deformed_image_raw);
end

author: tomdoel

repo: pulmonarytoolkit


text:
% Chapter 7 - Differential Equations.
% Program_7c - Solving Simple Differential Equations.
% Copyright Birkhauser 2013. Stephen Lynch.

% Chemical kinetics (Example 8).
% Note that INLINE will be removed in future releases.
% deqn=inline('.00713*(4-c(1))^2*(1-c(1)/2)','t','c');
k=0.00713;
deqn = @(t,c) [k*(4-c)^2*(1-c/2)];
[t,ca]=ode45(deqn,[0 400],0);
plot(t,ca(:,1))
axis([0 400 0 3])
fsize=15;
set(gca,'XTick',0:100:400,'FontSize',fsize)
set(gca,'YTick',0:1:3,'FontSize',fsize)
xlabel('t','FontSize',fsize)
ylabel('c(t)','FontSize',fsize)
hold off

% End of Program_7c.

author: Sable

repo: mcbench-benchmarks


text:
function geometry_test0415 ( )

%*****************************************************************************80
%
%% TEST0415 tests LINES_IMP_INT_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0415\n' );
  fprintf ( 1, '  LINES_IMP_INT_2D finds the intersection of\n' );
  fprintf ( 1, '  two lines written in implicit form.\n' );
  fprintf ( 1, '\n' );
%
%  x + 2y - 4 = 0
%
  a1 =  1.0;
  b1 =  2.0;
  c1 = -4.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Line 1 coefficients:  %8f  %8f  %8f\n', a1, b1, c1 );
%
%  x - y - 1 = 0
%
  a2 =  1.0;
  b2 = -1.0;
  c2 = -1.0;
  fprintf ( 1, '  Line 2 coefficients:  %8f  %8f  %8f\n', a2, b2, c2 );

  [ ival, p ] = lines_imp_int_2d ( a1, b1, c1, a2, b2, c2 );

  if ( ival == 1 )
    fprintf ( 1, '  Intersection at %8f  %8f\n', p(1:dim_num) )
  elseif ( ival == 0 )
    fprintf ( 1, '  Lines are parallel, no intersection.\n' );
  elseif ( ival == 2 )
    fprintf ( 1, '  Lines are coincident.\n' );
  else
    fprintf ( 1, '  Unknown return value of ival = %d\n', ival );
  end
%
%  2x + 4y - 1 = 0
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Line 1 coefficients:  %8f  %8f  %8f\n', a1, b1, c1 );

  a2 =  2.0;
  b2 = +4.0;
  c2 = -1.0;
  fprintf ( 1, '  Line 2 coefficients:  %8f  %8f  %8f\n', a2, b2, c2 );

  [ ival, p ] = lines_imp_int_2d (a1, b1, c1, a2, b2, c2 );

  if ( ival == 1 )
    fprintf ( 1, '  Intersection at %8f  %8f\n', p(1:dim_num) );
  elseif ( ival == 0 )
    fprintf ( 1, '  Lines are parallel, no intersection.\n' );
  elseif ( ival == 2 )
    fprintf ( 1, '  Lines are coincident.\n' );
  else
    fprintf ( 1, '  Unknown return value of ival = %d\n', ival );
  end
%
%  -3x - 6y +12 = 0
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Line 1 coefficients:  %8f  %8f  %8f\n', a1, b1, c1 );

  a2 =  -3.0;
  b2 =  -6.0;
  c2 = +12.0;
  fprintf ( 1, '  Line 2 coefficients:  %8f  %8f  %8f\n', a2, b2, c2 );

  [ ival, p ] = lines_imp_int_2d  ( a1, b1, c1, a2, b2, c2 );

  if ( ival == 1 )
    fprintf ( 1, '  Intersection at %8f  %8f\n', p(1:dim_num) );
  elseif ( ival == 0 )
    fprintf ( 1, '  Lines are parallel, no intersection.\n' );
  elseif ( ival == 2 )
    fprintf ( 1, '  Lines are coincident.\n' );
  else
    fprintf ( 1, '  Unknown return value of ival = %d\n', ival );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%WAVE_SIGNIF  Significance testing for the 1D Wavelet transform WAVELET
%
%   [SIGNIF,FFT_THEOR] = ...
%      wave_signif(Y,DT,SCALE,SIGTEST,LAG1,SIGLVL,DOF,MOTHER,PARAM)
%
% INPUTS:
%
%    Y = the time series, or, the VARIANCE of the time series.
%        (If this is a single number, it is assumed to be the variance...)
%    DT = amount of time between each Y value, i.e. the sampling time.
%    SCALE = the vector of scale indices, from previous call to WAVELET.
%
%
% OUTPUTS:
%
%    SIGNIF = significance levels as a function of SCALE
%    FFT_THEOR = output theoretical red-noise spectrum as fn of PERIOD
%
%
% OPTIONAL INPUTS:
% *** Note *** setting any of the following to -1 will cause the default
%               value to be used.
%
%    SIGTEST = 0, 1, or 2.    If omitted, then assume 0.
%
%         If 0 (the default), then just do a regular chi-square test,
%             i.e. Eqn (18) from Torrence & Compo.
%         If 1, then do a "time-average" test, i.e. Eqn (23).
%             In this case, DOF should be set to NA, the number
%             of local wavelet spectra that were averaged together.
%             For the Global Wavelet Spectrum, this would be NA=N,
%             where N is the number of points in your time series.
%         If 2, then do a "scale-average" test, i.e. Eqns (25)-(28).
%             In this case, DOF should be set to a
%             two-element vector [S1,S2], which gives the scale
%             range that was averaged together.
%             e.g. if one scale-averaged scales between 2 and 8,
%             then DOF=[2,8].
%
%    LAG1 = LAG 1 Autocorrelation, used for SIGNIF levels. Default is 0.0
%
%    SIGLVL = significance level to use. Default is 0.95
%
%    DOF = degrees-of-freedom for signif test.
%         IF SIGTEST=0, then (automatically) DOF = 2 (or 1 for MOTHER='DOG')
%         IF SIGTEST=1, then DOF = NA, the number of times averaged together.
%         IF SIGTEST=2, then DOF = [S1,S2], the range of scales averaged.
%
%       Note: IF SIGTEST=1, then DOF can be a vector (same length as SCALEs),
%            in which case NA is assumed to vary with SCALE.
%            This allows one to average different numbers of times
%            together at different scales, or to take into account
%            things like the Cone of Influence.
%            See discussion following Eqn (23) in Torrence & Compo.
%
%
%----------------------------------------------------------------------------
%   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
%   University of Colorado, Program in Atmospheric and Oceanic Sciences.
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.
%----------------------------------------------------------------------------
function [signif,fft_theor] = ...
	wave_signif(Y,dt,scale1,sigtest,lag1,siglvl,dof,mother,param);

if (nargin < 9), param = -1;, end
if (nargin < 8), mother = -1;, end
if (nargin < 7), dof = -1;, end
if (nargin < 6), siglvl = -1;, end
if (nargin < 5), lag1 = -1;, end
if (nargin < 4), sigtest = -1;, end
if (nargin < 3)
	error('Must input a vector Y, sampling time DT, and SCALE vector')
end

n1 = length(Y);
J1 = length(scale1) - 1;
scale(1:J1+1) = scale1;
s0 = min(scale);
dj = log(scale(2)/scale(1))/log(2.);

if (n1 == 1)
	variance = Y;
else
	variance = std(Y)^2;
end

if (sigtest == -1), sigtest = 0;, end
if (lag1 == -1), lag1 = 0.0;, end
if (siglvl == -1), siglvl = 0.95;, end
if (mother == -1), mother = 'MORLET';, end

mother = upper(mother);

% get the appropriate parameters [see Table(2)]
if (strcmp(mother,'MORLET'))  %----------------------------------  Morlet
	if (param == -1), param = 6.;, end
	k0 = param;
	fourier_factor = (4*pi)/(k0 + sqrt(2 + k0^2)); % Scale-->Fourier [Sec.3h]
	empir = [2.,-1,-1,-1];
	if (k0 == 6), empir(2:4)=[0.776,2.32,0.60];, end
elseif (strcmp(mother,'PAUL'))  %--------------------------------  Paul
	if (param == -1), param = 4.;, end
	m = param;
	fourier_factor = 4*pi/(2*m+1);
	empir = [2.,-1,-1,-1];
	if (m == 4), empir(2:4)=[1.132,1.17,1.5];, end
elseif (strcmp(mother,'DOG'))  %---------------------------------  DOG
	if (param == -1), param = 2.;, end
	m = param;
	fourier_factor = 2*pi*sqrt(2./(2*m+1));
	empir = [1.,-1,-1,-1];
	if (m == 2), empir(2:4) = [3.541,1.43,1.4];, end
	if (m == 6), empir(2:4) = [1.966,1.37,0.97];, end
else
	error('Mother must be one of MORLET,PAUL,DOG')
end

period = scale.*fourier_factor;
dofmin = empir(1);     % Degrees of freedom with no smoothing
Cdelta = empir(2);     % reconstruction factor
gamma_fac = empir(3);  % time-decorrelation factor
dj0 = empir(4);        % scale-decorrelation factor

freq = dt ./ period;   % normalized frequency
fft_theor = (1-lag1^2) ./ (1-2*lag1*cos(freq*2*pi)+lag1^2);  % [Eqn(16)]
fft_theor = variance*fft_theor;  % include time-series variance
signif = fft_theor;
if (dof == -1), dof = dofmin;, end

if (sigtest == 0)    % no smoothing, DOF=dofmin [Sec.4]
	dof = dofmin;
	chisquare = chisquare_inv(siglvl,dof)/dof;
	signif = fft_theor*chisquare ;  % [Eqn(18)]
elseif (sigtest == 1)  % time-averaged significance
	if (length(dof) == 1), dof=zeros(1,J1+1)+dof;, end
	truncate = find(dof < 1);
	dof(truncate) = ones(size(truncate));
	dof = dofmin*sqrt(1 + (dof*dt/gamma_fac ./ scale).^2 );   % [Eqn(23)]
	truncate = find(dof < dofmin);
	dof(truncate) = dofmin*ones(size(truncate));   % minimum DOF is dofmin
	for a1 = 1:J1+1
		chisquare = chisquare_inv(siglvl,dof(a1))/dof(a1);
		signif(a1) = fft_theor(a1)*chisquare;
	end
elseif (sigtest == 2)  % time-averaged significance
	if (length(dof) ~= 2)
		error('DOF must be set to [S1,S2], the range of scale-averages')
	end
	if (Cdelta == -1)
		error(['Cdelta & dj0 not defined for ',mother, ...
			' with param = ',num2str(param)])
	end
	s1 = dof(1);
	s2 = dof(2);
	avg = find((scale >= s1) & (scale <= s2));  % scales between S1 & S2
	navg = length(avg);
	if (navg == 0)
		error(['No valid scales between ',num2str(s1),' and ',num2str(s2)])
	end
	Savg = 1./sum(1 ./ scale(avg));       % [Eqn(25)]
	Smid = exp((log(s1)+log(s2))/2.);     % power-of-two midpoint
	dof = (dofmin*navg*Savg/Smid)*sqrt(1 + (navg*dj/dj0)^2);  % [Eqn(28)]
	fft_theor = Savg*sum(fft_theor(avg) ./ scale(avg));  % [Eqn(27)]
	chisquare = chisquare_inv(siglvl,dof)/dof;
	signif = (dj*dt/Cdelta/Savg)*fft_theor*chisquare;    % [Eqn(26)]
else
	error('sigtest must be either 0, 1, or 2')
end

return

author: grinsted

repo: wavelet-coherence


text:
function [EE,I] = conservative_edge_matching(E_orig,varargin)
  % CONSERVATIVE_EDGE_MATCHING  Find set of given edges so that no two edges
  % share a vertex.
  %
  % [EE,I] = conservative_edge_matching(E)
  %
  % Inputs:
  %   E  #E by 2 list of edges
  %   Optional:
  %      'Method'  followed by one of the following:
  %        'one-per-component'  simply take one edge per conneceted component
  %          of edges
  %        'random'  use a randomized algorithm to try to grab many per
  %          component
  % Outputs:
  %   EE  #EE by 2 list of edges
  %   I  #EE list of indices into E, so that EE = E(I,:)
  %

  if size(E_orig,1) == 0
    EE = E_orig;
    I = [];
    C = [];
    return;
  end

  % default values
  method = 'one-per-component';
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Method'}, ...
    {'method'});
  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

  % Do this in reduced graph of just vertices incident on E. Otherwise this is
  % O(max(E(:))) rather than O(size(E))
  [~,~,E] = unique(E_orig(:));
  E = reshape(E,size(E_orig));

  switch method
  case 'one-per-component'
    % Find one edge per connected component of edges (ideally we would find a
    % "perfect matching")
    ne = size(E,1);
    C = connected_components(E);
    C = C(E(:,1));
    % Use reverse order so we can take max
    E2C = sparse(1:ne,C,ne-(1:ne)+1,ne,max(C));
    I = max(E2C,[],1);
    I = I(I>0);
    % reverse
    I = ne - I + 1;
  case 'recursive'
    keep = 1:size(E,1);
    Ekeep = E(keep,:);
    J = 1:size(E,1);
    I = [];
    while true
      % This should use random for large meshes
      [~,Ikeep] = conservative_edge_matching(Ekeep,'Method','random');
      I = [I J(Ikeep)];
      % Remove edges connected to E
      keep = find(~any(ismember(Ekeep,Ekeep(Ikeep,:)),2));
      if isempty(keep)
        break;
      end
      % Indices in original edges
      J = J(keep);
      Ekeep = Ekeep(keep,:);
    end

  case 'random'
    ns = size(E,1);
    while true
      I = randperm(size(E,1));
      I = I(1:ns);
      % Reduced graph of randomly selected edges
      [~,~,EI] = unique(E(I,:));
      % Check if every vertex is valence 1
      V = full(sparse(EI,1,1));
      if max(V) == 1
        break;
      end
      ns = floor(ns/2);
      if ns <= 1
        I = 1;
        break;
      end
    end
  end

  % Output on original graph
  EE = E_orig(I,:);
end

author: alecjacobson

repo: gptoolbox


text:
function hammersley_test01 ( )

%*****************************************************************************80
%
%% TEST01 tests I4_TO_HAMMERSLEY_SEQUENCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  I4_TO_HAMMERSLEY_SEQUENCE computes N elements of\n' );
  fprintf ( 1, '  a Hammersley sequence on a single call.\n' );
  fprintf ( 1, '  All arguments are specified explicitly.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this example, we compute the first 10 elements\n' );
  fprintf ( 1, '  of a "classical" Hammersley sequence, and then\n' );
  fprintf ( 1, '  the "last" 10 elements.\n' );

  nmax = 1000;

  dim_num = 4;
  n = 10;
  step = 1;
  seed(1:dim_num) = 0;
  leap(1:dim_num) = 1;
  base(1) = -nmax;
  for i = 2 : dim_num
    base(i) = prime ( i - 1 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  DIM_NUM = %12d\n', dim_num );
  fprintf ( 1, '  N =    %12d\n', n );
  fprintf ( 1, '  STEP = %12d\n', step );
  i4vec_transpose_print ( dim_num, seed, '  SEED = ' );
  i4vec_transpose_print ( dim_num, leap, '  LEAP = ' );
  i4vec_transpose_print ( dim_num, base, '  BASE = ' );

  r = i4_to_hammersley_sequence ( dim_num, n, step, seed, leap, base );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    STEP   Hammersley\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  ', step+j-1 );
    for i = 1 : dim_num
      fprintf ( 1, '%12f  ', r(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  We can jump ahead in the sequence by changing STEP.\n' );

  step = nmax - n + 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  STEP = %12d\n', step );

  r = i4_to_hammersley_sequence ( dim_num, n, step, seed, leap, base );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    STEP   Hammersley\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  ', step+j-1 );
    for i = 1 : dim_num
      fprintf ( 1, '%12f  ', r(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function [W] = readDMAT(filename)
  % READDMAT  read a matrix from a dmat file.  first line is <# columns> <#
  % rows>, then values with columns running faster
  %
  % [W] = readDMAT(filename)
  %
  % Input:
  %   filename  name of .dmat file
  % Output:
  %   W  matrix read from file
  %
  % See also: writeDMAT
  %

  % open file
  fp = fopen(filename,'r');
  % read header
  size_W = fliplr(fscanf(fp,'%d %d',[1 2]));
  % read data
  W = fscanf(fp,'%g',size_W);

  if ~feof(fp)
    [size_B,c_B] = fscanf(fp,'%d %d',[1 2]);
    size_B = fliplr(size_B);
    if c_B==2 && all(size_B>0)
      assert(~any(size_W));
      assert(isempty(W));
      % Finish reading header: read '\n' char
      [lf,nlf] = fread(fp,1,'char*1');
      assert(nlf==1);
      assert(lf==int8(sprintf('\n')));
      % We're reading binary then
      size_W = size_B;
      [W,cW] = fread(fp,prod(size_W),'*double');
      assert(cW == prod(size_W));
      W = reshape(W,size_W);
    end
  end

  % close file
  fclose(fp);

  % size should match header
  if(~all(size_W == size(W)))
    error('Size in header (%d,%d) did not match size of data (%d,%d) in file',size_W,size(W));
  end

end


text:
function val = get(A, prop_name)
%
%  Get transformMatrix properties from specified object and return
%  the value.
%
%     val = get(A, prop_name);
%
%  Valid choices for prop_name are:
%    'transform', 'transpose'
%

%  J. Nagy  6/2/02

switch prop_name
case 'transform'
  val = A.transform;
case 'transpose'
  val = A.transpose;
otherwise
  error([prop_name, 'Is not valid transformMatrix property'])
end

repo: IRtools


text:
function b = r8vec_bracket6 ( nd, xd, ni, xi )

%*****************************************************************************80
%
%% R8VEC_BRACKET6 brackets data between successive entries of a sorted R8VEC.
%
%  Discussion:
%
%    We assume XD is sorted.
%
%    If XI(I) is contained in the interval [XD(1),XD(N)], then the value of
%    B(I) indicates that XI(I) is contained in [ XD(B(I)), XD(B(I)+1) ].
%
%    If XI(I) is not contained in the interval [XD(1),XD(N)], then B(I) = -1.
%
%    This code implements a version of binary search which is perhaps more
%    understandable than the usual ones.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ND, the number of data values.
%
%    Input, real XD(N,1), the sorted data.
%
%    Input, integer NI, the number of inquiry values.
%
%    Input, real XI(NI,1), the query values.
%
%    Output, integer B(NI,1), the bracket information.
%
  xd = xd(:);
  xi = xi(:);
  b = zeros ( ni, 1 );

  for i = 1 : ni

    if ( xi(i) < xd(1) || xd(nd) < xi(i) )

      b(i) = -1;

    else

      l = 1;
      r = nd;

      while ( l + 1 < r )
        m = floor ( ( l + r ) / 2 );
        if ( xi(i) < xd(m) )
          r = m;
        else
          l = m;
        end
      end

      b(i) = l;

    end

  end

  return
end


text:
function [refBody, scName] = orbitPanelGetOrbitFromKSPTOTConnectActiveVesselCallBack(hSMA, hECC, hINC, hRAAN, hARG, varargin)

    refBody = [];
    [orbit] = getSingularOrbitFromKSPTOTConnect([]);
    if(not(isempty(orbit)))
        if(not(isempty(hSMA)))
            set(hSMA, 'String', fullAccNum2Str(orbit{3}));
        end

        if(not(isempty(hECC)))
            set(hECC, 'String', fullAccNum2Str(orbit{4}));
        end

        if(not(isempty(hINC)))
            set(hINC, 'String', fullAccNum2Str(orbit{5}));
        end

        if(not(isempty(hRAAN)))
            set(hRAAN, 'String', fullAccNum2Str(orbit{6}));
        end

        if(not(isempty(hARG)))
            set(hARG, 'String', fullAccNum2Str(orbit{7}));
        end

        refBody = orbit{10};
        scName = orbit{1};

        if(nargin == 7)
            hMA = varargin{1};
            hEpoch = varargin{2};
            set(hMA, 'String', fullAccNum2Str(rad2deg(orbit{8})));
            set(hEpoch, 'String', fullAccNum2Str(orbit{9}));
        end
    end
end


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function [BB2 Conf Valid tld] = tldTrack_occlusion(tld,BB1,I,J)

GRID      = tld.tracker.grid;
BIG_FB    = tld.tracker.big_fb;
% OCCLUSION = tld.tracker.occlusion;

BB2    = [];
Conf   = [];
Valid  = 0;

if isempty(BB1) || ~bb_isdef(BB1),
    return;
end

% Forward
xFI    = bb_points(BB1,GRID,GRID,5);
xFJ    = lk(2,tld.img{I}.input,tld.img{J}.input,xFI,xFI,3);
idxFB  = xFJ(3,:) < 2;
idxNCC = xFJ(4,:) > 0.9;

% Occlusion detection
 xF1    = bb_points(tld.bb(:,1),GRID,GRID,5);
 xB1    = lk(2,tld.img{I}.input,tld.img{1}.input,xFI,xFI,1);
 idxOCC = xB1(3,:) < .1 & xB1(4,:) > 0.95 & vnormp(xF1(1:2,:)-xB1(1:2,:),2) > 10;

xFI    = xFI(:,~idxOCC & idxFB & idxNCC);
xFJ    = xFJ(:,~idxOCC & idxFB & idxNCC);

medFB  = median2(xFJ(3,:));
medNCC = median2(xFJ(4,:));
idxF   = xFJ(3,:) <= medFB & xFJ(4,:)>= medNCC;
BB2    = bb_predict(BB1,xFI(:,idxF),xFJ(1:2,idxF));

% Failure detection
if medFB > BIG_FB || ~bb_isdef(BB2) || bb_isout(BB2,tld.imgsize)
    disp('Tracker uncertain.');
    BB2 = [];
    return;
end

% Output
ex       = tldGetPattern(tld.img{J},BB2,tld.model.patchsize);
[~,Conf] = tldNN(ex,tld);
Valid    = tld.valid(I);

% if sum(idxOCC)/length(idxOCC) >= .20
%     tld.trackerfailure(J) = 1;
%     disp('Tracker occlusion.');
%     %disp('tracker failure');
%     Valid = 0;
% end

% Match error
% resF  = mean(vnormp(xFJ(1:2,idxF)-xBJ(:,idxF),2)); % residual of forward match
% resB  = mean(vnormp(xBI(1:2,idxB)-xFI(:,idxB),2)); % residual of backward match
% %BB_shift = vnormp(BB_shift,2);
% error = abs(resF-resB);

% if ~bb_isin(BB2,tld.imgsize) || error > OCCLUSION
%     Valid = 0;
% end

repo: OpenTLD


text:
function p = perms_m(n, m)
%perms_m(n, m): return the mth permutation of the digits 1:n
% where m is from 0 to factorial(n) - 1
%
% The ordering is fairly arbitrary, but there is a bijection between m to
% the n! permutations (asking for same m twice gets same permutation), and
% perms_m(n, 0) is 1:n in order.
%
% Example:
%   p = perms(1:4);
%   P = zeros(size(p));
%   for i = 1:size(P,1)
%       P(i, :) = perms_m(4, i-1)';
%   end
%   all(all(sortrows(P) == sortrows(p)))
%
% See also: perms, randperm, uperms, signs_m, nchoosek_m, sample_no_repl

% Copyright 2010 Ged Ridgway
% http://www.mathworks.com/matlabcentral/fileexchange/authors/27434

n = ceil(abs(n)); m = round(abs(m));

M = factorial(n) - 1; % max m
if m < 0 || m > M
    error('perms_m:m_range', 'Require 0 <= m <= %d', M);
end

p = zeros(n,1);
d = 1:n;
F = factorial(n); % (F in loop goes from factorial(n-1) to factorial(0))
for f = 1:n
    F = F / (n + 1 - f); % (here rather than end of loop to avoid div-by-0)
    x = floor(m / F);
    p(f) = d(f + x);
    d(f + x) = d(f);
    m = rem(m, F);
end


text:
function [tissueModel,coreRxnBool,coreMetBool,coreCtrsBool] = fastcore(model, coreRxnInd, epsilon, printLevel)
% Use the FASTCORE algorithm ('Vlassis et al, 2014') to extract a context
% specific model. FASTCORE algorithm defines one set of core
% reactions that is guaranteed to be active in the extracted model and find
% the minimum of reactions possible to support the core.
%
% USAGE:
%
%    tissueModel = fastcore(model, coreRxnInd)
%
% INPUTS:
%    model:             (the following fields are required - others can be supplied)
%                         * S  - `m x 1` Stoichiometric matrix
%                         * lb - `n x 1` Lower bounds
%                         * ub - `n x 1` Upper bounds
%                         * rxns   - `n x 1` cell array of reaction abbreviations
%
%   coreRxnInd:          indices of reactions in cobra model that are part of the
%                        core set of reactions (called 'C' in 'Vlassis et al,
%                        2014')
%
% OPTIONAL INPUTS:
%   epsilon:             smallest flux value that is considered nonzero
%                        (default getCobraSolverParams('LP', 'feasTol')*100)
%   printLevel:          0 = silent, 1 = summary, 2 = debug (default - 0)
%
% OUTPUT:
%
%   tissueModel:         extracted model
%
%   coreRxnBool:         n x 1 boolean vector indicating core reactions
%
% 'Vlassis, Pacheco, Sauter (2014). Fast reconstruction of compact
% context-specific metbolic network models. PLoS Comput. Biol. 10,
% e1003424.'
%
% .. Authors:
%       - Nikos Vlassis, Maria Pires Pacheco, Thomas Sauter, 2013 LCSB / LSRU, University of Luxembourg
%       - Ronan Fleming, commenting of code and inputs/outputs
%       - Anne Richelle, code adaptation to fit with createTissueSpecificModel

if nargin < 4 || ~exist('printLevel','var')
    printLevel = 0;
end
if nargin < 3 || isempty(epsilon)
    epsilon=getCobraSolverParams('LP', 'feasTol')*100;
end

if printLevel > 1
    tic
end

model_orig = model;

[nMets,nRxns] = size(model.S);

LPproblem = buildLPproblemFromModel(model);

%reactions irreversible in the reverse direction
Ir = model.lb < 0 & model.ub<=0;

%flip direction of reactions irreversible in the reverse direction
LPproblem.A(:,Ir) = -LPproblem.A(:,Ir);
tmp = LPproblem.ub(Ir);
LPproblem.ub(Ir) = -LPproblem.lb(Ir);
LPproblem.lb(Ir) = -tmp;

%Find irreversible reactions
irrevRxns = find(model.lb>=0);

A = [];
flipped = false;
singleton = false;

% Find irreversible core reactions
J = intersect(coreRxnInd, irrevRxns);

if printLevel > 0
    fprintf('|J|=%d  ', length(J));
end

%Find all the reactions that are not in the core
nbRxns = 1:nRxns;
% Non Core reactions (penalized)
P = setdiff(nbRxns, coreRxnInd);

% Find the minimum of set reactions from P that need to be included to
% support the irreversible core set of reactions
[Supp, basis] = findSparseMode(J, P, singleton, model, LPproblem, epsilon);

if ~isempty(setdiff(J, Supp))
    warning('fastcore.m Error: Global network is not flux consistent, ignoring the following irreversible core reactions:\n');
    model.rxns(setdiff(J, Supp))
    coreRxnInd = setdiff(coreRxnInd,setdiff(J, Supp));
    %error ('fastcore.m Error: Inconsistent irreversible core reactions.\n');
end

A = Supp;
if printLevel > 0
    fprintf('|A|=%d\n', length(A));
end

% J is the set of core reactions not already in the extracted model
J = setdiff(coreRxnInd, A);
if printLevel > 0
    fprintf('|J|=%d  ', length(J));
end

% Main loop that reduce at each iteration the number of reactions from P that need to be included to
% support the complete core set of reactions
while ~isempty(J)

    P = setdiff(P, A);

    %reuse the basis from the previous solve if it exists
    [Supp, basis] = findSparseMode(J, P, singleton, model, LPproblem, epsilon, basis);

    A = union(A, Supp);
    if printLevel > 0
        fprintf('|A|=%d\n', length(A));
    end

    if ~isempty( intersect(J, A))
        J = setdiff(J, A);
        if printLevel > 0
            fprintf('|J|=%d  ', length(J));
        end
        flipped = false;
    else
        if singleton
            JiRev = setdiff(J(1),irrevRxns);
        else
            JiRev = setdiff(J,irrevRxns);
        end
        if flipped || isempty(JiRev)
            if singleton
                warning('\n fastcore.m: Global network is not flux consistent, ignoring corresponding core reaction:\n');
                disp(model.rxns(J))
                J = [];
                %error('\n fastcore.m Error: Global network is not consistent.\n');
            else
                flipped = false;
                singleton = true;
            end
        else
            LPproblem.A(:,JiRev) = -LPproblem.A(:,JiRev);
            tmp = LPproblem.ub(JiRev);
            LPproblem.ub(JiRev) = -LPproblem.lb(JiRev);
            LPproblem.lb(JiRev) = -tmp;
            flipped = true;

            if printLevel > 0
                fprintf('(flip)  ');
            end
        end
    end
end
if printLevel > 0
    fprintf('|A|=%d\n', length(A)); % A : indices of reactions in the new model
end

if printLevel > 1
    toc
end

coreRxnBool=false(size(model.S,2),1);
coreRxnBool(A)=1;

rxnRemoveList = setdiff(model.rxns,model.rxns(A));

dummyMetBool = contains(model.mets,'dummy_Met_');
dummyRxnBool = contains(model.rxns,'dummy_Rxn_');
dummyRxnList = model.rxns(dummyRxnBool);
if any(dummyMetBool) || any(dummyRxnBool)
    model = destroyDummyModel(model,dummyMetBool,dummyRxnBool);
    rxnRemoveList = setdiff(rxnRemoveList,dummyRxnList);
end

%removes any infeasible coupling constraints also
[tissueModel, metRemoveList, ctrsRemoveList] = removeRxns(model, rxnRemoveList,'metRemoveMethod','exclusive','ctrsRemoveMethod','infeasible');

coreMetBool=~ismember(model_orig.mets,metRemoveList);
if isfield(model,'ctrs')
    coreCtrsBool = ~ismember(model_orig.ctrs,ctrsRemoveList);
else
    coreCtrsBool = ctrsRemoveList;
end

%coreGeneBool
tissueModel = removeUnusedGenes(tissueModel);

author: opencobra

repo: cobratoolbox


text:
function [masterOutput, masterTime] = TS_ComputeMasterLoop(x,x_z,masterCode,masterID,numMasterOps,howVocal,theTsID,iterNum)
% TS_ComputeMasterLoop     Used in a loop by TS_Compute to evaluate a given master function.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

if strcmp(howVocal,'full')
    beVocal = true;
else
    beVocal = false;
end

if beVocal
    % Display code name for error checking
    fprintf(1,'[TimeSeries_ID = %u, MasterOperation_ID = %u (%u/%u)] %s...',...
                theTsID, masterID, iterNum, numMasterOps, masterCode);
end

try
	masterTimer = tic;
    if beVocal
        % Any output text is printed to screen
    	masterOutput = BF_ParEval(x,x_z,masterCode,1);
    else
        % Output text stored in second output (could log this if you really want to)
        masterOutput = BF_ParEval(x,x_z,masterCode,0);
    end
	masterTime = toc(masterTimer);
    if beVocal
        fprintf(1,' evaluated (%s).\n',BF_TheTime(masterTime));
    end
	% For not-applicable/'real NaN', masterOutput is a NaN, otherwise a
	% structure with components to be called below by pointer operations.

catch emsg
    if beVocal
        fprintf(1,' error.\n'); % ,BF_TheTime(masterTime)
    end
	fprintf(1,'---Error evaluating %s:\n%s.\n',masterCode,emsg.message);
    masterOutput = {}; % Keep empty output
    masterTime = 0; % Set zero calculation time
	% Remains an empty cell entry.
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [l,n,k,k2] = get_dim(dat)
%   [numEx,vDim,oDim,numCls]=get_dim(data)
%    return (in given order) the number of examples, the dimension of vectors
%    the number of output dimensions and the number of classes.
%
%    Note:
%    Sometimes the number of classes differs from the output dimension
%    (e.g. in binary pattern recognition often oDim equals 1 and numCls
%    equals 2).

  l=length(dat.index);
  n=length(dat.findex);
  k=size(dat.Y,2);
  k2=k;

  if k2==1
    f1=find(dat.Y==1); f2=find(dat.Y==-1);
    if length(f1)>0 & length(f2)>0 & (length(f1)+length(f2))==l
      k2=2;  % -1s and +1s => two classes
    end
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef L2normalization < dagnn.ElementWise
    properties (Transient)
        numInputs
        SIZE_
    end

    methods
        function outputs = forward(obj, inputs, params)
            obj.numInputs = numel(inputs);
            assert(obj.numInputs == 1);

            X = inputs{1};
            obj.SIZE_ = size(X);
            if length(obj.SIZE_) < 3
                obj.SIZE_ = [obj.SIZE_, 1, 1];
            elseif length(obj.SIZE_) ==3
                obj.SIZE_ = [obj.SIZE_ 1];
            elseif length(obj.SIZE_) > 4
                error('length(SIZE_) > 4');
            end

            X = permute(X, [3,1,2,4]);
            X = reshape(X, [size(X,1), prod(obj.SIZE_)/obj.SIZE_(3)] );
            sumX = sqrt(sum(X.^2,1));
            %sumX = sumX + (sumX==0);
            sumX = sumX + 0.00001;
            X = bsxfun(@rdivide, X, sumX);

            X = reshape(X, [size(X,1), obj.SIZE_(1:2), obj.SIZE_(4)]);
            X = permute(X, [2,3,1,4]);
            outputs{1} = X;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            derInputs = cell(1, numel(inputs));

            dzdy = derOutputs{1};
            X = inputs{1};
%             obj.SIZE_ = size(X);
%             if length(obj.SIZE_) < 3
%                 obj.SIZE_ = [obj.SIZE_, 1, 1];
%             elseif length(obj.SIZE_) ==3
%                 obj.SIZE_ = [obj.SIZE_ 1];
%             elseif length(obj.SIZE_) > 4
%                 error('length(SIZE_) > 4');
%             end

            X = permute(X, [3,1,2,4]);
            X = reshape(X, [size(X,1), prod(obj.SIZE_)/obj.SIZE_(3)] );
            dzdy = permute(dzdy, [3,1,2,4]);
            dzdy = reshape(dzdy, [size(dzdy,1), prod(obj.SIZE_)/obj.SIZE_(3)] );

            lambda = 1./(sqrt(sum(X.^2, 1)) + 1e-10);
            dzdx = bsxfun(@times, lambda, dzdy) - bsxfun(@times, X, (lambda.^3) .* sum(X.*dzdy, 1));

            dzdx = reshape(dzdx, [size(dzdx,1), obj.SIZE_(1:2), obj.SIZE_(4)]);
            dzdx = permute(dzdx, [2,3,1,4]);

            derInputs{1} = dzdx;
            derParams = {} ;
        end

        function obj = Scale(varargin)
            obj.load(varargin) ;
        end
    end
end

author: aimerykong


text:
function p=imcorr(sys,par,dp)
%IMCORR corrects image coordinates, which are contaminated by radial
%and tangential distortion.
%
%Usage:
%   p=imcorr(name,par2,dp)
%
%where
%   name = string that is specific to the camera and the framegrabber.
%          This string must be defined in configc.m
%   par2 = camera intrinsic parameters for correcting the coordinates.
%          these parameters are computed by using invmodel.m.
%   dp   = distorted image coordinates in pixels (n x 2 matrix)
%   p    = corrected image coordinates

%   Version 3.0  10-17-00
%   Janne Heikkila, University of Oulu, Finland

NDX=sys(1); NDY=sys(2); Sx=sys(3); Sy=sys(4);
Asp=par(1); Foc=par(2);
Cpx=par(3); Cpy=par(4);
Rad1=par(5); Rad2=par(6);
Tan1=par(7); Tan2=par(8);

dx=(dp(:,1)-Cpx)*Sx/NDX/Asp;
dy=(dp(:,2)-Cpy)*Sy/NDY;

r2=dx.*dx+dy.*dy;
delta=Rad1*r2+Rad2*r2.*r2;

cx=dx.*(1+delta)+2*Tan1*dx.*dy+Tan2*(r2+2*dx.*dx);
cy=dy.*(1+delta)+Tan1*(r2+2*dy.*dy)+2*Tan2*dx.*dy;

p=NDX*Asp*cx/Sx+Cpx;
p(:,2)=NDY*cy/Sy+Cpy;

author: strawlab

repo: MultiCamSelfCal