text:
function inpu2()
%inpu2
    %  This scriptfile ask for several selection parameters that can be set
    %  for part of the analysis. The original input parameters can be reset using
    %  the  Reset button in the main figure map
    %  Changes ZG.newt2
    %
    % make the interface
    %
    ZG=ZmapGlobal.Data;
    report_this_filefun(mfilename('fullpath'));
    ZG.newt2 = catalog_overview(ZG.newt2);
end
%{
h=zmap_message_center();
message_window_candidates = findobj('Name','Message Window');

mess = max([message_window_candidates.Number]);
clear message_window_candidates
if isempty(mess) == 1; mess = 1 ; end
figure_w_normalized_uicontrolunits(mess)

mint = min(ZG.newt2.Date);
maxt = max(ZG.newt2.Date);
clf
set(gca,'visible','off');
set(gcf,'Units','pixel','pos',[20 fipo(4)-350 400 300])
set(gcf,'Name','Selection Parameters');

inp5t=uicontrol('Style','edit','Position',[.70 .70 .22 .06],...
    'Units','normalized','String',char(min(ZG.newt2.Date)),...
    'Callback','mint=str2double(inp5t.String); inp5t.String=num2str(mint,6);');

inp6t=uicontrol('Style','edit','Position',[.70 .60 .22 .06],...
    'Units','normalized','String',char(max(ZG.newt2.Date)),...
    'Callback','maxt=str2double(inp6t.String); inp6t.String=num2str(maxt,6);');

inp5=uicontrol('Style','edit','Position',[.70 .50 .22 .06],...
    'Units','normalized','String',num2str(minma2),...
    'Callback','minma2=str2double(inp5.String); inp5.String=num2str(minma2);');

inp6=uicontrol('Style','edit','Position',[.70 .40 .22 .06],...
    'Units','normalized','String',num2str(maxma2),...
    'Callback','maxma2=str2double(inp6.String); inp6.String=num2str(maxma2);');

inp7=uicontrol('Style','edit','Position',[.30 .20 .10 .06],...
    'Units','normalized','String',num2str(minde),...
    'Callback','minde=str2double(inp7.String); inp7.String=num2str(minde);');

inp8=uicontrol('Style','edit','Position',[.60 .20 .10 .06],...
    'Units','normalized','String',num2str(maxde),...
    'Callback','maxde=str2double(inp8.String); inp8.String=num2str(maxde);');

close_button=uicontrol('Style','Pushbutton',...
    'Position',[.65 .02 .20 .10 ],...
    'Units','normalized','Callback',@(~,~)zmap_Message_center(),'String','cancel');

go_button=uicontrol('Style','Pushbutton',...
    'Position',[.35 .02 .20 .10 ],...
    'Units','normalized',...
    'Callback','zmap_message_center(); magcumu',...
    'String','Go');

txt6 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.02 0.76 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.s ,...
    'FontWeight','bold',...
    'String','Minimum Time: ');

txt6 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.02 0.64 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.s ,...
    'FontWeight','bold',...
    'String','Maximum Time: ');

txt3 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.15 0.98 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'FontWeight','bold',...
    'String','Cut the catalog');

txt6 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.02 0.52 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.s ,...
    'FontWeight','bold',...
    'String','Minimum Magnitude: ');

txt6 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.02 0.39 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.s ,...
    'FontWeight','bold',...
    'String','Maximum Magnitude: ');

txt7 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0.02 0.22 0 ],...
    'Rotation',0 ,...
    'FontSize',ZmapGlobal.Data.fontsz.s ,...
    'FontWeight','bold',...
    'String','                 Min Depth                    Max Depth   ');

set(gcf,'visible','on')
%}

author: CelsoReyes


text:
function varargout = mitm_viewer( filename, varargin )
%MITM_VIEWER More-Image-Than-Memory complex image viewer
%   MITM_VIEWER (FILENAME, 'PropertyName', PropertyValue, ...) takes a
%   file of complex data (or set of files of identical sizes) and displays
%   them in a MATLAB IMSHOW/IMAGESC style viewer, even if the file is too
%   large to fit into memory.  Any file type handled by OPEN_READER can be
%   used.
%
%   Really just a convenient GUI wrapper for the HG_MITM_VIEWER component.
%
%       Property name     Description
%       figureSize        Size of the figure to be drawn (if a new figure
%                         is to be created).  If figureSize is not
%                         specified, the largest figure size that will fit
%                         in the current screen with an integer zoom factor
%                         is used.  This propert is ignored if the parent
%                         property is set.
%       initialFrame      Set the initial frame viewed.  Default is 1.
%       initialRemap      Set the initial remap function.  Default is
%                         'densityremap'.
%       mode              'navigate' (default), 'aoi', 'points', or
%                         'polygon'. Navigate just browses images.  'aoi'
%                         adds a button that allows the users to select a
%                         rectangular AOI which is returned as the output
%                         parameter.  'points' selects a single point or
%                         set of points in the image.  'polygon' actually
%                         behaves the same as 'points' with an indefinite
%                         number points, but draws lines between the
%                         points.
%       numPoints         If mode is 'points', this is the number of points
%                         to select.  Default is 1.  Setting to 0 or Inf
%                         will allow an indefinite number of points (same
%                         use as MATLAB getpts function).  If mode is not
%                         'points', this parameter is ignored.
%       selectCallback    Function handle for function that should be run
%                         after selection.  Callback is a function that
%                         takes a single argument, the selection
%                         parameters.  If selectCallback is set, multiple
%                         AOI/points can be chosen.  If not, only one
%                         AOI/point is chosen and returned in the output
%                         argument.
%       closeAfterSelect  Default is false.  Determines whether the figure
%                         should close after the AOI/point is chosen.  If
%                         mode is 'navigate', or if parent or
%                         selectCallback are set, then this property is
%                         ignored.
%       undulationFilename   The name/path of the file containing WGS84 to
%                            geoid undulations.  This is used in the data
%                            cursor (see geoid_undulation) to display
%                            the height above ellipsoid rather then height
%                            above mean sea level.  If the file doesn't
%                            exist no HOE measurement is shown.
%
% Author: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Parse input parameters
p = inputParser;
p.KeepUnmatched=true;
p.addParamValue('figureSize', [], @(x) isequal(size(x),[1 2])||isequal(size(x),[2 1]));
p.addParamValue('mode', 'navigate', @(x) any(strcmpi(x,{'navigate', 'aoi', 'points','polygon'})));
p.addParamValue('numPoints',1,@(x) isscalar(x));
p.addParamValue('initialFrame',1, @(x) isscalar(x)&&isnumeric(x)&&(x>0));
p.addParamValue('initialRemap',[], @(x) exist(x)>0);
p.addParamValue('closeAfterSelect', false);
p.addParamValue('selectCallback',[],@(x) isempty(x)||isa(x,'function_handle'))
p.addParamValue('undulationFilename', 'Und_min2.5x2.5_egm2008_isw=82_WGS84_TideFree_SE');
p.FunctionName = 'MITM_VIEWER';
p.parse(varargin{:});

%% Open files
if ((nargin<1)||isempty(filename)) % If no filename was give, use dialog box to ask for one
    % Recall last interactively selected path used
    if ispref('matlab_sar_toolbox','last_used_directory')
        pathname = getpref('matlab_sar_toolbox','last_used_directory');
        if ~ischar(pathname)||~exist(pathname,'dir')
            pathname = pwd;
        end
    else
        pathname = pwd;
    end
    [filename,pathname]=uigetfile(sar_file_extensions('complex'),...
        'Select Input File', pathname,'MultiSelect','on');
    if(iscell(filename)), filename=sort(filename); end;
    setpref('matlab_sar_toolbox','last_used_directory',pathname);
else
    pathname=[];
end
if(isinmem(filename)) % Filename is actually an array in memory
    fullfilename{1}=filename;
elseif(iscell(filename)) % Multiple files requested
    for i=1:length(filename)
        fullfilename{i}=[pathname filename{i}];
    end
elseif(filename) % Single filename
    fullfilename{1}=[pathname filename];
else % filename=0.  Cancel was pressed, instead of a file being chosen.
    if ~strcmpi(p.Results.mode,'navigate')
        varargout={[],1};
    end
    return;
end

%% Setup layout
fig_hand=figure('MenuBar','none','Toolbar','none');
uip_hand = uipanel('Parent',fig_hand,'Position',[0 0 1 1],'BorderType','none');
mitm_hand = hg_mitm_viewer(uip_hand);
% Save image button
saveButton=javaObjectEDT('javax.swing.JButton');
saveButton.setIcon(javaObjectEDT('javax.swing.ImageIcon',...
    fullfile(toolboxdir('matlab'),['icons' filesep 'file_save.png'])));
saveButton.setToolTipText('Save As...');
saveButton.setFocusable(false); % Make consistent with MATLAB toolbars
set(handle(saveButton,'callbackproperties'),'ActionPerformedCallback',...
    @(hbutton, eventStruct, hfig) save_image());
% Metadata button
metaButton=javaObjectEDT('javax.swing.JButton','Metadata');
set(handle(metaButton,'callbackproperties'),'ActionPerformedCallback',@(obj,eventdata) show_metadata());
% Geojump button
geojumpButton=javaObjectEDT('javax.swing.JButton','Geojump');
set(handle(geojumpButton,'callbackproperties'),'ActionPerformedCallback',@(obj,eventdata) geojump());
% AOI button
if ~strcmpi(p.Results.mode,'navigate')
    switch lower(p.Results.mode)
        case 'aoi'
            button_str = 'Select AOI';
        case 'points'
            button_str = 'Select point(s)';
        case 'polygon'
            button_str = 'Select polygon';
    end
    aoiButton=javaObjectEDT('javax.swing.JButton',button_str);
    set(handle(aoiButton,'callbackproperties'),'ActionPerformedCallback',@(obj,eventdata) select_AOI());
end

% Open file
mitm_hand.Frame = p.Results.initialFrame;
mitm_hand.openFile(fullfilename, true);
metadata = mitm_hand.Metadata;

% Add buttons to toolbar
main_toolbar = mitm_hand.main_toolbar;
main_toolbar.addSeparator;
main_toolbar.add(metaButton);
try
    if ~isempty(point_slant_to_ground([1; 1], metadata{mitm_hand.Frame}))
        main_toolbar.add(geojumpButton);
    end
end
if(~strcmpi(p.Results.mode,'navigate'))
    main_toolbar.addSeparator;
    main_toolbar.add(aoiButton);
end
main_toolbar.addSeparator;
main_toolbar.add(saveButton);
main_toolbar.repaint;
main_toolbar.revalidate;

%% Set figure size appropriate for data
datasize=double([metadata{mitm_hand.Frame}.ImageData.NumCols metadata{mitm_hand.Frame}.ImageData.NumRows]);
if numel(metadata)==1 % Ugly.  Hardcoded to match the toolbar height found in hg_mitm_viewer.
    toolbar_height = 25;
else
    toolbar_height = 50;
end
if isempty(p.Results.figureSize)
    screensize=get(0,'ScreenSize');
    % Image must fit in full screen size minus toolbars and border
    max_img_size=((screensize(3:4)-screensize(1:2))-[0 toolbar_height]) - [20 50];
else
    max_img_size=p.Results.figureSize;
end
if isfield(metadata{mitm_hand.Frame},'CollectionInfo')&&...
        isfield(metadata{mitm_hand.Frame}.CollectionInfo,'CoreName')
    fig_name_str = metadata{mitm_hand.Frame}.CollectionInfo.CoreName;
elseif (numel(metadata)>1)
    fig_name_str = ['Image #' num2str(mitm_hand.Frame)];
else
    fig_name_str = '';
end
scale = max(ceil(datasize./max_img_size)); % Initial scale is lowest integer decimation that will fit
% Calculation of desired size is very specific to internal workings of hg_mitm_viewer
new_size = ceil(datasize./scale) + [0 toolbar_height] + 1;
set(fig_hand,'Name',fig_name_str,'Position',[20 20 new_size]);
% The resize callback executes in parallel with the current code.  Make
% sure resize completes, before we continue, since the following zoom is
% dependent on the new size.
tid = tic;
while(~isequal(getpixelposition(fig_hand),[20 20 new_size]))&&(toc(tid)<1)
end
mitm_hand.Zoom = 'fit';
if ~isempty(p.Results.initialRemap)
    mitm_hand.Remap = p.Results.initialRemap;
end

%% Will selection parameters be returned in output argument?  If so, wait
if(~strcmpi(p.Results.mode,'navigate')&&isempty(p.Results.selectCallback))
    % Should disable metadata button here until AOI selection is finished
    waitfor(uip_hand,'UserData','AOI_finished');
    if(~exist('varargout','var'))
        varargout={[],1};
    end
end

%% Callbacks for figure toolbar button functions
    % All of the work here is temporarily removing the toolbar from the
    % figure, so it doesn't get saved.  We assume that we must have a
    % figure with nothing else other than the mitm_viewer.
    function save_image()
        old_resize = get(uip_hand,'ResizeFcn');
        old_visible = get(mitm_hand.movie_toolbar,'Visible');
        set(uip_hand,'ResizeFcn',[]);
        fig_pos = getpixelposition(fig_hand);
        axes_pos = getpixelposition(mitm_hand.AxesHandle);
        setpixelposition(fig_hand,[fig_pos(1:2) axes_pos(3:4)]);
        setpixelposition(mitm_hand.AxesHandle,axes_pos);
        filemenufcn(fig_hand,'FileSaveAs');
        set(mitm_hand.movie_toolbar,'Visible',old_visible); % Why does FileSaveAs change this?
        set(uip_hand,'ResizeFcn',old_resize);
        setpixelposition(fig_hand,fig_pos); % Replace toolbar
    end

    % Roughly the reverse of complexdatacursorfun.  Takes a lat/long and
    % places a datacursor-- as opposed to complexdatacursorfun which takes
    % a datacursor and computes a lat/long.
    function geojump()
        full_meta = mitm_hand.Metadata;
        meta = full_meta{mitm_hand.Frame}; % For easier referencing
        % Get lat/long from user
        if isfield(meta,'GeoData')&&isfield(meta.GeoData,'SCP')
            % Default elevation is elevation of scene center
            def_el = num2str(meta.GeoData.SCP.LLH.HAE);
        else
            def_el = '';
        end
        answers = inputdlg({'Latitude','Longitude','Elevation (HAE, m)'},'Geojump',1,{'','',def_el});
        if isempty(answers), return, end; % Cancel was pressed
        lla = [latlonnum(answers{1});...
            latlonnum(answers{2});...
            str2double(answers{3})];
        if any(isnan(lla))
            msgbox('Invalid Lat/Lon format!');
            return;
        end
        if ~mitm_hand.geojump(lla)
            msgbox('Specified point is not in image!');
        end
    end

    % Display image metadata
    function show_metadata
        % Just opens the MATLAB variable editor for the metadata.  Assumes
        % user is at the MATLAB commandline with no debugger.
        % assignin('base','metadata',metadata{mitm_hand.Frame});
        % openvar('metadata');
        % Brings up new figure with tree view of metadata
        % Doesn't handle cell arrays or display 2-D matrices as well as
        % the MATLAB variable editor, but cleaner because its standalone
        % and doesn't write to base MATLAB workspace.
        MetaViewer(metadata{mitm_hand.Frame});
    end

    function select_AOI
        % Get AOI
        if strcmpi(p.Results.mode,'aoi') % Rectangular AOI
            aoi_screen=[0 0 0 0];
            while(any(aoi_screen(3:4)==0)),
                aoi_screen=getrect(mitm_hand.AxesHandle);
            end
        elseif strcmpi(p.Results.mode,'polygon')
            [x,y]=getline(mitm_hand.AxesHandle);
        elseif (p.Results.numPoints>0)&&isfinite(p.Results.numPoints) % Predetermined number of points
            [x,y]=ginput(p.Results.numPoints);
        else % Indefinite number of points
            [x,y]=getpts(mitm_hand.AxesHandle);
        end
        % Draw AOI
        if strcmpi(p.Results.mode,'aoi')
            h=rectangle('Position',aoi_screen,'LineWidth',2,'EdgeColor','r');
        elseif strcmpi(p.Results.mode,'polygon')
            h=line([x; x(1)],[y; y(1)],zeros(size(x)+[1 0]),'Color', 'r', 'LineWidth', 2);
        else % Points
            % Emulate getpts-style markers
            h=line(x,y,zeros(size(x)),'Color', 'c', 'LineStyle', 'none', 'Marker', '+', 'MarkerSize', 9);
            h2=line(x,y,'Color', 'm', 'LineStyle', 'none', 'Marker', 'x', 'MarkerSize', 9);
        end
        % Confirm selection
        if ~strcmpi(p.Results.mode,'aoi')
            if p.Results.numPoints~=1
                question_string='Are these the points you wish to select?';
            else
                question_string='Is this the point you wish to select?';
            end
        else
            question_string='Is this the area of interest you wish to select?';
        end
        button=questdlg(question_string,'Confirm AOI','Yes','No','Yes');
        if(strcmp(button,'Yes'))
            varargout={[],1};
            % Convert to image coordinate space
            if strcmpi(p.Results.mode,'aoi')
                ulcorner=round(mitm_hand.axescoords2native(aoi_screen(1:2)));
                lrcorner=round(mitm_hand.axescoords2native(aoi_screen(1:2) + aoi_screen(3:4) - 1));
                % Constrain to be within image dimensions
                ulcorner=max(1,ulcorner);
                lrcorner=min(datasize,lrcorner);
                % Return values (UL corner and size) to calling function
                varargout{1}(1:2)=ulcorner;
                varargout{1}(3:4)=lrcorner-ulcorner+1;
            else % Series of points
                for i=1:length(x)
                    varargout{1}(i, 1:2)=round(mitm_hand.axescoords2native([x(i) y(i)]));
                end
            end
            varargout{2}=mitm_hand.FrameSubFileIndices{mitm_hand.Frame}; % Which frame was AOI selected in?
            if isa(p.Results.selectCallback,'function_handle')
                try
                    p.Results.selectCallback(varargout{:});
                catch
                    delete(h); if exist('h2','var'), delete(h2); end;
                    rethrow(lasterror);
                end
                delete(h); if exist('h2','var'), delete(h2); end;
            else
                if (p.Results.closeAfterSelect)
                    close(fig_hand);
                else
                    set(uip_hand,'UserData','AOI_finished'); % If you want to keep window, do this
                    drawnow; % If you don't do this, everything freezes when you disable button
                    aoiButton.setEnabled(false);
                end
            end
        else
            delete(h); if exist('h2','var'), delete(h2); end;
        end
    end
end

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

author: ngageoint


text:
function [nodesBF, rangeLoc, rangeOut] = treeBFranges(wt,varargin)
%TREEBFRANGES Tree nodes output ranges in BF order
%   Usage: [nodesBF, rangeLoc, rangeOut] = treeBFranges(wt);
%          [nodesBF, rangeLoc, rangeOut] = treeBFranges(wt,'rev');
%
%   Input parameters:
%         wt       : Filterbank tree struct.
%   Output parameters:
%         nodesBF  : All nodes in a breadth-first order
%         rangeLoc : Local ranges of unconnected (terminal) outputs
%         rangeOut : Global ranges of unconnected (terminal) outputs
%
%   `[nodesBF, rangeLoc, rangeOut] = treeBFranges(wt)` is a helper function
%   extracting all nodes of a tree in a BF order (root and low-pass first)
%   (numeric array of indexes `nodesBF`), and two cell arrays of ranges of
%   outputs. Each element of `rangeLoc` specifies range of unconnected
%   outputs of a node with at the corresponding position in `nodesBF`.
%   Elements `rangeOut` specify contain the resulting global indexes
%   (in the resulting coefficient cell array) of such unconnected nodes.
%
%   `[nodesBF, rangeLoc, rangeOut] = treeBFranges(wt,'rev')` does the same
%   but the arrays are reversed.

nodesBF = nodeBForder(0,wt);
do_rev = 0;
if ~isempty(varargin(strcmp('rev',varargin)));
   nodesBF = fliplr(nodesBF);
   do_rev = 1;
end

rangeLoc = nodesLocOutRange(nodesBF,wt);

rangeOut = treeOutRange(wt);
if do_rev
    %rangeOut = nodesOutRange(nodesBF,wt);
    rangeOut = rangeOut(end:-1:1);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [Population,RankSolution] = EnvironmentalSelection(Population,N,alpha)

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

% This function is written by Mengjun Ming

    %% Parameter
    popSize    = length(Population);
    NumSeq     = (1:popSize)';
    RankConvg  = zeros(popSize,1);
    RankDivs   = zeros(popSize,1);

    %% Modify the infeasible solutions
    PopObj = Population.objs;
    PopCon = Population.cons;
    z      = min(PopObj,[],1);
    n      = max(PopObj,[],1);
    M      = length(z);
    [W,~]  = UniformPoint(N,M);
    [~,Region]     = min(pdist2(PopObj-z,W,'cosine'),[],2);
    PopObj_2       = PopObj;
    Infeasible_all = any(PopCon>0,2);
    if sum(Infeasible_all) ~= 0
        [~,Region_Fmax] = min(pdist2(n-z,W,'cosine'),[],2);
        PopObj_2(Infeasible_all,:) = repmat(n,sum(Infeasible_all),1)+sum(max(0,PopCon(Infeasible_all,:)),2)*W(Region_Fmax,:)/norm(W(Region_Fmax,:));
    end

    %% Non-dominated sorting
    CV       = sum(max(0,PopCon),2);
    Dominate = false(popSize);
    for i = 1 : popSize-1
        for j = i+1 : popSize
            if CV(i) < CV(j)
                Dominate(i,j) = true;
            elseif CV(i) > CV(j)
                Dominate(j,i) = true;
            else
                k = any(PopObj(i,:)<PopObj(j,:)) - any(PopObj(i,:)>PopObj(j,:));
                if k == 1
                    Dominate(i,j) = true;
                elseif k == -1
                    Dominate(j,i) = true;
                end
            end
        end
    end

    %% Calculate S(i)
    S = sum(Dominate,2);

    %% Calculate R(i)
    R = zeros(1,popSize);
    for i = 1 : popSize
        R(i) = sum(S(Dominate(:,i)));
    end
    FrontNo = R + 1;

    %% Calculate the crowding distance of each solution
    PopObj   = Population.objs;
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Distance = sort(Distance,2);
    CrowdDis = Distance(:,floor(sqrt(popSize)));

    %% Add a middle column
    MiddleLevel = zeros(popSize,1);

    %% Environmental selection -- convergence
    Next = FrontNo == 1;
    if sum(Next) <= N
        [~,indx_Convg] = sortrows([FrontNo',-CrowdDis]);
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
        MiddleLevel(Temp(Del)) = 1;
        [~,indx_Convg] = sortrows([FrontNo',MiddleLevel,-CrowdDis]);
    end
    RankConvg(indx_Convg) = NumSeq;

    %% Environmental selection -- diversity
    FrontNo_D = ones(popSize,1);
    for i = 1:size(W,1)
        index = find(Region==i);
        if (~isempty(index))
            Objs_temp = PopObj_2(index,:);
            g_temp = sum((Objs_temp-z).*W(i,:),2);
            [~,index_FrontNo_D] = sort(g_temp);
            FrontNo_D(index(index_FrontNo_D)) = (1:length(g_temp))';
        end
    end
    FrontNo_D           = FrontNo_D+Infeasible_all*popSize;
    [~,indx_divs]       = sortrows([FrontNo_D,-CrowdDis]);
    RankDivs(indx_divs) = NumSeq;

    %% Population for next generation
    RankSolution = alpha*RankConvg+(1-alpha)*RankDivs;
    [~,Rank]     = sort(RankSolution);
    Population   = Population(Rank(1:N));
    RankSolution = 1 : N;
end
function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
		% Incarcarea matricei A (cu un patrat magic)
A=magic(3);
		% Salvarea continutului matricei A intr-un fisier text
save patrmag.dat A -ascii

author: Sable

repo: mcbench-benchmarks


text:
function [newL2] = bestMap(L1,L2)
%bestmap: permute labels of L2 match L1 as good as possible
%   [newL2] = bestMap(L1,L2);

%===========
L1 = L1(:);
L2 = L2(:);
if size(L1) ~= size(L2)
    error('size(L1) must == size(L2)');
end

Label1 = unique(L1);
nClass1 = length(Label1);
Label2 = unique(L2);
nClass2 = length(Label2);

nClass = max(nClass1,nClass2);
G = zeros(nClass);
for i=1:nClass1
    for j=1:nClass2
        G(i,j) = length(find(L1 == Label1(i) & L2 == Label2(j)));
    end
end
[c,t] = hungarian(-G);
newL2 = zeros(size(L2));
for i=1:nClass2
    newL2(L2 == Label2(i)) = Label1(c(i));
end

return;

%=======backup old===========

L1 = L1 - min(L1) + 1;      %   min (L1) <- 1;
L2 = L2 - min(L2) + 1;      %   min (L2) <- 1;
%===========    make bipartition graph  ============
nClass = max(max(L1), max(L2));
G = zeros(nClass);
for i=1:nClass
    for j=1:nClass
        G(i,j) = length(find(L1 == i & L2 == j));
    end
end
%===========    assign with hungarian method    ======
[c,t] = hungarian(-G);
newL2 = zeros(nClass,1);
for i=1:nClass
    newL2(L2 == i) = c(i);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function doy = convert_to_doy (epoch)
    doy = mydatedoy(epoch);
end

%!test
%! epoch = 1;
%! doy = mydatedoy(epoch);
%! doy2 = convert_to_doy(epoch);
%! myassert(doy2, doy)


text:
function [molCollectionReport, newMolFilesDir] = obtainMetStructures(model, metList, outputDir, sources)
% Obtain MDL MOL files from various databases, including KEGG, HMDB, ChEBI,
% and PubChem. Alternatively, openBabel can be used to convert InChI
% strings or SMILES in MDL MOL files.
%
% USAGE:
%
% molCollectionReport = obtainMetStructures(model, mets, sources, standardisationApproach)
%
% INPUTS:
%    model: COBRA model with following fields:
%
%        * .S - The m x n stoichiometric matrix for the metabolic network.
%        * .mets - An m x 1 array of metabolite identifiers.
%        * .metInChIString - An m x 1 array of metabolite identifiers.
%        * .metSmiles - An m x 1 array of metabolite identifiers.
%        * .metVMHID - An m x 1 array of metabolite identifiers.
%        * .metCHEBIID - An m x 1 array of metabolite identifiers.
%        * .metKEGGID - An m x 1 array of metabolite identifiers.
%        * .metPubChemID - An m x 1 array of metabolite identifiers.
%        * .metHMDBID - An m x 1 array of metabolite identifiers.
%        * .metDrugbankID - An m x 1 array of metabolite identifiers.
%        * .metLipidmassID - An m x 1 array of metabolite identifiers.
%
% OPTIONAL INPUTS:
%    mets: List of metabolites to be download (Default: All)
%    outputDir: Directory that will contain the obtained metabolite structures.
%    sources: Sources where the MOL files will be obtained (Default: all).
%             The sources supported are:
%
%        1.- 'inchi' (requires openBabel)
%        2.- 'smiles' (requires openBabel)
%        3.- 'kegg' (https://www.genome.jp/)
%        4.- 'hmdb' (https://hmdb.ca/)
%        5.- 'pubchem' (https://pubchem.ncbi.nlm.nih.gov/)
%        6.- 'chebi' (https://www.ebi.ac.uk/)
%        7.- 'drugbank' (https://go.drugbank.com/)
%        8.- 'lipidmass' (https://www.lipidmaps.org/)
%
% OUTPUTS:
%    molCollectionReport: Report of the obtained MDL MOL files
%
%        * .metList - List of metabolites to be download
%        * .sources -list of sources of metabolite structures
%        * .structuresObtained -Total of metabolite structures obtained.
%        * .structuresObtainedPerSource -Boolean table with the metabolite
%               structures obtained by source.
%        * .databaseCoverage -Table indicating the coverage of metabolites
%               obtained by each of the sources.
%        * .idsToCheck -Id source from which no molecular structures could
%               be obtained due to a webTimeout, conversion error, or
%               inconsistent id.
%    newMolFilesDir: directory where new mol files were written, subdirectories for each source.

if nargin < 2 || isempty(metList)
    metList = unique(regexprep(model.mets, '(\[\w\])', ''));
else
    metList = unique(regexprep(metList, '(\[\w\])', ''));
end
if nargin < 3 || isempty(outputDir)
    outputDir = [pwd filesep];
else
    % Make sure input path ends with directory separator
    outputDir = [regexprep(outputDir,'(/|\\)$',''), filesep];
end

% Set directory
newMolFilesDir  = [outputDir 'metStructures' filesep];
if exist(newMolFilesDir, 'dir') == 0
    mkdir(newMolFilesDir)
end

if nargin < 4 || isempty(sources)
    sources = {'chebi'; 'drugbank'; 'hmdb'; 'inchi'; 'kegg'; 'lipidmaps'; 'pubchem'; 'smiles'};
else
    sources = sort(sources);
end
allSources = {'chebi'; 'drugbank'; 'hmdb'; 'inchi'; 'kegg'; 'lipidmaps'; 'pubchem'; 'smiles'};

% Check openbabel installation
if isunix || ispc
    obabelCommand = 'obabel';
else
    obabelCommand = 'openbabel.obabel';
end
[oBabelInstalled, ~] = system(obabelCommand);
if oBabelInstalled ~= 1
    oBabelInstalled = false;
    disp('obabel is not installed, two features cannot be used: ')
    disp('1 - Generation of SMILES, InChI and InChIkey')
    disp('2 - MOL file standardisation')
end

% if ismac || ispc
%     [oBabelInstalled, results] = system('obabel');
%     if find(contains(results,'Usage:')) % it is actually installed but input file is missing
%         oBabelInstalled = 1;
%     end
% else
%     [oBabelInstalled, ~] = system('openbabel.obabel');
% end
% if oBabelInstalled == 127
%     oBabelInstalled = 0;
% end

webTimeout = weboptions('Timeout', 60);

%% Obtain met data

% Obtain ID's
fields = fieldnames(model);

% chebi
chebiFieldBool = ~cellfun(@isempty, regexpi(fields, 'chebi'));
if any(chebiFieldBool)
    chebiIDs = model.(fields{chebiFieldBool});
else
    chebiIDs = cell(size(model.mets));
end

% drugbank
drugbankFieldBool = ~cellfun(@isempty, regexpi(fields, 'drugbank'));
if any(drugbankFieldBool)
    drugbankIDs = model.(fields{drugbankFieldBool});
else
    drugbankIDs = cell(size(model.mets));
end

% HMDB
hmdbFieldBool = ~cellfun(@isempty, regexpi(fields, 'hmdb'));
if any(hmdbFieldBool)
    hmdbIDs = model.(fields{hmdbFieldBool});
end

% inchi
inchiFieldBool = ~cellfun(@isempty, regexpi(fields, 'inchi'));
if any(inchiFieldBool) && isstruct(model.(fields{ismember(fields, 'inchi')}))
    inchis = model.(fields{ismember(fields, 'inchi')}).standardWithStereoAndCharge;
elseif any(inchiFieldBool) && ~isstruct(model.(fields{ismember(fields, 'inchi')}))
    inchis = model.(fields{ismember(fields, 'inchi')});
else
    inchis = cell(size(model.mets));
end

% KEGG
keggFieldBool = ~cellfun(@isempty, regexpi(fields, 'kegg'));
if any(keggFieldBool)
    if sum(keggFieldBool) > 1
        metFieldBool = ~cellfun(@isempty, regexpi(fields, 'met'));
        keggFieldBool = keggFieldBool & metFieldBool;
    end
    keggIDs = model.(fields{keggFieldBool});
else
    keggIDs = cell(size(model.mets));
end

% lipidmaps
lipidmapsFieldBool = ~cellfun(@isempty, regexpi(fields, 'lipidmaps'));
if any(lipidmapsFieldBool)
    lipidmapsIDs = model.(fields{lipidmapsFieldBool});
else
    lipidmapsIDs = cell(size(model.mets));
end

% PubChem
PubChemFieldBool = ~cellfun(@isempty, regexpi(fields, 'PubChem'));
if any(PubChemFieldBool)
    PubChemIDs = model.(fields{PubChemFieldBool});
else
    PubChemIDs = cell(size(model.mets));
end

% SMILES
smilesFieldBool = ~cellfun(@isempty, regexpi(fields, 'smiles'));
if any(smilesFieldBool)
    smiles = model.(fields{smilesFieldBool});
else
    smiles = cell(size(model.mets));
end

%% Obtain met structures

% Unique metabolites idexes
mets = regexprep(model.mets, '(\[\w\])', '');

% Obtain MDL MOL files
idsToCheck = {};
[inchiMsg, smilesMsg] = deal(true);
idMatrix = false(length(metList), length(allSources));
for i = 1:length(metList)

    % identify met in model and start matrix counter
    idx = find(ismember(mets, metList{i}));
    matrixCounter = 0;

    % ChEBI
    if ~isempty(chebiIDs{idx(1)})  && ismember({'chebi'}, sources)
        saveFileDir = [newMolFilesDir 'chebi' filesep];
        try
            molFile = webread(['https://www.ebi.ac.uk/chebi/saveStructure.do?defaultImage=true&chebiId=' num2str(chebiIDs{idx(1)}) '&imageId=0'], webTimeout);
            if ~isempty(regexp(molFile, 'M  END'))
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                fid2 = fopen([newMolFilesDir 'chebi' filesep metList{i} '.mol'], 'w');
                fprintf(fid2, '%s\n', molFile);
                fclose(fid2);
                idMatrix(i, 1) = true;
            end
        catch ME
            disp(ME.message)
            idsToCheck{end + 1, 1} = ['chebi - ' chebiIDs{idx(1)}];
        end
    end

    % drugbank
    if ~isempty(drugbankIDs{idx(1)})  && ismember({'drugbank'}, sources)
        saveFileDir = [newMolFilesDir 'drugbank' filesep];
        try
            if contains(drugbankIDs{idx(1)}, 'MET')
                molFile = webread(['https://go.drugbank.com/structures/metabolites/' num2str(drugbankIDs{idx(1)}) '.mol'], webTimeout);
            else
                molFile = webread(['https://go.drugbank.com/structures/small_molecule_drugs/' num2str(drugbankIDs{idx(1)}) '.mol'], webTimeout);
            end
            if ~isempty(regexp(molFile, 'M  END'))
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                fid2 = fopen([newMolFilesDir 'drugbank' filesep metList{i} '.mol'], 'w');
                fprintf(fid2, '%s\n', molFile);
                fclose(fid2);
                idMatrix(i, 2) = true;
            end
        catch ME
            disp(ME.message)
            idsToCheck{end + 1, 1} = ['drugbank - ' drugbankIDs{idx(1)}];
        end
    end

    % HMDB
    if ~isempty(hmdbIDs{idx(1)})  && ismember({'hmdb'}, sources)
        saveFileDir = [newMolFilesDir 'hmdb' filesep];
        try
            numbersID = hmdbIDs{idx(1)}(5:end);
            if size(numbersID, 2) < 7
                numbersID = [repelem('0', 7 - size(numbersID, 2)) numbersID];
            end
            molFile = webread(['https://hmdb.ca/structures/metabolites/HMDB' numbersID '.mol'], webTimeout);
            if ~isempty(regexp(molFile, 'M  END'))
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                fid2 = fopen([newMolFilesDir 'hmdb' filesep metList{i} '.mol'], 'w');
                fprintf(fid2, '%s\n', molFile);
                fclose(fid2);
                idMatrix(i, 3) = true;
            end
        catch ME
            disp(ME.message)
            idsToCheck{end + 1, 1} = ['hmdb - ' hmdbIDs{idx(1)}];
        end
    end

    % InChI
    if ~isempty(inchis{idx(1)}) && ismember({'inchi'}, sources)
        if oBabelInstalled
            try
                saveFileDir = [newMolFilesDir 'inchi' filesep];
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                newFormat = openBabelConverter(inchis{idx(1)}, 'mol', [saveFileDir ...
                    metList{i} '.mol']);
                idMatrix(i, 4) = true;
            catch ME
                disp(ME.message)
                idsToCheck{end + 1, 1} = ['inchi - ' keggIDs{idx(1)}];
            end
        elseif inchiMsg && ~oBabelInstalled
            inchiMsg = false;
            display('OpenBabel is not installed to convert InChIs')
        end
    end

    % KEGG
    if ~isempty(keggIDs{idx(1)})  && ismember({'kegg'}, sources)
        saveFileDir = [newMolFilesDir 'kegg' filesep];
        try
            switch keggIDs{idx(1)}(1)
                case 'C'
                    molFile = webread(['https://www.genome.jp/dbget-bin/www_bget?-f+m+compound+' keggIDs{idx(1)}], webTimeout);
                case 'D'
                    molFile = webread(['https://www.kegg.jp/dbget-bin/www_bget?-f+m+drug+' keggIDs{idx(1)}], webTimeout);
            end
            if ~isempty(regexp(molFile, 'M  END'))
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                fid2 = fopen([newMolFilesDir 'kegg' filesep metList{i} '.mol'], 'w');
                fprintf(fid2, '%s\n', molFile);
                fclose(fid2);
                idMatrix(i, 5) = true;
            end
        catch ME
            disp(ME.message)
            idsToCheck{end + 1, 1} = ['kegg - ' keggIDs{idx(1)}];
        end
    end

    % lipidmaps
    if ~isempty(lipidmapsIDs{idx(1)})  && ismember({'lipidmaps'}, sources)
        saveFileDir = [newMolFilesDir 'lipidmaps' filesep];
        try
            molFile = webread(['https://www.lipidmaps.org/databases/lmsd/' num2str(lipidmapsIDs{idx(1)}) '?format=mdlmol'], webTimeout);
            if ~isempty(regexp(molFile, 'M  END'))
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                fid2 = fopen([newMolFilesDir 'lipidmaps' filesep metList{i} '.mol'], 'w');
                fprintf(fid2, '%s\n', molFile);
                fclose(fid2);
                idMatrix(i, 6) = true;
            end
        catch ME
            disp(ME.message)
            idsToCheck{end + 1, 1} = ['lipidmaps - ' lipidmapsIDs{idx(1)}];
        end
    end

    % PubChem
    if ~isempty(PubChemIDs{idx(1)})  && ismember({'pubchem'}, sources)
        saveFileDir = [newMolFilesDir 'pubchem' filesep];
        try
            molFile = webread(['https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/CID/'...
                num2str(PubChemIDs{idx(1)}) ...
                '/record/SDF/?record_type=2d&response_type=display'], webTimeout);
            if ~isempty(regexp(molFile, 'M  END'))
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                molFile(regexp(molFile, 'M  END') + 6:end) = [];
                fid2 = fopen([newMolFilesDir 'pubchem' filesep metList{i} '.mol'], 'w');
                fprintf(fid2, '%s\n', molFile);
                fclose(fid2);
                idMatrix(i, 7) = true;
            end
        catch ME
            disp(ME.message)
            idsToCheck{end + 1, 1} = ['pubchem - ' PubChemIDs{idx(1)}];
        end
    end

    % SMILES
    if ~isempty(inchis{idx(1)}) && ismember({'smiles'}, sources)
        if oBabelInstalled
            try
                saveFileDir = [newMolFilesDir 'smiles' filesep];
                if exist(saveFileDir, 'dir') == 0
                    mkdir(saveFileDir)
                end
                newFormat = openBabelConverter(smiles{idx(1)}, 'mol', [saveFileDir ...
                    metList{i} '.mol']);
                idMatrix(i, 8) = true;
            catch ME
                disp(ME.message)
                idsToCheck{end + 1, 1} = ['inchi - ' keggIDs{idx(1)}];
            end
        elseif smilesMsg && ~oBabelInstalled
            smilesMsg = false;
            display('OpenBabel is not isntalled to convert SMILES')
        end
    end
end

%% Report

% Delete databases not included
idMatrix(:, ~ismember(allSources, sources)) = [];

% molCollectionReport
molCollectionReport.metList = metList;
molCollectionReport.sources = sources;
molCollectionReport.structuresObtained = sum(any(idMatrix'));

% Structures obtained
nRows = size(idMatrix, 1);
varTypes = ['string', repmat({'logical'}, 1, size(idMatrix, 2))];
varNames = ['mets'; sources];
structuresObtainedPerSource = table('Size', [nRows length(varTypes)], 'VariableTypes', varTypes, 'VariableNames', varNames);
structuresObtainedPerSource.mets = metList;
for i = 1:length(sources)
    structuresObtainedPerSource.(sources{i}) = idMatrix(:, i);
end
molCollectionReport.structuresObtainedPerSource = structuresObtainedPerSource;

% Database coverage table
[nCols nRows] = size(idMatrix);
varTypes = {'string', 'double', 'double', 'double'};
varNames = {'sources', 'coverage', 'metsWithStructure', 'metsWithoutStructure'};
databaseCoverage = table('Size', [nRows length(varTypes)], 'VariableTypes', varTypes, 'VariableNames', varNames);
databaseCoverage.sources = sources;

%if size(idMatrix, 2)> 0
%    databaseCoverage.metsWithStructure = sum(idMatrix)';
%    databaseCoverage.metsWithoutStructure = sum(~idMatrix)';
%    databaseCoverage.coverage = round((databaseCoverage.metsWithStructure * 100) / size(idMatrix, 1), 2);
%    molCollectionReport.databaseCoverage = databaseCoverage;
%end

if nCols > 1
    databaseCoverage.metsWithStructure = sum(idMatrix)';
    databaseCoverage.metsWithoutStructure = sum(~idMatrix)';
else
    databaseCoverage.metsWithStructure = double(idMatrix)';
    databaseCoverage.metsWithoutStructure = double(~idMatrix)';
end
databaseCoverage.coverage = round((databaseCoverage.metsWithStructure * 100) / size(idMatrix, 1), 2);
molCollectionReport.databaseCoverage = databaseCoverage;

molCollectionReport.idsToCheck = idsToCheck;
if ~isempty(idsToCheck)
    disp('The following structures could not be obtained')
    disp(idsToCheck)
end
end

author: opencobra

repo: cobratoolbox


text:
%tRNA Aminoacylation
%
% @wholeCellModelID Process_tRNAAminoacylation
% @name             tRNA Aminoacylation
% @description
%   Biology
%   ===============
%   In biological systems tRNAs serve as mediators between the ribosome and
%   the amino acids which forms peptide polymers. In this process we
%   simulate the conjugation of amino acids to the tRNAs which deliver them
%   to the ribosome. Additionally we simulate the aminoacylation of the
%   tmRNA which similarly delivers the amino acid alanine to stalled
%   ribosomes.
%
%   Knowledge Base
%   ===============
%   As of 8/11/2010 the M. genitalium knowledge base contains 39 tRNA
%   aminoacylation reactions involving
%   - 37 aminoacylation reactions, each assuming a cost of 1 ATP per
%     aminoacylation
%   - 2 transfer reactions
%   - 36 tRNAs and 1 tmRNA
%   - 21 enzymes
%   - 20 amino acid substrates and 10 additional substrates
%
%   The reaction kinetics stored in the knowledge base were compiled
%   several sources including SABIO-RK [PUB_0100].
%
%   Note: Unlike knowledge base representation, for transfer reactions
%   reactionStoichiometryMatrix doesn't include the already conjugated
%   amino acid on the left-hand-side or the the ultimate conjugated amino
%   acid on the right-hand-side.
%
%   Representation
%   ===============
%   substrates represents the counts of free metabolites available for tRNA
%   aminoacylation. enzymes represents the counts of proteins available to
%   catalyze tRNA aminoacylation reactions. freeRNAs and aminoacylatedRNAs
%   represent the counts of free, unaminoacylated and aminocylated RNAs. We do
%   not represent intermediate states in the aminoacylation of tRNAs.  The
%   molecular weights of the unaminoacylated and aminocylated tRNAs are computed
%   by the knowledge base RNA classes.
%
%   reactionStoichiometryMatrix, reactionModificationMatrix,
%   reactionCatalysisMatrix, and enzymeBounds represent the tRNA aminoacylation
%   reactions. reactionStoichiometryMatrix represents the free metabolites
%   required for  each reaction. reactionModificationMatrix represents the tRNA
%   aminoacylated by each reaction. reactionCatalysis represents the enzyme
%   required to aminoacylate each tRNA. enzymeBounds represents the kcat of the
%   catalyzing enzyme of each reaction. Note: reactionStoichiometryMatrix here
%   is slightly modified from that computed by the super class: amino acids
%   conjugated to tRNAs before and after transfer reactions have been removed
%   from the left- and righ-hand-side of reactionStoichiometryMatrix.
%
%   Initialization
%   ===============
%   tRNAs are all initialized to the aminoacylated state.
%
%   Simulation
%   ===============
%   Uses greedy algorithm to simulate tRNA aminacylation (complexation of
%   tRNAs with the specific amino acids that aminoacylate them).
%   1. Deterministically activate tRNAs up to the minimum of free tRNAs and
%      amino acids, proportional to free tRNAs
%   2. Stoichastically activate residual tRNAs using residual amino acids
%      with probabilities proportional to remaining free tRNAs
%
%      While(true)
%        1. Calculate numbers of protein monomers that can be modified based on
%           substrate, enzyme, and unmodified protein monomer availability and
%           kinetics.
%        2. Randomly select protein monomer to modify weighted by limits
%           calculated in step (1).
%        3. Update substrates, enzymes, unmodified protein monomers
%        4. Repeat until insufficient resources to further modify protein
%           monomers
%      End
%
% Author: Markus Covert, mcovert@stanford.edu
% Author: Jayodita Sanghvi, jayodita@stanfod.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 8/11/2010

classdef tRNAAminoacylation < edu.stanford.covert.cell.sim.ReactionProcess
    %property annotations
    properties (Constant)
        optionNames__              = {}; %names of option properties
		fixedConstantNames__       = {   %names of fixed constant properties
			'monomerTRNACounts';
			};
        fittedConstantNames__      = {}; %names of fitted constant properties
        localStateNames__          = {   %names of simulation state properties redundant with timecourses in this or other processes or the simulation
            'freeRNAs';
            'aminoacylatedRNAs'};
    end

    %IDs, names, and local indices
    properties
        substrateIndexs_aminoAcids                   %indices within substrates of amino acids
        substrateIndexs_glutamate                    %index within substrates of glutamate
        substrateIndexs_glutamine                    %index within substrates of glutamine
        substrateIndexs_methionine                   %index within substrates of methionine
        substrateIndexs_fmethionine                  %index within substrates of formylmethionine
        substrateIndexs_atp                          %index within substrates of ATP
        substrateIndexs_amp                          %index within substrates of AMP
        substrateIndexs_adp                          %index within substrates of ADP
        substrateIndexs_diphosphate                  %index within substrates of diphosphate
        substrateIndexs_phosphate                    %index within substrates of inorganic phosphate
        substrateIndexs_water                        %index within substrates of water
        substrateIndexs_hydrogen                     %index within substrates of hydrogen
        substrateIndexs_fthf10                       %index within substrates of 10-forymltetrahydrofolate
        substrateIndexs_thf                          %index within substrates of tetrahydrofolate

        enzymeIndexs_tRNASynthetases                 %indices within enzymes of tRNA synthetases
        enzymeIndexs_tRNATransferases                %indices within enzymes of tRNA transferases
        enzymeIndexs_tRNAGlutamylSynthetase          %index within enzymes of tRNA glutamyl synthetase
        enzymeIndexs_tRNAMethionylSynthetase         %index within enzymes of tRNA methionyl synthetase
        enzymeIndexs_tRNAGlutamylAmidotransferase    %index within enzymes of tRNA glutamyl amidotransferase
        enzymeIndexs_tRNAMethionylFormyltransferase  %index within enzymes of tRNA methionyl formyltransferase
        enzymeIndexs_tRNAs                           %index within enzymes of tRNAs
        enzymeIndexs_tmRNA                           %index within enzymes of tmRNA

        reactionIndexs_aminoacylation                %indices within reactions of tRNA aminoacylation reactions
        reactionIndexs_transfer                      %indices within reactions of tRNA transfer reactions
        reactionIndexs_glutamylamidotransfer         %index within reactions of tRNA glutamylamidotransfer reaction
        reactionIndexs_methionylformyltransfer       %index within reactions of tRNA methionylformyltransfer reaction

        compartmentWholeCellModelIDs                 %whole cell model ids of compartments
        compartmentIndexs_cytosol                    %index within compartments of cytosol

        freeRNAWholeCellModelIDs                     %whole cell model ids of free RNAs
        aminoacylatedRNAWholeCellModelIDs            %whole cell model ids of aminoacylated RNAs

        aminoacylatedRNAGlobalIndexs                 %indices of aminoacylated RNAs within simulation.matureRNAIndexs
        aminoacylatedRNAtRNAIndexs                   %indices of aminoacylated tRNAs within aminoacylatedRNAGlobalIndexs

        speciesIndexs_enzymes                        %indexs of enzyme species within speciesReactantByproductMatrix, speciesReactantMatrix
    end

    %fixed biological constants
    properties
        monomerTRNACounts              %numbers of each tRNA species required to translate each protein monomer species

        speciesReactantByproductMatrix %stoichiometry of susbtrates, enzymes, RNA in modifications necessary to mature each modifying RNA: [RNA modifications] X [subtrates; enzymes; free RNAs]
        speciesReactantMatrix          %reactant stoichiometry of susbtrates, enzymes, RNA in modifications necessary to mature each modifying RNA: [RNA modifications] X [subtrates; enzymes; free RNAs]
    end

    %global state (stored locally for convenience)
    properties
        freeRNAs          %counts of free RNAs
        aminoacylatedRNAs %counts of aminoacylated RNAs
    end

    %dependent global state (implemented as dependent property for
    %convenience)
    properties (Dependent)
        tRNASynthetases   %counts of tRNA synthetases
        tRNATransferases  %counts of tRNA transferases
    end

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

    %communication between process/simulation
    methods
        %initialize constants
        function initializeConstants(this, knowledgeBase, simulation, varargin)
            %import classes
            import edu.stanford.covert.util.ComputationUtil;

            %call super class method to get knowledge base reactions
            this.initializeConstants@edu.stanford.covert.cell.sim.ReactionProcess(knowledgeBase, simulation, varargin{:});

            %compartments
            this.compartmentWholeCellModelIDs = this.compartment.wholeCellModelIDs(this.compartment.cytosolIndexs);
            this.compartmentIndexs_cytosol    = 1;

            %create representation for GLN which super class leaves off
            %because it cancels out in the knowledge base
            this.substrateWholeCellModelIDs  = [this.substrateWholeCellModelIDs; 'GLN'];
            this.reactionStoichiometryMatrix = [this.reactionStoichiometryMatrix; zeros(1, size(this.reactionStoichiometryMatrix, 2))];
            this.reactionCoenzymeMatrix      = [this.reactionCoenzymeMatrix zeros(size(this.reactionCoenzymeMatrix, 1), 1)];

            %substrate indices
            this.substrateIndexs_aminoAcids  = this.substrateIndexs({'ALA'; 'ARG'; 'ASN'; 'ASP'; 'CYS'; 'GLN'; 'GLU'; 'GLY'; 'HIS'; 'ILE'; 'LEU'; 'LYS'; 'MET'; 'PHE'; 'PRO'; 'SER'; 'THR'; 'TRP'; 'TYR'; 'VAL'});
            this.substrateIndexs_glutamate   = this.substrateIndexs({'GLU'});
            this.substrateIndexs_glutamine   = this.substrateIndexs({'GLN'});
            this.substrateIndexs_methionine  = this.substrateIndexs({'MET'});
            this.substrateIndexs_fmethionine = this.substrateIndexs({'FMET'});
            this.substrateIndexs_atp         = this.substrateIndexs({'ATP'});
            this.substrateIndexs_amp         = this.substrateIndexs({'AMP'});
            this.substrateIndexs_adp         = this.substrateIndexs({'ADP'});
            this.substrateIndexs_diphosphate = this.substrateIndexs({'PPI'});
            this.substrateIndexs_phosphate   = this.substrateIndexs({'PI'});
            this.substrateIndexs_water       = this.substrateIndexs({'H2O'});
            this.substrateIndexs_hydrogen    = this.substrateIndexs({'H'});
            this.substrateIndexs_fthf10      = this.substrateIndexs({'FTHF10'});
            this.substrateIndexs_thf         = this.substrateIndexs({'THF'});

            %include tRNAs as enzymes
            this.enzymeWholeCellModelIDs = [
                this.enzymeWholeCellModelIDs;
                this.rna.wholeCellModelIDs(this.rna.matureIndexs(this.rna.matureTRNAIndexs));
                'MG_0004'];

            %enzyme indices
            this.enzymeIndexs_tRNASynthetases = this.enzymeIndexs({
                'MG_292_TETRAMER';       %alanyl-tRNA synthetase
                'MG_378_MONOMER';        %arginyl-tRNA synthetase
                'MG_036_DIMER';          %aspartyl-tRNA synthetase
                'MG_113_DIMER';          %asparaginyl-tRNA synthetase
                'MG_253_MONOMER';        %cysteinyl-tRNA synthetase
                'MG_462_MONOMER';        %glutamyl-tRNA synthetase
                'MG_251_DIMER';          %glycyl-tRNA synthetase
                'MG_035_DIMER';          %histidyl-tRNA synthetase
                'MG_345_MONOMER';        %isoleucyl-tRNA synthetase
                'MG_266_MONOMER';        %leucyl-tRNA synthetase
                'MG_136_DIMER';          %lysyl-tRNA synthetase
                'MG_021_DIMER';          %methionyl-tRNA synthetase
                'MG_194_195_TETRAMER';   %phenylalanyl-tRNA synthetase
                'MG_283_DIMER';          %prolyl-tRNA synthetase
                'MG_005_DIMER';          %seryl-tRNA synthetase
                'MG_375_DIMER';          %threonyl-tRNA synthetase
                'MG_126_DIMER';          %tryptophanyl-tRNA synthetase
                'MG_455_DIMER';          %tyrosyl-tRNA synthetase
                'MG_334_MONOMER'});      %valyl-tRNA synthetase
            this.enzymeIndexs_tRNATransferases = this.enzymeIndexs({
                'MG_098_099_100_TRIMER'; %glutamyl-tRNA(Gln) and/or aspartyl-tRNA(Asn) amidotransferase
                'MG_365_MONOMER'});      %methionyl-tRNA formyltransferase
            this.enzymeIndexs_tRNAGlutamylSynthetase         = this.enzymeIndexs({'MG_462_MONOMER'});        %glutamyl-tRNA synthetase
            this.enzymeIndexs_tRNAMethionylSynthetase        = this.enzymeIndexs({'MG_021_DIMER'});          %methionyl-tRNA synthetase
            this.enzymeIndexs_tRNAGlutamylAmidotransferase   = this.enzymeIndexs({'MG_098_099_100_TRIMER'}); %glutamyl-tRNA(Gln) and/or aspartyl-tRNA(Asn) amidotransferase
            this.enzymeIndexs_tRNAMethionylFormyltransferase = this.enzymeIndexs({'MG_365_MONOMER'});        %methionyl-tRNA formyltransferase
            this.enzymeIndexs_tRNAs                          = this.enzymeIndexs(this.rna.wholeCellModelIDs(this.rna.matureIndexs(this.rna.matureTRNAIndexs)));
            this.enzymeIndexs_tmRNA                          = this.enzymeIndexs({'MG_0004'});

            %RNAs
            this.reactionModificationMatrix(isnan(this.reactionModificationMatrix)) = 1;

            aaGenes = any(this.reactionModificationMatrix, 1)';
            aaRNAs  = ComputationUtil.invertCompositionMatrix(this.rna.matureRNAGeneComposition) * aaGenes;
            this.aminoacylatedRNAGlobalIndexs  = find(aaRNAs);
            this.aminoacylatedRNAtRNAIndexs    = intersect(this.aminoacylatedRNAGlobalIndexs, this.rna.matureTRNAIndexs);

            this.freeRNAWholeCellModelIDs          = this.rna.wholeCellModelIDs(this.rna.matureIndexs(this.aminoacylatedRNAGlobalIndexs));
            this.aminoacylatedRNAWholeCellModelIDs = this.rna.wholeCellModelIDs(this.rna.aminoacylatedIndexs(this.aminoacylatedRNAGlobalIndexs));

            this.reactionModificationMatrix = this.reactionModificationMatrix(:, aaGenes);

            %reaction indices
            this.reactionIndexs_aminoacylation          = find(strcmp(this.reactionTypes, 'aminoacylation'));
            this.reactionIndexs_transfer                = find(strcmp(this.reactionTypes, 'transfer'));
            this.reactionIndexs_glutamylamidotransfer   = this.reactionIndexs({'MG502_Amidotransferase'});
            this.reactionIndexs_methionylformyltransfer = this.reactionIndexs({'MG488_Formyltransferase'});

            %for transfer reactions unlike knowledge base representation,
            %formulate reactionStoichiometryMatrix so that left-hand-side
            %doesn't include the already conjugated amino acid, and the
            %right-hand-side doesn't include the ultimate conjugated amino
            %acid
            this.reactionStoichiometryMatrix(this.substrateIndexs_glutamine,   this.reactionIndexs_glutamylamidotransfer)   = -1;
            this.reactionStoichiometryMatrix(this.substrateIndexs_glutamate,   this.reactionIndexs_glutamylamidotransfer)   =  1;
            this.reactionStoichiometryMatrix(this.substrateIndexs_methionine,  this.reactionIndexs_methionylformyltransfer) =  0;
            this.reactionStoichiometryMatrix(this.substrateIndexs_fmethionine, this.reactionIndexs_methionylformyltransfer) =  0;

            %re-run super class method to update substrate and enzyme
            %mappings to include GLN, the tRNAs, and the tmRNA
            this.initializeConstants@edu.stanford.covert.cell.sim.Process(knowledgeBase, simulation, varargin{:});

            this.initializeSpeciesNetwork();

            %protein monomers -- needed for fitting
            this.monomerTRNACounts = knowledgeBase.proteinMonomerTRNACounts;
        end

        %initialize reaction stoichiomety to be used in evolveState:
        %  [tRNA aminoacylations and transfers] X [subtrates; enzymes; free tRNAs]
        function initializeSpeciesNetwork(this)
            % formulate network: [tRNA aminoacylations and transfers] X [subtrates; enzymes; free tRNAs]
            numMetabolites = size(this.reactionStoichiometryMatrix, 1);
            numEnzymes     = size(this.reactionCatalysisMatrix, 2);
            numRNAs        = size(this.reactionModificationMatrix, 2);

            this.speciesReactantByproductMatrix = [this.reactionModificationMatrix' * [ ...
                -this.reactionStoichiometryMatrix' ...
                this.reactionCatalysisMatrix ./ repmat(this.enzymeBounds(:, 2) * this.stepSizeSec, 1, numEnzymes)] ...
                eye(numRNAs)];
            this.speciesReactantMatrix = [this.reactionModificationMatrix' * max(0, [ ...
                -this.reactionStoichiometryMatrix' ...
                this.reactionCatalysisMatrix ./ repmat(this.enzymeBounds(:, 2) * this.stepSizeSec, 1, numEnzymes)]) ...
                eye(numRNAs)];

            this.speciesIndexs_enzymes = (1:numEnzymes)' + numMetabolites;
        end

        %retrieve state from simulation
        function copyFromState(this)
            this.copyFromState@edu.stanford.covert.cell.sim.ReactionProcess();

            this.freeRNAs          = this.rna.counts(this.rna.matureIndexs(       this.aminoacylatedRNAGlobalIndexs), this.compartment.cytosolIndexs, :);
            this.aminoacylatedRNAs = this.rna.counts(this.rna.aminoacylatedIndexs(this.aminoacylatedRNAGlobalIndexs), this.compartment.cytosolIndexs, :);
        end

        %send state to simulation
        function copyToState(this)
            this.copyToState@edu.stanford.covert.cell.sim.ReactionProcess();

            this.rna.counts(this.rna.matureIndexs(       this.aminoacylatedRNAGlobalIndexs), this.compartment.cytosolIndexs, :) = this.freeRNAs;
            this.rna.counts(this.rna.aminoacylatedIndexs(this.aminoacylatedRNAGlobalIndexs), this.compartment.cytosolIndexs, :) = this.aminoacylatedRNAs;
        end
    end

    %allocate memory for state
    methods
        function allocateMemoryForState(this, numTimePoints)
            this.allocateMemoryForState@edu.stanford.covert.cell.sim.ReactionProcess(numTimePoints);

            this.freeRNAs          = zeros(length(this.aminoacylatedRNAGlobalIndexs), 1, numTimePoints);
            this.aminoacylatedRNAs = zeros(length(this.aminoacylatedRNAGlobalIndexs), 1, numTimePoints);
        end
    end

    %model
    methods
        %Calculate
        %- contribution to FBA objective
        %- minimum expression consistent with cell cycle length
        function [bmProd, byProd, minEnzExp, maxEnzExp] = calcResourceRequirements_LifeCycle(this, ~, states)
            %% substrate and byproducts
            %aminoacylations of each t(m)RNA
            [~, idxs] = ismember(this.aminoacylatedRNAtRNAIndexs, this.aminoacylatedRNAGlobalIndexs);
            aminoacylations = zeros(size(this.aminoacylatedRNAGlobalIndexs));
            aminoacylations(idxs) = aminoacylations(idxs) + ...
                this.monomerTRNACounts' * states.monomerProductions;

            %t(m)RNA aminoacylation reactions
            reactions = this.reactionModificationMatrix * aminoacylations;

            %metabolic load of t(m)RNA aminoacylation reactions
            bmProd = max(0, -this.reactionStoichiometryMatrix) * reactions;
            byProd = max(0,  this.reactionStoichiometryMatrix) * reactions;

            %% enzymes
            %initialize
            minEnzExp = zeros(size(this.enzymeWholeCellModelIDs));
            maxEnzExp = Inf(size(this.enzymeWholeCellModelIDs));

            %tRNAs
            minEnzExp(this.enzymeIndexs_tRNAs) = 1.95 * this.monomerTRNACounts' * states.monomerProductions0;

            %tmRNA
            minEnzExp(this.enzymeIndexs_tmRNA) = 0;

            %tRNA synthetases, transferases
            [~, idxs] = ismember(this.aminoacylatedRNAtRNAIndexs, this.aminoacylatedRNAGlobalIndexs);
            aminoacylations = states.rnaProductions0(this.aminoacylatedRNAGlobalIndexs);
            aminoacylations(idxs)                = max(aminoacylations(idxs),                minEnzExp(this.enzymeIndexs_tRNAs));
            aminoacylations(setdiff(1:end,idxs)) = max(aminoacylations(setdiff(1:end,idxs)), minEnzExp(this.enzymeIndexs_tmRNA));

            minEnzExp(setdiff(1:end, [this.enzymeIndexs_tRNAs; this.enzymeIndexs_tmRNA])) = ...
                + minEnzExp(setdiff(1:end, [this.enzymeIndexs_tRNAs; this.enzymeIndexs_tmRNA])) ...
                + 3 * this.reactionCatalysisMatrix' * ((this.reactionModificationMatrix * aminoacylations) ./ this.enzymeBounds(:, 2));
        end

        %initialization: RNAs initialized to mature state by simulation
        %initializeState method. This class initializes RNAs to
        %aminoacylated state.
        function initializeState(this)
            totalRNAs = this.aminoacylatedRNAs + this.freeRNAs;
            this.aminoacylatedRNAs = ceil(2 / 3 * totalRNAs);
            this.freeRNAs = totalRNAs - this.aminoacylatedRNAs;
        end

        %resource requirements
        function result = calcResourceRequirements_Current(this)
            result = max(0, -this.reactionStoichiometryMatrix) * min(...
                this.reactionModificationMatrix * (this.freeRNAs + this.aminoacylatedRNAs + 1), ...
                this.reactionCatalysisMatrix * this.enzymes(1:numel(this.speciesIndexs_enzymes)) * this.stepSizeSec);
        end

        %simulation
        function evolveState(this)
            %terminate early if no free RNAs
            if ~any(this.freeRNAs)
                return;
            end

            numRNAs = numel(this.freeRNAs);

            species = [
                this.substrates;
                this.enzymes(1:numel(this.speciesIndexs_enzymes));
                this.freeRNAs]';

            reactionFluxes = zeros(size(this.speciesReactantByproductMatrix, 1), 1);
            anyFlux = false;

            limits = species(ones(numRNAs, 1), :) ./ this.speciesReactantMatrix;
            limits(:, this.substrateIndexs_water) = NaN;
            limits(:, this.substrateIndexs_hydrogen) = NaN;
            reactionLimits = min(limits, [], 2)';
            reactionLimits(isinf(reactionLimits) | isnan(reactionLimits) | reactionLimits < 0) = 0;
            isReactionInactive = reactionLimits <= 0;

            while true
                %compute maximum number of each species that can be modified
                limits = species(ones(numRNAs, 1), :) ./ max(0, this.speciesReactantByproductMatrix);
                limits(:, this.substrateIndexs_water) = NaN;
                limits(:, this.substrateIndexs_hydrogen) = NaN;
                reactionLimits = min(...
                    this.randStream.stochasticRound(min(limits(:, this.speciesIndexs_enzymes), [], 2)), ...
                    min(limits(:, [1:this.speciesIndexs_enzymes(1)-1  this.speciesIndexs_enzymes(end)+1:end]), [], 2))';
                reactionLimits(isReactionInactive | isinf(reactionLimits) | isnan(reactionLimits) | reactionLimits < 1) = 0;

                %stop if no more tRNA can be aminoacylated
                if ~any(reactionLimits)
                    break;
                end
                anyFlux = true;

                edges = min([0 cumsum(reactionLimits(:)' / sum(reactionLimits))], 1);
                nRxns = min(reactionLimits(reactionLimits > 0));
                if nRxns <= 1
                    %pick substrate, reaction, enzymes
                    [~, selectedReaction] = histc(rand(this.randStream, 1, 1), edges);
                    reactionFluxes(selectedReaction) = ...
                        + reactionFluxes(selectedReaction) ...
                        + 1;

                    %decrement metabolites, enzymes, unmodified substrates; increment modified substrates
                    species = species - this.speciesReactantByproductMatrix(selectedReaction, :);
                else
                    %pick substrate, reaction, enzymes
                    edges(end) = 1.1;
                    selectedReactions = histc(rand(this.randStream, nRxns, 1), edges);
                    selectedReactions = selectedReactions(1:end-1);
                    selectedReactions = selectedReactions(:);
                    reactionFluxes = ...
                        + reactionFluxes ...
                        + selectedReactions;

                    %decrement metabolites, enzymes, unmodified substrates; increment modified substrates
                    species = species - selectedReactions' * this.speciesReactantByproductMatrix;
                end
            end

            %stop if no reaction can proceed
            if ~anyFlux
                return;
            end

            % update substrates
            this.substrates = this.substrates + ...
                this.reactionStoichiometryMatrix * this.reactionModificationMatrix * reactionFluxes;

            % update RNAs
            this.freeRNAs          = this.freeRNAs          - reactionFluxes;
            this.aminoacylatedRNAs = this.aminoacylatedRNAs + reactionFluxes;
        end
    end

    %get methods of dependent global state
    methods
        function value = get.tRNASynthetases(this)
            value = this.enzymes(this.enzymeIndexs_tRNASynthetases,:,:);
        end

        function value = get.tRNATransferases(this)
            value = this.enzymes(this.enzymeIndexs_tRNATransferases,:,:);
        end
    end

    %get methods of dependent local state
    methods
        function value = getDryWeight(this)
            if size(this.freeRNAs, 3) == 1
                value = this.getDryWeight@edu.stanford.covert.cell.sim.ReactionProcess() + (...
                    this.rna.molecularWeights(this.rna.matureIndexs(this.aminoacylatedRNAGlobalIndexs))'        * this.freeRNAs         + ...
                    this.rna.molecularWeights(this.rna.aminoacylatedIndexs(this.aminoacylatedRNAGlobalIndexs))' * this.aminoacylatedRNAs) ...
                    / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            else
                value = this.getDryWeight@edu.stanford.covert.cell.sim.ReactionProcess() + (...
                    permute(this.rna.molecularWeights(this.rna.matureIndexs(this.aminoacylatedRNAGlobalIndexs))'        * permute(this.freeRNAs,         [1 3 2]),[1 3 2]) + ...
                    permute(this.rna.molecularWeights(this.rna.aminoacylatedIndexs(this.aminoacylatedRNAGlobalIndexs))' * permute(this.aminoacylatedRNAs,[1 3 2]),[1 3 2])) ...
                    / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function vl_test_ttlayers(gpu, tests)
% VL_TEST_TTLAYERS Test the TT-layer with numeric differentiation
%    VL_TEST_TTLAYERS(0) Test the CPU implementation.
%    VL_TEST_TTLAYERS(1) Test the GPU implementation.

range = 100;

if nargin < 1, gpu = false ; end
if gpu
  grandn = @(varargin) range * gpuArray.randn(varargin{:});
  grand = @(varargin) range * gpuArray.rand(varargin{:});
else
  grandn = @(varargin) range * randn(varargin{:});
  grand = @(varargin) range * rand(varargin{:});
end

switch gpu
  case 0,
    fprintf('testing the CPU code\n');
  case 1
    fprintf('testing the GPU code\n');
end

rng(1);

if nargin < 2
  tests = 1:3;
end

function y = vl_nntt_forward_weights(layer, in, out, iGroup, values)
  layer.weights{iGroup} = values;
  outIn = vl_nntt_forward(layer, in, out);
  y = outIn.x;
end

function y = vl_nntt_forward_x(layer, in, out, x)
  in.x = x;
  outIn = vl_nntt_forward(layer, in, out);
  y = outIn.x;
end

for l = tests
  fprintf('test number %d\n', l)
  % resets random number generator to obtain reproducible results
  if gpu
    parallel.gpu.rng(0, 'combRecursive');
  else
    rng(0, 'combRecursive');
  end
  switch l
    case 1
      disp('Testing vl_nntt_* with the identity TT-matrix.');

      in.x = grandn(8, 32, 3, 4, 'single');
      W = tt_ones([4, 4, 4, 4, 3]);
      W.core = single(W.core);
      W = diag(W);
      layer.W = W;
      layer.weights{1} = W.core;
      if gpu
        layer.weights{1} = gpuArray(layer.weights{1});
      end
      layer.outHeight = 8;
      layer.outWidth = 32;
      layer.outChannels = 3;
      layer.weights{2} = grandn(8 * 32 * 3, 1, 'single');
      out = [];
      out = vl_nntt_forward(layer, in, out);
      y = out.x;
      out.dzdx = grandn(size(y), 'single');
      in = vl_nntt_backward(layer, in, out);
      for iGroup = 1:numel(layer.weights)
          vl_testder(@(w) vl_nntt_forward_weights(layer, in, out, iGroup, w), layer.weights{iGroup}, out.dzdx, in.dzdw{iGroup}, range * 1e-2);
      end
      vl_testder(@(x) vl_nntt_forward_x(layer, in, out, x), in.x, out.dzdx, in.dzdx, range * 1e-2);

    case 2
      disp('Testing vl_nntt_* with a random square TT-matrix.');
      % Shape for the input and output tensors.
      tensorShape = [4, 4, 4, 4, 3];
      batchSize = 10;
      ranks = [1, 4, 6, 10, 5, 1];
      W = tt_rand(tensorShape.^2, 5, ranks);
      W.core = single(W.core);
      W = tt_matrix(W, tensorShape, tensorShape);
      layer.W = W;
      layer.weights{1} = W.core;
      if gpu
        layer.weights{1} = gpuArray(layer.weights{1});
      end
      layer.outHeight = 8;
      layer.outWidth = 32;
      layer.outChannels = 3;
      in.x = grandn(8, 32, 3, batchSize, 'single');
      layer.weights{2} = grandn(8 * 32 * 3, 1, 'single');
      out = [];
      out = vl_nntt_forward(layer, in, out);
      y = out.x;
      exactY = full(W) * reshape(in.x, [], batchSize);
      exactY = bsxfun(@plus, exactY, layer.weights{2});
      vl_testsim(y, reshape(exactY, 8, 32, 3, batchSize));
      out.dzdx = grandn(size(y), 'single');
      in = vl_nntt_backward(layer, in, out);
      for iGroup = 1:numel(layer.weights)
          vl_testder(@(w) vl_nntt_forward_weights(layer, in, out, iGroup, w), layer.weights{iGroup}, out.dzdx, in.dzdw{iGroup}, range * 1e-2);
      end
      vl_testder(@(x) vl_nntt_forward_x(layer, in, out, x), in.x, out.dzdx, in.dzdx, range * 1e-2);

    case 3
      disp('Testing vl_nntt_* with random rectangular TT-matrices.');
      for bias = [false true]
        for batchSize = [1 3]
          inputTensorShape = [3, 6, 4, 5];
          outputTensorShape = [4, 11, 7, 13];
          layer.outHeight = 2 * 11;
          layer.outWidth = 7 * 13;
          layer.outChannels = 2;
          ranks = [1, 5, 9, 5, 1];
          W = tt_rand(outputTensorShape .* inputTensorShape, length(inputTensorShape), ranks, []);
          W.core = single(W.core);
          W = tt_matrix(W, outputTensorShape, inputTensorShape);
          layer.W = W;
          layer.weights{1} = W.core;
          if gpu
            layer.weights{1} = gpuArray(layer.weights{1});
          end
          if bias
            layer.weights{2} = grandn(prod(outputTensorShape), 1, 'single');
          else
            layer.weights{2} = [];
          end
          in.x = grandn(9, 8, 5, batchSize, 'single');
          out = [];
          out = vl_nntt_forward(layer, in, out);
          y = out.x;
          exactY = full(W) * reshape(in.x, [], batchSize);
          if bias
            exactY = bsxfun(@plus, exactY, layer.weights{2});
          end
          vl_testsim(y, reshape(exactY, layer.outHeight, layer.outWidth, layer.outChannels, batchSize));
          out.dzdx = grandn(size(y), 'single');
          in = vl_nntt_backward(layer, in, out);
          for iGroup = 1:numel(layer.weights)
              vl_testder(@(w) vl_nntt_forward_weights(layer, in, out, iGroup, w), layer.weights{iGroup}, out.dzdx, in.dzdw{iGroup}, range * 1e-2);
          end
          vl_testder(@(x) vl_nntt_forward_x(layer, in, out, x), in.x, out.dzdx, in.dzdx, range * 1e-2);
        end
      end
    end
  end
end

author: Bihaqo

repo: TensorNet

save_path: github-repos/MATLAB/Bihaqo-TensorNet


text:
function jed = epoch_to_jed_saka ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_SAKA returns the epoch of the Saka calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1749994.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [min_val max_val] = GradDescLimits(modelname)
%
% function [min_val max_val] = GradDescLimits(modelname)
%
% Returns maximum and minimum settings for the parameters of different
% models to use during direct fitting.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

LARGE = 100000;

D_MIN = 0.001;
D_MAX = 3;
D_PERP_MIN = 0.001;

ANGLEMAX=100;

R_MIN = 0.1;
R_MAX = 5;

% the values need to account for scaling defined
% in GetScalingFactors.m
K_MIN = 0;
K_MAX = 6.4;

B_MIN = 0;
B_MAX = 3.2;

parameterStrings = GetParameterStrings(modelname);
min_val = zeros(1, length(parameterStrings));
max_val = zeros(1, length(parameterStrings));

for i=1:length(parameterStrings)
    if (strcmp(parameterStrings(i), 'ficvf') ||...
        strcmp(parameterStrings(i), 'fiso') ||...
        strcmp(parameterStrings(i), 'irfrac'))
        min_val(i) = 0;
        max_val(i) = LARGE;
    elseif (strcmp(parameterStrings(i), 'di') ||...
            strcmp(parameterStrings(i), 'diso'))
        min_val(i) = sqrt(D_MIN);
        max_val(i) = sqrt(D_MAX);
    elseif (strcmp(parameterStrings(i), 'dh'))
        min_val(i) = sqrt(D_PERP_MIN);
        max_val(i) = LARGE;
    elseif (strcmp(parameterStrings(i), 'rad'))
        min_val(i) = sqrt(R_MIN);
        max_val(i) = sqrt(R_MAX);
    elseif (strcmp(parameterStrings(i), 'kappa'))
        min_val(i) = sqrt(K_MIN);
        max_val(i) = sqrt(K_MAX);
    elseif (strcmp(parameterStrings(i), 'beta'))
        min_val(i) = sqrt(B_MIN);
        max_val(i) = sqrt(B_MAX);
    elseif (strcmp(parameterStrings(i), 'theta') ||...
            strcmp(parameterStrings(i), 'phi') ||...
            strcmp(parameterStrings(i), 'psi'))
        min_val(i) = -ANGLEMAX;
        max_val(i) = ANGLEMAX;
    elseif (strcmp(parameterStrings(i), 'b0'))
        min_val(i) = 0.001;
        max_val(i) = LARGE;
    elseif (strcmp(parameterStrings(i), 't1'))
        min_val(i) = 0.1;
        max_val(i) = 1;
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function aux_parkfield(vResult, hParentFigure)

hSets = gui_result2('GetDatasetsHandle', hParentFigure, [], guidata(hParentFigure));
sList = get(hSets, 'String');
nLen = length(sList);
vResults = gui_result2('GetResults', hParentFigure, [], guidata(hParentFigure));

nCS=[[-154.01 60.88 -150.50 60.37];
    [-153.45 61.68 -149.81 61.21]];

for nCnt = 1:size(nCS,1)  % nLen

    vResults(nLen+nCnt)=aux_cs2(vResults(nCnt),nCS(nCnt,:));

    hPlot = gui_result2('GetFigureHandle', hParentFigure, [], guidata(hParentFigure));
    %   exportfig(hPlot, [num2str(nCnt) '.eps'], 'Color', 'cmyk');

end

save vResults_CS.mat vResults -mat

disp('Calc of CS finished: vResults_CS.mat')


text:
function F = Gaus_filter(I,sigma)

R = ceil(3*sigma);
for i = -R:R,
    for j = -R:R,
        M(i+ R+1,j+R+1) = exp(-(i*i+j*j)/2/sigma/sigma)/(2*pi*sigma*sigma);
    end
end
M = M/sum(sum(M));  % normalize
F = xconv2(I,M);


text:
function [thrs, fars, frrs] = slhistroc(hist_a, hist_r, sepvals, op)
%SLHISTROC Computes the ROC curve from value histogram
%
% $ Syntax $
%   - [thrs, fas, frs] = slhistroc(hist_a, hist_r, sepvals, op)
%
% $ Arguments $
%   - hist_a:       the histogram of the values that should be accepted
%   - hist_r:       the histogram of the values that should be rejected
%   - sepvals:      the separation values of the histograms
%   - op:           the option of the attributes of the values
%   - thrs:         the sampled threshold values
%   - fars:         the false accept rates at the sampled thresholds
%   - frrs:         the false reject rates at the sampled thresholds
%
% $ Description $
%   - [thrs, fas, frs] = slhistroc(hist_a, hist_r, sepvals, op) computes
%     the ROC curve based on histograms. The hist_a and hist_r should have
%     corresponding bins, and the bins are sorted in ascending order of
%     the values they represent. Suppose the number of bins be n, then
%     sepvals are the values on the bin edges, and the length of sepvals
%     should be n. Note that the histtograms should be created following
%     the rules as in histc, so the last element of histogram is the
%     number of elements that exactly match the last sep value.
%     The op can be either 'low' or 'high', when it is 'low', the value
%     lower than threshold would be accepted, otherwise the value higher
%     than threshold would be accepted.
%     The output thrs is just sepvals, while fars and frrs are
%     corresponding false accept rates and false reject rates, with
%     n = nbins + 1 elements.
%
% $ History $
%   - Created by Dahua Lin, on Aug 8th, 2006
%

%% parse and verify arguments

if nargin < 4
    raise_lackinput('slhistroc', 4);
end

nbins = length(hist_a);
if length(hist_r) ~= nbins
    error('sltoolbox:sizmismatch', ...
        'The sizes of hist_a and hist_r are inconsistent');
end
if length(sepvals) ~= nbins
    error('sltoolbox:sizmismatch', ...
        'The length of sepvals should be number of bins');
end

%% Compute ROC

% preprocess
hist_r(end-1) = hist_r(end-1) + hist_r(end);
hist_a(end-1) = hist_a(end-1) + hist_a(end);
hist_r = hist_r(1:end-1);
hist_a = hist_a(1:end-1);

hist_r = hist_r(:);
hist_a = hist_a(:);
thrs = sepvals(:);

nr = sum(hist_r);
na = sum(hist_a);

switch op
    case 'low'
        fars = [0; cumsum(hist_r)] / nr;
        frrs = [na; na - cumsum(hist_a)] / na;

    case 'high'
        fars = [nr; nr - cumsum(hist_r)] / nr;
        frrs = [0; cumsum(hist_a)] / na;

    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid option %s for roc', op);
end

author: lmthang


text:
function y = tfocs_round(x)
if iscell( x ),
    y = cellfun( @tfocs_round, x, 'UniformOutput', false );
else
    y = round( x );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [rbin]=rdfcalc(C,xlim,ylim)

%This file takes a file containing the co-ordinates of the centers as
%input, and finds the rdf

%xlim=512; % The maximum x co-ordinate of the image
%ylim=512; % The maximum y co-ordinate of the image

n=size(C,1);

rmin=1;
rmax=100;
dr=2;

rvecarray=rmin:dr:rmax;

m=size(rvecarray,2);

rbin=zeros(m,2);
rbin(:,1)=rvecarray;

for j=1:m

    rvec=rvecarray(j);

for i=1:n
    x1=C(i,1);
    y1=C(i,2);
    bini=0;

    for k=1:n
        rk= sqrt((C(k,1)-C(i,1))^2 +(C(k,2)-C(i,2))^2);
        if rk>= rvec && rk<=(rvec+dr)
            bini=bini+1;
        end
    end

    theta1=checkquad(x1,y1,rvec,1,xlim,ylim);
    theta2=checkquad(x1,y1,rvec,2,xlim,ylim);
    theta3=checkquad(x1,y1,rvec,3,xlim,ylim);
    theta4=checkquad(x1,y1,rvec,4,xlim,ylim);

    theta=theta1+theta2+theta3+theta4;

    area1=rvec*theta*dr;

    bini=bini/area1;
    rbin(j,2)=rbin(j,2)+bini;

end
end

% rbin=rbin/n; % This is to divide the result by the number of particles.
%
% numden=(n/(xlim*ylim)); %This is the number density
%
% rbin=rbin/numden;
%
% plot(rvecarray,rbin,'ro','MarkerSize',10)


text:
% predicting numerical values using Linear Regression
clear all;
format long
disp('===== Linear Regression ====');
disp('Reading featur vector');

for feat = 1:3
    featurs = csvread('data\forWeka_featuresonly.csv');
    num_data = size(featurs,1); %5000;
    disp(sprintf('Number of datapoints %d',num_data))

    possiblefeaturizations =  {'bernouli', 'tfidf','multinomial'};
    %featurization = 'bernouli'%'tfidf'%'tfidf'%'multinomial'%'tfidf' %'multinomial'; % 'bernouli', 'tfidf'
    featurization  = possiblefeaturizations{feat}

    featurs = featurs(:,2:size(featurs,2));
    if strcmp(featurization,'multinomial')
        %just pass
    elseif strcmp(featurization,'bernouli')
        featurs = bernoulli(featurs);
    elseif strcmp(featurization,'tfidf')
        featurs = tfidf(featurs);
    end

    size_training = floor(.8*num_data);

    trainingset = featurs(1:size_training,:);
    testset = featurs((size_training+1):num_data,:);

    disp('Splitting up data into training/test sets');
    [num,txt,raw] = xlsread('data\final104.xls');

    % reading the description of each shoe
    descriptions = raw(2:size(raw,1),2);
    style_ratings = num(1:size(num,1),1);
    comfort_ratings = num(1:size(num,1),4);
    overal_ratings = num(1:size(num,1),5);

    % only take m data points
    m=num_data;
    descriptions = descriptions(1:m);
    style_ratings = style_ratings(1:m);
    comfort_ratings = comfort_ratings(1:m);
    overal_ratings = overal_ratings(1:m);

    responsevals = [style_ratings, comfort_ratings, overal_ratings];

    responsevals_training = responsevals(1:size_training,:);
    responsevals_test = responsevals((size_training+1):num_data,:);

    disp('Discreminate Analsys');
    %

    tic;

    predictions = [];
    actual = [];
    for i =1:3
        a = responsevals_training(:,i);
        b = responsevals_test(:,i)';
        %regresscoeff = regress(a, trainingset);
        %C2 = (regresscoeff'*(testset'));
        class = classify(testset,trainingset,a)
        predictions = [predictions, class];
        actual = [actual, b'];
    end

    MSE = mean(sum(((predictions-actual).^2)'))
    toc;

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [ x, y, z ] = image_fun ( )

%*****************************************************************************80
%
%% IMAGE_FUN processes an image using SPMD.
%
%  Discussion:
%
%    Only 3 SPMD workers should be allocated, one each for the
%    R, G, and B components of the image.
%
%    The image "balloons.tif" is read in by the client, and distributed.
%    This RGB image is stored as a (:,:,3) array, and by default, is
%    distributed by the last dimension.  That means the first three workers
%    get the R, G and B arrays and if there are more workers they get nothing.
%
%    Each worker applies the MEDFILT2 operation to its data.
%
%    The client then assembles the filtered R, G, and B data back into
%    a (:,:,3) array and returns that as the output argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, uint8 X(:,:,3), the RGB data for the original image.
%
%    Output, uint8 Y(:,:,3), the RGB data for the noisy image.
%
%    Output, uint8 Z(:,:,3), the RGB data for the filtered image.
%

%
%  Read an image X.
%  This happens to be a color image, and is stored as 480x640x3 array.
%
  x = imread ( 'balloons.tif' );
%
%  Create an image Y by adding "salt and pepper" noise to X.
%
  y = imnoise ( x, 'salt & pepper', 0.30 );
%
%  Make YD, a distributed version of Y, by copying a separate
%  portion of Y to each worker.
%
  yd = distributed ( y );
%
%  Each worker creates YL, a name for its portion of the array.
%  It then applies the median filter to YL, using a 3x3 block of
%  data around each pixel.
%
  spmd
    yl = getLocalPart ( yd );
    yl = medfilt2 ( yl, [ 3 3 ] );
  end
%
%  The client retrieves the data from each worker.
%  Assuming three workers were used, and the data was divided by
%  the last dimension (R,G,B), then the three pieces are
%  tacked together as follows:
%
  z        = yl{1};
  z(:,:,2) = yl{2};
  z(:,:,3) = yl{3};

  return
end


text:
function tnet = vl_simplenn_tidy(net)
%VL_SIMPLENN_TIDY  Fix an incomplete or outdated SimpleNN network.
%   NET = VL_SIMPLENN_TIDY(NET) takes the NET object and upgrades
%   it to the current version of MatConvNet. This is necessary in
%   order to allow MatConvNet to evolve, while maintaining the NET
%   objects clean.
%
%   The function is also generally useful to fill in missing default
%   values in NET.
%
%   See also: VL_SIMPLENN().

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

tnet = struct('layers', {{}}, 'meta', struct()) ;

% copy meta information in net.meta subfield
if isfield(net, 'classes')
  tnet.meta.classes = net.classes ;
end

if isfield(net, 'normalization')
  tnet.meta.normalization = net.normalization ;
end

if isfield(net, 'meta')
  tnet.meta = net.meta ;
end

% copy layers
for l = 1:numel(net.layers)
  defaults = {'precious', false};
  layer = net.layers{l} ;

  % check weights format
  switch layer.type
    case {'conv', 'convt', 'bnorm'}
      if ~isfield(layer, 'weights')
        layer.weights = {...
          layer.filters, ...
          layer.biases} ;
        layer = rmfield(layer, 'filters') ;
        layer = rmfield(layer, 'biases') ;
      end
  end

  % check that weights inlcude moments in batch normalization
  if strcmp(layer.type, 'bnorm')
    if numel(layer.weights) < 3
      layer.weights{3} = ....
        zeros(numel(layer.weights{1}),2,'single') ;
    end
  end

  % fill in missing values
  switch layer.type
    case {'conv', 'pool'}
      defaults = [ defaults {...
        'pad', 0, ...
        'stride', 1, ...
        'opts', {}}] ;

    case 'convt'
      defaults = [ defaults {...
        'crop', 0, ...
        'upsample', 1, ...
        'numGroups', 1, ...
        'opts', {}}] ;

    case {'pool'}
      defaults = [ defaults {...
        'method', 'max', ...
        'pad', 0, ...
        'stride', 1, ...
        'opts', {}}] ;

    case 'relu'
      defaults = [ defaults {...
        'leak', 0}] ;

    case 'dropout'
      defaults = [ defaults {...
        'rate', 0.5}] ;

    case {'normalize', 'lrn'}
      defaults = [ defaults {...
        'param', [5 1 0.0001/5 0.75]}] ;

    case {'pdist'}
      defaults = [ defaults {...
        'noRoot', false, ...
        'aggregate', false, ...
        'p', 2, ...
        'epsilon', 1e-3} ];
  end

  for i = 1:2:numel(defaults)
    if ~isfield(layer, defaults{i})
      layer.(defaults{i}) = defaults{i+1} ;
    end
  end

  % save back
  tnet.layers{l} = layer ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function w2 = waveform_remove_empty(w);
%WAVEFORM_REMOVE_EMPTY remove empty waveform objects from a vector of waveform objects
%   Note that the opposite function, of expanding empty waveform objects can largely
%   be achieved with waveform/pad

%   Glenn Thompson
    e = 1;
    for c=1:length(w)
            nsamp = get(w(c), 'data_length');
            if nsamp > 1 % 0 means blank waveform object, 1 means corrupt waveform object
                    w2(e) = w(c);
                    e = e + 1;
            end
    end
    if e==1
            w2 = [];
    end
end

author: geoscience-community-codes

repo: GISMO


text:
s = tf([1 0],1);
Fi = (s-1)/(s^2+3*s+2);
Fo = 1/(s+1);

F = Fi*Fo;
step(F)
bode(F)
% clear all
% t = 0:.001:100;
% thetastar(:,1) = t;
% thetastar(:,2) = .01*exp(.01*t);
% fstar(:,1) = t;
% fstar(:,2) = [(0:.001:10)*0 (10.001:.001:100)*0+.01];

author: dynamicslab