text:
function varargout = rfDesignPanel(varargin)
% RFDESIGNPANEL M-file for rfDesignPanel.fig
%      RFDESIGNPANEL, by itself, creates a new RFDESIGNPANEL or raises the existing
%      singleton*.
%
%      H = RFDESIGNPANEL returns the handle to a new RFDESIGNPANEL or the handle to
%      the existing singleton*.
%
%      RFDESIGNPANEL('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in RFDESIGNPANEL.M with the given input arguments.
%
%      RFDESIGNPANEL('Property','Value',...) creates a new RFDESIGNPANEL or
%      raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before rfDesignPanel_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to rfDesignPanel_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above Gz_text to modify the response to help rfDesignPanel

% Last Modified by GUIDE v2.5 20-Jan-2014 15:46:34

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @rfDesignPanel_OpeningFcn, ...
    'gui_OutputFcn',  @rfDesignPanel_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before rfDesignPanel is made visible.
function rfDesignPanel_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to rfDesignPanel (see VARARGIN)

handles.Simuh=varargin{1};
% Load pulse element list
try
    handles.SeqElemListStruct=DoParseXML([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'SeqElem.xml']);
    guidata(hObject,handles);
    Root=DoConvStruct2Tree(handles.SeqElemListStruct);
    handles.SeqElemListTree=uitree('v0','Root',Root,'SelectionChangeFcn', {@ChkSeqElemListTreeNode,handles});
    set(handles.SeqElemListTree.getUIContainer,'Units','normalized');
    set(handles.SeqElemListTree.getUIContainer,'Position',[0.0,0.2,0.125,0.8]);
    handles.Attrh0=[];
    handles.Attrh1=[];
    handles.Attrh2=[];
    handles.Attrh3=[];
    handles.Attrh4=[];
catch ME
    errordlg('SeqElem.xml file is missing or can not be loaded!');
    close(hObject);
    return;
end

handles.SeqAttri_tabgroup=uitabgroup(handles.SeqAttri_uipanel);
handles.rf_tab=uitab( handles.SeqAttri_tabgroup, 'title', 'rf');
handles.Gz_tab=uitab( handles.SeqAttri_tabgroup, 'title', 'Gz');
handles.Gy_tab=uitab( handles.SeqAttri_tabgroup, 'title', 'Gy');
handles.Gx_tab=uitab( handles.SeqAttri_tabgroup, 'title', 'Gx');

%Load tabs for Spin Property & Gradient
try
    handles.rfSimuAttrStruct=DoParseXML([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'rfSimuAttr.xml']);
catch ME
    errordlg('rfSimuAttr.xml file is missing or can not be loaded!');
    close(hObject);
    return;
end

handles.Setting_tabgroup=uitabgroup(handles.Setting_uipanel);
guidata(handles.rfDesignPanel_figure,handles);
for i=1:length(handles.rfSimuAttrStruct.Children)
    eval(['handles.' handles.rfSimuAttrStruct.Children(i).Name '_tab=uitab( handles.Setting_tabgroup,' '''title'',' '''' handles.rfSimuAttrStruct.Children(i).Name ''',''Units'',''normalized'');']);
    eval(['DoEditValue(handles,handles.' handles.rfSimuAttrStruct.Children(i).Name '_tab,handles.rfSimuAttrStruct.Children(' num2str(i) ').Attributes,0,[0.25,0.25,0.0,0.08,0.08,0.02]);']);
    handles=guidata(handles.rfDesignPanel_figure);
end

%Create tabs for Slice Profile
handles.SliProfile_tabgroup=uitabgroup(handles.SliProfile_uipanel);
handles.MxMyMz_tab=uitab(handles.SliProfile_tabgroup,'title','Mx My Mz','Units','normalized');
handles.MxyMz_tab=uitab(handles.SliProfile_tabgroup,'title','|Mxy| Mz','Units','normalized');
handles.MP_tab=uitab(handles.SliProfile_tabgroup,'title','Mg Pe','Units','normalized');
handles.MxMyMz_axes=axes('parent', handles.MxMyMz_tab,'Position', [0.06 0.1 0.88 0.85],'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);
handles.MxyMz_axes=axes('parent', handles.MxyMz_tab,'Position', [0.06 0.1 0.88 0.85],'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);
handles.MP_axes=axes('parent', handles.MP_tab,'Position', [0.06 0.1 0.88 0.85],'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);

%Create tabs for XYZ Gradient
handles.G_tabgroup=uitabgroup(handles.Gradient_uipanel);
handles.Gradz_tab=uitab(handles.G_tabgroup,'title','Gz (G/cm)','Units','normalized');
handles.Grady_tab=uitab(handles.G_tabgroup,'title','Gy (G/cm)','Units','normalized');
handles.Gradx_tab=uitab(handles.G_tabgroup,'title','Gx (G/cm)','Units','normalized');
handles.Gz_axes=axes('parent', handles.Gradz_tab,'Position', [0.06 0.1 0.88 0.85],'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);
handles.Gy_axes=axes('parent', handles.Grady_tab,'Position', [0.06 0.1 0.88 0.85],'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);
handles.Gx_axes=axes('parent', handles.Gradx_tab,'Position', [0.06 0.1 0.88 0.85],'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);

% Choose default command line output for rfDesignPanel
handles.output=hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes rfDesignPanel wait for user response (see UIRESUME)
% uiwait(handles.rfDesignPanel_figure);

function ChkSeqElemListTreeNode(tree, value, handles)

handles=guidata(handles.rfDesignPanel_figure);
SelNode=tree.SelectedNodes;
SelNode=SelNode(1);
Level=SelNode.getValue;
if ~ischar(Level)
    LevelChg=diff(Level,1);
    Level([LevelChg; -1]==1)=[];
end

Node='handles.SeqElemListStruct';
if Level=='0'
    eval(['handles.SeqElemNode=' Node ';']);
else
    for i=1:length(Level)
        Node=[Node '.Children(' num2str(Level(i)) ')'];
    end
    eval(['handles.SeqElemNode=' Node ';']);
end

set(handles.SeqAttri_uipanel,'Title',['Checking ''' char(SelNode.getName) '''']);
if SelNode.isLeaf | SelNode.isRoot

    if strcmp(handles.SeqElemNode.Name(1:2),'rf') & length(handles.SeqElemNode.Name)>2

        if ~isempty(handles.Attrh1)
            Attrh1name=fieldnames(handles.Attrh1);
            for i=1:length(Attrh1name)
                delete(handles.Attrh1.(Attrh1name{i}));
            end
            handles.Attrh1=[];
        end
        handles.rfNode=handles.SeqElemNode;
        handles.rfNodeLvl=Node;
        set(handles.rf_tab,'Title',['rf:' handles.SeqElemNode.Name]);
        DoEditValue(handles,handles.rf_tab,handles.SeqElemNode.Attributes,1,[0.50,0.50,0.0,0.05,0.05,0.005]);

        %----load rf memo
        fid=fopen([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'rf' filesep handles.SeqElemNode.Name '_Memo.txt'],'r');
        if fid==-1
            set(handles.rfMemo_edit,'String','rf Memo is not available!');
        else
            tline = fgetl(fid);
            i=1;
            while ischar(tline)
                Memo{i,1}=tline;
                tline = fgetl(fid);
                i=i+1;
            end
            if i==1
                set(handles.rfMemo_edit,'String','rf Memo is empty!');
            else
                set(handles.rfMemo_edit,'String',Memo);
            end
            fclose(fid);
        end
        %----end
        set(handles.UpdateAttrXML_pushbutton,'Enable','on');
     elseif strcmp(handles.SeqElemNode.Name(1:2),'Gz') & length(handles.SeqElemNode.Name)>4

        if ~isempty(handles.Attrh2)
            Attrh2name=fieldnames(handles.Attrh2);
            for i=1:length(Attrh2name)
                delete(handles.Attrh2.(Attrh2name{i}));
            end
            handles.Attrh2=[];
        end
        handles.GzNode=handles.SeqElemNode;
        handles.GzNodeLvl=Node;
        set(handles.Gz_tab,'Title',['Gz:' handles.SeqElemNode.Name]);
        DoEditValue(handles,handles.Gz_tab,handles.SeqElemNode.Attributes,2,[0.50,0.50,0.0,0.05,0.05,0.015]);
        set(handles.UpdateAttrXML_pushbutton,'Enable','on');
     elseif strcmp(handles.SeqElemNode.Name(1:2),'Gy') & length(handles.SeqElemNode.Name)>4

        if ~isempty(handles.Attrh3)
            Attrh3name=fieldnames(handles.Attrh3);
            for i=1:length(Attrh3name)
                delete(handles.Attrh3.(Attrh3name{i}));
            end
            handles.Attrh3=[];
        end
        handles.GyNode=handles.SeqElemNode;
        handles.GyNodeLvl=Node;
        set(handles.Gy_tab,'Title',['Gy:' handles.SeqElemNode.Name]);
        DoEditValue(handles,handles.Gy_tab,handles.SeqElemNode.Attributes,3,[0.50,0.50,0.0,0.05,0.05,0.015]);
        set(handles.UpdateAttrXML_pushbutton,'Enable','on');
     elseif strcmp(handles.SeqElemNode.Name(1:2),'Gx') & length(handles.SeqElemNode.Name)>4

        if ~isempty(handles.Attrh4)
            Attrh4name=fieldnames(handles.Attrh4);
            for i=1:length(Attrh4name)
                delete(handles.Attrh4.(Attrh4name{i}));
            end
            handles.Attrh4=[];
        end
        handles.GxNode=handles.SeqElemNode;
        handles.GxNodeLvl=Node;
        set(handles.Gx_tab,'Title',['Gx:' handles.SeqElemNode.Name]);
        DoEditValue(handles,handles.Gx_tab,handles.SeqElemNode.Attributes,4,[0.50,0.50,0.0,0.05,0.05,0.015]);
        set(handles.UpdateAttrXML_pushbutton,'Enable','on');
    end
else
    guidata(handles.rfDesignPanel_figure, handles);
end

% --- Outputs from this function are returned to the command line.
function varargout = rfDesignPanel_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = hObject;

% --- Executes on slider movement.
function rf_slider_Callback(hObject, eventdata, handles)
% hObject    handle to rf_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% set(handles.rf_line,'XData',[get(hObject,'Value'),get(hObject,'Value')]','YData',[min(handles.rfAmp) max(handles.rfAmp)]');

rf_step=round(get(hObject,'Value'));
[az,el] = view(handles.SpinRot_axes);

if strcmp(handles.Freq_Flag,'off') & strcmp(handles.Spat_Flag,'off')

    Mx=handles.Mx(:,:,:,:,:,rf_step);
    My=handles.My(:,:,:,:,:,rf_step);
    Mz=handles.Mz(:,:,:,:,:,rf_step);
    Mxy=sqrt(Mx.^2+My.^2);
    Pe=angle(Mx+1i.*My);
    for j=1:str2double(get(handles.Attrh0.TypeNum,'String'))
        spinver=quiver3(handles.SpinRot_axes,handles.Gxgrid,handles.Gygrid,handles.Gzgrid,Mx(:,:,:,:,j),My(:,:,:,:,j),Mz(:,:,:,:,j));
        set(spinver,'AutoScale','off');
        hold(handles.SpinRot_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mx(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MxMyMz_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(My(:,:,:,:,j)),'Color',[0 1-(j-1)*0.25 0],'linewidth',2);
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MxyMz_axes,'on');
        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MP_axes,'on');
        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Pe(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);
    end
    hold(handles.SpinRot_axes,'off');
    hold(handles.MxMyMz_axes,'off');
    hold(handles.MxyMz_axes,'off');
    hold(handles.MP_axes,'off');

    set(handles.SpinRot_axes,'YLim',[-1 1],'XLim',[-1 1],'ZLim',[min(handles.Gzgrid,[],3)-1 max(handles.Gzgrid,[],3)+1]);
    xlabel(handles.SpinRot_axes,'X axis');
    ylabel(handles.SpinRot_axes,'Y axis');
    zlabel(handles.SpinRot_axes,'Z axis (Gradient Direction)');
    set(handles.SpinRot_axes,'view',[az,el]);

    set(handles.MxMyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxMyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxMyMz_axes,'Mx','My','Mz');

    set(handles.MxyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxyMz_axes,'|Mxy|','Mz');

    set(handles.MP_axes,'YLim',[-pi pi],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MP_axes,'YGrid','on','XGrid','on');
    legend(handles.MP_axes,'Mg','Pe');

else
    updateSS(handles, rf_step);
end

plot(handles.rf_axes,handles.rfTime,handles.rfAmp,...
    [handles.Muts(rf_step) handles.Muts(rf_step)],[min(handles.rfAmp) max(handles.rfAmp)],'g-','linewidth',2);
set(handles.rf_axes,'YLim',[min(handles.rfAmp) max(handles.rfAmp)],'XLim',[min(handles.rfTime) max(handles.rfTime)]);
set(handles.rf_axes,'YGrid','on','XGrid','on');
%set(handles.SpinRot_axes,'ZTickLabel',get(handles.SpinRot_axes,'ZTick')*10);
%set(handles.MxMyMz_axes,'XTickLabel',get(handles.MxMyMz_axes,'XTick')*100);
set(handles.rf_axes,'XTickLabel',get(handles.rf_axes,'XTick')*1000);
set(handles.rf_axes,'YTickLabel',get(handles.rf_axes,'YTick')*10000);
set(handles.rfFreq_axes,'YTickLabel',get(handles.rfFreq_axes,'YTick')/1000);

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

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

[az,el] = view(handles.SpinRot_axes);
Starti=round(get(handles.rf_slider,'Value'));
for i=Starti+1:max(size(handles.Muts))
    handles=guidata(hObject);
    if handles.PauseFlag==1
        return;
    end
    set(handles.rf_slider,'Value',i);
    rf_step=round(get(handles.rf_slider,'Value'));

    if strcmp(handles.Freq_Flag,'off') & strcmp(handles.Spat_Flag,'off')

        Mx=handles.Mx(:,:,:,:,:,rf_step);
        My=handles.My(:,:,:,:,:,rf_step);
        Mz=handles.Mz(:,:,:,:,:,rf_step);
        Mxy=sqrt(Mx.^2+My.^2);
        Pe=angle(Mx+1i.*My);
        for j=1:str2double(get(handles.Attrh0.TypeNum,'String'))
            spinver=quiver3(handles.SpinRot_axes,handles.Gxgrid,handles.Gygrid,handles.Gzgrid,Mx(:,:,:,:,j),My(:,:,:,:,j),Mz(:,:,:,:,j));
            set(spinver,'AutoScale','off');
            hold(handles.SpinRot_axes,'on');
            plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mx(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
            hold(handles.MxMyMz_axes,'on');
            plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(My(:,:,:,:,j)),'Color',[0 1-(j-1)*0.25 0],'linewidth',2);
            plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

            plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
            hold(handles.MxyMz_axes,'on');
            plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

            plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
            hold(handles.MP_axes,'on');
            plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Pe(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);
        end
        hold(handles.SpinRot_axes,'off');
        hold(handles.MxMyMz_axes,'off');
        hold(handles.MxyMz_axes,'off');
        hold(handles.MP_axes,'off');

        set(handles.SpinRot_axes,'YLim',[-1 1],'XLim',[-1 1],'ZLim',[min(handles.Gzgrid,[],3)-1 max(handles.Gzgrid,[],3)+1]);
        xlabel(handles.SpinRot_axes,'X axis');
        ylabel(handles.SpinRot_axes,'Y axis');
        zlabel(handles.SpinRot_axes,'Z axis (Gradient Direction)');
        set(handles.SpinRot_axes,'view',[az,el]);

        set(handles.MxMyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
        set(handles.MxMyMz_axes,'YGrid','on','XGrid','on');
        legend(handles.MxMyMz_axes,'Mx','My','Mz');

        set(handles.MxyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
        set(handles.MxyMz_axes,'YGrid','on','XGrid','on');
        legend(handles.MxyMz_axes,'|Mxy|','Mz');

        set(handles.MP_axes,'YLim',[-pi pi],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
        set(handles.MP_axes,'YGrid','on','XGrid','on');
        legend(handles.MP_axes,'Mg','Pe');

    else
        updateSS(handles, rf_step);
    end

    plot(handles.rf_axes,handles.rfTime,handles.rfAmp,...
        [handles.Muts(rf_step) handles.Muts(rf_step)],[min(handles.rfAmp) max(handles.rfAmp)],'g-','linewidth',2);
    set(handles.rf_axes,'YLim',[min(handles.rfAmp) max(handles.rfAmp)],'XLim',[min(handles.rfTime) max(handles.rfTime)]);
    set(handles.rf_axes,'YGrid','on','XGrid','on');
    %set(handles.SpinRot_axes,'ZTickLabel',get(handles.SpinRot_axes,'ZTick')*10);
    %set(handles.MxMyMz_axes,'XTickLabel',get(handles.MxMyMz_axes,'XTick')*100);
    set(handles.rf_axes,'XTickLabel',get(handles.rf_axes,'XTick')*1000);
    set(handles.rf_axes,'YTickLabel',get(handles.rf_axes,'YTick')*10000);
    set(handles.rfFreq_axes,'YTickLabel',get(handles.rfFreq_axes,'YTick')/1000);

    pause(0.1);
end

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

[az,el] = view(handles.SpinRot_axes);
Starti=round(get(handles.rf_slider,'Value'));
for i=Starti+1:max(size(handles.Muts))
    handles=guidata(hObject);
    if handles.PauseFlag==1
        return;
    end
    set(handles.rf_slider,'Value',i);
    rf_step=round(get(handles.rf_slider,'Value'));

    if strcmp(handles.Freq_Flag,'off') & strcmp(handles.Spat_Flag,'off')

        Mx=handles.Mx(:,:,:,:,:,rf_step);
        My=handles.My(:,:,:,:,:,rf_step);
        Mz=handles.Mz(:,:,:,:,:,rf_step);
        Mxy=sqrt(Mx.^2+My.^2);
        Pe=angle(Mx+1i.*My);
        for j=1:str2double(get(handles.Attrh0.TypeNum,'String'))
            spinver=quiver3(handles.SpinRot_axes,handles.Gxgrid,handles.Gygrid,handles.Gzgrid,Mx(:,:,:,:,j),My(:,:,:,:,j),Mz(:,:,:,:,j));
            set(spinver,'AutoScale','off');
            hold(handles.SpinRot_axes,'on');
            plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mx(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
            hold(handles.MxMyMz_axes,'on');
            plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(My(:,:,:,:,j)),'Color',[0 1-(j-1)*0.25 0],'linewidth',2);
            plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

            plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
            hold(handles.MxyMz_axes,'on');
            plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

            plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
            hold(handles.MP_axes,'on');
            plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Pe(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);
        end
        hold(handles.SpinRot_axes,'off');
        hold(handles.MxMyMz_axes,'off');
        hold(handles.MxyMz_axes,'off');
        hold(handles.MP_axes,'off');

        set(handles.SpinRot_axes,'YLim',[-1 1],'XLim',[-1 1],'ZLim',[min(handles.Gzgrid,[],3)-1 max(handles.Gzgrid,[],3)+1]);
        xlabel(handles.SpinRot_axes,'X axis');
        ylabel(handles.SpinRot_axes,'Y axis');
        zlabel(handles.SpinRot_axes,'Z axis (Gradient Direction)');
        set(handles.SpinRot_axes,'view',[az,el]);

        set(handles.MxMyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
        set(handles.MxMyMz_axes,'YGrid','on','XGrid','on');
        legend(handles.MxMyMz_axes,'Mx','My','Mz');

        set(handles.MxyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
        set(handles.MxyMz_axes,'YGrid','on','XGrid','on');
        legend(handles.MxyMz_axes,'|Mxy|','Mz');

        set(handles.MP_axes,'YLim',[-pi pi],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
        set(handles.MP_axes,'YGrid','on','XGrid','on');
        legend(handles.MP_axes,'Mg','Pe');

    else
        updateSS(handles, rf_step);
    end

    plot(handles.rf_axes,handles.rfTime,handles.rfAmp,...
        [handles.Muts(rf_step) handles.Muts(rf_step)],[min(handles.rfAmp) max(handles.rfAmp)],'g-','linewidth',2);
    set(handles.rf_axes,'YLim',[min(handles.rfAmp) max(handles.rfAmp)],'XLim',[min(handles.rfTime) max(handles.rfTime)]);
    set(handles.rf_axes,'YGrid','on','XGrid','on');
    %set(handles.SpinRot_axes,'ZTickLabel',get(handles.SpinRot_axes,'ZTick')*10);
    %set(handles.MxMyMz_axes,'XTickLabel',get(handles.MxMyMz_axes,'XTick')*100);
    set(handles.rf_axes,'XTickLabel',get(handles.rf_axes,'XTick')*1000);
    set(handles.rf_axes,'YTickLabel',get(handles.rf_axes,'YTick')*10000);
    set(handles.rfFreq_axes,'YTickLabel',get(handles.rfFreq_axes,'YTick')/1000);

    pause(0.001);
end

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

set(handles.rf_slider,'Value',1);
rf_step=round(get(handles.rf_slider,'Value'));
[az,el] = view(handles.SpinRot_axes);

if strcmp(handles.Freq_Flag,'off') & strcmp(handles.Spat_Flag,'off')

    Mx=handles.Mx(:,:,:,:,:,rf_step);
    My=handles.My(:,:,:,:,:,rf_step);
    Mz=handles.Mz(:,:,:,:,:,rf_step);
    Mxy=sqrt(Mx.^2+My.^2);
    Pe=angle(Mx+1i.*My);
    for j=1:str2double(get(handles.Attrh0.TypeNum,'String'))
        spinver=quiver3(handles.SpinRot_axes,handles.Gxgrid,handles.Gygrid,handles.Gzgrid,Mx(:,:,:,:,j),My(:,:,:,:,j),Mz(:,:,:,:,j));
        set(spinver,'AutoScale','off');
        hold(handles.SpinRot_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mx(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MxMyMz_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(My(:,:,:,:,j)),'Color',[0 1-(j-1)*0.25 0],'linewidth',2);
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MxyMz_axes,'on');
        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MP_axes,'on');
        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Pe(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);
    end
    hold(handles.SpinRot_axes,'off');
    hold(handles.MxMyMz_axes,'off');
    hold(handles.MxyMz_axes,'off');
    hold(handles.MP_axes,'off');

    set(handles.SpinRot_axes,'YLim',[-1 1],'XLim',[-1 1],'ZLim',[min(handles.Gzgrid,[],3)-1 max(handles.Gzgrid,[],3)+1]);
    xlabel(handles.SpinRot_axes,'X axis');
    ylabel(handles.SpinRot_axes,'Y axis');
    zlabel(handles.SpinRot_axes,'Z axis (Gradient Direction)');
    set(handles.SpinRot_axes,'view',[az,el]);

    set(handles.MxMyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxMyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxMyMz_axes,'Mx','My','Mz');

    set(handles.MxyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxyMz_axes,'|Mxy|','Mz');

    set(handles.MP_axes,'YLim',[-pi pi],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MP_axes,'YGrid','on','XGrid','on');
    legend(handles.MP_axes,'Mg','Pe');

else
    updateSS(handles, rf_step);
end

plot(handles.rf_axes,handles.rfTime,handles.rfAmp,...
    [handles.Muts(rf_step) handles.Muts(rf_step)],[min(handles.rfAmp) max(handles.rfAmp)],'g-','linewidth',2);
set(handles.rf_axes,'YLim',[min(handles.rfAmp) max(handles.rfAmp)],'XLim',[min(handles.rfTime) max(handles.rfTime)]);
set(handles.rf_axes,'YGrid','on','XGrid','on');
%set(handles.SpinRot_axes,'ZTickLabel',get(handles.SpinRot_axes,'ZTick')*10);
%set(handles.MxMyMz_axes,'XTickLabel',get(handles.MxMyMz_axes,'XTick')*100);
set(handles.rf_axes,'XTickLabel',get(handles.rf_axes,'XTick')*1000);
set(handles.rf_axes,'YTickLabel',get(handles.rf_axes,'YTick')*10000);
set(handles.rfFreq_axes,'YTickLabel',get(handles.rfFreq_axes,'YTick')/1000);

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

set(handles.rf_slider,'Value',max(size(handles.Muts)));
rf_step=round(get(handles.rf_slider,'Value'));
[az,el] = view(handles.SpinRot_axes);

if strcmp(handles.Freq_Flag,'off') & strcmp(handles.Spat_Flag,'off')

    Mx=handles.Mx(:,:,:,:,:,rf_step);
    My=handles.My(:,:,:,:,:,rf_step);
    Mz=handles.Mz(:,:,:,:,:,rf_step);
    Mxy=sqrt(Mx.^2+My.^2);
    Pe=angle(Mx+1i.*My);
    for j=1:str2double(get(handles.Attrh0.TypeNum,'String'))
        spinver=quiver3(handles.SpinRot_axes,handles.Gxgrid,handles.Gygrid,handles.Gzgrid,Mx(:,:,:,:,j),My(:,:,:,:,j),Mz(:,:,:,:,j));
        set(spinver,'AutoScale','off');
        hold(handles.SpinRot_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mx(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MxMyMz_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(My(:,:,:,:,j)),'Color',[0 1-(j-1)*0.25 0],'linewidth',2);
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MxyMz_axes,'on');
        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);

        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,j)),'Color',[0 0 1-(j-1)*0.25],'linewidth',2);
        hold(handles.MP_axes,'on');
        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Pe(:,:,:,:,j)),'Color',[1-(j-1)*0.25 0 0],'linewidth',2);
    end
    hold(handles.SpinRot_axes,'off');
    hold(handles.MxMyMz_axes,'off');
    hold(handles.MxyMz_axes,'off');
    hold(handles.MP_axes,'off');

    set(handles.SpinRot_axes,'YLim',[-1 1],'XLim',[-1 1],'ZLim',[min(handles.Gzgrid,[],3)-1 max(handles.Gzgrid,[],3)+1]);
    xlabel(handles.SpinRot_axes,'X axis');
    ylabel(handles.SpinRot_axes,'Y axis');
    zlabel(handles.SpinRot_axes,'Z axis (Gradient Direction)');
    set(handles.SpinRot_axes,'view',[az,el]);

    set(handles.MxMyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxMyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxMyMz_axes,'Mx','My','Mz');

    set(handles.MxyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxyMz_axes,'|Mxy|','Mz');

    set(handles.MP_axes,'YLim',[-pi pi],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MP_axes,'YGrid','on','XGrid','on');
    legend(handles.MP_axes,'Mg','Pe');

else
    updateSS(handles, rf_step);
end

plot(handles.rf_axes,handles.rfTime,handles.rfAmp,...
    [handles.Muts(rf_step) handles.Muts(rf_step)],[min(handles.rfAmp) max(handles.rfAmp)],'g-','linewidth',2);
set(handles.rf_axes,'YLim',[min(handles.rfAmp) max(handles.rfAmp)],'XLim',[min(handles.rfTime) max(handles.rfTime)]);
set(handles.rf_axes,'YGrid','on','XGrid','on');
set(handles.rf_axes,'XTickLabel',get(handles.rf_axes,'XTick')*1000);
set(handles.rf_axes,'YTickLabel',get(handles.rf_axes,'YTick')*10000);
set(handles.rfFreq_axes,'YTickLabel',get(handles.rfFreq_axes,'YTick')/1000);

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

set(hObject,'Enable','off');
set(hObject,'String','Calc...');
pause(0.01);
ExecFlag=DorfExec(handles);
set(handles.Simuh.TimeWait_text,'String', ['Est. Time Left :  ' '~' ' : ' '~' ' : ' '~']);
handles=guidata(handles.rfDesignPanel_figure);
if ExecFlag==0
    set(hObject,'Enable','on');
    set(hObject,'String','Execute');
    return;
end
handles.PauseFlag=0;
set(handles.Leftend_pushbutton,'Enable','on');
set(handles.Rightend_pushbutton,'Enable','on');
set(handles.Pause_pushbutton,'Enable','on');
set(handles.Right_pushbutton,'Enable','on');
set(handles.DoubleRight_pushbutton,'Enable','on');
set(handles.rf_slider,'Enable','on');
%% Display
% rf Response & Spin Rotation
if strcmp(handles.Freq_Flag,'off') & strcmp(handles.Spat_Flag,'off')
    if ~isfield(handles,'MxMyMz_tab')
        %----clear tabs
        tabs=get(handles.SliProfile_tabgroup,'Children');
        for i=1:length(tabs)
            delete(get(tabs(i),'Children'));
        end
        delete(tabs);
        delete(handles.SliProfile_tabgroup);
        handles=rmfield(handles,'Mx_tab');
        %----end
        %----recreate tabs
        handles.SliProfile_tabgroup=uitabgroup(handles.SliProfile_uipanel);
        handles.MxMyMz_tab=uitab(handles.SliProfile_tabgroup,'title','Mx My Mz','Units','normalized');
        handles.MxyMz_tab=uitab(handles.SliProfile_tabgroup,'title','|Mxy| Mz','Units','normalized');
        handles.MP_tab=uitab(handles.SliProfile_tabgroup,'title','Mg Pe','Units','normalized');
        handles.MxMyMz_axes=axes('parent', handles.MxMyMz_tab,'Position', [0.06 0.1 0.88 0.85]);
        handles.MxyMz_axes=axes('parent', handles.MxyMz_tab,'Position', [0.06 0.1 0.88 0.85]);
        handles.MP_axes=axes('parent', handles.MP_tab,'Position', [0.06 0.1 0.88 0.85]);
        %----end
    end
    Mx=handles.Mx(:,:,:,:,:,end);
    My=handles.My(:,:,:,:,:,end);
    Mz=handles.Mz(:,:,:,:,:,end);
    Mxy=sqrt(Mx.^2+My.^2);
    Pe=angle(Mx+1i.*My);
    for i=1:str2double(get(handles.Attrh0.TypeNum,'String'))
        spinver=quiver3(handles.SpinRot_axes,handles.Gxgrid,handles.Gygrid,handles.Gzgrid,Mx(:,:,:,:,i),My(:,:,:,:,i),Mz(:,:,:,:,i));
        set(spinver,'AutoScale','off');
        hold(handles.SpinRot_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mx(:,:,:,:,i)),'Color',[0 0 1-(i-1)*0.25],'linewidth',2);
        hold(handles.MxMyMz_axes,'on');
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(My(:,:,:,:,i)),'Color',[0 1-(i-1)*0.25 0],'linewidth',2);
        plot(handles.MxMyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,i)),'Color',[1-(i-1)*0.25 0 0],'linewidth',2);

        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,i)),'Color',[0 0 1-(i-1)*0.25],'linewidth',2);
        hold(handles.MxyMz_axes,'on');
        plot(handles.MxyMz_axes,squeeze(handles.Gzgrid),squeeze(Mz(:,:,:,:,i)),'Color',[1-(i-1)*0.25 0 0],'linewidth',2);

        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Mxy(:,:,:,:,i)),'Color',[0 0 1-(i-1)*0.25],'linewidth',2);
        hold(handles.MP_axes,'on');
        plot(handles.MP_axes,squeeze(handles.Gzgrid),squeeze(Pe(:,:,:,:,i)),'Color',[1-(i-1)*0.25 0 0],'linewidth',2);
    end
    hold(handles.SpinRot_axes,'off');
    hold(handles.MxMyMz_axes,'off');
    hold(handles.MxyMz_axes,'off');
    hold(handles.MP_axes,'off');
    set(handles.SpinRot_axes,'YLim',[-1 1],'XLim',[-1 1],'ZLim',[min(handles.Gzgrid,[],3)-1 max(handles.Gzgrid,[],3)+1]);
    xlabel(handles.SpinRot_axes,'X axis');
    ylabel(handles.SpinRot_axes,'Y axis');
    zlabel(handles.SpinRot_axes,'Z axis (Gradient Direction)');
    set(handles.SpinRot_axes,'view',[134,24]);

    set(handles.MxMyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    legend(handles.MxMyMz_axes,'Mx','My','Mz');
    set(handles.MxMyMz_axes,'YGrid','on','XGrid','on');

    set(handles.MxyMz_axes,'YLim',[-1 1],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MxyMz_axes,'YGrid','on','XGrid','on');
    legend(handles.MxyMz_axes,'|Mxy|','Mz');

    set(handles.MP_axes,'YLim',[-pi pi],'XLim',[min(handles.Gzgrid,[],3) max(handles.Gzgrid,[],3)]);
    set(handles.MP_axes,'YGrid','on','XGrid','on');
    legend(handles.MP_axes,'Mg','Pe');
else

    if ~isfield(handles,'Mx_tab')
        %----clear tabs
        tabs=get(handles.SliProfile_tabgroup,'Children');
        for i=1:length(tabs)
            delete(get(tabs(i),'Children'));
        end
        delete(tabs);
        delete(handles.SliProfile_tabgroup);
        handles=rmfield(handles,'MxMyMz_tab');
        %----end
        %----recreate tab
        handles.SliProfile_tabgroup=uitabgroup(handles.SliProfile_uipanel);
        handles.Mx_tab=uitab(handles.SliProfile_tabgroup,'title','Mx','Units','normalized');
        handles.My_tab=uitab(handles.SliProfile_tabgroup,'title','My','Units','normalized');
        handles.Mz_tab=uitab(handles.SliProfile_tabgroup,'title','Mz','Units','normalized');
        handles.M_tab=uitab(handles.SliProfile_tabgroup,'title','Mag','Units','normalized');
        handles.P_tab=uitab(handles.SliProfile_tabgroup,'title','Ph','Units','normalized');

        handles.Mx_axes=axes('parent', handles.Mx_tab,'Position', [0.06 0.1 0.88 0.85]);
        handles.My_axes=axes('parent', handles.My_tab,'Position', [0.06 0.1 0.88 0.85]);
        handles.Mz_axes=axes('parent', handles.Mz_tab,'Position', [0.06 0.1 0.88 0.85]);
        handles.M_axes=axes('parent', handles.M_tab,'Position', [0.06 0.1 0.88 0.85]);
        handles.P_axes=axes('parent', handles.P_tab,'Position', [0.06 0.1 0.88 0.85]);
        %----end
    end

    updateSS(handles, size(handles.Mx, 3));
    set(handles.SpinRot_axes,'view',[134,24]);
end

set(handles.rf_slider,'Min',1);
set(handles.rf_slider,'Max',max(size(handles.Muts)));
set(handles.rf_slider,'Value',max(size(handles.Muts)));

plot(handles.rf_axes,handles.rfTime,handles.rfAmp,[handles.Muts(get(handles.rf_slider,'Value')) handles.Muts(get(handles.rf_slider,'Value'))],[min(handles.rfAmp) max(handles.rfAmp)],'g-','linewidth',2);
plot(handles.rfPhase_axes,handles.rfTime,handles.rfPhase,'linewidth',2);
plot(handles.rfFreq_axes,handles.rfTime,handles.rfFreq,'linewidth',2);

if strcmp(handles.Spat_Flag,'off')
    plot(handles.Gz_axes,handles.GzTime,handles.GzAmp,'linewidth',2);
    set(handles.Gz_axes,'YGrid','on','XGrid','off','XTick',[]);
    set(handles.Gz_axes,'YTickLabel',get(handles.Gz_axes,'YTick')*100);
    linkaxes([handles.rf_axes handles.rfPhase_axes handles.rfFreq_axes handles.Gz_axes],'x');
else
    plot(handles.Gy_axes,handles.GyTime,handles.GyAmp,'linewidth',2);
    set(handles.Gy_axes,'YGrid','on','XGrid','off','XTick',[]);
    set(handles.Gy_axes,'YTickLabel',get(handles.Gy_axes,'YTick')*100);
    plot(handles.Gx_axes,handles.GxTime,handles.GxAmp,'linewidth',2);
    set(handles.Gx_axes,'YGrid','on','XGrid','off','XTick',[]);
    set(handles.Gx_axes,'YTickLabel',get(handles.Gx_axes,'YTick')*100);
    linkaxes([handles.rf_axes handles.rfPhase_axes handles.rfFreq_axes handles.Gy_axes handles.Gx_axes ],'x');
end

set(handles.rf_axes,'YLim',[min(handles.rfAmp) max(handles.rfAmp)],'XLim',[min(handles.rfTime) max(handles.rfTime)]);
set(handles.rf_axes,'YGrid','on','XGrid','on');
set(handles.rfPhase_axes,'YGrid','on','XGrid','off','XTick',[]);
set(handles.rfFreq_axes,'YGrid','on','XGrid','off','XTick',[]);

set(handles.rf_axes,'XTickLabel',get(handles.rf_axes,'XTick')*1000);
set(handles.rf_axes,'YTickLabel',get(handles.rf_axes,'YTick')*10000);
set(handles.rfFreq_axes,'YTickLabel',get(handles.rfFreq_axes,'YTick')/1000);

set(hObject,'Enable','on');
set(hObject,'String','Execute');
guidata(handles.rfDesignPanel_figure, handles);

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

% update SeqElem.xml
if verLessThan('matlab','8.5')
    % Code to run in MATLAB R2015a and earlier here
    idx = get(handles.SeqAttri_tabgroup,'SelectedIndex');
else
    % Code to run in MATLAB R2015a and later here
    tabtitle=get(get(handles.SeqAttri_tabgroup,'SelectedTab'),'Title');
    if tabtitle(1:2) == 'rf'
        idx = 1;
    elseif tabtitle(1:2) == 'Gz'
        idx = 2;
    end
end

if idx==1 % update rf
    if ~isempty(handles.Attrh1)
        Attrh1name=fieldnames(handles.Attrh1);
        for i=1:2:length(Attrh1name)

            if ~iscell(get(handles.Attrh1.(Attrh1name{i+1}),'String'))
                handles.rfNode.Attributes((i+1)/2).Value=get(handles.Attrh1.(Attrh1name{i+1}),'String');
            else
                handles.rfNode.Attributes((i+1)/2).Value(2)=num2str(get(handles.Attrh1.(Attrh1name{i+1}),'Value'));
            end
        end
        eval([handles.rfNodeLvl '=handles.rfNode;']);
        DoWriteXML(handles.SeqElemListStruct,[handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'SeqElem.xml']);
    end

    fid=fopen([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'rf' filesep handles.SeqElemNode.Name '_Memo.txt'],'wt+');
    Memo=get(handles.rfMemo_edit,'String');
    if ~isempty(Memo)
        if ischar(Memo)
            Memo=cellstr(Memo);
        end
        for i=1:length(Memo)
            fprintf(fid,'%s\n', Memo{i,1});
        end
    end
    fclose(fid);

elseif idx==2 % update GzSS
    if ~isempty(handles.Attrh2)
        Attrh2name=fieldnames(handles.Attrh2);
        for i=1:2:length(Attrh2name)

            if ~iscell(get(handles.Attrh2.(Attrh2name{i+1}),'String'))
                handles.GzNode.Attributes((i+1)/2).Value=get(handles.Attrh2.(Attrh2name{i+1}),'String');
            else
                handles.GzNode.Attributes((i+1)/2).Value(2)=num2str(get(handles.Attrh2.(Attrh2name{i+1}),'Value'));
            end
        end
        eval([handles.GzNodeLvl '=handles.GzNode;']);
        DoWriteXML(handles.SeqElemListStruct,[handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'SeqElem.xml']);
    end
end

% update rfSimuAttr.xml
tabs=get(handles.Setting_tabgroup,'Children');
for j=1:length(handles.rfSimuAttrStruct.Children)
    Attrh0=get(tabs(j),'Children');
    for i=1:2:length(Attrh0)
        if ~iscell(get(Attrh0(end-i),'String'))
            handles.rfSimuAttrStruct.Children(j).Attributes((i+1)/2).Value=get(Attrh0(end-i),'String');
        else
            handles.rfSimuAttrStruct.Children(j).Attributes((i+1)/2).Value(2)=num2str(get(Attrh0(end-i),'Value'));
        end
    end
end
DoWriteXML(handles.rfSimuAttrStruct,[handles.Simuh.MRiLabPath filesep 'Macro' filesep 'SeqElem' filesep 'rfSimuAttr.xml']);

guidata(hObject, handles);

% --------------------------------------------------------------------
function SavePanel_uipushtool_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to SavePanel_uipushtool (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
DoSaveSnapshot(handles.rfDesignPanel_figure);

% --- Executes on button press in Pause_pushbutton.
function Pause_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Pause_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if handles.PauseFlag==1
    handles.PauseFlag=0;
    set(hObject,'String','X');
else
    handles.PauseFlag=1;
    set(hObject,'String','O');
end

guidata(hObject, handles);

% --- Executes when user attempts to close rfDesignPanel_figure.
function rfDesignPanel_figure_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to rfDesignPanel_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure

if isfield(handles,'PauseFlag')
    if handles.PauseFlag==1
        delete(hObject);
    else
        warndlg('Please pause rf animation before closing the window (i.e. press X).');
    end
else
    delete(hObject);
end

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

% Hints: get(hObject,'String') returns contents of rfMemo_edit as Gz_text
%        str2double(get(hObject,'String')) returns contents of rfMemo_edit as a double

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

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

% update spin response for spatial-spectral analysis
function updateSS(handles, rf_step)

Mx=handles.Mx(:,:,rf_step);
My=handles.My(:,:,rf_step);
Mz=handles.Mz(:,:,rf_step);
Mag=sqrt(Mx.^2+My.^2);
Ph=angle(Mx+1i.*My);

axes(handles.Mx_axes);
imagesc(Mx);
colormap gray
colorbar
axes(handles.My_axes);
imagesc(My);
colorbar
axes(handles.Mz_axes);
imagesc(Mz);
colorbar
axes(handles.M_axes);
imagesc(Mag);
colorbar
axes(handles.P_axes);
imagesc(Ph);
colorbar

axes(handles.SpinRot_axes);
[az,el] = view(handles.SpinRot_axes);
cla(handles.SpinRot_axes);

if verLessThan('matlab','8.5')
    % Code to run in MATLAB R2015a and earlier here
    idx = get(handles.SliProfile_tabgroup,'SelectedIndex');
    switch idx
        case 1 % Mx
            surf(double(Mx));
            title('Mx');
        case 2 % My
            surf(double(My));
            title('My');
        case 3 % Mz
            surf(double(Mz));
            title('Mz');
        case 4 % Mag
            surf(double(Mag));
            title('Mag');
        case 5 % Ph
            surf(double(Ph));
            title('Ph');
    end
else
    % Code to run in MATLAB R2015a and later here
    tabtitle=get(get(handles.SliProfile_tabgroup,'SelectedTab'),'Title');
    switch tabtitle
        case 'Mx' % Mx
            surf(double(Mx));
            title('Mx');
        case 'My' % My
            surf(double(My));
            title('My');
        case 'Mz' % Mz
            surf(double(Mz));
            title('Mz');
        case 'Mag' % Mag
            surf(double(Mag));
            title('Mag');
        case 'Ph' % Ph
            surf(double(Ph));
            title('Ph');
    end
end

if strcmp(handles.Freq_Flag,'on')
    xlabel('Frequency');
    ylabel('Distance');
else
    xlabel('X-Distance');
    ylabel('Y-Distance');
end
shading interp;
set(handles.SpinRot_axes,'view',[az,el],'xdir','reverse');

% --- Executes on mouse press over figure background, over a disabled or
% --- inactive control, or over an axes background.
function rfDesignPanel_figure_WindowButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to rfDesignPanel_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

tempaxes = gca;
try
    if strcmp(get(gcf,'Selectiontype'),'alt')
        figure 			% Create a new figure
        ax_new = axes;
        copyaxes(tempaxes, ax_new);
    end
catch me
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [model] = integrate_cVector_into_model(model)
% sintegrate_cVector_into_model incorporates the c vector directly into the
% model as the reaction 'obj_fun_rxn'.

% USAGE:
%
%    [model] = integrate_cVector_into_model(model)
%
% INPUTS:
%   model:  COBRA model structure with minimal fields:

% OUTPUT:
%   model:   Augmented COBRA model
%
%
% .. Authors: - Bronson R. Weston   2022

model=addMetabolite(model, 'obj_fun_met');
model.S(end,:)=model.c;
model = addMultipleReactions(model, {'obj_fun_rxn'}, {'obj_fun_met'}, [-1], 'lb', [-1e7], 'ub', [1e7]);

model.c=zeros(length(model.c),1);
model.c(end)=1;

end

author: opencobra

repo: cobratoolbox


text:
function varargout = gridmeshToQuadmesh(x, y, varargin)
%GRIDMESHTOQUADMESH Create a quad mesh from a grid mesh
%
%   Deprecated: replaced by surfToMesh (2012.10.25)
%
%   [V F] = gridmeshToQuadmesh(X, Y)
%   [V F] = gridmeshToQuadmesh(X, Y, Z)
%   Converts the surface grid given by two or three coordinate arrays into
%   a face-vertex quad mesh.
%
%   Example
%     % transform a surface into a mesh
%     [X,Y] = meshgrid(-2:.2:2, -2:.2:2);
%     Z = X .* exp(-X.^2 - Y.^2);
%     [V F] = gridmeshToQuadmesh(X, Y, Z);
%     figure;
%     drawMesh(V, F);
%
%     % Transform surface of a cylinder as a mesh
%     [x y z] = cylinder(5*ones(1, 10));
%     [v f] = gridmeshToQuadmesh(x, y, z);
%     figure;
%     drawMesh(v, f);
%     view(3); axis equal;
%
%   See also
%     meshes3d, meshgrid, drawMesh
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-12-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

warning('MatGeom:deprecated', ...
    'function "gridmeshToQuadmesh" is deprecated, and was replaced by "surfToMesh".');

% number of vertices
dim = size(x);
nv = prod(dim);

% vertex indices in grid
inds = reshape(1:nv, dim);

% create vertices
if isempty(varargin)
    vertices = [x(:) y(:)];
else
    z = varargin{1};
    vertices = [x(:) y(:) z(:)];
end

% create faces
v1 = inds(1:end-1, 1:end-1);
v2 = inds(1:end-1, 2:end);
v3 = inds(2:end, 2:end);
v4 = inds(2:end, 1:end-1);

faces = [v1(:) v2(:) v3(:) v4(:)];

% format output
if nargout <= 1
    % concatenate into one structure
    mesh.vertices = vertices;
    mesh.edges = edges;
    mesh.faces = faces;
    varargout = {mesh};

elseif nargout == 2
    % returns as separate arguments
    varargout = {vertices, faces};

elseif nargout == 3
    % also return vertex indices
    varargout = {vertices, faces, inds};

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% changeFileExtensions
% Below is a demonstration of the features of the |changeFileExtensions| function

%%
clear; close all; clc;

%% Syntax
% |changeFileExtensions(pathName,extOld,extNew);|

%% Description
% The |changeFileExtensions| function changes the file extensions for all
% files in the folder pathName, and with the extension extOld, to have the
% extension extNew.

%% Examples
%

%%
%Create example files in the data/temp directory. Here a set of files with
%the txt extension are created. Later these are changed to have a csv
%file extension.

%Create .txt files
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
pathName=fullfile(defaultFolder,'data','temp','renameCheck');
mkdir(pathName)

n=5;
for q=1:1:n %Create n files
    fileID=fopen(fullfile(pathName,['temp',num2str(q),'.txt']),'w');
    fprintf(fileID,'%d\n',pi);
    fclose(fileID);
end

% Add one additional file with a different extension to show this file will
% not be effected.
fileID=fopen(fullfile(pathName,'temp.txp'),'w');
fprintf(fileID,'%d\n',pi);
fclose(fileID);

%%
% Show current folder content
disp('Old folder content:')
ls(pathName)

%%
% Change file extensions
extOld='txt'; %Old extension
extNew='csv'; %New extension
changeFileExtensions(pathName,extOld,extNew)

%%
% Show current folder content
disp('New folder content:')
ls(pathName)

%%
% remove the temporary folder created for this example
rmdir(pathName,'s')

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function fem2d_pack_test22 ( )

%*****************************************************************************80
%
%% TEST22 tests SPHERE_GRID_T3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  element_order = 3;
  nelemx = 8;
  nelemy = 8;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST22\n' );
  fprintf ( 1, '  SPHERE_GRID_T3_ELEMENT sets up a grid of T3 triangles\n' );
  fprintf ( 1, '    on a sphere.\n' );
  fprintf ( 1, '  SPHERE_GRID_T3_ELEMENT_NUM returns the number\n' );
  fprintf ( 1, '    of elements in the grid\n' );
  fprintf ( 1, '  SPHERE_GRID_T3_NODE_NUM returns the number\n' );
  fprintf ( 1, '    of nodes in the grid.\n' );
  fprintf ( 1, '  SPHERE_GRID_T3_NODE_XYZ returns the coordinates\n' );
  fprintf ( 1, '    of nodes in the grid.\n' );

  element_num = sphere_grid_t3_element_num ( nelemx, nelemy );
  node_num = sphere_grid_t3_node_num ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected number of nodes =    %d\n', node_num );
  fprintf ( 1, '  Expected number of elements = %d\n', element_num );

  element_node = sphere_grid_t3_element ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The elements and their nodes:\n' );
  fprintf ( 1, '\n' );

  for element = 1 : element_num
    fprintf ( 1, '%4d  ', element );
    for order = 1 : element_order
      fprintf ( 1, '%4d', element_node(order,element) );
    end
    fprintf ( 1, '\n' );
  end

  node_xyz = sphere_grid_t3_node_xyz ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The node coordinates:\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num
    fprintf ( 1, '  %4d  %12f  %12f  %12f\n', node, node_xyz(1:3,node) );
  end
%
%  Write the elements and nodes to files.
%
  r8mat_write ( 'sphere_t3_nodes.txt', 3, node_num, node_xyz );

  i4mat_write ( 'sphere_t3_elements.txt', element_order, element_num, ...
    element_node );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function orderedNoiseDisplay(noise, spacing)

% ORDEREDNOISEDISPLAY Display parameters of the ORDERED noise.
% FORMAT
% DESC displays the parameters of the ordered categorical
% noise and the noise type to the console.
% ARG noise : the noise to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG noise : the noise to display.
% ARG spacing : how many spaces to indent the display of the noise by.
%
% SEEALSO : orderedNoiseParamInit, modelDisplay, noiseDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
for i = 1:noise.numProcess
  fprintf(spacing);
  fprintf('Ordered bias on process %d: %2.4f\n', i, noise.bias(i))
end
for i = 1:noise.C-2
  fprintf('Ordered noise model width %d: %2.4f\n', i, noise.widths(i))
end
fprintf(spacing);
fprintf(spacing);
fprintf('Ordered Sigma2: %2.4f\n', noise.variance);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [document, scores] = report_difficulty(context, experiment, trackers, sequences, varargin)
% report_difficulty Generate a difficulty report for tags or sequences
%
% Performs A-R ranking analysis and generates a report that shows the difficulty of individual
% tags or sequences.
%
% Input:
% - context (structure): Report context structure.
% - experiment (struct): An experiment structure.
% - trackers (cell): An array of tracker structures.
% - sequences (cell): An array of sequence structures.
% - varargin[UsePractical] (boolean): Use practical difference.
% - varargin[UseTags] (boolean): Rank according to tags (otherwise rank according to sequences).
% - varargin[Alpha] (boolean): Statistical significance parameter.
%
% Output:
% - document (structure): Resulting document structure.
%

usetags = true;
usepractical = true;
alpha = 0.05;
adaptation = get_global_variable('report_ranking_adaptation', 'mean');

scores = [];

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'usepractical'
            usepractical = varargin{i+1};
        case 'usetags'
            usetags = varargin{i+1};
        case 'alpha'
            alpha = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i}, '!']) ;
    end
end

document = create_document(context, 'difficulty', 'title', 'Difficulty');

trackers_hash = md5hash(strjoin((cellfun(@(x) x.identifier, trackers, 'UniformOutput', false)), '-'), 'Char', 'hex');
parameters_hash = md5hash(sprintf('%f-%d-%d-%s', alpha, usetags, usepractical, adaptation));

if ~strcmp(experiment.type, 'supervised')
   error('Difficulty analysis only suitable for supervised experiments!');
end

document.section('Experiment %s', experiment.name);

tags = {};

if usetags && isfield(experiment, 'tags')
    tags = experiment.tags;

    sequences_hash = md5hash(strjoin(tags, '-'), 'Char', 'hex');
else
    sequences_hash = md5hash(strjoin((cellfun(@(x) x.name, sequences, 'UniformOutput', false)), '-'), 'Char', 'hex');
end;

cache_identifier = sprintf('ranking_%s_%s_%s_%s', experiment.name, trackers_hash, sequences_hash, parameters_hash);

result = report_cache(context, cache_identifier, @analyze_ranks, experiment, trackers, ...
    sequences, 'tags', tags, 'usepractical', usepractical, ...
    'alpha', alpha, 'adaptation', adaptation);

selector_tags = result.tags;

median_accuracy = nanmedian(result.accuracy.values, 2);
median_robustness = nanmedian(result.robustness.normalized, 2);

median_diff_accuracy = nanmedian(bsxfun(@minus, result.accuracy.values, mean(result.accuracy.values, 1)), 2);
median_diff_robustness = nanmedian(bsxfun(@minus, result.robustness.normalized, mean(result.robustness.normalized, 1)), 2);

table_data = [median_accuracy, median_robustness * 100, median_diff_accuracy, median_diff_robustness * 100];

table_data = highlight_best_rows(num2cell(table_data), {'ascending', 'descending', 'ascending', 'descending'})';
row_tags = {'Accuracy', 'Robustness', 'Accuracy difference', 'Robustness difference'};

title = sprintf('Difficulty for experiment %s', experiment.name);

document.table(table_data, 'columnLabels', selector_tags, 'rowLabels', row_tags', 'title', title);

document.raw('
\n');

hf = figure('Visible', 'off');
hold on;
for i = 1:numel(selector_tags)
    if median_diff_accuracy(i) > 0
        rectangle('Position',[i-0.3, 0.5, 0.6, median_diff_accuracy(i)], 'FaceColor', 'green');
    elseif median_diff_accuracy(i) < 0
        rectangle('Position',[i-0.3, 0.5+median_diff_accuracy(i), 0.6, -median_diff_accuracy(i)], 'FaceColor', 'red');
    end;

    plot(ones(size(result.accuracy.values(i, :))) * i, result.accuracy.values(i, :), 'bo');
end
hold off;
set(gca, 'YLim', [0, 1], 'XLim', [0.5, numel(selector_tags) + 0.5], 'XTick', 1:numel(selector_tags), 'XTickLabel', selector_tags);

document.figure(hf, sprintf('difficulty_%s_accuracy_scatter', experiment.name), ...
    sprintf('Accuracy scatter'));

hf = figure('Visible', 'off');
hold on;
origin = mean(median_robustness);
for i = 1:numel(selector_tags)
    if median_diff_robustness(i) > 0
        rectangle('Position',[i-0.3, origin, 0.6, median_diff_robustness(i)], 'FaceColor', 'red');
    elseif median_diff_robustness(i) < 0
        rectangle('Position',[i-0.3, origin+median_diff_robustness(i), 0.6, -median_diff_robustness(i)], 'FaceColor', 'green');
    end;

    plot(ones(size(result.robustness.normalized(i, :))) * i, result.robustness.normalized(i, :), 'bo');
end
hold off;
set(gca, 'YLim', [0, max(result.robustness.normalized(:))], ...
    'XLim', [0.5, numel(selector_tags) + 0.5], 'XTick', 1:numel(selector_tags), 'XTickLabel', selector_tags);

document.figure(hf, sprintf('difficulty_%s_robustness_scatter', experiment.name), ...
    sprintf('Robustness scatter'));

document.raw('
\n');

document.write();

end

author: votchallenge

repo: toolkit-legacy


text:
function [imageCha, ssim_map, metrics] = algo_BFCSA_proxA_2D_real( obj,input )
% 2D real-valued variant of FCSA with an outer Bregman Iteration
% based on Huang et al. paper on FCSA
%
% input:
% obj           CS reconstruction object (holding all parameters)
% input         struct containing recon parameters and image
%
% output:
% imageCha      reconstructed channel individual image
% ssim_map      structural similarity map
% metrics       evaluation metrics
%
% (c) Marc Fischer, Thomas Kuestner, May 2015
% -------------------------------------------------------------------------

%%
timer_proxA = tic;

%% variables:
% internal flags
flag_wavetree = true;
flag_fast = true;
flag_extendImage = true;

% internal variables:
L = 1.5; % > 1 suggested, otherwise it may diverge depending on "b{1,j} - FFT2D_mask(z{1,j},mask)"
t_old = 1;
NLTV_struct.kernelratio = 3;
NLTV_struct.windowratio = 6;
NLTV_struct.nThreads = 1;  % mind this option if used with -singleCompThread on BWCluster
itrNLTV = obj.maxitr - 5;
% chambolle tv:
parsin.MAXITER=100; parsin.tv='iso'; % 'iso' or 'l1'

% from obj:
% maxitr = obj.maxitr;
maxitr = obj.iNINNER;
% maxitrOUTER = obj.maxitrOUTER;
maxitrOUTER = obj.iNOUTER;
maxitrINNER = ceil( maxitr / maxitrOUTER );
% n1 = obj.measPara.dim(1);
% n2 = obj.measPara.dim(2);
n1 = input.n1;
n2 = input.n2;
% nSlices = obj.measPara.dim(3);
nCha = obj.measPara.dim(5);
lambdaWave = obj.lambda;
lambdaTV = obj.lambdaTV;
lambdaGroup = obj.lambdaGroup;
NLTV_struct.filterstrength = obj.lambdaNLTV_h; % 0.03 % converted from NLTV h = 0.01 %old: used: 3e-10
lambdaNLTV = obj.lambdaNLTV;
lambdaNLTV_h = obj.lambdaNLTV_h;
regularizerWeights = obj.regularizerWeights;
flagTV = obj.flagTV;
flagTV_iso = obj.flagTV_iso;
flagWave = obj.flagWave;
flagGroup = obj.flagGroup;
flagNLTV = obj.flagNLTV;
flagSBNLTV = obj.flagSBNLTV;
waveletStages = obj.trafo.waveletStages;
waveletFilterName_l1 = obj.trafo.waveletFilterName_l1;
waveletFilterName_l12 = obj.trafo.waveletFilterName_l12;

% from input:
b=input.b;
mask = input.mask;
G_prox = input.G_prox;
Gt_prox = input.Gt_prox;
groupnorm_index = input.groupnorm_index;
waveS_l1 = input.waveS_l1;
waveS_l12 = input.waveS_l12;
waveS_l12proxA = input.waveS_l12proxA;
proxA_extend_y = waveS_l12proxA(waveletStages+2,1) - waveS_l12(waveletStages+2,1);
proxA_extend_x = waveS_l12proxA(waveletStages+2,2) - waveS_l12(waveletStages+2,2);
z_proxA = cell(1,nCha);

% im_ref = input.im_ref;
% im_ref_full = zeros(n1,n2);
% for j = 1:nCha
%     im_ref_full = im_ref_full + abs(im_ref{1,j}).^2;
% end;
% im_ref_full = sqrt(im_ref_full);
clear input

% initialize cells/vectors
for j=1:nCha
    FTb{1,j} = real(iFFT2D(b{1,j}));
    % y{1,j} = real(FTb{1,j});
    % y{1,j+nCha} = imag(FTb{1,j});
end;
z = FTb; % starting point
b_outer = b;

x_wave = cell(1,nCha);
x_helper = x_wave;
x_wave_helper = x_wave;
x_tv = x_wave;
x_nltv = x_wave;
for j=1:nCha
    x_nltv{1,j} = 0;
end;
x_g = x_wave;
x_g_helper = x_wave;
x_g_proxA = x_wave;
y = x_wave;

%% MAD dependent lambdas:
flag_MAD = true;
if flag_MAD
    x_wavedec = cell(1,nCha);
    threshold = zeros(1:nCha);
    for j=1:nCha % 2*nCha
        x_wavedec{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
        x_wave_fine_scale = size(x_wavedec{1,j},2) - (3*waveS_l1(waveletStages+1,1)*waveS_l1(waveletStages+1,2));
        threshold(j) = mad(x_wavedec{1,j}(x_wave_fine_scale:end),1);
    end;
    clear x_wavedec
else
    threshold(1:nCha) = 1;
end;

threshold_wave(j) = lambdaWave * threshold(j) * 2/L;
threshold_TV(j) = lambdaTV * threshold(j) * 2/L;
threshold_group(j) = lambdaGroup * threshold(j) * 2/L;
threshold_NLTV(j) = lambdaNLTV; % * threshold(j) * 2/L;
threshold_NLTV_h(j) = lambdaNLTV_h; %*threshold(j); % adjust carefully or NLTV won't find a solution. lambdaNLTV_h should be < 0.01

%% initialize metrics:
itr = 0;
metrics.xtime(itr+1)= 0;
% [metrics, ssim_map{1,1}] = get_metrics_itr( im_ref, im_ref_full, z, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );
ssim_map = [];

%% recon

dispProgress('Proximal Average', 0, maxitrOUTER*maxitrINNER);
itr = 0;
for itrOUTER = 1:maxitrOUTER  % total iter counter
    for itrINNER = 1:maxitrINNER

        t_new = (1+sqrt(1+4*t_old^2))/2;
        t_old = t_new;

        y_old = z;   % y_old = y for complex case

        %% landweber step
        for j = 1:nCha
            x_helper{1,j} = real(iFFT2D(FFT2D_mask(z{1,j},mask))) -FTb{1,j} ;
            z{1,j} = z{1,j} - real(x_helper{1,j})/L;
        end;

        %% l1-Wavelet
        if flagWave
            for j = 1:nCha % 2*nCha
                x_wave_helper{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l1);
                x_wave_helper{1,j} = softthresh_real(x_wave_helper{1,j},threshold_wave(j));
                x_wave{1,j} = waverec2(x_wave_helper{1,j},waveS_l1,waveletFilterName_l1);
            end;
        end;

        %% TV
        if flagTV
            if ~flagTV_iso
                for j = 1:nCha
                    x_tv{1,j} = MTV_2D(z{1,j},threshold_TV(j),n1,n2);
                end;
            else
                for j = 1:nCha
                    if (itr < 2)
                        [x_tv{1,j}, P]=denoise_TV_One((z{1,j}), threshold_TV(j),-inf,inf,[],parsin);
                    else
                        [x_tv{1,j}, P]=denoise_TV_One((z{1,j}), threshold_TV(j),-inf,inf,P,parsin);
                    end;
                end;
            end;
        end;

        %% NLTV
        if flagNLTV
            if itr >= itrNLTV
                if flagSBNLTV
                    for j = 1:nCha
                        x_nltv{1,j} = SB_NLTVfunc_slim_rescale(z{1,j},n1,n2, threshold_NLTV(j), threshold_NLTV_h(j) );
                    end;
                else
                    for j = 1:nCha
                        % if mod(itr,5) == 0 || itr == 1
                        x_nltv{1,j} = NLMF(z{1,j},NLTV_struct);
                        x_nltv{1,j} = (L.*z{1,j} + 2*threshold_NLTV(j)*x_nltv{1,j})./(L+2*threshold_NLTV(j));
                        % end;
                    end;
                end;
            end;
        end;

        %% l12-Wavelet
        if flagGroup
            for j = 1:nCha
                if flag_extendImage
                    z_proxA{1,j} = extend_image(z{1,j}, waveS_l12proxA, waveletStages, proxA_extend_y, proxA_extend_x);
                    x_g_helper{1,j} = wavedec2(z_proxA{1,j},waveletStages,waveletFilterName_l12);
                else
                    x_g_helper{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l12);
                end;

                if flag_wavetree
                    x_g_helper{2,j} = (G_prox*x_g_helper{1,j}')';
                else
                    x_g_helper{2,j} = zeros(1,size(x_g_helper{1,j},2));
                end;
            end;

            x_g_helper(1,:) = softthresh_proxA_cha(x_g_helper,threshold_group(j),nCha,waveS_l12proxA,groupnorm_index);

            for j = 1:nCha
                x_g_proxA{1,j} = waverec2(x_g_helper{1,j},waveS_l12proxA,waveletFilterName_l12);
                x_g{1,j} =  x_g_proxA{1,j}(1:end-proxA_extend_y,1:end-proxA_extend_x);
            end;
        end;

        %% add prox(.)
        for j = 1:nCha
            y{1,j} = zeros(n1,n2);
            if flagWave y{1,j} = y{1,j} + x_wave{1,j}.*regularizerWeights(1); end;
            if flagTV y{1,j} = y{1,j} + x_tv{1,j}.*regularizerWeights(2); end;
            if flagGroup y{1,j} = y{1,j} + x_g{1,j}.*regularizerWeights(3); end;
            if flagNLTV y{1,j} = y{1,j} + x_nltv{1,j}.*regularizerWeights(4); end;

            if itr < itrNLTV
                y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagWave.*regularizerWeights(1) + flagGroup.*regularizerWeights(3));
            else
                y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagNLTV.*regularizerWeights(4) + flagWave.*regularizerWeights(1) + flagGroup.*regularizerWeights(3));
            end;

            if flag_fast
                y{1,j}=y{1,j}+((t_old-1)/t_new).*(y{1,j}-y_old{1,j});
            end;
        end;

        %% metrics of current itr:
        itr = itr + 1;
%         disp(itr);
        dispProgress('Proximal Average', itr/(maxitrOUTER*maxitrINNER));

        metrics.xtime(itr+1)= toc(timer_proxA);
        for j = 1:nCha
            z{1,j} = y{1,j};
        end;
%         [metrics, ssim_map{1,2}] = get_metrics_itr( im_ref, im_ref_full, z, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );

        if itr >= maxitr break; end;
    end;
    if itr >= maxitr break; end;
    for j=1:nCha
        b_outer{1,j} = b_outer{1,j} + b{1,j} - FFT2D_mask(z{1,j},mask);
        FTb{1,j} = real(iFFT2D(b_outer{1,j}));
    end;
end;
dispProgress('Proximal Average', 'Close');

imageCha = z;
for j = 1:nCha
    imageCha{1,j} = turn_image( imageCha{1,j} );
end;
% for j = 1:nCha+1
%     ssim_map{1,1}{1,j} = turn_image( ssim_map{1,1}{1,j} );
%     ssim_map{1,2}{1,j} = turn_image( ssim_map{1,2}{1,j} );
% end;

end

author: thomaskuestner


text:
function [P, X] = Get_Parameters()

%This function is used to define all parameters that a user can set. It
%returns a struct P with parameters, grouped into subfields, and a struct X
%with the initial guess at the state that is to be passed to fmincon.

%% Dynamics:
P.dyn.gravity = 9.81;  % (m/s^2)   Acceleration of gravity
P.dyn.length = 0.25;  % (m)   Length of Handle
P.dyn.mass = 0.2;   % (kg)  Hammer Point Mass
P.dyn.coeffRestitution = 0.2;  % () Coefficient of restitution

%% Optimization Algorithm:
% Type: optimset('fmincon') at command line to see all available options
P.opt = optimset(...
    'TolFun', 1e-8,...
    'TolX', 1e-12,...
    'TolCon', 1e-8,...
    'MaxIter',500,...
    'Display', 'iter-detailed',... % [ off | iter | iter-detailed | notify | notify-detailed | final | final-detailed ]
    'Algorithm', 'sqp',... % [ active-set | interior-point | interior-point-convex | levenberg-marquardt | sqp | trust-region-dogleg | trust-region-reflective ]
    'MaxFunEvals', 1e6,...
    'UseParallel', 'never'); % [ always | {never} ]

%% Discritization:
% This is the number of gridpoints along with trajectory
P.nGridPts = 20;

%% Cost Function
P.cost.Method = 0;
  % 0 = Time Integral of Torque^2
  % 1 = Time Integral of Torque*Rate

  % IF P.cost.Method == 1
        P.cost.Count_Negative_Work = 0.2;
            %Negative work is when the system is backdriving the motor.
             % 1 => negative work equal to positive work,
             % 0 => ignore negative work,
             %-1 => full regeneration
        P.cost.Motor_Cost_Smoothing = 1;
            %It turns out that smoothing is very important for convergence.
            %Here I use a modified form hyperbolic tangent smoothing:
            % 0    ==>  Sharp
            % inf  ==>  Smooth
            % 0.1  -->  feasible, non-optimal, hit iteration limit
            % 0.5  -->  close, but with some numerical artifacts
            % 1.0  -->  finds correct solution after 121 iterations
            % 5.0  -->  heavy smoothing - correct solution after 158
            %           iterations. Note that the cost function increases
            %           nearly linearly. This is because the optimization
            %           has found the minimum value of the smoothed cost
            %           function, which is nearly flat at this point.
  % END

%% Initial Guess:
%Initialize the trajectory by guessing a few rough points.
  P.init.durationGuess = 1;
  P.init.angleGuess = (pi/180)*[90, 0, 90];
  P.init.rateGuess = 2*pi*linspace(-1, 2, length(P.init.angleGuess));

  %P.init.torqueGuess = linspace(-1, 1, length(P.init.angleGuess));    %Start with values in toruq
  P.init.torqueGuess = zeros(size(P.init.angleGuess));    %start with empty torque

%% Bounds:
%States
    P.bnd.angle = (pi/180)*[-25,91];   %90 = Table, 0 = Vertical
    P.bnd.rate = 2*pi*[-2,2];
%Motor
    P.bnd.torque = [-1,1];
%Duration
    P.bnd.duration = [0.8,1.2];

%% Display:
% As this optimization runs, you can choose to have it display it's
% progress towards the solution. There are also options to automatically
% save the figures and generate the LaTeX code to put the figures into a
% document.

%Size the text on the figures nicely.
    P.disp.TitleFontSize = 16;
    P.disp.LabelFontSize = 16;
    P.disp.AxisFontSize = 12;
    P.disp.LineWidth = 3;
    P.disp.DefectLineWidth = 4;

%Decide if the plots should be shown as the code runs:
    P.disp.intermediatePlot = true;   %Use to toggle the plotting while running
    %IF P.disp.intermediatePlot, THEN
        P.disp.showIterationNum = ...%Creates intermediate plots at these iterations:
            [0,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987];
        P.disp.saveIntermediatePlots = false;   %Saves each auto-generated plot
        P.disp.createTex = true;   %Creates the skeleton of the LaTeX code to call figures
            P.disp.TexFileName = 'TexFigureCode.txt';
    %END

%% Constraints:
    P.cst.strikeAngle = pi/2;  %Hammer strikes exactly at this angle
    P.cst.strikeRate = max(P.bnd.rate);  %Reaches maximum speed at impact

%% Automatic:
% This section of code consists of things that are automatically done by
% the code, and are not parameters to be adjusted by the user.

% Later on in the code, it is necessary to pass the 'state' (all of the
% decision variables in the optimization problem) as a vector. This format
% is difficuly to work with, so there is a function to convert between a
% vector and a struct. Part of the conversion requires knowledge of the
% structure of the state struct. Here we create the initial guess at the
% state, and then use this as a template to get the information about the
% structure of this struct and store is in P.structDat.
    X = Initial_State(P);
    [~,P.structDat] = Convert_State(X);

% PhysicsIntegration uses persistent variables to keep track of some
% things. Before running the optimization, we need to reset these
% variables. One way of doing this is to clear the function from memory.
% Another way (not used here) is to pass a flag to the function (such as
% calling it with no arguments) and it can then reset the persistent
% variables itself.
    clear PhysicsIntegration;   %Reset all of the persistent variables

% Matlab takes a bit of time to generate a function handle with extra
% arguments (in this case P), so we do it here to save some time. NOTE -
% this means that any future changes to P will not be passed to the
% dynamics function. This is fine here, but could produce problems if the
% code was written differently.
    P.dyn.UserFunc = @(Z,u)HammerDynamics([],Z,u,P);

% This next block of code ensures that the LaTeX code fragmant starts with
% an empty file.
    if P.disp.createTex
        if exist(P.disp.TexFileName,'file')
            fid = fopen(P.disp.TexFileName,'w');
            fclose(fid);
        end
    end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [ x, obj_val ] = QCQP_UB( H_wave,y_wave,N,l,u,x1)
%UNTITLED9 Summary of this function goes here
%   Detailed explanation goes here
Q = 2*H_wave*H_wave.';
f = -2*H_wave*y_wave;
c = y_wave.'*y_wave;
A = zeros(N,2*N);
for ii = 1:N
    A(ii,ii) = cos((l(ii)+u(ii))/2);
    A(ii,ii+N) = sin((l(ii)+u(ii))/2); %Linear Constraint
end

function [y,grady] = quadobj(x,Q,f,c) %Objective Function
    y = 1/2*x.'*Q*x+f.'*x+c;
    if nargout>1
        grady = Q*x + f;
    end
end

function [y_con,yeq,grady,gradyeq] = quadconstr(x)    %Nonlinear Constraint
    for nn = 1:N
        yeq(nn) = x(nn)^2+x(nn+N)^2 - 1;
    end
    y_con = [];
    if nargout>2
       gradyeq = zeros(length(x),N);
       for nn = 1:N
           gradyeq(nn,nn) = 2*x(nn);
           gradyeq(nn+N,nn) = 2*x(nn+N);
       end
    end
    grady = [];
end

Acon = -A;

for ii = 1:N
    b(ii) = -cos((u(ii)-l(ii))/2);
end

% options = optimoptions(@fmincon,'Algorithm','sqp',...
%     'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
%     'Display','off','OptimalityTolerance',1e-6,'MaxIterations',400);

options = optimoptions('fmincon','Algorithm','sqp', 'GradObj','on');

fun = @(x)quadobj(x,Q,f,c);
nonlconstr = @(x)quadconstr(x);
% x1 = ones(2*N,1);
x = fmincon(fun,x1,Acon,b,[],[],[],[],nonlconstr,options);
obj_val = objval_func(x,H_wave,y_wave);
end

author: yuanhao-cui

repo: Must-Reading-on-ISAC

path: github-repos/MATLAB/yuanhao-cui-Must-Reading-on-ISAC/Must-Reading-on-ISAC-34cd6615c52ebca121428a979e756c608b195040/Codes/Fan2018TSP/Codes for DFRC Waveform Design/Constant Modulus/QCQP_UB.m


text:
classdef EdgeFunctionInterpolator < handle

    properties (Access = private)
        edgeMesh
        fInEdge
    end

    properties (Access = private)
        mesh
        fNodes
    end

    methods (Access = public)

        function obj = EdgeFunctionInterpolator(cParams)
            obj.init(cParams)
        end

        function fE = compute(obj)
            obj.createP1FunctionInEdges();
            fE = obj.interpolateInMiddleEdge();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.edgeMesh = cParams.edgeMesh;
            obj.fNodes   = cParams.fNodes;
        end

        function f = createP1FunctionInEdges(obj)
            s.mesh    = obj.edgeMesh;
            s.fValues = obj.fNodes;
            s.functionType = 'P1';
            f = P1Function(s);
            obj.fInEdge = f;
        end

        function fE = interpolateInMiddleEdge(obj)
            m = obj.edgeMesh;
            q = Quadrature.set(m.type);
            q.computeQuadrature('CONSTANT');
            xV = q.posgp;
            f  = obj.fInEdge;
            fE = squeeze(f.evaluate(xV));
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function checkreg = spm_cfg_checkreg
% SPM Configuration file for Check Reg
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_checkreg.m 6952 2016-11-25 16:03:13Z guillaume $

%--------------------------------------------------------------------------
% data Images to Display
%--------------------------------------------------------------------------
data         = cfg_files;
data.tag     = 'data';
data.name    = 'Images to Display';
data.help    = {'Images to display.'};
data.filter  = 'image';
data.ufilter = '.*';
data.num     = [1 24];

%--------------------------------------------------------------------------
% checkreg Check Registration
%--------------------------------------------------------------------------
checkreg      = cfg_exbranch;
checkreg.tag  = 'checkreg';
checkreg.name = 'Check Registration';
checkreg.val  = {data};
checkreg.help = {
    'Display of orthogonal views of one or more images.'
    ''
    'Clicking in any image moves the centre of the orthogonal views. Images are shown in orientations relative to that of the first selected image. The first specified image is shown at the top-left, and the last at the bottom right.  The fastest increment is in the left-to-right direction (the same as you are reading this).'
    ''
    'If you have put your images in the correct file format, then (possibly after specifying some rigid-body rotations):'
    '    The top-left image is coronal with the top (superior) of the head displayed at the top and the left shown on the left. This is as if the subject is viewed from behind.'
    '    The bottom-left image is axial with the front (anterior) of the head at the top and the left shown on the left. This is as if the subject is viewed from above.'
    '    The top-right image is sagittal with the front (anterior) of the head at the left and the top of the head shown at the top. This is as if the subject is viewed from the left.'
    }';
checkreg.prog = @check_reg;

%==========================================================================
function check_reg(job)
spm_check_registration(char(job.data));

author: spm


text:
function [out]=woody(x,varargin)
%
% [out]=woody(x,tol,max_it,est_mthd,xcorr_mthd)
%
% Weighted average using Woody average for a signal
% with jitter. Parameters:
%
% x             Signal measurements. Each COLUMN represents
%               and independent measure of the signal (or channel).
% tol           Tolerance paremeter to stop average (default is 0.1)
% max_it        Maximum number of iterations done on the average (default is 100).
% est_mthd      Estimation method to use. Options are:
%               'woody'     : classical approach (default)
%               'thornton'  : implements the Thornton approach that is also useful for different noise sources.
% xcorr_mthd    Determines what estimation method to use for the estimating the correlaation function using the
%               XCORR function. Options are:
%               'biased'   - scales the raw cross-correlation by 1/M.
%               'unbiased' - scales the raw correlation by 1/(M-abs(lags)). (Default)
% out           Final averaged waveform (time aligned).
%
%
%
% Written by Ikaro Silva
%
% Since 0.9.5
%
% %%%Example 1 %%%%
% t=[0:1/1000:1];
% N=1001;
% x=sin(2*pi*t)+sin(4*pi*t)+sin(8*pi*t);
% y=exp(0.01*[-1*[500:-1:1] 0 -1*[1:500]]);
% s=x.*y;
% sig1=0;
% sig2=0.1;
% M=100;
% S=zeros(N,M);
% center=501;
% TAU=round((rand(1,M)-0.5)*160);
% for i=1:M,
%     tau=TAU(i);
%
%     if(tau<0)
%         S(:,i)=[s(-1*tau:end)'; zeros(-1*(tau+1),1)];
%     else
%         S(:,i)=[zeros(tau,1);s(1:N-tau)'; ];
%     end
%     if(i<50)
%        S(:,i)=S(:,i) + randn(N,1).*sig1;
%     else
%         S(:,i)=S(:,i) + randn(N,1).*sig2;
%     end
% end
%
% [wood]=woody(S,[],[],'woody','biased');
% [thor]=woody(S,[],[],'thornton','biased');
% figure;
% subplot(211)
% plot(s,'b','LineWidth',2);grid on;hold on;plot(S,'r');plot(s,'b','LineWidth',2)
% legend('Signal','Measurements')
% subplot(212)
% plot(s);hold on;plot(mean(S,2),'r');plot(wood,'g');plot(thor,'k')
% legend('Signal','Normal Ave','Woody Ave','Thornton Ave');grid on

%endOfHelp
%Default parameter values
tol= 0.1;
max_it=100;
est_mthd='woody';
xcorr_mthd='unbiased';
thornton_sub=3;         %number of subaverages to use in the thornton procedure

if(nargin>1)
    if(~isempty(varargin{1}))
        tol=varargin{1};
    end
    if(nargin>2)
        if(~isempty(varargin{2}))
            max_it=varargin{2};
        end
        if(nargin>3)
            if(~isempty(varargin{3}))
                est_mthd=varargin{3};
            end
            if(nargin>4)
                if(~isempty(varargin{4}))
                    xcorr_mthd=varargin{4};
                end
            end
        end
    end
end

%Call repective averaging technique
switch est_mthd

    case 'woody'
        out=woody_core(x,tol,max_it,xcorr_mthd);

    case 'thornton'
        %Implement procedure from Thornton 2008
        [N,M]=size(x);
        K=floor(M/thornton_sub);

        %Call woody several times implementing the subaverages
        for k=1:K

            sub=thornton_sub*k;
            ind=round(linspace(1,M,sub+1));

            if((length(ind)-2) > (M/2))
                %Number of subaverages is equal to or just less than
                %half the number of trials, move to the final stage
                %and exit loop
                [out,est_lags]=woody_core(x,tol,max_it,xcorr_mthd);
                break
            end

            %Get woody average from the subaverages
            %procedure converges when there is no lag changes
            y=gen_subave(x,ind); %Generate sub averages
            y_old=y;
            err=1;
            while(err)
                [trash,est_lags]=woody_core(y,tol,max_it,xcorr_mthd);
                x=shift_data(x,est_lags,ind,N,M);
                y=gen_subave(x,ind); %Re-generate sub averages
                err=sum(abs(y(:)-y_old(:)));
                y_old=y;
            end

        end

    otherwise
        error(['Invalid option for est_mthd parameter: ' xcorr_mthd ' valid options are: woody, weighted, and thornton']);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%End of Maing Function%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%Helper Functions%%%%%%%%%%%%

function    x=shift_data(x,est_lags,ind,N,M)

%Shifts individual trials within each subaverage
K=length(est_lags);
for k=1:K
    lag=est_lags(k);
    if(lag)
        if(k~=K)
            sel_ind=[ind(k):ind(k+1)-1];
        else
            sel_ind=[ind(k):M];
        end
        pad=length(sel_ind);
        if(lag>0)
            x(:,sel_ind)=[zeros(lag-1,pad); x(lag:end,sel_ind)];
            %x(:,sel_ind)=[randn(lag-1,pad).*mean(std(x(:,sel_ind))).*0.001; x(lag:end,sel_ind)];
        elseif(lag<0)
            x(:,sel_ind)=[x(1:N+lag,sel_ind); zeros(lag*-1,pad)];
            %x(:,sel_ind)=[x(1:N+lag,sel_ind); randn(lag*-1,pad).*mean(std(x(:,sel_ind))).*0.001];
        end
    end
end

function [out,varargout]=woody_core(x,tol,max_it,xcorr_mthd)
[N,M]=size(x);
mx=mean(x,2);
p=zeros(N,1);
conv=1;
run=0;
sig_x=diag(sqrt(x'*x));
X=xcorr(mx);
ref=length(X)/2;
if(mod(ref,2))
    ref=ceil(ref);
else
    ref=floor(ref);
end

if(nargout>1)
    %In this case we output the lag of the trials as well
    lag_data=zeros(1,M);
end

while(conv*(run<max_it))

    z=zeros(N,1);
    w=ones(N,1);
    for i=1:M,

        y=x(:,i);
        xy=xcorr(mx,y,xcorr_mthd);
        [val,ind]=max(xy);
        if(ind>ref)
            lag=ref-ind-1;
        else
            lag=ref-ind;
        end
        if(lag>0)
            num=w(lag:end)-1;
            z(1:N-lag+1)=( z(1:N-lag+1).*num + y(lag:end))./w(lag:end);
            w(lag:end)=w(lag:end)+1;
        elseif(lag<0)
            num=w(lag*(-1)+1:end)-1;
            z(lag*(-1)+1:end)=( z(lag*(-1)+1:end).*num + y(1:N+lag) )./w(lag*(-1)+1:end);
            w(lag*(-1)+1:end)=w(lag*(-1)+1:end)+1;
        else
            z=z.*(w-1)./w + y./w;
            w=w+1;
        end
        if(exist('lag_data','var'))
            lag_data(i)=lag;
        end

    end

    old_mx=mx;
    mx=z;
    p_old=p;
    p=mx'*x./(sqrt(mx'*mx).*sig_x');
    p=sum(p)./M;
    err=abs(p-p_old);
    if(err<tol)
        conv=0;
    end
    run=run+1;

end

out=mx;

if(exist('lag_data','var'))
    varargout(1)={lag_data};
end

function [y]=gen_subave(x,ind)

[N,M]=size(x);
T=length(ind)-1;
y=zeros(N,T);

%Generate Subaverages
for i=1:T-1
    y(:,i)=mean(x(:,ind(i):ind(i+1)-1),2);
end

y(:,end)=mean(x(:,ind(T):end),2);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function [ni,canXform] = niftiApplyCannonicalXform(ni, canXform, phaseDir)
% Reorient NIFTI data and metadata fields into RAS space
%
%    [ni,canXform] = ...
%       niftiApplyCannonicalXform(ni, ...
%          [canXform=mrAnatComputeCannonicalXformFromDicomXform], phaseDir=[])
%
% Inputs:
%       ni:  Nifti-1 data structure.  (Should allow a filename in future).
% canXform:  Canonical transform from current form to RAS (Right-left,
%            Anterior-posterior, Superior-inferior).  If not passed in,
%            then the transform is estimated
% phaseDir:  ????  Phase encode direction.  (Columns, rows ... ???)
%
%  ** PLEASE EXPLAIN WHY WE NEED TO KNOW PHASE ENCODE DIRECTION ***
%
% Returns:
%       ni:  New nifti data set
% canXform:  Canonical transform
%
% Given a NIFTI-1 file ni- reorients the data and adjusts relevant metadata
% fields (like quaternion, pixfdim, freq/phase/slice dims, etc.) to be
% consistent.
%
% The transformation is specified by canXform (canonical transform). If the
% canXform isn't provided, it is estimated using
%
%     mrAnatComputeCannonicalXformFromDicomXform
%
% You can also pass the phase dir to overwrite a bad phase_dir field in
% your nifti header. Note that you should specify the phase_dir of the
% input data *before* the cannonical xform is applied. The value should be
% 1, 2, or 3 to specify phase-encoding along the first, second or third
% image dimension.  ** PLEASE EXPLAIN WHICH FUNCTION USES THIS INFORMATION
%
% See also: mrAnatComputeCannonicalXformFromDicomXform, applyCannonicalXform
%
% HISTORY:
% 2007.05.17 RFD: wrote it.
% 2008.08.15 RFD: fixed dim order bug.
%
% (c) Stanford VISTALAB

if(nargin<1), help(mfile); end

% Do a sanity-check on the nifti transform
% ni = niftiCheckQto(ni);
ni = niftiSet(ni,'Check Qto',''); %The final value is not necessary here

if(~exist('canXform','var') || isempty(canXform))
    canXform = mrAnatComputeCannonicalXformFromDicomXform(ni.qto_xyz, ni.dim(1:3));
end
%canXform(:,4) = [0 0 0 1]';
if(exist('phaseDir','var') && ~isempty(phaseDir))
    ni.phase_dim = phaseDir;
end

if(all(all(canXform == eye(4))))
 fprintf('[%s]: Data are already in canonical orientation.\n', mfilename)
else

    % Apply the transform
    [ni.data,newPixdim,dimOrder] = ...
        applyCannonicalXform(ni.data, canXform, ni.pixdim(1:3), false);

    % Fill the NIFTI image slots
    ni = niftiSet(ni,'dim',size(niftiGet(ni,'data')));
    ni.pixdim(1:3) = newPixdim;

    % Fill the transform slots
    ni = niftiSet(ni, 'qto', inv(canXform * niftiGet(ni,'qto_ijk')));
    if(any(ni.sto_xyz(:)>0))
        ni.sto_ijk = canXform*ni.sto_ijk;
        ni.sto_xyz = inv(ni.sto_ijk);
    end

    if(ni.freq_dim>0 && ni.freq_dim<4)
        ni.freq_dim = dimOrder(ni.freq_dim);
    else
        disp('freq_dim not set correctly in NIFTI header.');
    end

    if(ni.phase_dim>0 && ni.phase_dim<4)
        ni.phase_dim = dimOrder(ni.phase_dim);
    else
        disp('phase_dim not set correctly in NIFTI header.');
    end

    if(ni.slice_dim>0 && ni.slice_dim<4)
        ni.slice_dim = dimOrder(ni.slice_dim);
    end
end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f = langermannsfcn(x)

if strcmpi(x,'init')
    f.PopInitRange = [0;10] ;
    f.PopulationSize = 200 ;
elseif size(x,2) == 2
    x = reshape(x,1,[]) ;
    a = [3,5,2,1,7; 5,2,1,4,9] ;
    c = [1,2,5,2,3] ;
%     xa = (repmat(x,size(a,2),1)'-a)'*(repmat(x,size(a,2),1)'-a) ;
%     % xa is now 5x2
%     f = exp(-1/pi*xa)*cos(pi*xa);
%     f = 0 ;
    f = sum(c.*exp(-sum((repmat(x',1,size(a,2)) - a).^2,1)/pi).*...
        cos(pi*sum((repmat(x',1,size(a,2)) - a).^2,1)),2) ;
else
    error('Two dimensional input required for Langermann''s Function')
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [y,t,optw,W,C,confb95,yb] = sskernel(x,tin,W)
% [y,t,optw,W,C,confb95,yb] = sskernel(x,t,W)
%
% Function `sskernel' returns an optimized kernel density estimate
% using a Gauss kernel function.
%
% Examples:
% >> x = 0.5-0.5*log(rand(1,1e3)); t = linspace(0,3,1000);
% >> [y,t,optw] = sskernel(x,t);
% This example produces a vector of kernel density estimates, y, at points
% specified in a vector t, using an optimized bandwidth, optw (a standard
% deviation of a normal density function).
%
% >> sskernel(x);
% By calling the function without output arguments, the estimated density
% is displayed along with 95% bootstrap confidence intervals.
%
% Input arguments:
% x:    Sample data vector.
% tin (optinal):
%       Points at which estimation are computed. Please use fine resolution
%       to obtain a correct optimal bandwidth.
% W (optinal):
%       A vector of kernel bandwidths.
%       If W is provided, the optimal bandwidth is selected from the
%       elements of W.
%       * Do not search bandwidths smaller than a sampling resolution of data.
%       If W is not provided, the program searches the optimal bandwidth
%       using a golden section search method.
%
% Output arguments:
% y:    Estimated density
% t:    Points at which estimation was computed.
%       The same as tin if tin is provided.
%       (If the sampling resolution of tin is smaller than the sampling
%       resolution of the data, x, the estimation was done at smaller
%       number of points than t. The results, t and y, are obtained by
%       interpolating the low resolution sampling points.)
% optw: Optimal kernel bandwidth.
% W:    Kernel bandwidths examined.
% C:    Cost functions of W.
% conf95:
%       Bootstrap confidence intervals.
% yb:   Booststrap samples.
%
%
% Usage:
% >> [y,t,optw] = sskernel(x);
% When t is not given in the input arguments, i.e., the output argument t
% is generated automatically.
%
% >> W = linspace(0.01,1,20);
% >> [y,t,optw] = sskernel(x,t,W);
% The optimal bandwidth is selected from the elements of W.
%
% >> [y,t,optw] = sskernel(x,t,0.1);
% If the density estimate with a given bandwidth, simply put a scalar value
% as W. The computation is faster than the built-in function, ksdensity.
%
% >> [y,t,optw,confb95,yb] = sskernel(x);
% This additionally computes 95% bootstrap confidence intervals, confb95.
% The bootstrap samples are provided as yb.
%
%
% Optimization principle:
% The optimal bandwidth is obtained as a minimizer of the formula,
% sum_{i,j} \int k(x - x_i) k(x - x_j) dx  -  2 sum_{i~=j} k(x_i - x_j),
% where k(x) is the kernel function, according to
%
% Hideaki Shimazaki and Shigeru Shinomoto
% Kernel Bandwidth Optimization in Spike Rate Estimation
% Journal of Computational Neuroscience 2010
% http://dx.doi.org/10.1007/s10827-009-0180-4
%
% The above optimization is based on a principle of minimizing
% expected L2 loss function between the kernel estimate and an unknown
% underlying density function. An assumption is merely that samples
% are drawn from the density independently each other.
%
% For more information, please visit
% http://2000.jukuin.keio.ac.jp/shimazaki/res/kernel.html
%
% See also SSVKERNEL, SSHIST
%
% Bug fix
% 131004 fixed a problem for large values
%
% Hideaki Shimazaki
% http://2000.jukuin.keio.ac.jp/shimazaki

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters Settings
x = reshape(x,1,numel(x));

if nargin == 1
    T = max(x) - min(x);
    [mbuf,nbuf,dt_samp] = find( sort(diff(sort(x))),1,'first');
    tin = linspace(min(x),max(x), min(ceil(T/dt_samp),1e3));
    t = tin;
    x_ab = x( logical((x >= min(tin)) .*(x <= max(tin))) ) ;
else
    T = max(tin) - min(tin);
    x_ab = x( logical((x >= min(tin)) .*(x <= max(tin))) ) ;
    [mbuf,nbuf,dt_samp] = find( sort(diff(sort(x_ab))),1,'first');

    if dt_samp > min(diff(tin))
        t = linspace(min(tin),max(tin), min(ceil(T/dt_samp),1e3));
    else
        t = tin;
    end
end

dt = min(diff(t));

% Create a finest histogram
y_hist = histc(x_ab,t-dt/2);
L = length(y_hist);
N = sum(y_hist);
y_hist = y_hist/N/dt;   %density

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute a Cost Function

if nargin >= 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global search
C = zeros(1,length(W));
C_min = Inf;

for k = 1: length(W)
	w = W(k);
    [C(k) yh] = CostFunction(y_hist,N,w,dt);

    if C(k) < C_min
        C_min = C(k);
        optw = w;
        y = yh;
    end
end

else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Golden section search on a log-exp scale
% Initialize
Wmin = 2*dt; Wmax = 1*(max(x) - min(x));

tol = 10^-5;
phi = (sqrt(5) + 1)/2;        %golden ratio
%logexp = @(x) log(1+exp(x));
%ilogexp = @(x) log(exp(x)-1);

%a = Wmin; b = Wmax;
a = ilogexp(Wmin); b = ilogexp(Wmax);

c1 = (phi-1)*a + (2-phi)*b;
c2 = (2-phi)*a + (phi-1)*b;

f1 = CostFunction(y_hist,N,logexp(c1),dt);
f2 = CostFunction(y_hist,N,logexp(c2),dt);

k = 1;
while abs(b-a) > tol*(abs(c1)+abs(c2)) && k <= 20
	if (f1 < f2)
        b = c2;
        c2 = c1;

        c1 = (phi - 1)*a + (2 - phi)*b;

        f2 = f1;
        [f1 yh1] = CostFunction(y_hist,N,logexp(c1),dt);

        W(k) = logexp(c1);
        C(k) = f1;
        optw = logexp(c1);
        y = yh1./sum(yh1.*dt);  %make the final output a density
    else
        a = c1;
        c1 = c2;

        c2 = (2 - phi)*a + (phi - 1)*b;

        f1 = f2;
        [f2 yh2] = CostFunction(y_hist,N,logexp(c2),dt);

        W(k) = logexp(c2);
        C(k) = f2;
        optw = logexp(c2);
        y = yh2./sum(yh2.*dt);
    end

    k = k + 1;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bootstrap Confidence Intervals
if nargout == 0 || nargout >= 6
    nbs = 1e3;        %number of bootstrap samples
    yb = zeros(nbs,length(tin));

    for i = 1: nbs,
        %y_histb = poissrnd(y_hist*dt*N)/dt/N;

        idx = ceil(rand(1,N)*N);
        xb = x_ab(idx);
        y_histb = histc(xb,t-dt/2)/dt/N;

        yb_buf = fftkernel(y_histb,optw/dt);
        yb_buf = yb_buf / sum(yb_buf*dt);

        yb(i,:) = interp1(t,yb_buf,tin);
    end

    ybsort = sort(yb);
    y95b = ybsort(floor(0.05*nbs),:);
    y95u = ybsort(floor(0.95*nbs),:);

    confb95 = [y95b; y95u];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Return results
y = interp1(t,y,tin);
t = tin;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results
if nargout == 0
            hold on;

            line([t; t],[y95b; y95u]...
                    ,'Color',[7 7 7]/8,'LineWidth',1 );
            plot(t,y95b,'Color',[7 7 7]/9,'LineWidth',1);
            plot(t,y95u,'Color',[7 7 7]/9,'LineWidth',1);

            plot(t,y,'Color',[0.9 0.2 0.2],'LineWidth',2);

            grid on;
            ylabel('density');
            set(gca,'TickDir','out');
else
    if nargin >= 4
        if strcmp(option,'Visible')
            hold on;

            if nargout >= 6
                line([t; t],[y95b; y95u]...
                    ,'Color',[7 7 7]/8,'LineWidth',1 );
                plot(t,y95b,'Color',[7 7 7]/9,'LineWidth',1);
                plot(t,y95u,'Color',[7 7 7]/9,'LineWidth',1);
            end

            plot(t,y,'Color',[0.9 0.2 0.2],'LineWidth',1);

            grid on;
            ylabel('density');
            set(gca,'TickDir','out');

        end
    end
end

function [C yh] = CostFunction(y_hist,N,w,dt)
yh = fftkernel(y_hist,w/dt);  %density

%formula for density
C = sum(yh.^2)*dt - 2* sum(yh.*y_hist)*dt...
        + 2*1/sqrt(2*pi)/w/N;
C = C * N* N;

%formula for rate
%C = dt*sum( yh.^2 - 2*yh.*y_hist + 2/sqrt(2*pi)/w*y_hist );

function y = fftkernel(x,w)
% y = fftkernel(x,w)
%
% Function `fftkernel' applies the Gauss kernel smoother to
% an input signal using FFT algorithm.
%
% Input argument
% x:    Sample signal vector.
% w: 	Kernel bandwidth (the standard deviation) in unit of
%       the sampling resolution of x.
%
% Output argument
% y: 	Smoothed signal.
%
% MAY 5/23, 2012 Author Hideaki Shimazaki
% RIKEN Brain Science Insitute
% http://2000.jukuin.keio.ac.jp/shimazaki

L = length(x);
Lmax = max(1:L+3*w);
n = 2^(nextpow2(Lmax));

X = fft(x,n);

f = (0:n-1)/n;
f = [-f(1:n/2+1) f(n/2:-1:2)];

K = exp(-0.5*(w*2*pi*f).^2);

y = ifft(X.*K,n);

y = y(1:L);

function y = logexp(x)
if x<1e2
    y = log(1+exp(x));
else
    y = x;
end

function y = ilogexp(x)
%ilogexp = @(x) log(exp(x)-1);
if x<1e2
    y = log(exp(x)-1);
else
    y = x;
end

repo: REST