text:
function varargout = ARMAX_viewer(varargin)
% ARMAX_viewer M-file for ARMAX_viewer.fig
%      ARMAX_viewer, by itself, creates a new ARMAX_viewer or raises the existing
%      singleton*.
%
%      H = ARMAX_viewer returns the handle to a new ARMAX_viewer or the handle to
%      the existing singleton*.
%
%      ARMAX_viewer('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ARMAX_viewer.M with the given input arguments.
%
%      ARMAX_viewer('Property','Value',...) creates a new ARMAX_viewer or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ARMAX_viewer_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ARMAX_viewer_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

% TODO : Add parameter display for hyperparapmeter for distributions
% FIXME: Fix the behavior for very large models

% Edit the above text to modify the response to help ARMAX_viewer

% Last Modified by GUIDE v2.5 09-Sep-2006 14:45:09

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @ARMAX_viewer_OpeningFcn, ...
    'gui_OutputFcn',  @ARMAX_viewer_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 ARMAX_viewer is made visible.
function ARMAX_viewer_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 ARMAX_viewer (see VARARGIN)

% Choose default command line output for ARMAX_viewer
handles.output  = hObject;
handles.Results = varargin{1};
if length(varargin)==2
    handles.DisplayNumber = varargin{2};
else
    handles.DisplayNumber = length(handles.Results);
end

% Update handles structure
guidata(hObject, handles);
popup_cell=[];
for i=1:length(handles.Results);
    popup_cell{i}=handles.Results{i}.StringID;
end
set(handles.model_popupmenu,'String',popup_cell)
set(handles.model_popupmenu,'Value',length(handles.Results))
set(gcf,'Position',[10 10 1000 706])
%
display_model(handles,hObject);

% UIWAIT makes ARMAX_viewer wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = ARMAX_viewer_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} = handles.output;

% --- Executes on mouse press over axes background.
function model_axes_ButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to model_axes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function generate_text(handles,hObject);
% Get rid of results not needed
number = handles.DisplayNumber;
results=handles.Results{number};
% Compute the needed strings
[ParameterString,VariableString]= generate_strings(results);
% Layout the model
[model_layout_info]=model_layout(ParameterString,VariableString,handles);
% Layout the model statistics
[model_stats_layout_info]=model_stats_layout(results,handles);
% Resize the Parameter Estimates Window
[estimate_window]=resize_parameter_estimate_uipanel(handles);
% Layout the column titles and return their position
[title_positions]=title_layout(handles);
% Compute the text elements for the main display and their positions
[text_positions]=compute_text_positions(handles,ParameterString,title_positions);
% Update the results with a flag that the text layout is done and with all
% of the relevant information, also set the handles
handles.Results{number}.ParameterString=ParameterString;
handles.Results{number}.VariableString=VariableString;
handles.Results{number}.model_layout_info=model_layout_info;
handles.Results{number}.model_stats_layout_info=model_stats_layout_info;
handles.Results{number}.estimate_window=estimate_window;
handles.Results{number}.title_position=title_positions;
handles.Results{number}.text_position=text_positions;
handles.Results{number}.TextFinished=true;
guidata(hObject, handles);

function [ParameterString,VariableString]= generate_strings(results);
% Not finished
% TO DO: ARARCH, EGARCH and AGARCH
index=1;

ParameterString{index}='y_t =';
c{index}='';
index=index+1;

if results.constant
    ParameterString{index}='\phi_0';
    VariableString{index}='';
    index=index+1;
end

if ~isempty(results.ARlags)
    if max(results.ARlags)>0
        for i=1:length(results.ARlags)
            ParameterString{index}=['\phi_{' num2str(results.ARlags(i)) '}'];
            VariableString{index}=['y_{t-' num2str(results.ARlags(i)) '}'];
            index=index+1;
        end
    end
end

if ~isempty(results.MAlags)
    if max(results.MAlags)>0
        for i=1:length(results.MAlags)
            ParameterString{index}=['\theta_{' num2str(results.MAlags(i)) '}'];
            VariableString{index}=['\epsilon_{t-' num2str(results.MAlags(i)) '}'];
            index=index+1;
        end
    end
end

ParameterString{index}='';
VariableString{index}='\epsilon_t';

function [layout_info]=model_layout(ParameterString,VariableString,handles)

axes(handles.model_axes);
%set(gcf,'Position',[ -3          35        1024         664],'Units','pixels')
%set(ax,'Units','pixels','Position',[134.1200   500.0400  793.6000  120])
line=0.75;
count=1;

result = handles.Results{handles.DisplayNumber};
ModelString = result.StringID;
% Add on the model specific string
final_str{count}=ModelString;
count=count+1;

%FIX ME, need to worry about no constant!
if result.constant
    final_str{count}=[ParameterString{1} ParameterString{2}];
else
    final_str{count}=[ParameterString{1} ParameterString{2} VariableString{2}];
end

xpos=.05;
for i=3:length(ParameterString)
    t=text(xpos,line,['$' final_str{count} '+' ParameterString{i} VariableString{i} '$' ],'Interpreter','latex','FontSize',12);
    ex=get(t,'Extent');
    delete(t);
    if ex(1)+ex(3)<1
        % No problems, add and try next
        final_str{count}=[final_str{count} '+' ParameterString{i} VariableString{i}];
    else
        %Problem, no not add and start a new line
        xpos=.1;
        count=count+1;
        final_str{count}=['+' ParameterString{i} VariableString{i}];
    end
end

t=[];
line=0.75;
for i=1:count;
    if i<3
        xpos=.02;
    else
        xpos=.07;
    end
    if i>1
        t{i}=text(xpos,line,['$' final_str{i} '$'] ,'Interpreter','latex','FontSize',12);
    else
        t{i}=text(xpos,line,['$' final_str{i} '$'] ,'FontSize',12);
    end
    line=line-0.25;
end

% Now that the text is layed out, the axes and uipanel need to be resized
% accordingly and the text strings need to be vertically repositioned
% axes height is 33* no lines, y position is always 9
% uipanel height is 33*no lines+24, y position is 650-uiheight
no_lines=length(final_str);
uiHeight = no_lines*33+24;
uiY = 670 - uiHeight;
uiP=get(handles.model_uipanel,'Position');
uiP(2)=uiY;
uiP(4)=uiHeight;
set(handles.model_uipanel,'Position',uiP);

axHeight = 33*no_lines;
axY = 9;
axP=get(handles.model_axes,'Position');
axP(2)=axY;
axP(4)=axHeight;
set(handles.model_axes,'Position',axP);

% Finally uniformly distribute the t{i} between at
% [1:no_lines]/(no_lines+1)
positions = (no_lines:-1:1)/(no_lines+1);
for i=1:no_lines;
    tP=get(t{i},'Position');
    tP(2)=positions(i);
    set(t{i},'Position',tP);
end

% Need to return everything needed to make relayout easy:
%   uipanle position
%   axes position
%   text strings
%   text string position
layout_info.model_uipanel_position = get(handles.model_uipanel,'Position');
layout_info.model_axes_position = get(handles.model_axes,'Position');
layout_info.model_text_strings = final_str;
for i=1:length(final_str);
    layout_info.model_text_strings_position{i}=get(t{i},'Position');
end
% Clear the axis
cla

function [layout_info]=model_stats_layout(results,handles)
% 1. Fix the position of and size of the UI and axes
% 2. Layout the important stats and their labels

axes(handles.stats_axes);
set(handles.stats_axes,'Position',[15     9   920   60])
model_uipanel_position=get(handles.model_uipanel ,'Position');
set(handles.model_stats_uipanel,'Position',[25  model_uipanel_position(2)-84    950   84])

titles={'Log-likelihood','No. of Params',' Akaike IC','Schwartz/Bayesian IC','$\hat{\sigma}^2$'};
positions=linspace(.02,.9,length(titles));
total_space=0;
for i=1:length(titles)
    if i==5
        t(i)=text(positions(i),.66,titles{i},'FontName','Tahoma','FontSize',10,'Interpreter','latex');
    else
        t(i)=text(positions(i),.66,titles{i},'FontName','Tahoma','FontSize',10);
    end
    ex(i,:)=get(t(i),'Extent');
    if i>1
        total_space=total_space+ex(i,1)-(ex(i-1,3)+ex(i-1,1));
    end
end
total_space=total_space+.98-(ex(5,3)+ex(5,1));
% Compute the total amount of space between the objects and uniformly
% respace them
spaces=[4 4 4 4]*total_space/22;
for i=2:length(titles)
    ex=get(t(i-1),'Extent');
    set(t(i),'Position',[ex(1)+ex(3)+spaces(i-1) .66 0]);
end

% Set up the text strings.  I need to limit the width
data{1}=sprintf('%0.7g',results.ll);
data{2}=num2str(results.K);
data{3}=sprintf('%0.7g',results.AIC);
data{4}=sprintf('%0.7g',results.BIC);
data{5}=sprintf('%0.7g',results.seregression);

for i=1:5
    d(i)=text(positions(i),.33,data{i},'FontName','Tahoma','FontSize',10);
    % Now move it so that it is right aligned
    p1 = get(d(i),'Extent');
    p2 = get(t(i),'Extent');
    set(d(i),'Position',[p2(1)+p2(3)-p1(3) .33 0]);
end
AX=axis;
hold on;
for i=1:5
    % Now move it so that it is right aligned
    p2 = get(t(i),'Extent');
    axis(AX);
    x_line_data{i}=[p2(1) p2(1)+p2(3)]';
    y_line_data{i}=[p2(2) p2(2)];
    plot([p2(1) p2(1)+p2(3)],[p2(2) p2(2)],'k');
end
hold off;

% Need to return all of teh strings along with all of their positions
layout_info.Titles=titles;
layout_info.Data=data;
layout_info.TitlePosition=[];
layout_info.DataPosition=[];

for i=1:5;
    layout_info.TitleHandles{i}=t(i);
    layout_info.TitlePosition{i}=get(t(i),'Position');
    layout_info.DataHandle{i}=d(i);
    layout_info.DataPosition{i}=get(d(i),'Position');
end
layout_info.model_stats_uipanel_position = get(handles.model_stats_uipanel,'Position');
layout_info.stats_axes_position = get(handles.stats_axes,'Position');
layout_info.x_line_data=x_line_data;
layout_info.y_line_data=y_line_data;
% Clear the axis
cla

% --- Executes on selection change in model_popupmenu.
function model_popupmenu_Callback(hObject, eventdata, handles)
% hObject    handle to model_popupmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns model_popupmenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from model_popupmenu
handles.DisplayNumber = get(hObject,'Value');
guidata(hObject, handles);
display_model(handles,hObject);

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

% Hint: popupmenu 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

function page_number_edit_Callback(hObject, eventdata, handles)
% hObject    handle to page_number_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 page_number_edit as text
%        str2double(get(hObject,'String')) returns contents of page_number_edit as a double
val = get(hObject,'String');
result = handles.Results{handles.DisplayNumber};
LastPage = handles.LastPage;
NumberOfPages = size(result.text_position.x_pos,3);
[CurrentPage, OK] = bounded_integer_validate(val,LastPage,1,NumberOfPages);
if ~OK
    set(hObject,'String',num2str(CurrentPage));
end
draw_parameter_table(result,CurrentPage,handles)
handles.LastPage=CurrentPage;
% Update handles structure
guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function page_number_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to page_number_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

function num_pages_text_Callback(hObject, eventdata, handles)
% hObject    handle to num_pages_text (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 num_pages_text as text
%        str2double(get(hObject,'String')) returns contents of num_pages_text as a double

% --- Executes during object creation, after setting all properties.
function num_pages_text_CreateFcn(hObject, eventdata, handles)
% hObject    handle to num_pages_text (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

% --- Executes on button press in increase_page_pushbutton.
function increase_page_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to increase_page_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
result = handles.Results{handles.DisplayNumber};
NumberOfPages = size(result.text_position.x_pos,3);
CurrentPage = str2double(get(handles.page_number_edit,'String'));
if CurrentPage < NumberOfPages
    CurrentPage = CurrentPage +1;
    set(handles.page_number_edit,'String',num2str(CurrentPage));

    draw_parameter_table(result,CurrentPage,handles)
    handles.LastPage=CurrentPage;
    % Update handles structure
    guidata(hObject, handles);

end

% --- Executes on button press in decrease_page_pushbutton.
function decrease_page_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to decrease_page_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
result = handles.Results{handles.DisplayNumber};
CurrentPage = str2double(get(handles.page_number_edit,'String'));
if CurrentPage >1
    CurrentPage = CurrentPage - 1;
    set(handles.page_number_edit,'String',num2str(CurrentPage));

    draw_parameter_table(result,CurrentPage,handles)
    handles.LastPage=CurrentPage;
    % Update handles structure
    guidata(hObject, handles);
end

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

function [window_stats]=resize_parameter_estimate_uipanel(handles)
model_uipanel_pos=get(handles.model_stats_uipanel,'Position');
estimates_uipanel_pos=get(handles.estimates_uipanel ,'Position');
estimates_uipanel_pos(4)=estimates_uipanel_pos(4)+model_uipanel_pos(2)-sum(estimates_uipanel_pos([2 4]));
set(handles.estimates_uipanel ,'Position',estimates_uipanel_pos);

estimates_axes_pos=get(handles.estimates_axes,'Position');
estimates_axes_pos(2)=12;
set(handles.estimates_axes,'Position',estimates_axes_pos)

estimates_uipanel_pos=get(handles.estimates_uipanel ,'Position');
estimates_axes_pos=get(handles.estimates_axes,'Position');
estimates_axes_pos(4)=estimates_uipanel_pos(4)-24;
set(handles.estimates_axes,'Position',estimates_axes_pos);
estimates_axes_pos=get(handles.estimates_axes,'Position');
% Note, I want ot make sure that the estimates window is at least 330
% pixels high

% FIXME: Fix the behavior for very large models
% ScreenSize = get(0,'ScreenSize')
% set(handles.figure1,'Position',[1 ScreenSize-HeightofWindow-33 WidthofWindow HeightofWindow])
% Return the position of the uipanel, the axes and the figure!
window_stats.axes_pos=estimates_axes_pos;
window_stats.uipanel_pos=estimates_uipanel_pos;
window_stats.figure_pos=get(handles.figure1,'Position');

function [title_position]=title_layout(handles)
titles={'Parameter','Estimate','Std. Error','T-stat','P-val'};
axes(handles.estimates_axes);
t(1)=text(.02,.5,titles{1},'FontName','Tahoma','FontSize',10);
axes_pos=get(handles.estimates_axes,'Position');
set(t(1),'Units','pixels');
t_pos=get(t(1),'Position');
t_ex = get(t(1),'Extent');
y=floor(axes_pos(4)-t_ex(4)/2-10);
t_pos(2)=y;
set(t(1),'Position',t_pos)
position{1}=t_pos;
for i=2:length(titles);
    t(i)=text((i-1)/6+.02,.5,titles{i},'FontName','Tahoma','FontSize',10);
    set(t(i),'Units','pixels') ;
    t_pos=get(t(i),'Position');
    t_pos(2)=y;
    set(t(i),'Position',t_pos);
    position{i}=t_pos;
end
AX=axis;
hold on;
for i=1:length(titles)
    % Now move it so that it is right aligned
    set(t(i),'Units','normalized');
    p2 = get(t(i),'Extent');
    % Need to get the right points of each label to do right alignment
    right_points{i}=p2(1)+p2(3);
    axis(AX);
    x_plot{i}=[p2(1) p2(1)+p2(3)];
    y_plot{i}=[p2(2) p2(2)]';
    plot(x_plot{i},y_plot{i},'k');
end

hold off;
tt=text(.5,.5,'$\alpha \beta \gamma \lambda \omega \nu \mu$','Interpreter','latex','FontSize',14);
row_height=get(tt,'Extent');
delete(tt);
row_height=row_height(4);
title_position.Titles=titles;
title_position.TitlePosition=position;
title_position.x_line_data=x_plot;
title_position.y_line_data=y_plot;
title_position.RowHeight=row_height;
title_position.NumRows=floor((1-0.5*row_height)/row_height);
title_position.RightPoints=right_points;
% Clear the axis
cla

function [text_positions]=compute_text_positions(handles,str,title_position)
% First set up the texts

results = handles.Results{handles.DisplayNumber};
right_points = title_position.RightPoints;
actual_rows=title_position.NumRows-1;
for i=2:length(str)-1
    % first compute the page
    page = ceil((i-1)/actual_rows);

    if mod(i-1,actual_rows)==0
        table_position = actual_rows;
    else
        table_position = mod(i-1,actual_rows);
    end
    for j=1:5
        switch j
            case 1
                % Latex string
                window_text{table_position,j,page} = ['$' str{i} '$'];
            case 2
                % Estimate
                temp = results.parameters(i-1);
            case 3
                % Estimate
                temp = results.StdErr(i-1);
            case 4
                % Estimate
                temp = results.Tstat(i-1);
            case 5
                % Estimate
                temp = results.Pval(i-1);
        end
        if j~=1
            if temp<.000001
                window_text{table_position,j,page} = sprintf('%0.3g',temp);
            else
                window_text{table_position,j,page} = sprintf('%0.5f',temp);
            end
        end
    end
    if mod(i-1,actual_rows)==0
        actual_position = actual_rows+1;
    else
        actual_position = mod(i-1,actual_rows)+1;
    end
    row_height=title_position.RowHeight;
    for j=1:5
        y_pos{table_position,j,page} = 1-(actual_position)*row_height;
        if j==1
            temp_text =text(.5,y_pos{table_position,j,page},window_text{table_position,j,page},'Interpreter','latex','FontSize',12);
        else
            temp_text =text(.5,y_pos{table_position,j,page},window_text{table_position,j,page},'FontName','Tahoma','FontSize',10);
        end
        ex=get(temp_text,'Extent');
        delete(temp_text);
        x_pos{table_position,j,page}=right_points{j}-ex(3);
    end
end
% Set up the return data
text_positions.x_pos = x_pos;
text_positions.y_pos = y_pos;
text_positions.window_text = window_text;

function draw_data(result,handles,hObject)
% Resize the figure
set(handles.figure1,'Position',result.estimate_window.figure_pos);
% Resize the uipanels
set(handles.model_uipanel,'Position',result.model_layout_info.model_uipanel_position);
set(handles.model_stats_uipanel,'Position',result.model_stats_layout_info.model_stats_uipanel_position);
set(handles.estimates_uipanel,'Position',result.estimate_window.uipanel_pos)
% Resize the axes
set(handles.model_axes,'Position',result.model_layout_info.model_axes_position);
set(handles.stats_axes ,'Position',result.model_stats_layout_info.stats_axes_position);
set(handles.estimates_axes,'Position',result.estimate_window.axes_pos);
% Draw the model
axes(handles.model_axes)
% Make sure itis clear
cla
str = result.model_layout_info.model_text_strings;
pos = result.model_layout_info.model_text_strings_position;
N = length(str);
for i=1:N
    x=pos{i}(1);
    y=pos{i}(2);
    if i>1
        text(x,y,['$' str{i} '$'],'Unit','normalized','interpreter','latex','FontSize',12)
    else
        text(x,y,[str{i}],'Unit','normalized','FontSize',12)
    end
end
% Draw the model statistic titles, lines, and values
axes(handles.stats_axes)
% Make sure itis clear
cla
str = result.model_stats_layout_info.Titles;
pos = result.model_stats_layout_info.TitlePosition;
N = length(str);
for i=1:N
    x=pos{i}(1);
    y=pos{i}(2);
    if i~=5
        text(x,y,str{i},'Unit','normalized','FontSize',10,'FontName','Tahoma')
    else
        text(x,y,str{i},'Unit','normalized','FontSize',10,'FontName','Tahoma','Interpreter','latex')
    end
end
hold on;
for i=1:N
    % Now move it so that it is right aligned
    x=result.model_stats_layout_info.x_line_data{i};
    y=result.model_stats_layout_info.y_line_data{i};
    plot(x,y,'k');
end
hold off;
str = result.model_stats_layout_info.Data;
pos = result.model_stats_layout_info.DataPosition;
N = length(str);
for i=1:N
    x=pos{i}(1);
    y=pos{i}(2);
    text(x,y,str{i},'Unit','normalized','FontSize',10,'FontName','Tahoma')
end
% Call the function to draw the correct page, presumably one
draw_parameter_table(result,1,handles) % In the initial draw we will always draw page 1

handles.LastPage=1;
% Update handles structure
guidata(hObject, handles);

% Finally update the UI elements that have to do with the pages
NumberOfPages = size(result.text_position.x_pos,3);
if NumberOfPages==1
    % Easy case
    set(handles.page_number_edit,'BackgroundColor',[0.5 0.5 0.5],'String','1')
    set(handles.increase_page_pushbutton,'Enable','off')
    set(handles.decrease_page_pushbutton,'Enable','off')
    set(handles.number_of_pages_text,'String',num2str(NumberOfPages));
else
    set(handles.page_number_edit,'BackgroundColor',[1 1 1],'String','1')
    set(handles.increase_page_pushbutton,'Enable','on')
    set(handles.decrease_page_pushbutton,'Enable','on')
    set(handles.number_of_pages_text,'String',num2str(NumberOfPages));
end

function draw_parameter_table(result,page,handles)
axes(handles.estimates_axes)
% Make sure it is clear
cla
% Draw the column titles
str = result.title_position.Titles;
pos = result.title_position.TitlePosition;
N = length(str);
for i=1:N
    x=pos{i}(1);
    y=pos{i}(2);
    text(x,y,str{i},'Unit','pixel','FontSize',10,'FontName','Tahoma')
end
% Draw the lines
hold on;
for i=1:N
    % Now move it so that it is right aligned
    x=result.title_position.x_line_data{i};
    y=result.title_position.y_line_data{i};
    plot(x,y,'k');
end
hold off;
x = result.text_position.x_pos(:,:,page);
y = result.text_position.y_pos(:,:,page);
str = result.text_position.window_text(:,:,page);

N = size(x,1);

for i=1:N
    for j=1:5
        % Draw the data on the correct page
        if j==1
            text(x{i,j},y{i,j},str{i,j},'Interpreter','latex','FontSize',12);
        else
            text(x{i,j},y{i,j},str{i,j},'FontName','Tahoma','FontSize',10);
        end
    end
end

function display_model(handles,hObject)
% Check to see if the moodel has been processed.  If it has, go striaght to
% the draw model, else compute te strings then draw the model
results = handles.Results{handles.DisplayNumber};
if isfield(results,'TextFinished')
    draw_data(results,handles,hObject)
else
    % Generate the text labels
    generate_text(handles,hObject);
    % update handles
    handles = guidata(hObject);
    results = handles.Results{handles.DisplayNumber};
    % Draw the labels
    draw_data(results,handles,hObject)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to do integer vallidation for ACF Lags
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [val, OK] = bounded_integer_validate(val,orig_val,LB,UB)
val = str2double(val);
% Test cases
%   Nan: return old value and  OK=0
%   <0 return 0 and OK=0
%   + but not integer, floor and OK=0
%   +, integer, OK =1
OK = false;
if isnan(val)
    val = orig_val;
elseif val<LB
    val = LB;
elseif val>UB;
    val = UB;
elseif floor(val)~=val
    val=floor(val);
else
    OK = true;
end

author: bashtage

repo: mfe-toolbox


text:
function test_failed = test_libltfat_dgtreal_fb(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

[~,~,enuminfo]=libltfatprotofile;
phaseconv = enuminfo.ltfat_phaseconvention;

fftwflags = struct('FFTW_MEASURE',0,'FFTW_ESTIMATE',64,'FFTW_PATIENT',32,'FFTW_DESTROY_INPUT',1,...
    'FFTW_UNALIGNED',2,'FFTW_EXHAUSTIVE',8,'FFTW_PRESERVE_INPUT',16);

Larr  = [350 350   9   1];
glarr = [ 20  10   9   1];
aarr  = [ 10  10   9   1];
Marr  = [ 35  35   3   1];
Warr  = [  1   3   3   1];

for idx = 1:numel(Larr)
    L = Larr(idx);
    W = Warr(idx);
    a = aarr(idx);
    M = Marr(idx);
    M2 = floor(M/2) + 1;
    gl = glarr(idx);

    N = L/a;

    g = randn(gl,1,flags.complexity);
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    f = randn(L,W,flags.complexity);

    fPtr = libpointer(dataPtr,f);
    gPtr = libpointer(dataPtr,g);
    coutPtr = libpointer(dataPtr,cout);

    truec = dgtreal(f,g,a,M);

    funname = makelibraryname('dgtreal_fb',flags.complexity,0);
    status = calllib('libltfat',funname,fPtr,gPtr,L,gl,W,a,M,phaseconv.LTFAT_FREQINV,coutPtr);

    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
    fprintf(['DGTREAL FREQINV    L:%3i, gl:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,gl,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    truec = dgtreal(f,g,a,M,'timeinv');
    status = calllib('libltfat',funname,fPtr,gPtr,L,gl,W,a,M,phaseconv.LTFAT_TIMEINV,coutPtr);

    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+status,test_failed);
    fprintf(['DGTREAL TIMEINV    L:%3i, gl:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,gl,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    % With plan
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('dgtreal_fb_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,gPtr,gl,a,M,phaseconv.LTFAT_FREQINV,fftwflags.FFTW_MEASURE,plan);

    funname = makelibraryname('dgtreal_fb_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan, fPtr,L,W,coutPtr);

    funname = makelibraryname('dgtreal_fb_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    truec = dgtreal(f,g,a,M);
    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['DGTREAL FREQINV WP L:%3i, gl:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,gl,W,a,M,flags.complexity,ltfatstatusstring(status),fail);

    %%%%%%
    c = cast(randn(M2,N,W)+1i*randn(M2,N,W),flags.complexity);
    cout = complex2interleaved(c);
    coutPtr = libpointer(dataPtr,cout);

    plan = libpointer();
    funname = makelibraryname('dgtreal_fb_init',flags.complexity,0);
    statusInit = calllib('libltfat',funname,gPtr,gl,a,M,phaseconv.LTFAT_TIMEINV,fftwflags.FFTW_MEASURE,plan);

    funname = makelibraryname('dgtreal_fb_execute',flags.complexity,0);
    statusExecute = calllib('libltfat',funname,plan, fPtr,L,W,coutPtr);

    funname = makelibraryname('dgtreal_fb_done',flags.complexity,0);
    statusDone = calllib('libltfat',funname,plan);

    truec = dgtreal(f,g,a,M,'timeinv');
    res = norm(reshape(truec,M2,N*W) - interleaved2complex(coutPtr.Value),'fro');
    [test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    fprintf(['DGTREAL TIMEINV WP L:%3i, gl:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],L,gl,W,a,M,flags.complexity,ltfatstatusstring(status),fail);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function judge = TestIfSO3(mat)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes mat: A 3x3 matrix.
% Check if mat is close to or on the manifold SO(3).
% Example Inputs:
%
% clear; clc;
% mat = [1.0, 0.0,   0.0;
%        0.0, 0.1, -0.95;
%        0.0, 1.0,   0.1];
% judge = TestIfSO3(mat)
%
% Output:
% dudge =
%     0

judge = norm(DistanceToSO3(mat)) < 1e-3;
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function [ a, more ] = vec_colex_next2 ( dim_num, base, a, more )

%*****************************************************************************80
%
%% VEC_COLEX_NEXT2 generates vectors in colex order.
%
%  Discussion:
%
%    The vectors are produced in colexical order, starting with
%    (0,0,...,0),
%    (1,0,...,0),
%    ...
%    (BASE(1)-1,BASE(2)-1,...,BASE(DIM_NUM)-1).
%
%  Example:
%
%    DIM_NUM = 2,
%    BASE = [ 3, 3]
%
%    0   0
%    1   0
%    2   0
%    0   1
%    1   1
%    2   1
%    0   2
%    1   2
%    2   2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dennis Stanton, Dennis White,
%    Constructive Combinatorics,
%    Springer, 1986,
%    ISBN: 0387963472,
%    LC: QA164.S79.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer BASE(DIM_NUM), the base to be used in each dimension.
%
%    Input, integer A(DIM_NUM), except on the first call, this should
%    be the output value of A on the last call.
%
%    Input, logical MORE, should be FALSE on the first call, and
%    thereafter should be the output value of MORE from the previous call.
%
%    Output, integer A(DIM_NUM), the next vector.
%
%    Output, logical MORE, is TRUE if another vector was computed.
%    If MORE is FALSE on return, then ignore the output value A, and
%    stop calling the routine.
%
  if ( ~more )

    a(1:dim_num) = 0;
    more = 1;

  else

    for i = 1 : dim_num

      a(i) = a(i) + 1;

      if ( a(i) < base(i) )
        return
      end

      a(i) = 0;

    end

    more = 0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [index]=elementdof(node,nnel,ndof)
%----------------------------------------------------------
%  Purpose:
%     Compute system dofs associated with each element
%
%  Synopsis:
%     [index]=elementdof(node,nnel,ndof)
%
%  Variable Description:
%     index - system dof vector associated with element "iel"
%     iel - element number whose system dofs are to be determined
%     node - nodes associated with the element "iel"
%     nnel - number of nodes per element
%     ndof - number of dofs per node
%-----------------------------------------------------------

   k=0;
   for i=1:nnel
     start = (node(i)-1)*ndof;
       for j=1:ndof
         k=k+1;
         index(k)=start+j;
       end
   end

author: Sable

repo: mcbench-benchmarks


text:
function [optimState,vp,gp,timer] = ...
    activesample_vbmc(optimState,Ns,funwrapper,vp,vp_old,gp,stats,timer,options)
%ACTIVESAMPLE_VBMC Actively sample points iteratively based on acquisition function.

NSsearch = options.NSsearch;    % Number of points for acquisition fcn
t_func = 0;

time_active = tic;

if isempty(gp)

    % No GP yet, just use provided points or sample from plausible box
    [optimState,t_func] = ...
        initdesign_vbmc(optimState,Ns,funwrapper,t_func,options);

else                    % Active uncertainty sampling

    % One or multiple acquisition functions (most commonly only one)
    SearchAcqFcn = options.SearchAcqFcn;

    % Use "hedge" strategy to propose an acquisition function? (unused)
    if options.AcqHedge && numel(SearchAcqFcn) > 1
        % Choose acquisition function via hedge strategy
        optimState.hedge = acqhedge_vbmc('acq',optimState.hedge,[],options);
        idxAcq = optimState.hedge.chosen;
    end

    % Compute time cost (used by some acquisition functions)
    if optimState.iter > 2
        deltaNeff = max(1,stats.Neff(optimState.iter-1) - stats.Neff(optimState.iter-2));
    else
        deltaNeff = stats.Neff(1);
    end
    timer_iter = stats.timer(optimState.iter-1);
    % t_base = timer_iter.activeSampling + timer_iter.variationalFit + timer_iter.finalize;
    gpTrain_vec = [stats.timer.gpTrain];

    if options.ActiveVariationalSamples > 0 % Unused
        options_activevar = options;
        options_activevar.TolWeight = 0;
        options_activevar.NSentFine = options.NSent;
        options_activevar.ELCBOmidpoint = false;
        Ns_activevar = options.ActiveVariationalSamples;
    end

    % Perform GP (and possibly variational) update after each active sample
    ActiveSampleFullUpdate_flag = ...
        (options.ActiveSampleVPUpdate || options.ActiveSampleGPUpdate) && ...
        (((optimState.iter - options.ActiveSampleFullUpdatePastWarmup) <= optimState.LastWarmup) ...
        || (stats.rindex(end) > options.ActiveSampleFullUpdateThreshold));

    if ActiveSampleFullUpdate_flag && Ns > 1
        % Temporarily change options for local updates

        RecomputeVarPost_old = optimState.RecomputeVarPost;
        entropy_alpha_old = optimState.entropy_alpha;

        % optimState.RecomputeVarPost = false;
        options_update = options;
        % options_update.GPRetrainThreshold = Inf;
        % options_update.GPSampleThin = 1;
        options_update.GPTolOpt = options.GPTolOptActive;
        options_update.GPTolOptMCMC = options.GPTolOptMCMCActive;
        options_update.TolWeight = 0;
        options_update.NSent = options.NSentActive;
        options_update.NSentFast = options.NSentFastActive;
        options_update.NSentFine = options.NSentFineActive;
%        options_update.ELCBOWeight = -options.OptimisticVariationalBound;

%        options_update.TolFunStochastic = 3*options.TolFunStochastic;
%        options_update.DetEntTolOpt = 3*options.DetEntTolOpt;
%        options_update.NSgpMaxMain = 3;
%        options_update.NSgpMaxWarmup = 3;
        hypstruct = [];
        vp0 = vp;
    end

    % if Ns > 1; vp_old = vp; end

    %% Active sampling loop (sequentially acquire Ns new points)
    for is = 1:Ns

        optimState.N = optimState.Xn;  % Number of training inputs
        optimState.Neff = sum(optimState.nevals(optimState.X_flag));

        if options.ActiveVariationalSamples > 0 % Unused
            [vp,~,output] = vpsample_vbmc(Ns_activevar,0,vp,gp,optimState,options_activevar,options.ScaleLowerBound);
            if isfield(output,'stepsize'); optimState.mcmc_stepsize = output.stepsize; end
%            output.stepsize
%            ELCBOWeight = sqrt(0.2*2*log(vp.D*optimState.Neff^2*pi^2/(6*0.1)));
%            options_temp.ELCBOWeight = -log(rand())*ELCBOWeight;
%            vp = vpoptimize_vbmc(0,1,vp,gp,vp.K,optimState,options_temp,0);
%            vbmc_iterplot(vp,gp,optimState,stats,stats.elbo(end));
%            vbmc_iterplot(vp_old,gp,optimState,stats,stats.elbo(end));
        end

        Nextra = evaloption_vbmc(options.SampleExtraVPMeans,vp.K);
        if Nextra > 0   % Unused
            vp_base = vp;
            NsFromGP = 4e3;
            Nextra = evaloption_vbmc(options.SampleExtraVPMeans,vp.K);
            gpbeta = -options.OptimisticVariationalBound;

            xrange = max(gp.X) - min(gp.X);
            LB_extra = min(gp.X) - 0.1*xrange; UB_extra = max(gp.X) + 0.1*xrange;
            if 0
                x0 = vbmc_rnd(vp,1,0);
                xx = gplite_sample(gp,NsFromGP,x0,[],[],gpbeta,Inf,[],[],[LB_extra;UB_extra]);
            elseif 0
                K = 2*(vp.D+1);
                x0 = vbmc_rnd(vp,K,0);
                xx = gplite_sample(gp,NsFromGP,x0,'parallel',[],gpbeta,Inf,[],[],[LB_extra;UB_extra]);
            else
                [X_hpd,y_hpd] = gethpd_vbmc(gp.X,gp.y,options.HPDFrac);
                Nextra = min(size(X_hpd,1),Nextra);
                xx = gp.X(randperm(size(X_hpd,1),Nextra),:);
                xx = gp.X;
                Nextra = size(xx,1);
            end
            OptimizeMu = vp.optimize_mu;
            OptimizeWeights = vp.optimize_weights;
            vp.optimize_mu = false;
            vp.optimize_weights = true;
            vp.K = vp.K + Nextra;
            vp.mu = [vp.mu,xx(round(linspace(1,size(xx,1),Nextra)),:)'];
            vp.sigma = [vp.sigma,exp(mean(log(vp.sigma)))*ones(1,Nextra)];
            vp.w = [vp.w,exp(mean(log(vp.w)))*ones(1,Nextra)];
            vp.w = vp.w(:)'/sum(vp.w);
            options_vp = options;
            options_vp.NSent = 0;
            options_vp.NSentFast = 0;
            options_vp.NSentFine = 0;
            options_vp.TolWeight = 0;
%            options_vp.ELCBOWeight = 0;
%            vp2 = vpoptimize_vbmc(0,1,vp,gp,[],optimState,options_vp,0);
%            xrnd = vbmc_rnd(vp2,1e5); cornerplot(xrnd);
            options_vp.ELCBOWeight = -options.OptimisticVariationalBound;
            vp = vpoptimize_vbmc(0,1,vp,gp,[],optimState,options_vp,0);
%            xrnd = vbmc_rnd(vp,1e5); cornerplot(xrnd);
            vp.optimize_mu = OptimizeMu;
            vp.optimize_weights = OptimizeWeights;
        end

        if ~options.AcqHedge
            % If multiple acquisition functions are provided and not
            % following a "hedge" strategy, pick one at random
            idxAcq = randi(numel(SearchAcqFcn));
        end

        %% Pre-computations for acquisition functions

        % Re-evaluate variance of the log joint if requested
        if isfield(optimState.acqInfo{idxAcq},'compute_varlogjoint') ...
                && optimState.acqInfo{idxAcq}.compute_varlogjoint
            [~,~,varF] = gplogjoint(vp,gp,0,0,0,1);
            optimState.varlogjoint_samples = varF;
        end

        % Evaluate noise at each training point
        Ns_gp = numel(gp.post);
        sn2new = zeros(size(gp.X,1),Ns_gp);
        for s = 1:Ns_gp
            hyp_noise = gp.post(s).hyp(gp.Ncov+1:gp.Ncov+gp.Nnoise); % Get noise hyperparameters
            if isfield(optimState,'S')
                s2 = (optimState.S(optimState.X_flag).^2).*optimState.nevals(optimState.X_flag);
            else
                s2 = [];
            end
            if options.NoiseShaping
                s2 = noiseshaping_vbmc(s2,gp.y,options);
            end
            sn2new(:,s) = gplite_noisefun(hyp_noise,gp.X,gp.noisefun,gp.y,s2);
        end
        gp.sn2new = mean(sn2new,2);

        % Evaluate GP input length scale (use geometric mean)
        D = size(gp.X,2);
        ln_ell = zeros(D,Ns_gp);
        for s = 1:Ns_gp; ln_ell(:,s) = gp.post(s).hyp(1:D); end
        optimState.gplengthscale = exp(mean(ln_ell,2))';

        % Rescale GP training inputs by GP length scale
        gp.X_rescaled = bsxfun(@rdivide,gp.X,optimState.gplengthscale);

        % Algorithmic time per iteration (from last iteration)
        t_base = timer_iter.activeSampling + timer_iter.variationalFit + timer_iter.finalize + timer_iter.gpTrain;

        % Estimated increase in cost for a new training input
        if optimState.iter > 3
            len = 10;
            xx = log(stats.N(max(end-len,ceil(end/2)):end));
            yy = log(gpTrain_vec(max(end-len,ceil(end/2)):end));
            if numel(unique(xx)) > 1
                p = polyfit(xx,yy,1);
                gpTrain_diff = diff(exp(polyval(p,log([stats.N(end),stats.N(end)+1]))));
            else
                gpTrain_diff = 0;
            end
        else
            gpTrain_diff = 0;
        end

        % Algorithmic cost per function evaluation
        optimState.t_algoperfuneval = t_base/deltaNeff + max(0,gpTrain_diff);
        % [t_base/deltaNeff,gpTrain_diff]

        % Prepare for importance sampling based acquisition function
        if isfield(optimState.acqInfo{idxAcq},'importance_sampling') ...
                && optimState.acqInfo{idxAcq}.importance_sampling
            optimState.ActiveImportanceSampling = activeimportancesampling_vbmc(vp,gp,SearchAcqFcn{idxAcq},optimState.acqInfo{idxAcq},options);
        end

        %% Start active search

        optimState.acqrand = rand();    % Seed for random acquisition fcn

        % Create fast search set from cache and randomly generated
        [Xsearch,idx_cache] = getSearchPoints(NSsearch,optimState,vp,gp,options);
        Xsearch = real2int_vbmc(Xsearch,vp.trinfo,optimState.integervars);

%        acqEval = @(Xs_,vp_,gp_,optimState_,transpose_flag_) ...
%            SearchAcqFcn{idxAcq}(Xs_,vp_,gp_,optimState_,transpose_flag_);
        acqEval = @(Xs_,vp_,gp_,optimState_,transpose_flag_) ...
            acqwrapper_vbmc(Xs_,vp_,gp_,optimState_,transpose_flag_,...
            SearchAcqFcn{idxAcq},optimState.acqInfo{idxAcq});

        % Evaluate acquisition function
        acq_fast = acqEval(Xsearch,vp,gp,optimState,0);

        if options.SearchCacheFrac > 0
            [~,ord] = sort(acq_fast,'ascend');
            optimState.SearchCache = Xsearch(ord,:);
            idx = ord(1);
        else
            [~,idx] = min(acq_fast);
        end
        % idx/numel(acq_fast)
        Xacq = Xsearch(idx,:);
        idx_cache_acq = idx_cache(idx);

        % Remove selected points from search set
        Xsearch(idx,:) = []; idx_cache(idx) = [];
        % [size(Xacq,1),size(Xacq2,1)]

        % Additional search via optimization
        if ~strcmpi(options.SearchOptimizer,'none')
            fval_old = acqEval(Xacq(1,:),vp,gp,optimState,0);
            x0 = real2int_vbmc(Xacq(1,:),vp.trinfo,optimState.integervars);
            if all(isfinite([optimState.LB_search,optimState.UB_search]))
                LB = min([x0;optimState.LB_search]);
                UB = max([x0;optimState.UB_search]);
            else
                xrange = max(gp.X) - min(gp.X);
                LB = min([gp.X;x0]) - 0.1*xrange; UB = max([gp.X;x0]) + 0.1*xrange;
            end

            if isfield(optimState.acqInfo{idxAcq},'log_flag') ...
                    && optimState.acqInfo{idxAcq}.log_flag
                TolFun = 1e-2;
            else
                TolFun = max(1e-12,abs(fval_old*1e-3));
            end

            try
                switch lower(options.SearchOptimizer)
                    case 'cmaes'
                        if options.SearchCMAESVPInit
                            [~,Sigma] = vbmc_moments(vp,0);
                            %[~,idx_nearest] = min(sum(bsxfun(@minus,vp.mu,x0(:)).^2,1));
                            %Sigma = diag(vp.sigma(idx_nearest)^2.*vp.lambda.^2);
                        else
                            X_hpd = gethpd_vbmc(gp.X,gp.y,options.HPDFrac);
                            Sigma = cov(X_hpd,1);
                        end
                        insigma = sqrt(diag(Sigma));
                        cmaes_opts = options.CMAESopts;
                        cmaes_opts.TolFun = TolFun;
                        cmaes_opts.MaxFunEvals = options.SearchMaxFunEvals;
                        cmaes_opts.LBounds = LB(:);
                        cmaes_opts.UBounds = UB(:);
%                        [xsearch_optim,fval_optim,~,~,out_optim,bestever] = cmaes_modded(func2str(SearchAcqFcn{idxAcq}),x0(:),insigma,cmaes_opts,vp,gp,optimState,1);
                        [xsearch_optim,fval_optim,~,~,out_optim,bestever] = ...
                            cmaes_modded('acqwrapper_vbmc',x0(:),insigma,cmaes_opts,vp,gp,optimState,1,SearchAcqFcn{idxAcq},optimState.acqInfo{idxAcq});
                        nevals = out_optim.evals;
                        if options.SearchCMAESbest
                            xsearch_optim = bestever.x;
                            fval_optim = bestever.f;
                        end
                        % out_optim.evals
                        xsearch_optim = xsearch_optim';
                    case 'fmincon'
                        fmincon_opts.Display = 'off';
                        fmincon_opts.TolFun = TolFun;
                        fmincon_opts.MaxFunEvals = options.SearchMaxFunEvals;
                        [xsearch_optim,fval_optim,~,out_optim] = fmincon(@(x) acqEval(x,vp,gp,optimState,0),x0,[],[],[],[],LB,UB,[],fmincon_opts);
                        nevals = out_optim.funcCount;
                        % out_optim
                    case 'bads'
                        bads_opts.Display = 'off';
                        bads_opts.TolFun = TolFun;
                        bads_opts.MaxFunEvals = options.SearchMaxFunEvals;
                        [xsearch_optim,fval_optim,~,out_optim] = bads(@(x) acqEval(x,vp,gp,optimState,0),x0,LB,UB,LB,UB,[],bads_opts);
                    case 'slicesample'
                        samplefun = @(x) -acqEval(x,vp,gp,optimState,0);
                        sampleopts.Thin = 1;
                        sampleopts.Burnin = 10+D;
                        sampleopts.Display = 'off';
                        sampleopts.Diagnostics = false;
                        Widths = (UB - LB)/sqrt(12);
                        [xsearch_optim,fval_optim,~,out_optim] = ...
                            slicesamplebnd(samplefun,x0,1,Widths,LB,UB,sampleopts);
                        [x0;xsearch_optim]

                    otherwise
                        error('vbmc:UnknownOptimizer','Unknown acquisition function search optimization method.');
                end
            catch
                fprintf('Active search failed.\n');
                fval_optim = Inf;
            end

            % [fval_old, fval_optim,nevals]

            if fval_optim < fval_old
                Xacq(1,:) = real2int_vbmc(xsearch_optim,vp.trinfo,optimState.integervars);
                idx_cache_acq(1) = 0;
                % idx_cache = [idx_cache(:); 0];
                % Double check if the cache indexing is correct
            end
        end

        % [acq,tr_p] = acqEval(Xacq(1,:),vp,gp,optimState,0)

        if options.UncertaintyHandling && options.MaxRepeatedObservations > 0
            if optimState.RepeatedObservationsStreak >= options.MaxRepeatedObservations
                % Maximum number of consecutive repeated observations
                % (to prevent getting stuck in a wrong belief state)
                optimState.RepeatedObservationsStreak = 0;
            else
                % Re-evaluate acquisition function on training set
                X_train = get_traindata_vbmc(optimState,options);

                % Disable variance-based regularization first
                oldflag = optimState.VarianceRegularizedAcqFcn;
                optimState.VarianceRegularizedAcqFcn = false;
                % Use current cost of GP instead of future cost
                old_t_algoperfuneval = optimState.t_algoperfuneval;
                optimState.t_algoperfuneval = t_base/deltaNeff;
                acq_train = acqEval(X_train,vp,gp,optimState,0);
                optimState.VarianceRegularizedAcqFcn = oldflag;
                optimState.t_algoperfuneval = old_t_algoperfuneval;
                [acq_train,idx_train] = min(acq_train);

                acq_now = acqEval(Xacq(1,:),vp,gp,optimState,0);

                % [acq_train,acq_now]

                if acq_train < options.RepeatedAcqDiscount*acq_now
                    Xacq(1,:) = X_train(idx_train,:);
                    optimState.RepeatedObservationsStreak = optimState.RepeatedObservationsStreak + 1;
                else
                    optimState.RepeatedObservationsStreak = 0;
                end
            end
        end

        y_orig = [NaN; optimState.Cache.y_orig(:)]; % First position is NaN (not from cache)
        yacq = y_orig(idx_cache_acq+1);
        idx_nn = ~isnan(yacq);
        if any(idx_nn)
            yacq(idx_nn) = yacq(idx_nn) + warpvars_vbmc(Xacq(idx_nn,:),'logp',optimState.trinfo);
        end

        xnew = Xacq(1,:);
        idxnew = 1;

        % See if chosen point comes from starting cache
        idx = idx_cache_acq(idxnew);
        if idx > 0; y_orig = optimState.Cache.y_orig(idx); else; y_orig = NaN; end
        timer_func = tic;
        if isnan(y_orig)    % Function value is not available, evaluate
            try
                [ynew,optimState,idx_new] = funlogger_vbmc(funwrapper,xnew,optimState,'iter');
            catch func_error
                % pause
            end
        else
            [ynew,optimState,idx_new] = funlogger_vbmc(funwrapper,xnew,optimState,'add',y_orig);
            % Remove point from starting cache
            optimState.Cache.X_orig(idx,:) = [];
            optimState.Cache.y_orig(idx) = [];
        end
        t_func = t_func + toc(timer_func);

        % ynew = outputwarp(ynew,optimState,options);
        if isfield(optimState,'S')
            s2new = optimState.S(idx_new)^2;
        else
            s2new = [];
        end
        tnew = optimState.funevaltime(idx_new);

        if 1
            % Store acquisition information (mostly for debug)
            if ~isfield(optimState,'acqtable'); optimState.acqtable = []; end
            [~,~,fmu,fs2] = gplite_pred(gp,xnew);
            v = [idxAcq,ynew,fmu,sqrt(fs2)];
            optimState.acqtable = [optimState.acqtable; v];
        end

        if Nextra > 0
            vp = vp_base;
        end

        if is < Ns
            % If not the last sample, update GP and possibly other things
            % (no need perform updates after the last sample)

%             w_orig = optimState.UB_orig - optimState.LB_orig;
%             x_min = warpvars_vbmc(optimState.LB_orig + 1e-6*w_orig,'d',vp.trinfo);
%             x_max = warpvars_vbmc(optimState.UB_orig - 1e-6*w_orig,'d',vp.trinfo);
%             if is == 1
%                 [optimState.LB_search;x_min]
%                 [optimState.UB_search;x_max]
%             end
%             optimState.LB_search = max(optimState.LB_search,x_min);
%             optimState.UB_search = min(optimState.UB_search,x_max);

            if ActiveSampleFullUpdate_flag
                % If performing full updates with active sampling, the GP
                % hyperparameters are updated after each acquisition

                % Quick GP update
                if isempty(hypstruct); hypstruct = optimState.hypstruct; end

                fESS_thresh = options.ActiveSamplefESSThresh;
                gptmp = [];
                if fESS_thresh < 1
                    gptmp = gpreupdate(gp,optimState,options);
                    fESS = fess_vbmc(vp,gptmp,100);
                else
                    fESS = 0;
                end

                if fESS <= fESS_thresh
                    if options.ActiveSampleGPUpdate
                        t = tic;
                        [gp,hypstruct,Ns_gp,optimState] = ...
                            gptrain_vbmc(hypstruct,optimState,stats,options_update);
                        timer.gpTrain = timer.gpTrain + toc(t);
                    else
                        if isempty(gptmp)
                            gp = gpreupdate(gp,optimState,options);
                        else
                            gp = gptmp;
                        end
                    end

                    if options.ActiveSampleVPUpdate
                        % Quick variational optimization
                        t = tic;
                        % Decide number of fast optimizations
                        Nfastopts = ceil(options_update.NSelboIncr * evaloption_vbmc(options_update.NSelbo,vp.K));
                        if options.UpdateRandomAlpha
                            optimState.entropy_alpha = 1 - sqrt(rand());
                        end
                        vp = vpoptimize_vbmc(Nfastopts,1,vp,gp,[],optimState,options_update,0);
                        optimState.vp_repo{end+1} = get_vptheta(vp);
                        timer.variationalFit = timer.variationalFit + toc(t);
                    end
                else
                    gp = gptmp;
                end

            else
                % If NOT performing full updates with active sampling, only
                % the GP posterior is updated (but not the hyperparameters)

                % Perform simple rank-1 update if no noise and first sample
                t = tic;
                update1 = (isempty(s2new) || optimState.nevals(idx_new) == 1) && ~options.NoiseShaping;
                if update1
                    gp = gplite_post(gp,xnew,ynew,[],[],[],s2new,1);
                    gp.t(end+1) = tnew;
                else
                    gp = gpreupdate(gp,optimState,options);
                end
                timer.gpTrain = timer.gpTrain + toc(t);
            end
        end

        % Check if active search bounds need to be expanded
        delta_search = 0.05*(optimState.UB_search - optimState.LB_search);

        % ADD DIFFERENT CHECKS FOR INTEGER VARIABLES!
        idx = abs(xnew - optimState.LB_search) < delta_search;
        optimState.LB_search(idx) = max(optimState.LB(idx), ...
            optimState.LB_search(idx) - delta_search(idx));
        idx = abs(xnew - optimState.UB_search) < delta_search;
        optimState.UB_search(idx) = min(optimState.UB(idx), ...
            optimState.UB_search(idx) + delta_search(idx));

        % Hard lower/upper bounds on search
        prange = optimState.PUB - optimState.PLB;
        LB_searchmin = max(optimState.PLB - 2*prange*options.ActiveSearchBound,optimState.LB);
        UB_searchmin = min(optimState.PUB + 2*prange*options.ActiveSearchBound,optimState.UB);
%         optimState.LB_search = max(optimState.LB_search,LB_searchmin);
%         optimState.UB_search = min(optimState.UB_search,UB_searchmin);

        % [[optimState.LB_search, min(gp.X)]; [optimState.UB_search, max(gp.X)]]

    end

    if ActiveSampleFullUpdate_flag && Ns > 1
        % Reset OPTIMSTATE
        optimState.RecomputeVarPost = RecomputeVarPost_old;
        optimState.entropy_alpha = entropy_alpha_old;
        optimState.hypstruct = hypstruct;

        % If variational posterior has changed, check if old variational
        % posterior is better than current
        theta0 = get_vptheta(vp0);
        theta = get_vptheta(vp);
        if numel(theta0) ~= numel(theta) || any(theta0 ~= theta)
            NSentFineK = ceil(evaloption_vbmc(options.NSentFineActive,vp0.K)/vp0.K);
            elbo0 = -negelcbo_vbmc(theta0,0,vp0,gp,NSentFineK,0,1);
            if elbo0 > vp.stats.elbo(end); vp = vp0; end
            % [elbo0,vp.stats.elbo]
        end
    end
end

timer.activeSampling = timer.activeSampling + toc(time_active) ...
    - t_func - timer.gpTrain - timer.variationalFit;
timer.funEvals = timer.funEvals + t_func;

% Remove temporary fields (unnecessary here because GP is not returned)
% if isfield(gp,'sn2new'); gp = rmfield(gp,'sn2new'); end
% if isfield(gp,'X_rescaled'); gp = rmfield(gp,'X_rescaled'); end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Xsearch,idx_cache] = getSearchPoints(NSsearch,optimState,vp,gp,options)
%GETSEARCHPOINTS Get search points from starting cache or randomly generated.

% Take some points from starting cache, if not empty
x0 = optimState.Cache.X_orig;

if ~isempty(x0)
    cacheFrac = options.CacheFrac;  % Fraction of points from cache (if nonempty)
    Ncache = ceil(NSsearch*cacheFrac);
    idx_cache = randperm(size(x0,1),min(Ncache,size(x0,1)));
    Xsearch = warpvars_vbmc(x0(idx_cache,:),'d',optimState.trinfo);
else
    Xsearch = []; idx_cache = [];
end

% Randomly sample remaining points
if size(Xsearch,1) < NSsearch
    Nrnd = NSsearch-size(Xsearch,1);

    Xrnd = [];
    Nsearchcache = round(options.SearchCacheFrac*Nrnd);
    if Nsearchcache > 0 % Take points from search cache
        Xrnd = [Xrnd; optimState.SearchCache(1:min(end,Nsearchcache),:)];
    end
    Nheavy = round(options.HeavyTailSearchFrac*Nrnd);
    if Nheavy > 0
        Xrnd = [Xrnd; vbmc_rnd(vp,Nheavy,0,1,3)];
    end
    Nmvn = round(options.MVNSearchFrac*Nrnd);
    if Nmvn > 0
        [mubar,Sigmabar] = vbmc_moments(vp,0);
        Xrnd = [Xrnd; mvnrnd(mubar,Sigmabar,Nmvn)];
    end
    Nhpd = round(options.HPDSearchFrac*Nrnd);
    if Nhpd > 0
        hpd_min = options.HPDFrac/8;
        hpd_max = options.HPDFrac;
        hpdfracs = sort([rand(1,4)*(hpd_max-hpd_min) + hpd_min,hpd_min,hpd_max]);
        Nhpd_vec = diff(round(linspace(0,Nhpd,numel(hpdfracs)+1)));
        X = optimState.X(optimState.X_flag,:);
        y = optimState.y(optimState.X_flag);
        D = size(X,2);
        for ii = 1:numel(hpdfracs)
            if Nhpd_vec(ii) == 0; continue; end
            X_hpd = gethpd_vbmc(X,y,hpdfracs(ii));
            if isempty(X_hpd)
                [~,idxmax] = max(y);
                mubar = X(idxmax,:);
                Sigmabar = cov(X);
            else
                mubar = mean(X_hpd,1);
                Sigmabar = cov(X_hpd,1);
            end
            if isscalar(Sigmabar); Sigmabar = Sigmabar*ones(D,D); end
            %[~,idxmax] = max(y);
            %x0 = optimState.X(idxmax,:);
            %[Sigmabar,mubar] = covcma(X,y,x0,[],hpdfracs(ii));
            Xrnd = [Xrnd; mvnrnd(mubar,Sigmabar,Nhpd_vec(ii))];
        end
    end
    Nbox = round(options.BoxSearchFrac*Nrnd);
    if Nbox > 0
        X = optimState.X(optimState.X_flag,:);
        D = size(X,2);
        X_diam = max(X) - min(X);
        plb = optimState.PLB;
        pub = optimState.PUB;
        if all(isfinite([optimState.LB_search,optimState.UB_search]))
            box_lb = optimState.LB_search;
            box_ub = optimState.UB_search;
            box_lb = max(min(X) - 0.5*X_diam,box_lb);
            box_ub = min(max(X) + 0.5*X_diam,box_ub);
        else
            box_lb = plb - 3*(pub - plb);
            box_ub = pub + 3*(pub - plb);
            box_lb = max(min(X) - 0.5*X_diam,box_lb);
            box_ub = min(max(X) + 0.5*X_diam,box_ub);
        end
        Xrnd = [Xrnd; ...
            bsxfun(@plus,bsxfun(@times,rand(Nbox,D),box_ub-box_lb),box_lb)];
    end

    Nvp = max(0,Nrnd-Nsearchcache-Nheavy-Nmvn-Nhpd-Nbox);
    if Nvp > 0
        Xrnd = [Xrnd; vbmc_rnd(vp,Nvp,0,1)];
    end

    Xsearch = [Xsearch; Xrnd];
    idx_cache = [idx_cache(:); zeros(Nrnd,1)];
end

% Apply search bounds
Xsearch = bsxfun(@min,bsxfun(@max,Xsearch,optimState.LB_search),optimState.UB_search);

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function varargout = coeffs2( f, varargin)
%COEFFS2   componentwise Fourier--Chebyshev coefficients of a DISKFUNV.
%
%   [X, Y] = COEFFS2( F ) returns the coefficients for each component of
%   the DISKFUNV F in the Fourier--Chebyshev bases, where the columns
%   are the Chebyshev coefficients and the rows are Fourier coefficients.
%
%
%   [X, Y] = COEFFS2(F, M, N) returns bivariate coefficients with N Chebyshev
%   coefficients in the radial direction and M Fourier coefficients in
%   angular direction.

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

%empty check
if isempty(f)
    varargout = {};
    return
end

F1 = f.components{1};
F2 = f.components{2};
if isempty(varargin)
    varargout{1} = coeffs2(F1);
    varargout{2} = coeffs2(F2);
else
    varargout{1} = coeffs2(F1, varargin{:});
    varargout{2} = coeffs2(F2, varargin{:});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y= fcnNormalize0_1(a)
% normalizes matrix between 100 and 0; input, output: matrix
% higheest value = 0; lowest = 100;
% use y= fcnNormalize0_1(-a) for high=100; low=0;

b=[];
for i=1:size(a)
    b = [b a(i,:)]; % converts matrix to a vector
end

n=(b-(min(b)))/max(b-min(b)); % normalizes between 0 and 1
n = abs(n-1);                 % changes the normaliation between 1 and 0
n= n.*100;                      % range = 0 :100

y = vec2mat ( n,(size(a,2))); % converts vector to same size matrix

% *************************************************************************


text:
function [yi, ypi, yppi] = chermite(x, y, yp, xi, c)

% CHERMITE 1-D piecewise cubic Hermite spline
%    CHERMITE(X,Y,YP,XI,C) interpolates to find YI, the values of the
%    underlying function Y at the points in the array XI, using
%    piecewise cubic Hermite splines.  X and Y must be vectors
%    of length N.
%
%    C specifies how tangents are calculated when YP is not specified.
%    C can be:
%       0 : Finite difference (default)
%       1 : Catmull-Rom spline
%       2 : Monotone interpolation
%       3 : Monotone with Lam harmonic mean
%
%    [YI,YPI,YPPI] = CHERMITE() also returns the interpolated
%    quadratic derivative and linear second derivative of the
%    underlying function Y at points XI.

% Joe Henning - Fall 2011

if (nargin < 5)
   c = 0;
end

if (~isempty(yp))
   c = -1;
end

n = length(x);

if (c == 0)
   % precompute finite difference derivative
   yp = lfindiff (x, y);
elseif (c == 3)
   % precompute harmonic mean
   yp = zeros(1,n);
   d = zeros(1,n-1);

   for i = 1:n-1
      h = x(i+1) - x(i);
      dy = y(i+1) - y(i);
      d(i) = dy/h;
   end

   for i = 2:n-1
      if (d(i-1)*d(i) > 0)
         yp(i) = 2*d(i-1)*d(i)/(d(i-1) + d(i));
      else
         yp(i) = 0;
      end
   end

   if (d(1)*(2*d(1)-yp(2)) > 0)
      yp(1) = 2*d(1) - yp(2);
   else
      yp(1) = 0;
   end

   if (d(n-1)*(2*d(n-1)-yp(n-1)) > 0)
      yp(n) = 2*d(n-1) - yp(n-1);
   else
      yp(n) = 0;
   end
end

for i = 1:length(xi)
   % Find the right place in the table by means of a bisection.
   klo = 1;
   khi = n;
   while (khi-klo > 1)
      k = fix((khi+klo)/2.0);
      if (x(k) > xi(i))
         khi = k;
      else
         klo = k;
      end
   end

   h = x(khi) - x(klo);
   if (h == 0.0)
      fprintf('??? Bad x input to chermite ==> x values must be distinct\n');
      yi(i) = NaN;
      ypi(i) = NaN;
      yppi(i) = NaN;
      continue;
   end

   if (c == -1)
      a = yp(klo);
      b = yp(khi);
   elseif (c == 0)
      % Finite difference
      a = yp(klo);
      b = yp(khi);
   elseif (c == 1)
      % Catmull-Rom spline
      if (klo == 1)
         a = (y(khi) - y(klo))/h;
         b = (y(khi+1) - y(klo))/(x(khi+1) - x(klo));
      elseif (khi == n)
         a = (y(khi) - y(klo-1))/(x(khi) - x(klo-1));
         b = (y(khi) - y(klo))/h;
      else
         a = (y(khi) - y(klo-1))/(x(khi) - x(klo-1));
         b = (y(khi+1) - y(klo))/(x(khi+1) - x(klo));
      end
   else
      % Monotone interpolation
      if (c == 2)
         if (klo == 1)
            a = (y(khi) - y(klo))/h;
            b = ( (y(khi) - y(klo))/h + (y(khi+1) - y(khi))/(x(khi+1) - x(khi)) )/2;
         elseif (khi == n)
            a = ( (y(klo) - y(klo-1))/(x(klo) - x(klo-1)) + (y(khi) - y(klo))/h )/2;
            b = (y(khi) - y(klo))/h;
         else
            a = ( (y(klo) - y(klo-1))/(x(klo) - x(klo-1)) + (y(khi) - y(klo))/h )/2;
            b = ( (y(khi) - y(klo))/h + (y(khi+1) - y(khi))/(x(khi+1) - x(khi)) )/2;
         end
      else
         % Harmonic mean
         a = yp(klo);
         b = yp(khi);
      end

      if ( y(khi) == y(klo) )
         a = 0;
         b = 0;
      else
         alpha = a/( (y(khi) - y(klo))/h );
         beta = b/( (y(khi) - y(klo))/h );
         if (alpha < 0 || beta < 0)
            a = 0;
            b = 0;
         else
            if ( (alpha*alpha + beta*beta) > 9 )
               tau = 3/sqrt(alpha*alpha + beta*beta);
               a = tau*alpha*(y(khi) - y(klo))/h;
               b = tau*beta*(y(khi) - y(klo))/h;
            end
         end
      end
   end

   % Evaluate cubic Hermite polynomial
   t = (xi(i) - x(klo))/h;
   t2 = t*t;
   t3 = t2*t;
   h2 = h*h;
   h00 = 2*t3 - 3*t2 + 1;
   h10 = t3 - 2*t2 + t;
   h01 = -2*t3 + 3*t2;
   h11 = t3 - t2;

   yi(i) = h00*y(klo) + h10*h*a + h01*y(khi) + h11*h*b;

   % Differentiate to find the second-order interpolant
   h00 = 6*t2 - 6*t;
   h10 = 3*t2 - 4*t + 1;
   h01 = -h00;
   h11 = 3*t2 - 2*t;

   ypi(i) = h00*y(klo)/h + h10*a + h01*y(khi)/h + h11*b;

   % Differentiate to find the first-order interpolant
   h00 = 12*t - 6;
   h10 = 6*t - 4;
   h01 = -h00;
   h11 = 6*t - 2;

   yppi(i) = h00*y(klo)/h2 + h10*a/h + h01*y(khi)/h2 + h11*b/h;
end

% 3-pt finite difference
function [yp] = lfindiff(x, y)
% Finite Difference Formulae for Unequal Sub-Intervals Using Lagrange's Interpolation Formula
% Singh, Ashok K. and B. S. Bhadauria
% Int. Journal of Math. Analysis, Vol. 3, 2009, no. 17, 815-827

n = length(x);

if (n ~= length(y))
   fprintf('Error, len(x) != len(y), returning\n');
   yp = NaN;
   return;
end

for k = 1:n
   if (k == 1)
      h1 = x(k+1) - x(k);
      h2 = x(k+2) - x(k+1);
      yp(k) = -(2*h1+h2)/(h1*(h1+h2))*y(k) + (h1+h2)/(h1*h2)*y(k+1) - h1/(h2*(h1+h2))*y(k+2);
   elseif (k == n)
      h1 = x(k-1) - x(k-2);
      h2 = x(k) - x(k-1);
      yp(k) = h2/(h1*(h1+h2))*y(k-2) - (h1+h2)/(h1*h2)*y(k-1) + (h1+2*h2)/(h2*(h1+h2))*y(k);
   else
      h1 = x(k) - x(k-1);
      h2 = x(k+1) - x(k);
      yp(k) = -h2/(h1*(h1+h2))*y(k-1) - (h1-h2)/(h1*h2)*y(k) + h1/(h2*(h1+h2))*y(k+1);
   end
end


text:
%DEMO_OFDM  Demo of Gabor systems used for OFDM
%
%   This demo shows how to use a Gabor Riesz basis for OFDM.
%
%   We want to transmit a signal consisting of 0's and 1's through a
%   noisy communication channel. This is accomplished in the following
%   steps in the demo:
%
%     1) Convert this digital signal into complex valued coefficients by
%        QAM modulation.
%
%     2) Construct the signal to be transmitted by an inverse Gabor
%        transform of the complex coefficients
%
%     3) "Transmit" the signal by applying a spreading operator to the
%        signal and adding white noise
%
%     4) Convert the received signal into noisy coefficients by a Gabor
%        transform
%
%     5) Convert the noisy coefficients into bits by inverse QAM.
%
%   Some simplifications used to make this demo simple:
%
%     * We assume that the whole spectrum is available for transmission.
%
%     * The window and its dual have full length support. This is not
%       practical, because all data would have to be processed at once.
%       Instead, an FIR should be used, with both the window and its dual
%       having a short length.
%
%     * The window is periodic. The data at the very end interferes with
%       the data at the very beginning. A simple way to solve this is to
%       transmit zeros at the beginning and at the end, to flush the system
%       properly.
%
%   .. figure::
%
%      Received coefficients.
%
%      This figure shows the distribution in the complex plane of the
%      received coefficients. If the channel was perfect, all the points
%      should appear at the complex roots of unity (1,i,-1 and -i). This
%      demo is random, so everytime it is run it produces a new plot, and
%      the error rate may vary.
%

disp('Type "help demo_ofdm" to see a description of how this demo works.');
%% ----------- setup of signal and transmission system --------------------

% Number of channels to use
M=20;

% Time-distance between succesive transmission. This must be
% larger than M, otherwise the symbols will interfere.
a=24;

% Number of bits to transmit, must be divisable by 2*M
nbits=16000;

% Length (in samples) of transmitted signal.
L=nbits/(2*M)*a;

% We choose an orthonormal window.
g=gabtight(a,M,L);

%% ----------- Setup of communication channel ---------------------------

% Larger means more random
howrandom=.3;

% Rate of decay away from (1,1). Larger means smaller spread (faster decay).
spreaddecay=1.2;

% Noiselevel for the channel.
noiselevel=0.05;

% Define the symbol of the spreading operator
symbol=sparse(L,L);
for ii=1:3
  for jj=1:3
    symbol(ii,jj)=(1-abs(randn(1)*howrandom))*exp(-(ii+jj-1)*spreaddecay);
  end;
end;

% Make the symbol conserve real signals.
symbol=(symbol+involute(symbol))/2;

% Make sure that energy is conserved
symbol=symbol/sum(abs(symbol(:)));

%% ------------ Convert input data into analog signal -------------------

% Create a random stream of bits.
inputdata=round(rand(nbits,1));

% QAM modulate it
transmitdata=qam4(inputdata);

% Create the signal to be tranmitted
f=idgt(reshape(transmitdata,M,[]),g,a);

% --- transmission of signal - influence of the channel ----------

% Apply the underspread operator.
f=spreadop(f,symbol);

% add white noise.
noise = ((randn(size(f))-.5)+i*(randn(size(f))-.5));
f=f+noise*noiselevel/norm(noise)*norm(f);

% --- reconstruction of received signal ------------------------

% Obtain the noisy coefficients from the transmitted signal
receivedcoefficients = dgt(f,g,a,M);

% Convert the analog signal to the digital coefficients by inverse QAM
receivedbits=iqam4(receivedcoefficients(:));

%% --- visualization and print output -------------------------

% Plot the coefficients in the complex plane.
figure(1);
plot(receivedcoefficients(:),'.');
axis([-1 1 -1 1]);

% Test for errors.

disp(' ');
disp('Number of faulty bits:');
faulty=sum(abs(receivedbits-inputdata))

disp(' ');
disp('Error rate:');
faulty/nbits


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                                                                  %%%%%
%%%%    IEEE PES Power Grid Library - Optimal Power Flow - v21.07     %%%%%
%%%%          (https://github.com/power-grid-lib/pglib-opf)           %%%%%
%%%%             Benchmark Group - Small Angle Difference             %%%%%
%%%%                         29 - July - 2021                         %%%%%
%%%%                                                                  %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mpc = pglib_opf_case5_pjm__sad
mpc.version = '2';
mpc.baseMVA = 100.0;

%% area data
%	area	refbus
mpc.areas = [
	1	 4;
];

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
	2	 1	 300.0	 98.61	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
	3	 2	 300.0	 98.61	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
	4	 3	 400.0	 131.47	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
	5	 2	 0.0	 0.0	 0.0	 0.0	 1	    1.00000	    0.00000	 230.0	 1	    1.10000	    0.90000;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin
mpc.gen = [
	1	 20.0	 0.0	 30.0	 -30.0	 1.0	 100.0	 1	 40.0	 0.0;
	1	 85.0	 0.0	 127.5	 -127.5	 1.0	 100.0	 1	 170.0	 0.0;
	3	 260.0	 0.0	 390.0	 -390.0	 1.0	 100.0	 1	 520.0	 0.0;
	4	 100.0	 0.0	 150.0	 -150.0	 1.0	 100.0	 1	 200.0	 0.0;
	5	 300.0	 0.0	 450.0	 -450.0	 1.0	 100.0	 1	 600.0	 0.0;
];

%% generator cost data
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	 0.0	 0.0	 3	   0.000000	  14.000000	   0.000000;
	2	 0.0	 0.0	 3	   0.000000	  15.000000	   0.000000;
	2	 0.0	 0.0	 3	   0.000000	  30.000000	   0.000000;
	2	 0.0	 0.0	 3	   0.000000	  40.000000	   0.000000;
	2	 0.0	 0.0	 3	   0.000000	  10.000000	   0.000000;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	 2	 0.00281	 0.0281	 0.00712	 400.0	 400.0	 400.0	 0.0	 0.0	 1	 -1.33164584752	 1.33164584752;
	1	 4	 0.00304	 0.0304	 0.00658	 426.0	 426.0	 426.0	 0.0	 0.0	 1	 -1.33164584752	 1.33164584752;
	1	 5	 0.00064	 0.0064	 0.03126	 426.0	 426.0	 426.0	 0.0	 0.0	 1	 -1.33164584752	 1.33164584752;
	2	 3	 0.00108	 0.0108	 0.01852	 426.0	 426.0	 426.0	 0.0	 0.0	 1	 -1.33164584752	 1.33164584752;
	3	 4	 0.00297	 0.0297	 0.00674	 426.0	 426.0	 426.0	 0.0	 0.0	 1	 -1.33164584752	 1.33164584752;
	4	 5	 0.00297	 0.0297	 0.00674	 240.0	 240.0	 240.0	 0.0	 0.0	 1	 -1.33164584752	 1.33164584752;
];

% INFO    : === Translation Options ===
% INFO    : Phase Angle Bound:           1.33164584752 (deg.)
% INFO    :
% INFO    : === Generator Bounds Update Notes ===
% INFO    :
% INFO    : === Base KV Replacement Notes ===
% INFO    :
% INFO    : === Transformer Setting Replacement Notes ===
% INFO    :
% INFO    : === Line Capacity Monotonicity Notes ===
% INFO    :
% INFO    : === Writing Matpower Case File Notes ===

author: power-grid-lib

repo: pglib-opf


text:
function [mappedX, mapping] = npe(X, no_dims, k, eig_impl)
%NPE Perform the Neighborhood Preserving Embedding algorithm
%
%       [mappedX, mapping] = npe(X, no_dims, k)
%       [mappedX, mapping] = npe(X, no_dims, k, eig_impl)
%
% Runs the Neighborhood Preserving Embedding algorithm on dataset X to
% reduce it to dimensionality no_dims. The number of neighbors that is used
% by LPP is specified by k (default = 12).
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if size(X, 2) > size(X, 1)
        error('Number of samples should be higher than number of dimensions.');
    end
    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('k', 'var')
        k = 12;
    end
    if ~exist('eig_impl', 'var')
        eig_impl = 'Matlab';
    end

    % Get dimensionality and number of dimensions
    [n, d] = size(X);
    mapping.mean = mean(X, 1);

    % Compute pairwise distances and find nearest neighbours (vectorized implementation)
    disp('Finding nearest neighbors...');
    [distance, neighborhood] = find_nn(X, k);
    max_k = size(neighborhood, 2);
    if nargout > 1
        mapping.nbhd = distance;
    end
    X = X';
    neighborhood = neighborhood';

    % Find reconstruction weights for all points by solving the MSE problem
    % of reconstructing a point from each neighbours. A used constraint is
    % that the sum of the reconstruction weights for a point should be 1.
    disp('Compute reconstruction weights...');
    if k > d
        tol = 1e-5;
    else
        tol = 0;
    end

    % Construct reconstruction weight matrix
    W = zeros(max_k, n);
    for i=1:n
        nbhd = neighborhood(:,i);
        if ischar(k)
           nbhd = nbhd(nbhd ~= 0);
        end
        kt = numel(nbhd);
        z = X(:,nbhd) - repmat(X(:,i), 1, kt);                  % Shift point to origin
        C = z' * z;												% Compute local covariance
        C = C + eye(kt, kt) * tol * trace(C);					% Regularization of covariance (if K > D)
        wi = C \ ones(kt, 1);                                   % Solve linear system
        wi = wi / sum(wi);                                      % Make sure that sum is 1
        W(:,i) = [wi; nan(max_k - kt, 1)];
    end

    % Now that we have the reconstruction weights matrix, we define the
    % sparse cost matrix M = (I-W)'*(I-W).
    M = sparse(1:n, 1:n, ones(1, n), n, n, 4 * max_k * n);
    for i=1:n
       w = W(:,i);
       w(isnan(w)) = 0;
       j = neighborhood(:,i);
       w = w(j ~= 0);
       j = j(j ~= 0);
       M(i, j) = M(i, j) - w';
       M(j, i) = M(j, i) - w;
       M(j, j) = M(j, j) + w * w';
    end

	% For sparse datasets, we might end up with NaNs or Infs in M. We just set them to zero for now...
	M(isnan(M)) = 0;
	M(isinf(M)) = 0;

    % Compute XWX and XX and make sure these are symmetric
    X = X';
    WP = X' * M * X;
    DP = X' * X;
    DP = (DP + DP') / 2;
    WP = (WP + WP') / 2;

    % Solve generalized eigenproblem
    if size(X, 1) > 1500 && no_dims < (size(X, 1) / 10)
        if strcmp(eig_impl, 'JDQR')
            options.Disp = 0;
            options.LSolver = 'bicgstab';
            [eigvector, eigvalue] = jdqz(WP, DP, no_dims, 'SA', options);
        else
            options.disp = 0;
            options.issym = 1;
            options.isreal = 0;
            [eigvector, eigvalue] = eigs(WP, DP, no_dims, 'SA', options);
        end
    else
        [eigvector, eigvalue] = eig(WP, DP);
    end

    % Sort eigenvalues in descending order and get smallest eigenvectors
    [eigvalue, ind] = sort(diag(eigvalue), 'ascend');
    eigvector = eigvector(:,ind(1:no_dims));

    % Compute final linear basis and map data
    mappedX = X * eigvector;
    mapping.M = eigvector;

repo: cluster


text:
% function acc = LRAccuracy(GroundTruth, Predictions) compares the
% vector of predictions with the vector of ground truth values,
% and returns the accuracy (fraction of predictions that are correct).
%
% Input:
% GroundTruth    (numInstances x 1 vector)
% Predictions    (numInstances x 1 vector)
%
% Output:
% err            (scalar between 0 and 1 inclusive)
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function acc = LRAccuracy(GroundTruth, Predictions)

    GroundTruth = GroundTruth(:);
    Predictions = Predictions(:);
    assert(all(size(GroundTruth) == size(Predictions)));

    acc = mean(GroundTruth == Predictions);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
%% Filter the best QRS detections according to the *m* quality criterion
% This is an auxiliar function for the ECGtask_QRS_detections_post_process.
% Lead selection strategy based on the best (quality metric) *m*.
% This is a lead selection algorithm, so this algorithm filters the
% first ranked lead of all the algorithms outputs included in the
% payload struct. In contrast with other postprocessing algorithm, this
% doesn't create a new detection serie, it is just a filter.
%
%   all_detections = best_m_lead(struct_in, ECG_header, start_end_this_segment)
%
% Arguments:
%
%      + struct_in: the structure which results from loading the result of
%                   the QRS detection task:
%
%                   cached_filenames = ECG_w.GetCahchedFileName('QRS_detection');
%                   struct_in = load(cached_filenames{1});
%
%      +ECG_header: [struct] OPTIONAL.
%
%             Description of the ECG typically available in the
%             ECG_header. Structure with fields:
%
%               -freq: Sampling rate in Hz. (1)
%
%               -nsig: Number of ECG leads. (size(ECG,2))
%
%               -nsamp: Number of ECG samples. (size(ECG,1))
%
%      + start_end_this_segment: an array with the first and last sample
%                                indexes.
%
% Output:
%
%      + all_detections: struct with the detections  of the best lead
%                        according to the *m* quality metric.
%
%
% Example
%
%     ECG_w.ECGtaskHandle = 'QRS_detections_post_process';
%     % Mix QRS detections strategy function
%     ECG_w.ECGtaskHandle.post_proc_func = 'calculate_artificial_QRS_detections';
%     ECG_w.ECGtaskHandle.payload = load(cached_filenames{1});
%     ECG_w.ECGtaskHandle.CalculatePerformance = true;
%     ECG_w.Run;
%
% See also ECGtask_QRS_detections_post_process, mixartif, wavedetmix, qrs_detection_and_correction
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Birthdate: 12/05/2018
% Last update: 12/05/2018
% Copyright 2008-2018
%
function all_detections = best_m_lead(struct_in, ECG_header, start_end_this_segment)

    all_detections = [];

    if( nargin < 3 || isempty(start_end_this_segment) )
        start_end_this_segment = [1 ECG_header.nsamp];
    end

    if( isfield(struct_in, 'series_quality') )
        % choose the best ranked automatic detection
        best_detections_idx = max_index(struct_in.series_quality.ratios);
        best_det_name = struct_in.series_quality.AnnNames{best_detections_idx,1};
        all_detections.(best_det_name) = struct_in.(best_det_name);

    else
        disp_string_framed(2, 'Could not identify QRS detections in the input payload.');
    end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function d = simulateDataError(d)
global ERROR_ODDS_CONST
global logger

logger = InitLogger(logger);

if nargin == 0
    d = [];
    return
end
if isempty(d)
    return
end
if isempty(ERROR_ODDS_CONST)
    return
end
if ERROR_ODDS_CONST == 0
    return
end

% Set error odds in percent units
r = round(ERROR_ODDS_CONST*rand());
if r==1
    percentErrs = 7;
    n = length(d(:,1));
    idxs = floor(n * rand(1, uint32((n*percentErrs)/100)));
    k = find(idxs==0);
    idxs(k) = 1;
    d0 = d;
    d(idxs,1) = d(idxs,1) + d(idxs,1).*rand(length(idxs),1);
    logger.Write('\n**** Simulated %d errors:    d(%d,1) orig:  %0.4g,  d(%d,1) new:  %0.4g,  error size = %0.4g ****\n\n',  ...
          length(idxs), idxs(1), d0(idxs(1),1), idxs(1), d(idxs(1),1), abs(d0(idxs(1),1)-d(idxs(1),1)));
end

author: BUNPC


text:
function x = testMaps()
%tests the functionality of the functions in testMaps especially drawLine,
%   drawConc, drawFlux
%
%   testMaps
%
%
%

oriFolder = pwd;

mFilePath = mfilename('fullpath');
cd(mFilePath(1:end-length(mfilename)));

mapCoordinateFilename='ecoli_core_map.txt';

map = readCbMap(mapCoordinateFilename);
load('ecoli_core_model.mat');
sol = optimizeCbModel(model);
global CB_MAP_OUTPUT
if ~strcmp(CB_MAP_OUTPUT,'svg')
    f = figure;
end
drawCbMap(map);

%tests drawConc with zero concentrations of metabolites
disp('Testing drawConc with zero concentrations');
conc = zeros(size(model.metNames,1), 1);
drawConc(map, model, conc);
disp(' ');
disp(' ');

disp('Testing drawConc with non-zero concentrations');
conc = rand(size(model.metNames,1),1);
drawConc(map, model, conc);
disp(' ');
disp(' ');

%tests drawFlux with zero fluxes
disp('Testing drawFlux with zero flux');
flux = zeros(size(model.rxns,1), 1);
drawFlux(map, model, flux);
disp(' ');
disp(' ');

%tests drawFlux with non-zero fluxes w/ reaction with zero flux set to
%width of 1
disp('Testing drawFlux with non-zero flux');
flux = sol.x;
drawFlux(map, model, flux, [],'zeroFluxWidth',1);
disp(' ');
disp(' ');

%tests drawFluxVariablity with enlarged arrowheads to denote flux
%directionality
load('testMapsData.mat','minFlux','maxFlux');
display('Testing drawFluxVariability')
drawFluxVariability(map,model,minFlux,maxFlux,[],'rxnDirMultiplier',2.5);
display(' ');
display(' ');

x=1;
if ~strcmp(CB_MAP_OUTPUT,'svg')
    close(f);
end
cd(oriFolder);

end

author: opencobra

repo: cobratoolbox


text:
%%************************************************************
%% Prod3: compute the entries of Q = A*B*C specified in
%%        nzlistQ.
%%
%% Q = Prod3(blk,A,B,C,sym,nzlistQ)
%% Important: (a) A is assumed to be symmetric if nzlistQ
%%                has 2 columns (since mexProd2nz computes A'*B).
%%            (b) The 2nd column of nzlistQ must be sorted in
%%                ascending order.
%%
%% (optional) sym = 1, if Q is symmetric.
%%                = 0, otherwise.
%% (optional) nzlistQ = list of non-zero elements of Q to be
%%                      computed.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function Q = Prod3(blk,A,B,C,sym,nzlistQ)

  if (nargin<5); sym = 0; end;
  checkcell = [iscell(A) iscell(B) iscell(C)];
  if (nargin==6)
     checkcell(1,4) = iscell(nzlistQ);
  else
     nzlistQ = inf;
  end
%%
  if any(checkcell-1)
     if (size(blk,1) > 1)
        error('Prod3: blk and A,B,C are not compatible');
     end
     if strcmp(blk{1},'s')
        [len,len2] = size(nzlistQ);
        if (len == 0); nzlistQ = inf; len2 = 1; end;
        if (len2 == 1) & (nzlistQ == inf)
           tmp = Prod2(blk,A,B,0);
           Q = Prod2(blk,tmp,C,sym);
        else
           tmp = Prod2(blk,B,C,0);
           Q = mexProd2nz(blk,A,tmp,nzlistQ);
           if sym; Q = 0.5*(Q+Q'); end;
        end
     elseif strcmp(blk{1},'q') | strcmp(blk{1},'l') | strcmp(blk{1},'u')
        Q = A.*B.*C;
     end
  else
     error('Prod3: A,B,C,nzlistQ must all be matrices');
  end
%%************************************************************

author: yu-jiang

repo: radpbook

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


text:
function im = displaySegments(smap, segim, grayim, dodisp, varargin)
% im = displaySegments(smap, segim, grayim, dodisp, options)
% options: 'color': [r g b]
%          'width': line width

if ~exist('dodisp', 'var')
    dodisp = 1;
end

dispColor = [1 0 0];
scale = round(max(size(grayim))/200);
for k = 1:2:numel(varargin)
    if strcmp(varargin{k}, 'color')
        dispColor = varargin{k+1};
    elseif strcmp(varargin{k}, 'width')
        scale = varargin{k+1};
    end
end

 [h,s,v] = rgb2hsv(label2rgb(smap(segim)));
%
 imagesc(hsv2rgb(h,s,im2double(grayim)))

[gx, gy] = gradient(double(smap(segim)));
g = gx.^2 + gy.^2;

g = conv2(g, ones(scale), 'same');
edgepix = find(g>0);

npix = prod(size(segim));

im = repmat(grayim, [1 1 3]);
%im = 0.5*im2double(label2rgb(smap(segim))) + 0.5*im;
% size(im)
% size(segim)
% max(edgepix)
for b = 1:3
    im((b-1)*npix+edgepix) = dispColor(b);
end

if dodisp
    imagesc(im), axis image
end

author: zouchuhang

repo: LayoutNet