text:
classdef clusteringGUI < handle

    properties
        currentCluster = 1
        page = 1
        thumbnail_size = [200 200]
        clustAssign
        clusters
        rejected
        ClusteringData
        minfreq
        maxfreq
        fig
        image_axes = gobjects()
        handle_image = gobjects()
        ColorData
        totalCount
        count
        clusterName
        pagenumber
        finished
        call_id_text
        txtbox
    end

    methods
        function [obj, NewclusterName, NewRejected, NewFinished, NewClustAssign] = clusteringGUI(clustAssign, ClusteringData)

            obj.clustAssign = clustAssign;
            %Image, Lower freq, delta time, Time points, Freq points, File path, Call ID in file, power, RelBox
            obj.ClusteringData = ClusteringData;
            obj.rejected = zeros(1,length(obj.clustAssign));

            obj.minfreq = prctile(ClusteringData.MinFreq, 5);
            obj.maxfreq = prctile(ClusteringData.MinFreq + ClusteringData.Bandwidth, 95);
            obj.ColorData = jet(256); % Color by mean frequency
            % obj.ColorData = HSLuv_to_RGB(256, 'H',  [270 0], 'S', 100, 'L', 75, 'type', 'HSL'); % Make a color map for each category
            obj.ColorData = reshape(obj.ColorData,size(obj.ColorData,1),1,size(obj.ColorData,2));

            if iscategorical(obj.clustAssign)
                obj.clusterName =unique(obj.clustAssign);
                obj.clusters = unique(obj.clustAssign);
            else
                obj.clusterName = categorical(unique(obj.clustAssign(~isnan(obj.clustAssign))));
                obj.clusters = (unique(obj.clustAssign(~isnan(obj.clustAssign))));
            end

            obj.fig = dialog('Visible','off','Position',[360,500,600,600],'WindowStyle','Normal','resize', 'on','WindowState','maximized' );
            obj.fig.CloseRequestFcn = @(src,event) finished_Callback(obj, src, event);
            set(obj.fig,'color',[.1, .1, .1]);

            movegui(obj.fig,'center');
            %             set(obj.fig,'WindowButtonMotionFcn', @(hObject, eventdata) mouse_over_Callback(obj, hObject, eventdata));

            txt = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.1 .1 .1],...
                'ForegroundColor','w',...
                'Style','text',...
                'Position',[120 565 80 30],...
                'String','Name:');

            obj.txtbox = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.149 .251 .251],...
                'ForegroundColor','w',...
                'Style','edit',...
                'String','',...
                'Position',[120 550 80 30],...
                'Callback',@(src,event) txtbox_Callback(obj,src,event));

            obj.totalCount = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.1 .1 .1],...
                'ForegroundColor','w',...
                'Style','text',...
                'String','',...
                'Position',[330 542.5 200 30],...
                'HorizontalAlignment','left');

            back = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.149 .251 .251],...
                'ForegroundColor','w',...
                'Position',[20 550 80 30],...
                'String','Back',...
                'Callback',@(src,event) back_Callback(obj, src, event));

            next = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.149 .251 .251],...
                'ForegroundColor','w',...
                'Position',[220 550 80 30],...
                'String','Next',...
                'Callback',@(src,event) next_Callback(obj, src, event));

            apply = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.149 .251 .251],...
                'ForegroundColor','w',...
                'Position',[440 550 60 30],...
                'String','Save',...
                'Callback',@(src,event)  finished_Callback(obj, src, event));

            if nargin == 2
                redo = uicontrol('Parent',obj.fig,...
                    'BackgroundColor',[.149 .251 .251],...
                    'ForegroundColor','w',...
                    'Position',[510 550 60 30],...
                    'String','Redo',...
                    'Callback',@(src,event) finished_Callback(obj, src, event));
            else
                redo = uicontrol('Parent',obj.fig,...
                    'BackgroundColor',[.149 .251 .251],...
                    'ForegroundColor','w',...
                    'Position',[510 550 60 30],...
                    'String','Cancel',...
                    'Callback',@(src,event) finished_Callback(obj, src, event));
            end
            %% Paging
            nextpage = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.149 .251 .251],...
                'ForegroundColor','w',...
                'Position',[220 517 80 30],...
                'String','Next Page',...
                'Callback',@(src,event) nextpage_Callback(obj, src, event));

            backpage = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.149 .251 .251],...
                'ForegroundColor','w',...
                'Position',[20 517 80 30],...
                'String','Previous Page',...
                'Callback',@(src,event, h) backpage_Callback(obj, src, event));

            obj.pagenumber = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.1 .1 .1],...
                'ForegroundColor','w',...
                'Style','text',...
                'String','',...
                'Position',[118 509 80 30],...
                'HorizontalAlignment','center');

            obj.call_id_text = uicontrol('Parent',obj.fig,...
                'BackgroundColor',[.1 .1 .1],...
                'ForegroundColor','w',...
                'Style','text',...
                'String','',...
                'FontSize',12,...
                'Position',[100 470 400 30],...
                'HorizontalAlignment','center');

            obj.render_GUI();

            % Wait for d to close before running to completion
            set( findall(obj.fig, '-property', 'Units' ), 'Units', 'Normalized');
            obj.fig.Visible = 'on';

            % Enable pointer management for the figure for mouse hover over
            iptPointerManager(obj.fig, 'enable');

            uiwait(obj.fig);
            NewclusterName = obj.clusterName;
            NewRejected = obj.rejected;
            NewFinished = obj.finished;
            NewClustAssign = obj.clustAssign;

        end

        function render_GUI(obj)

            %% Colormap
            xdata = obj.minfreq:obj.maxfreq;
            caxis = axes(obj.fig,'Units','Normalized','Position',[.88 .05 .04 .8]);
            image(1,xdata,obj.ColorData,'parent',caxis)
            caxis.YDir = 'normal';
            set(caxis,'YColor','w','box','off','YAxisLocation','right');
            ylabel(caxis, 'Frequency (kHz)')

            %% Make the axes
            aspectRatio = median(cellfun(@(im) size(im,1) ./ size(im,2), obj.ClusteringData.Spectrogram));

            % Choose a number of rows and columns to fill the space with
            % the average call aspect ratio
            % nFrames = 10;
            % figureAspectRatio = 1;
            % x_grids = sqrt(aspectRatio * figureAspectRatio * nFrames);
            % x_grids = ceil(x_grids);
            % y_grids = ceil(nFrames / x_grids);

            obj.thumbnail_size = round(sqrt(20000 .* [aspectRatio, 1/aspectRatio]));

            axes_spacing = .70; % Relative width of each image
            y_range = [.05, .75]; % [Start, End] of the grid
            x_range = [.05, .85];
            x_grids = 9; % Number of x grids
            y_grids = 3; % Number of y grids

            ypos = linspace(y_range(1), y_range(2) - axes_spacing * range(y_range) / y_grids, y_grids );
            xpos = linspace(x_range(1), x_range(2) - axes_spacing * range(x_range) / x_grids, x_grids );
            xpos = fliplr(xpos);

            pos = [];
            for i = 1:length(ypos)
                for j = 1:length(xpos)
                    pos(end+1,:) = [xpos(j), ypos(i), (xpos(1)-xpos(2)) * axes_spacing, (ypos(2)-ypos(1)) * axes_spacing];
                end
            end
            pos = flipud(pos);
            for i = 1 : length(ypos) * length(xpos)
                    im = zeros([obj.thumbnail_size, 3]);
                    obj.image_axes(i) = axes(obj.fig,'Units','Normalized','Position',pos(i,:));
                    obj.handle_image(i) = image(im,'parent',obj.image_axes(i));
                    set(obj.image_axes(i),'Visible','off')
                    set(get(obj.image_axes(i),'children'),'Visible','off');
            end
            plotimages(obj);
        end

        function [colorIM, rel_x, rel_y] = create_thumbnail(obj, ClusteringData,clustIndex,callID)
            % Resize the image while maintaining the aspect ratio by
            % padding with zeros
            im_size = size(ClusteringData.Spectrogram{clustIndex(callID)}) ;
            new_size = floor(im_size .* min(obj.thumbnail_size ./ im_size));
            im = double(imresize(ClusteringData.Spectrogram{clustIndex(callID)}, new_size));
            pad = (obj.thumbnail_size - size(im)) / 2;
            im = padarray(im, floor(pad), 'pre');
            im = padarray(im, ceil(pad), 'post');

            % Relative offsets for setting the tick values
            rel_size = pad ./ obj.thumbnail_size;
            rel_x = [rel_size(2), 1-rel_size(2)];
            rel_y = [rel_size(1), 1-rel_size(1)];

            % Apply color to the greyscale images
            freqRange = [ClusteringData.MinFreq(clustIndex(callID)),...
                ClusteringData.MinFreq(clustIndex(callID)) + ClusteringData.Bandwidth(clustIndex(callID))];
            % Account for any padding on the y axis
            freqRange = freqRange + range(freqRange) .* rel_y(1) .* [-1, 1];

            freqdata = linspace(freqRange(2) ,freqRange(1), obj.thumbnail_size(1));
            colorMask = interp1(linspace(obj.minfreq, obj.maxfreq, size(obj.ColorData,1)), obj.ColorData, freqdata, 'nearest', 'extrap');
            colorIM = im .* colorMask ./ 255;
        end

        function obj = config_axis(obj, axis_handles,i, rel_x, rel_y)
            set(axis_handles,'xcolor','w');
            set(axis_handles,'ycolor','w');

            x_lim = xlim(axis_handles);
            x_span = x_lim(2) - x_lim(1);
            xtick_positions = linspace(x_span*rel_x(1)+x_lim(1), x_span*rel_x(2)+x_lim(1),4);
            x_ticks = linspace(0,obj.ClusteringData.Duration(i),4);
            x_ticks = arrayfun(@(x) sprintf('%.3f',x),x_ticks(2:end),'UniformOutput',false);

            y_lim = ylim(axis_handles);
            y_span = y_lim(2) - y_lim(1);
            ytick_positions = linspace(y_span*rel_y(1)+y_lim(1), y_span*rel_y(2)+y_lim(1),3);

            y_ticks = linspace(obj.ClusteringData.MinFreq(i),obj.ClusteringData.MinFreq(i)+obj.ClusteringData.Bandwidth(i),3);
            y_ticks = arrayfun(@(x) sprintf('%.1f',x),y_ticks(1:end),'UniformOutput',false);
            y_ticks = flip(y_ticks);

            yticks(axis_handles,ytick_positions);
            xticks(axis_handles,xtick_positions(2:end));
            xticklabels(axis_handles,x_ticks);
            yticklabels(axis_handles,y_ticks);
            xlabel(axis_handles,'Time (s)');
            ylabel(axis_handles,'Frequency (kHz)');
        end

        function obj = plotimages(obj)
            % Number of calls in each cluster
            for cl = 1:length(obj.clusterName)
                obj.count(cl) = sum(obj.clustAssign==obj.clusters(cl));
            end

            clustIndex = find(obj.clustAssign==obj.clusters(obj.currentCluster));

            for i=1:length(obj.image_axes)
                if i <= length(clustIndex) - (obj.page - 1)*length(obj.image_axes)
                    % set(image_axes(i),'Visible','off')

                    set(get(obj.image_axes(i),'children'),'Visible','on');

                    callID = i + (obj.page - 1)*length(obj.image_axes);
                    [colorIM, rel_x, rel_y] = obj.create_thumbnail(obj.ClusteringData,clustIndex,callID);
                    set(obj.handle_image(i), 'ButtonDownFcn',@(src,event) clicked(obj,src,event,clustIndex(callID),i,callID));
                    obj.add_cluster_context_menu(obj.handle_image(i),clustIndex(callID));

                    % Display the file ID and call number on mouse hover
                    [~,call_file,~] = fileparts(obj.ClusteringData.Filename(clustIndex(callID)));
                    call_id = sprintf('Call: %u', obj.ClusteringData.callID(clustIndex(callID)));
                    pointerBehavior.enterFcn = @(~,~) set(obj.call_id_text, 'string', {call_id, call_file});
                    pointerBehavior.traverseFcn = [];
                    pointerBehavior.exitFcn = @(~,~) set(obj.call_id_text, 'string', '');
                    iptSetPointerBehavior(obj.handle_image(i), pointerBehavior);

                    % Make the image red if the call is rejected
                    if obj.rejected(clustIndex(callID))
                        colorIM(:,:,1) = colorIM(:,:,1) + .5;
                    end

                    set(obj.handle_image(i),'CData',colorIM, 'XData', []);

                    obj.config_axis(obj.image_axes(i),clustIndex(callID), rel_x, rel_y);

                    set(obj.image_axes(i),'Visible','on')

                else
                    set(obj.image_axes(i),'Visible','off')
                    set(get(obj.image_axes(i),'children'),'Visible','off');
                end

            end

            % Update text
            obj.pagenumber.String = sprintf('Page %u of %u', obj.page, ceil(obj.count(obj.currentCluster) / length(obj.image_axes)));
            obj.txtbox.String = string(obj.clusterName(obj.currentCluster));
            obj.totalCount.String = sprintf('total count: %u', obj.count(obj.currentCluster));
            obj.fig.Name = sprintf('Cluster %u of %u', obj.currentCluster, length(obj.count));

        end

        function obj = add_cluster_context_menu(obj, hObject, i)
            unique_clusters = unique(obj.clusterName);

            c = uicontextmenu(obj.fig);
            for ci=1:length(unique_clusters)
                uimenu(c,'text',string(obj.clusterName(ci)),'Callback',@(src,event) assign_cluster(obj, src, event,i,unique_clusters(ci)));
            end

            set(hObject, 'UIContextMenu',c);
        end

        function obj = assign_cluster(obj, hObject,eventdata,i, clusterLabel)
            obj.clustAssign(i) = clusterLabel;
            obj.plotimages();
        end

        function obj = clicked(obj, hObject,eventdata,i,plotI,callID)
            if( eventdata.Button ~= 1 ) % Return if not left clicked
                return
            end

            clustIndex = find(obj.clustAssign == obj.clusters(obj.currentCluster));

            obj.rejected(i) = ~obj.rejected(i);

            [colorIM, ~, ~] = obj.create_thumbnail(obj.ClusteringData,clustIndex,callID);

            if obj.rejected(i)
                colorIM(:,:,1) = colorIM(:,:,1) + .5;
            end
            set(obj.handle_image(plotI),'CData',colorIM);
        end

        function obj = next_Callback(obj, hObject, eventdata)
            obj.clusterName(obj.currentCluster) = get(obj.txtbox,'String');
            if obj.currentCluster < length(obj.clusterName)
                obj.currentCluster = obj.currentCluster + 1;
                obj.page = 1;
                obj.plotimages();
            end
        end

        function obj = back_Callback(obj, hObject, eventdata)
            obj.clusterName(obj.currentCluster) = get(obj.txtbox,'String');
            if obj.currentCluster > 1
                obj.currentCluster = obj.currentCluster-1;
                obj.page = 1;
                obj.plotimages();
            end
        end

        function obj = nextpage_Callback(obj, hObject, eventdata)
            if obj.page < ceil(obj.count(obj.currentCluster) / length(obj.image_axes))
                obj.page = obj.page + 1;
                obj.plotimages();
            end
        end

        function obj = backpage_Callback(obj, hObject, eventdata)
            if obj.page > 1
                obj.page = obj.page - 1;
                obj.plotimages();
            end
        end

        function obj = txtbox_Callback(obj, hObject, eventdata)
            obj.clusterName(obj.currentCluster) = get(hObject,'String');
        end

        function obj = finished_Callback(obj, hObject, eventdata)
            % If window is closed, finished = 2
            % If clicked apply, finished = 1
            % If clicked redo, finished = 0
            switch eventdata.EventName
                case 'Close'
                    obj.finished = 2;
                otherwise
                    switch hObject.String
                        case 'Save'
                            obj.finished = 1;
                        case 'Redo'
                            obj.finished = 0;
                    end
            end
            set(obj.fig,  'closerequestfcn', '');
            delete(obj.fig);
            obj.fig = [];
        end

    end
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function out = full(fa)
% Convert to numeric form
% FORMAT full(fa)
% fa - a file_array
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

%
% $Id: full.m 7501 2018-11-30 12:16:58Z guillaume $

[vo{1:ndims(fa)}] = deal(':');
out = subsref(fa,struct('type','()','subs',{vo}));

author: spm


text:
function hdrdat = fmri_lddat2(datfile)
%
% Load information and parameters from a data file into an hdrdata
% structure. To replace fmri_lddat()
%
% datfile - (string) name of data file
%
% $Id: fmri_lddat2.m,v 1.1 2005/08/24 22:19:57 sayres Exp $

if(nargin ~= 1)
  msg = 'Usage: hdrdat = fmri_lddat2(datfile)'
  qoe(msg); error(msg);
end

fid=fopen(deblank(datfile),'r');
if( fid == -1 )
  msg = sprintf('Could not open dof file %s\n',datfile);
  qoe(msg); error(msg);
end

hdrdat = fmri_hdrdatstruct;

%% --- selavg stuff --- %%
hdrdat.TR      = fscanf2(fid,'%f');
hdrdat.TimeWindow      = fscanf2(fid,'%f');
hdrdat.TPreStim     = fscanf2(fid,'%f');
nPS     = floor(hdrdat.TPreStim/hdrdat.TR);
hdrdat.Nc   = fscanf2(fid,'%d');
hdrdat.Nnnc    = hdrdat.Nc - 1;
hdrdat.Nh   = fscanf2(fid,'%d'); % nperevent

%% --- selxavg stuff --- %%
hdrdat.Version = fscanf2(fid,'%d');

if(isempty(hdrdat.Version))
  %fprintf('    INFO: Data file in selavg format\n');
  hdrdat.Version = 0;
  hdrdat.TER     = hdrdat.TR;
  hdrdat.DOF     = [];
  return;
end

hdrdat.TER   = fscanf2(fid,'%f');
hdrdat.DOF   = fscanf2(fid,'%f');

dummy = fscanf(fid,'%s',1);
hdrdat.Npercond = fscanf(fid,'%d',hdrdat.Nc);

hdrdat.Nruns = fscanf2(fid,'%d');
hdrdat.Ntp   = fscanf2(fid,'%d');
hdrdat.Nrows = fscanf2(fid,'%d');
hdrdat.Ncols = fscanf2(fid,'%d');
hdrdat.Nskip = fscanf2(fid,'%d');
hdrdat.DTOrder  = fscanf2(fid,'%d');
hdrdat.RescaleFactor  = fscanf2(fid,'%f');
hdrdat.HanningRadius   = fscanf2(fid,'%f');
hdrdat.nNoiseAC = fscanf2(fid,'%d');
hdrdat.BrainAirSeg    = fscanf2(fid,'%d');

%%--------- Gamma Fit -------------%%%
hdrdat.GammaFit   = fscanf2(fid,'%d');
if(hdrdat.GammaFit == 0)
  hdrdat.gfDelta = [];
  hdrdat.gfTau = [];
else
  dummy = fscanf(fid,'%s',1);
  hdrdat.gfDelta  = fscanf(fid,'%f',hdrdat.GammaFit);
  dummy = fscanf(fid,'%s',1);
  hdrdat.gfTau    = fscanf(fid,'%f',hdrdat.GammaFit);
end

hdrdat.NullCondId = fscanf2(fid,'%d');

if(~hdrdat.GammaFit)  Nch = hdrdat.Nnnc*hdrdat.Nh;
else                  Nch = hdrdat.Nnnc*length(hdrdat.gfDelta);
end

dummy      = fscanf(fid,'%s',1);
hdrdat.SumXtX     = fscanf(fid,'%f',Nch*Nch);

nSumXtX = prod(size(hdrdat.SumXtX));
if(nSumXtX ~= Nch*Nch )
  msg = sprintf('Inconsistent Dimensions: nSumXtX=%d, Nch^2 = %d\n',...
                nSumXtX,Nch*Nch);
  qoe(msg);error(msg);
end
hdrdat.SumXtX     = reshape(hdrdat.SumXtX,[Nch Nch]);

%--------- Version 2 -------------%
if(hdrdat.Version > 1)
  dummy          = fscanf(fid,'%s',1);
  [hdrdat.hCovMtx c] = fscanf(fid,'%f',Nch*Nch);
  if(c ~= Nch*Nch)
    msg = sprintf('Could not read hCovMtx from %s\n',datfile);
    hdrdat = [];
    qoe(msg); error(msg);
  end
  hdrdat.hCovMtx = reshape(hdrdat.hCovMtx,[Nch Nch]);

  dummy            = fscanf(fid,'%s',1);
  hdrdat.CondIdMap = fscanf(fid,'%f',hdrdat.Nc);
else
  hdrdat.hCovMtx   = hdrdat.SumXtX;
  hdrdat.CondIdMap = [0:hdrdat.Nc];
end

fclose(fid);

return

%%%%%%%%%%%%%%------------------------%%%%%%%%%%%%%%%%%%
%- reads a string followed by a value -----%
function x = fscanf2(fid,fmt)
  fscanf(fid,'%s',1);
  [x c] = fscanf(fid,fmt,1);
  if(c ~= 1)
    x = []; return;
  end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function N = tangentspacefactory(M, x)
% Returns a manifold structure representing the tangent space to M at x.
%
% N = tangentspacefactory(M, x)
%
% N defines a (linear) manifold that is the tangent space to M at x. Points
% are represented as tangent vectors to M at x. Tangent vectors are also
% represented as tangent vectors to M at x.
%
% This is chiefly useful to solve optimization problems involving tangent
% vectors to M at x, which notably comes up when solving linear systems
% involving, for example, the Hessian of the cost on M at x. The Riemannian
% (actually, Euclidean) structure on N is that of the tangent space to M,
% that is, the inner product is inherited.
%
% See also: preconhessiansolve

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 9, 2015.
% Contributors:
% Change log:

    % N is the manifold we build. y will be a point on N, thus also a
    % tangent vector to M at x. This is a typical Euclidean space, hence it
    % will be easy to describe in terms of the tools available for M.
    N = struct();

    % u, u1 and u2 will be tangent vectors to N at y. The tangent space to
    % N at y is the tangent space to M at x, thus u, u1 and u2 are also
    % tangent vectors to M at x.

    N.dim   = @() M.dim();
    N.inner = @(y, u1, u2) M.inner(x, u1, u2);
    N.norm  = @(y, u) M.norm(x, u);
    N.proj  = M.proj;
    N.typicaldist = @() N.dim();
    N.tangent = @(y, u) u;
    N.egrad2rgrad = @(x, g) g;
    N.ehess2rhess = @(x, eg, eh, d) eh;
    N.exp = @exponential;
    N.retr = @exponential;
    N.log = @(y1, y2) M.lincomb(x, 1, y2, -1, y1);
    N.pairmean = @(y1, y2) M.lincomb(x, 0.5, y1, 0.5, y2);
    N.rand = @() M.randvec(x);
    N.randvec = @(y) M.randvec(x);
    N.zerovec = M.zerovec;
    N.lincomb = M.lincomb;
    N.transp = @(y1, y2, u) u;
    N.hash = @(y) ['z' hashmd5(M.vec(x, y))];

    % In a Euclidean space, the exponential is merely the sum: y + tu.
    function yy = exponential(y, u, t)
        if nargin == 2
            t = 1;
        end
        yy = M.lincomb(x, 1, y, t, u);
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [D, Nrows, Ncols] = addimagesize(D, HOMEIMAGES)
%
% [D, nrows, ncols] = addimagesize(D, HOMEIMAGES)

n = length(D);
Nrows = zeros(n,1);
Ncols = zeros(n,1);

for i = 1:n
    if mod(i,10)==1; disp(n-i); end
    disp(fullfile(D(i).annotation.folder, D(i).annotation.filename))

    if ~isfield(D(i).annotation,'imagesize')
        try
            info = imfinfo(strrep(strrep(fullfile(HOMEIMAGES, D(i).annotation.folder, D(i).annotation.filename), '\', '/'), ' ', '%20'));
            nrows = info.Height;
            ncols = info.Width;
        catch
            [ncols, nrows] = getaproximagesize(D(i).annotation);
        end

        D(i).annotation.imagesize.nrows = (nrows);
        D(i).annotation.imagesize.ncols = (ncols);

    else
        if ischar(D(i).annotation.imagesize.nrows)
            D(i).annotation.imagesize.nrows = str2num(D(i).annotation.imagesize.nrows);
            D(i).annotation.imagesize.ncols = str2num(D(i).annotation.imagesize.ncols);
        end
    end

    Nrows(i) = D(i).annotation.imagesize.nrows;
    Ncols(i) = D(i).annotation.imagesize.ncols;
end

author: CSAILVision

repo: LabelMeToolbox


text:
check_active_images;

for kk =  ind_active,

   if ~exist(['x_' num2str(kk)]),

      fprintf(1,'WARNING: Need to extract grid corners on image %d\n',kk);

      active_images(kk) = 0;

      eval(['dX_' num2str(kk) ' = NaN;']);
      eval(['dY_' num2str(kk) ' = NaN;']);

      eval(['wintx_' num2str(kk) ' = NaN;']);
      eval(['winty_' num2str(kk) ' = NaN;']);

      eval(['x_' num2str(kk) ' = NaN*ones(2,1);']);
      eval(['X_' num2str(kk) ' = NaN*ones(3,1);']);

      eval(['n_sq_x_' num2str(kk) ' = NaN;']);
      eval(['n_sq_y_' num2str(kk) ' = NaN;']);

   else

      eval(['xkk = x_' num2str(kk) ';']);

      if isnan(xkk(1)),

	 fprintf(1,'WARNING: Need to extract grid corners on image %d - This image is now set inactive\n',kk);

	 active_images(kk) = 0;

      end;

   end;

end;

repo: ekfmonoslam


text:
% wls_grpr_test.m
% test WLS gradient projection method for various preconditionners
% Copyright Dec. 2000, Jeff Fessler, The University of Michigan

%
% generate data
%
if ~isvar('yi') || 1
%	xtrue = [1 2]';
%	a = -0.5;
%	H = [1 a; a 1];
%	b = [-1 1]';

	G = [1 -1; 0 1];
	G = [1 3; 1 1];
	W = 1;
	H = G' * W * G

	b = H * [-1 1]';
	yi = G' \ b;

%	yi = [-2 1]';
	b = G' * W * yi
end

%
% GP
%
if ~isvar('xgp') || 1
	f.niter = 18;
	xinit = [0.5 0]';
	xinit = [0.5 2]';

	D = 0.3 * diag([1 1]);
	D = diag([0.3 0.1]);
	D = 1.5 / norm(H);
	D = diag([1.3 0.3]) / norm(H);
%	D = inv(H);
	xmin = 0;
	xgp = wls_grpr(xinit, G, W, yi, D, xmin, inf, f.niter);
end

if 1
	xmin = H \ b;
	xcon = [0; b(2) / H(2,2)];

	n1 = 51;
	n2 = 53;
	x1 = linspace(-2, 1, n1);
	x2 = linspace(-1, 2, n2);
	[xx1, xx2] = ndgrid(x1, x2);

	x = [xx1(:) xx2(:)]';
	f = sum(x .* (H * x))/2 - b' * x;
	f = reshape(f, n1, n2);
	f = f';

	if im
		clf
		contour(x1, x2, sqrt(f-min(f(:))))
		hold on
		plot([0 0], minmax(x2), 'w-')
		plot(minmax(x1), [0 0], 'w-')
		plot(xmin(1), xmin(2), 'y*')
		plot(xcon(1), xcon(2), 'yx')
		hold off
		grid
	end
end

if im
	hold on
	plot(xgp(1,:), xgp(2,:), 'g-o')
	hold off
	xlabel x1, ylabel x2
end

% ir_savefig 'fig_wls_grpr'

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function Qd = calculateQd(tau,alpha,h)
%calculates Qd for noise simulation algorithm from paper DISCRETE SIMULATION OF POWER LAW NOISE
%equation: Qd = h / (2*(2pi)^alpha * tau^(alpha-1))

Qd = h / (2*(2*pi)^alpha * tau^(alpha-1)); %calculate Qd

author: Sable

repo: mcbench-benchmarks


text:
function im = imNameToNum(imname)
	im = str2num(imname(5:8));
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
% PutSegyTrace(segyid,tracedata,SegyTraceHeader,SegyHeader)
% Write a SegyTrace to a filehandle 'segyid'
%
% (C) 2001-2012, Thomas Mejer Hansen, thomas.mejer.hansen@gmail.com
%

%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function PutSegyTrace(segyid,tracedata,SegyTraceHeader,SegyHeader);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% WRITE SegyTraceHeader

  if ~exist('TraceStart')==1, TraceStart=ftell(segyid);end

  WRITETRACEHEADER=1;
  if WRITETRACEHEADER==1,
    fseek(segyid,0,'bof');
    fseek(segyid,TraceStart,'bof');
    fwrite(segyid,SegyTraceHeader.TraceSequenceLine,'int32');    % 0
    fwrite(segyid,SegyTraceHeader.TraceSequenceFile,'int32');    % 4
    fwrite(segyid,SegyTraceHeader.FieldRecord,'int32');          % 8
    fwrite(segyid,SegyTraceHeader.TraceNumber,'int32');          % 12
    fwrite(segyid,SegyTraceHeader.EnergySourcePoint,'int32');    % 16
    fwrite(segyid,SegyTraceHeader.cdp,'int32');                  % 20
    fwrite(segyid,SegyTraceHeader.cdpTrace,'int32');             % 24
    fwrite(segyid,SegyTraceHeader.TraceIdentifactionCode,'int16'); % 28
    fwrite(segyid,SegyTraceHeader.NSummedTraces,'int16'); % 30
    fwrite(segyid,SegyTraceHeader.NStackedTraces,'int16'); % 32
    fwrite(segyid,SegyTraceHeader.DataUse,'int16'); % 34
    fwrite(segyid,SegyTraceHeader.offset,'int32');             %36
    fwrite(segyid,SegyTraceHeader.ReceiverGroupElevation,'int32');             %40
    fwrite(segyid,SegyTraceHeader.SourceSurfaceElevation,'int32');             %44
    fwrite(segyid,SegyTraceHeader.SourceDepth,'int32');             %48
    fwrite(segyid,SegyTraceHeader.ReceiverDatumElevation,'int32');             %52
    fwrite(segyid,SegyTraceHeader.SourceDatumElevation,'int32');             %56
    fwrite(segyid,SegyTraceHeader.SourceWaterDepth,'int32');  %60
    fwrite(segyid,SegyTraceHeader.GroupWaterDepth,'int32');  %64
    fwrite(segyid,SegyTraceHeader.ElevationScalar,'int16');  %68
                                                             % Multiply/divide next number for following 4 values
    fwrite(segyid,SegyTraceHeader.SourceGroupScalar,'int16');  %70
    fwrite(segyid,SegyTraceHeader.SourceX,'int32');  %72
    fwrite(segyid,SegyTraceHeader.SourceY,'int32');  %76
    fwrite(segyid,SegyTraceHeader.GroupX,'int32');  %80
    fwrite(segyid,SegyTraceHeader.GroupY,'int32');  %84
    fwrite(segyid,SegyTraceHeader.CoordinateUnits,'int16');  %88
    fwrite(segyid,SegyTraceHeader.WeatheringVelocity,'int16');  %90
    fwrite(segyid,SegyTraceHeader.SubWeatheringVelocity,'int16');  %92
    fwrite(segyid,SegyTraceHeader.SourceUpholeTime,'int16');  %94
    fwrite(segyid,SegyTraceHeader.GroupUpholeTime,'int16');  %96
    fwrite(segyid,SegyTraceHeader.SourceStaticCorrection,'int16');  %98
    fwrite(segyid,SegyTraceHeader.GroupStaticCorrection,'int16');  %100
    fwrite(segyid,SegyTraceHeader.TotalStaticApplied,'int16');  %102
    fwrite(segyid,SegyTraceHeader.LagTimeA,'int16');  %104
    fwrite(segyid,SegyTraceHeader.LagTimeB,'int16');  %106
    fwrite(segyid,SegyTraceHeader.DelayRecordingTime,'int16');  %108
    fwrite(segyid,SegyTraceHeader.MuteTimeStart,'int16');  %110
    fwrite(segyid,SegyTraceHeader.MuteTimeEND,'int16');  %112
    fwrite(segyid,SegyTraceHeader.ns,'uint16');  %114
    fwrite(segyid,SegyTraceHeader.dt,'uint16');  %116
    fwrite(segyid,SegyTraceHeader.GainType,'int16');  %118
    fwrite(segyid,SegyTraceHeader.InstrumentGainConstant,'int16');  %120
    fwrite(segyid,SegyTraceHeader.InstrumentInitialGain,'int16');  %%122
    fwrite(segyid,SegyTraceHeader.Correlated,'int16');  %124
    fwrite(segyid,SegyTraceHeader.SweepFrequenceStart,'int16');  %126
    fwrite(segyid,SegyTraceHeader.SweepFrequenceEnd,'int16');  %128
    fwrite(segyid,SegyTraceHeader.SweepLength,'int16');  %130
    fwrite(segyid,SegyTraceHeader.SweepType,'int16');  %132
    fwrite(segyid,SegyTraceHeader.SweepTraceTaperLengthStart,'int16');  %134
    fwrite(segyid,SegyTraceHeader.SweepTraceTaperLengthEnd,'int16');  %136
    fwrite(segyid,SegyTraceHeader.TaperType,'int16');  %138
    fwrite(segyid,SegyTraceHeader.AliasFilterFrequency,'int16');  %140
    fwrite(segyid,SegyTraceHeader.AliasFilterSlope,'int16');  %142
    fwrite(segyid,SegyTraceHeader.NotchFilterFrequency,'int16');  %144
    fwrite(segyid,SegyTraceHeader.NotchFilterSlope,'int16');  %146
    fwrite(segyid,SegyTraceHeader.LowCutFrequency,'int16');  %148
    fwrite(segyid,SegyTraceHeader.HighCutFrequency,'int16');  %150
    fwrite(segyid,SegyTraceHeader.LowCutSlope,'int16');  %152
    fwrite(segyid,SegyTraceHeader.HighCutSlope,'int16');  %154
    fwrite(segyid,SegyTraceHeader.YearDataRecorded,'int16');  %156
    fwrite(segyid,SegyTraceHeader.DayOfYear,'int16');  %158
    fwrite(segyid,SegyTraceHeader.HourOfDay,'int16');  %160
    fwrite(segyid,SegyTraceHeader.MinuteOfHour,'int16');  %162
    fwrite(segyid,SegyTraceHeader.SecondOfMinute,'int16');  %164
    fwrite(segyid,SegyTraceHeader.TimeBaseCode,'int16');  %166
    fwrite(segyid,SegyTraceHeader.TraceWeightningFactor,'int16');  %170
    fwrite(segyid,SegyTraceHeader.GeophoneGroupNumberRoll1,'int16');  %172
    fwrite(segyid,SegyTraceHeader.GeophoneGroupNumberFirstTraceOrigField,'int16');  %174
    fwrite(segyid,SegyTraceHeader.GeophoneGroupNumberLastTraceOrigField,'int16');  %176
    fwrite(segyid,SegyTraceHeader.GapSize,'int16');  %178
    fwrite(segyid,SegyTraceHeader.OverTravel,'int16');  %178
    fwrite(segyid,SegyTraceHeader.cdpX,'int32');  %180
    fwrite(segyid,SegyTraceHeader.cdpY,'int32');  %184
    fwrite(segyid,SegyTraceHeader.Inline3D,'int32');  %188
    fwrite(segyid,SegyTraceHeader.Crossline3D,'int32');  %192
    fwrite(segyid,SegyTraceHeader.ShotPoint,'int32');  %196
    fwrite(segyid,SegyTraceHeader.ShotPointScalar,'int16');  %200
    fwrite(segyid,SegyTraceHeader.TraceValueMeasurementUnit,'int16');  %202
    fwrite(segyid,SegyTraceHeader.TransductionConstantMantissa,'int32');  %204
    fwrite(segyid,SegyTraceHeader.TransductionConstantPower,'int16'); %208
    fwrite(segyid,SegyTraceHeader.TransductionUnit,'int16');  %210
    fwrite(segyid,SegyTraceHeader.TraceIdentifier,'int16');  %212
    fwrite(segyid,SegyTraceHeader.ScalarTraceHeader,'int16');  %214
    fwrite(segyid,SegyTraceHeader.SourceType,'int16');  %216
    fwrite(segyid,SegyTraceHeader.SourceEnergyDirectionMantissa,'int32');  %218
    fwrite(segyid,SegyTraceHeader.SourceEnergyDirectionExponent,'int16');  %222
    fwrite(segyid,SegyTraceHeader.SourceMeasurementMantissa,'int32');  %224
    fwrite(segyid,SegyTraceHeader.SourceMeasurementExponent,'int16');  %228
    fwrite(segyid,SegyTraceHeader.SourceMeasurementUnit,'int16');  %230
                                                                   % WRITE UNASSIGNED CHARACTERS FOR THE REST
    fwrite(segyid,SegyTraceHeader.UnassignedInt1,'int32');  %232
    fwrite(segyid,SegyTraceHeader.UnassignedInt2,'int32');  %236

  end

  Revision=SegyHeader.SegyFormatRevisionNumber;
  if Revision>0, Revision=1; end
  Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).format;
  BPS=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).bps;

  SegymatVerbose([mfilename,' SegyRevision ',sprintf('%0.4g',Revision),', ',Format],3);
  %% WRITE TRACE DATA;
  if (strcmp(Format,'uint32')==1)|(strcmp(Format,'uint16')==1), % IBM FLOATING POINT
                                                                % CONVERT FROM FLOATING POINT
    SegymatVerbose([mfilename,'Converting from IBM, DataFormat :',SegyHeader.DataSampleFormat],2);
    tracedata=double(num2ibm(tracedata));
  end;
  if SegyHeader.FixedLengthTraceFlag
      fwrite(segyid,tracedata,Format);
  else
      fwrite(segyid,tracedata(1:SegyTraceHeader.ns),Format);
  end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function [D,Tc,res,dD,dTc,dRes,d2Phi] = distance(varargin)
% NP, JH 2006/11/25
%
% the distance function is phrased as
% D(y) = D(R,T,y) = phi(res(y))
%
% hence
% D   = phi(res(y))
% dD  = dPhi(res(y)) * dRes(y)
% d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider
%
% example SSD:
%   phi = hd/2 res'*res,  dPhi = hd * res', d2Phi = hd
%   res = T(y) - R,       dRes = dT

D     = [];
Tc    = [];
res   = [];
dD    = [];
dTc   = [];
dRes  = [];
d2Phi = [];
%
% -----------------------------------------------------------------------------
persistent PARA

if nargout>0 & nargin == 0,            % return PARA
    D  = PARA;
    return;
end;

% clear PARA
if isstr(varargin{1}) & strcmp(varargin{1},'clear'),
    PARA = [];
    return;
end;
% set PARA
if isstr(varargin{1}) & strcmp(varargin{1},'set'),
    PARA = setPARA(PARA,varargin{2:end});
    return;
end;
% -----------------------------------------------------------------------------

% -----------------------------------------------------------------------------
% do the work
% -----------------------------------------------------------------------------

Rc     = varargin{1};
TD     = varargin{2};
Omega  = varargin{3};
m      = varargin{4};
Y      = varargin{5};

doDerivative = (nargout > 3);

h   = Omega./m;
hd  = prod(h);
n   = prod(m);
%X   = getGrid(Omega,m);
%Rc  = interpolation(RD,Omega,X);
[Tc,dTc] = interpolation(TD,Omega,Y,doDerivative);

switch PARA.MODE,

    case 'SSD',
        % D   = phi(res(y))
        % dD  = dPhi(res(y)) * dRes(y)
        % d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider
        %
        % example SSD:
        %   phi = hd/2 res'*res,  dPhi = hd * res, d2Phi = hd
        %   res = T(y) - R,       dRes = dT

        res   = (Tc-Rc);
        D     = hd/2 * res' * res;

        if ~doDerivative, return; end;

        dRes  = dTc;
        dD    = hd * res' * dRes;
        d2Phi = hd;

    case 'SSD_W',
        % D   = phi(res(y))
        % dD  = dPhi(res(y)) * dRes(y)
        % d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider
        %
        % example SSD:
        %   phi = hd/2 res'*res,  dPhi = hd * res, d2Phi = hd
        %   res = T(y) - R,       dRes = dT

        mask   = varargin{6};

        res    = (Tc-Rc).*mask;
        D      = hd/2 * res' * res;

        if ~doDerivative, return; end;

        dRes  = sdiag(mask) * dTc;
        dD    = hd * res' * dRes;
        d2Phi = hd;

    case 'SSD_WP',
        % D   = phi(res(y))
        % dD  = dPhi(res(y)) * dRes(y)
        % d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider
        %
        % example SSD:
        %   phi = hd/2 res'*res,  dPhi = hd * res, d2Phi = hd
        %   res = T(y) - R,       dRes = dT
        period = PARA.period;
        period2 = period/2;

        mask   = varargin{6};

        res1   = Tc - Rc;

        p      = zeros(size(res1));
        p(res1 > period2)  = -period;
        p(res1 < -period2) =  period;

        res1   = res1 + p;
        res    = res1.*mask;
        D      = hd/2 * res' * res;

        if ~doDerivative, return; end;

        dRes  = sdiag(mask) * dTc;
        dD    = hd * res' * dRes;
        d2Phi = hd;

    case 'NGF',
        % D   = phi(res(y))
        % dD  = dPhi(res(y)) * dRes(y)
        % d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider
        %
        % example NGF:
        %   phi = hd/2 res'*res,  dPhi = hd * res, d2Phi = hd
        %   dR_i  = \nabla R_i / |\nabla R_i|_edge
        %   res = (\nabla T(y_i)' * dR_i) / ( |\nabla T(y_i)|_edge )
        %   dRes = complicated, see below

        edge = PARA.edge;
        [G1,G2,G3] = getGrad('c',Omega,m);
        if isempty(G3), G3 = 0; end;

        d1R = G1 * Rc(:);    d2R = G2 * Rc(:);    d3R = G3 * Rc(:);
        d1T = G1 * Tc(:);    d2T = G2 * Tc(:);    d3T = G3 * Tc(:);

        ndR = sqrt(d1R.^2 + d2R.^2 + d3R.^2 + edge^2);
        ndT = sqrt(d1T.^2 + d2T.^2 + d3T.^2 + edge^2);

        nd1R = d1R./ndR;      %nd1T = d1T./ndT;
        nd2R = d2R./ndR;      %nd2T = d2T./ndT;
        nd3R = d3R./ndR;      %nd3T = d3T./ndT;

        res1 = (nd1R.*d1T + nd2R.*d2T + nd3R.*d3T);
        res2 = 1./ndT;
        res  = res1 .* res2;

        D    = -hd/2 * res' * res;

        if ~doDerivative, return; end;

        dRes1 = sdiag(nd1R)*G1 + sdiag(nd2R)*G2 + sdiag(nd3R)*G3;
        dRes2 = -sdiag(1./ndT.^3)*(sdiag(d1T)*G1 + sdiag(d2T)*G2 + sdiag(d3T)*G3);
        dRes  = (sdiag(res2)*dRes1 + sdiag(res1)*dRes2) * dTc;

        dD    = -hd * res' * dRes;
        d2Phi = hd; % note the missing minus sign is not a bug!

    case 'MI'
        % D   = phi(res(y))
        % dD  = dPhi(res(y)) * dRes(y)
        % d2D = res(y)' * d2Phi(res(y)) * dRes(y) + stuff we don't consider
        %
        % example MI:
        %   phi   = res' * log(res + tol) + ...
        %   dPhi  = log(res + tol) + res./(res + tol) + ...
        %   d2Phi = (res + 2*tol)./(res + tol)^2 + ...
        %   res   = rho(T(y),R)
        %   dRes  = drho, see pdfestimate
        tol         = PARA.entropyTol;
        [rho,drho]  = pdfestimate(Rc,Tc,PARA,doDerivative);
        [n1,n2]     = size(rho);

        rhoR = sum(rho,2);
        rhoT = sum(rho,1)';
        rho  = rho(:);

        res  = rho;
        D    = rhoR'*log(rhoR+tol)+rhoT'*log(rhoT+tol) - rho'*log(rho+tol);

        if ~doDerivative, return; end;

        SR    = sparse(kron(ones(1,n2),speye(n1,n1)));
        ST    = sparse(kron(speye(n2,n2),ones(1,n1)));

        dPhi  = ...
            (log(rhoR+tol)+rhoR./(rhoR+tol))'*SR ...
            +(log(rhoT+tol)+rhoT./(rhoT+tol))'*ST ...
            -(log(rho +tol)+rho ./(rho +tol))';

        dRes  = drho*dTc;
        dD    = dPhi * dRes;

        d2Phi = ...
            SR'*sdiag((rhoR + 2*tol)./(rhoR+tol).^2)*SR ...
            +ST'*sdiag((rhoT + 2*tol)./(rhoT+tol).^2)*ST ...
            -sdiag((rho + 2*tol)./(rho+tol).^2);

        a = 1/sqrt(PARA.ngvR*PARA.ngvT);
        a = 1e0;
        d2Phi = - a * d2Phi;

    otherwise, error(sprintf('MODE = %s not implemented',MODE));
end;

return;

% =============================================================================
function d = sdiag(d);
d = spdiags(d,0,length(d),length(d));
% =============================================================================
function testMe(varargin)

fprintf('test %s\n',mfilename)

if nargin == 0,
    TD = zeros(100,100);
    TD(31:60,31:60) = 100;
    RD = TD + 1e1*randn(size(TD));
    Omega = [1,1];
    m     = [32,32];
    X  = getGrid(Omega,m);
    interpolation('set','MODE','linear-smooth');
    Rc = interpolation(RD,Omega,X);

    X = X+0.1;
end;

fctn = @(X) distance(Rc,TD,Omega,m,X);
testDerivative(fctn,X);

return;

% =============================================================================
function testDerivative(fctn,xc)

fprintf('derivative(%s)\n',mfilename)

[f,Tc,res,df]  = feval(fctn,xc);
v       = randn(size(xc));
dfv     = df*v;
hh      = logspace(0,-10,11);
fprintf('%12s %12s %12s \n','h','|fc-ft|','|fc+vdfc-ft|');
for j=1:length(hh),
    xt = xc + hh(j)*v;
    ft = feval(fctn,xt);
    n1 = norm(f(:)-ft(:));
    n2 = norm(f(:)+hh(j)*dfv-ft(:));
    fprintf('%12.4e %12.4e %12.4e\n',hh(j),n1,n2);
end;
return;
% =============================================================================

% ------------------------------------------------------------------------------
% set persistent parameter
% -----------------------------------------------------------------------------
function PARA = setPARA(PARA,varargin)

%disp(mfilename)
%varargin{:}

if ~isfield(PARA,'MODE') | isempty(getfield(PARA,'MODE')),
    PARA = setfield(PARA,'MODE','SSD');
end;

for j=1:length(varargin),
    if strcmp(varargin{j},'MODE'),
        PARA.MODE = varargin{j+1};
        varargin([j,j+1]) = [];
        break;
    end;
end;

fprintf('set distance MODE to %s, ',PARA.MODE);
switch PARA.MODE,
    case 'SSD',
        fprintf('no additional parameter needed\n');
    case 'SSD_W',
        fprintf('no additional parameter needed\n');
    case 'SSD_WP',
        if ~isfield(PARA,'period'), PARA.period = [];  end;
        if isempty(PARA.period),    PARA.period = 10;  end;
        for k=1:length(varargin)/2,
            str = sprintf('PARA=setfield(PARA,''%s'',varargin{%d});',...
                varargin{2*k-1},2*k);
            %disp(str);
            eval(str);
        end;
        fprintf('period=%s\n',num2str(PARA.period));
    case 'NGF',
        if ~isfield(PARA,'edge'), PARA.edge = [];  end;
        if isempty(PARA.edge),    PARA.edge = 10;  end;
        for k=1:length(varargin)/2,
            str = sprintf('PARA=setfield(PARA,''%s'',varargin{%d});',...
                varargin{2*k-1},2*k);
            %disp(str);
            eval(str);
        end;
        fprintf('edge=%s\n',num2str(PARA.edge));
    case 'MI',
        PARA = setMIpara(PARA,varargin{:});
        fprintf('setMIpara used\n');
    otherwise,
        error(sprintf('MODE = %s',MODE))
end;

return;
% -----------------------------------------------------------------------------
% -----------------------------------------------------------------------------


text:
function a = in( x, y )

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '==' );
if nargout, a = b; end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function cPlot(x,y,varargin)
mode = 3;  % 1 = magnitude and phase,
           % 2 = real and imag
           % 3 = real and imag in one window
doHolding = false;
plotLegend = 0;
deleteInds = zeros(length(varargin),1);
parent = [];
if nargin == 1
    y = x;
    x = 1 : length(y);
end
if nargin > 2
    for i = 1 : floor(length(varargin)/2)
        option       = varargin{i*2-1};
        option_value = varargin{i*2};

        switch lower(option)
            case 'hold'
                doHolding  = option_value;
                deleteInds(i*2-1 : i*2) = [1;1];
            case 'mode'
                mode       = option_value;
                deleteInds(i*2-1 : i*2) = [1;1];
            case 'parent'
                parent       = option_value;
                deleteInds(i*2-1 : i*2) = [1;1];
        end
    end
end
varargin(deleteInds ~= 0) = [];

if isempty(parent)
    parent = gca;
end

if doHolding
    hold on
end

switch mode
    case 1
        subplot(parent,2,1,1)
        plot(x,abs(y),varargin{:});
        title('abs');
        subplot(parent,2,1,2)
        plot(x,angle(y),varargin{:});
        title('phase');
    case 2
        subplot(parent,2,1,1)
        plot(x,real(y),varargin{:});
        title('real');
        subplot(parent,2,1,2)
        plot(x,imag(y),varargin{:});
        title('imag');
    case 3
        % put parent option at the end of varargin vector to allow for a
        % preceeding cursor style option (e.g. "x-");
        varargin = [varargin,{'parent'},{parent}];

        plot(x,[real(y),imag(y)],varargin{:});
        if plotLegend
            legend('real','imag');
        end

end
hold off
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function str = info(i, game)

data = get(game.Players(i), 'UserData');
str = {game.names{i}; ...
       ['Troops: ' num2str(data(1))]; ...
       ['Conquered: ' num2str(data(2))]};


text:
% ijcvLocationFigureScript

outdir = '/IUS/vmr20/dhoiem/data/ijcv06/results';

nb = 50; % nbins

vimages = ones(nb, nb, 3);
himages = ones(nb, nb, 5);

for f = 1:numel(imsegs)

    segimage = imsegs(f).segimage;
    im = im2double(imread([imdir '/' imsegs(f).imname]));
    [hue, sat, val] = rgb2hsv(im);
    hue = max(ceil(hue*nb),1);
    sat = max(ceil((1-sat)*nb),1);

    vlab = imsegs(f).vert_labels(imsegs(f).segimage);
    hlab = imsegs(f).horz_labels(imsegs(f).segimage);

    for k = 1:numel(segimage)
        if vlab(k)~=0
            vimages(sat(k), hue(k), vlab(k)) = vimages(sat(k), hue(k), vlab(k)) + 1;
        end
        if hlab(k)~=0
            himages(sat(k), hue(k), hlab(k)) = himages(sat(k), hue(k), hlab(k)) + 1;
        end
    end

    if mod(f, 50)==0
        disp(num2str(f))
    end
end

probcolor = sum(vimages, 3);
probcolor = probcolor / sum(probcolor(:));
probcolor = probcolor / max(probcolor(:));

vimages = vimages ./ repmat(sum(vimages, 3), [1 1 3]);
himages = himages ./ repmat(sum(himages, 3), [1 1 5]);

hueim = repmat([1:nb]/nb, [nb 1]);
satim = repmat(1-[1:nb]'/nb, [1 nb]);

imwrite(hsv2rgb(imresize(cat(3, hueim, satim, probcolor),1)), [outdir '/color_prob.jpg'], 'Quality', 100);
imwrite(hsv2rgb(imresize(cat(3, hueim, satim, ones(nb,nb)),1)), [outdir '/color_full.jpg'], 'Quality', 100);
for v = 1:3
    imwrite(hsv2rgb(imresize(cat(3, hueim, satim, vimages(:, :, v)),1)), ...
        [outdir '/vcol' num2str(v) '.jpg'], 'Quality', 100);
end

for h = 1:5
    imwrite(hsv2rgb(imresize(cat(3, hueim, satim, himages(:, :, h)),1)), ...
        [outdir '/hcol' num2str(h) '.jpg'], 'Quality', 100);
end

% [tmp, bestv] = max(vimages, [], 3);
% [tmp, besth] = max(himages, [], 3);
% colimseg.nseg = 7;
% colimseg.segimage = (bestv==1) + (bestv==2).*(1+besth) + 7*(bestv==3);
% for k = 1:7
%     colimseg.npixels(k) = sum(colimseg.segimage(:)==k);
% end
% lim = APPgetLabeledImage(ones(200,200, 3), colimseg, ...
%     {'000', '090', '090', '090', '090', '090', 'sky'}, ones(7, 1), ...
%     {'---', '045', '090', '135', 'por', 'sol', '---'}, ones(7, 1));
% imwrite(lim, [outdir '/labeledcol.jpg'], 'Quality', 100);
% lim2 = (lim>0).*vimages(:, :, [2 1 3]);
% imwrite(lim2, [outdir '/labeledcol.jpg'], 'Quality', 100);

author: Cloud-CV

repo: object-proposals


text:
function PlotStates(timeAll, stateAll, P)

clf;

    %break out the data:
        time = timeAll;
        xPos = stateAll(1,:);
        yPos = stateAll(2,:);
        xVel = stateAll(3,:);
        yVel = stateAll(4,:);

    subplot(2,2,1)
        plot(time,xPos,'LineWidth',P.CurveLineWidth);
        xlabel('time (s)','FontSize',P.LabelFontSize)
        ylabel('horizontal position (m)','FontSize',P.LabelFontSize)
        set(gca,'fontsize',P.AxisFontSize);
    subplot(2,2,2)
        plot(time,yPos,'LineWidth',P.CurveLineWidth);
        xlabel('time (s)','FontSize',P.LabelFontSize)
        ylabel('vertical position (m)','FontSize',P.LabelFontSize)
        set(gca,'fontsize',P.AxisFontSize);
    subplot(2,2,3)
        plot(time,xVel,'LineWidth',P.CurveLineWidth);
        xlabel('time (s)','FontSize',P.LabelFontSize)
        ylabel('horizontal speed (m/s)','FontSize',P.LabelFontSize)
        set(gca,'fontsize',P.AxisFontSize);
    subplot(2,2,4)
        plot(time,yVel,'LineWidth',P.CurveLineWidth);
        xlabel('time (s)','FontSize',P.LabelFontSize)
        ylabel('vertical speed (m/s)','FontSize',P.LabelFontSize)
        set(gca,'fontsize',P.AxisFontSize);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% D training using ImageNet model and RPN proposals
% --------------------------------------------------------
% D&T implementation
% Modified from MATLAB  R-FCN (https://github.com/daijifeng001/R-FCN/)
% and Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2017, Christoph Feichtenhofer
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------
clc;
clear mex;
clear is_valid_handle;
run(fullfile(fileparts(fileparts(mfilename('fullpath'))), 'startup'));

%% -------------------- CONFIG --------------------
opts.caffe_version          = 'caffe-rfcn';
opts.gpu_id                 = auto_select_gpu;

active_caffe_mex(opts.gpu_id, opts.caffe_version);

res50 = 0 ;
with_ohem = true;
concat = false;
modality = 'rgb';
with_resNext = 1;
root_path = get_root_path;

if res50
  model                       = Model.ResNet50_for_RFCN_ILSVRCvid_dilation();
  opts.cache_name             = 'rfcn_ILSVRC_ResNet50_OHEM_rpn_dilation';
  conf                        = rfcn_config_ohem('image_means', model.mean_image);
else
  model                       = Model.ResNet101_for_RFCN_ILSVRCvid_dilation();
  opts.cache_name             = 'rfcn_ILSVRC_ResNet101_OHEM_rpn_dilation';
  conf                        = rfcn_config_ohem('image_means', model.mean_image);
  if with_resNext
    model                       = Model.ResNetX101_for_RFCN_ILSVRCvid_dilation();
    opts.cache_name             = 'rfcn_ILSVRC_ResNeXt101_OHEM_rpn_dilation';
    conf                        = rfcn_config_ohem('image_means', model.mean_image);
  end
end

% check base model existance
check_dl_model(model.net_file);

conf.root_path = root_path ;
fprintf('Loading dataset...')
dataset = [];
dataset = Dataset.ilsvrc15(dataset, 'VIDval16lite3frames', false, root_path);
dataset = Dataset.ilsvrc15(dataset, 'VID+DETtrain16sample', false, root_path);
% do validation after val_interval iters, or not
opts.do_val                 = true;

%% -------------------- TRAINING --------------------
caffe.reset_all();
conf.input_modality = 'rgb';
conf.ims_per_batch = 2 ;
conf.max_size = [1000 ];

conf.sample_vid = false;
conf.use_flipped = true;
opts.cache_name_  = [opts.cache_name '_' conf.input_modality '300RPN_batch4_1e3LR_160kiter' ];

conf.num_proposals = 300;

% model.net_file = fullfile(root_path, 'output', 'rfcn_cachedir', opts.cache_name_, imdbs_name , ...
%   'iter_140000');

opts.rfcn_model        = rfcn_train(conf, dataset.imdb_train, dataset.roidb_train, ...
                                'offline_roidb',           true, ...
                                'output_dir',     root_path, ...
                                'do_val',           opts.do_val, ...
                                'imdb_val',         dataset.imdb_test, ...
                                'roidb_val',        dataset.roidb_test, ...
                                'solver_def_file',  model.solver_def_file, ...
                                'net_file',         model.net_file, ...
                                'cache_name',       opts.cache_name_, ...
                                'caffe_version',    opts.caffe_version, ...
                                'resume_iter', 0, ...
                                'val_iters', 500, ...
                                'val_interval', 5000, ...
                                'visualize_interval', 0);

%% -------------------- TESTING --------------------
%%
caffe.reset_all();
imdbs_name = cell2mat(cellfun(@(x) x.name, dataset.imdb_train, 'UniformOutput', false));
opts.rfcn_model = fullfile(root_path, 'output', 'rfcn_cachedir', opts.cache_name_, imdbs_name , ...
  'iter_140000');
proposal_path = fullfile(root_path, 'Data/VID/val/RPN_proposals/');

conf.root_path = root_path ;
imdb_val_lite = Dataset.ilsvrc15([], 'VIDval16', false, root_path);
% imdb_val_lite = Dataset.ilsvrc15([], 'VIDval16liteStride10', false, root_path);
% imdb_val_lite = Dataset.ilsvrc15([], 'VIDval16lite3frames', false, root_path);
test_classes = [];

conf.num_proposals = 300;
rfcn_test(conf, imdb_val_lite.imdb_test, imdb_val_lite.roidb_test, ...
      'net_def_file',     model.test_net_def_file, ...
      'net_file',         opts.rfcn_model, ...
      'cache_name',       opts.cache_name_,...
      'ignore_cache',     false, ...
      'proposal_path', proposal_path, ...
      'visualize_every', 0, ...
      'test_classes', test_classes, ...
      'suffix', '');

%% proposal test
recall_per_cls = compute_recall_ilsvrc15(proposal_path, 300, imdb_val_lite.imdb_test);
mean_recall = mean([recall_per_cls.recall]);
fprintf('mean rec:: %.2f\n\n', 100*mean_recall);

author: feichtenhofer

repo: Detect-Track


text:
function uvTformDataT = sr_src_domain_tform(uvPlaneID, modelPlane, uvTformAData, srcPos, trgPos)
% SR_SRC_DOMAIN_TFORM
%
% The function estimates the source patch domain transformation based on
% the source and target position, the plane parameters and apply the affine
% transformation
%
% Input:
%     - uvPlaneID:    [numUvPix] x [1], the plane label of uv pixels
%     - modelPlane:   plane model
%     - uvTformAData: [numUvPix] x 4, affine transformation
%     - srcPos:       [numUvPix] x [2], the center position of source patch
%     - trgPos:       [numUvPix] x [2], the center position of target patch
% Output:
%     - uvTformDataT: [numUvPix] X [9]
%
% The function of computing homography transformation induced by plane
% is modified from the following paper:
%
% Jia-Bin Huang, Sing Bing Kang, Narendra Ahuja, and Johannes Kopf,
% Image Completion using Planar Structure Guidance,
% ACM Transactions on Graphics (Proceedings of SIGGRAPH 2014), 33(4), 2014
% =========================================================================

% TO-DO: Make sure that the predicted transformation can point it
% the desired source patch position

% =========================================================================
% Compute homography transformation induced by plane
% =========================================================================

numUvPix = size(srcPos, 1);
uvTformDataH = zeros(numUvPix, 9, 'single');
I = eye(3);

for indPlane = 1: modelPlane.numPlane
    % The rectifying transformation for the plane
    rectMat = modelPlane.rectMat{indPlane};
    h7 = rectMat(3,1);
    h8 = rectMat(3,2);

    % Retrieve the uv pixels that have the current plane label
    uvPlaneIndCur = uvPlaneID == indPlane;
    numPlanePixCur = sum(uvPlaneIndCur);

    if(numPlanePixCur)
        % Target patch center position in the rectified domain
        trgPosCur = trgPos(uvPlaneIndCur, :) - 1;
        trgPosCurR = sr_apply_tform_H(trgPosCur, h7, h8);

        % Source patch center position in the rectified domain
        srcPosCur = srcPos(uvPlaneIndCur, :) - 1;
        srcPosCurR = sr_apply_tform_H(srcPosCur, h7, h8);

        % Displacement vector from target to source position
        dRect = srcPosCurR - trgPosCurR;

        % Compute the transformation that maps from target to source
        % (See Eqn 8 in the paper [Huang et al. TOG 2014])
        uvTformCur = zeros(numPlanePixCur, 9, 'single');
        uvTformCur(:,[1,4,7]) = bsxfun(@times, dRect(:,1), [h7, h8, 1]);
        uvTformCur(:,[2,5,8]) = bsxfun(@times, dRect(:,2), [h7, h8, 1]);

        dTemp = dRect*[h7;h8]; % dTemp = dx*h7 + dy*h8
        uvTformCur(:,[3,6,9]) = bsxfun(@times, dTemp, -[h7, h8, 1]);
        uvTformCur = bsxfun(@plus, uvTformCur, I(:)');

        % Apply the offset to cancel out the dependency of the target position
        % (See Eqn 9 in the paper [Huang et al. TOG 2014])
        uvTformDataH(uvPlaneIndCur, :)   = sr_trans_tform(uvTformCur, trgPosCur);
        uvTformDataH(uvPlaneIndCur, 7:8) = uvTformDataH(uvPlaneIndCur, 7:8) + 1;
    end
end

% =========================================================================
% Apply the similarity transformation
% =========================================================================
% The transformation T first map from original reference points to affine
% transformation coordinate using A, and then map to the desired target
% position using H
% -> T = H * A
uvTformDataT = uvTformDataH;
uvTformDataT(:,1:3) = bsxfun(@times, uvTformDataH(:,1:3), uvTformAData(:,1)) + ...
    bsxfun(@times, uvTformDataH(:,4:6), uvTformAData(:,2));
uvTformDataT(:,4:6) = bsxfun(@times, uvTformDataH(:,4:6), uvTformAData(:,3)) + ...
    bsxfun(@times, uvTformDataH(:,4:6), uvTformAData(:,4));

uvTformDataT = bsxfun(@rdivide, uvTformDataT, uvTformDataT(:,9));

end

function x = sr_apply_tform_H(x, h7, h8)

% Apply homography H with third row [h7, h8, 1] to 2D points x

y = x(:,1)*h7 + x(:,2)*h8 + 1;
x = bsxfun(@rdivide, x(:,1:2), y + eps);

% numPix = size(x, 1);
% x = cat(2, x, ones(numPix, 1, 'single'));
% x = x*H';
% x = bsxfun(@rdivide, x(:,1:2), x(:,3) + eps);

end

repo: SelfExSR


text:
function btv_test04 ( )

%*****************************************************************************80
%
%% BTV_TEST04 tests BURGERS_TIME_VISCOUS with the shock initial condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTV_TEST04\n' );
  fprintf ( 1, '  Test BURGERS_TIME_VISCOUS with the shock initial condition.\n' );
  fprintf ( 1, '  Use periodic boundaries.\n' );

  nx = 81;
  nt = 300;
  t_max = 3.0;
  nu = 0.01;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial condition: shock\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Viscosity = %g\n', nu );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_viscous ( @ic_shock, nx, nt, t_max, nu, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 4 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers equation solutions over time, initial condition shock' )

  filename = 'btv_test04.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% WriteSu : writes seismic data in the Seismic Unix format
%
% EX
% WriteSu('datacube.su',data,'dt',.004,'Inline3D',Inline,'Crossline3D',Crossline,'cdpX',X,'cdpY',Y);
%
% to use a specific SEG revision use :
%
% to use a specific Data Sampling Format use :
% WriteSu('test.su',seisdata,'dsf',1); % IBM FLAOTING POINT
%
% WriteSu('test.su',seisdata,'dsf');
%

%
% (C) 2001-2004, Thomas Mejer Hansen, thomas@cultpenguin.com
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function WriteSu(filename,data,varargin);

[ns,ntraces]=size(data);

if exist('dt')~=1
  dt=0.004;
  SegymatVerbose([mfilename,' - No dt set. Using dt=',num2str(dt)])
end

for i=1:2:length(varargin)
  var=varargin{i};
  val=varargin{i+1};
  eval([var,'=[',num2str(val),'];']);
end

SegyHeader.SegyFormatRevisionNumber=256;
SegyHeader.DataSampleFormat=5;
SegyHeader.Rev=GetSegyHeaderBasics;

% OPEN SEGY FILE HANDLE
segyid = fopen(filename,'w');
% segyid = fopen(filename,'w','b'); % BIG ENDIAN

SegyTraceHeader=InitSegyTraceHeader(ns,dt*1e+6);
for i=1:ntraces;
    if (i/100)==round(i/100),SegymatVerbose(['writing trace ',num2str(i),' of ',num2str(ntraces)],1),end
    % Basic TraceHeader information
    %SegyTraceHeader.ns=ns;
    %SegyTraceHeader.dt=dt.*1e+6;
    % Update TraceHeader information if available
    if exist('cdpX')==1,SegyTraceHeader.cdpX = cdpX(i);end
    if exist('offset')==1,SegyTraceHeader.offset = offset(i);end
    if exist('cdpY')==1,SegyTraceHeader.cdpY = cdpY(i);end
    if exist('Inline3D')==1,SegyTraceHeader.Inline3D = Inline3D(i);end
    if exist('Crossline3D')==1,SegyTraceHeader.Crossline3D=Crossline3D(i);end
    % Write the Trace
    PutSegyTrace(segyid,data(:,i),SegyTraceHeader,SegyHeader);
end

% CLOSE SEGY FILE HANDLE

fclose(segyid);


text:
function bezier_surface_node_print ( node_num, node_xyz )

%*****************************************************************************80
%
%% BEZIER_SURFACE_NODE_PRINT prints nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the coordinates of the
%    nodes.
%
  dim_num = 3;

  r8mat_transpose_print ( dim_num, node_num, node_xyz, ...
    '  Bezier Surface Nodes:' );

  return
end


text:
function checkxy_N(x, y, function_name, x_var_name, y_var_name, x_pos, y_pos)
%CHECKXY_N Check validity of map x and y vectors
%
%   CHECKXY(X, Y, FUNCTION_NAME, X_VAR_NAME, Y_VAR_NAME, X_POS, X_POS)
%   ensures that X and Y are real vectors of matching size and equal NaN
%   locations.

% Copyright 2006-2011 The MathWorks, Inc.
% $Revision: 1.1.6.3 $  $Date: 2011/05/17 02:12:40 $

% Input arguments are not checked for validity.

if ~isempty(x) || ~isempty(y)
    % Check numeric, 2d, and real.
    validateattributes(x, ...
        {'numeric'}, {'real','2d','vector'}, function_name, x_var_name, x_pos);
    validateattributes(y, ...
        {'numeric'}, {'real','2d','vector'}, function_name, y_var_name, y_pos);

    if ~isequal(isnan(x), isnan(y))
        error(sprintf('map:%s:inconsistentXY', function_name), ...
            'Function %s expected its %s and %s input arguments, %s and %s, to match in size or NaN locations.', ...
            upper(function_name), num2ordinal(x_pos), num2ordinal(y_pos), ...
            x_var_name, y_var_name)
    end
end

author: zouchuhang

repo: LayoutNet


text:
function [kWh] = cal2kWh(cal)
% Convert energy or work from calories to kilowatt-hours.
% Note: these calories are different from the capitalized Calories found on
% American cereal boxes.  I suppose we'd be afraid to eat a candy bar with
% 250,000 calories in it.  Because then we'd feel like fatties.
% Chad A. Greene 2012
kWh = cal*0.000001163;


text:
RadarName = 'staddon_clutter';
RadarCoords.lat = 50.346069;
RadarCoords.lon = -4.113670;

RadarCoords_lat = 50.339026;
RadarCoords_lon = -4.126307;

% Load dataset
load(strcat(RadarName,'.mat'));
filename = 'smc_phd_clutter_fakellr.mat';
load(filename);
% load('smc_phd_no_clutter.mat');
% load('smc_phd_clutter_fakellr.mat');

% lon_lim = [-4.2637   -4.1024];
% lat_lim = [50.2699   50.3853];

lon_lim = [-4.176846509809199,-4.103030078951917];
lat_lim = [50.297532151330770,50.337022151330770];
% radar_meas_convert;

% Map plot
figure('units','normalized','outerposition',[0 0 .5 0.8])
ax(1) = gca;
plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
plots = [];
hold on;

N = size(MeasurementScans,2); % Simulation length

Tracks = [TrackList, DeletedTracks];

%% PLot detections
% for k=2:N
%
%     MeasurementList = MeasurementScans(k);
%     if(MeasurementList.NumMeasurements>0)
%         % Convert measurements to LLA and plot them
%     %     [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
%     %                                MeasurementList.Vectors(1,:),...
%     %                                0,...
%     %                                RadarCoords.lat,...
%     %                                RadarCoords.lon,...
%     %                                0,...
%     %                                referenceEllipsoid('wgs84'));
%         latlon = [MeasurementList.Metadata.LatLon];
%         lat = latlon(1,:);
%         lon = latlon(2,:);
%         plots(end+1) = plot(ax(1), lon,lat,'y.','MarkerSize', 4);
%     end
% end

%% Plot all existing tracks

for j=1:numel(Tracks)

    track = Tracks{j};
    if(strcmp(filename, 'smc_phd_clutter_fakellr_0.mat'))
        t = [183531,328072];
        if(find(t==track.Tag.ID,1))
            continue
        end
    else
    end
    % Convert track trajectory to LLA and plot it
    means = [track.Trajectory.Mean];
    [lat,lon,~] = ned2geodetic(means(3,:),...
                               means(1,:),...
                               0,...
                               RadarCoords.lat,...
                               RadarCoords.lon,...
                               0,...
                               referenceEllipsoid('wgs84'));
    traj_length = size(lon,2);

    plots(end+1) = plot(ax(1), lon(:,1),lat(:,1),'go','MarkerSize',15,'LineWidth',2);
    t = [766616, 128226, 592651, 237840];
    if (find(t == track.Tag.ID,1))
        plots(end+1) = plot(ax(1), lon, lat,'-','LineWidth',4,'Color', [0.9856, 0.7372, 0.2537]);
        plots(end+1) = plot(ax(1), lon, lat,'-.w','LineWidth',1);
    else
        plots(end+1) = plot(ax(1), lon, lat,'-.w','LineWidth',2);
    end
%     plots(end+1) = plot(ax(1), lon, lat,'-.w','LineWidth',2);
    plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'rx','MarkerSize',15,'LineWidth',2);
%     text(ax(1), lon(:,end), lat(:,end), num2str(track.Tag.ID),'FontSize',18,'Color','g');
end

% Plot extra track
means = [-4.12069846817424,50.3000236150299;-4.12101955629663,50.2999210786118;-4.12152794582374,50.2999381680148;-4.12192930597672,50.2999381680148;-4.12227715144264,50.2999210786118;-4.12259823956503,50.2999039892088;-4.12310662909214,50.2999210786118;-4.12348123190159,50.2999381680148;-4.12484585642174,50.3000065256269;-4.12588939281949,50.3000919726419;-4.12698644390432,50.3000748832389;-4.12741456140083,50.3000748832389;-4.12768213483616,50.3000236150299;-4.12800322295854,50.3000065256269;-4.12811025233267,50.3000065256269;-4.12843134045506,50.3000236150299;-4.12859188451625,50.3000236150299;-4.12979596497520,50.2999894362238;-4.13033111184585,50.2999552574178;-4.13078598668589,50.2999381680148;-4.13094653074709,50.2999723468208;-4.13142816293067,50.3000065256269;-4.13150843496126,50.2999894362238;-4.13169573636599,50.2999552574178;-4.13207033917544,50.2999723468208;-4.13241818464136,50.2999894362238;-4.13273927276375,50.2999723468208;-4.13289981682494,50.2999552574178;-4.13314063291673,50.2999039892088;-4.13343496369559,50.2999381680148;-4.13383632384857,50.2999894362238;-4.13410389728389,50.2999894362238;-4.13455877212394,50.3000065256269;-4.13466580149807,50.2999723468208;-4.13501364696399,50.2999723468208;-4.13525446305578,50.2999381680148;-4.13576285258289,50.2998868998058;-4.13616421273587,50.2999723468208;-4.13664584491946,50.3000236150299;-4.13702044772891,50.3000065256269;-4.13768938131721,50.3000065256269;-4.13833155756199,50.3000748832389;-4.13883994708910,50.3000748832389;-4.13913427786795,50.3000748832389;-4.13926806458562,50.3001090620449;-4.13942860864681,50.3000919726419;-4.13996375551745,50.3000919726419;-4.14028484363984,50.3000919726419;-4.14060593176223,50.3001432408509;-4.14098053457168,50.3000748832389;-4.14138189472466,50.3000577938359;-4.14178325487764,50.3000577938359;-4.14253246049655,50.3000236150299;-4.14376329829903,50.3000065256269;-4.14411114376495,50.2998868998058;-4.14429844516968,50.2998698104028;-4.14453926126147,50.2998014527907;-4.14496737875798,50.2998185421937;-4.14520819484977,50.2998868998058;-4.14582361375101,50.2999381680148;-4.14595740046867,50.3000065256269;-4.14622497390400,50.2999552574178;-4.14646578999579,50.2999381680148;-4.14689390749230,50.2999381680148;-4.14726851030175,50.2999210786118;-4.14745581170648,50.2999210786118;-4.14780365717240,50.2999381680148;-4.14804447326419,50.2999723468208;-4.14823177466892,50.3000065256269;-4.14857962013483,50.2999723468208;-4.14879367888309,50.2999723468208;-4.14906125231841,50.2999894362238;-4.14940909778433,50.3000407044329;-4.14964991387612,50.3000407044329;-4.14991748731145,50.2999723468208;-4.15055966355622,50.3001603302539;-4.15104129573980,50.3000748832389;-4.15144265589278,50.3002628666720;-4.15200456010696,50.3004166712991;-4.15283403775646,50.3003654030901;-4.15334242728357,50.3003654030901;-4.15371703009302,50.3003824924931;-4.15414514758954,50.3004337607021;-4.15438596368133,50.3004337607021]';
plots(end+1) = plot(ax(1), means(1,1),means(2,1),'go','MarkerSize',15,'LineWidth',2);
plots(end+1) = plot(ax(1), means(1,:), means(2,:),'-','LineWidth',4,'Color', [0.9856, 0.7372, 0.2537]);
plots(end+1) = plot(ax(1), means(1,:), means(2,:),'-.w','LineWidth',1);
plots(end+1) = plot(ax(1), means(1,end),means(2,end),'rx','MarkerSize',15,'LineWidth',2);
means = [-4.15651854566730,50.3007257705160;-4.15696081580833,50.3007386100374;-4.15728246681998,50.3007771286017;-4.15760411783164,50.3007642890803;-4.15776494333747,50.3008156471659;-4.15846855492547,50.3008926842945;-4.15877010274890,50.3009183633373;-4.15921237288993,50.3009312028587;-4.15963453984273,50.3009697214230;-4.15999639723084,50.3009825609444;-4.16029794505427,50.3010467585515;-4.16045877056010,50.3010724375943;-4.16074021519530,50.3010595980729;-4.16132320765393,50.3010467585515;-4.16146392997153,50.3010724375943;-4.16232836706536,50.3011366352014;-4.16279074039461,50.3011879932871;-4.16323301053564,50.3012393513728;-4.16345414560616,50.3012136723300;-4.16453971777050,50.3013677465870;-4.16530363892319,50.3013677465870;-4.16578611544067,50.3015346603655;-4.16606756007587,50.3015860184511]';
plots(end+1) = plot(ax(1), means(1,1),means(2,1),'go','MarkerSize',15,'LineWidth',2);
plots(end+1) = plot(ax(1), means(1,:), means(2,:),'-','LineWidth',4,'Color', [0.9856, 0.7372, 0.2537]);
plots(end+1) = plot(ax(1), means(1,:), means(2,:),'-.w','LineWidth',1);
plots(end+1) = plot(ax(1), means(1,end),means(2,end),'rx','MarkerSize',15,'LineWidth',2);

%% Plot region of low PD
x = [-3422.85261310741,-3191.58441229017,-2993.55608122595,-3228.08558459284, -3422.85261310741];
y = [1.472966334316646e+03,1.123217741935625e+03,1.243952901078573e+03,1.579667558371468e+03, 1.472966334316646e+03];
[y,x,~] = ned2geodetic(y,...
                       x,...
                       0,...
                       RadarCoords.lat,...
                       RadarCoords.lon,...
                       0,...
                       referenceEllipsoid('wgs84'));
plots(end+1) = plot(ax(1), x, y, 'm-', 'LineWidth', 2);

%% Plot region of high clutter
x = [-469.710602736631,-2572.92295458513,-2516.31234127506,-1320.74163107899,-835.307363807233, -469.710602736631];
y = [-4915.41679013513,-4065.50122858976,-2858.60918927224,-1380.19032924249,-1261.08679763585, -4915.41679013513];
%         x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
%         y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
[y,x,~] = ned2geodetic(y,...
                       x,...
                       0,...
                       RadarCoords.lat,...
                       RadarCoords.lon,...
                       0,...
                       referenceEllipsoid('wgs84'));
plots(end+1) = plot(ax(1), x, y, 'r-', 'LineWidth', 2);

plot(ax(1), RadarCoords_lon,RadarCoords_lat,...
                               '-s','MarkerSize',10,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);
%% Draw fragmented tracks
pos = [-4.166866227065740,50.308021682580770,0.003389529988345,0.002102476851853];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
pos = [-4.162648145302467,50.302582666377070,0.002937592656566,0.001919652777772];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
pos = [-4.157676834652895,50.299840305265950,0.004293404651904,0.001508298611114];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'c');
%% PLot false tracks
pos = [-4.129054136973541,50.326578326099290,0.002937592656566,0.001873946759261];
rectangle(ax(1), 'Position', pos, 'LineWidth', 2, 'EdgeColor', 'm');
%% Legend
h1 = plot(NaN, NaN, 'o','MarkerSize', 5, 'MarkerFaceColor', 'y', 'MarkerEdgeColor','k');
h2 = plot(NaN, NaN, 's','MarkerSize',10,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);
h3 = plot(NaN, NaN, 'ms', 'MarkerSize',15,'LineWidth', 2);
h4 = plot(NaN, NaN, 'r-', 'LineWidth', 2);
h5 = plot(NaN, NaN, 'go','MarkerSize',15,'LineWidth',2);
h6 = plot(NaN, NaN, 'rx','MarkerSize',15,'LineWidth',2);
h7 = plot(NaN, NaN, 'cs','MarkerSize',15,'LineWidth',2);
h8 = plot(NaN, NaN, 'b-','MarkerSize',15,'LineWidth',4);
h9 = plot(NaN, NaN, '-','LineWidth',4,'Color', [0.9856, 0.7372, 0.2537]);
legend([h4, h5, h6, h3, h7, h9], 'High-clutter Region', 'Track Birth', 'Track Death','False Tracks', 'Fragmentation Examples', 'Bonus Tracks','Location','northeast');
% legend([h1, h3, h5, h6], 'Detections', 'Obscured Region', 'Track Birth', 'Track Death','Location','southeast');
% legend([h1, h2, h3, h4], 'Detections', 'Radar Position', 'Obscured Region', 'High-clutter region');
% legend([h1, h3, h5, h6, h7, h9], 'Detections', 'Obscured Region', 'Track Birth', 'Track Death', 'Fragmentation Examples', 'Bonus Tracks','Location','southeast');
% legend([h1, h2, h3, h4], 'Detections', 'Radar Position', 'Obscured Region', 'High-clutter region');
% xlabel('Longitude');
% ylabel('Latitude');
set(ax(1),'YTickLabel',[]);
set(ax(1),'XTickLabel',[]);

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');
% Number of measurements
M = 64;
% Signal space
N = 2 * M;
% Dirac-Fourier Dictionary
Phi = spx.dict.simple.dirac_fourier_dict(M);

% Number of signals
Ss = [1, 2, 4, 8, 16, 32];
% Sparsity levels
Ks = 2:2:50;
num_trials = 100;
num_ks = length(Ks);
num_ss = length(Ss);
bp_success_with_k = zeros(num_ss, num_ks);

snr_threshold = 100;

for ns=1:num_ss
    % Current sparsity level
    S = Ss(ns);
    for nk=1:num_ks
        K = Ks(nk);
        num_bp_successes = 0;
        for nt=1:num_trials
            X = model_3_data(N, K, S);
            % Measurement vectors
            Y = Phi * X;

            % BP MMV solver instance
            bp_solver = spx.pursuit.joint.BasisPursuit(Phi);
            % Solve the sparse recovery problem
            result = bp_solver.solve_l2_l1_complex(Y);
            % Solution vectors
            X_BP = result.Z;
            % Comparison
            cs = spx.commons.SparseSignalsComparison(X, X_BP, K);
            snr = cs.cum_signal_to_noise_ratio;
            bp_success = snr > snr_threshold;
            num_bp_successes = num_bp_successes + bp_success;
            fprintf('S: %d, K=%d, trial=%d, residual bp: %e, SNR: %f dB\n'...
                , S, K, nt, cs.cum_difference_norm, snr);
        end
        bp_success_with_k(ns, nk) = num_bp_successes / num_trials;
    end
end

save ('bin/figure_2_dirac_fourier_dict_model_3_bp_success_with_k.mat');

author: indigits

repo: sparse-plex


text:
function sphere = createSphere(varargin)
%CREATESPHERE Create a sphere containing 4 points.
%
%   s = createSphere(p1, p2, p3, p4);
%   return in s the sphere common to the 4 pointsp1, p2, p3 and p4.
%
%   Ref: P. Bourke
%   http://astronomy.swin.edu.au/~pbourke/geometry/spherefrom4/
%
%   See also
%   spheres, circles3d

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-03-22
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

if length(varargin)==4
    pts = [varargin{1};varargin{2};varargin{3};varargin{4}];
elseif length(varargin)==1
    pts = varargin{1};
else
    error('wrong number of arguments in createSphere');
end

m1 = det([pts ones(4,1)]);
s2 = sum(pts.*pts, 2);
m2 = det([s2 pts(:,2) pts(:,3) ones(4,1)]);
m3 = det([pts(:,1) s2 pts(:,3) ones(4,1)]);
m4 = det([pts(:,1) pts(:,2) s2 ones(4,1)]);

m5 = det([s2 pts]);

x0 = m2*.5/m1;
y0 = m3*.5/m1;
z0 = m4*.5/m1;
r  = sqrt(x0*x0 + y0*y0 + z0*z0 - m5/m1);

sphere = [x0 y0 z0 r];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Compare the speeds of various inference engines on the water DBN
seed = 0;
rand('state', seed);
randn('state', seed);

%bnet = mk_water_dbn;
bnet = mk_orig_water_dbn;

T = 3;
engine = {};
%engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));
%engine{end+1} = smoother_engine(hmm_2TBN_inf_engine(bnet));
%engine{end+1} = jtree_dbn_inf_engine(bnet);
engine{end+1} = jtree_unrolled_dbn_inf_engine(bnet, T);
engine{end+1} = cbk_inf_engine(bnet, 'clusters', {[1],[2],[3],[4],[5],[6],[7],[8]}); %ff
engine{end+1} = cbk_inf_engine(bnet, 'clusters', {[1 2],[3 4 5 6],[7 8]}); %manually designed marginally independent by BK
engine{end+1} = cbk_inf_engine(bnet, 'clusters', {[1:5], [3:7], [7:8]}); %manually designed conditionally independent by BK
engine{end+1} = cbk_inf_engine(bnet, 'clusters', {[1 3], [2 3 7], [3 5], [3 4 7], [6 7 8]}); %automatically found using TJTs offline
engine{end+1} = cbk_inf_engine(bnet, 'clusters', {[1 3 5], [2 3 5 7], [3 4 7], [4 6 7], [6 7 8]}); %automatically found using TJTs offline
engine{end+1} = cbk_inf_engine(bnet, 'clusters', {[1 3 4 5], [2 3 4 7 8], [4 6 7 8]}); %automatically found using TJTs offline

% bk_inf_engine yields exactly the same results for the marginally independent cases.
%engine{end+1} = bk_inf_engine(bnet, 'clusters', 'ff');
%engine{end+1} = bk_inf_engine(bnet, 'clusters', { [1 2], [3 4 5 6], [7 8] });

inf_time = cmp_inference_dbn(bnet, engine, T, 'exact', 1)
learning_time = cmp_learning_dbn(bnet, engine, T, 'exact', 1)

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% SUBSASGN Subscripted assignment.
%    A(I) = B assigns the values of B into the elements of A specified by
%    the subscript vector I.  B must have the same number of elements as I
%    or be a scalar.
%
%    A(I,J) = B assigns the values of B into the elements of the rectangular
%    submatrix of A specified by the subscript vectors I and J.  A colon used as
%    a subscript, as in A(I,:) = B, indicates all columns of those rows
%    indicated by vector I. Similarly, A(:,J) = B means all rows of columns J.
%
%    A(I,J,K,...) = B assigns the values of B to the submatrix of A specified
%    by the subscript vectors I, J, K, etc. A colon used as a subscript, as in
%    A(I,:,K) = B, indicates the entire dimension.
%
%    For both A(I,J) = B and the more general multi-dimensional
%    A(I,J,K,...) = B, B must be LENGTH(I)-by-LENGTH(J)-by-LENGTH(K)-... , or
%    be shiftable to that size by adding or removing singleton dimensions, or
%    contain a scalar, in which case its value is replicated to form a matrix
%    of that size.
%
%    A{I} = B when A is a cell array and I is a scalar places a copy of
%    the array B into the specified cell of A.  If I has more than one
%    element, this expression is an error.  Use [A{I}] = DEAL(B) to place
%    copies of B into multiple cells of A.  Multiple subscripts that
%    specify a scalar element, as in A{3,4} = magic(3), also work.
%
%    A(I).field = B when A is a structure array and I is a scalar places
%    a copy of the array B into the field with the name 'field'.  If I
%    has more than one element, this expression is an error. Use
%    [A(I).field] = DEAL(B) to place copies of B into multiple elements
%    of A.  If A is a 1-by-1 structure array, then the subscript can be
%    dropped.  In this case, A.field = B is the same as A(1).field = B.
%
%    A = SUBSASGN(A,S,B) is called for the syntax A(I)=B, A{I}=B, or
%    A.I=B when A is an object.  S is a structure array with the fields:
%        type -- character vector or string containing '()', '{}', or '.'
%                specifying the subscript type.
%        subs -- cell array, character vector, or string containing the
%                actual subscripts.
%
%    For instance, the syntax A(1:2,:)=B calls A=SUBSASGN(A,S,B) where
%    S is a 1-by-1 structure with S.type='()' and S.subs = {1:2,':'}. A
%    colon used as a subscript is passed as the character vector ':'.
%
%    Similarly, the syntax A{1:2}=B invokes A=SUBSASGN(A,S,B) where
%    S.type='{}' and the syntax A.field=B invokes SUBSASGN(A,S,B) where
%    S.type='.' and S.subs='field'.
%
%    These simple calls are combined in a straightforward way for
%    more complicated subscripting expressions.  In such cases
%    length(S) is the number of subscripting levels.  For instance,
%    A(1,2).name(3:5)=B invokes A=SUBSASGN(A,S,B) where S is 3-by-1
%    structure array with the following values:
%        S(1).type='()'       S(2).type='.'        S(3).type='()'
%        S(1).subs={1,2}      S(2).subs='name'     S(3).subs={3:5}
%
%    See also SUBSREF, SUBSTRUCT, PAREN, SUBSINDEX, LISTS,
%             NUMARGUMENTSFROMSUBSCRIPT.
%
%    Reference page in Doc Center
%       doc subsasgn
%
%    Other functions named subsasgn
%
%       axisobj/subsasgn             gbm/subsasgn
%       axistext/subsasgn            gpuArray/subsasgn
%       bm/subsasgn                  graph/subsasgn
%       calendarDuration/subsasgn    heston/subsasgn
%       categorical/subsasgn         hgbin/subsasgn
%       cev/subsasgn                 hwv/subsasgn
%       cir/subsasgn                 instrument/subsasgn
%       classregtree/subsasgn        scribehandle/subsasgn
%       codistributed/subsasgn       scribehgobj/subsasgn
%       Composite/subsasgn           sde/subsasgn
%       dataset/subsasgn             sdeddo/subsasgn
%       datetime/subsasgn            sdeld/subsasgn
%       diffusion/subsasgn           sdemrd/subsasgn
%       digraph/subsasgn             serial/subsasgn
%       distributed/subsasgn         sym/subsasgn
%       drift/subsasgn               tabular/subsasgn
%       duration/subsasgn            tall/subsasgn
%       fighandle/subsasgn           ts/subsasgn
%       figobj/subsasgn              tscollection/subsasgn
%       fints/subsasgn
%

author: jmaih


text:
function matirc()
% MATLAB_IRC - A small proof of concept IRC client written in MATLAB
%
% DESCRIPTION.  This is a pretty minimal proof-of-concept client.  There
% is not much error checking, or robustness.
%
% Supported commands:
% /msg  - send private message:  /msg <user> <message>
%         example: /msg nickserv identify hunter123
%
% /quit - disconnect and quit the client
%
%
% REQUIREMENTS:
%  * TCP/UDP/IP Toolbox from:
%  https://www.mathworks.com/matlabcentral/fileexchange/345-tcpudpip-toolbox-2-0-6
%
% xixor
% Feb. 13, 2012
%
% MATIRC 0.1
%
% website: http://matirc.xixor.net
% email: xixor@xixor.net

%% OPTIONS
host = 'irc.freenode.net';
port = 6667;
opt.nick       = ['matirc_user_' num2str(randi(1000,1))];  % Doesn't handle it if someone already has your nick
opt.channel    = '##matlab';
opt.username   = 'matirc';
opt.hostname   = 'www.xixor.net';
opt.servername = 'matirc';
opt.realname   = 'im too lame to edit matirc.m';
timer_period   = 0.1;

%% CONNECT
con=pnet('tcpconnect',host,port);
pnet(con,'setwritetimeout',5);
pnet(con,'setreadtimeout',5);

if con==-1
	error 'Bad url or server down.....';
end

disp(['MATLAB IRC Client connected to: ' host]);

%% setup the data buffer callback event to check if messages have arrived
t = timer('StartDelay', 0, 'Period', timer_period,  ...
        	'ExecutionMode', 'fixedRate');
t.TimerFcn = {@tcpip_data_callback, con};

start(t);

c = onCleanup(@()on_exit_event(t,con));

%% Connection Things

pause(2);  % TODO: need a better way to wait for the connection to be initialized
register(con,opt);
pause(2);
register(con,opt);

msg_text = opt.channel;
while 1
	user_msg = input([msg_text ' >> '],'s');
	if isempty(user_msg), continue, end;
	user_msg = strtrim(user_msg);

	if numel(user_msg) >=4 && strcmp(user_msg(1:4),'/msg')
		send_priv_msg(user_msg,con);
		continue;
	end

	switch user_msg
		case '/quit'
			break;
		case '/help'
			display_help();
		otherwise
			send_msg = ['PRIVMSG ' opt.channel ' : ' user_msg '\r\n'];
			pnet(con,'printf',send_msg);
	end
	end

function send_priv_msg(msg,con)

[~,msg] = strtok(msg,' ');
[nick,msg] = strtok(msg,' ');
msg = strtrim(deblank(msg));
nick = strtrim(deblank(nick));

send_msg = ['PRIVMSG ' nick ' : ' msg '\r\n'];
pnet(con,'printf',send_msg);
fprintf('\n    --> %s : %s\n', nick, msg);

function display_help()
fprintf('------------------------------------------\r\n');
fprintf('MATLAB IRC Client version 0.0.69-pre-alpha\r\n');
fprintf('Just type your message and it will be sent to the channel\r\n');
fprintf('To send a private message: /msg <user> <msg>\r\n');
fprintf('To quit type: /quit\r\n');
fprintf('------------------------------------------\r\n');

function on_exit_event(t,con)
% ON_EXIT_EVENT - callback for onCleanup
% This function deletes the timer and closes the TCP/IP connection object
%

fprintf('\r\n Quit(): cleaning timer and TCP/IP connection.\r\n');

if isvalid(t)
	stop(t);
	delete(t);
end
pnet(con,'close');

function tcpip_data_callback(~, ~, con)
% TCPIP CALLBACK - callback to check the TCP/IP data queue
% This function is executed by a timerand is called automatically
% INPUTS:
% con - connection object
% OUTPUTS:
% none
% no error checking or printable output.
%

data = pnet(con,'read',4096,'noblock');
if isempty(data), return, end

d = regexp(data,'\n','split');   % break the message into a cell array

for i=1:size(d,1)
	process_message(con,d{i});
end

function process_message(con,msg)
% PROCESS_MESSAGE - process irc messages received from the server
%
% INPUTS:
% con: connection object
% msg: the messgae to parse
%
% OUTPUTS: none.  This function calls other functions based on the message.
%
% TODO: need to parse irc messages in a general way

if ~isempty(strfind(msg,'PRIVMSG'))
	parse_privmsg(msg)
	return;
end

prefix = find(msg ==':',1,'first');
msg_name = strtrim(msg(1:prefix-1));
msg_body = strtrim(msg(prefix+1:end));

switch msg_name
	case 'NOTICE AUTH'
		% server messages
		fprintf('%s\n',msg_body);
	case 'PING'
		% ping requests from the server
		pnet(con,'printf',['PONG ' msg_body '\r\n']);
		fprintf('PING? PONG!\n');
	otherwise
		% TODO: parse other messages
end

function register(con,opt)
%REGISTER - set nick, user info
% INPUTS:
% con - connection object
% opt - options struct, members: username, hostname, servername, realname,
%       nick, channel
%
% OUTPUTS:
% no outputs or error checking
%

pnet(con,'printf',['NICK ' opt.nick '\r\n']);
pnet(con,'printf',['USER ' opt.username ' ' opt.hostname ' ' opt.servername ' :' opt.realname '\r\n']);
pnet(con,'printf',['JOIN ' opt.channel '\r\n']);

function parse_privmsg(msg,con,opt)
% PARSE_PRIVMSG - parse a PRIVMSG, print to screen
%
% INPUTS
% msg : a privmsg string from the irc server
% con : connection object
% opt : options structure
% (con and opt unused, left for future implementation)
%
% OUTPUT:
% prints the <nick> msg to the command window
%
% No error checking or confirmation
%

% TODO: need a much more general and robust way of parsing irc messages
% fprintf('\n MSG: %s\n',msg); // debug statement
sent_nick = msg(find(msg==':',1,'first')+1:find(msg=='!',1,'first')-1);

fprintf('<%s> %s\n',sent_nick,strtrim(msg(find(msg ==':',1,'last')+1:end)));


text:
function [Subcomponents,Population] = DividingDistanceVariables(Problem,NIA,DiverIndexes,ConverIndexes)
% Dividing distance variables

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

    %% Generate the initial population
    PopDec = zeros(Problem.N,Problem.D);
    % Generate the diverse variables by uniformly sampling method
    if sum(DiverIndexes) == 1
        PopDec(:,DiverIndexes) = (0:Problem.N-1)/(Problem.N-1);
    elseif sum(DiverIndexes) > 4
        PopDec(:,DiverIndexes) = rand(Problem.N,sum(DiverIndexes));
    else
        PopDec(:,DiverIndexes) = UDall(Problem.N,sum(DiverIndexes));
    end
    % Randomly generate the distance variables
    PopDec(:,ConverIndexes) = rand(Problem.N,sum(ConverIndexes));
    % Generate the initial population
    PopDec     = PopDec.*repmat(Problem.upper-Problem.lower,Problem.N,1) + repmat(Problem.lower,Problem.N,1);
    Population = Problem.Evaluation(PopDec);

    %% Interdependence analysis
    interaction = false(Problem.D);
    interaction(logical(eye(Problem.D))) = true;
    for i = 1 : Problem.D-1
        for j = i+1 : Problem.D
            drawnow('limitrate');
            for time2try = 1 : NIA
                % Detect whether the i-th and j-th decision variables are
                % interacting
                x    = randi(Problem.N);
                a2   = rand*(Problem.upper(i)-Problem.lower(i)) + Problem.lower(i);
                b2   = rand*(Problem.upper(j)-Problem.lower(j)) + Problem.lower(j);
                Decs = repmat(Population(x).dec,3,1);
                Decs(1,i) = a2;
                Decs(2,j) = b2;
                Decs(3,[i,j]) = [a2,b2];
                F = Problem.Evaluation(Decs);
                delta1 = F(1).obj - Population(x).obj;
                delta2 = F(3).obj - F(2).obj;
                interaction(i,j) = interaction(i,j) | any(delta1.*delta2<0);
                interaction(j,i) = interaction(i,j);
                % Update the solution
                if ConverIndexes(j) && all(F(2).obj <=Population(x).obj)
                    Population(x) = F(2);
                end
                if ConverIndexes(i) && all(F(1).obj <=Population(x).obj)
                    Population(x) = F(1);
                end
                if all(ConverIndexes([i,j])) && all(F(3).obj <=Population(x).obj)
                    Population(x) = F(3);
                end
            end
        end
    end

    %% Dividing distance variables
    Subcomponents = {};
    divided = false(1,Problem.D);
    while ~all(divided(ConverIndexes))
        x = find(~divided & ConverIndexes,1);
        while sum(any(interaction(x,ConverIndexes),1)) > length(x)
            x = find(any(interaction(x,:),1) & ConverIndexes);
        end
        Subcomponents = [Subcomponents,x];
        divided(x) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%MetabolicReaction
%
% @wholeCellModelID State_MetabolicReaction
% @name             Metabolic reaction
% @description
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 1/5/2011
classdef MetabolicReaction < edu.stanford.covert.cell.sim.CellState
    %property annotations
    properties (Constant)
        optionNames             = {   %names of properties that are options
            'verbosity';
            'seed';
            };
        fixedConstantNames      = {}; %names of process properties that are considered fixed constants
        fittedConstantNames     = {}; %names of process properties that are considered fitted constants, and should be stored with the simulation as such
        stateNames              = {   %names of properties which are part of the simulation's state
            'growth';
            'fluxs';
            };
        dependentStateNames     = {   %names of properties which can be calculated from the simulation's state
            'doublingTime'
            };
    end

    %constants
    properties
        reactionWholeCellModelIDs  %reaction whole cell model ids
        growth0                    %target growth rate
        initialGrowthFilterWidth   %tolerance of initial growth rate
        meanInitialGrowthRate      %mean initial growth rate
    end

    %state
    properties
        growth  %growth rate (cell/s)
        fluxs   %flux of each reaction, size: [length(reactionWholeCellModelIDs) X 1 X time] (rxn/s)
    end

    %dependent state
    properties (Dependent = true, SetAccess = protected)
        doublingTime  %doubling time (s)
    end

    %dependent state
    properties (Constant)
        dryWeight = 0; %dry weight of this class' state properties
    end

    %constructor
    methods
        function this = MetabolicReaction(wholeCellModelID, name)
            this = this@edu.stanford.covert.cell.sim.CellState(wholeCellModelID, name);
        end
    end

    methods
        function initializeConstants(this, knowledgeBase, simulation)
            this.initializeConstants@edu.stanford.covert.cell.sim.CellState(knowledgeBase, simulation);

            m = findobj(knowledgeBase.processes, 'wholeCellModelID', 'Process_Metabolism');
            this.reactionWholeCellModelIDs = {m.reactions.wholeCellModelID}';
        end
    end

    %allocate memory for state
    methods
        function allocateMemory(this, numTimePoints)
            this.growth = zeros(1, 1, numTimePoints);
            this.fluxs  = zeros(numel(this.reactionWholeCellModelIDs), 1, numTimePoints);
        end
    end

    %initialize state
    methods
        function initialize(~)
        end
    end

    %getters
    methods
        %doubling time (seconds/cell)
        function value = get.doublingTime(this)
            value = 1 ./ this.growth;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell