text:
%% Execute this file to collect training data for model identification

%% Parameters: Model
n = 5;
options = odeset('RelTol',1e-10,'AbsTol',1e-10*ones(1,n));
SIM_DURATION = 10;  % Physical time simulation is integrated before integrator is switched, adjust if necessary, depends on length of time series to be computed
dt = 1/12;       % time step of data

run_HIV_params

% Steady-state of model corresponding to progressive infection
xSTEADY1 = zeros(1,5);
xSTEADY1(2) = ((c2*(lambda1-d*q)-b2*alpha1) - ...
    sqrt((c2*(lambda1-d*q)-b2*alpha1)^2 - 4*alpha1*c2*q*d*b2))/(2*alpha1*c2*q);
xSTEADY1(1) = lambda1/(d+alpha1*xSTEADY1(2));
xSTEADY1(4) = 0;
xSTEADY1(5) = (xSTEADY1(2)*c2*(alpha1*q-a) + b2*alpha1)/(c2*p2*xSTEADY1(2));
xSTEADY1(3) = h*xSTEADY1(5)/(c2*q*xSTEADY1(2));

% Steady-state of model corresponding to successful immune response
xSTEADY2 = zeros(1,5);
xSTEADY2(1) = (lambda1*c1)/(d*c1 + b1*alpha1);
xSTEADY2(2) = b1/c1;
xSTEADY2(3) = 0;
xSTEADY2(4) = (alpha1*xSTEADY2(1)-a)/p1;
xSTEADY2(5) = 0;

xref = xSTEADY1; % recovered healthy steady state as reference, might be used in system identification
%% Collect data
if DATA_ENSEMBLE == 1

    Nic = 32; %568
    if exist(fullfile(datapath,['DATA_',SystemModel,'_TRAINING-ENSEMBLE_',InputSignalType,'_N',num2str(Nic),'.mat']))==2
        load(fullfile(datapath,['DATA_',SystemModel,'_TRAINING-ENSEMBLE_',InputSignalType,'_N',num2str(Nic),'.mat']))
    else
        tspan=[0:dt:20];
        Ntrain = (length(tspan)-1)/2+1;
        Nt = length(tspan);

        if strcmp(ModelName,'SINDYc')
            [x10, x20, x30, x40, x50] = ndgrid([1000],[0,10], [0,10], [0.1,1,10], [0.1,1,10]);
        elseif strcmp(ModelName,'NARX') || strcmp(ModelName,'DMDc')
            [x10, x20, x30, x40, x50] = ndgrid([10,10],[0.1,0.1], [0.1,0.1],[0.1,0.1], [0.1,0.1]);

        end
        [N1,N2,N3,N4,N5] = size(x10);
        x0_ensemble = [reshape(x10,[N1*N2*N3*N4*N5,1]), reshape(x20,[N1*N2*N3*N4*N5,1]), reshape(x30,[N1*N2*N3*N4*N5,1]), reshape(x40,[N1*N2*N3*N4*N5,1]), reshape(x50,[N1*N2*N3*N4*N5,1])];
        Nic = size(x0_ensemble,1)

        xensemble = zeros(Nt,n,Nic); % Init ensemble data
        switch InputSignalType

            case 'sine2'

                rng(1,'twister')
                Nrand = [1*rand(Nic,5)];
                A = 2;
                forcing = @(x,t) [(0.4*(sin(2*pi*0.2*t).*sin(2*pi*0.05*t)))+0.3];
                for iIC = 1:Nic
                    tic
                    try
                        forcing = @(x,t) [(Nrand(iIC,1)*A* (sin(Nrand(iIC,1)*0.7*t).*sin(Nrand(iIC,1)*.1*t).*sin(Nrand(iIC,1)*.2*t).*sin(Nrand(iIC,1)*.05*t)) ).^2]; %[(A*(sin(0.01*t)+sin(.1*t))).^2];
                        [t,x]=ode45(@(t,x) HIVsys_KWON(t,x,forcing(x,t)),tspan,x0_ensemble(iIC,:),options);
                        xensemble(:,:,iIC) = x;
                        for i = 1:length(tspan)
                            u(i,:,iIC) = forcing(0,tspan(i));
                        end
                    catch
                        disp(['ERROR: Simulation failed for i=',num2str(iIC)])
                        xensemble(:,:,iIC) = nan(Nt,n);
                        u(:,:,iIC) = zeros(Nt,length(A));
                    end
                    tend = toc;
                    disp(['Time for ',num2str(iIC), ' of ',num2str(Nic),': ', num2str(tend)])
                end
                figure,plot(t,xensemble(:,:,iIC))
                figure,plot(t,squeeze(xensemble(:,5,:)))

            case 'prbs'
                A = 1;
                taulim = [0.2 8];
                states = [0,0:0.25:1,0,0,0]
                Nswitch = 200;

                rng(1,'twister');
                seed = randi(Nic,Nic,1); % vary actuation in each trajectory
                forcing = @(x,t,seedval) [A(1)*prbs(taulim, Nswitch, states, t,0,seedval)];
                u = zeros(length(tspan),length(A),Nic);

                for iIC = 1:Nic
                    tic
                    try
                        % IF simulation takes longer than SIM_DURATION,
                        % stop simulation and switch solver
                        [t,x,isterminal] = integrateODE(tspan, x0_ensemble(iIC,:), 'ode45', SIM_DURATION,@(x,t) forcing(x,t,seed(iIC)));
                        if isterminal == 1 % try different solver
                            [t,x,isterminal] = integrateODE(tspan, x0_ensemble(iIC,:), 'ode15s', SIM_DURATION,@(x,t) forcing(x,t,seed(iIC)));
                        end
                        xensemble(:,:,iIC) = x;
                        for i = 1:length(tspan)
                            u(i,:,iIC) = forcing(0,tspan(i),seed(iIC));
                        end
                    catch
                        disp(['ERROR: Simulation failed for i=',num2str(iIC)])
                        xensemble(:,:,iIC) = nan(Nt,n);
                        u(:,:,iIC) = zeros(Nt,length(A));
                    end
                    tend = toc;
                    disp(['Time for ',num2str(iIC), ' of ',num2str(Nic),': ', num2str(tend)])
                end

        end

        %% Clean up data
        % from unstable simulations or which failed
        IXnaninf = zeros(Nic,1);
        for i = 1:Nic
            IXnan = isnan(squeeze((xensemble(:,:,i))));
            IXinf = isinf(squeeze(abs(xensemble(:,:,i))));
            if any(IXnan(:)==1) || any(IXinf(:)==1)
                IXnaninf(i) = 1;
            end
        end
        xensemble(:,:,logical(IXnaninf)) = [];
        u(:,:,logical(IXnaninf)) = [];
        Nic = size(xensemble,3);

        %% Split into training and validation data set
        Ntrain = ceil(length(tspan)/2);
        xv = xensemble(Ntrain+1:end,:,:);
        x = xensemble(1:Ntrain,:,:);

        uv = u(Ntrain+1:end,:,:);
        u = u(1:Ntrain,:,:);

        tv = t(Ntrain+1:end);
        t = t(1:Ntrain);

        tspanv = tspan(Ntrain+1:end);
        tspan = tspan(1:Ntrain);

        %% Show data

        figure;
        subplot(6,1,1)
        plot(tspan,squeeze(x(:,1,:)),'LineWidth',1.5)
        ylabel('x1')
        set(gca,'LineWidth',1, 'FontSize',14)

        subplot(6,1,2)
        plot(tspan,squeeze(x(:,2,:)),'LineWidth',1.5)
        ylabel('x2')
        set(gca,'LineWidth',1, 'FontSize',14)

        subplot(6,1,3)
        plot(tspan,squeeze(x(:,3,:)),'LineWidth',1.5)
        ylabel('x3')
        set(gca,'LineWidth',1, 'FontSize',14)

        subplot(6,1,4)
        plot(tspan,squeeze(x(:,4,:)),'LineWidth',1.5)
        ylabel('x4')
        set(gca,'LineWidth',1, 'FontSize',14)

        subplot(6,1,5)
        plot(tspan,squeeze(x(:,5,:)),'LineWidth',1.5)
        ylabel('x5')
        set(gca,'LineWidth',1, 'FontSize',14)

        subplot(6,1,6)
        plot(tspan,squeeze(u(:,1,:)),'LineWidth',1.5)
        ylabel('u')
        set(gca,'LineWidth',1, 'FontSize',14)
        %%
        T = length(tspan);
        N = length(tspan);
        save(fullfile(datapath,['DATA_',SystemModel,'_TRAINING-ENSEMBLE_',InputSignalType,'_N',num2str(Nic),'.mat']))
    end
elseif DATA_ENSEMBLE == 0
    rng(3,'twister')

    % Initial Condition
    x0 = [10, 0.1, 0.1, 1, 0.1];
    tspan=[0:dt:400];
    Ntrain = (length(tspan)-1)/2+1;

    switch InputSignalType
        case 'unforced'
            forcing = @(x,t) [0];
            [t,x]=ode45(@(t,x) HIVsys_ZURAKOWSKI(t,x,forcing(x,t)),tspan,x0,options);
            u = zeros(length(tspan),1);
            for i = 1:length(tspan)
                u(i,:) = forcing(0,tspan(i));
            end
        case 'sine2'
            A = 1.5;
            Nic = 1; iIC = 1;
            Nrand = [rand(Nic,1),rand(Nic,1),randn(Nic,1),rand(Nic,1),randn(Nic,1),rand(Nic,1),randn(Nic,1),rand(Nic,1),randn(Nic,1)];
            forcing = @(x,t) [mod((A* (sin(Nrand(iIC,2)*0.7*(t-Nrand(iIC,3))).* ...
                sin(Nrand(iIC,4)*.1*(t-Nrand(iIC,5))).*sin(Nrand(iIC,6)*.2*(t-Nrand(iIC,7))).*sin(Nrand(iIC,8)*.05*(t-Nrand(iIC,9)))) ).^2,1)];
            [t,x]=ode45(@(t,x) HIVsys_ZURAKOWSKI(t,x,forcing(x,t)),tspan,x0,options);
            u = forcing(0,tspan)';

        case 'prbs'
            A = 1; % PAPER
            taulim = [0.2 10];
            states = [0:0.25:1,0,0,0,0];
            Nswitch = 200;
            forcing = @(x,t) [A(1)*prbs(taulim, Nswitch, states, t,0)];

            [t,x]=ode45(@(t,x) HIVsys_ZURAKOWSKI(t,x,forcing(x,t)),tspan,x0,options);

            u = zeros(length(tspan),1);
            for i = 1:length(tspan)
                u(i,:) = forcing(0,tspan(i));
            end
            figure,plot(tspan,u)
            size(u)

    end

    %% Split into training and validation data set
    Ntrain = ceil(length(tspan)/2);
    xv = x(Ntrain+1:end,:);
    x = x(1:Ntrain,:);

    uv = u(Ntrain+1:end,:);
    u = u(1:Ntrain,:);

    tv = t(Ntrain+1:end);
    t = t(1:Ntrain);

    tspanv = tspan(Ntrain+1:end);
    tspan = tspan(1:Ntrain);

    %% Show results
    betta_eff = alpha1*(1-eta*u(1));
    betta_eff<c1*(c2*b1*(lambda1-q*d)-b2*c1*d)/(b1*(c2*b1*q+b2*c1))
    xSTEADY = xSTEADY1;
    figure;
    subplot(5,1,1),hold on
    % plot(t,xSTEADY(1)*ones(size(t)),'--k')
    plot(t,x(:,1),'LineWidth',1.5)
    ylabel('x1')
    set(gca,'LineWidth',1, 'FontSize',14)
    axis tight

    subplot(5,1,2),hold on
    % plot(t,xSTEADY(2)*ones(size(t)),'--k')
    plot(t,x(:,2),'LineWidth',1.5)
    ylabel('x2')
    set(gca,'LineWidth',1, 'FontSize',14)
    axis tight

    subplot(5,1,3),hold on
    % plot(t,xSTEADY(3)*ones(size(t)),'--k')
    plot(t,x(:,3),'LineWidth',1.5)
    ylabel('x3')
    set(gca,'LineWidth',1, 'FontSize',14)
    axis tight

    subplot(5,1,4),hold on
    % plot(t,xSTEADY(4)*ones(size(t)),'--k')
    plot(t,x(:,4),'LineWidth',1.5)
    ylabel('x4')
    set(gca,'LineWidth',1, 'FontSize',14)
    axis tight

    subplot(5,1,5),hold on
    % plot(t,xSTEADY(5)*ones(size(t)),'--k')
    plot(t,x(:,5),'LineWidth',1.5)
    ylabel('x5')
    set(gca,'LineWidth',1, 'FontSize',14)
    axis tight

    T = length(tspan);
    N = length(tspan);
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function [g]= spm_gx_adem_write(x,v,a,P)
% returns the prediction for a two-joint arm (proprioception and vision)
% FORMAT [g]= spm_gx_adem_write(x,v,a,P)
%
% x    - hidden states:
%   x(1) - joint angle
%   x(2) - joint angle
%   x(3) - angular velocity
%   x(4) - angular velocity
% v    - causal states{
%   v(1) - exogenous force (x)
%   v(2) - exogenous force (y)
% a    - action
% P    - parameters
%
% g    - sensations:
%   g(1) - joint angle (proprioception)
%   g(2) - joint angle (proprioception)
%   g(3) - arm location (visual)
%   g(4) - arm location (visual)
%
% As for spm_dem_reach but with no visual target
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_gx_adem_write.m 3901 2010-05-27 16:14:36Z karl $

% evaluate positions
%--------------------------------------------------------------------------
J  = spm_dem_reach_x2J(x);

% stretch (angular) and visual (positional) information about motor plant
%==========================================================================
g  = [x; J{1}; J{1} + J{2}];

author: spm


text:
function manageMiji(varargin)
	% Controls starting and stopping of Miji in a way to reduce issues
	% Biafra Ahanonu
	% started: 2019.04.08
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.06.20 [00:20:38] - Add support for setting up ImageJ along with closing all windows to future proof any changes to those calls.
        % 2021.07.16 [13:34:57] - Check that ImageJ already in java path to prevent duplicate loading and wasting time.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2021.12.19 [10:40:47] - Added `forceLoad` option.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Str:
		% 'start' - start Miji'.
		% 'setupImageJ' - Sets up ImageJ by pointing Java path to the jar files.
		% 'closeAllWindows' - Closes all open windows but leaves ImageJ running.
		% 'exit' - close Miji.
	options.startStop = 'start';
    % Binary: 1 = force loading Miji.
    options.forceLoad = 0;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		startStop = options.startStop;
		switch startStop
			case 'start'
				try
					% If MIJ class not loaded, load Miji.m.
					if exist('MIJ','class')~=8
						% Load Miji so paths added to javaclasspath('-dynamic')
						currP = pwd;
						% Miji;
						Miji(false);
						cd(currP);
						MIJ.exit;
					end
				catch err
					disp(repmat('@',1,7))
					disp(getReport(err,'extended','hyperlinks','on'));
					disp(repmat('@',1,7))

					disp('Reset Java class path and Miji then try again');
					resetMiji

					% Try again after resetting Miji.
					try
						MIJ.start;
					catch err
						disp(repmat('@',1,7))
						disp('Apparently Miji hates your computer, sorry!')
						disp(getReport(err,'extended','hyperlinks','on'));
						disp(repmat('@',1,7))
						return;
					end
				end

				% If Miji.m not in path, ask user
				if exist('MIJ','class')~=8
				else
					if exist('Miji.m','file')~=2
						modelAddOutsideDependencies('miji');
					end
				end

				% First attempt to open Miji
				try
					MIJ.start;
				catch err
					disp(repmat('@',1,7))
					disp(getReport(err,'extended','hyperlinks','on'));
					disp(repmat('@',1,7))

					disp('Reset Java class path and Miji then try again');
					resetMiji

					% Try again after resetting Miji.
					try
						MIJ.start;
					catch err
						disp(repmat('@',1,7))
						disp('Apparently Miji hates your computer, sorry!')
						disp(getReport(err,'extended','hyperlinks','on'));
						disp(repmat('@',1,7))
					end
				end
			case 'setupImageJ'
				% Sets up ImageJ by pointing Java path to the jar files.

				imagejPath = [ciapkg.getDirExternalPrograms filesep 'imagej'];
				pathsToAdd = {[imagejPath filesep 'mij.jar'],[imagejPath filesep 'ij.jar']};
				nPaths = length(pathsToAdd);

                javaDyna = javaclasspath('-dynamic');
                % matchIdx = cellfun(@isempty,regexpi(javaDyna,pathsToAdd));
                matchIdx = ismember(pathsToAdd,javaDyna);
                % cellfun(@(x) javarmpath(x),javaDyna(matchIdx));
                % javaDynaPathStr = join(javaDyna(matchIdx),''',''');
                if options.forceLoad==1||any(~matchIdx)
                    for i = 1:nPaths
                        if matchIdx(i)==1
                            continue;
                        end
                        thisPath = pathsToAdd{i};
                        disp('Loading MIJI + ImageJ.')
                        fprintf('Adding to Java path: %s\n',thisPath);
                        javaaddpath(thisPath);
                    end
                    % eval(sprintf('javarmpath(''%s'');',javaDynaPathStr{1}))
                else
                    disp('ImageJ already loaded!')
                end
			case 'closeAllWindows'
				% Closes all open windows but leaves ImageJ running.

				% MIJ.run('Close All Without Saving');
				% MIJ.closeAllWindows;
				allClosedFlag = 0
				nImagesToBreak = 20;
				imageNo = 1
				while allClosedFlag==0
					try
						MIJ.getListImages
						MIJ.run('Close')
					catch
						allClosedFlag = 1;
					end
					imageNo = imageNo + 1;
				end
				MIJ.closeAllWindows;
				% MIJ.exit;
			case 'exit'
				try
					MIJ.exit;
				catch err
					disp(repmat('@',1,7))
					disp('Miji likely already closed, throwing Java error.')
					disp(getReport(err,'extended','hyperlinks','on'));
					disp(repmat('@',1,7))
				end
			otherwise
				% Body
		end
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function view3d(arg,arg2)

% view3d  Interactively rotate, zoom and pan the view of a 3-D plot
% --------------------------------------------------------------------
%
% VIEW3D ROT turns on mouse-based 3-D rotation
% VIEW3D ZOOM turns on mouse-based 3-D zoom and pan
% VIEW3D OFF turns it off
%
% VIEW3D(FIG,...) works on the figure FIG
%
% Double click to restore the original view
%
% hit "z" key over the figure to switch from ROT to ZOOM
% hit "r" key over the figure to switch from ZOOM to ROT
%
% in ROT mode:
% press and hold left mouse button to rotate about screen xy axis
% press and hold middle mouse button to rotate about screen z axis
% in ZOOM mode:
% press and hold left mouse button to zoom in and out
% press and hold middle mouse button to move the plot
%
% --------------------------------------------------------------------
% inspired from rotate3d by The MathWorks, Inc.
%
% Torsten Vogel 09.04.1999
% tv.volke@bmw.de
% tested under Matlab 5.2
% --------------------------------------------------------------------

% ---------------------------------------------- inputs --------------
if nargin == 0
   error('not enough inputs')
elseif nargin == 1
   if ishandle(arg)
      error('not enough inputs')
      return
   else
      switch(lower(arg))
      case 'rot'
         viewact(gcf,'rot')
      case 'zoom'
         viewact(gcf,'zoom')
      case 'off'
         viewact(gcf,'off')
      case 'down'
         view3dDownFcn
      case 'up'
         view3dUpFcn
      case 'keypress'
         view3dkeypressFcn
      case 'view_xy' % rotate via screen xy axis
         view3dxyFcn
      case 'view_z'  % rotate via screen z axis
         view3dzFcn
      case 'view_zoom' % zoom in and out
         view3dzoomFcn
      case 'view_pan'  % move the plot
         view3dpanFcn
      otherwise
         error('misspelled command argument')
      end
   end
elseif nargin==2
   if ~ishandle(arg)
      error('bad figure handle')
   end
   switch(lower(arg2))
   case 'rot'
      viewact(arg,'rot')
   case 'zoom'
      viewact(arg,'zoom')
   case 'off'
      viewact(arg,'off')
   otherwise
      error('misspelled command argument')
   end
end

% ---------------------------------------------- activation ----------
function viewact(fig,what)

% de-/activates view3d for the given figure

view3dObj = findobj(allchild(fig),'Tag','view3dObj');

if strcmp(what,'rot')
   if isempty(view3dObj)
      view3dObj = makeview3dObj(fig); %the small text box at the lower left corner
   end
   vdata = get(view3dObj,'UserData');
   vdata.what = 'rot';
   set(view3dObj,'UserData',vdata);
elseif strcmp(what,'zoom')
   if isempty(view3dObj)
      view3dObj = makeview3dObj(fig); %the small text box at the lower left corner
   end
   vdata = get(view3dObj,'UserData');
   vdata.what = 'zoom';
   set(view3dObj,'UserData',vdata);
elseif strcmp(what,'off')
   if isempty(view3dObj)
      return
   end
   vdata = get(view3dObj,'UserData');
   uirestore(vdata.uistate);
   set(fig,'KeyPressFcn',vdata.oldkeypressfcn)
   delete(view3dObj);
end

% ---------------------------------------------- view3dDownFcn -------
function view3dDownFcn

view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
mouseclick = get(gcf,'SelectionType');
if isempty(view3dObj)
   return
end
vdata = get(view3dObj,'UserData');
vdata.oldunits = get(gcf,'Units');
set(gcf,'Units','pixels');
vdata.old_pt = get(0,'PointerLocation');
%  ----------------- store or restore previous view
ViewData = get(get(gca,'zlabel'),'UserData');
if isempty(ViewData)
   ViewData = manageViewData('get_from_axes');
   set(get(gca,'zlabel'),'UserData',ViewData)
end
if strcmp(mouseclick,'open')
   manageViewData('set_axes',ViewData);
   set(gcf,'Units',vdata.oldunits)
   return
end
%  ----------------- display text box
fig_color = get(gcf,'Color');
c = sum([.3 .6 .1].*fig_color);
set(vdata.textbox,'BackgroundColor',fig_color);
if(c > .5)
   set(vdata.textbox,'ForegroundColor',[0 0 0]);
else
   set(vdata.textbox,'ForegroundColor',[1 1 1]);
end
%  ----------------- what to do?
if strcmp(vdata.what,'rot')
   if strcmp(mouseclick,'normal')
      set(vdata.textbox,'string','Screen XY Rotation');
      set(gcf,'WindowButtonMotionFcn','view3d(''view_xy'')');
      set(gcf,'Pointer','custom','pointershapecdata',pointershapes('rot'));
   elseif strcmp(mouseclick,'extend')
      set(vdata.textbox,'string','Screen Z Rotation');
      set(gcf,'WindowButtonMotionFcn','view3d(''view_z'')');
      set(gcf,'Pointer','custom','pointershapecdata',pointershapes('rot'));
   end
else
   if strcmp(mouseclick,'normal')
      set(vdata.textbox,'string','Zoom');
      set(gcf,'WindowButtonMotionFcn','view3d(''view_zoom'')');
      set(gcf,'Pointer','custom','pointershapecdata',pointershapes('zoom'));
   elseif strcmp(mouseclick,'extend')
      set(vdata.textbox,'string','Pan');
      set(gcf,'WindowButtonMotionFcn','view3d(''view_pan'')');
      set(gcf,'Pointer','custom','pointershapecdata',pointershapes('pan'));
   end
end
set(view3dObj,'UserData',vdata)
set(vdata.textbox,'visi','on')

% ---------------------------------------------- view3dUpFcn ---------
function view3dUpFcn

view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
if isempty(view3dObj)
   return
end
vdata = get(view3dObj,'UserData');
set(gcf,'WindowButtonMotionFcn','','Units',vdata.oldunits,'pointer','arrow')
set(view3dObj,'visi','off')

% ---------------------------------------------- view3dkeypressFcn ---
function view3dkeypressFcn

view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
if isempty(view3dObj)
   return
end
vdata = get(view3dObj,'UserData');
currchar = lower(get(gcf,'currentchar'));
if strcmp(currchar,'r')
   vdata.what = 'rot';
elseif strcmp(currchar,'z')
   vdata.what = 'zoom';
end
set(view3dObj,'UserData',vdata)

% ---------------------------------------------- view3dxyFcn ---------
function view3dxyFcn

view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
vdata = get(view3dObj,'UserData');
new_pt = get(0,'PointerLocation');
old_pt = vdata.old_pt;
dx = (new_pt(1) - old_pt(1))*.5;
dy = (new_pt(2) - old_pt(2))*.5;
direction = [0 0 1];
coordsys  = 'camera';
pos  = get(gca,'cameraposition' );
targ = get(gca,'cameratarget'   );
dar  = get(gca,'dataaspectratio');
up   = get(gca,'cameraupvector' );
[newPos newUp] = camrotate(pos,targ,dar,up,-dx,-dy,coordsys,direction);
set(gca,'cameraposition', newPos, 'cameraupvector', newUp);
vdata.old_pt = new_pt;
set(view3dObj,'UserData',vdata)

% ---------------------------------------------- view3dzFcn ----------
function view3dzFcn

view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
vdata = get(view3dObj,'UserData');
new_pt = get(0,'PointerLocation');
old_pt = vdata.old_pt;
dy = (new_pt(2) - old_pt(2))*.5;
camroll(gca,-dy)
vdata.old_pt = new_pt;
set(view3dObj,'UserData',vdata)

% ---------------------------------------------- view3dzoomFcn -------
function view3dzoomFcn
view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
vdata = get(view3dObj,'UserData');
new_pt = get(0,'PointerLocation');
old_pt = vdata.old_pt;
dy = (new_pt(2) - old_pt(2))/abs(old_pt(2));
camzoom(gca,1-dy)
vdata.old_pt = new_pt;
set(view3dObj,'UserData',vdata)

% ---------------------------------------------- view3dpanFcn --------
function view3dpanFcn

view3dObj  = findobj(allchild(gcf),'Tag','view3dObj');
vdata = get(view3dObj,'UserData');
new_pt = get(0,'PointerLocation');
old_pt = vdata.old_pt;
dx = (new_pt(1) - old_pt(1))/old_pt(1)*4;
dy = (new_pt(2) - old_pt(2))/old_pt(2)*4;
campan(gca,-dx,-dy,'camera')
vdata.old_pt = new_pt;
set(view3dObj,'UserData',vdata)

% ---------------------------------------------- make view3dObj ------
function view3dObj = makeview3dObj(fig)

% save the previous state of the figure window
vdata.uistate  = uisuspend(fig);
% the data structure
vdata.what     = [];
vdata.olp_pt   = [];
vdata.textbox  = [];
vdata.oldunits = [];
vdata.oldkeypressfcn = get(fig,'KeyPressFcn');
% view3dObj
view3dObj = uicontrol('style','text','parent',fig,'Units','Pixels',...
                      'Position',[2 2 130 20],'Visible','off', ...
                      'HandleVisibility','off','tag','view3dObj');
vdata.textbox  = view3dObj;
% store current view
ViewData = manageViewData('get_from_axes');
set(get(gca,'zlabel'),'UserData',ViewData);
% functions
set(fig,'WindowButtonDownFcn','view3d(''down'')');
set(fig,'WindowButtonUpFcn','view3d(''up'')');
set(fig,'WindowButtonMotionFcn','');
set(fig,'ButtonDownFcn','');
set(fig,'KeyPressFcn','view3d(''keypress'')');

set(view3dObj,'UserData',vdata);
% ---------------------------------------------- manage ViewData -----
function ViewData = manageViewData(how,data)

if nargin == 1 ; data = [];end
props = {
   'DataAspectRatio'
   'DataAspectRatioMode'
   'CameraPosition'
   'CameraPositionMode'
   'CameraTarget'
   'CameraTargetMode'
   'CameraUpVector'
   'CameraUpVectorMode'
   'CameraViewAngle'
   'CameraViewAngleMode'
   'PlotBoxAspectRatio'
   'PlotBoxAspectRatioMode'
   'Units'
   'Position'
   'View'
   'Projection'
};
if strcmp(how,'get_from_axes')
   ViewData = get(gca,props);
elseif strcmp(how,'get_stored')
   ViewData = get(get(gca,'zlabel'),'UserData');
elseif strcmp(how,'set_axes')
   set(gca,props,data)
   ViewData = [];
end
% -------------------------------------------------------------------------
% get some pointer shapes
function shape = pointershapes(arg)

if strcmp(arg,'zoom')
% -- zoom
shape=[ 2   2   2   2   2   2   2   2   2   2 NaN NaN NaN NaN NaN NaN  ;
        2   1   1   1   1   1   1   1   1   2 NaN NaN NaN NaN NaN NaN  ;
        2   1   2   2   2   2   2   2   2   2 NaN NaN NaN NaN NaN NaN  ;
        2   1   2   1   1   1   1   1   1   2 NaN NaN NaN NaN NaN NaN  ;
        2   1   2   1   1   1   1   1   2 NaN NaN NaN NaN NaN NaN NaN  ;
        2   1   2   1   1   1   1   2 NaN NaN NaN NaN NaN NaN NaN NaN  ;
        2   1   2   1   1   1   1   1   2 NaN NaN NaN   2   2   2   2  ;
        2   1   2   1   1   2   1   1   1   2 NaN   2   1   2   1   2  ;
        2   1   2   1   2 NaN   2   1   1   1   2   1   1   2   1   2  ;
        2   2   2   2 NaN NaN NaN   2   1   1   1   1   1   2   1   2  ;
      NaN NaN NaN NaN NaN NaN NaN NaN   2   1   1   1   1   2   1   2  ;
      NaN NaN NaN NaN NaN NaN NaN   2   1   1   1   1   1   2   1   2  ;
      NaN NaN NaN NaN NaN NaN   2   1   1   1   1   1   1   2   1   2  ;
      NaN NaN NaN NaN NaN NaN   2   2   2   2   2   2   2   2   1   2  ;
      NaN NaN NaN NaN NaN NaN   2   1   1   1   1   1   1   1   1   2  ;
      NaN NaN NaN NaN NaN NaN   2   2   2   2   2   2   2   2   2   2  ];
elseif strcmp(arg,'pan')
% -- pan
shape=[ NaN NaN NaN NaN NaN NaN NaN   2   2 NaN NaN NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN NaN   2   1   1   2 NaN NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN   2   1   1   1   1   2 NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN   1   1   1   1   1   1 NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN NaN   2   1   1   2 NaN NaN NaN NaN NaN NaN ;
        NaN NaN   2   1 NaN NaN   2   1   1   2 NaN NaN   1   2 NaN NaN ;
        NaN   2   1   1   2   2   2   1   1   2   2   2   1   1   2 NaN ;
          2   1   1   1   1   1   1   1   1   1   1   1   1   1   1   2 ;
          2   1   1   1   1   1   1   1   1   1   1   1   1   1   1   2 ;
        NaN   2   1   1   2   2   2   1   1   2   2   2   1   1   2 NaN ;
        NaN NaN   2   1 NaN NaN   2   1   1   2 NaN NaN   1   2 NaN NaN ;
        NaN NaN NaN NaN NaN NaN   2   1   1   2 NaN NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN   1   1   1   1   1   1 NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN   2   1   1   1   1   2 NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN NaN   2   1   1   2 NaN NaN NaN NaN NaN NaN ;
        NaN NaN NaN NaN NaN NaN NaN   2   2 NaN NaN NaN NaN NaN NaN NaN ];
elseif strcmp(arg,'rot')
% -- rot
shape=[ NaN NaN NaN   2   2   2   2   2 NaN   2   2 NaN NaN NaN NaN NaN ;
        NaN NaN NaN   1   1   1   1   1   2   1   1   2 NaN NaN NaN NaN ;
        NaN NaN NaN   2   1   1   1   1   2   1   1   1   2 NaN NaN NaN ;
        NaN NaN   2   1   1   1   1   1   2   2   1   1   1   2 NaN NaN ;
        NaN   2   1   1   1   2   1   1   2 NaN NaN   2   1   1   2 NaN ;
        NaN   2   1   1   2 NaN   2   1   2 NaN NaN   2   1   1   2 NaN ;
          2   1   1   2 NaN NaN NaN NaN NaN NaN NaN NaN   2   1   1   2 ;
          2   1   1   2 NaN NaN NaN NaN NaN NaN NaN NaN   2   1   1   2 ;
          2   1   1   2 NaN NaN NaN NaN NaN NaN NaN NaN   2   1   1   2 ;
          2   1   1   2 NaN NaN NaN NaN NaN NaN NaN NaN   2   1   1   2 ;
        NaN   2   1   1   2 NaN NaN   2   1   2 NaN   2   1   1   2 NaN ;
        NaN   2   1   1   2 NaN NaN   2   1   1   2   1   1   1   2 NaN ;
        NaN NaN   2   1   1   1   2   2   1   1   1   1   1   2 NaN NaN ;
        NaN NaN NaN   2   1   1   1   2   1   1   1   1   2 NaN NaN NaN ;
        NaN NaN NaN NaN   2   1   1   2   1   1   1   1   1 NaN NaN NaN ;
        NaN NaN NaN NaN NaN   2   2 NaN   2   2   2   2   2 NaN NaN NaN ];

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function gray = image_threshold3 ( gray, a, b )

%*****************************************************************************80
%
%% IMAGE_THRESHOLD3 resets gray pixels to white based on a threshold.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, uint8 GRAY(:,:), the original data.
%
%    Input, uint8 A, B, the threshold values.  Pixels of value less
%    than A or greater than B are reset to white.  0 <= A <= B <= 255
%
%    Output, uint8 GRAY(:,:), the thresholded image data.
%
  i = find ( gray < a );
  j = find ( b < gray );

  gray(i) = 255;
  gray(j) = 255;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function results = lowrank_sparse_fista( data, params, funs )
% FISTA implementation for solving
%   min 0.5*||AY - T|| + \lambda_*\sum_i||Xi,(i)||_* + \lambda1*||E||_1
% where Z = [ X1,...,X_N,E ]', Y is the FISTA auxiliary variable,
% A is a linear operator, T is given observation.

% tic;
% function handles
eval_f = funs.f;
eval_grad_f = funs.grad_f;
Aprod = funs.Aprod;
Atprod = funs.Atprod;
solve_prox_grad = funs.prox_grad_solver;

T = data.T;
if params.IsTC
    T = data.b;
end
N = length( size(params.X0) );
Y = cell(1,N+1);
for i = 1:N
    Y{i} = params.X0;
end
Y{N+1} = params.E0;
Z = Y;

use_cont = isfield( params, 'use_cont' ) && params.use_cont;
lambda1 = params.lambda;
lambdaS = params.lambdaS;
if use_cont
    if lambdaS > 0
        lambdaS_min = lambdaS;
        r = lambda1 / lambdaS;      lambdaS = 0.99*norm(T);     lambda1 = lambdaS*r;
    else
        r = params.rRatio /sqrt( max(size(data.T)) );
        lambdaS = 0.99*norm(T);     lambda1 = lambdaS*r;
        lambdaS_min = lambdaS * 1e-5;
    end
    eta = 0.97;
end

mu = params.mu0;
beta = 0.5;
t = 0;

fZ = eval_f( Z, Aprod, T );
trnorm = 0;
oneNorm = 0;

for iter = 1:params.max_iter

    % compute gradient at Y
    grad_f = eval_grad_f( Y, Aprod, Atprod, T );

    Zp = Z;     trnorm_p = trnorm;    fZp = fZ;     oneNorm_p = oneNorm;
    [ Z, trnorm ] = solve_prox_grad( Y, mu, grad_f, lambdaS, lambda1 );
    % backtrack linesearch
    fY = eval_f( Y, Aprod, T );
    ZYdiff = tensor_array_diff( Z, Y );
    approx = fY + tensor_array_innerprod( grad_f, ZYdiff ) + tensor_array_norm(ZYdiff)^2 / (2*mu);
    fZ = eval_f( Z, Aprod, T );
    while fZ > approx
        mu = mu * beta;
        [ Z, trnorm ] = solve_prox_grad( Y, mu, grad_f, lambdaS, lambda1 );
        ZYdiff = tensor_array_diff( Z, Y );
        approx = fY + tensor_array_innerprod( grad_f, ZYdiff ) + tensor_array_norm(ZYdiff)^2 / (2*mu);
        fZ = eval_f( Z, Aprod, T );
    end
    oneNorm = tensor_1norm(Z{N+1});

    % FISTA acceleration step
    tp = t;
    t = (1 + sqrt(1+4*tp^2))/2;
    Zdiff = tensor_array_diff( Z, Zp );
    Y = tensor_array_add( Zp, tensor_array_scale(Zdiff,(tp-1)/t) );

    % compute optimality stats
    Fp = fZp + lambdaS*trnorm_p + lambda1*oneNorm_p;
    F = fZ + lambdaS*trnorm + lambda1*oneNorm;
    rel_F = abs(F-Fp)/Fp;
    denom = tensor_array_norm(Zp);
    rel_X = tensor_array_norm(Zdiff) / denom;   if denom == 0; rel_X = 1; end
    rel_err = norm(ten_sum_all(Z(1:N))-data.X) / norm(data.X);

    % print
    if params.verbose && rem( iter, 50 ) == 0
        fprintf('Iter: %d,   fit: %3.2e,   rel_F: %3.2e,    rel_X: %3.2e,   rel_err: %3.2e, S: %3.2e\n', ...
            iter, fZ, rel_F, rel_X, rel_err, lambdaS);
    end

    if max(rel_F,rel_X) < params.opt_tol
        break;
    end

    % update lambda's
    if use_cont
        lambdaS = max( lambdaS*eta, lambdaS_min );
        lambda1 = lambdaS*r;
    end
end

results.vars = Z;
results.V = params.V0;
results.T = data.T;
results.iter = iter;
% results.cpu = toc;
results.mu = mu;
results.lambda1 = lambda1;
results.lambdaS = lambdaS;
results.IsTC = params.IsTC;

end

author: andrewssobral

repo: lrslibrary


text:
% LINSOLVE_LDLCHOL - Solves a matrix-vector equation when the matrix is a
%                    sparse symmetric positive definite matrix.
%
% Solves K*Y = X. Usage:
%
%   Y = LINSOLVE_LDLCHOL(LD,X)
%
% where LD = LDLCHOL(K).
%
% The function is just a simple wrapper in accordance with the naming
% conventions.
%
% See also LDLCHOL, LDLSOLVE.

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

function x = linsolve_ldlchol(LD,x,q)

%fprintf('Sparsity in linsolve_ldlchol: %f\n', sparsity(LD));

if nargin < 3
  x = ldlsolve(LD,x);
else
  x(q(:),:) = ldlsolve(LD,x(q(:),:));
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [ml] = l2ml(l)
% Convert volume from liters to milliliters.
% Chad Greene 2012
ml = l*1000;

author: Sable

repo: mcbench-benchmarks


text:
function options = cplexoptimset (varargin)
%%
% Purpose
% Create or edit optimization options structure.
%
% Syntax
%    options = cplexoptimset('param1',value1,'param2',value2,...)
%    cplexoptimset
%    options = cplexoptimset
%    options = cplexoptimset('cplex')
%    options = cplexoptimset(oldopts,'param1',value1,...)
%    options = cplexoptimset(oldopts,newopts)
%
% Description
% options = cplexoptimset('param1',value1,'param2',value2,...) creates an
% optimization options structure called options, in which the specified
% options (param) have specified values. The specified options correspond
% to the options in the MATLAB(R) Optimization Toolbox.Any unspecified
% options are set to [] (options with value [] indicate that the default
% value for that option should be used when options is passed to the
% optimization function).
%
% cplexoptimset with no input or output arguments displays a complete list
% of options with their valid values.
%
% options = cplexoptimset (with no input arguments) creates an option
% structure options where all fields are set to [], which instructs
% CPLEX(R) to use all default parameter values.
%
% options = cplexoptimset ('cplex') creates a structure options, which
% contains all of the CPLEX parameters.  Use this method when you need to
% set parameters that do not correspond to the options in the MATLAB
% Optimization Toolbox.
%
% options = cplexoptimset(oldopts,'param1',value1,...) creates a copy of
% oldopts, modifying the specified options with the specified values.
%
% options = cplexoptimset(oldopts,newopts) combines an existing options
% structure, oldopts, with a new options structure, newopts. Any options in
% newopts with nonempty values overwrite the corresponding old options in
% oldopts.
%
% Options corresponding to the MATLAB Optimization Toolbox
% Diagnostics          'on' | {'off'}
% Display              'off' | 'iter' | 'final' | 'notify'
% MaxIter               refer to cplex.Param.simplex.limits.iterations
% Simplex               refer to cplex.Param.lpmethod
%                       refer to cplex.Param.qpmethod
%                       refer to cplex.Param.mip.strategy.startalgorithm
% BranchStrategy        refer to cplex.Param.mip.strategy.variableselect
% MaxNodes              refer to cplex.Param.mip.limits.nodes
% MaxTime               refer to cplex.Param.timelimit
% NodeDisplayInterval   refer to cplex.Param.mip.interval
% NodeSearchStrategy    refer to cplex.Param.mip.strategy.nodeselect
% TolFun                refer to cplex.Param.simplex.tolerances.optimality
% TolXInteger           refer to cplex.Param.mip.tolerances.integrality
% TolRLPFun             refer to cplex.Param.simplex.tolerances.optimality
%
% For other CPLEX options, refer to the CPLEX Parameters Reference
% Manual.
%
% Examples
% To turn on the optimizer output and set a node limit of 400:
% options = cplexoptimset('Diagnostics', 'on', 'MaxNodes', 400);
%
% To set a node limit of 400 and instruct CPLEX to use traditional branch
% and cut style search:
% opt = cplexoptimset('cplex');
% opt.mip.limits.nodes=400;
% opt.mip.strategy.search=1;
%
%  See also cplexoptimget
%

% ---------------------------------------------------------------------------
% File: cplexoptimset.m
% ---------------------------------------------------------------------------
% Licensed Materials - Property of IBM
% 5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55
% Copyright IBM Corporation 2008, 2010. All Rights Reserved.
%
% US Government Users Restricted Rights - Use, duplication or
% disclosure restricted by GSA ADP Schedule Contract with
% IBM Corp.
% ---------------------------------------------------------------------------

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [params, names] = mixturePriorExtractParam(prior)

% MIXTUREPRIOREXTRACTPARAM Extract params from mixture prior structure.

% COPYRIGHT : Antti Honkela, 2013

% SHEFFIELDML

params = zeros(1, prior.nParams);
if nargout > 1
  names = cell(1, prior.nParams);
end

begindex = 1;
for k=1:length(prior.comp),
  myindex = begindex:(begindex+prior.comp{k}.nParams-1);
  if nargout > 1,
    [params(myindex), names(myindex)] = priorExtractParam(prior.comp{k});
  else
    params(myindex) = priorExtractParam(prior.comp{k});
  end
  begindex = begindex + prior.comp{k}.nParams;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [datapt, unitStr] = lvd_ExtremaTasks(stateLogEntry, subTask, extremum)
%lvd_ExtremaTasks Summary of this function goes here
%   Detailed explanation goes here

    switch subTask
        case 'extremumValue'
            extremaStates = stateLogEntry.getAllExtremaStates();
            extremaState = extremaStates([extremaStates.extrema] == extremum);

            if(not(isempty(extremaState)))
                extremaState = extremaState(1);
                datapt = extremaState.value;

                if(isempty(datapt))
                    datapt = -1;
                end

                unitStr = extremaState.extrema.unitStr;
            else
                datapt = -1;
            end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function c = assemble_mass ( node_num, node_x, element_num, element_node, ...
  quad_num )

%*****************************************************************************80
%
%% ASSEMBLE_MASS assembles the finite element mass matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_X(NODE_NUM), the coordinates of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(2,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer QUAD_NUM, the number of quadrature points used in assembly.
%
%    Output, sparse real C(NODE_NUM,NODE_NUM), the finite element mass matrix.
%
%  Local parameters:
%
%    Local, real BI, DBIDX, the value of some basis function
%    and its first derivative at a quadrature point.
%
%    Local, real BJ, DBJDX, the value of another basis
%    function and its first derivative at a quadrature point.
%

%
%  Initialize the arrays.
%
  c = [];
  c = sparse ( node_num, node_num );
%
%  Get the quadrature weights and nodes.
%
  [ reference_w, reference_q ] = quadrature_set ( quad_num );
%
%  Consider each ELEMENT.
%
  for element = 1 : element_num

    element_x(1:2,1) = node_x(element_node(1:2,element));

    element_q(1:quad_num,1) = reference_to_physical ( element, element_node, ...
      node_x, quad_num, reference_q );

    element_area = element_x(2) - element_x(1);

    element_w(1:quad_num,1) = ( element_area / 2.0 ) * reference_w(1:quad_num);
%
%  Consider the QUAD-th quadrature point in the element.
%
    for quad = 1 : quad_num
%
%  Consider the TEST-th test function.
%
%  We generate an integral for every node associated with an unknown.
%
      for i = 1 : 2

        test = element_node(i,element);

        [ bi, dbidx ] = basis_function ( test, element, node_x, 1, element_q(quad) );
%
%  Consider the BASIS-th basis function, which is used to form the
%  value of the solution function.
%
        for j = 1 : 2

          basis = element_node(j,element);

          [ bj, dbjdx ] = basis_function ( basis, element, node_x, 1, element_q(quad) );

          c(test,basis) = c(test,basis) + element_w(quad) * bi * bj;

        end

      end

    end

  end

  return
end


text:
function extractLabels(jabfile,outfile)

J = load(jabfile,'-mat');
oo = fopen(outfile,'w');

for exp = 1:numel(J.x.labels)
  for ndx = 1:numel(J.x.labels(exp).t0s)
    for bndx = 1:numel(J.x.labels(exp).t0s{ndx})
      fprintf(oo,'%s %d %d %d %s\n',J.x.expDirNames{exp},...
        J.x.labels(exp).flies(ndx),...
        J.x.labels(exp).t0s{ndx}(bndx),...
        J.x.labels(exp).t1s{ndx}(bndx)-1,...
        J.x.labels(exp).names{ndx}{bndx}...
        );
    end

  end

end
fclose(oo);

author: kristinbranson

repo: JAABA


text:
function [y] = spm_phase_shuffle(x,n)
% phase-shuffling of a vector
% FORMAT [y] = spm_phase_shuffle(x,[n])
% x   - data matrix (time-series in columns)
% n   - optional window length for windowed shuffling
%__________________________________________________________________________
% Copyright (C) 2007-2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_phase_shuffle.m 6654 2015-12-22 12:55:36Z spm $

try

    % randomise phase - WFT
    %----------------------------------------------------------------------
    k     = 1:fix(n/2);
    for i = 1:size(x,2);
        C      = spm_wft(x(:,i),k,n);
        W      = abs(C).*exp(1i*angle(C(randperm(size(C,1)),:)));
        y(:,i) = spm_iwft(W,k,n)';
    end

catch

    % randomise phase - FT
    %----------------------------------------------------------------------
    n              = size(x,1);
    s              = fft(x);
    i              = 2:ceil(n/2);
    r              = rand(length(i),size(x,2))*2*pi - pi;
    p              = zeros(n,size(x,2));
    p(i,:)         =  r;
    p(n - i + 2,:) = -r;
    s              = abs(s).*exp(1i*p);
    y              = real(ifft(s));

end


text:
function [SDR,ISR,SIR,SAR]=bss_mcrit(s_true,e_spat,e_interf,e_artif)

% BSS_MCRIT Computation of some separation performance criteria given the
% decomposition of an estimated source image into four components
% representing respectively the true source image, spatial (or filtering)
% distortion, interference and artifacts.
%
% [SDR,ISR,SIR,SAR]=bss_mcrit(s_true,e_spat,e_interf,e_artif)
%
% Inputs:
% s_true: I x T matrix containing the true source image (one row per channel)
% e_spat: I x T matrix containing the spatial (or filtering) distortion component
% e_interf: I x T matrix containing the interference component
% e_artif: I x T matrix containing the artifacts component
%
% Outputs:
% SDR: Signal to Distortion Ratio
% ISR: source Image to Spatial distortion Ratio
% SIR: Source to Interference Ratio
% SAR: Sources to Artifacts Ratio

%%% Errors %%%
if nargin<4, error('Not enough input arguments.'); end
[It,Tt]=size(s_true);
[Is,Ts]=size(e_spat);
[Ii,Ti]=size(e_interf);
[Ia,Ta]=size(e_artif);
if ~((It==Is)&&(It==Ii)&&(It==Ia)), error('All the components must have the same number of channels.'); end
if ~((Tt==Ts)&&(Tt==Ti)&&(Tt==Ta)), error('All the components must have the same duration.'); end

%%% Energy ratios %%%
% SDR
SDR=10*log10(sum(sum(s_true.^2))/sum(sum((e_spat+e_interf+e_artif).^2)));
% ISR
ISR=10*log10(sum(sum(s_true.^2))/sum(sum(e_spat.^2)));
% SIR
SIR=10*log10(sum(sum((s_true+e_spat).^2))/sum(sum(e_interf.^2)));
% SAR
SAR=10*log10(sum(sum((s_true+e_spat+e_interf).^2))/sum(sum(e_artif.^2)));

return;

repo: Voice-Separation-and-Enhancement


text:
function spm_mask(P1, P2, thresh)
% Mask images
% FORMAT spm_mask(P1, P2, thresh)
% P1     - matrix of input image filenames from which
%          to compute the mask.
% P2     - matrix of input image filenames on which
%          to apply the mask.
% thresh - optional threshold(s) for defining the mask.
% The masked images are prepended with the prefix `m'.
%
% If any voxel in the series of images is zero (for data types without
% a floating point representation) or does not have a finite value (for
% floating point and double precision images), then that voxel is set to
% NaN or zero in all the images.  If a threshold, or vector of
% thresholds is passed, then the masking is based on voxels whos
% values are above all the thresholds.
%
% Images sampled in different orientations and positions can be passed
% to the routine.
%__________________________________________________________________________
% Copyright (C) 1999-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_mask.m 4419 2011-08-03 18:42:35Z guillaume $

persistent runonce
if isempty(runonce)
    warning('ImCalc should be preferred to spm_mask whenever possible.');
    runonce = 1;
end

SVNid = '$Rev: 4419 $';

%-Say hello
%--------------------------------------------------------------------------
SPMid = spm('FnBanner',mfilename,SVNid);

%-Parameters & Arguments
%--------------------------------------------------------------------------
if ~nargin
    [P1, sts] = spm_select(Inf,'image','Images to compute mask from');
    if ~sts, return; end
    [P2, sts] = spm_select(Inf,'image','Images to apply mask to');
    if ~sts, return; end
end

if nargin==1
    P2 = P1;
end

V1 = spm_vol(P1);
V2 = spm_vol(P2);

if nargin==3
    if numel(thresh)==1
        thresh = repmat(thresh,numel(V1),1);
    else
        if numel(V1) ~= numel(thresh)
            error('Input argument ''thresh'' has wrong size.');
        end
        thresh = thresh(:);
    end
end

m1 = numel(V1);
m2 = numel(V2);

%-Create headers
%--------------------------------------------------------------------------
VO = V2;
for i=1:m2
    [pth,nm,ext,num] = spm_fileparts(deblank(VO(i).fname));
    VO(i).fname      = fullfile(pth,['m', nm, ext, num]);
    VO(i).descrip    = 'Masked';
    VO(i).mat        = VO(1).mat;
    VO(i).dim(1:3)   = VO(1).dim(1:3);
end
VO  = spm_create_vol(VO);
M   = VO(1).mat;
dim = VO(1).dim(1:3);

%-Compute masked images
%--------------------------------------------------------------------------
spm_progress_bar('Init',VO(1).dim(3),'Masking','planes completed');

for j=1:dim(3)

    msk = true(dim(1:2));
    Mi  = spm_matrix([0 0 j]);

    % Load slice j from all images
    for i=1:m1
        M1  = M\V1(i).mat\Mi;
        %if sum((M1(:)-Mi(:)).^2<eps) M1 = Mi; end;

        img = spm_slice_vol(V1(i),M1,dim(1:2),[0 NaN]);

        msk = msk & isfinite(img);

        if ~spm_type(V1(i).dt(1),'nanrep')
            msk = msk & (img ~= 0);
        end

        if nargin == 3
            msk = msk & (img >= thresh(i));
        end
    end

    % Write the images
    for i=1:m2
        M1        = M\V2(i).mat\Mi;
        img       = spm_slice_vol(V2(i),M1,dim(1:2),[1 0]);
        img(~msk) = NaN;
        VO(i)     = spm_write_plane(VO(i),img,j);
    end

    spm_progress_bar('Set',j);
end

spm_progress_bar('Clear');


text:
function xyzankurAnimCompare(X, X2, fps)

% XYZANKURANIMCOMPARE Animate a prediction and ground truth for stick man from Agarwal & Triggs dataset.
% FORMAT
% DESC animates a matrix of x,y,z point clound positions representing the
% motion of the figure used to generate the silhouttes for Agarwal &
% Triggs silhouette data.
% ARG y : the data to animate.
% ARG y2 : other data to animate.
% ARG fps : the number of frames per second to animate (defaults to 24).
%
% SEEALSO : xyzankurVisualise, xyzankurModify
%
% COPYRIGHT : Carl Henrik Ek and Neil D. Lawrence, 2008

% MOCAP

  if(nargin<3)
    fps = 24;
    if(nargin<2)
      fid = 1;
      if(nargin<1)
        error('Too few arguments');
      end
    end
  end
  clf
  for(i = 1:1:size(X,1))
    if(i==1)
      subplot(1, 2, 1)
      handle = xyzankurVisualise(X(i,:), 1);
      subplot(1, 2, 2)
      handle2 = xyzankurVisualise(X2(i,:), 1);
    else
      xyzankurModify(handle, X(i,:));
      xyzankurModify(handle2, X2(i,:));
    end
    pause(1/fps);
  end
end


function [x, y] = trans_cam2fisheye(xcam, ycam, zcam, M, D)
fx = M(1,1); fy = M(2,2);
cx = M(1,3); cy = M(2,3);
if D ~= 0
    k1 = D(1); k2 = D(2); k3 = D(3); k4 = D(4);
end

rcam = sqrt(xcam .^ 2 + ycam .^ 2 + zcam .^ 2) ;
alpha = acos(zcam ./ rcam) ;

if D ~= 0
    alpha_2 = alpha .^ 2;
    k_radial =  1 + k1 * alpha_2 + k2 * alpha_2.^2 + k3 * alpha_2.^3 + k4 * alpha_2.^4;
    alpha_d = alpha .* k_radial;
else
    alpha_d = alpha;
end

rcam_xy = sqrt(xcam.^2 + ycam.^2);
alpha_d_x = xcam ./ rcam_xy .* alpha_d;
alpha_d_y = ycam ./ rcam_xy .* alpha_d;

x = alpha_d_x .* fx + cx ;
y = alpha_d_y .* fy + cy ;

end


text:
%OUTEROP calculate an outer operation on two vectors
%   function y=outerop(A,B,OPERATOR)
%
%   Calculates resultant matrix when the OPERATOR is applied
%   to all combinations of the elements of vector A and the
%   elements of vector B e.g. the outer product of A and B
%   is outerop(A,B,'*'), the outer sum of A and B
%   is outerop(A,B,'+')
%
%   If OPERATOR is omitted '+' is assumed
%
%   This function is equivalent to the
%   APL language's circle.dot operation. e.g.
%   in APL Ao.*B is the outer product
%   of A and B % and Ao.+B is the outer sum.
%   Ideally it would work on matrices but in practice
%   I've only ever needed to use it with vectors.
%
%   Copyright Murphy O'Brien 2005
%   all rights unreserved
%
function y=outerop(a,b,operator)

if nargin<3
    operator='+';                       % for only two arguments assume outerproduct
end

if isequal(operator,'*')                % common operator
    y=a(:)*b(:)';
else
  outera=a(:)*ones(1,length(b));        % these are the matrices that
  outerb=ones(length(a),1)*b(:).';      % meshgrid(A,B) would generate
  functionHandle=str2func(operator);  % new R14 functionality
  y=functionHandle(outera,outerb);    % allows faster/neater method
end


text:
function f = p02_f ( m, n, x )

%*****************************************************************************80
%
%% P02_F evaluates the objective function for problem 02.
%
%  Discussion:
%
%    This function is also known as the weighted sphere model.
%
%    The function is continuous, convex, and unimodal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marcin Molga, Czeslaw Smutnicki,
%    Test functions for optimization needs.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of arguments.
%
%    Input, real X(M,N), the arguments.
%
%    Output, real F(N), the function evaluated at the arguments.
%
  f = zeros ( n, 1 );

  y = r8vec_indicator ( m );

  for j = 1 : n
    f(j) = sum ( y(1:m) .* x(1:m,j).^2 );
  end

  return
end


text:
%% Angular correlation between spherical harmonics profiles
%% Reference:
%%     Adam W. Anderson.
%%     Measurement of Fiber Orientation Distributions Using High Angular Resolution Diffusion Imaging.
%%     Magn. Reson. Med., 54(5):1194-1206, 2005.  (Eq.(71))
%% Example:
%%     [corr] = angular_correlation(coeff_u,coeff_v);
function [corr] = angular_correlation(coeff_u,coeff_v)
%%=============================================================
%% Project:   Spherical Harmonics
%% Module:    $RCSfile: angular_correlation.m,v $
%% Language:  MATLAB
%% Author:    $Author: bjian $
%% Date:      $Date: 2007/12/27 06:23:35 $
%% Version:   $Revision: 1.3 $
%%=============================================================

corr = (coeff_u'*coeff_v)/(norm(coeff_u)*norm(coeff_v));


text:
function [min] = s2min(s)
% Convert time from seconds to minutes.
% Chad A. Greene 2012
min = s/60;


text:
function display(item)

% function display(item)
% Display a configuration object
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: display.m 4166 2011-01-17 15:06:41Z volkmar $

rev = '$Rev: 4166 $'; %#ok

disp(' ');
disp([inputname(1),' = ']);
disp(' ');
disp(item);
disp(' ');


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z) .* (1 .- sigmoid(z));

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

end

author: ecmadao

repo: Coding-Guide


text:
function i4row_sum_test ( )

%*****************************************************************************80
%
%% I4ROW_SUM_TEST tests I4ROW_SUM;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4ROW_SUM_TEST\n' );
  fprintf ( 1, '  I4ROW_SUM computes row sums;\n' );

  k = 0;
  for i = 1 : m
    for j = 1 : n
      k = k + 1;
      a(i,j) = k;
    end
  end

  i4mat_print ( m, n, a, '  The matrix:' );

  rowsum(1:m) = i4row_sum ( m, n, a );

  mean(1:m) = i4row_mean ( m, n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sum, mean:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    fprintf ( 1, '  %6d  %6d  %9f\n', i, rowsum(i), mean(i) );
  end

  return
end


text:
function [estimate] = ft_inverse_beamformer_dics(leadfield, Cf, varargin)

% FT_INVERSE_BEAMFORMER_DICS estimates the source power or source
% coherence according to the Dynamic Imaging of Coherent Sources
% method.
%
% Use as
%   estimate = ft_inverse_beamformer_dics(leadfield, Cf, ...)
% where
%   leadfield  = leadfield of the source of interest or a cell-array with leadfields for multiple sources
%   Cf         = cross-spectral density matrix of the data
% and
%   estimate   = structure with the estimated source parameters
%
% Additional options should be specified in key-value pairs and can be
%  'Pr'               = power of the external reference channel
%  'Cr'               = cross spectral density between all data channels and the external reference channel
%  'refdip'           = location of dipole with which coherence is computed
%  'lambda'           = regularisation parameter
%  'powmethod'        = can be 'trace' or 'lambda1'
%  'feedback'         = give ft_progress indication, can be 'text', 'gui' or 'none'
%  'fixedori'         = use fixed or free orientation,                 can be 'yes' or 'no'
%  'projectnoise'     = project noise estimate through filter,         can be 'yes' or 'no'
%  'realfilter'       = construct a real-valued filter,                can be 'yes' or 'no'
%  'keepfilter'       = remember the beamformer filter,                can be 'yes' or 'no'
%  'keepleadfield'    = remember the forward computation,              can be 'yes' or 'no'
%  'keepcsd'          = remember the estimated cross-spectral density, can be 'yes' or 'no'
%
% This implements Joachim Gross et al. 2001

% Copyright (C) 2003-2010, Robert Oostenveld

% these optional settings do not have defaults
Pr             = keyval('Pr',            varargin);
Cr             = keyval('Cr',            varargin);
refdip         = keyval('refdip',        varargin);
powmethod      = keyval('powmethod',     varargin); % the default for this is set below
realfilter     = keyval('realfilter',    varargin); % the default for this is set below
% these optional settings have defaults
feedback       = keyval('feedback',      varargin); if isempty(feedback),      feedback = 'text';            end
keepcsd        = keyval('keepcsd',       varargin); if isempty(keepcsd),       keepcsd = 'no';               end
keepfilter     = keyval('keepfilter',    varargin); if isempty(keepfilter),    keepfilter = 'no';            end
keepleadfield  = keyval('keepleadfield', varargin); if isempty(keepleadfield), keepleadfield = 'no';         end
lambda         = keyval('lambda',        varargin); if isempty(lambda  ),      lambda = 0;                   end
projectnoise   = keyval('projectnoise',  varargin); if isempty(projectnoise),  projectnoise = 'yes';         end
fixedori       = keyval('fixedori',      varargin); if isempty(fixedori),      fixedori = 'no';              end

% convert the yes/no arguments to the corresponding logical values
% FIXME use istrue
keepcsd        = strcmp(keepcsd,       'yes');
keepfilter     = strcmp(keepfilter,    'yes');
keepleadfield  = strcmp(keepleadfield, 'yes');
projectnoise   = strcmp(projectnoise,  'yes');
fixedori       = strcmp(fixedori,      'yes');

% FIXME besides regular/complex lambda1, also implement a real version

% default is to use the largest singular value of the csd matrix, see Gross 2001
if isempty(powmethod)
  powmethod = 'lambda1';
end

% default is to be consistent with the original description of DICS in Gross 2001
if isempty(realfilter)
  realfilter = 'no';
end

% use these two logical flags instead of doing the string comparisons each time again
powtrace   = strcmp(powmethod, 'trace');
powlambda1 = strcmp(powmethod, 'lambda1');

% dics has the following sub-methods, which depend on the additional input arguments
if ~isempty(Cr) && ~isempty(Pr) && isempty(refdip)
  % compute cortico-muscular coherence, using reference cross spectral density
  submethod = 'dics_refchan';
elseif isempty(Cr) && isempty(Pr) && ~isempty(refdip)
  % compute cortico-cortical coherence with a dipole at the reference position
  submethod = 'dics_refdip';
elseif isempty(Cr) && isempty(Pr) && isempty(refdip)
  % only compute power of a dipole at the grid positions
  submethod = 'dics_power';
else
  error('invalid combination of input arguments for dics');
end

if ~iscell(leadfield)
  % the leadfield specifies a single source
  leadfield = {leadfield};
end
ndipoles = length(leadfield);

if ~isempty(Cr)
  % ensure that the cross-spectral density with the reference signal is a column matrix
  Cr = Cr(:);
end

isrankdeficient = (rank(Cf)<size(Cf,1));
if isrankdeficient
  warning('cross-spectral density matrix is rank deficient')
end

% it is difficult to give a quantitative estimate of lambda, therefore also
% support relative (percentage) measure that can be specified as string (e.g. '10%')
if ~isempty(lambda) && ischar(lambda) && lambda(end)=='%'
  ratio = sscanf(lambda, '%f%%');
  ratio = ratio/100;
  lambda = ratio * trace(Cf)/size(Cf,1);
end

if projectnoise
  % estimate the noise power, which is further assumed to be equal and uncorrelated over channels
  if isrankdeficient
    % estimated noise floor is equal to or higher than lambda
    noise = lambda;
  else
    % estimate the noise level in the covariance matrix by the smallest singular value
    noise = svd(Cf);
    noise = noise(end);
    % estimated noise floor is equal to or higher than lambda
    noise = max(noise, lambda);
  end
end

% the inverse only has to be computed once for all dipoles
if strcmp(realfilter, 'yes')
  % the filter is computed using only the leadfield and the inverse covariance or CSD matrix
  % therefore using the real-valued part of the CSD matrix here ensures a real-valued filter
  invCf = pinv(real(Cf) + lambda * eye(size(Cf)));
else
  invCf = pinv(Cf + lambda * eye(size(Cf)));
end

% start the scanning with the proper metric
ft_progress('init', feedback, 'scanning grid');
switch submethod

  case 'dics_power'
    % only compute power of a dipole at the grid positions
    for i=1:ndipoles
      lf = leadfield{i};

      if fixedori
        % compute the leadfield for the optimal dipole orientation
        % subsequently the leadfield for only that dipole orientation will be used for the final filter computation
        filt = pinv(lf' * invCf * lf) * lf' * invCf;
        [u, s, v] = svd(real(filt * Cf * ctranspose(filt)));
        maxpowori = u(:,1);
        eta = s(1,1)./s(2,2);
        lf  = lf * maxpowori;
        estimate.ori{i} = maxpowori;
        estimate.eta{i} = eta;
      end

      % construct the spatial filter
      filt = pinv(lf' * invCf * lf) * lf' * invCf;                % Gross eqn. 3, use PINV/SVD to cover rank deficient leadfield
      csd = filt * Cf * ctranspose(filt);                         % Gross eqn. 4 and 5

      % assign the output values
      if powlambda1
        estimate.pow(i) = lambda1(csd);                             % compute the power at the dipole location, Gross eqn. 8
      elseif powtrace
        estimate.pow(i) = real(trace(csd));                         % compute the power at the dipole location
      end
      if keepcsd
        estimate.csd{i} = csd;
      end
      if projectnoise
        if powlambda1
          estimate.noise(i) = noise * lambda1(filt * ctranspose(filt));
        elseif powtrace
          estimate.noise(i) = noise * real(trace(filt * ctranspose(filt)));
        end
        if keepcsd
          estimate.noisecsd{i} = noise * filt * ctranspose(filt);
        end
      end
      if keepfilter
        estimate.filter{i} = filt;
      end
      if keepleadfield
        estimate.leadfield{i} = lf;
      end
      ft_progress(i/ndipoles, 'scanning grid %d/%d\n', i, ndipoles);
    end

  case 'dics_refchan'
    % compute cortico-muscular coherence, using reference cross spectral density
    for i=1:ndipoles

      % get the leadfield for this source
      lf = leadfield{i};

      if fixedori
        % compute the leadfield for the optimal dipole orientation
        % subsequently the leadfield for only that dipole orientation will be used for the final filter computation
        filt = pinv(lf' * invCf * lf) * lf' * invCf;
        [u, s, v] = svd(real(filt * Cf * ctranspose(filt)));
        maxpowori = u(:,1);
        lf  = lf * maxpowori;
        estimate.ori{i} = maxpowori;
      end

      % construct the spatial filter
      filt = pinv(lf' * invCf * lf) * lf' * invCf;                     % use PINV/SVD to cover rank deficient leadfield

      if powlambda1
        [pow, ori] = lambda1(filt * Cf * ctranspose(filt));            % compute the power and orientation at the dipole location, Gross eqn. 4, 5 and 8
      elseif powtrace
        pow = real(trace(filt * Cf * ctranspose(filt)));               % compute the power at the dipole location
      end
      csd = filt*Cr;                                                   % Gross eqn. 6
      if powlambda1
        % FIXME this should use the dipole orientation with maximum power
        coh = lambda1(csd)^2 / (pow * Pr);                             % Gross eqn. 9
      elseif powtrace
        coh = norm(csd)^2 / (pow * Pr);
      end

      estimate.pow(i) = pow;
      estimate.coh(i) = coh;
      if keepcsd
        estimate.csd{i} = csd;
      end
      if projectnoise
        if powlambda1
          estimate.noise(i) = noise * lambda1(filt * ctranspose(filt));
        elseif powtrace
          estimate.noise(i) = noise * real(trace(filt * ctranspose(filt)));
        end
        if keepcsd
          estimate.noisecsd{i} = noise * filt * ctranspose(filt);
        end
      end
      if keepfilter
        estimate.filter{i} = filt;
      end
      ft_progress(i/ndipoles, 'scanning grid %d/%d\n', i, ndipoles);
    end

  case 'dics_refdip'

    if fixedori
      error('fixed orientations are not supported for beaming cortico-cortical coherence');
    end

    % get the leadfield of the reference source
    lf1 = refdip;

    % construct the spatial filter for the first (reference) dipole location
    filt1 = pinv(lf1' * invCf * lf1) * lf1' * invCf;       % use PINV/SVD to cover rank deficient leadfield
    if powlambda1
      Pref = lambda1(filt1 * Cf * ctranspose(filt1));      % compute the power at the first dipole location, Gross eqn. 8
    elseif powtrace
      Pref = real(trace(filt1 * Cf * ctranspose(filt1)));  % compute the power at the first dipole location
    end

    for i=1:ndipoles

      % get the leadfield for the second source, i.e. the one that is being scanned
      lf2 = leadfield{i};

      % construct the spatial filter for the second source
      filt2 = pinv(lf2' * invCf * lf2) * lf2' * invCf;     % use PINV/SVD to cover rank deficient leadfield
      csd = filt1 * Cf * ctranspose(filt2);                % compute the cross spectral density between the two dipoles, Gross eqn. 4

      if powlambda1
        pow = lambda1(filt2 * Cf * ctranspose(filt2));     % compute the power at the second dipole location, Gross eqn. 8
      elseif powtrace
        pow = real(trace(filt2 * Cf * ctranspose(filt2))); % compute the power at the second dipole location
      end
      if powlambda1
        coh = lambda1(csd)^2 / (pow * Pref);               % compute the coherence between the first and second dipole
      elseif powtrace
        coh = real(trace((csd)))^2 / (pow * Pref);         % compute the coherence between the first and second dipole
      end

      estimate.pow(i) = pow;
      estimate.coh(i) = coh;
      if keepcsd
        estimate.csd{i} = csd;
      end
      if projectnoise
        if powlambda1
          estimate.noise(i) = noise * lambda1(filt2 * ctranspose(filt2));
        elseif powtrace
          estimate.noise(i) = noise * real(trace(filt2 * ctranspose(filt2)));
        end
        if keepcsd
          estimate.noisecsd{i} = noise * filt2 * ctranspose(filt2);
        end
      end
      ft_progress(i/ndipoles, 'scanning grid %d/%d\n', i, ndipoles);
    end

end % switch submethod

ft_progress('close');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to obtain the largest singular value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [s, ori] = lambda1(x)
% determine the largest singular value, which corresponds to the power along the dominant direction
[u, s, v] = svd(x);
s   = s(1);
ori = u(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function to compute the pseudo inverse. This is the same as the
% standard Matlab function, except that the default tolerance is twice as
% high.
%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 2710 $  $Date: 2009/06/17 13:40:37 $
%   default tolerance increased by factor 2 (Robert Oostenveld, 7 Feb 2004)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = pinv(A,varargin)
[m,n] = size(A);
if n > m
  X = pinv(A',varargin{:})';
else
  [U,S,V] = svd(A,0);
  if m > 1, s = diag(S);
  elseif m == 1, s = S(1);
  else s = 0;
  end
  if nargin == 2
    tol = varargin{1};
  else
    tol = 10 * max(m,n) * max(s) * eps;
  end
  r = sum(s > tol);
  if (r == 0)
    X = zeros(size(A'),class(A));
  else
    s = diag(ones(r,1)./s(1:r));
    X = V(:,1:r)*s*U(:,1:r)';
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function [img1_rot, rot, err] = AlignLLPanoramas(img1, img2, bVisualization)
%
%
%     [img1_rot, rot, err] = AlignLLPanoramas(img1, img2, bVisualization)
%
%     This function finds the rotation around Y-axis in pixel for aligning
%     the panorma img1 (in longitude-latitude format) to the panorma
%     img2 (in longitude-latitude format).
%
%     Note that img1 and img2 have to share the same Y-axis in order to
%     produce a meaningful result!
%
%
%     Input:
%       -img1: unaligned image
%       -img2: reference panorama for alignment
%       -bVisualization: if it set to 1 this will show the result of
%                        minimization
%
%     Output:
%       -img1_rot: img1 rotated to be aligned to img2
%       -rot: rotation in pixel. Img1 needs to be shifted of rot pixels in
%             order to be aligned to img2. For the rotation use imShiftWrap.m
%       -err: matching error
%
%     Copyright (C) 2012-16  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

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

c = size(img1, 2);

%Calculate the descriptor
Lwa_ext = (logMean(lum(img1)) + logMean(lum(img2))) / 2.0;

img1_tmo = ReinhardTMO(img1, 0.15, 1e9, 'global', [], Lwa_ext);
img2_tmo = ReinhardTMO(img2, 0.15, 1e9, 'global', [], Lwa_ext);

line1 = LLDescriptor(img1_tmo, 1)';
line2 = LLDescriptor(img2_tmo, 1)';

%minimization of the rotation
rot = 0;
err = sum((line1 - line2).^2);

for i=1:(c - 1)
    line1 = circshift(line1, 1);
    tmpErr = sum((line1 - line2).^2);

    if(tmpErr < err)
        rot = i;
        err = tmpErr;
    end
end

%plotting the result
if(bVisualization)
    line1 = circshift(line1, rot);
    plot(1:c, line1, 1:c, line2);
end

img1_rot = imShiftWrap(img1, rot);

end

author: banterle


text:
function p = binary_search(A,t)
%% Binary Search
% This function binary searches target value (t) in sorted (in increasing order) array A.
% Binary search compares the target value to the middle element of the
% array. If they are not equal, it determines the half part of array in which the taget might be existed,
% then changes the search range from left or right to half array range and
% repeat searching for this new range.
% If target can be found in array, this function returns its index.
% If target can not be found in array, it displays "target is not found in
% array"

array_length = length(A);
counter = 0;                                      % number of iteration in searching algorithm
L_SearchRange = 1;                                %initial search range
R_SearchRange = array_length;

while counter <= floor(log2(array_length))+1       %maximum iteration needed to find the target
mid = (L_SearchRange + R_SearchRange)/2;

if t == A(floor(mid))
    p = floor(mid);
    break
else if t > A(floor(mid))
        L_SearchRange = floor(mid)+1;
    else
        R_SearchRange = floor(mid)-1;
        if R_SearchRange == 0            %to stop searching when t is less than the minimum value of array
           counter = counter +1;
        end
    end
 counter = counter+1;
end
end
if counter > floor(log2(array_length))+1
    disp('target is not found in aray')
end
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function S = globMatrixVIFE3DStiff(fun,mesh,fem1,fem2)

%% USAGE: generate stiffness global matrix on a 3D mesh
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
%
% OUTPUTS:
% [IN JN XN] --- triplets of the sparse matrix from regular elements.
% [II JI XI] --- triplets of the sparse matrix from interface elements.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initializaiton
if strcmp(fem1.type,'P1')||strcmp(fem1.type,'DGP1')||strcmp(fem1.type,'CR')
    feEvalBas1 = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem1.type,'DGP2')
    feEvalBas1 = @evalP2Bas3D;
end

if strcmp(fem2.type,'P1')||strcmp(fem2.type,'DGP1')||strcmp(fem2.type,'CR')
    feEvalBas2 = @evalP1Bas3D;
elseif strcmp(fem2.type,'P2')||strcmp(fem2.type,'DGP2')
    feEvalBas2 = @evalP2Bas3D;
end

%% 1. Matrix on noninterface elements
dof1 = fem1.ldof; dof2 = fem2.ldof; nloc = dof1*dof2;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);
AN = fem1.area(ntID);
gxN = fem1.gx(ntID,:); gyN = fem1.gy(ntID,:); gzN = fem1.gz(ntID,:); gw = fem1.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(fun,gxN,gyN,gzN);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);
for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem1.bas(ntID,:,i), gxN, gyN, gzN, [1,0,0]);
    Ibasy{i} = feEvalBas1(fem1.bas(ntID,:,i), gxN, gyN, gzN, [0,1,0]);
    Ibasz{i} = feEvalBas1(fem1.bas(ntID,:,i), gxN, gyN, gzN, [0,0,1]);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem2.bas(ntID,:,j), gxN, gyN, gzN, [1,0,0]);
    Jbasy{j} = feEvalBas2(fem2.bas(ntID,:,j), gxN, gyN, gzN, [0,1,0]);
    Jbasz{j} = feEvalBas2(fem2.bas(ntID,:,j), gxN, gyN, gzN, [0,0,1]);
end

IN = reshape(repmat(fem1.t(ntID,:),4,1),nloc*ntN,1);
JN = repmat(reshape(fem2.t(ntID,:),dof2*ntN,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibasx{i}.*(coefN.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coefN.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coefN.*Jbasz{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
Ndof = size([mesh.p;mesh.eIntP],1);
S = sparse(IN(ID),JN(ID),XN(ID),Ndof,Ndof);

repo: ifem


text:
function rVectECEF = getrVectEcefFromLatLongAlt(lat, long, alt, bodyInfo)
    r = bodyInfo.radius + alt;

    x = r.*cos(lat).*cos(long);
    y = r.*cos(lat).*sin(long);
    z = r.*sin(lat);

    rVectECEF = [x(:)';y(:)';z(:)'];
end


text:
function fem1d_bvp_linear_test09 ( )

%*****************************************************************************80
%
%% FEM1D_BVP_LINEAR_TEST09 carries out test case #9.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/fem1d_bvp_linear/fem1d_bvp_linear_test09.m
%
%  Discussion:
%
%    Use A9, C9, F9, EXACT9, EXACT_UX9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne O'Leary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%    LC: QA401.O44.
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_BVP_LINEAR_TEST09\n' );
  fprintf ( 1, '  Solve -( A(x) U''(x) )'' + C(x) U(x) = F(x)\n' );
  fprintf ( 1, '  for 0 < x < 1, with U(0) = U(1) = 0.\n' );
  fprintf ( 1, '  A9(X)  = 1.0\n' );
  fprintf ( 1, '  C9(X)  = 0.0\n' );
  fprintf ( 1, '  F9(X)  = X * ( X + 3 ) * exp ( X ),   X <= 2/3\n' );
  fprintf ( 1, '         = 2 * exp ( 2/3),                   2/3 < X\n' );
  fprintf ( 1, '  U9(X)  = X * ( 1 - X ) * exp ( X ),   X <= 2/3\n' );
  fprintf ( 1, '         = X * ( 1 - X ),                    2/3 < X\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', n );
%
%  Geometry definitions.
%
  x_first = 0.0;
  x_last = 1.0;
  x = linspace ( x_first, x_last, n );
  x = x(:);

  u = fem1d_bvp_linear ( n, @a9, @c9, @f9, x );

  uexact = exact9 ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    X         U         Uexact    Error\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %4d  %8f  %8f  %8f  %8e\n', ...
      i, x(i), u(i), uexact(i), abs ( u(i) - uexact(i) ) );
  end

  e1 = l1_error ( n, x, u, @exact9 );
  e2 = l2_error_linear ( n, x, u, @exact9 );
  h1s = h1s_error_linear ( n, x, u, @exact_ux9 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 norm of error  = %g\n', e1 );
  fprintf ( 1, '  L2 norm of error  = %g\n', e2 );
  fprintf ( 1, '  Seminorm of error = %g\n', h1s  );

  return
end
function value = a9 ( x )

%*****************************************************************************80
%
%% A9 evaluates A function #9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c9 ( x )

%*****************************************************************************80
%
%% C9 evaluates C function #9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 0.0;

  return
end
function value = exact9 ( x )

%*****************************************************************************80
%
%% EXACT9 evaluates exact solution #9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of U(X).
%
  value = x .* ( 1.0 - x ) .* exp ( x )         .* ( x <= 2.0 / 3.0 ) ...
        + x .* ( 1.0 - x )                      .* ( 2.0 / 3.0 < x );

  return
end
function value = exact_ux9 ( x )

%*****************************************************************************80
%
%% EXACT_UX9 evaluates the derivative of exact solution #9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of dUdX(X).
%
  value = ( 1.0 - x - x .* x ) .* exp ( x )         .* ( x <= 2.0 / 3.0 ) ...
        + ( 1.0 - 2.0 * x )                         .* (      2.0 / 3.0 < x );

  return
end
function value = f9 ( x )

%*****************************************************************************80
%
%% F9 evaluates right hand side function #9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = x .* ( x + 3.0 ) .* exp ( x )         .* ( x <= 2.0 / 3.0 ) ...
        +               2.0                     .* (      2.0 / 3.0 < x );

  return
end


text:
% RM_histogram Computes the frequency histogram of the row vector x.
%   [RESULT,DESCRIPTOR] = HISTOGRAM(X) or
%   [RESULT,DESCRIPTOR] = HISTOGRAM(X,DESCRIPTOR) or
% where
%   DESCRIPTOR = [LOWER,UPPER,NCELL]
%
%   RESULT    : A row vector containing the histogram
%   DESCRIPTOR: The used descriptor
%
%   X         : The row vector be analyzed
%   DESCRIPTOR: The descriptor of the histogram
%     LOWER   : The lowerbound of the histogram
%     UPPER   : The upperbound of the histogram
%     NCELL   : The number of cells of the histogram
%
%   See also: http://www.cs.rug.nl/~rudy/matlab/
%
%   R. Moddemeijer
%   Copyright (c) by R. Moddemeijer
%   $Revision: 1.2 $  $Date: 2001/02/05 09:54:29 $

function [result,descriptor] = RM_histogram(x,descriptor)

if nargin < 1
   disp('Usage: RESULT = HISTOGRAM(X)')
   disp('       RESULT = HISTOGRAM(X,DESCRIPTOR)')
   disp('Where: DESCRIPTOR = [LOWER,UPPER,NCELL]')
   return
end

% Some initial tests on the input arguments

[NRowX, NColX] = size(x);

if NRowX ~= 1
  error('Invalid dimension of X');
end;

if nargin > 2
  error('Too many arguments');
end;

if nargin==1
  minx=min(x);
  maxx=max(x);
  delta=(maxx-minx)/(length(x)-1);
  ncell=ceil(sqrt(length(x)));
  descriptor=[minx-delta/2,maxx+delta/2,ncell];
end;

lower=descriptor(1);
upper=descriptor(2);
ncell=descriptor(3);

if ncell<1
  error('Invalid number of cells')
end;

if upper<=lower
  error('Invalid bounds')
end;

result(1:ncell)=0;

y=round( (x-lower)/(upper-lower)*ncell + 1/2 );
for n=1:NColX
  index=y(n);
  if index >= 1 & index<=ncell
    result(index)=result(index)+1;
  end;
end;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function d = duplicated(x)
%DUPLICATED   Find duplicated rows.
% DUPLICATED(x) returns a vector d such that d(i) = 1 if x(i,:) is a
% duplicate of an earlier row.
%
% Examples:
% duplicated([2 7 8 7 1 2 8]') = [0 0 0 1 0 1 1]'
% duplicated([0 0 1 1 0; 0 1 0 1 1]') = [0 0 0 0 1]'
% duplicated(eye(100))
% duplicated(kron((1:3)',ones(3)))
%
% You can simulate unique(x) or unique(x,'rows') by x(~duplicated(x)).
% The difference is that the latter form will not sort the contents of x,
% as unique will.
%
% See also UNIQUE.

% (c) Microsoft Corporation. All rights reserved.

% This function is not well optimized.
% In particular, it is slower than unique.

[nr,nc] = size(x);
if nc == 1
  d = duplicated1(x);
  return;
end
if nr == 1
  d = 0;
  return;
end
hash = x*rand(nc,1);
[dummy,ord] = sort(hash);
xo = x(ord,:);
%d = [0 all(diff(xo)==0,2)'];
d = [0 all(xo(1:end-1,:)==xo(2:end,:),2)'];
dd = diff([d 0]);
dstart = find(dd > 0);
dend = find(dd < 0);
% loop each run of duplicated columns
for i = 1:length(dstart)
  % place the zero at the first element in the original order
  d(dstart(i)) = 1;
  d(dstart(i)-1 + argmin(ord(dstart(i):dend(i)))) = 0;
end
d(ord) = d;
d = d';
%d = duplicated1(hash);

function d = duplicated1(x)
% special case where x is a column vector.

[s,ord] = sort(x);
d = zeros(size(x));
d(ord) = [0; s(1:end-1)==s(2:end)];
%d(ord) = [0; diff(s)==0];


text:
function [ NITF_meta ] = read_sicd_nitf_offsets( filename )
%READ_SICD_NITF_OFFSETS Read metadata from Sensor Independent Complex Data
%(SICD) file, version 0.3 and above
%
% SICD (version >= 0.3) is stored in a NITF container.  NITF is a
% complicated format that involves lots of fields and configurations
% possibilities. Fortunately, SICD only really uses a small, specific
% portion of the NITF format.  This function extracts only the few parts of
% the NITF metadata necessary for reading a SICD NITF file.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Open file
fid = fopen(filename);

%% Read NITF file header
if ~strcmp(fread(fid,9,'uint8=>char')','NITF02.10') % Check format
    fclose(fid);
    error('READ_SICD_META:INVALID_FILE_FORMAT','Only NITF version 2.1 recognized.');
end
fseek(fid,354,'bof'); % Offset to first field of interest
HL = str2double(fread(fid,6,'uint8=>char')'); % File header length
NUMI = str2double(fread(fid,3,'uint8=>char')'); % Number of image segments
if NUMI > 0
    ImgSubhdrLngths=zeros(NUMI,1); % Length of the image segment subheader
    ImgLngths=zeros(NUMI,1); % Length of the image segment data
    NITF_meta.minimal.imageSegmentOffsets=zeros(NUMI,1); % Offset to image data from beginning of file (in bytes)
    NITF_meta.minimal.imageSegmentRows=zeros(NUMI,1); % Number of rows in each image segment (in case data is spread across multiple image segments)
    NITF_meta.minimal.imageSegmentColumns=zeros(NUMI,1); % Number of rows in each image segment (in case data is spread across multiple image segments)
    for i = 1:NUMI
        ImgSubhdrLngths(i) = str2double(fread(fid,6,'uint8=>char')');
        NITF_meta.minimal.imageSegmentOffsets(i)=HL+sum(ImgSubhdrLngths)+sum(ImgLngths);
        ImgLngths(i) = str2double(fread(fid,10,'uint8=>char')');
    end
else
    NITF_meta.minimal.imageSegmentOffsets=[];
    NITF_meta.minimal.imageSegmentRows=[];
    NITF_meta.minimal.imageSegmentColumns=[];
end
for i=1:2 % Jump over segments not used by SICD to DES description
    segment_length = str2double(fread(fid,3,'uint8=>char')');
    if(isfinite(segment_length)&&(segment_length>0))
        error('READ_SICD_META:INVALID_NITF_SEGMENT','SICD does not allow for graphics or reserved extension segments.');
    end
end
% Jump over text segments
NUMT = str2double(fread(fid,3,'uint8=>char')');
LTSH = zeros(NUMT,1); LT = zeros(NUMT,1);
for lp = 1:NUMT
    LTSH(lp) = str2double(fread(fid,4,'uint8=>char')'); % Length of text segment subheader
    LT(lp) = str2double(fread(fid,5,'uint8=>char')'); % Length of text segment data
end
NUMDES = str2double(fread(fid,3,'uint8=>char')');  % Number of data extension segments
LDSH = zeros(NUMDES,1); LD = zeros(NUMDES,1);
for lp = 1:NUMDES
    LDSH(lp) = str2double(fread(fid,4,'uint8=>char')'); % Length of data extension segment subheader
    LD(lp) = str2double(fread(fid,9,'uint8=>char')'); % Length of data extension segment data
end
NITF_meta.minimal.desLengths = LD;
NITF_meta.minimal.desOffsets = HL + sum(ImgSubhdrLngths) + sum(ImgLngths) + ...
    sum(LTSH) + sum(LT) + cumsum(LDSH) + [0; cumsum(LD(1:(end-1)))];

%% Get number of rows for each image segment from image segment headers
nextImgSubhdroffset=HL;
for i = 1:NUMI
    fseek(fid,nextImgSubhdroffset,'bof'); % Jump to ith image segment
    fseek(fid,333,'cof'); % Jump to number of rows field
    NITF_meta.minimal.imageSegmentRows(i)=str2double(fread(fid,8,'uint8=>char')');
    NITF_meta.minimal.imageSegmentColumns(i)=str2double(fread(fid,8,'uint8=>char')');
    nextImgSubhdroffset=nextImgSubhdroffset+ImgSubhdrLngths(i)+ImgLngths(i);
end

%% Close file
fclose(fid);

% If you want the full NITF header data, you can use the NITF functions
% from any other NITF toolbox, but its not really necessary for anything in
% SICD.
% NITF_meta.IPT = nitfinfo(filename); % MATLAB Image Processing Toolbox
% NITF_meta.matlab_sar = read_nitf_meta(filename); % NITF reading in the MATLAB SAR Toolbox

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function dtiXformRoiToMrVistaVolume(handles, roiNum)
%
% mrVistaFibers = dtiXformRoiToMrVistaVolume(handles, [roiNum])
%
% Uses the xformVAnatToAcpc (see dtiXformVanatCompute) to convert the specified
% ROI to mrVista vAnatomy coords.
%
% HISTORY:
% 2004.05.19 RFD (bob@white.stanford.edu) wrote it.

if(~isfield(handles, 'xformVAnatToAcpc') | isempty(handles.xformVAnatToAcpc))
    error('xformVAnatToAcpc not computed yet- run "Compute mrVista Xform" from the "Xform" menu.');
end
if(~exist('roiNum','var') | isempty(roiNum))
    roiNum = handles.curRoi;
end

view = getSelectedVolume;

for(ii=1:length(roiNum))
    roi = handles.rois(roiNum(ii));
    coords = unique(round(mrAnatXformCoords(inv(handles.xformVAnatToAcpc),roi.coords)),'rows')';
    view = newROI(view, ['dti_' roi.name], 1);
    view.ROIs(view.selectedROI).coords = coords;
end

% Is there a better way to do this?
mrGlobals;
eval([view.name '=view;']);
return;


text:
function Y=EWMASTD(X,d)

%   Y=EWMASTD(X,d) returns the EWMA (Exponentially Weighted Moving Average)
%   standard deviation using the historical returns in vector X and a decay
%   factor, d.
% % ======================================================================
%
%   Author: Lorenzo Brancali
%   E-mail: lbrancali@gmail.com
%   Date:   20th Febryary 2012
%
% % ======================================================================

t=length(X);

for i=1:length(X)-1

    F(i,:)=((d^(i-1))*(X(t-i,:)-mean(X))^2);

end

Y=sqrt((1-d)*sum(F));


text:
function [Ymat Xmat N n m p T k q h]=panel5prelim(data_endo,data_exo,const,lags)

% first compute N, the number of units, as the dimension of the data_endo matrix
N=size(data_endo,3);

% compute p, the number of lags in the model
p=lags;

% then compute n, the number of endogenous variables in the model; it is simply the number of columns in the matrix 'data_endo'
n=size(data_endo,2);

% if the constant has been selected, augment the matrix of exogenous with a column of ones (number of rows equal to number of rows in data_endo)
if const==1
data_exo=[ones(size(data_endo,1),1) data_exo];
% if no constant was included, do nothing
else
end

% compute m, the number of exogenous variables in the model
% if data_exo is empty, set m=0
if isempty(data_exo)==1
m=0;
% if data_exo is not empty, count the number of exogenous variables that will be included in the model
else
m=size(data_exo,2);
% Also, trim a number initial rows equal to the number of lags, as they will be suppressed from the endogenous as well to create initial conditions
data_exo=data_exo(p+1:end,:);
end

% determine k, the number of parameters to estimate in each equation; it is equal to np+m
k=N*n*p+m;

% determine q, the total number of VAR parameters for each unit
q=n*k;

% determine h, the total number of VAR parameters for the whole model
h=N*q;

% obtain Ymat and Xmat
temp=[];
% stack the matrices of endogenous variables to obtain a temporary matrix
for ii=1:N
temp=[temp data_endo(:,:,ii)];
end
% use the lagx function on this matrix
temp=bear.lagx(temp,lags);

% set Ymat as the first Nn columns of the result
Ymat=temp(:,1:N*n);

% to build Xmat, take off the Nn initial columns of temp, and concatenate the exogenous
Xmat=[temp(:,N*n+1:end) data_exo];

% Define T, the number of periods of the model, as the number of rows of X
T=size(Xmat,1);

author: european-central-bank

repo: BEAR-toolbox


text:
classdef MetaDataTagsClass  < FileLoadSaveClass

    properties
        tags
    end

    methods

        % -------------------------------------------------------
        function obj = MetaDataTagsClass(varargin)
            % Set class properties not part of the SNIRF format
            obj.SetFileFormat('hdf5');
            obj.tags.SubjectID = 'default';
            obj.tags.MeasurementDate = datestr(now,29);
            obj.tags.MeasurementTime = datestr(now,'hh:mm:ss');
            obj.tags.LengthUnit = 'mm';
            obj.tags.TimeUnit = 'unknown';
            obj.tags.FrequencyUnit = 'unknown';
            obj.tags.AppName  = 'homer3-DataTree';

            if nargin==1 && ~isempty(varargin{1})
                obj.SetFilename(varargin{1});
                obj.Load();
            end
            if nargin==2
                obj.tags.LengthUnit = varargin{2};
            end
        end

        % -------------------------------------------------------
        function err = LoadHdf5(obj, fileobj, location)
            err = 0;

            % Arg 1
            if ~exist('fileobj','var') || (ischar(fileobj) && ~exist(fileobj,'file'))
                fileobj = '';
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/metaDataTags';
            elseif location(1)~='/'
                location = ['/',location];
            end

            % Error checking
            if ~isempty(fileobj) && ischar(fileobj)
                obj.SetFilename(fileobj);
            elseif isempty(fileobj)
                fileobj = obj.GetFilename();
            end
            if isempty(fileobj)
               err = -1;
               return;
            end

            %%%%%%%%%%%% Ready to load from file

            try
                % Reset tags
                obj.tags = struct();

                % Open group
                [gid, fid] = HDF5_GroupOpen(fileobj, location);

                metaDataStruct = h5loadgroup(gid);
                tags = fieldnames(metaDataStruct); %#ok<*PROPLC>
                for ii=1:length(tags)
                    eval(sprintf('obj.tags.%s = metaDataStruct.%s;', tags{ii}, tags{ii}));
                end

                HDF5_GroupClose(fileobj, gid, fid);

                % Detect old or invalid metadatatag data
                assert(obj.IsValid());

            catch

                err = -1;

            end
            obj.SetError(err);

        end

        % -------------------------------------------------------
        function SaveHdf5(obj, fileobj, location) %#ok<*INUSD>
            % Arg 1
            if ~exist('fileobj', 'var') || isempty(fileobj)
                error('Unable to save file. No file name given.')
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/metaDataTags';
            elseif location(1)~='/'
                location = ['/',location]; %#ok<*NASGU>
            end

            if ~exist(fileobj, 'file')
                fid = H5F.create(fileobj, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
                H5F.close(fid);
            end
            props = propnames(obj.tags);
            for ii = 1:length(props)
                eval(sprintf('hdf5write_safe(fileobj, [location, ''/%s''], obj.tags.%s);', props{ii}, props{ii}));
            end
        end

        % -------------------------------------------------------
        function b = IsValid(obj)
            b = false;

            % Use latest required fields to determine if we're loading old
            % metaDataTag format version of SNIRF spec
            if ~isproperty(obj.tags, 'FrequencyUnit')
                return;
            end

            b = true;
        end

        % -------------------------------------------------------
        function B = eq(obj, obj2)
            B = false;
            props1 = propnames(obj.tags);
            props2 = propnames(obj2.tags);
            for ii=1:length(props1)
                if ~isproperty(obj2.tags, props1{ii})
                    return;
                end
                if eval(sprintf('~strcmp(obj.tags.%s, obj2.tags.%s)', props1{ii}, props1{ii}))
                    return;
                end
            end
            for ii=1:length(props2)
                if ~isproperty(obj.tags, props2{ii})
                    return;
                end
                if eval(sprintf('~strcmp(obj.tags.%s, obj2.tags.%s)', props2{ii}, props2{ii}))
                    return;
                end
            end
            B = true;
        end

        % -------------------------------------------------------
        function Add(obj, key, value) %#ok<INUSL>
            key(key==' ') = '';
            eval(sprintf('obj.tags.%s = value', key));
        end

        % ----------------------------------------------------------------------------------
        function val = Get(obj, name)
            val = [];
            if ~exist('name', 'var')
                return;
            end
            if isfield(obj.tags, name)
                val = eval( sprintf('obj.tags.%s;', name) );
            end
        end

        % ----------------------------------------------------------------------------------
        function Set(obj, name, value)
            if ~exist('name', 'var') || ~exist('value', 'var')
                retrun
            end
            eval( sprintf('obj.tags.%s = ''%s'';', name, value) )
        end

        % ----------------------------------------------------------------------------------
        function SetLengthUnit(obj, unit)
            if isempty(obj)
                return
            end
            obj.tags.LengthUnit = unit;
        end

        % ----------------------------------------------------------------------------------
        function val = GetLengthUnit(obj)
            val = '';
            if isempty(obj)
                return
            end
            val = obj.tags.LengthUnit;
        end

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            fields = propnames(obj.tags);
            for ii=1:length(fields)
                nbytes = nbytes + eval(sprintf('sizeof(obj.tags.%s)', fields{ii}));
            end
        end

    end
end

author: BUNPC


text:
function inspect_bug645

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_event

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug645/Events_ECI_TCPIP_55513.xml');

% there is nothing in this script to be tested automatically
% it only serves to determine the timings once

cd /Volumes/Data/roboos/matlab/fieldtrip/fileio/private
tic; ver1 = xml2struct(filename); toc
% xml2struct reading /home/common/matlab/fieldtrip/data/test/bug645/Events_ECI_TCPIP_55513.xml
% Elapsed time is 22.311995 seconds.

cd /Volumes/Data/roboos/matlab/fieldtrip/external/xml4mat
tic; ver2 = xml2struct(filename); toc
% Elapsed time is 404.807848 seconds.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cvx_optval,P,q,r,X,lambda] = cheb(A,b,Sigma) %#ok

% Computes Chebyshev lower bounds on probability vectors
%
% Calculates a lower bound on the probability that a random vector
% x with mean zero and covariance Sigma satisfies A x <= b
%
% Sigma must be positive definite
%
% output arguments:
% - prob: lower bound on probability
% - P,q,r: x'*P*x + 2*q'*x + r is a quadratic function
%   that majorizes the 0-1 indicator function of the complement
%   of the polyhedron,
% - X, lambda:  a discrete distribution with mean zero, covariance
%   Sigma and Prob(X not in C)  >= 1-prob

%
% maximize  1 - Tr Sigma*P - r
% s.t.      [ P  q     ]             [ 0      a_i/2 ]
%           [ q' r - 1 ] >= tau(i) * [ a_i'/2  -b_i ], i=1,...,m
%           taui >= 0
%           [ P q  ]
%           [ q' r ] >= 0
%
% variables P in Sn, q in Rn, r in R
%

[ m, n ] = size( A ); %#ok
cvx_begin sdp quiet
    variable P(n,n) symmetric
    variables q(n) r tau(m)
    dual variables Z{m}
    maximize( 1 - trace( Sigma * P ) - r )
    subject to
        for i = 1 : m,
            qadj = q - 0.5 * tau(i) * A(i,:)';
            radj = r - 1 + tau(i) * b(i);
            [ P, qadj ; qadj', radj ] >= 0 : Z{i}; %#ok
        end
        [ P, q ; q', r ] >= 0; %#ok
        tau >= 0; %#ok
cvx_end

if nargout < 4,
    return
end

X = [];
lambda = [];
for i=1:m
   Zi = Z{i};
   if (abs(Zi(3,3)) > 1e-4)
      lambda = [lambda; Zi(3,3)]; %#ok
      X = [X Zi(1:2,3)/Zi(3,3)]; %#ok
   end;
end;
mu = 1-sum(lambda);
if (mu>1e-5)
   w = (-X*lambda)/mu;
   W = (Sigma - X*diag(lambda)*X')/mu;
   [v,d] = eig(W-w*w');
   d = diag(d);
   s = sum(d>1e-5);
   if (d(1) > 1e-5)
      X = [X w+sqrt(s)*sqrt(d(1))*v(:,1) ...
            w-sqrt(s)*sqrt(d(1))*v(:,1)];
      lambda = [lambda; mu/(2*s); mu/(2*s)];
   elseif (d(2) > 1e-5)
      X = [X w+sqrt(s)*sqrt(d(2))*v(:,2) ...
            w-sqrt(s)*sqrt(d(2))*v(:,2)];
      lambda = [lambda; mu/(2*s); mu/(2*s)];
   else
      X = [X w];
      lambda = [lambda; mu];
   end;
end;

author: yu-jiang

repo: radpbook

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