text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the HEART-TUBE-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function HeartTube()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  128;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  128;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 5.0;        % Length of Eulerian Grid in x-Direction
Ly = 5.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
ds = Lx/(2*Nx);  % Lagrangian Spacing
d = 1.0;         % Diameter of the tube
L = 3.0;         % Length of the heart-tube
struct_name = 'HeartTube'; % Name for .vertex, .spring, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,L,d,Lx,Ly);

% Plot Geometry to test
plot(xLag(1:end/2),yLag(1:end/2),'r-'); hold on;
plot(xLag(end/2+1:end),yLag(end/2+1:end),'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Ly]);

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
k_Spring = 1e7;
print_Lagrangian_Springs(xLag,k_Spring,ds,struct_name);

% Prints .muscle file! [ a_f * Fmax *exp( -( (Q-1)/SK )^2 ) * (1/P0)*(b*P0-a*v)/(v+b); Q = LF/LFO ]
LFO = d; SK = 0.3; a = 0.25; b = 4.0; Fmax = 1e5;
print_Lagrangian_Muscles(xLag,LFO,SK,a,b,Fmax,struct_name)

% Prints .beam file!
k_Beam = 7.5e7; C = 0.0;
print_Lagrangian_Beams(xLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1e6;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag); % Total Number of Lagrangian Pts
    num = 1;          % Number of target points on each end

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', 4*num );

    %Loops over all Lagrangian Pts.
    %for s = 1:N

    %Left Bottom Target Points
    for s=1:num
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %Right Bottom Target Points
    for s=(N/2)-(num-1):N/2
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %Left Top Target Points
    for s=(N/2+1):(N/2+1)+(num-1)
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %Right Top Target Points
    for s=N-(num-1):N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N - 4 );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N/2-1
                % Bottom of tube
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            elseif ( ( s >= N/2+2 ) && ( s <= N-1 ) )
                % Top of Tube
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints MUSCLE points to a file called "struct_name".muscle
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Muscles(xLag,LFO,SK,a,b,Fmax,struct_name)

    N = length(xLag); %Number of Lagrangian Pts. Total

    muscle_fid = fopen([struct_name '.muscle'], 'w');

    fprintf(muscle_fid, '%d\n', N/2-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %MUSCLES BETWEEN VERTICES
    for s = 2:N/2-1
            fprintf(muscle_fid, '%d %d %1.16e %1.16e %1.16e %1.16e %1.16e\n', s, s+N/2, LFO, SK, a,b,Fmax);
    end
    fclose(muscle_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag); %Number of Lagrangian Pts. Total

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N/2
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif ( ( s >= N/2+1 ) && ( s < N ) )
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,L,d,Lx,Ly)

% The immsersed structure is a straight heart-tube %
x1 = [-L/2:ds:L/2 L/2];             % Constructs x-Values for bottom of tube
y1 = -d/2*ones(1,length(x1)); % Constructs y-Values for bottom of tube

x2 = x1;                      % Constructs x-Values for top of tube
y2 = -y1;                     % Constructs y-Values for top of tube

x1 = x1 + Lx/2;               % Shift into correct box
x2 = x2 + Lx/2;                % Shift into correct box

y1 = y1 + Ly/2;               % Shift into correct box
y2 = y2 + Ly/2;               % Shift into correct box

xLag = [x1 x2];
yLag = [y1 y2];

author: nickabattista


text:
%Navigation Navigation superclass
%
% An abstract superclass for implementing planar grid-based navigation classes.
%
% Methods::
%   plot        Display the occupancy grid
%   visualize   Display the occupancy grid (deprecated)
%   plan        Plan a path to goal
%   path        Return/animate a path from start to goal
%   display     Display the parameters in human readable form
%   char        Convert to string
%
%   rand        Uniformly distributed random number
%   randn       Normally distributed random number
%   randi       Uniformly distributed random integer
%
% Properties (read only)::
%   occgrid   Occupancy grid representing the navigation environment
%   goal      Goal coordinate
%   seed0     Random number state
%
% Methods that must be provided in subclass::
%   plan      Generate a plan for motion to goal
%   next      Returns coordinate of next point along path
%
% Methods that may be overriden in a subclass::
%   goal_set        The goal has been changed by nav.goal = (a,b)
%   navigate_init   Start of path planning.
%
% Notes::
% - Subclasses the MATLAB handle class which means that pass by reference semantics
%   apply.
% - A grid world is assumed and vehicle position is quantized to grid cells.
% - Vehicle orientation is not considered.
% - The initial random number state is captured as seed0 to allow rerunning an
%   experiment with an interesting outcome.
%
% See also Dstar, Dxform, PRM, RRT.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% Peter Corke 8/2009.

% TODO
%  keep dimensions of workspace in this object, have a setaxes() method
%  which transfers the dimensions to the current axes.

classdef Navigation < handle

    properties
        occgrid     % occupancy grid
        goal        % goal coordinate

        navhook     % function handle, called on each navigation iteration
        verbose     % verbosity
        seed            % current random seed
        spincount

        randstream
        seed0
    end

    % next() should be protected and abstract, but this doesnt work
    % properly
    methods (Abstract)
        plan(obj)
        n = next(obj)
    end % method Abstract

    methods

        % TODO fix up set methods for goal
        % setup argument callback like features, can we inherit from that.
        % occ grid should be an option

        % constructor

        function nav = Navigation(varargin)
        %Navigation.Navigation Create a Navigation object
        %
        % N = Navigation(OCCGRID, OPTIONS) is a Navigation object that holds an
        % occupancy grid OCCGRID.  A number of options can be be passed.
        %
        % Options::
        % 'navhook',F   Specify a function to be called at every step of path
        % 'goal',G      Specify the goal point (2x1)
        % 'verbose'     Display debugging information
        % 'inflate',K   Inflate all obstacles by K cells.
        % 'private'     Use private random number stream.
        % 'reset'       Reset random number stream.
        % 'seed',S      Set the initial state of the random number stream.  S must
        %               be a proper random number generator state such as saved in
        %               the seed0 property of an earlier run.
        %
        % Notes::
        % - In the occupancy grid a value of zero means free space and non-zero means
        %   occupied (not driveable).
        % - Obstacle inflation is performed with a round structuring element (kcircle)
        %   with radius given by the 'inflate' option.
        % - The 'private' option creates a private random number stream for the methods
        %   rand, randn and randi.  If not given the global stream is used.

            if nargin >= 1 && isnumeric(varargin{1}) && ~isscalar(varargin{1})
                nav.occgrid = varargin{1};
                varargin = varargin(2:end);
            end

            % default values of options
            opt.goal = [];
            opt.inflate = 0;
            opt.navhook = [];
            opt.private = false;
            opt.reset = false;
            opt.seed = [];

            [opt,args] = tb_optparse(opt, varargin);

            % optionally inflate the obstacles
            if opt.inflate > 0
                nav.occgrid = idilate(nav.occgrid, kcircle(opt.inflate));
            end

            % copy other options into the object
            nav.verbose = opt.verbose;
            nav.navhook = opt.navhook;
            if ~isempty(opt.goal)
                nav.goal = opt.goal(:)';
            end

            % create a private random number stream if required
            if opt.private
                nav.randstream = RandStream.create('mt19937ar');
            else
                nav.randstream = RandStream.getGlobalStream();
            end

            % reset the random number stream if required
            if opt.reset
                nav.randstream.reset();
            end

            % return the random number stream to known state if required
            if ~isempty(opt.seed)
                set(nav.randstream.set(opt.seed));
            end

            % save the current state in case it later turns out to give interesting results
            nav.seed0 = nav.randstream.State;

            nav.spincount = 0;
        end

        function r = rand(nav, varargin)
        %Navigation.rand Uniformly distributed random number
        %
        % R = N.rand() return a uniformly distributed random number from
        % a private random number stream.
        %
        % R = N.rand(M) as above but return a matrix (MxM) of random numbers.
        %
        % R = N.rand(L,M) as above but return a matrix (LxM) of random numbers.
        %
        % Notes::
        % - Accepts the same arguments as rand().
        % - Seed is provided to Navigation constructor.
        % - Provides an independent sequence of random numbers that does not
        %   interfere with any other randomised algorithms that might be used.
        %
        % See also Navigation.randi, Navigation.randn, rand, RandStream.
            r = nav.randstream.rand(varargin{:});
        end

        function r = randn(nav, varargin)
        %Navigation.randn Normally distributed random number
        %
        % R = N.randn() returns a normally distributed random number from
        % a private random number stream.
        %
        % R = N.randn(M) as above but returns a matrix (MxM) of random numbers.
        %
        % R = N.randn(L,M) as above but returns a matrix (LxM) of random numbers.
        %
        %
        % Notes::
        % - Accepts the same arguments as randn().
        % - Seed is provided to Navigation constructor.
        % - Provides an independent sequence of random numbers that does not
        %   interfere with any other randomised algorithms that might be used.
        %
        % See also Navigation.rand, Navigation.randi, randn, RandStream.
            r = nav.randstream.randn(varargin{:});
        end

        function r = randi(nav, varargin)
        %Navigation.randi Integer random number
        %
        % I = N.randi(RM) returns a uniformly distributed random integer in the
        % range 1 to RM from a private random number stream.
        %
        % I = N.randi(RM, M) as above but returns a matrix (MxM) of random integers.
        %
        % I = N.randn(RM, L,M) as above but returns a matrix (LxM) of random integers.
        %
        %
        % Notes::
        % - Accepts the same arguments as randn().
        % - Seed is provided to Navigation constructor.
        % - Provides an independent sequence of random numbers that does not
        %   interfere with any other randomised algorithms that might be used.
        %
        % See also Navigation.rand, Navigation.randn, randi, RandStream.
            r = nav.randstream.randi(varargin{:});
        end

        % invoked whenever the goal is set
        function set.goal(nav, goal)

            if ~isempty(nav.occgrid) && nav.occgrid( goal(2), goal(1)) == 1
                error('Navigation: cant set goal inside obstacle');
            end

            goal = goal(:);
            if ~(all(size(goal) == size(nav.goal)) && all(goal == nav.goal))
                % goal has changed
                nav.goal = goal(:);
                nav.goal_change();
            end
        end

        function goal_change(nav)
            %Navigation.goal_change Notify change of goal
            %
            % Invoked when the goal property of the object is changed.  Typically this
            % is overriden in a subclass to take particular action such as invalidating
            % a costmap.
        end

        function pp = path(nav, start)
            %Navigation.path Follow path from start to goal
            %
            % N.path(START) animates the robot moving from START (2x1) to the goal (which is a
            % property of the object).
            %
            % N.path() as above but first displays the occupancy grid, and prompts the user to
            % click a start location.
            % the object).
            %
            % X = N.path(START) returns the path (2xM) from START to the goal (which is a property of
            % the object).
            %
            % The method performs the following steps:
            %
            %  - Get start position interactively if not given
            %  - Initialize navigation, invoke method N.navigate_init()
            %  - Visualize the environment, invoke method N.plot()
            %  - Iterate on the next() method of the subclass until the goal is
            %    achieved.
            %
            % See also Navigation.plot, Navigation.goal.

            % if no start point given, display the map, and prompt the user to select
            % a start point
            if nargin < 2
                % display the world
                nav.plot();

                % prompt the user to click a goal point
                fprintf('** click a starting point ');
                [x,y] = ginput(1);
                fprintf('\n');
                start = round([x;y]);
            end
            start = start(:);

            % if no output arguments given, then display the world
            if nargout == 0
                % render the world
                nav.plot();
                hold on
            end

            nav.navigate_init(start);

            p = [];
            % robot is a column vector
%             if nav.backProp()==1
%                 % subclass algorithm calls for back propagation
%                 robot = nav.goal(:);
%             else
%                 robot = start;
%             end
robot = start;

            % iterate using the next() method until we reach the goal
            while true
                if nargout == 0
                    plot(robot(1), robot(2), 'g.', 'MarkerSize', 12);
                    drawnow
                end

                % move to next point on path
                robot = nav.next(robot);

                % are we there yet?
                if isempty(robot)
                    % yes, exit the loop
                    break
                else
                    % no, append it to the path
                    p = [p; robot(:)'];
                end

                % invoke the navhook function
                if isa(nav.navhook, 'function_handle')
                    nav.navhook(nav, robot(1), robot(2));
                end
            end

            % only return the path if required
            if nargout > 0
                pp = p;
            end
        end

        function visualize(nav, varargin)
            warning('visualize method deprecated for Navigation classes, use plot instead');
            nav.plot(varargin{:});
        end

        function plot(nav, varargin)
        %Navigation.plot  Visualize navigation environment
        %
        % N.plot() displays the occupancy grid in a new figure.
        %
        % N.plot(P) as above but overlays the points along the path (2xM) matrix.
        %
        % Options::
        %  'goal'         Superimpose the goal position if set
        %  'distance',D   Display a distance field D behind the obstacle map.  D is
        %                 a matrix of the same size as the occupancy grid.
        %
        % Notes::
        % - The distance field at a point encodes its distance from the goal, small
        %   distance is dark, a large distance is bright.  Obstacles are encoded as
        %   red.

            opt.goal = false;
            opt.distance = [];

            [opt,args] = tb_optparse(opt, varargin);

            clf
            if isempty(opt.distance)
                % create color map for free space + obstacle:
                %   free space, color index = 1, white,
                %   obstacle, color index = 2, red
                cmap = [1 1 1; 1 0 0];  % non obstacles are white
                image(nav.occgrid+1, 'CDataMapping', 'direct');
                colormap(cmap)

            else
                % create color map for distance field + obstacle:
                %   obstacle, color index = 1, red
                %   free space, color index > 1, greyscale

                % find maximum distance, ignore infinite values in
                % obstacles
                d = opt.distance(isfinite(opt.distance));
                maxdist = max(d(:)) + 1;

                % create the color map
                cmap = [1 0 0; gray(maxdist)];

                % ensure obstacles appear as red pixels
                opt.distance(nav.occgrid > 0) = 0;

                % display it with colorbar
                image(opt.distance+1, 'CDataMapping', 'direct');
                set(gcf, 'Renderer', 'Zbuffer')
                colormap(cmap)
                colorbar
            end

            % label the grid
            set(gca, 'Ydir', 'normal');
            xlabel('x');
            ylabel('y');
            grid on
            hold on

            if ~isempty(args)
                p = args{1};
                if numcols(p) ~= 2
                    error('expecting Nx2 matrix of points');
                end
                plot(p(:,1), p(:,2), 'g.', 'MarkerSize', 12);
            end

            if ~isempty(nav.goal) && opt.goal
                plot(nav.goal(1), nav.goal(2), 'bd', 'MarkerFaceColor', 'b');
            end
            hold off
        end

        function navigate_init(nav, start)
            %Navigation.navigate_init Notify start of path
            %
            % N.navigate_init(start) is called when the path() method is invoked.
            % Typically overriden in a subclass to take particular action such as
            % computing some path parameters. start is the initial position for this
            % path, and nav.goal is the final position.
        end

        function display(nav)
        %Navigation.display Display status of navigation object
        %
        % N.display() displays the state of the navigation object in
        % human-readable form.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a Navigation object and the command has no trailing
        %   semicolon.
        %
        % See also Navigation.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( nav.char() );
        end % display()

        function s = char(nav)
        %Navigation.char Convert to string
        %
        % N.char() is a string representing the state of the navigation
        % object in human-readable form.
            s = [class(nav) ' navigation class:'];
            s = char(s, sprintf('  occupancy grid: %dx%d', size(nav.occgrid)));
            if ~isempty(nav.goal)
                if length(nav.goal) == 2
                    s = char(s, sprintf('  goal: (%d,%d)', nav.goal) );
                else
                    s = char(s, sprintf('  goal: (%g,%g, %g)', nav.goal) );

                end
            end
        end

        function verbosity(nav, v)
        %Navigation.verbosity Set verbosity
        %
        % N.verbosity(V) set verbosity to V, where 0 is silent and greater
        % values display more information.
            nav.verbose = v;
        end

        % called at each point on the path as
        %   navhook(nav, robot)
        %
        % can be used for logging data, animation, etc.
        function navhook_set(nav, navhook)
            nav.navhook = navhook;
        end

        function message(nav, varargin)
        %Navigation.message Print debug message
        %
        % N.message(S) displays the string S if the verbose property is true.
        %
        % N.message(FMT, ARGS) as above but accepts printf() like semantics.
            if nav.verbose
                fprintf([class(nav) ' debug:: ' sprintf(varargin{:}) '\n']);
            end
        end

        function spinner(nav)
        %Navigation.spinner Update progress spinner
        %
        % N.spinner() displays a simple ASCII progress spinner, a rotating bar.
            spinchars = '-\|/';
            nav.spincount = nav.spincount + 1;
            fprintf('\b%c', spinchars( mod(nav.spincount, length(spinchars))+1 ) );
        end

    end % method

end % classde

author: Allopart

repo: rbpf-gmapping


text:
function[h]=letterlabels(arg1,arg2,arg3)
%LETTERLABELS	For automatically putting letter labels on subplots.
%
%   LETTERLABELS puts letter labels '(a)','(b)','(c)' etc., on the
%   subplots of the current figure, in the upper left-hand corner of
%   each subplot window.
%
%   LETTERLABELS(I) specifies the labelling position, clockwise from
%   top left:
%
%		I=1	Top left
%		I=2	Top right
%		I=3	Bottom right
%		I=4	Bottom left
%
%   LETTERLABELS(H,I), where H is a vector of handles to subplots,
%   puts labels on the subplots indicated by H.
%
%   LETTERLABELS(H,I,'d'), begins the labelling with letter 'd'.
%
%   LETTERLABELS will not put the letters on the right subplots if you
%   click on them before hand (because this changes the order of the
%   subplot handles).  Also if you want to label the subplots in a
%   nonstandard order, do this by reordering H.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2021 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(arg1,'--t')
    return
end

py=0.08;
ar=get(gca,'plotboxaspectratio');
xstretch=ar(1)./ar(2);
px=py/xstretch*1.8;

firstlet=real('a');
i=1;
axhand=flipud(get(gcf,'children'));

bool=false(size(axhand));
for i=1:length(bool)
    bool(i)=strcmp(axhand(i).Type,'axes');
end
axhand=axhand(bool);

for j=1:nargin
    xx=eval(['arg' int2str(j)]);
    isaxhandle=false;
    if ishandle(xx)
        if strcmp(get(xx,'Type'),'axes')
            isaxhandle=true;
        end
    end
    if ischar(xx)&~isaxhandle
        firstlet=real(xx);
    elseif (length(xx)==1)&&~isaxhandle
        i=xx;
    else
        axhand=xx;
    end
end
nax=length(axhand);
fact=1/100;

for j=1:length(axhand)
    axes(axhand(j))
    ax=axis;
    isrevx=strcmpi(get(gca,'xdir'),'reverse');
    isrevy=strcmpi(get(gca,'ydir'),'reverse');

    if ~isrevx
        x1=ax(1);
        x2=ax(2);
    else
        x1=ax(2);
        x2=ax(1);
    end
    if ~isrevy
        y1=ax(3);
        y2=ax(4);
    else
        y1=ax(4);
        y2=ax(3);
    end

    if i==1
        t=y2-(y2-y1)*fact;
        b=y2-(y2-y1)*py;
        l=x1+(x2-x1)*fact;
        r=x1+(x2-x1)*px;
    elseif i==2
        t=y2-(y2-y1)*fact;
        b=y2-(y2-y1)*py;
        l=x2-(x2-x1)*px;
        r=x2-(x2-x1)*fact;
    elseif i==3
        t=y1+(y2-y1)*py;
        b=y1+(y2-y1)*fact;
        l=x2-(x2-x1)*px;
        r=x2-(x2-x1)*fact;
    elseif i==4
        t=y1+(y2-y1)*py;
        b=y1+(y2-y1)*fact;
        l=x1+(x2-x1)*fact;
        r=x1+(x2-x1)*px;
    end

    %Tried setting log to lin before adding label but didn't help
    %if ~strcmpi(get(gca,'yscale'),'log')
    %    h(j)=patch([l l r r],[b t t b],'w');
    %    set(h(j),'edgecolor',[1 1 1])
    %end
    cx=l/2+r/2;
    cy=t/2+b/2;
    text(l+(r-l)*0.2,cy,['(' char(firstlet-1+j),')']);
    axis(ax)
end

if nargout==0
    clear h
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% op_ampScale.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_ampScale(in,A);
%
% DESCRIPTION:
% Scale the amplitude of a spectrum by factor A.
%
% INPUTS:
% in    = input data in matlab structure format
% A     = Amplitude scaling factor.
%
% OUTPUTS:
% out   = Output following amplitude scaling.

function out=op_ampScale(in,A);

out=in;
out.specs=in.specs*A;
out.fids=in.fids*A;

author: CIC-methods

repo: FID-A

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


text:
classdef HaloOrbitApproximator < matlab.mixin.SetGet
    %HaloOrbitApproximator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        gm1(1,1) double
        gm2(1,1) double
        sma2(1,1) double
        radius1(1,1) double
        radius2(1,1) double
        AzUnscaled(1,1) double
        LPt(1,1) string = "L1"
        side(1,1) string = "northern"

        %diff corr options
        diffCorAlpha(1,1) double = 0.01;
        adjustIndEnum(1,1) HaloCoordAdjustEnum = HaloCoordAdjustEnum.Z;
    end

    properties(Access=public, Dependent)
        T2(1,1) double
        LU(1,1) double
        TU(1,1) double
        VU(1,1) double
        mu(1,1) double
        gL(1,1) double
        AzScaled(1,1) double
    end

    properties(Transient)
        corrSol
        manifoldLine(1,1)
        rVectCR3BP_Prim_Manifold(3,1) double
        vVectCR3BP_Prim_Manifold(3,1) double
    end

    events
        UpdateCalcWaitbar
    end

    methods
        function obj = HaloOrbitApproximator(gm1, gm2, sma2, radius1, radius2, AzUnscaled, LPt, side)
            obj.gm1 = gm1;
            obj.gm2 = gm2;
            obj.sma2 = sma2;
            obj.radius1 = radius1;
            obj.radius2 = radius2;
            obj.AzUnscaled = AzUnscaled;
            obj.LPt = string(upper(LPt));
            obj.side = string(lower(side));
        end

        function val = getT2(obj)
            val = computePeriod(obj.sma2, obj.gm1);
        end

        function val = get.LU(obj)
            val = obj.sma2;
        end

        function val = get.TU(obj)
            val = 1/sqrt((obj.gm1 + obj.gm2)/obj.sma2^3);
        end

        function val = get.VU(obj)
            val = obj.LU / obj.TU;
        end

        function val = get.mu(obj)
            val = obj.gm2/(obj.gm1 + obj.gm2);
        end

        function val = get.gL(obj)
            LP = HaloOrbitApproximator.lagrangePoints(obj.mu);

            switch obj.LPt
                case "L1"
                    val = abs(1 - obj.mu - LP(1,1));
                case "L2"
                    val = abs(LP(2,1) - 1 + obj.mu);
                otherwise
                    error('Unknown Lagrange point.');
            end
        end

        function val = get.AzScaled(obj)
            val = (obj.AzUnscaled/obj.LU)/obj.gL;
        end

        function computeManifold(obj, tFrac, whichManifold, propDurDimensioned, hAx)
            sol = obj.corrSol;
            if(not(isempty(sol)))
                T = sol.x';
                Y = sol.y';

                stateF_Monodromy = reshape(Y(end,7:end),6,6);
                [V,D] = eig(stateF_Monodromy);
                eigvalues = diag(D,0);
                [~, IMaxEig] = max(real(eigvalues));
                [~, IMinEig] = min(real(eigvalues));

                Vs = V(:, IMinEig);
                Vu = V(:, IMaxEig);

                tFracActual = tFrac * (T(end) - T(1));
                Ynd_corr_Tfrac = deval(sol,tFracActual)';
                stateF_STM = reshape(Ynd_corr_Tfrac(1,7:end),6,6);

                Vs_i = stateF_STM * Vs;
                Vu_i = stateF_STM * Vu;

                epsilson = [0; 0; 0; 1E-4; 1E-4; 1E-4];
                Xi = Ynd_corr_Tfrac(1,1:6);
                XiS = Xi(:) + epsilson .* normVector(Vs_i);
                XiU = Xi(:) + epsilson .* normVector(Vu_i);

                propDur = propDurDimensioned / obj.TU;
                if(propDur > 0)
                    switch whichManifold
                        case HaloArriveDepartManifoldEnum.Arrive
                            [T_manifold, Y_manifold] = obj.simulateOrbit(XiU(:)', propDur, false);

                        case HaloArriveDepartManifoldEnum.Depart
                            [T_manifold, Y_manifold] = obj.simulateOrbit(XiS(:)', -propDur, false);

                        otherwise
                            error('Unknown manifold type.');
                    end

                else
                    T_manifold = 0;
                    Y_manifold = XiS(:)';
                end

                Y_manifold_Pos_Dim = Y_manifold(:,1:3) * obj.LU;
                if(isempty(obj.manifoldLine) || not(isgraphics(obj.manifoldLine, 'Line')))
                    obj.manifoldLine = plot3(hAx, Y_manifold_Pos_Dim(:,1), Y_manifold_Pos_Dim(:,2), Y_manifold_Pos_Dim(:,3));
                    obj.manifoldLine.DisplayName = 'Arrival/Departure Orbit';
                else
                    obj.manifoldLine.XData = Y_manifold_Pos_Dim(:,1);
                    obj.manifoldLine.YData = Y_manifold_Pos_Dim(:,2);
                    obj.manifoldLine.ZData = Y_manifold_Pos_Dim(:,3);
                end

                primary = [-obj.mu, 0, 0];

                obj.rVectCR3BP_Prim_Manifold = (Y_manifold(1,1:3)-primary)' * obj.LU;
                obj.vVectCR3BP_Prim_Manifold = Y_manifold(1,4:6)' * obj.VU;
            end
        end

        function [rVectCR3BP_Prim, vVectCR3BP_Prim, period, approxAzAchieved, corrAzAchieved] = computeCorrectedInitialState(obj, hAx, dispPrimBody, dispSecBody)
            [~,~,~,~,~,~, statesND_approx, ~, Tapprox] = obj.getApproxHaloOrbitInitState(0);
            state0ND_diffCor = obj.differentialCorrector(statesND_approx(1,:), Tapprox);

            [Tnd_uncorr,Ynd_uncorr] = obj.simulateOrbit(statesND_approx(1,:), Tapprox, false);
            [Tnd_corr,Ynd_corr, ~,~,~, sol_corr] = obj.simulateOrbit(state0ND_diffCor, 2*Tapprox, false);
            obj.corrSol = sol_corr;

            approxAzAchieved = max(abs(statesND_approx(:,3))*obj.LU);
            corrAzAchieved = max(abs(Ynd_corr(:,3))*obj.LU);

            fprintf('Approx Az Achieved: %0.3f km\n', approxAzAchieved);
            fprintf('Corrected Az Achieved: %0.3f km\n', corrAzAchieved);

            if(isempty(hAx))
                hAx = axes(figure());
            else
                cla(hAx);
            end

            [hTrajApprox, hLPts] = obj.plotHaloOrbit(hAx, statesND_approx);
            [hTrajUnCorr, ~] = obj.plotHaloOrbit(hAx, Ynd_uncorr);
            [hTrajCorr,   ~] = obj.plotHaloOrbit(hAx, Ynd_corr);

            legendObjs = [hTrajApprox, hTrajUnCorr, hTrajCorr];

            hTrajApprox.DisplayName = 'Approx Halo';
            hTrajUnCorr.DisplayName = 'Uncorrected Halo';
            hTrajCorr.DisplayName   = 'Corrected Halo';
            hLPts.DisplayName = 'L Points';

            axis(hAx,'equal');
            hold(hAx,'on');

            %Plot secondary
            if(dispSecBody)
                dRad = obj.radius2;
                [X,Y,Z] = sphere(30);
                X = dRad*X + (1-obj.mu)*obj.LU;
                CData = getCDataForSphereWithColormap(Z, 'gray');
                legendObjs(end+1) = surf(hAx, X,dRad*Y,dRad*Z, 'CData',CData, 'BackFaceLighting','lit', 'FaceLighting','gouraud', 'EdgeLighting','gouraud', 'LineWidth',0.1, 'EdgeAlpha',0.1, 'DisplayName','Primary Body');
            end

            %Plot primary
            if(dispPrimBody)
                dRad = obj.radius1;
                [X,Y,Z] = sphere(30);
                X = dRad*X + (-obj.mu)*obj.LU;
                CData = getCDataForSphereWithColormap(Z, 'gray');
                legendObjs(end+1) = surf(hAx, X,dRad*Y,dRad*Z, 'CData',CData, 'BackFaceLighting','lit', 'FaceLighting','gouraud', 'EdgeLighting','gouraud', 'LineWidth',0.1, 'EdgeAlpha',0.1, 'DisplayName','Primary Body');
            end

%             {'Approx Halo', 'Uncorrected Halo', 'Corrected Halo', 'Primary Body', 'Secondary Body', 'L Points'};
            legend(legendObjs, 'Location','best');

            primary = [-obj.mu, 0, 0];

            rVectCR3BP_Prim = (Ynd_corr(1,1:3)-primary)' * obj.LU;
            vVectCR3BP_Prim = Ynd_corr(1,4:6)' * obj.VU;
            period = Tapprox * obj.TU;
        end

        function state0ND_approx = differentialCorrector(obj, state0ND_approx, Tapprox)
            MU = obj.mu;
            adjustInd = obj.adjustIndEnum.adjustInd;

            deltaXDot = Inf;
            deltaZDot = Inf;
            tol = 1E-12;
            while(abs(deltaXDot) > tol || ...
                    abs(deltaZDot) > tol)
                % 1. Propagate the initial state vector using an ODE solver, such as ode45 in Matlab,
                % until the position in y is equal to zero again (T/2).
                [Tnd_uncorr,Ynd_uncorr, te,ye,ie] = obj.simulateOrbit(state0ND_approx, 10*Tapprox, true);
                stateF = Ynd_uncorr(end,:);

                % 2. Find the error in x and z velocities at T/2 (x_ and z_).
                deltaXDot = -stateF(4);
                deltaZDot = -stateF(6);
                deltaVect = [deltaXDot;
                    deltaZDot];

                % 3. Calculate the change in the initial state required to reduce the error. It is only
                % necessary to change two of the initial states. Since all the zero terms are correct
                % those are left alone, and x and y_ are the only terms that will need to be altered.
                % The initial position in z will remain xed. x0 and y_0 are the values that will
                %be used to change the initial state.
                stateF_kinematic = stateF(1:6);
                stateF_STM = stateF(7:end);
                stateF_STM = reshape(stateF_STM,6,6);

                A = [stateF_STM(4,adjustInd), stateF_STM(4,5);
                    stateF_STM(6,adjustInd), stateF_STM(6,5)];

                b = [stateF_STM(2,adjustInd), stateF_STM(2,5)];

                accelVec = HaloOrbitApproximator.cr3bpODE(Tnd_uncorr(end),stateF_kinematic,MU);
                accelVec = [accelVec(4); accelVec(6)];

                vy = stateF_kinematic(5);

                RHS = (A - (1/vy)*accelVec*b);
                deltas = RHS \ deltaVect;
                deltas = deltas * obj.diffCorAlpha;

                % 4. Now, a new initial state is calculated by adding x0 and y_0 to the original
                %initial state as in Equation 3.59 and Equation 3.60.
                r0_new = state0ND_approx(adjustInd) + deltas(1);
                ydot0_new = state0ND_approx(5) + deltas(2);

                % 5. The new estimate for the initial state vector is now:
                state0ND_approx(adjustInd) = r0_new;
                state0ND_approx(5) = ydot0_new;

                curVal = log10(max(abs(deltaVect)));
                evtData = matlab.ui.eventdata.ValueChangedData(curVal, log10(tol));
                notify(obj, 'UpdateCalcWaitbar', evtData);
            end
        end

        function [x,y,z,xdot,ydot,zdot, statesND, statesDim, Tapprox] = getApproxHaloOrbitInitState(obj, phi0)
            [k, wp, nu, ...
                a21, a22, a23, a24, a31, a32, ...
                b21, b22, b31, b32, ...
                d21, d31, d32, Ax, Az] = obj.getRichardsonConstants();

            if(obj.side == "northern")
                m = 1;
            elseif(obj.side == "southern")
                m = 3;
            else
                error('Unknown side parameter');
            end

            Tapprox = (2*pi)/(wp*nu);
            tau = linspace(0,Tapprox,100);

            deltaM = 2 - m;
            tau1 = wp .* tau + phi0;

            g = obj.gL;

            x = g*(a21*Ax^2 + a22*Az - Ax*cos(tau1) + (a23*Ax^2 - a24*Az^2)*cos(2*tau1) + ...
                (a31*Ax^3 - a32*Ax*(Az^2)*cos(3*tau1)));

            y = g*(k*Ax*sin(tau1) + (b21*Ax^2 - b22*Az^2)*sin(2*tau1) + (b31*Ax^3 - b32*Ax*Az^2)*sin(3*tau1));

            z = g*(deltaM*Az*cos(tau1) + deltaM*d21*Ax*Az*(cos(2*tau1) - 3) + deltaM*(d32*Az*Ax^2 - d31*Az^3)*cos(3*tau1));

            xdot = g*(wp*nu*Ax*sin(tau1) - 2*wp*nu*(a23*Ax^2 - a24*Az^2)*sin(2*tau1) - 3*wp*nu*(a31*Ax^3 - a32*Ax*Az^2)*sin(3*tau1));

            ydot = g*(wp*nu*k*Ax*cos(tau1) + 2*wp*nu*(b21*Ax^2 - b22*Az^2)*cos(2*tau1) + 3*wp*nu*(b31*Ax^3 - b32*Ax*Az^2)*cos(3*tau1));

            zdot = g*(-wp*nu*deltaM*Az*sin(tau1) - 2*wp*nu*deltaM*d21*Ax*Az*sin(2*tau1) - 3*wp*nu*deltaM*(d32*Az*Ax^2 - d31*Az^2)*sin(3*tau1));

            LP = HaloOrbitApproximator.lagrangePoints(obj.mu);
            switch obj.LPt
                case "L1"
                    deltaX = LP(1,1);
                case "L2"
                    deltaX = LP(2,1);
                otherwise
                    error('Unknown Lagrange point.');
            end

            x = x + deltaX; %need to move the x-coords back to the barycenter for integration in the synodic frame

            statesND = [x(:), y(:), z(:), xdot(:), ydot(:), zdot(:)];
            statesDim = [x(:)*obj.LU, ...
                y(:)*obj.LU, ...
                z(:)*obj.LU, ...
                xdot(:)*obj.VU, ...
                ydot(:)*obj.VU, ...
                zdot(:)*obj.VU];
        end

        function [hTraj, hLPts] = plotHaloOrbit(obj, hAx, statesND)
            LP = HaloOrbitApproximator.lagrangePoints(obj.mu);
            LP = LP * obj.LU;

            statesDim = statesND * obj.LU;

            hold(hAx, 'on');
            hTraj = plot3(hAx, statesDim(:,1), statesDim(:,2), statesDim(:,3));
            grid(hAx, 'off');
            grid(hAx, 'minor');
            %             plot3(-obj.mu, 0, 0, 'ko');
            %             hSecBody = plot3(hAx, (1-obj.mu)*obj.LU, 0, 0, 'ko');
            %             hSecBody = obj.plotBody(hAx);
            hLPts = plot3(hAx, LP(1:2,1), LP(1:2,2), LP(1:2,3), 'bo');
            %             axis(hAx,'equal');
            hold(hAx, 'off');
            xlabel(hAx,'X [km]');
            ylabel(hAx,'Y [km]');
            zlabel(hAx,'Z [km]');
        end

        function hCBodySurf = plotBody(obj, hAx)
            dRad = obj.radius2 / obj.LU;
            [X,Y,Z] = sphere(50);

            CData = getCDataForSphereWithColormap(Z, 'gray');
            hCBodySurf = surf(hAx, dRad*X,dRad*Y,dRad*Z, 'CData',CData, 'BackFaceLighting','lit', 'FaceLighting','gouraud', 'EdgeLighting','gouraud', 'LineWidth',0.1, 'EdgeAlpha',0.1);
            material(hCBodySurf,'dull');
        end

        function [T,Y, te,ye,ie, sol] = simulateOrbit(obj, y0, tf, stopOnYCrossing)
            tspan=[0 tf];

            if(numel(y0) == 6)
                y0 = [y0(:)', reshape(eye(6),1,36)];
            end

            odefun = @(t,y) HaloOrbitApproximator.stateTransMatrixODE(t,y,obj.mu);

            evtFcn = @(t,y) HaloOrbitApproximator.odeEvents(t,y, stopOnYCrossing);
            options=odeset('RelTol',1e-12,'AbsTol',1e-12,'Events',evtFcn);

            sol = ode113(odefun, tspan, y0, options);

            T = sol.x';
            Y = sol.y';
            te = sol.xe';
            ye = sol.ye';
            ie = sol.ie';
        end
    end

    methods(Access=private)
        function f = correctedHaloOrbitObjFun(obj, x, state0ND_approx, Tapprox)
            xNew = x(1);
            zNew = x(2);
            ydotNew = x(3);
            state0ND_approx([1,3,5]) = [xNew, zNew, ydotNew];

            [~,Ynd, ~,~,~] = obj.simulateOrbit(state0ND_approx, 10*Tapprox, true);
            aZAchieved = max(abs(Ynd(:,3))*obj.LU);

            f = abs((aZAchieved - obj.AzUnscaled)/abs(obj.AzUnscaled));
        end

        function [c, ceq] = correctedHaloOrbitNonlcon(obj, x, state0ND_approx, Tapprox)
            xNew = x(1);
            zNew = x(2);
            ydotNew = x(3);
            state0ND_approx([1,3,5]) = [xNew, zNew, ydotNew];

            [~,Ynd, ~,~,~] = obj.simulateOrbit(state0ND_approx, 10*Tapprox, true);

            c = [];
            ceq(1) = Ynd(end,2); %y=0
            ceq(2) = Ynd(end,3); %xdot = 0
            ceq(3) = Ynd(end,5); %zdot = 0
        end

        function [k, wp, nu, ...
                a21, a22, a23, a24, a31, a32, ...
                b21, b22, b31, b32, ...
                d21, d31, d32, Ax, Az] = getRichardsonConstants(obj)

            %             c2 = (1/obj.gL^3) * (obj.mu + ((1-obj.mu)*obj.gL^3)/(1-obj.gL)^3);
            %             c3 = (1/obj.gL^3) * (obj.mu - ((1-obj.mu)*obj.gL^4)/(1-obj.gL)^4);
            %             c4 = (1/obj.gL^3) * (obj.mu + ((1-obj.mu)*obj.gL^5)/(1-obj.gL)^5);

            [~, c2, c3, c4] = obj.computeCCoeff();

            wp = sqrt(2 - c2 + ((9*c2^2 - 8*c2)/2)^(1/2));
            k = (wp^2 + 1 + 2*c2)/(2*wp);

            d1 = ((3*wp^2)/k) * (k*(6*wp^2-1)-2*wp);
            d2 = ((8*wp^2)/k) * (k*(11*wp^2-1)-2*wp);

            a21 = (3*c3*(k^2-2))/(4*(1+2*c2));
            a22 = (3*c3)/(4*(1+2*c2));
            a23 = ((-3*c3*wp)/(4*k*d1))*(3*(k^3)*wp - 6*k*(k-wp) + 4);
            a24 = ((-3*c3*wp)/(4*k*d1))*(2+3*k*wp);

            d21 = -c3/(2*wp^2);
            d31 = (3/(64*wp^2)) * (4*c3*a24 + c4);
            d32 = (3/(64*wp^2)) * (4*c3*(a23 - d21) + c4*(4+k^2));

            b21 = ((-3*c3*wp)/(2*d1))*(3*k*wp - 4);
            b22 = (3*c3*wp)/d1;
            b31 = (3/(8*d2)) * (8*wp*(3*c3*(k*b21 - 2*a23) - c4*(2+3*k^2)) + ...
                (9*wp^2 + 1 + 2*c2)*(4*c3*(k*a23-b21)+(k*c4*(4+k^2))));
            b32 = (1/d2) * (9*wp*(c3*(k*b22 + d21 - 2*a24) - c4) + ...
                (3/8)*(9*wp^2 + 1 + 2*c2)*(4*c3*(k*a24-b22)+(k*c4)));

            a31 = (-9*wp/(4*d2)) * (4*c3*(k*a23 - b21) + k*c4*(4+k^2)) + ...
                ((9*wp^2 + 1 - c2)/(2*d2))*(3*c3*(2*a23-k*b21)+(c4*(2+3*k^2)));
            a32 = (-1/d2) * ((9/4)*wp*(4*c3*(k*a24 - b22) + k*c4) + ...
                (3/2)*(9*wp^2 + 1 - c2)*(c3*(k*b22 + d21 - 2*a24) - c4));

            s1 = (1/(2*wp*(wp*(1+k^2)-2*k))) * ((3/2)*c3*(2*a21*(k^2-2) - a23*(k^2+2) - 2*k*b21) - ...
                (3/8)*c4*(3*k^4 - 8*k^2 + 8));
            s2 = (1/(2*wp*(wp*(1+k^2)-2*k))) * ((3/2)*c3*(2*a22*(k^2-2) + a24*(k^2+2) + 2*k*b22 + 5*d21) + ...
                (3/8)*c4*(12 - k^2));

            a1 = (-3/2)*c3*(2*a21 + a23 + 5*d21) - (3/8)*c4*(12-k^2);
            a2 = (3/2)*c3*(a24 - 2*a22) + (9/8)*c4;

            l1 = a1 + 2*(wp^2)*s1;
            l2 = a2 + 2*(wp^2)*s2;

            delta = wp^2 - c2;

            Az = obj.AzScaled;
            Ax = sqrt((-l2 * Az^2 - delta)/l1);

            nu = 1 + s1*Ax^2 + s2*Az^2;
        end

        function [c1, c2, c3, c4] = computeCCoeff(obj)
            g = obj.gL;
            m = obj.mu;

            switch obj.LPt
                case "L1"
                    Ceqn = @(n) (1/((g^3)*(1-g^(n+1)))) * (m + ((-1)^n)*(1 - m)*g^(n+1));
                case "L2"
                    Ceqn = @(n) (1/((g^3)*(1+g^(n+1)))) * (((-1)^n)*m + ((-1)^n)*(1 - m)*g^(n+1));
                otherwise
                    error('Unknown Lagrange point.');
            end

            c = NaN(4,1);
            for(n=1:4)
                c(n) = Ceqn(n);
            end

            c1 = c(1);
            c2 = c(2);
            c3 = c(3);
            c4 = c(4);
        end
    end

    methods(Static)
        function [value,isterminal,direction] = odeEvents(~,y, stopOnYCrossing)
            value = y(2); %when y=0
            direction = 0;

            if(stopOnYCrossing)
                isterminal = 1;
            else
                isterminal = 0;
            end
        end

        function LP = lagrangePoints(mu)
            %For a given value of mu, this function computes the location
            %of all five libration points for the circular restricted three
            %body problem.  It returns them as equilibrium points in R^3
            %space.  Then the output is five points each with three components.
            %Each point is a column in a matrix with three rows.  The first column is
            %L1, the second L2, and so on to the last which is L5.

            %Compute the location of the libration points
            l=1-mu;

            LP = zeros(5,3);

            %L1
            p_L1=[1, 2*(mu-l), l^2-4*l*mu+mu^2, 2*mu*l*(l-mu)+mu-l, mu^2*l^2+2*(l^2+mu^2), mu^3-l^3];
            L1roots=roots(p_L1);
            %initialize L1 for loop
            L1=0;
            for i=1:5
                if (L1roots(i) > -mu) && (L1roots(i) < l)
                    L1=L1roots(i);
                end
            end
            LP(1,1) = L1;

            %L2
            p_L2=[1, 2*(mu-l), l^2-4*l*mu+mu^2, 2*mu*l*(l-mu)-(mu+l), mu^2*l^2+2*(l^2-mu^2), -(mu^3+l^3)];
            L2roots=roots(p_L2);
            %initialize L2 for loop
            L2=0;
            for i=1:5
                if (L2roots(i) > -mu) && (L2roots(i) > l)
                    L2=L2roots(i);
                end
            end
            LP(2,1) = L2;

            %L3
            p_L3=[1, 2*(mu-l), l^2-4*mu*l+mu^2, 2*mu*l*(l-mu)+(l+mu), mu^2*l^2+2*(mu^2-l^2), l^3+mu^3];
            L3roots=roots(p_L3);
            %initialize L3 for loop
            L3=0;
            for i=1:5
                if L3roots(i) < -mu
                    L3=L3roots(i);
                end
            end
            LP(3,1) = L3;

            %L4
            LP(4,1) = 0.5 - mu;
            LP(4,2) = sqrt(3)/2;

            %L5
            LP(5,1) = 0.5 - mu;
            LP(5,2) = -sqrt(3)/2;
        end

        function xdot = stateTransMatrixODE(t,xinput,mu)
            xstate = xinput(1:6);
            xstatedot = HaloOrbitApproximator.cr3bpODE(t,xstate,mu);

            PHI = reshape(xinput(7:6*6+6),6,6);

            x=xstate(1);
            y=xstate(2);
            z=xstate(3);

            x1 = -mu;
            x2 = 1-mu;

            r1 = sqrt((x-x1)^2 + y^2 + z^2);
            r2 = sqrt((x-x2)^2 + y^2 + z^2);

            OMEGA_XX = 1 + (1-mu)*((3*(x-x1)^2/r1^5) - 1/r1^3) ...
                + mu*((3*(x-x2)^2/r2^5) - 1/r2^3);

            OMEGA_YY = 1 + (1-mu)*((3*y^2/r1^5) - 1/r1^3) ...
                + mu*((3*y^2/r2^5) - 1/r2^3);

            OMEGA_XY = 3*(1-mu)*(x-x1)*y/r1^5 + 3*mu*(x-x2)*y/r2^5;

            OMEGA_XZ = 3*(1-mu)*z*(x-x1)/r1^5 + 3*mu*z*(x-x2)/r2^5;

            OMEGA_YZ = 3*(1-mu)*y*z/r1^5 + 3*mu*y*z/r2^5;

            OMEGA_ZZ = 3*(1-mu)*z^2/r1^5 - (1-mu)/r1^3 + 3*mu*z^2/r2^5 - mu/r2^3;

            A  = [0 0 0 1 0 0;
                0 0 0 0 1 0;
                0 0 0 0 0 1;
                OMEGA_XX OMEGA_XY OMEGA_XZ 0 2 0;
                OMEGA_XY OMEGA_YY OMEGA_YZ -2 0 0;
                OMEGA_XZ OMEGA_YZ OMEGA_ZZ 0 0 0];
            PHIdot = A * PHI;

            xdot = [xstatedot;reshape(PHIdot,6*6,1)];
        end

        function ydot = cr3bpODE(t,y, mu)
            %Define the distances from 1 and 2 to the s/c
            r1 = sqrt((y(1)+mu)^2 + y(2)^2 + y(3)^2);
            r2 = sqrt((y(1)-(1-mu))^2 + y(2)^2 + y(3)^2);

            %the deriative functions here
            %velocities
            ydot(1) = y(4);
            ydot(2) = y(5);
            ydot(3) = y(6);

            %accelerations
            ydot(4) = 2*y(5) + y(1) - ((1-mu)*(y(1)+mu))/r1^3 - (mu*(y(1)-1+mu))/r2^3;
            ydot(5) = -2*y(4) + y(2) - ((1-mu)*y(2))/r1^3 - mu*y(2)/r2^3;
            ydot(6) = (-(1-mu)*y(3))/r1^3 - (mu*y(3))/r2^3;

            ydot=ydot';
        end
    end
end

%         function [k, lambda, omega, ...
%                   a21, a22, a23, a24, a31, a32, ...
%                   b21, b22, b31, b32, ...
%                   d21, d31, d32, Ax, Ay, Az] = getRichardsonConstants(obj)
%
%             c2 = (1/obj.gL^3) * (obj.mu + ((1-obj.mu)*obj.gL^3)/(1-obj.gL)^3);
%             c3 = (1/obj.gL^3) * (obj.mu - ((1-obj.mu)*obj.gL^4)/(1-obj.gL)^4);
%             c4 = (1/obj.gL^3) * (obj.mu + ((1-obj.mu)*obj.gL^5)/(1-obj.gL)^5);
%
%             p = [1, 0, (c2-2), 0, -(c2-1)*(1+2*c2)];
%             p_roots = roots(p);
%             lambda = real(p_roots(real(p_roots) > 0 & imag(p_roots) == 0));
%
%             k = 2*lambda/(lambda^2 + 1 - c2);
%
%             d1 = ((3*lambda^2)/k) * (k*(6*lambda^2-1)-2*lambda);
%             d2 = ((8*lambda^2)/k) * (k*(11*lambda^2-1)-2*lambda);
%
%             a21 = (3*c3*(k^2-2))/(4*(1+2*c2));
%             a22 = (3*c3)/(4*(1+2*c2));
%             a23 = ((-3*c3*lambda)/(4*k*d1))*(3*(k^3)*lambda - 6*k*(k-lambda) + 4);
%             a24 = ((-3*c3*lambda)/(4*k*d1))*(2+3*k*lambda);
%
%             d21 = -c3/(2*lambda^2);
%             d31 = (3/(64*lambda^2)) * (4*c3*a24 + c4);
%             d32 = (3/(64*lambda^2)) * (4*c3*(a23 - d21) + c4*(4+k^2));
%
%             b21 = ((-3*c3*lambda)/(2*d1))*(3*k*lambda - 4);
%             b22 = (3*c3*lambda)/d1;
%             b31 = (3/(8*d2)) * (8*lambda*(3*c3*(k*b21 - 2*a23) - c4*(2+3*k^2)) + ...
%                                          (9*lambda^2 + 1 + 2*c2)*(4*c3*(k*a23-b21)+(k*c4*(4+k^2))));
%             b32 = (1/d2) * (9*lambda*(c3*(k*b22 + d21 - 2*a24) - c4) + ...
%                                          (3/8)*(9*lambda^2 + 1 + 2*c2)*(4*c3*(k*a24-b22)+(k*c4)));
%
%             a31 = (-9*lambda/(4*d2)) * (4*c3*(k*a23 - b21) + k*c4*(4+k^2)) + ...
%                                          ((9*lambda^2 + 1 - c2)/(2*d2))*(3*c3*(2*a23-k*b21)+(c4*(2+3*k^2)));
%             a32 = (-1/d2) * ((9/4)*lambda*(4*c3*(k*a24 - b22) + k*c4) + ...
%                                           (3/2)*(9*lambda^2 + 1 - c2)*(c3*(k*b22 + d21 - 2*a24) - c4));
%
%             s1 = (1/(2*lambda*(lambda*(1+k^2)-2*k))) * ((3/2)*c3*(2*a21*(k^2-2) - a23*(k^2+2) - 2*k*b21) - ...
%                                                         (3/8)*c4*(3*k^4 - 8*k^2 + 8));
%             s2 = (1/(2*lambda*(lambda*(1+k^2)-2*k))) * ((3/2)*c3*(2*a22*(k^2-2) + a24*(k^2+2) + 2*k*b22 + 5*d21) + ...
%                                                         (3/8)*c4*(12 - k^2));
%
%             a1 = (-3/2)*c3*(2*a21 + a23 + 5*d21) - (3/8)*c4*(12-k^2);
%             a2 = (3/2)*c3*(a24 - 2*a22) + (9/8)*c4;
%
%             l1 = a1 + 2*(lambda^2)*s1;
%             l2 = a2 + 2*(lambda^2)*s2;
%
%             delta = lambda^2 - c2;
%
%             Az = (obj.AzUnscaled/obj.LU)*obj.gL;
%             Ay = Az/obj.gL;
%             Ax = sqrt((-l2 * Az^2 - delta)/l1);
%
%             omega = 1 + s1*Ax^2 + s2*Az^2;
%         end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ DetectionReview ] = DetectionReview(obj,event )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
%obj = findobj('tag','EventExplorerMaster');
FO = guidata(obj);
%% Select the time windows to look at
%User input for this
numwins = 30; %number of windows to look at. determine to maximize sampling or have user input (FO.uinput.field?)
%Selecting from events:
%randevents = randsample(FO.EventTimes,numevents);
%Selecting from random times (RestrictInts takes way too long...)
set(findobj(FO.fig,'Type','uicontrol'),'Enable','off');
drawnow;
if ~isempty(FO.detectionints)
    [status,interval,index] = InIntervals(FO.data.lfp.timestamps,double(FO.detectionints));
    restrictedtimes = FO.data.lfp.timestamps(status);
else
    restrictedtimes = FO.data.lfp.timestamps;
end

%use InIntervals
set(findobj(FO.fig,'Type','uicontrol'),'Enable','on');
randevents = randsample(restrictedtimes,numwins);
%Find any events that are within winsize of another event to remove them?
%Also should look at windows in which no events were detected?  Maybe just
%random windows in the detectionwin
closeevents = abs(diff(randevents))<FO.winsize;
% numtries =0;
% while any(closeevents) && numtries<50
%     randevents(closeevents)
% end
%Problem:overlapping ints
%[ ints ] = MergeSeparatedInts( ints,minseparation )

%%
%Store the current user action
FO.currentuseraction = 'MarkEvents';
FO.markedevents = []; %clear any previously marked events
FO.viewmode = 'timepoint';
guidata(FO.fig, FO); %store the data in the figure

counter = numwins;
miss=[];hit=[];falsealarm=[];
lookedatwins = [];

%UI panel for event review
FO.EventPanel = uipanel('Title','Detection Review','FontSize',12,...
        'Position',[.65 .05 0.25 0.32]);
%Buttons for next loop or finishing
nextbtn = uicontrol('Parent',FO.EventPanel,...
    'Position',[200 20 150 40],'String','Next Window (Return)',...
          'Callback','uiresume(gcbf)');
qutbtn = uicontrol('Parent',FO.EventPanel,...
    'Position',[200 65 150 40],'String','Quit Early (Esc)',...
         'Callback','QUITLOOP=true;uiresume(gcbf)');
%Instruction text
instructtext = uicontrol('Parent',FO.EventPanel,...
    'Position',[25 100 400 100],'style','text',...
    'string',{'Instructions:',...
    '   - LEFT click missed events (o)','   - RIGHT click False Alarms (x)'},...
    'HorizontalAlignment','left');
%Display counter text
correcttext = uicontrol('Parent',FO.EventPanel,...
    'Position',[25 60 125 15],'style','text',...
    'string','Correct: 0','HorizontalAlignment','left');
misstext = uicontrol('Parent',FO.EventPanel,...
    'Position',[25 40 125 15],'style','text',...
    'string','Miss: 0','HorizontalAlignment','left');
FAtext = uicontrol('Parent',FO.EventPanel,...
    'Position',[25 20 125 15],'style','text',...
    'string','FA: 0','HorizontalAlignment','left');
countetext = uicontrol('Parent',FO.EventPanel,...
    'Position',[25 100 175 15],'style','text',...
    'string',['Windows Remaining: ',num2str(counter)],'HorizontalAlignment','left');

%The Event Review Loop
QUITLOOP = false;
while counter>0 && QUITLOOP~=true
    %EventUI(FO)   %function that will run the user interface
    thiseventtime = randevents(counter);
    FO.currevent = thiseventtime; guidata(FO.fig, FO);
    viewinfo = EventVewPlot;
    uiwait(FO.fig)  %Wait here until the user clicks quit or next

    lookedatwins = [lookedatwins; viewinfo.thiseventwin];
    inwinevents = viewinfo.inwinevents;

    %Tally the user selections for that window
    obj = findobj('tag','EventExplorerMaster');  FO = guidata(obj);
    if ~isempty(FO.markedevents)
        FO.markedevents(isnan(FO.markedevents(:,1)))=[]; %remove nans - bug fix
        miss = [miss; FO.markedevents(FO.markedevents(:,3)==1,1)];
        %This is messy to account for interp errors with 0-1 reference points
        if isempty(FO.markedevents(FO.markedevents(:,3)==3,1))
            newfalsealarms = [];
        elseif length(FO.markedevents(FO.markedevents(:,3)==3,1))==1 && length(inwinevents)==1
            newfalsealarms = inwinevents;
        else
            newfalsealarms = interp1(inwinevents,inwinevents,...
                FO.markedevents(FO.markedevents(:,3)==3,1),'nearest');
        end
        falsealarm = [falsealarm; newfalsealarms];
    end
    hit = [hit; inwinevents(~ismember(inwinevents,falsealarm))];
    FO.markedevents = [];  guidata(FO.fig, FO); %reset the marked events

    counter = counter-1;
    set(correcttext, 'String',['Correct: ',num2str(length(hit))]);
    set(misstext, 'String', ['Miss: ',num2str(length(miss))]);
    set(FAtext, 'String', ['FA: ',num2str(length(falsealarm))]);
    set(countetext, 'String', ['Windows Remaining: ',num2str(counter)]);
end
%Calculate total number of miss,hit,FA
numMiss = length(miss);
numHit = length(hit);
numFA = length(falsealarm);

%Calculate total amount of time/percentage of detection time (detectionintervals) looked at

%Put things in the output structure
DetectionReview.lookedatwins = lookedatwins;
DetectionReview.miss = miss;
DetectionReview.hit = hit;
DetectionReview.falsealarm = falsealarm;
DetectionReview.estMissperc = numMiss./(numHit+numMiss);
DetectionReview.estFAperc = numFA./(numHit+numFA);
DetectionReview.ReviewDate = today;
DetectionReview.EventsType = FO.EventName;

%UI: Done!  Would you like to save the results to (eventsfilename?)
%Make function that does this: SaveResults(FO,EEoutput)
if isfield(FO,'eventsfilename')
    button = questdlg(['Event Review Complete! Would you like to add the results to ',...
        FO.eventsfilename,'?'],'Good Job!');
    switch button
        case 'Yes'
            %Load the events file, add the field, save the events file
            try %Only do this if the correct named structure lives in the file
                eventsfile = load(FO.eventsfilename,FO.EventName);
                eventsfile.(FO.EventName).EventExplorer.DetectionReview = DetectionReview;
                save(FO.eventsfilename,'-struct','eventsfile',FO.EventName,'-append')
            catch
                warndlg({' Save failed... ',[FO.eventsfilename,' may not ',...
                    'contain a structure titled ',FO.EventName,'.'],...
                    'Or you may not have sudo priviliges...?'},'Oh No!')
            end
    end
end

%% Return to explorer mode
FO.viewmode = 'events';
FO.currentuseraction = 'none';
FO.currevent = 1;
set(FO.EventPanel,'Visible','off')
FO.DetectionReview = DetectionReview;
guidata(FO.fig, FO); %Save the detection review back to GUI data
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function gui_popup( jPopup, jParent, x, y )
% GUI_POPUP: Show a JPopupMenu at the current location of the mouse cursor.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2014

% Parse inputs
if (nargin < 2) || isempty(jParent)
    jParent = jPopup;
end
if (nargin < 4) || isempty(x) || isempty(y)
    mouseInfo = java.awt.MouseInfo.getPointerInfo();
    mouseLoc  = mouseInfo.getLocation();
else
    mouseLoc = java.awt.Point(x,y);
end
jPopup.pack();

% If parent is a figure handle
if ~isjava(jParent) && ~isempty(jParent)
    % Fix stupid bug in some older Matlab versions
    if any(bst_get('MatlabVersion') == [709, 710, 711])
        jParent = jPopup;
    % Matlab >= 2019b deprecated the JavaFrame property, impossible to attach JPopupMenu to a figure
    elseif (bst_get('MatlabVersion') >= 907)
        jParent = jPopup;
    else
        % Get the handle of the java figure
        jParent = bst_get('JFrame', jParent);
        % If the call crashed
        if isempty(jParent)
            jParent = jPopup;
        end
    end
end
% Show popup
jPopup.setLocation(mouseLoc);
jPopup.setInvoker(jPopup);
jPopup.setVisible(1);
drawnow;

author: brainstorm-tools


text:
function [ekg_RSlinB__ELF_Cto_SFu] = fuse_rr_is(ekg_RSlinB_am_ELF_Cto, ekg_RSlinB_bw_ELF_Cto, ekg_RSlinB_fm_ELF_Cto, up)
%UNTITLED5 Summary of this function goes here
%   Detailed explanation goes here
rel_sub_comps = up.al.sub_components.fus_mod;
% Identify relevant data
feat_mods = {'bw', 'am', 'fm'};
for mod_no = 1:3
    mod = feat_mods(mod_no);
    eval(['rel_data.' mod{1,1} ' = ekg_RSlinB_', mod{1,1}, '_ELF_Cto' ';']);
    %eval(['rel_data.' mod{1,1} ' = rrEsts.' temp.bef{fusable_rr_est_no}, mod{1,1}, temp.aft{fusable_rr_est_no} ';']);
    %eval(['rel_data.' mod{1,1} ' = rrEsts.' temp.bef{fusable_rr_est_no}, mod{1,1}, temp.aft{fusable_rr_est_no} ';']);
end

for sub_comp_no = 1 : length(rel_sub_comps)
    % Do fusion
    temp_rr = feval(rel_sub_comps{sub_comp_no}, rel_data, up); clear rel_data
    % store this series of rrs:
    save_name = ['ekg_RSlinB__ELF_Cto_SFu'];
    eval([save_name, ' = temp_rr;']);
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function result = runScriptFile(fileName)
% This function runs the test in fileName
% It can distinguish between skipped and Failed tests. A test is considered
% to be skipped if it throws a COBRA:RequirementsNotMet error.
%
% OUTPUTS:
%
%    result:    A structure array with the following fields:
%
%                - `.passed`: true if the test passed otherwise false
%                - `.skipped`: true if the test was skipped otherwise false
%                - `.failed`: true if the test failed, or was skipped, otherwise false
%                - `.status`: a string representing the status of the test ('failed','skipped' or'passed')
%                - `.fileName`: the fileName of the test
%                - `.time`: the duration of the test (if passed otherwise NaN)
%                - `.statusMessage`: Informative string about potential problems.
%                - `.Error`: Error message received from a failed or skipped test
%
% .. Author:   - Thomas Pfau Jan 2018.

global CBT_MISSING_REQUIREMENTS_ERROR_ID

COBRA_TESTSUITE_TESTFILE = fileName;

% get the timinig (and hope these values are not overwritten.
COBRA_TESTSUITE_STARTTIME = clock();

try
    % run the file
    executefile(fileName);
catch ME
    % vatch errors and interpret them
    clearvars -except ME COBRA_TESTSUITE_STARTTIME COBRA_TESTSUITE_TESTFILE CBT_MISSING_REQUIREMENTS_ERROR_ID
    scriptTime = etime(clock(), COBRA_TESTSUITE_STARTTIME);
    result = struct('status', 'failed', 'failed', true, 'passed', false, 'skipped', false, 'fileName', ...
                    COBRA_TESTSUITE_TESTFILE, 'time', scriptTime, 'statusMessage', 'fail', 'Error', ME);
    if strcmp(ME.identifier, CBT_MISSING_REQUIREMENTS_ERROR_ID)
        % requirement missing, so the test was skipped.
        result.status = 'skipped';
        result.skipped = true;
        result.failed = false;
        result.statusMessage = ME.message;
    else
        % actual error in the test.
        result.skipped = false;
        result.status = 'failed';
        result.statusMessage = ME.message;
    end
    return
end

% get the timinig.
scriptTime = etime(clock(), COBRA_TESTSUITE_STARTTIME);

result = struct('status', 'passed', 'failed', false, 'passed', true, 'skipped', false, 'fileName', ...
                COBRA_TESTSUITE_TESTFILE, 'time', scriptTime, 'statusMessage', 'success', 'Error', MException('', ''));

end

function executefile(fileName)
% runs a script file (used to separate workspaces)
    run(fileName)
end

author: opencobra

repo: cobratoolbox


text:
function [x, infos] = div_mu_partial_nmf(V, rank, in_options)
%
% This file is part of NMFLibrary
%
% Created by H.Kasai on Feb. 16, 2017
%
% Change log:
%
%       June 16, 2022 (Hiroyuki Kasai): Initial version.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.alg           = 'mu';
    local_options.norm_h        = 0;
    local_options.norm_w        = 1;
    local_options.alpha         = 2;
    local_options.delta         = 0.1;
    local_options.metric_type   = 'kl-div'; % 'kl-div' (default)
    local_options.d_alpha       = -1; % for alpha divergence
    local_options.d_beta        = 0; % for beta divergence
    local_options.myeps         = 1e-16;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    if ~strcmp(options.alg, 'mu')
        fprintf('Invalid algorithm: %s. Therfore, we use mu (i.e., multiplicative update).\n', options.alg);
        options.alg = 'mu';
    end

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = sprintf('MU-Partial (%s:%s)', options.alg, options.metric);
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % store initial info
    clear infos;
%     if strcmp(options.metric, 'alpha-div')
%         metric_param = options.d_alpha;
%     elseif strcmp(options.metric, 'beta-div')
%         metric_param = options.d_beta;
%     end
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('MU-Partial (%s:%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', options.alg, options.metric, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        if strcmp(options.alg, 'mu')
            if strcmp(options.metric, 'euc')

                % update H
                H = H .* (W' * V) ./ (W' * W * H);
                H = H + (H<options.myeps) .* options.myeps;

                % update W
                W(:, options.updateW) = W(:, options.updateW) .* (V * H(options.updateW, :)') ...
                    ./ (W * (H * H(options.updateW, :)'));
                W = W + (W<options.myeps) .* options.myeps;

            elseif strcmp(options.metric, 'kl-div')

                % update W
                W(:, options.updateW) = W(:, options.updateW) .* ...
                    ((V./(W*H + options.myeps))*H(options.updateW, :)')./(ones(m,1)*sum(H(options.updateW, :)'));
                if options.norm_w ~= 0
                    W = normalize_W(W, options.norm_w);
                end

                % update H
                H = H .* (W'*(V./(W*H + options.myeps)))./(sum(W)'*ones(1,n));
                if options.norm_h ~= 0
                    H = normalize_H(H, options.norm_h);
                end

            elseif strcmp(options.metric, 'alpha-div')

                % update W
                W(:, options.updateW) = W(:, options.updateW) .* ...
                    ( ((V+options.myeps) ./ (W*H+options.myeps)).^options.d_alpha * H(options.updateW, :)').^(1/options.d_alpha);
                if options.norm_w ~= 0
                    W = normalize_W(W, options.norm_w);
                end
                W = max(W, options.myeps);

                % update H
                H = H .* ( (W'*((V+options.myeps)./(W*H+options.myeps)).^options.d_alpha) ).^(1/options.d_alpha);
                if options.norm_h ~= 0
                    H = normalize_H(H, options.norm_h);
                end
                H = max(H, options.myeps);

            elseif strcmp(options.metric, 'beta-div')

                WH = W * H;

                % update W
                W(:, options.updateW) = W(:, options.updateW) .* ...
                    ( ((WH.^(options.d_beta-2) .* V)*H(options.updateW, :)') ./ ...
                    max(WH.^(options.d_beta-1)*H(options.updateW, :)', options.myeps) );

                if options.norm_w ~= 0
                    W = normalize_W(W, options.norm_w);
                end

                WH = W * H;

                % update H
                H = H .* ( (W'*(WH.^(options.d_beta-2) .* V)) ./ max(W'*WH.^(options.d_beta-1), options.myeps) );
                if options.norm_h ~= 0
                    H = normalize_H(H, options.norm_h);
                end
            else
                error('Invalid metric.')
            end

        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [model, solutionThermoRecon, solutionRecon, model1] = secondPassDirectionalityAssignment(model)
% Driver to call model specific code to manually generate a physiological model (if first pass does not result in a physiological model).
%
% The second pass directionality assignment needs careful manual curation
% since the adjustments necessary to get one organism to grow will not
% necessarily be the same as the ones which will get another
% organism to grow. There's no avioding manual debugging at this stage.
%
% USAGE:
%
%    [model, solutionThermoRecon, solutionRecon, model1] = secondPassDirectionalityAssignment(model)
%
% INPUTS:
%    model:
%
% OUTPUTS:
%    model:
%    solutionThermoRecon:
%    solutionRecon:
%    model1:
%
% NOTE:
%
%    This is the code  used for a number of organisms in order to point out
%    the kind of issues that arise. This is NOT supposed to work in the
%    general case.
%
% .. Author: - Ronan M. T. Fleming

global CBTLPSOLVER
switch model.description
    case 'iAF1260'
        % now assign reaction directions based on the P(\Delta_{r}G^{\primeo}<0)
        % doing so may prevent the model from growing, therefore at this stage it
        % is necessary to manually adjust some of the reaction directionalities
        % such that the model can grow, and grow at a similar rate as observed in
        % vivo. Therefore this script cannot be made model invariant as there is an
        % essential manual debugging stage. The script below is for E. coli iAF1260
        fprintf('%s\n%s\n%s\n','The second pass assignment of reaction directionality should',...
            'be a compromise between quantitative and qualitative assignment',...
            'this step requires manual imput as it is specific to each organism.');

        fprintf('%s\n',['The second pass assignment of reaction directionality is using the ' CBTLPSOLVER ' LP solver.']);
        fprintf('\n%s\n','...setThermoReactionDirectionality (based on the P(\Delta_{r}G^{\primeo}<0))');
        % this will update the model.lb_reconThermo & model.ub_reconThermo
        [model,solutionThermoRecon,solutionRecon,model1]=setThermoReactionDirectionalityiAF1260(model);
    otherwise
        fprintf('No manually generated .m file second pass directionality assignment is available for this model.')
end


text:
function p = prior_invt(varargin)
%PRIOR_INVT  Student-t prior structure for the inverse of the parameter
%
%  Description
%    P = PRIOR_INVT('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates for the inverse of the parameter Student's
%    t-distribution prior structure in which the named parameters
%    have the specified values. Any unspecified parameters are set
%    to default values.
%
%    P = PRIOR_INVT(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parameterisation is done as in Bayesian Data Analysis,
%    second edition, Gelman et.al 2004.
%
%    Parameters for Student-t prior [default]
%      mu       - location [0]
%      s2       - scale [1]
%      nu       - degrees of freedom [4]
%      mu_prior - prior for mu [prior_fixed]
%      s2_prior - prior for s2 [prior_fixed]
%      nu_prior - prior for nu [prior_fixed]
%
%  See also
%    PRIOR_T, PRIOR_*

% Copyright (c) 2000-2001,2010,2012 Aki Vehtari
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010 Jaakko Riihim�ki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_INVT';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('mu',0, @(x) isscalar(x));
  ip.addParamValue('mu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('s2',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('s2_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('nu',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('nu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Inv-t';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Inv-t')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('mu',ip.UsingDefaults)
    p.mu = ip.Results.mu;
  end
  if init || ~ismember('s2',ip.UsingDefaults)
    p.s2 = ip.Results.s2;
  end
  if init || ~ismember('nu',ip.UsingDefaults)
    p.nu = ip.Results.nu;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('mu_prior',ip.UsingDefaults)
    p.p.mu=ip.Results.mu_prior;
  end
  if init || ~ismember('s2_prior',ip.UsingDefaults)
    p.p.s2=ip.Results.s2_prior;
  end
  if init || ~ismember('nu_prior',ip.UsingDefaults)
    p.p.nu=ip.Results.nu_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_invt_pak;
    p.fh.unpak = @prior_invt_unpak;
    p.fh.lp = @prior_invt_lp;
    p.fh.lpg = @prior_invt_lpg;
    p.fh.recappend = @prior_invt_recappend;
  end

end

function [w, s] = prior_invt_pak(p)

  w=[];
  s={};
  if ~isempty(p.p.mu)
    w = p.mu;
    s=[s; 'Inv-t.mu'];
  end
  if ~isempty(p.p.s2)
    w = [w log(p.s2)];
    s=[s; 'log(Inv-t.s2)'];
  end
  if ~isempty(p.p.nu)
    w = [w log(p.nu)];
    s=[s; 'log(Inv-t.nu)'];
  end
end

function [p, w] = prior_invt_unpak(p, w)

  if ~isempty(p.p.mu)
    i1=1;
    p.mu = w(i1);
    w = w(i1+1:end);
  end
  if ~isempty(p.p.s2)
    i1=1;
    p.s2 = exp(w(i1));
    w = w(i1+1:end);
  end
  if ~isempty(p.p.nu)
    i1=1;
    p.nu = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_invt_lp(x, p)

  lJ = -log(x)*2;   % log(1/x^2) log(|J|) of transformation
  xt  = 1./x;       % transformation
  lp = sum(gammaln((p.nu+1)./2) -gammaln(p.nu./2) -0.5*log(p.nu.*pi.*p.s2) -(p.nu+1)./2.*log(1+(xt-p.mu).^2./p.nu./p.s2) + lJ);

  if ~isempty(p.p.mu)
    lp = lp + p.p.mu.fh.lp(p.mu, p.p.mu);
  end
  if ~isempty(p.p.s2)
    lp = lp + p.p.s2.fh.lp(p.s2, p.p.s2) +log(p.s2);
  end
  if ~isempty(p.p.nu)
    lp = lp + p.p.nu.fh.lp(p.nu, p.p.nu) +log(p.nu);
  end
end

function lpg = prior_invt_lpg(x, p)

  lJg = -2./x;      % gradient of log(|J|) of transformation
  xt  = 1./x;       % transformation
  xtg = -1./x.^2;   % derivative of transformation
  lpg = xtg.*(-(p.nu+1).* (xt-p.mu) ./ (p.nu.*p.s2 + (xt-p.mu).^2)) + lJg;

  if ~isempty(p.p.mu)
    lpgmu = sum( (p.nu+1).* (xt-p.mu) ./ (p.nu.*p.s2 + (xt-p.mu).^2) ) + p.p.mu.fh.lpg(p.mu, p.p.mu);
    lpg = [lpg lpgmu];
  end
  if ~isempty(p.p.s2)
    lpgs2 = (sum( -1./(2.*p.s2) +((p.nu + 1)*(p.mu - xt)^2)./(2*p.s2*((p.mu-xt)^2 + p.nu*p.s2))) + p.p.s2.fh.lpg(p.s2, p.p.s2)).*p.s2 + 1;
    lpg = [lpg lpgs2];
  end
  if ~isempty(p.p.nu)
    lpgnu = (0.5*sum( digamma1((p.nu+1)./2)-digamma1(p.nu./2)-1./p.nu-log(1+(xt-p.mu).^2./p.nu./p.s2)+(p.nu+1)./(1+(xt-p.mu).^2./p.nu./p.s2).*(xt-p.mu).^2./p.s2./p.nu.^2) + p.p.nu.fh.lpg(p.nu, p.p.nu)).*p.nu + 1;
    lpg = [lpg lpgnu];
  end
end

function rec = prior_invt_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.mu)
    rec.mu(ri,:) = p.mu;
  end
  if ~isempty(p.p.s2)
    rec.s2(ri,:) = p.s2;
  end
  if ~isempty(p.p.nu)
    rec.nu(ri,:) = p.nu;
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Let Cpn_est=(I-S(e))Cpn. Then e=[Egyro Eplat]*[gyro error;Platform error]
%platform err def:Cbp_est=(I-S(p))*Cbp
%gyro err def: gyro_meas=gyro+gyro_error
function [Cnp, Egyro, Eplat]=alingc_hd(gyro, Lat, Cbp)

WIE_E=7292115e-11;
nrot=WIE_E*cos(Lat);

%normalize gyro;
gyro=gyro/norm(gyro)*WIE_E;

%Compute the earth rotation in platform frame
if isempty(Cbp)
    Cbp=eye(3);
end
gyro_p=Cbp*gyro;

%Transformation
Cnp=zeros(3);
Cnp(1,1)=gyro_p(1)/nrot;
Cnp(2,2)=gyro_p(1)/nrot;
Cnp(2,1)=gyro_p(2)/nrot;
Cnp(1,2)=-Cnp(2,1);
Cnp(3,3)=1;

%Error matrices
Egyro=[-gyro_p(2)/nrot^2 gyro_p(1)/nrot^2 0]*Cbp;
Eplat=[-gyro_p(2)/nrot^2 gyro_p(1)/nrot^2 0]*skew(gyro_p);

repo: ekfmonoslam


text:
function planC = warp_dose(deformS,doseCreationScanNum,movDoseNum,movPlanC,planC,tmpDirPath,inverseFlag)
% function planC = warp_dose(deformS,doseCreationScanNum,movPlanC,planC)
%
% APA, 07/19/2012

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

if ~exist('tmpDirPath','var') || isempty(tmpDirPath)
    tmpDirPath = fullfile(getCERRPath,'ImageRegistration','tmpFiles');
end

indexMovS = movPlanC{end};
indexS = planC{end};

optS = getCERROptions();
plmFlag = 0; elxFlag = 0; antsFlag = 0;
% Create b-spline coefficients file
if isstruct(deformS)
    if ~inverseFlag
        baseScanUID = deformS.baseScanUID;
        movScanUID  = deformS.movScanUID;
    else
        baseScanUID = deformS.movScanUID;
        movScanUID  = deformS.baseScanUID;
    end
%     algorithm = deformS.algorithm;
    registration_tool = deformS.registrationTool;
    switch upper(registration_tool)
        case 'PLASTIMATCH'
            plmFlag = 1;
            disp('Plastimatch selected');
            bspFileName = fullfile(tmpDirPath,['bsp_coeffs_',baseScanUID,'_',movScanUID,'.txt']);
            success = write_bspline_coeff_file(bspFileName,deformS.algorithmParamsS);
        case 'ELASTIX'
            elxFlag = 1;
            disp('Elastix selected');
        case 'ANTS'
            antsFlag = 1;
            disp('ANTs selected');
    end

else
    bspFileName = deformS;
    indexS = planC{end};
    indexMovS = movPlanC{end};
    movScanNum = getDoseAssociatedScan(movDoseNum,movPlanC);
    movScanUID = movPlanC{indexMovS.scan}(movScanNum).scanUID;
    baseScanUID = planC{indexS.scan}(doseCreationScanNum).scanUID;
end

% Convert dose to .mha
movDoseUID = movPlanC{indexMovS.dose}(movDoseNum).doseUID;
randPart = floor(rand*1000);
movDoseUniqName = [movDoseUID,num2str(randPart)];
movDoseFileName = fullfile(tmpDirPath,['movDose_',movDoseUniqName,'.mha']);

% Write .mha file for this dose
success = createMhaDosesFromCERR(movDoseNum, movDoseFileName, movPlanC);

%write target space reference mha scan
refScanNum = findScanByUID(planC,baseScanUID);
[refScanUniqName, ~] = genScanUniqName(planC,refScanNum);
refScanFileName = fullfile(tmpDirPath,['baseScan_',refScanUniqName,'.mha']);
success = createMhaScansFromCERR(refScanNum, refScanFileName, planC);

% Generate name for the output .mha file
warpedMhaFileName = fullfile(tmpDirPath,['warped_dose_',baseScanUID,'_',movScanUID,'.mha']);

if plmFlag
    % Switch to plastimatch directory if it exists
    prevDir = pwd;
    plmCommand = 'plastimatch warp ';
    optName = fullfile(getCERRPath,'CERROptions.json');
    optS = opts4Exe(optName);
    if exist(optS.plastimatch_build_dir,'dir') && isunix
        cd(stateS.optS.plastimatch_build_dir)
        plmCommand = ['./',plmCommand];
    end

    % Issue plastimatch warp command with nearest neighbor interpolation
    %fail = system([plmCommand, '--input ', movDoseFileName, ' --output-img ',
    %warpedMhaFileName, ' --xf ', bspFileName, ' --algorithm itk']); %
    %ITK-based warping, consider adding it an an option.
    fail = system([plmCommand, '--input ', movDoseFileName, ' --output-img ',...
        warpedMhaFileName, ' --xf ', bspFileName, ' --fixed ', refScanFileName]);
    if fail % try escaping slashes
        system([plmCommand, '--input ', escapeSlashes(movDoseFileName),...
            ' --output-img ', escapeSlashes(warpedMhaFileName),...
            ' --xf ', escapeSlashes(bspFileName),...
            ' --fixed ', escapeSlashes(refScanFileName)])
    end
end

if antsFlag
            % ANTs path setup
        if ~exist(optS.antspath_dir,'dir')
            error(['ANTSPATH ' optS.antspath_dir ' not found on filesystem. Please review CERROptions.']);
        end
        antspath = fullfile(optS.antspath_dir,'bin');
        setenv('ANTSPATH', antspath);
        antsScriptPath = fullfile(optS.antspath_dir, 'Scripts');
        antsCERRScriptPath = fullfile(getCERRPath,'CERR_core','ImageRegistration','antsScripts');
        if isunix
            setenv('PATH',[antspath ':' antsScriptPath ':' antsCERRScriptPath ':' getenv('PATH')]);
            antsCommand = '';
        else
            setenv('PATH',[antspath ';' antsScriptPath ';' antsCERRScriptPath ';' getenv('PATH')]);
            antsCommand = '';
        end
        % generate transform command
        transformParams = buildAntsTransform(deformS.algorithmParamsS.antsWarpProducts);
        antsCommand = [antsCommand ' antsApplyTransforms -d 3 -r ' refScanFileName ' -i ' movDoseFileName ' -o ' warpedMhaFileName ' ' transformParams]; % ' ' interpParams];
        disp(antsCommand);
        system(antsCommand);
        disp('Warp complete, importing warped scan to planC...');

end

% Read the warped output .mha file within CERR
infoS  = mha_read_header(warpedMhaFileName);
data3M = mha_read_volume(infoS);
% [data3M,infoS] = readmha(warpedMhaFileName);
doseName = movPlanC{indexMovS.dose}(movDoseNum).fractionGroupID;
assocScanUID = planC{indexS.scan}(doseCreationScanNum).scanUID;
planC = dose2CERR(flipdim(permute(data3M,[2,1,3]),3),[],...
    ['Warped_',doseName],[],[],'UniformCT',[],'no',...
    assocScanUID,planC);

% Cleanup
try
    delete(movDoseFileName)
    delete(warpedMhaFileName)
    %delete(bspFileName)
end
try
    if isstruct(deformS)
        delete(bspFileName)
    end
end

% Switch back to the previous directory
% cd(prevDir)

% Get DICOMHeader for the original dose
dcmHeaderS = [];
if isfield(movPlanC{indexMovS.dose}(movDoseNum),'DICOMHeaders')
    dcmHeaderS = movPlanC{indexMovS.dose}(movDoseNum).DICOMHeaders;
end
if ~isempty(dcmHeaderS)
    dcmHeaderS.PixelSpacing = [];
    dcmHeaderS.ImagePositionPatient = [];
    dcmHeaderS.ImageOrientationPatient = [];
    dcmHeaderS.GridFrameOffsetVector = [];
    dcmHeaderS.Rows = [];
    dcmHeaderS.Columns = [];
    dcmHeaderS.SliceThickness = [];
end
planC{indexS.dose}(end).DICOMHeaders = dcmHeaderS;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function disp(F)
%DISP   Display a CHEBFUN2V.
%
% See also DISPLAY.

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

loose = strcmp(get(0,'FormatSpacing'),'loose');

% Compact version
if ( isempty( F ) )
    fprintf('empty chebfun2v\n')
    return
end

if ( F.isTransposed )
    tString = 'Row vector';
else
    tString = 'Column vector';
end

disp(['   chebfun2v object ' '(' tString ') containing:' ])
if ( loose )
    fprintf('\n');
end

% Display its two CHEBFUN2 halves.
for j = 1:F.nComponents
    disp( F.components{j} );
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [I] = corrindex(XData, YData, YDataM, par_number)
%
% Function for calculation of the correlation index
%
% Inputs:   XData:  x data
%           YData:  y data
%           YDataM: model data
%           par_number: number of model parameters
% Output:   I: correlation index
%
% Authors:
% Ivo Petras (ivo.petras@tuke.sk)
% Dagmar Bednarova (dagmar.bednarova@tuke.sk)
%
% Date: 21/002/2007
%
kx=length(XData);
ky=length(YData);
kym=length(YDataM);
if kx ~= ky
    disp('Incompatible X and Y data.');
    close all;
end
n=kym;
sey=(sum((YData-YDataM).^2))/(n-par_number);
sy=var(YData);
I=(1-(sey./sy))^0.5;
%

author: Sable

repo: mcbench-benchmarks


text:
function gd=gabprojdual(gm,g,a,M,varargin);
%GABPROJDUAL   Gabor Dual window by projection
%   Usage:  gd=gabprojdual(gm,g,a,M)
%           gd=gabprojdual(gm,g,a,M,L)
%
%   Input parameters:
%         gm    : Window to project.
%         g     : Window function.
%         a     : Length of time shift.
%         M     : Number of modulations.
%         L     : Length of transform to consider
%   Output parameters:
%         gd    : Dual window.
%
%   `gabprojdual(gm,g,a,M)` calculates the dual window of the Gabor frame given
%   by *g*, *a* and *M* closest to *gm* measured in the $l^2$ norm. The
%   function projects the suggested window *gm* onto the subspace of
%   admissable dual windows, hence the name of the function.
%
%   `gabprojdual(gm,g,a,M,L)` first extends the windows *g* and *gm* to
%   length *L*.
%
%   `gabprojdual(...,'lt',lt)` does the same for a non-separable lattice
%   specified by *lt*. Please see the help of |matrix2latticetype| for a
%   precise description of the parameter *lt*.
%
%   See also:  gabdual, gabtight, gabdualnorm, fir2long

%   AUTHOR : Peter L. Søndergaard

if nargin<4
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
definput.flags.phase={'freqinv','timeinv'};
[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);

%% ------ step 2: Verify a, M and L
if isempty(L)
    % Minimum transform length by default.
    Ls=1;

    % Use the window lengths, if any of them are numerical
    if isnumeric(g)
        Ls=max(length(g),Ls);
    end;

    if isnumeric(gm)
        Ls=max(length(gm),Ls);
    end;

    % ----- step 2b : Verify a, M and get L from the window length ----------
    L=dgtlength(Ls,a,M,kv.lt);

else

    % ----- step 2a : Verify a, M and get L

    Luser=dgtlength(L,a,M,kv.lt);
    if Luser~=L
        error(['%s: Incorrect transform length L=%i specified. Next valid length ' ...
               'is L=%i. See the help of DGTLENGTH for the requirements.'],...
              upper(mfilename),L,Luser)
    end;

end;

[g, info_g]  = gabwin(g, a,M,L,kv.lt,'callfun',upper(mfilename));
[gm,info_gm] = gabwin(gm,a,M,L,kv.lt,'callfun',upper(mfilename));

% gm must have the correct length, otherwise dgt will zero-extend it
% incorrectly using postpad instead of fir2long
gm=fir2long(gm,L);

% Calculate the canonical dual.
gamma0=gabdual(g,a,M,'lt',kv.lt);

% Get the residual
gres=gm-gamma0;

% Calculate parts that lives in span of adjoint lattice.
if isreal(gres) && isreal(gamma0) && isreal(g) && kv.lt(2)<=2
    gk=idgtreal(dgtreal(gres,gamma0,M,a,'lt',kv.lt),g,M,a,'lt',kv.lt)*M/a;
else
    gk=idgt(dgt(gres,gamma0,M,a,'lt',kv.lt),g,M,'lt',kv.lt)*M/a;
end;

% Construct dual window
gd=gamma0+(gres-gk);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function mv = mv_analysisMenu(mv, hfig);
%
% mv = mv_analysisMenu(mv, hfig);
%
% Add menus with analysis options for
% the MultiVoxel UI.
%
%
% ras,  04/05
if ieNotDefined('hfig')
    hfig = gcf;
end

if ieNotDefined('mv')
    mv = get(hfig, 'UserData');
end

mv.ui.analMenu = uimenu('ForegroundColor',  'k',  'Label',  'Analysis',  'Separator',  'on');

% split-half option (odd v even trials)
uimenu(mv.ui.analMenu, 'Label', 'Split-Half analysis', ...
       'Separator', 'off', ...
       'Callback', 'mv_splitHalf;');

% split-half option (odd v even trials), remove means
uimenu(mv.ui.analMenu, 'Label', 'Split-Half analysis (remove means)', ...
       'Separator', 'off', ...
       'Callback', 'mv_splitHalf([], 1);');

% % add bootstrapping option (odd v even trials)
% uimenu(mv.ui.analMenu, 'Label', 'Reliability Analysis (odd vs even)', ...
%        'Separator', 'off', ...
%        'Callback', 'results = mv_reliability;');

% % add bootstrapping option (choose subsets)
% uimenu(mv.ui.analMenu, 'Label', 'Reliability Analysis (select subsets)', ...
%        'Separator', 'off', ...
%        'Callback', 'results = mv_reliability([], ''user'');');

% % sort by voxR option
% uimenu(mv.ui.analMenu, 'Label', 'Sort Voxels by Voxel Reliability', ...
%    'Separator', 'on', 'Callback', 'mv_sortByVoxR;');
%
% % sort by omniR option
% uimenu(mv.ui.analMenu, 'Label', 'Sort Voxels by Omnibus Correlation', ...
%    'Separator', 'off', 'Callback', 'mv_sortByOmniR;');
%
% % sort by MI option
% uimenu(mv.ui.analMenu, 'Label', 'Sort Voxels by Mutual Information', ...
%    'Separator', 'off', 'Callback', 'mv_sortByMutualInf;')

% sort by voxel amplitude
uimenu(mv.ui.analMenu, 'Label', 'Sort Voxel Order (sortrows)', ...
   'Separator', 'on', 'Callback', 'mv_sortVoxels([],  ''sortrows'',  ''dialog'');');

% sort by dprime
uimenu(mv.ui.analMenu, 'Label', 'Sort Voxel Order (d'')', ...
   'Separator', 'off', 'Callback', 'mv_sortVoxels([],  ''dprime'');');

% sort by GLM variance explained
uimenu(mv.ui.analMenu, 'Label', 'Sort Voxel Order (GLM variance explained)', ...
   'Separator', 'off', 'Callback', 'mv_sortVoxels([],  ''varexplained'');');

% clustering analysis option
uimenu(mv.ui.analMenu, 'Label', 'Cluster Voxels', ...
   'Separator', 'off', 'Callback', 'mv_basicCluster;');

% blur time courses option
uimenu(mv.ui.analMenu, 'Label', 'Temporal Blur Each Voxel''s Time Course', ...
   'Separator', 'on', ...
   'Callback', 'mv_blurTimeCourse;');

% sub-select voxels option
uimenu(mv.ui.analMenu,  'Label',  'Select Voxel Subset Based on Criteria...',  ...
    'Separator',  'on',  'Callback',  'mv_selectSubset([],  [],  ''dialog''); multiVoxelUI; ');

% remove outliers option
uimenu(mv.ui.analMenu, 'Label', 'Remove Outlier Voxels', ...
   'Separator', 'off', ...
   'Callback', 'mv_removeOutliers;');

% pick a subset option
uimenu(mv.ui.analMenu, 'Label', 'Click on voxels for a new UI', ...
   'Separator', 'off', 'Callback', 'mv_pickVoxels;');

% make a time course UI option
uimenu(mv.ui.analMenu, 'Label', 'Time Course UI for these voxels', ...
   'Separator', 'on', 'Callback', 'mv_makeTCUI;');

% % dump data to workspace option
% uimenu(mv.ui.analMenu, 'Label', 'Dump Data to Workspace', ...
%    'Separator', 'on', 'Callback', 'tc_dumpDataToWorkspace;');
%

set(hfig, 'UserData', mv);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pstruct = tapas_ehgf_jget_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2017-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

l = length(pvec)/8;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

pstruct.mux_0     = pvec(1:l);
pstruct.sax_0     = pvec(l+1:2*l);
pstruct.mua_0     = pvec(2*l+1:3*l);
pstruct.saa_0     = pvec(3*l+1:4*l);
pstruct.kau       = pvec(4*l+1);
pstruct.kax       = pvec(4*l+2:5*l);
pstruct.kaa       = pvec(5*l+1:6*l-1);
pstruct.omu       = pvec(6*l);
pstruct.omx       = pvec(6*l+1:7*l);
pstruct.oma       = pvec(7*l+1:8*l);

end

author: translationalneuromodeling

repo: tapas


text:
function [m, pot] = marginal_nodes(engine, query)
% MARGINAL_NODES Compute the marginal on the specified set of nodes (global_joint)
% [m, pot] = marginal_nodes(engine, query)

pot = marginalize_pot(engine.jpot, query);
m = pot_to_marginal(pot);
%[m.T, lik] = normalize(m.T);
%loglik = log(lik);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt