text:
function [K] = ku1u1(x, y, xp, yp, hyp, i)

logsigma = hyp(1);
logthetax = hyp(2);
logthetay = hyp(3);

n_x = size(x,1);
n_y = size(y,1);
n_xp = size(xp,1);
n_yp = size(yp,1);

x = repmat(x,1,n_xp);
y = repmat(y,1,n_yp);
xp = repmat(xp',n_x,1);
yp = repmat(yp',n_y,1);

switch i

case 0

K=exp(1).^(logsigma+(-2).*logthetay+(-1/2).*exp(1).^((-1).*logthetax).*(x+ ...
  (-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2);

case 1 % logsigma

K=exp(1).^(logsigma+(-2).*logthetay+(-1/2).*exp(1).^((-1).*logthetax).*(x+ ...
  (-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2);

case 2 % logthetax

K=(1/2).*exp(1).^(logsigma+(-1).*logthetax+(-2).*logthetay+(-1/2).*exp(1) ...
  .^((-1).*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).* ...
  (y+(-1).*yp).^2).*(x+(-1).*xp).^2.*(exp(1).^logthetay+(-1).*(y+(-1).*yp) ...
  .^2);

case 3 % logthetay

K=(-1/2).*exp(1).^(logsigma+(-3).*logthetay+(-1/2).*exp(1).^((-1).* ...
  logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).* ...
  yp).^2).*(2.*exp(1).^(2.*logthetay)+(-5).*exp(1).^logthetay.*(y+(-1).* ...
  yp).^2+(y+(-1).*yp).^4);

otherwise

        K = zeros(n_x, n_xp);
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


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

% Last Modified by GUIDE v2.5 22-Oct-2020 04:46:05

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

% Choose default command line output for spaceTimeOptions
handles.output = hObject;

[mov_file,trx_file,prev_st] = myparse(varargin,'mov_file','','trx_file','','prev_st',[]);
handles.mov_file = mov_file;
handles.trx_file = trx_file;

[rfn,nframes,fid,headerinfo] = get_readframe_fcn(mov_file);
handles.rfn = rfn;
handles.nframes = nframes;
ex_frame = rfn(1);
fr_sz = size(ex_frame);
handles.width = size(ex_frame,2);
handles.height = size(ex_frame,1);
handles.fid = fid;
handles.headerinfo = headerinfo;

if ~isempty(trx_file)
    trx = load_tracks(trx_file);
    handles.trx = trx;
    handles.n_trx = length(trx);
    handles.cur_fly = 1;
    handles.frame_num = randsample(trx(handles.cur_fly).nframes,1) - trx(handles.cur_fly).off;
else
    % if not trx then dummy tracks at the center of the frame.
    trx = struct('x',repmat(fr_sz(2)/2,[1,nframes]),...
                'y',repmat(fr_sz(1)/2,[1,nframes]),...
                'theta',repmat(90,[1,nframes]),...
                'firstframe',1,...
                'endframe',nframes,...
                'off',0,...
                'nframes',nframes);

    handles.trx = trx;
    handles.n_trx = 1;
    handles.cur_fly = 1;
    handles.frame_num = randsample(nframes,1);
end

if isempty(prev_st)
    prev_st = getSTParams();
end
if ~isfield(prev_st,'x_center')
  prev_st.x_center = 0;
  prev_st.y_center = 0;
  prev_st.imsz_display = prev_st.psize*max(prev_st.npatches_x,prev_st.npatches_y);
end
handles.params = prev_st;

updateGUI(handles);
updatePlots(handles);

handles.is_canceled = true;

% Update handles structure
guidata(hObject, handles);

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

function updateGUI(handles)
params = handles.params;
set(handles.edit_patch,'String',sprintf('%d',params.psize));
set(handles.edit_horz,'String',sprintf('%d',params.npatches_x));
set(handles.edit_vert,'String',sprintf('%d',params.npatches_y));
set(handles.edit_frame,'String',sprintf('%d',handles.frame_num));
set(handles.edit_trx,'String',sprintf('%d',handles.cur_fly));
set(handles.edit_x_center,'String',sprintf('%d',params.x_center));
set(handles.edit_y_center,'String',sprintf('%d',params.y_center));
method_ndx = find(strcmp(handles.params.methods,handles.params.cur_method));
set(handles.menu_methods,'String',handles.params.method_names);
set(handles.menu_methods,'Value',method_ndx);

function handles = updateimsz(handles)
params = handles.params;
grid_sz_x = params.npatches_x*params.psize;
grid_sz_y = params.npatches_y*params.psize;
x_c = params.x_center; y_c = params.y_center;
max_sz_x = max(x_c+grid_sz_x/2,-x_c+grid_sz_x/2);
max_sz_y = max(y_c+grid_sz_y/2,-y_c+grid_sz_y/2);
max_sz = max(max_sz_x,max_sz_y)*2;
params.imsz_display = max(params.imsz_display,max_sz);
handles.params = params;
set(handles.edit_imsz,'String',sprintf('%d',params.imsz_display));

function handles = updatePlots(handles)

set(handles.figure1, 'pointer', 'watch')
handles = updateimsz(handles);
trx = handles.trx;
cur_fly = handles.cur_fly;
frame_num = handles.frame_num;
rfn = handles.rfn;
trx_fn = frame_num + trx(cur_fly).off;
cur_frame = rfn(frame_num);
x = trx(cur_fly).x(trx_fn);
y = trx(cur_fly).y(trx_fn);
theta = trx(cur_fly).theta(trx_fn);
psize = handles.params.psize;
nbins = handles.params.nbins;

ncells_x = handles.params.npatches_x;
ncells_y = handles.params.npatches_y;
x_c = handles.params.x_center;
y_c = handles.params.y_center;
winwidth = handles.params.imsz_display;
winheight = handles.params.imsz_display;
hsz = handles.params.imsz_display/2;
grid_hsz_x = ncells_x*psize/2;
grid_hsz_y = ncells_y*psize/2;
handles.cur_im = CropImAroundTrx(cur_frame,x,y,theta-pi/2,...
    (winwidth-1)/2,(winheight-1)/2);
cla(handles.axes_hog);
cla(handles.axes_hof);

ftrs = genSTFeatures(handles.rfn,handles.headerinfo,frame_num,frame_num+1,...
    handles.trx,handles.params.is_stationary,handles.params.cur_method,handles.params);

% hog ftrs
scale = 6;
im1 = handles.cur_im;

axes(handles.axes_hog);
him = imshow(imresize(im1,scale));
axis image;
colormap gray;
hold on;
axis off;

colors = hsv(nbins);
colors = colors([ (end/2+1):end 1:end/2],:);

[nr,nc,~] = size(im1);
% maxv2 = max(H(:));
maxv2 = 0.03;
wd = 0.5;
hogpatch = [wd wd -wd -wd wd;-psize psize psize -psize -psize]/2;
h = [];
H = ftrs.hogftrs{cur_fly}(:,:,:,1);
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);
dt = mean(diff(bincenters));
binedges = [bincenters(1)-dt/2,(bincenters(1:end-1)+bincenters(2:end))/2,bincenters(end)+dt/2];

for xi = 1:ncells_x
  cx = (psize/2 + (xi-1)*psize - grid_hsz_x + x_c + hsz )*scale+ 1 ;
  if cx+psize/2 > nc*scale
    break;
  end
  for yi = 1:ncells_y
    cy = (psize/2 + (yi-1)*psize - grid_hsz_y + y_c + hsz)*scale+ 1 ;
    if cy+psize/2 > nr*scale
      break;
    end

    for bini = 1:nbins
      tmp = bincenters(bini);
      curpatch = [cos(tmp) -sin(tmp); sin(tmp) cos(tmp)]*hogpatch;
      xcurr = cx + curpatch(1,:)*scale;
      ycurr = cy + curpatch(2,:)*scale;
      h(yi,xi,bini) = patch(xcurr,ycurr,colors(bini,:),...
          'LineStyle','none',...
          'FaceAlpha',min(1,H(yi,xi,bini)/maxv2));
    end
  end
end

for xi = 0:ncells_x
    plot([(xi*psize+x_c + hsz- grid_hsz_x)*scale,(xi*psize+x_c+hsz-grid_hsz_x)*scale],...
      [ (y_c+hsz-grid_hsz_y)*scale,(y_c+ncells_y*psize+hsz-grid_hsz_y)*scale],...
      'Color',[0.3,0.3,0.3,0.2]);
end
for yi = 0:ncells_y
    plot([(x_c+hsz-grid_hsz_x)*scale,(x_c+ncells_x*psize+hsz-grid_hsz_x)*scale],...
      [ (yi*psize+y_c+hsz-grid_hsz_y)*scale,(yi*psize+y_c+hsz-grid_hsz_y)*scale],...
      'Color',[0.3,0.3,0.3,0.2]);
end
scatter(hsz*scale,hsz*scale,'marker','*');
scatter((x_c+hsz)*scale,(y_c+hsz)*scale,'marker','o');

% Flow (hof) ftrs
if handles.params.is_stationary && frame_num < handles.nframes,
  ny = trx(cur_fly).y(trx_fn+1);
  nx = trx(cur_fly).x(trx_fn+1);
  ntheta = trx(cur_fly).theta(trx_fn+1);
  n_fr = rfn(frame_num+1);
else
  if frame_num == handles.nframes
      n_fr = cur_frame;
  else
      n_fr = rfn(frame_num+1);
  end
  ny = y; nx = x; ntheta = theta;
end

im2 = CropImAroundTrx(n_fr,nx,ny,ntheta-pi/2,(winwidth-1)/2,(winheight-1)/2);
axes(handles.axes_hof);
him = imshowpair(imresize(im1,scale),imresize(im2,scale),'Scaling','none');
axis image;
colormap gray;
hold on;

colors = hsv(nbins);
[nr,nc,~] = size(im1);
% maxv2 = max(F(:));
maxv2 = 1;

F = ftrs.flowftrs{cur_fly}(:,:,:,1);
h = [];
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);

bincenters2 = bincenters*2;
dt = mean(diff(bincenters2));
binedges2 = [bincenters2(1)-dt/2,(bincenters2(1:end-1)+bincenters2(2:end))/2,bincenters2(end)+dt/2];

for xi = 1:ncells_x
  cx = (psize/2  + (xi-1)*psize  + x_c +hsz-grid_hsz_x)*scale+ 1;
  if cx+psize/2 > nc*scale
    break;
  end
  for yi = 1:ncells_y
    cy = (psize/2 + (yi-1)*psize + y_c+hsz-grid_hsz_y)*scale+ 1;
    if cy+psize/2 > nr*scale
      break;
    end

    for bini = 1:nbins
      tmp = linspace(binedges2(bini),binedges2(bini+1),20);
      xcurr = cx + [0,psize/2*cos(tmp),0]*scale;
      ycurr = cy + [0,psize/2*sin(tmp),0]*scale;
      h(yi,xi,bini) = patch(xcurr,ycurr,colors(bini,:),'LineStyle','none','FaceAlpha',min(1,F(yi,xi,bini)/maxv2));
    end

  end
end
for xi = 0:ncells_x
    plot([(xi*psize+x_c+hsz-grid_hsz_x)*scale,(xi*psize+x_c+hsz-grid_hsz_x)*scale],...
      [ (y_c+hsz-grid_hsz_y)*scale,(y_c+ncells_y*psize+hsz-grid_hsz_y)*scale],...
      'Color',[0.3,0.3,0.3,0.2]);
end
for yi = 0:ncells_y
    plot([(x_c+hsz-grid_hsz_x)*scale,(x_c+ncells_x*psize+hsz-grid_hsz_x)*scale],...
      [ (yi*psize+y_c+hsz-grid_hsz_y)*scale,(yi*psize+y_c+hsz-grid_hsz_y)*scale],...
      'Color',[0.3,0.3,0.3,0.2]);
end
scatter(hsz*scale,hsz*scale,'marker','*');
scatter((x_c+hsz)*scale,(y_c+hsz)*scale,'marker','o');

set(handles.figure1, 'pointer', 'arrow')

% --- Outputs from this function are returned to the command line.
function varargout = spaceTimeOptions_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.is_canceled;
varargout{2} = handles.params;
figure1_CloseRequestFcn(hObject,eventdata,handles);

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

patch_sz = str2double(get(hObject,'String'));
prev_patch_sz = handles.params.psize;
if isempty(patch_sz) || isnan(patch_sz) || round(patch_sz)~=patch_sz || patch_sz < 1
    warndlg('Cell size must be a positive integer');
    set(handles.edit_patch,'String',sprintf('%d',prev_patch_sz));
    return;
end
handles.params.psize = patch_sz;
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_patch_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_patch (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 edit_horz_Callback(hObject, eventdata, handles)
% hObject    handle to edit_horz (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 edit_horz as text
%        str2double(get(hObject,'String')) returns contents of edit_horz as a double
ncells = str2double(get(hObject,'String'));
prev_cells = handles.params.npatches_x;
if isempty(ncells) || isnan(ncells) || round(ncells)~=ncells || ncells < 1
    warndlg('Cell size must be a positive integer');
    set(handles.edit_patch,'String',sprintf('%d',prev_cells));
    return;
end
handles.params.npatches_x = ncells;
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_horz_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_horz (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 edit_frame_Callback(hObject, eventdata, handles)
% hObject    handle to edit_frame (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 edit_frame as text
%        str2double(get(hObject,'String')) returns contents of edit_frame as a double

frame_num = str2double(get(hObject,'String'));
if isempty(frame_num) || isnan(frame_num) || round(frame_num)~=frame_num || frame_num < 1
    warndlg('Frame number must be a positive integer');
    set(handles.edit_frame,'String',sprintf('%d',handles.frame_num));
    return;
end
min_frame = handles.trx(handles.cur_fly).firstframe;
max_frame = handles.trx(handles.cur_fly).endframe;
if frame_num > max_frame || frame_num < min_frame
    warndlg(sprintf('Frame number must be between %d and %d for the current animal',min_frame,max_frame));
    set(handles.edit_frame,'String',sprintf('%d',handles.frame_num));
    return;
end
handles.frame_num = frame_num;
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_frame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_frame (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 edit_trx_Callback(hObject, eventdata, handles)
% hObject    handle to edit_trx (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 edit_trx as text
%        str2double(get(hObject,'String')) returns contents of edit_trx as a double

fly_num = str2double(get(hObject,'String'));
if isempty(fly_num) || isnan(fly_num) || round(fly_num)~=fly_num || fly_num < 1
    warndlg('Animal number must be a positive integer');
    set(handles.edit_trx,'String',sprintf('%d',handles.cur_fly));
    return;
end
if fly_num > handles.n_trx
    warndlg(sprintf('Animal number must be less than %d',handles.n_trx));
    set(handles.edit_trx,'String',sprintf('%d',handles.cur_fly));
    return;
end
handles.cur_fly = fly_num;
if handles.frame_num < handles.trx(fly_num).firstframe
    handles.frame_num = handles.trx(fly_num).firstframe;
end
if handles.frame_num > (handles.trx(fly_num).endframe-1)
    handles.frame_num = handles.trx(fly_num).endframe-1;
end
updateGUI(handles);
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_trx_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_trx (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 pushbutton_cancel.
function pushbutton_cancel_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.is_canceled=true;
guidata(hObject,handles);
uiresume(handles.figure1);

% --- Executes on button press in pushbutton_done.
function pushbutton_done_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_done (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.is_canceled=false;
guidata(hObject,handles);
uiresume(handles.figure1);

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

% Hint: get(hObject,'Value') returns toggle state of toggle_align
if get(hObject,'Value')
    handles.params.is_stationary= true;
else
    handles.params.is_stationary= false;
end
handles = updatePlots(handles);
guidata(hObject,handles);

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

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (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
try
    close(handles.fid);
catch ME

end
delete(hObject);

% --- Executes during object deletion, before destroying properties.
function figure1_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function edit_vert_Callback(hObject, eventdata, handles)
% hObject    handle to edit_vert (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 edit_vert as text
%        str2double(get(hObject,'String')) returns contents of edit_vert as a double
ncells = str2double(get(hObject,'String'));
prev_cells = handles.params.npatches_y;
if isempty(ncells) || isnan(ncells) || round(ncells)~=ncells || ncells < 1
    warndlg('Cell size must be a positive integer');
    set(handles.edit_patch,'String',sprintf('%d',prev_cells));
    return;
end
handles.params.npatches_y = ncells;
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_vert_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_vert (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 selection change in menu_methods.
function menu_methods_Callback(hObject, eventdata, handles)
% hObject    handle to menu_methods (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns menu_methods contents as cell array
%        contents{get(hObject,'Value')} returns selected item from menu_methods
method_ndx = get(hObject,'Value');
handles.params.cur_method = handles.params.methods{method_ndx};
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function menu_methods_CreateFcn(hObject, eventdata, handles)
% hObject    handle to menu_methods (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 edit_x_center_Callback(hObject, eventdata, handles)
% hObject    handle to edit_x_center (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 edit_x_center as text
%        str2double(get(hObject,'String')) returns contents of edit_x_center as a double
ncells = str2double(get(hObject,'String'));
prev_cells = handles.params.x_center;
if isempty(ncells) || isnan(ncells) || round(ncells)~=ncells
    warndlg('Cell size must be a integer');
    set(handles.edit_x_center,'String',sprintf('%d',prev_cells));
    return;
end
handles.params.x_center = ncells;
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_x_center_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_x_center (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 edit_y_center_Callback(hObject, eventdata, handles)
% hObject    handle to edit_y_center (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 edit_y_center as text
%        str2double(get(hObject,'String')) returns contents of edit_y_center as a double
ncells = str2double(get(hObject,'String'));
prev_cells = handles.params.y_center;
if isempty(ncells) || isnan(ncells) || round(ncells)~=ncells
    warndlg('Cell size must be a integer');
    set(handles.edit_y_center,'String',sprintf('%d',prev_cells));
    return;
end
handles.params.y_center = ncells;
handles = updatePlots(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_y_center_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_y_center (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 edit_imsz_Callback(hObject, eventdata, handles)
% hObject    handle to edit_imsz (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 edit_imsz as text
%        str2double(get(hObject,'String')) returns contents of edit_imsz as a double
ncells = str2double(get(hObject,'String'));
prev_cells = handles.params.imsz_display;
if isempty(ncells) || isnan(ncells) || round(ncells)~=ncells || ncells < 1
    warndlg('Cell size must be a positive integer');
    set(handles.edit_imsz,'String',sprintf('%d',prev_cells));
    return;
end
params = handles.params;
grid_sz_x = params.npatches_x*params.psize;
grid_sz_y = params.npatches_y*params.psize;
x_c = params.x_center; y_c = params.y_center;
max_sz_x = max(x_c+grid_sz_x/2,-x_c+grid_sz_x/2);
max_sz_y = max(y_c+grid_sz_y/2,-y_c+grid_sz_y/2);
max_sz = max(max_sz_x,max_sz_y)*2;
prev_sz = params.imsz_display;
if ncells < max_sz
  warndlg('Input size cannot include the whole grid. Select a size large enough to include the whole grid');
  set(handles.edit_imsz,'String',sprintf('%d',prev_sz));
  return;
end

handles.params.imsz_display = ncells;
handles = updatePlots(handles);
guidata(hObject,handles);

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

author: kristinbranson

repo: JAABA


function vl_simplenn_diagnose(net, res)
% VL_SIMPLENN_DIAGNOSE  Plot diagnostic information
%   VL_SIMPLENN_DIAGNOSE(NET, RES) plots in the current window
%   the average, maximum, and miminum element for all the filters
%   and biases in the network NET. If RES is also provided, it will
%   plot the average, minimum, and maximum element for all the
%   intermediate responses and deriviatives stored in RES as well.
%
%   This function can be used to rapidly glance at the evolution
%   of the paramters during training.

n = numel(net.layers) ;
fmu = NaN + zeros(1, n) ;
fmi = fmu ;
fmx = fmu ;
bmu = fmu ;
bmi = fmu ;
bmx = fmu ;
xmu = fmu ;
xmi = fmi ;
xmx = fmx ;
dxmu = fmu ;
dxmi = fmi ;
dxmx = fmx ;
dfmu = fmu ;
dfmi = fmu ;
dfmx = fmu ;
dbmu = fmu ;
dbmi = fmu ;
dbmx = fmu ;

for i=1:numel(net.layers)
  ly = net.layers{i} ;
  if strcmp(ly.type, 'conv') && numel(ly.filters) > 0
    x = gather(ly.filters) ;
    fmu(i) = mean(x(:)) ;
    fmi(i) = min(x(:)) ;
    fmx(i) = max(x(:)) ;
  end
  if strcmp(ly.type, 'conv') && numel(ly.biases) > 0
    x = gather(ly.biases) ;
    bmu(i) = mean(x(:)) ;
    bmi(i) = min(x(:)) ;
    bmx(i) = max(x(:)) ;
  end
  if nargin > 1
    if numel(res(i).x) > 1
      x = gather(res(i).x) ;
      xmu(i) = mean(x(:)) ;
      xmi(i) = min(x(:)) ;
      xmx(i) = max(x(:)) ;
    end
    if numel(res(i).dzdx) > 1
      x = gather(res(i).dzdx);
      dxmu(i) = mean(x(:)) ;
      dxmi(i) = min(x(:)) ;
      dxmx(i) = max(x(:)) ;
    end
    if strcmp(ly.type, 'conv') && numel(res(i).dzdw{1}) > 0
      x = gather(res(i).dzdw{1}) ;
      dfmu(i) = mean(x(:)) ;
      dfmi(i) = min(x(:)) ;
      dfmx(i) = max(x(:)) ;
    end
    if strcmp(ly.type, 'conv') && numel(res(i).dzdw{2}) > 0
      x = gather(res(i).dzdw{2}) ;
      dbmu(i) = mean(x(:)) ;
      dbmi(i) = min(x(:)) ;
      dbmx(i) = max(x(:)) ;
    end
  end
end

if nargin > 1
  np = 6 ;
else
  np = 2 ;
end

clf ; subplot(np,1,1) ;
errorbar(1:n, fmu, fmi, fmx, 'bo') ;
grid on ;
xlabel('layer') ;
ylabel('filters') ;
title('coefficient ranges') ;

subplot(np,1,2) ;
errorbar(1:n, bmu, bmi, bmx, 'bo') ;
grid on ;
xlabel('layer') ;
ylabel('biases') ;

if nargin > 1
  subplot(np,1,3) ;
  errorbar(1:n, xmu, xmi, xmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('x') ;

  subplot(np,1,4) ;
  errorbar(1:n, dxmu, dxmi, dxmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('dzdx') ;

  subplot(np,1,5) ;
  errorbar(1:n, dfmu, dfmi, dfmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('dfilters') ;

  subplot(np,1,6) ;
  errorbar(1:n, dbmu, dbmi, dbmx, 'bo') ;
  grid on ;
  xlabel('layer') ;
  ylabel('dbiases') ;
end

drawnow ;


text:
%MIRT3D_MEXINTERP  Fast 3D linear interpolation
%
% Output_image = mirt3D_mexinterp(Input_image, XI,YI,ZI) interpolates the 3D image 'Input_image' at
%    the points with coordinates X,Y,Z. Input_image is assumed to  be defined at a regular grid 1:N, 1:M, 1:K,
%    where [M,N,K]=size(Input_images). Points outside the boudary return NaNs.
%    This is equivalent (but much faster) to Matlab's:
%    Output_image = interp3(Input_image,XI,YI,ZI,'linear',NaN);
%
% Output_images = mirt3D_mexinterp(Input_images, XI,YI,ZI). Input_images can be a stack of many 3D images (4D).
%   The function interpolates each of the 3D images at X,Y,Z coordinates and return a stack of corresponding
%   interpolated images. This is equivalent to Matlab's
%
%   Input_images(:,:,:,1)=Input_image1;
%   Input_images(:,:,:,2)=Input_image2;
%   Input_images(:,:,:,3)=Input_image3;
%   Input_images(:,:,:,4)=Input_image4;
%
%   Output_images(:,:,:,1) = interp3(Input_image1,XI,YI,ZI,'linear',NaN);
%   Output_images(:,:,:,2) = interp3(Input_image2,XI,YI,ZI,'linear',NaN);
%   Output_images(:,:,:,3) = interp3(Input_image3,XI,YI,ZI,'linear',NaN);
%   Output_images(:,:,:,4) = interp3(Input_image4,XI,YI,ZI,'linear',NaN);
%
%  This is especially usefull fpr vector valued 3D images, RGB images, to interpolate the whole 3D video at the same coordinates
%  or to interpolate image and its gradients at the same time (in image registration).
%  The speed gain is also from the precomputation of nearest points for interpolation, which are the same for all images in a stack.
%
%  Andriy Myronenko, Feb 2008, email: myron@csee.ogi.edu,
%  homepage: http://www.bme.ogi.edu/~myron/

% The function below compiles the mirt3D_mexinterp.cpp file if you haven't done it yet.
% It will be executed only once at the very first run.
function Output_images = mirt3D_mexinterp(Input_images, XI,YI,ZI)

pathtofile=which('mirt3D_mexinterp.cpp');
pathstr = fileparts(pathtofile);
mex(pathtofile,'-outdir',pathstr);

Output_images = mirt3D_mexinterp(Input_images, XI,YI,ZI);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function L = evaluate_log_posterior(this, uv)
%EVALUATE_LOG_POSTERIOR computes the log-posterior (negative energy) of the
%   flow fields UV
%   Actually only proportional to the log posterior since the variance of neither the
%   spatial nor the data terms is considered
%
%   This is a member function of the class 'hs_optical_flow'.
%
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%   $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

% spatial term

S = {[1 -1], [1; -1]};
p = 0;
for i = 1:length(S)
    u_  = conv2(uv(:,:,1), S{i}, 'valid');
    v_  = conv2(uv(:,:,2), S{i}, 'valid');
    p   = p - sum(u_(:).^2) - sum(v_(:).^2);
end;

% data term
It  = partial_deriv(this.images, uv, this.interpolation_method, this.deriv_filter);
l   = -sum(It(:).^2);

L = this.lambda/this.sigmaS2*p + l/this.sigmaD2;

if this.display
    fprintf('spatial\t%3.2e\tdata\t%3.2e\n', this.lambda*p/this.sigmaS2, l/this.sigmaD2);
end;


text:
function cavity_flow_movie ( )

%*****************************************************************************80
%
%% CAVITY_FLOW_MOVIE creates an animation from flow data.
%
%  Discussion:
%
%    This MATLAB script file reads the Cavity flow data:
%
%      geometry (XY values at nodes, assumed to be in 'xy.txt')
%      flow (UV values at nodes, in a sequence of files starting with 'up001.txt')
%
%    and plots the velocity vectors (U,V)(X,Y), saving each plot and
%    making an animation.
%
%    The file plots either the velocity vector field, or the velocity
%    direction field, depending on the value of the internal logical
%    parameter "normalized".
%
%    The MATLAB routine quiver internally scales the vectors, but this
%    can be adjusted by using a value of SCALE that is not 1.
%
%    For the unnormalized case, the velocity field is computed to a fixed
%    scale determined by finding the largest velocity vector over time.
%
%    This routine requires some auxiliary routines in order to "increment"
%    the name of the current velocity data file to get the name of the next
%    one, with the main such routine being called FILE_NAME_INC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2003
%
%  Author:
%
%    John Burkardt
%
  FALSE = 0;
  TRUE = 1;
  scale = 1.0;
  normalized = TRUE;

  nframes = 500;
%
%  Get the XY coordinates of the nodes.
%
  load xy.txt;

  x = xy(:,1);
  y = xy(:,2);
%
%  Thin the data.
%
  thin_factor = 2;
  thin_dex = thin_index ( x, y, thin_factor );
  thin_num = length ( thin_dex );
  x = x(thin_dex);
  y = y(thin_dex);
%
%  For unnormalized plots, you need to make sure that a fixed scale is preserved from
%  step to step.  The only way I can see to do this requires that I determine the
%  maximum velocity over all time, and then append one extract node and velocity
%  to the data structure.
%
  if ( normalized == FALSE )

    upfile = 'up000.txt';
    vnorm_max = 0.0;

    for ( i = 1 : nframes )
      i
      upfile = file_name_inc ( upfile );
      uv = load ( upfile );
      u = uv(thin_dex,1);
      v = uv(thin_dex,2);
      norm = sqrt ( u.^2 + v.^2 );
      vnorm_max = max ( vnorm_max, max ( norm ) );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Maximum visible velocity magnitude is %f\n', vnorm_max );

  end
%
%  Set the coordinates of the boundary lines, and an extra set of
%  lines that will be invisible, but which make some space for the
%  file name to be displayed within the plot area.
%
  bx1 = [ 0.00, 0.00, 1.00, 1.00, 0.99, 0.99, 0.01, 0.01, 0.00 ];
  by1 = [ 1.00, 0.00, 0.00, 1.00, 1.00, 0.01, 0.01, 1.00, 1.00 ];

  bx2 = [ -0.10,  1.10, 1.10, -0.10, -0.10 ];
  by2 = [ -0.10, -0.10, 1.10,  1.10, -0.10 ];
%
%  Set the name of the 0th (nonexistent) velocity file.
%  All the velocity files with have this format, with the
%  numeric part of the name incremented to get the next one.
%  In particular, the first velocity file is called UP001.TXT.
%
  upfile = 'up000.txt';

  for ( i = 1 : nframes )
    i
    upfile = file_name_inc ( upfile );
    uv = load ( upfile );
    u = uv(thin_dex,1);
    v = uv(thin_dex,2);

    if ( normalized == TRUE )
      norm = sqrt ( u.^2 + v.^2 );
      nonzero = find ( norm ~= 0.0 );
      u(nonzero) = u(nonzero) ./ norm(nonzero);
      v(nonzero) = v(nonzero) ./ norm(nonzero);
      vnorm_max = 1.0;
    end

    x(thin_num+1) = 0.00;
    y(thin_num+1) = 1.05;
    u(thin_num+1) = vnorm_max;
    v(thin_num+1) = 0.0;

    quiver ( x, y, u, v, scale )
%
%  Draw the boundary lines, and the invisible bounding line.
%
    line ( bx1, by1, 'color', 'r' )
    line ( bx2, by2, 'color', 'w' )
    axis equal
    if ( normalized == TRUE )
      title ( 'Cavity Direction Field' )
    else
      title ( 'Cavity Flow Field' )
    end
    text ( 0.420, 1.03, upfile )

    my_frames(:,i) = getframe;

  end
%
%  Display the movie a given number of times.
%
  movie ( my_frames, 2 )

  return
end
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
%% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the string to be incremented.
%
%    Output, string FILE_NAME, the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end
function thin_dex = thin_index ( x, y, thin_factor )

%*****************************************************************************80
%
%%  THIN_INDEX determines thinning indices for a X, Y data.
%
%  Discussion:
%
%    A set of X, Y data is given, that is presumably, not too far off
%    from being on a rectangular grid.
%
%    The input value of THIN_FACTOR indicates by how much the data should
%    be thinned.
%
%    The X and Y ranges are computed, and only those data points are
%    retained for which both X and Y lie in an appropriate subrange.
%
%    For instance, a THIN_FACTOR of 2 would essentially save data
%    that lay in the black squares of a checkerboard.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(NODE_NUM), Y(NODE_NUM), the X and Y coordinates
%    of the nodes.
%
%    Input, integer THIN_FACTOR, the thinning factor.
%
%    Output, integer THIN_DEX(NODE_NUM), contains in (1:THIN_NUM) the
%    indices into X and Y of the vectors to be retained after thinning.
%
  TRUE = 1;
  FALSE = 0;

  node_num = length ( x );

  x_unique_num = 0;

  for ( i = 1 : node_num )

    unique = TRUE;

    for ( j = 1 : x_unique_num )
      if ( x(i) == x_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      x_unique_num = x_unique_num + 1;
      x_unique(x_unique_num) = x(i);
    end

  end

  sort ( x_unique );

  y_unique_num = 0;

  for ( i = 1 : node_num )

    unique = TRUE;

    for ( j = 1 : y_unique_num )
      if ( y(i) == y_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      y_unique_num = y_unique_num + 1;
      y_unique(y_unique_num) = y(i);
    end

  end

  sort ( y_unique );

  thin_num = 0;

  for ( i = 1 : node_num )

    for ( j = 1 : x_unique_num-1 )
      if ( x_unique(j) <= x(i) & x(i) <= x_unique(j+1) )
        x_bin = j;
        break;
      end
    end

    for ( j = 1 : y_unique_num-1 )
      if ( y_unique(j) <= y(i) & y(i) <= y_unique(j+1) )
        y_bin = j;
        break;
      end
    end

    if ( mod ( y_bin, thin_factor ) == thin_factor / 2 & ...
         mod ( x_bin, thin_factor ) == thin_factor / 2 )

      thin_num = thin_num + 1;
      thin_dex(thin_num) = i;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function T = imageStats5(f)
%imageStats5 Sample function used in Chapter 2.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt
P = size(f,3);
% Initialize table variables to hold information for P images.
ArraySize = zeros(P,2);
GlobalMean = zeros(P,1);
RowMeans = zeros(P,size(f,1));
ColumnMeans = zeros(P,size(f,2));
for k = 1:P
   fk = f(:,:,k);
   ArraySize(k,:) = size(fk);
   GlobalMean(k) = mean2(fk);
   RowMeans(k,:) = mean(fk,2)'; % Transpose to store row means as a row
                                % vector
   ColumnMeans(k,:) = mean(fk,1);
end
% Create the table from the individual variables.
T = table(ArraySize,GlobalMean,RowMeans,ColumnMeans);

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function F = real(F)
%REAL   Complex real part of a CHEBFUN.
%   REAL(F) is the real part of F.
%
% See also IMAG, ISREAL.

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

% Handle the empty case:
if ( isempty(F) )
    return
end

for j = 1:numel(F)
    % Take real part of the pointValues:
    F(j).pointValues = real(F(j).pointValues);

    % Take real part of the FUNs:
    for k = 1:numel(F(j).funs)
        F(j).funs{k} = real(F(j).funs{k});
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [est_C_b_e_new,est_v_eb_e_new,est_r_eb_e_new,est_IMU_bias_new,...
            P_matrix_new, corrections, Phi_matrix, Q_matrix] = ...
         LC_KF_Epoch(GNSS_epoch, GNSS_r_eb_e,GNSS_v_eb_e,tor_s, est_C_b_e_old,...
                     est_v_eb_e_old,est_r_eb_e_old,est_IMU_bias_old,...
                     P_matrix_old,meas_f_ib_b,est_L_b_old,LC_KF_config,...
                     Q_matrix, R_matrix, meas_omega_ib_b)
%LC_KF_Epoch - Implements one cycle of the loosely coupled INS/GNSS
% Kalman filter plus closed-loop correction of all inertial states
%
% Software for use with "Principles of GNSS, Inertial, and Multisensor
% Integrated Navigation Systems," Second Edition.
%
% This function created 12/4/2012 by Paul Groves
%
% Inputs:
%   GNSS_r_eb_e           GNSS estimated ECEF user position (m)
%   GNSS_v_eb_e           GNSS estimated ECEF user velocity (m/s)
%   tor_s                 propagation interval (s)
%   est_C_b_e_old         prior estimated body to ECEF coordinate
%                         transformation matrix
%   est_v_eb_e_old        prior estimated ECEF user velocity (m/s)
%   est_r_eb_e_old        prior estimated ECEF user position (m)
%   est_IMU_bias_old      prior estimated IMU biases (body axes)
%   P_matrix_old          previous Kalman filter error covariance matrix
%   meas_f_ib_b           measured specific force
%   est_L_b_old           previous latitude solution
%   LC_KF_config
%     .gyro_noise_PSD     Gyro noise PSD (rad^2/s)
%     .accel_noise_PSD    Accelerometer noise PSD (m^2 s^-3)
%     .accel_bias_PSD     Accelerometer bias random walk PSD (m^2 s^-5)
%     .gyro_bias_PSD      Gyro bias random walk PSD (rad^2 s^-3)
%     .pos_meas_SD            Position measurement noise SD per axis (m)
%     .vel_meas_SD            Velocity measurement noise SD per axis (m/s)
%
% Outputs:
%   est_C_b_e_new     updated estimated body to ECEF coordinate
%                      transformation matrix
%   est_v_eb_e_new    updated estimated ECEF user velocity (m/s)
%   est_r_eb_e_new    updated estimated ECEF user position (m)
%   est_IMU_bias_new  updated estimated IMU biases
%     Rows 1-3          estimated accelerometer biases (m/s^2)
%     Rows 4-6          estimated gyro biases (rad/s)
%   P_matrix_new      updated Kalman filter error covariance matrix

% Copyright 2012, Paul Groves
% License: BSD; see GrovesCode/license.txt for details
%
% Edits made by Adam Werries where noted as (Werries), Copyright 2015
% GNSS_epoch
% GNSS_r_eb_e
% GNSS_v_eb_e
% tor_s
% est_C_b_e_old
% est_v_eb_e_old
% est_r_eb_e_old
% est_IMU_bias_old
% P_matrix_old
% meas_f_ib_b
% est_L_b_old
% Q_matrix
% R_matrix
% meas_omega_ib_b
% Constants (sone of these could be changed to inputs at a later date)
c = 299792458; % Speed of light in m/s
omega_ie = 7.292115E-5;  % Earth rotation rate in rad/s
R_0 = 6378137; %WGS84 Equatorial radius in meters
e = 0.0818191908425; %WGS84 eccentricity
Lab = LC_KF_config.lever_arm;
% Begins

% Skew symmetric matrix of Earth rate
Omega_ie = Skew_symmetric([0,0,omega_ie]);

% SYSTEM PROPAGATION PHASE

% 1. Determine transition matrix using (14.50) (first-order approx)
F_21 = -Skew_symmetric(est_C_b_e_old * meas_f_ib_b);
Phi_matrix = eye(15);
Phi_matrix(1:3,1:3) = Phi_matrix(1:3,1:3) - Omega_ie * tor_s;
Phi_matrix(1:3,13:15) = est_C_b_e_old * tor_s;
Phi_matrix(4:6,1:3) = tor_s * F_21;
Phi_matrix(4:6,4:6) = Phi_matrix(4:6,4:6) - 2 * Omega_ie * tor_s;
geocentric_radius = R_0 / sqrt(1 - (e * sin(est_L_b_old))^2) *...
    sqrt(cos(est_L_b_old)^2 + (1 - e^2)^2 * sin(est_L_b_old)^2); % from (2.137)
Phi_matrix(4:6,7:9) = -tor_s * 2 * Gravity_ECEF(est_r_eb_e_old) /...
    geocentric_radius * est_r_eb_e_old' / sqrt (est_r_eb_e_old' *...
    est_r_eb_e_old);
Phi_matrix(4:6,10:12) = est_C_b_e_old * tor_s;
Phi_matrix(7:9,4:6) = eye(3) * tor_s;

% If we were supposed to adapt Q_matrix last iteration, don't replace it.
if GNSS_epoch <= LC_KF_config.n + 1
    % 2. Determine approximate system noise covariance matrix using (14.82)
    Srg = LC_KF_config.gyro_noise_PSD;
    Sra = LC_KF_config.accel_noise_PSD;
    Sbad = LC_KF_config.accel_bias_PSD;
    Sbgd = LC_KF_config.gyro_bias_PSD;
    Q11 = eye(3) * (Srg * tor_s + Sbgd*tor_s^3/3);
    Q21 = (Srg*tor_s^2 / 2 + Sbgd*tor_s^4 / 4) * F_21;
    Q31 = (Srg*tor_s^3 / 3 + Sbgd*tor_s^5 / 5) * F_21;
    Q15 = Sbgd*tor_s^2*est_C_b_e_old / 2;
    Q22 = (Sra*tor_s + Sbad*tor_s^3/3)*eye(3) + (Srg*tor_s^3/3 + Sbgd*tor_s^5/5)*(F_21*F_21');
    Q32 = (Sra*tor_s^2/2 + Sbad*tor_s^4/4)*eye(3) + (Srg*tor_s^4/4 + Sbgd*tor_s^6/6)*(F_21*F_21');
    Q24 = Sbad*tor_s^2*est_C_b_e_old / 2;
    Q25 = Sbgd*tor_s^3*F_21*est_C_b_e_old / 3;
    Q33 = (Sra*tor_s^3/3 + Sbad*tor_s^5/5)*eye(3) + (Srg*tor_s^5/5 + Sbgd*tor_s^7/7)*(F_21*F_21');
    Q34 = Sbad*tor_s^3*est_C_b_e_old/3;
    Q35 = Sbgd*tor_s^4*F_21*est_C_b_e_old/4;
    Q44 = Sbad*tor_s*eye(3);
    Q55 = Sbgd*tor_s*eye(3);
    Q_matrix = [Q11      Q21' Q31' zeros(3) Q15;
                Q21      Q22  Q32' Q24      Q25;
                Q31      Q32  Q33  Q34      Q35;
                zeros(3) Q24' Q34' Q44      zeros(3);
                Q15'     Q25' Q35' zeros(3) Q55;];
else
    Q_matrix = tor_s*Q_matrix;
end
[~,p] = chol(Q_matrix);
if p ~= 0
    Q_matrix = nearestSPD(Q_matrix);
end
% 3. Propagate state estimates using (3.14) noting that all states are zero
% due to closed-loop correction.
x_est_propagated(1:15,1) = 0;

% 4. Propagate state estimation error covariance matrix using (3.46)
P_matrix_propagated = Phi_matrix * P_matrix_old * Phi_matrix' + Q_matrix;

% MEASUREMENT UPDATE PHASE
% 5. Set-up measurement matrix using (14.115)
H_matrix = zeros(6,15);
H_matrix(1:3,1:3) = Skew_symmetric(est_C_b_e_old*Lab);
H_matrix(4:6,1:3) = Skew_symmetric(est_C_b_e_old*cross(meas_omega_ib_b,Lab) ...
                                   - Omega_ie*est_C_b_e_old*Lab);
H_matrix(1:3,7:9) = -eye(3);
H_matrix(4:6,4:6) = -eye(3);
H_matrix(4:6,13:15) = est_C_b_e_old*Skew_symmetric(Lab);

% 6. Set-up measurement noise covariance matrix assuming all components of
% GNSS position and velocity are independent and have equal variance.

% 7. Calculate Kalman gain using (3.21)
K_matrix = P_matrix_propagated * H_matrix' / (H_matrix * P_matrix_propagated * H_matrix' + R_matrix);

% 8. Formulate measurement innovations using (14.102)
% Lever arm from INS to GPS
delta_z(1:3,1) = GNSS_r_eb_e - est_r_eb_e_old - est_C_b_e_old*Lab - est_C_b_e_old*LC_KF_config.gps_correction;
delta_z(4:6,1) = GNSS_v_eb_e - est_v_eb_e_old ...
                  - est_C_b_e_old*cross(meas_omega_ib_b,Lab) ...
                  + Omega_ie*est_C_b_e_old*Lab;

% 9. Update state estimates using (3.24)
x_est_new = x_est_propagated + K_matrix * delta_z;
corrections = x_est_new;

% 10. Update state estimation error covariance matrix using (3.25)
P_matrix_new = (eye(15) - K_matrix * H_matrix) * P_matrix_propagated;
[~,p] = chol(P_matrix_new);
if p ~= 0
    P_matrix_new = nearestSPD(P_matrix_new);
end
% CLOSED-LOOP CORRECTION

% Correct attitude, velocity, and position using (14.7-9)
est_C_b_e_new = (eye(3) - Skew_symmetric(x_est_new(1:3))) * est_C_b_e_old;
est_v_eb_e_new = est_v_eb_e_old - x_est_new(4:6);
est_r_eb_e_new = est_r_eb_e_old - x_est_new(7:9);

% Return residual error
% residual = [GNSS_r_eb_e; GNSS_v_eb_e] - [est_r_eb_e_new; est_v_eb_e_new];
% residual = [est_r_eb_e_old; est_v_eb_e_old] - [est_r_eb_e_new; est_v_eb_e_new];
% Update IMU bias estimates
est_IMU_bias_new = est_IMU_bias_old + x_est_new(10:15);

% Ends

author: awerries

repo: kalman-localization


text:
function linplus_test5705 ( )

%*****************************************************************************80
%
%% TEST5705 tests R8SP_IJ_TO_K.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 7;
  n = 5;
  nz_num = 10;

  col = [ 2, 5, 1, 5, 1, 2, 3, 4, 4, 1 ];
  row = [ 1, 1, 2, 2, 4, 4, 4, 5, 6, 7 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST5705\n' );
  fprintf ( 1, '  R8SP_IJ_TO_K returns the R8SP index of (I,J).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %8d\n', m );
  fprintf ( 1, '  Matrix columns N = %8d\n', n );
  fprintf ( 1, '  Matrix nonzeros =  %8d\n', nz_num );

  check = r8sp_check ( m, n, nz_num, row, col );

  if ( ~check )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8SP_CHECK - Error!\n' );
    fprintf ( 1, '  The matrix is not in the proper sorted format.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '         I         J         K\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    for j = 1 : n
      k = r8sp_ij_to_k ( nz_num, row, col, i, j );
      fprintf ( 1, '  %8d  %8d  %8d\n', i, j, k );
    end
  end

  return
end


text:
function sched = par_beginschedule(tasks,varargin)
% Begin the scheduling of some set of tasks across a pool of (possibly remote) workers.
% Id = par_beginschedule(Tasks,Options...)
%
% Returns a scheduler handle to wait for and obtain results upon completion.
%
% In:
%   Tasks : cell array of tasks; each cell can be formatted as
%           * (evaluatable) string
%           * {function_handle, arg1,arg2,arg3, ...}
%           * struct('head',function_handle, 'parts',{{arg1,arg2,arg3, ...}})
%
%   Options...: optional name-value pairs, with possible names:
%               'engine': parallelization engine to use, can be one of:
%                         'global': select the global setting (tracking.parallel.engine)
%                         'local': do all computations locally, skipping serialization
%                         'BLS': use the BCILAB Scheduler (uses the resources specified in the pool argument) (default)
%                         'Reference': local reference implementation for testing BLS (using the same task serialization mechanism)
%                         'ParallelComputingToolbox': use the Mathworks Parallel Computing Toolbox (tm); uses resources allocated via the matlabpool command or a configuration file
%
%               'pool': pool of workers to consider for the BLS scheduler (default: 'global')
%                       if 'global', the global setting (tracking.parallel.pool) will be chosen
%                       (with the BLS engine, an empty pool implies local computation)
%
%               'policy': name of the scheduling policy function for the BLS (default: 'global')
%                         if 'global', the global setting (tracking.parallel.policy) will be chosen
%
%               'pushscope' : whether to "push" the current symbol scope (see hlp_scope and hlp_resolve)
%                             over the network (default: true)
%
%               'keep': keep this scheduler alive for later re-use (default: false)
%                       if false, the scheduler will be destroyed after use, and re-created during the next run
%
% Out:
%   Id : output handle; used to collect the results
%
% See also:
%   par_endschedule, par_worker, par_schedule
%
% Example:
%   % schedule two computations across a pool of some IP:port addresses (assuming that MATLAB is running there,
%   % and is executing the par_worker() function
%   id = par_beginschedule({'sin(randn(10))','exp(randn(10))'}, 'pool',{'192.168.1.1:23547','192.168.1.2:23547','192.168.1.2:23548','192.168.1.3:23547'});
%   ... optionally do something in the meantime
%   results = par_endschedule(id);
%
% Expert note:
%  The 'keep' option is easy to use with the wrapper function par_schedule; otherwise, the following holds:
%   * if passed as true to par_beginschedule, it *must* also be passed as true to par_endschedule
%   * nested schedules are not allowed if they use the same worker pool
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-08-29

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

% read options
opts = arg_define('allow-unlisted-names',varargin, ...
    arg_norep({'tasks','Tasks'},[],[],'Tasks to execute.','type','expression'), ...
    arg({'engine','ParallelEngine'},'global',{'global','local','BLS','Reference','ParallelComputingToolbox'}, 'Parallel engine to use. This can either be one of the supported parallel engines (BLS for BCILAB Scheduler, Reference for a local reference implementation, and ParallelComputingToolbox for a PCT-based implementation), or local to skip parallelization altogether, or global to select the currently globally selected setting (in the global tracking variable).'), ...
    arg({'pool','WorkerPool'},'global',[], 'Worker pool to use. This is typically a cell array, but can also be the string ''gobal'', which stands for the currently globally set up worker pool (see global tracking variable).','type','expression'), ...
    arg({'policy','ReschedulingPolicy'},'global',[], 'Rescheduling policy. This is the name of the rescheduling policy function that controls if and when tasks are being rescheduled. If set to global, the current global setting will be used.'), ...
    arg({'verbosity','Verbosity'},'global',{'global','0','1','2','3'}, 'Verbosity level. Verbosity level for the scheduling process, or global to use the global setting.'), ...
    arg({'receiver_backlog','ReceiverBacklog'},5,uint32([0 1 7 10]), 'Max backlogged connections at receiver. Maximum number of connections that can be in the backlog of the result receiver. '), ...
    arg({'receiver_timeout','ReceiverTimeout'},1,[0 1 60 3600], 'Receiver timeout. In seconds. When a worker that has previously announced that it is ready to transmit results does not start transmitting results within this time period after pinged, the connection will be discarded and the worker has to re-connect.'), ...
    arg({'reschedule_interval','RescheduleInterval'},5,[0 1 15 3600], 'Rescheduling interval. In seconds. Time between periodic checks if jobs have to be re-scheduled.'), ...
    arg({'keep','RetainScheduler'},false,[], 'Retain scheduler instance. If true, the scheduler instance will be reused on the next call to par_beginschedule. Can be more efficient, but does not allow for nested calls to par_beginschedule.'), ...
    arg({'pushscope','PushScope'},true,[], 'Push variable scope. If true, the current variable scope (established via hlp_scope) will be pushed over the network; necessary for some features to run seamlessly.'));

% max length of java arrays that we construct (technically 2^31-1 minus some slack but playing it safe)
max_java_array_len = 2^30;

% look up global settings, if requested
if strcmp(opts.engine,'global')
    opts.engine = par_globalsetting('engine'); end
if strcmp(opts.pool,'global')
    opts.pool = par_globalsetting('pool'); end
if strcmp(opts.policy,'global')
    opts.policy = par_globalsetting('policy'); end
if strcmp(opts.verbosity,'global')
    opts.verbosity = par_globalsetting('verbosity'); end
if ischar(opts.verbosity)
    opts.verbosity = str2num(opts.verbosity); end
if isa(opts.policy,'function_handle')
    opts.policy = char(opts.policy); end
if strcmp(opts.engine,'BLS') && isempty(opts.pool)
    opts.engine = 'local'; end

% validate options
if ~iscellstr(opts.pool)
    error('The given worker pool must be a cell array of strings, but was: %s',hlp_tostring(opts.pool)); end
if ~ischar(opts.policy)
    error('The given rescheduling policy must be a string, but was: %s',hlp_tostring(opts.policy)); end
if ~exist(opts.policy,'file')
    error('The given rescheduling policy refers to a non-existing file: %s',opts.policy); end

% canonlicalize task format to {function-handle,arguments...}
for t=1:length(tasks)
    if isfield(tasks{t},{'head','parts'})
        % task given as Mathematica-style expression struct (see expressions/exp_*)
        tasks{t} = [{tasks{t}.head} tasks{t}.parts];
    elseif ischar(tasks{t})
        % task given as string
        tasks{t} = {@eval,tasks{t}};
    elseif ~(iscell(tasks{t}) && ~isempty(tasks{t}) && isa(tasks{t}{1},'function_handle'))
        % incorrect task format...
        error('Unsupported task format, please see documentation: %s',hlp_tostring(tasks{t},10000));
    end
end

% push current symbol context over the network
if opts.pushscope && ~strcmp(opts.engine,'local')
    % get the current scope
    scope = hlp_resolveall;
    % and wrap a hlp_scope() around the task
    for t=1:length(tasks)
        tasks{t} = [{@hlp_scope, scope} tasks{t}]; end
end

% create a scheduler (Java code, see dependencies/Scheduling-*)
if strcmp(opts.engine,'BLS')
    if opts.verbosity>0 %#ok<*ST2NM>
        fprintf('Creating scheduler...\n'); end
    % trim the 'pid@' part from opts.pool (the scheduler expects the 'host:port' format)
    pool = opts.pool;
    for p=1:length(pool)
        pos = pool{p}=='@';
        if any(pos)
            pool{p} = pool{p}(find(pos,1)+1:end); end
    end
    if opts.keep
        tmp = hlp_microcache('schedulers',@(varargin)Scheduler(varargin{:}),pool,opts.policy,'par_accept_results',opts.receiver_backlog,round(1000*opts.receiver_timeout),round(1000*opts.reschedule_interval),opts.verbosity,length(tasks));
    else
        tmp = Scheduler(pool,opts.policy,'par_accept_results',opts.receiver_backlog,round(1000*opts.receiver_timeout),round(1000*opts.reschedule_interval),opts.verbosity,length(tasks));
    end
    sched = struct('sched',{tmp},'finisher',{onCleanup(@()tmp.clear())});
end

% serialize tasks for network transmission (and prepend order id)
if any(strcmp(opts.engine,{'BLS','Reference'}))
    if opts.verbosity>0
        t0=tic; fprintf('Serializing %d tasks...\n', length(tasks)); end
    for t=1:length(tasks)
        tasks{t} = hlp_serialize([{t} tasks{t}]); end
    if opts.verbosity>0
        fprintf('Tasks serialized (%.1f seconds).\n', toc(t0)); end
end

% submit tasks for execution
switch opts.engine
    case 'local'
        sched = cell(1,length(tasks));
        for t=1:length(tasks)
            sched(t) = {{t,tasks{t}{1}(tasks{t}{2:end})}}; end
    case 'ParallelComputingToolbox'
        sched = {};
        parfor t=1:length(tasks)
            sched(t) = {{t,tasks{t}{1}(tasks{t}{2:end})}}; end
    case 'BLS'
        % over the scheduler
        fprintf('Converting tasks for scheduler...');
        bam = ByteArrayMaker();
        for t=1:length(tasks)
            % we need to efficiently split the task into chunks no larger than max_java_array_len
            % (since the data can easily exceed the max java array length)
            if length(tasks{t}) > max_java_array_len
                num_blocks = ceil(length(tasks{t})/max_java_array_len);
                rest_size = mod(length(tasks{t}),max_java_array_len);
                sizes = [max_java_array_len*ones(1,num_blocks-1) rest_size];
            else
                % we split into at least 2 blocks to not get confused with the gnarly
                % construction/marshaling of Java objects in MATLAB
                sizes = [floor(length(tasks{t})/2), length(tasks{t})-floor(length(tasks{t})/2)];
            end
            blocks = cell(length(sizes),1);
            [blocks{:}] = chopdeal(tasks{t},sizes);
            for b=1:length(blocks)
                blocks{b} = bam.makeByteArray(blocks{b}); end
            tasks{t} = [blocks{:}];
        end
        tasklist = java.util.ArrayList;
        for t=1:length(tasks)
            tasklist.add(tasks{t}); end
        fprintf('done.\n');
        if opts.verbosity>0
            fprintf('Submitting tasks to scheduler...\n'); end
        sched.sched.submit(tasklist);
    case 'Reference'
        % evaluate locally, but go through the same evaluation function as the BLS workers
        unencoded_results = cell(1,length(tasks)); % holds unencoded results for debugging
        for t=1:length(tasks)
            [tasks{t},unencoded_results{t}] = par_evaluate(tasks{t}); end
        % return the collected result in sched
        sched = struct('ReferenceResults',{tasks},'ReferenceResultsUnencoded',{unencoded_results});
    otherwise
        error('Unsupported parallelization engine selected: %s',hlp_tostring(opts.engine));
end

repo: REST


text:
function fx = p40_fun ( n, x )

%*****************************************************************************80
%
%% P40_FUN evaluates the integrand for problem 40.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P40_PARAM_SET.
%
%    The integrand has a singularity at an internal point ( x = pi/4 ).
%
%    The suggested range for ALPHA is -0.8 through 2.1.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    ( abs ( x - pi/4 ) )^alpha
%
%  Exact Integral:
%
%    ( (1-pi/4)^(alpha+1) + (pi/4)^(alpha+1) ) / ( alpha + 1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 83.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  alpha = p40_param_get ( );

  fx = ( abs ( x - 0.25 * pi ) ).^alpha;

  return
end


text:
 function [Cone,EndPlate1,EndPlate2] = Cone(hAxis,X1,X2,R,n,cyl_color,closed,lines)
%
% This function constructs a cylinder connecting two center points
%
% Usage :
% [Cone,EndPlate1,EndPlate2] = Cone(X1,X2,R,n,cyl_color,closed,lines)
%
% Cone-------Handle of the cone
% EndPlate1------Handle of the Starting End plate
% EndPlate2------Handle of the Ending End plate
% X1 and X2 are the 3x1 vectors of the two points
% R is the radius of the cylinder/cone R(1) = start radius, R(2) = end radius
% n is the no. of elements on the cylinder circumference (more--> refined)
% cyl_color is the color definition like 'r','b',[0.52 0.52 0.52]
% closed=1 for closed cylinder or 0 for hollow open cylinder
% lines=1 for displaying the line segments on the cylinder 0 for only
% surface
%
% Typical Inputs
% X1=[10 10 10];
% X2=[35 20 40];
% r=[1 5];
% n=20;
% cyl_color='b';
% closed=1;
%
% NOTE: There is a MATLAB function "cylinder" to revolve a curve about an
% axis. This "Cylinder" provides more customization like direction and etc

% Calculating the length of the Cone
length_cyl=norm(X2-X1);

% Creating 2 circles in the YZ plane
t=linspace(0,2*pi,n)';
xa2=R(1)*cos(t);
xa3=R(1)*sin(t);
xb2=R(2)*cos(t);
xb3=R(2)*sin(t);

% Creating the points in the X-Direction
x1=[0 length_cyl];

% Creating (Extruding) the cylinder points in the X-Directions
xx1=repmat(x1,length(xa2),1);
xx2=[xa2 xb2];%xx2=repmat(x2,1,2);
xx3=[xa3 xb3];%xx3=repmat(x3,1,2);

% Drawing two filled cirlces to close the cylinder
if closed==1
    EndPlate1=fill3(xx1(:,1),xx2(:,1),xx3(:,1),'r', 'Parent', hAxis);
    EndPlate2=fill3(xx1(:,2),xx2(:,2),xx3(:,2),'r', 'Parent', hAxis);
end

% Plotting the cylinder along the X-Direction with required length starting
% from Origin
Cone=mesh(hAxis,xx1,xx2,xx3, 'Parent', hAxis);

% Defining Unit vector along the X-direction
unit_Vx=[1 0 0];

% Calulating the angle between the x direction and the required direction
% of Cone through dot product
angle_X1X2=acos( dot( unit_Vx,(X2-X1) )/( norm(unit_Vx)*norm(X2-X1)) )*180/pi;

% Finding the axis of rotation (single rotation) to roate the Cone in
% X-direction to the required arbitrary direction through cross product
axis_rot=cross([1 0 0],(X2-X1) );

% Rotating the plotted Cone and the end plate circles to the required
% angles
if angle_X1X2~=0 % Rotation is not needed if required direction is along X
    rotate(Cone,axis_rot,angle_X1X2,[0 0 0])
    if closed==1
        rotate(EndPlate1,axis_rot,angle_X1X2,[0 0 0])
        rotate(EndPlate2,axis_rot,angle_X1X2,[0 0 0])
    end
end

% Till now Cone has only been aligned with the required direction, but
% position starts from the origin. so it will now be shifted to the right
% position
if closed==1
    set(EndPlate1,'XData',get(EndPlate1,'XData')+X1(1))
    set(EndPlate1,'YData',get(EndPlate1,'YData')+X1(2))
    set(EndPlate1,'ZData',get(EndPlate1,'ZData')+X1(3))

    set(EndPlate2,'XData',get(EndPlate2,'XData')+X1(1))
    set(EndPlate2,'YData',get(EndPlate2,'YData')+X1(2))
    set(EndPlate2,'ZData',get(EndPlate2,'ZData')+X1(3))
end
set(Cone,'XData',get(Cone,'XData')+X1(1))
set(Cone,'YData',get(Cone,'YData')+X1(2))
set(Cone,'ZData',get(Cone,'ZData')+X1(3))

% Setting the color to the Cone and the end plates
set(Cone,'AmbientStrength',1,'FaceColor',cyl_color,'FaceLighting','gouraud', 'Parent', hAxis);%,'EdgeColor','none')
if closed==1
    set([EndPlate1 EndPlate2],'AmbientStrength',1,'FaceColor',cyl_color,'FaceLighting','gouraud', 'Parent', hAxis);%,'EdgeColor','none')
else
    EndPlate1=[];
    EndPlate2=[];
end

% If lines are not needed making it disapear
if lines==0
    set(Cone,'EdgeAlpha',0)
end

%shading faceted % faceted flat interp;
%camlight;
%light;
%lighting gouraud; %flat gouraud phong none
% material shiny; %shiny dull metal
%colormap(bone)

%camlight  headlight;
%light('Style','local','Position',[720 0 500]);
%light('Style','local','Position',[0 480 500]);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef MappingComputer < handle

    properties (Access = private)
        LHS
        RHS
    end

    properties (Access = private)
        meshDisc
        mesh
        orientation
        interp
        interpolator
    end

    methods (Access = public)

        function obj = MappingComputer(cParams)
            obj.init(cParams);
        end

        function uF = compute(obj)
            obj.computeLHS();
            obj.computeRHS();
            uC = obj.solveSystem();
            In = obj.interpolator;
            u  = In*uC;
            uV(1,:,:) = reshape(u,obj.mesh.nnodeElem,[]);
            s.mesh    = obj.mesh;
            s.fValues = uV;
            uF = P1DiscontinuousFunction(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh        = cParams.mesh;
            obj.orientation  = cParams.orientation;
            obj.interpolator = cParams.interpolator;
            obj.meshDisc     = obj.mesh.createDiscontinuousMesh();
        end

        function computeLHS(obj)
            K = obj.computeStiffnessMatrix();
            In = obj.interpolator;
            Kn = In'*K*In;
            obj.LHS = Kn;
        end

        function K = computeStiffnessMatrix(obj)
            s.mesh = obj.mesh;
            s.type = 'StiffnessMatrix';
            s.fun  = P1DiscontinuousFunction.create(obj.mesh,1);
            lhs2 = LHSintegrator.create(s);
            K = lhs2.compute();
        end

        function computeRHS(obj)
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature('QUADRATIC');
%             fG = obj.orientation.evaluate(q.posgp);
            s.mesh  = obj.meshDisc;
            s.type = 'ShapeDerivative';
            s.quadratureOrder = q.order;
            rhs  = RHSintegrator.create(s);
            rhsF = rhs.compute(obj.orientation);
            In = obj.interpolator;
            rhsV = In'*rhsF.fValues;
            obj.RHS = rhsV;
        end

        function u = solveSystem(obj)
            a.type = 'DIRECT';
            s = Solver.create(a);
            u = s.solve(obj.LHS,obj.RHS);
            u = u(1:end);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Unpacks an object of any complexity (structure, cell array ect.)
% from a  serialized form.
%
% mex-function, written by Dima Mihaylov (dmih@in-solve.ru)
%
% See also SERIALIZE

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef SplitterInConnectedComponents < handle

    properties (Access = private)
        edges
        adjancyMatrix
        graphV
        connectedComp
    end

    properties (Access = private)
        faces
    end

    methods (Access = public)

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

        function c = split(obj)
            obj.computeEdges();
            obj.computeAdjacencyMatrix();
            obj.computeGraph();
            obj.computeConnectedComponents();
            c = obj.connectedComp;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.faces = cParams.faces;
        end

        function computeEdges(obj)
            s.nodesByElem = obj.faces;
            edge = EdgesConnectivitiesComputer(s);
            edge.compute();
            obj.edges = edge.nodesInEdges;
        end

        function [A] = computeAdjacencyMatrix(obj)
            nodeA = obj.edges(:,1);
            nodeB = obj.edges(:,2);
            A = sparse(nodeA,nodeB,1);
            obj.adjancyMatrix = A+A';
        end

        function computeGraph(obj)
            A = obj.adjancyMatrix;
            g = graph(A,'omitselfloops');
            obj.graphV = g;
        end

        function computeConnectedComponents(obj)
            g = obj.graphV;
            comp = conncomp(g);
            obj.connectedComp = comp;
        end

    end

end


text:
function phiFaceAverage = upwindMean3D(phi, u)
% This function gets the value of the field variable phi defined
% over the MeshStructure and calculates the upwind average on
% the cell faces, based on the direction of the velocity vector for a uniform mesh.
%
% SYNOPSIS:
%   phiFaceAverage = upwindMean3D(phivar, u)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;
uy = u.yvalue;
uz = u.zvalue;

% check the size of the variable and the mesh dimension
Nxyz = phi.domain.dims;
Nx = Nxyz(1); Ny = Nxyz(2); Nz = Nxyz(3);

% assign to a temp variable for boundary corrections
phi_tmp = phi.value;

% correct the value of phi at the boundary (calculation trick)
% assign the value of the left boundary to the left ghost cells
phi_tmp(1,:,:) = (phi.value(1,:,:)+phi.value(2,:,:))/2;
% assign the value of the right boundary to the right ghost cells
phi_tmp(end,:,:) = (phi.value(end,:,:)+phi.value(end-1,:,:))/2;
% assign the value of the bottom boundary to the bottom ghost cells
phi_tmp(:,1,:) = (phi.value(:,1,:)+phi.value(:,2,:))/2;
% assign the value of the top boundary to the top ghost cells
phi_tmp(:,end,:) = (phi.value(:,end,:)+phi.value(:,end-1,:))/2;
% assign the value of the back boundary to the back ghost cells
phi_tmp(:,:,1) = (phi.value(:,:,1)+phi.value(:,:,2))/2;
% assign the value of the front boundary to the front ghost cells
phi_tmp(:,:,end) = (phi.value(:,:,end)+phi.value(:,:,end-1))/2;

% calculate the average value
xvalue = (ux>0).*phi_tmp(1:Nx+1,2:Ny+1,2:Nz+1)+ ...
                        (ux<0).*phi_tmp(2:Nx+2,2:Ny+1,2:Nz+1)+ ...
                        0.5*(ux==0).*(phi.value(1:Nx+1,2:Ny+1,2:Nz+1)+phi.value(2:Nx+2,2:Ny+1,2:Nz+1));
yvalue = (uy>0).*phi_tmp(2:Nx+1,1:Ny+1,2:Nz+1)+ ...
                        (uy<0).*phi_tmp(2:Nx+1,2:Ny+2,2:Nz+1)+ ...
                        0.5*(uy==0).*(phi.value(2:Nx+1,1:Ny+1,2:Nz+1)+phi.value(2:Nx+1,2:Ny+2,2:Nz+1));
zvalue = (uz>0).*phi_tmp(2:Nx+1,2:Ny+1,1:Nz+1)+ ...
                        (uz<0).*phi_tmp(2:Nx+1,2:Ny+1,2:Nz+2)+ ...
                        0.5*(uz==0).*(phi.value(2:Nx+1,2:Ny+1,1:Nz+1)+phi.value(2:Nx+1,2:Ny+1,2:Nz+2));
phiFaceAverage=FaceVariable(phi.domain, xvalue, yvalue, zvalue);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%ADDPREPROC Add commands in the PREPROC field of a datafile
%
%   A = ADDPREPROC(A,PREPROC,PARS,OUTSIZE)
%
% INPUT
%   A       - Datafile
%   PREPROC - String with preprocessing command
%   PARS    - Cell array with parameters (default empty)
%   OUTSIZE - Output size of objects
%
% OUTPUT
%   A       - Datafile
%
% DESCRIPTION
% Extends the structure array of preprocessing commands, A.PREPROC( )
% with one element, having two fields:
% A.PREPROC(N).PREPROC = PREPROC
% A.PREPROC(N).PARS = PARS
%
% Preprocessing can only be defined for raw datafiles that do not contain
% datasets. It is an error to define preprocessing for a datafile that
% points to MAT files. In that case datasets are expected and preprocessing
% is superfluous.
%
% The first command in the PREPROC field should always be a file read
% command. The DATAFILE constructor stores by default IMREAD. Use
% SETPREPROC to clear the PREPROC field and replace it by another
% command.
%
% PRTools needs to know the size of the output objects. If not supplied in
% the call, it is found by processing the first object. This may take some
% time.
%
% SEE ALSO (PRTools Guide)
% DATAFILES, SETPREPROC, ADDPOSTPROC, DATASETS

author: marianux

repo: ecg-kit

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


text:
classdef TestBitwiseOr
    %TestBitwiseOr

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_rgb_images
            img1 = cv.imread(TestBitwiseOr.im, 'ReduceScale',2);
            img2 = cv.cvtColor(img1, 'RGB2HSV');

            % rectangular mask
            [h,w,~] = size(img1);
            mask = false([h w]);
            mask(100:h-100,100:w-100) = true;

            out = cv.bitwise_or(img1, img2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_or(img1, img2);
            assert(isequal(out, expected));

            out = cv.bitwise_or(img1, img2, 'Mask',mask);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_or(img1, img2, mask);
            assert(isequal(out, expected));

            out = cv.bitwise_or(img1, img2, 'Mask',mask, 'Dest',img1);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_or(img1, img2, mask, img1);
            assert(isequal(out, expected));
        end

        function test_float_images
            img1 = cv.imread(TestBitwiseOr.im, 'ReduceScale',2);
            img1 = single(img1) ./ 255;
            img2 = cv.cvtColor(img1, 'RGB2HSV');
            img2(:,:,1) = img2(:,:,1) ./ 360;

            % circular mask
            [h,w,~] = size(img1);
            [X,Y] = meshgrid(1:w,1:h);
            c = fix([w h]/2); r = 50;
            mask = ((X-c(1)).^2 + (Y-c(2)).^2) < r^2;

            out = cv.bitwise_or(img1, img2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_or(img1, img2);
            assert(isequaln(out, expected));

            out = cv.bitwise_or(img1, img2, 'Mask',mask);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_or(img1, img2, mask);
            assert(isequaln(out, expected));

            out = cv.bitwise_or(img1, img2, 'Mask',mask, 'Dest',img1);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_or(img1, img2, mask, img1);
            assert(isequaln(out, expected));
        end

        function test_vectors
            A = uint8([250 253 0]);
            B = uint8([1 2 255]);
            mask = [true false true];

            % uint8([251 255 255])
            out = cv.bitwise_or(A, B);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitor(A, B);
            assert(isequal(out, expected));

            % uint8([251 0 255])
            out = cv.bitwise_or(A, B, 'Mask',mask);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitor(A, B);
            expected(~mask) = 0;
            assert(isequal(out, expected));

            % uint8([251 253 255])
            out = cv.bitwise_or(A, B, 'Mask',mask, 'Dest',A);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitor(A, B);
            expected(~mask) = A(~mask);
            assert(isequal(out, expected));
        end

        function test_error_argnum
            try
                cv.bitwise_or();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

function out = my_bitwise_or(src1, src2, mask, dst)
    %MY_BITWISE_OR  Similar to cv.bitwise_or using core MATLAB functions

    if nargin < 3, mask = true(size(src1,1), size(src1,2)); end
    if nargin < 4, dst = zeros(size(src1), class(src1)); end

    % calculate bitwise OR
    if isinteger(src1)
        out = bitor(src1, src2);
    elseif isfloat(src1)
        if isa(src1, 'double')
            klass = 'uint64';
        elseif isa(src1, 'single')
            klass = 'uint32';
        end
        out = bitor(typecast(src1(:), klass), typecast(src2(:), klass));
        out = reshape(typecast(out, class(src1)), size(src1));
    end

    % apply masking
    for k=1:size(out,3)
        out_slice = out(:,:,k);
        dst_slice = dst(:,:,k);
        out_slice(~mask) = dst_slice(~mask);
        out(:,:,k) = out_slice;
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function rez = datashift2(rez, do_correction)

NrankPC = 6;
[wTEMP, wPCA]    = extractTemplatesfromSnippets(rez, NrankPC);
rez.wTEMP = gather(wTEMP);
rez.wPCA  = gather(wPCA);

ops = rez.ops;

% The min and max of the y and x ranges of the channels
ymin = min(rez.yc);
ymax = max(rez.yc);
xmin = min(rez.xc);
xmax = max(rez.xc);

dmin = median(diff(unique(rez.yc)));
fprintf('vertical pitch size is %d \n', dmin)
rez.ops.dmin = dmin;
rez.ops.yup = ymin:dmin/2:ymax; % centers of the upsampled y positions

% dminx = median(diff(unique(rez.xc)));
yunq = unique(rez.yc);
mxc = zeros(numel(yunq), 1);
for j = 1:numel(yunq)
    xc = rez.xc(rez.yc==yunq(j));
    if numel(xc)>1
       mxc(j) = median(diff(sort(xc)));
    end
end
dminx = max(5, median(mxc));
fprintf('horizontal pitch size is %d \n', dminx)

rez.ops.dminx = dminx;
nx = round((xmax-xmin) / (dminx/2)) + 1;
rez.ops.xup = linspace(xmin, xmax, nx); % centers of the upsampled x positions
disp(rez.ops.xup)

if  getOr(rez.ops, 'nblocks', 1)==0
    rez.iorig = 1:rez.temp.Nbatch;
    return;
end

% binning width across Y (um)
dd = 5;
% min and max for the range of depths
dmin = ymin - 1;
dmax  = 1 + ceil((ymax-dmin)/dd);
disp(dmax)

spkTh = 10; % same as the usual "template amplitude", but for the generic templates

% Extract all the spikes across the recording that are captured by the
% generic templates. Very few real spikes are missed in this way.
[st3, rez] = standalone_detector(rez, spkTh);
%%

% detected depths
% dep = st3(:,2);
% dep = dep - dmin;

Nbatches      = rez.temp.Nbatch;
% which batch each spike is coming from
batch_id = st3(:,5); %ceil(st3(:,1)/dt);

% preallocate matrix of counts with 20 bins, spaced logarithmically
F = zeros(dmax, 20, Nbatches);
for t = 1:Nbatches
    % find spikes in this batch
    ix = find(batch_id==t);

    % subtract offset
    dep = st3(ix,2) - dmin;

    % amplitude bin relative to the minimum possible value
    amp = log10(min(99, st3(ix,3))) - log10(spkTh);

    % normalization by maximum possible value
    amp = amp / (log10(100) - log10(spkTh));

    % multiply by 20 to distribute a [0,1] variable into 20 bins
    % sparse is very useful here to do this binning quickly
    M = sparse(ceil(dep/dd), ceil(1e-5 + amp * 20), ones(numel(ix), 1), dmax, 20);

    % the counts themselves are taken on a logarithmic scale (some neurons
    % fire too much!)
    F(:, :, t) = log2(1+M);
end

%%
% determine registration offsets
ysamp = dmin + dd * [1:dmax] - dd/2;
[imin,yblk, F0, F0m] = align_block2(F, ysamp, ops.nblocks);

if isfield(rez, 'F0')
    d0 = align_pairs(rez.F0, F0);
    % concatenate the shifts
    imin = imin - d0;
end

%%
if getOr(ops, 'fig', 1)
    figure;
    set(gcf, 'Color', 'w')

    % plot the shift trace in um
    plot(imin * dd)
    box off
    xlabel('batch number')
    ylabel('drift (um)')
    title('Estimated drift traces')
    drawnow

    figure;
    set(gcf, 'Color', 'w')
    % raster plot of all spikes at their original depths
    st_shift = st3(:,2); %+ imin(batch_id)' * dd;
    for j = spkTh:100
        % for each amplitude bin, plot all the spikes of that size in the
        % same shade of gray
        ix = st3(:, 3)==j; % the amplitudes are rounded to integers
        plot(st3(ix, 1)/ops.fs, st_shift(ix), '.', 'color', [1 1 1] * max(0, 1-j/40)) % the marker color here has been carefully tuned
        hold on
    end
    axis tight
    box off

    xlabel('time (sec)')
    ylabel('spike position (um)')
    title('Drift map')

end
%%
% convert to um
dshift = imin * dd;

% this is not really used any more, should get taken out eventually
[~, rez.iorig] = sort(mean(dshift, 2));

if do_correction
    % sigma for the Gaussian process smoothing
    sig = rez.ops.sig;
    % register the data batch by batch
    dprev = gpuArray.zeros(ops.ntbuff,ops.Nchan, 'single');
    for ibatch = 1:Nbatches
        dprev = shift_batch_on_disk2(rez, ibatch, dshift(ibatch, :), yblk, sig, dprev);
    end
    fprintf('time %2.2f, Shifted up/down %d batches. \n', toc, Nbatches)
else
    fprintf('time %2.2f, Skipped shifting %d batches. \n', toc, Nbatches)
end
% keep track of dshift
rez.dshift = dshift;
% keep track of original spikes
rez.st0 = st3;

rez.F = F;
rez.F0 = F0;
rez.F0m = F0m;

% next, we can just run a normal spike sorter, like Kilosort1, and forget about the transformation that has happened in here

%%

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function gmt_xy(hLine, sOutput, sPrefix)

% Get attributes
sLineStyle = lower(get(hLine, 'LineStyle'));
sMarker = get(hLine, 'Marker');

if ~((strcmp(sLineStyle, 'none'))  &&  (strcmp(sMarker, 'none')))

  % Get data
  mXData = get(hLine, 'XData');
  mYData = get(hLine, 'YData');

  % Transform data and save it
  [nRow, nCnt] = size(mXData);
  if nRow == 1
    mXData = mXData';
  end
  [nRow, nCnt] = size(mYData);
  if nRow == 1
    mYData = mYData';
  end
  mXYValues = [mXData mYData];
  save([sPrefix '_xy.tmp'], 'mXYValues', '-ascii');

  hFile = fopen(sOutput, 'a');
  fprintf(hFile, '# Plot xy\n');

  unix(['sed -f xy.sed ' sPrefix '_xy.tmp > ' sPrefix '_xy.dat']);
  unix(['rm -f ' sPrefix '_xy.tmp']);

  % Get attributes
  sLineStyle = lower(get(hLine, 'LineStyle'));
  fLineWidth = get(hLine, 'LineWidth');
  sMarker = get(hLine, 'Marker');
  fMarkerSize = get(hLine, 'MarkerSize');
  vColor = get(hLine, 'Color');

  if ~strcmp(sLineStyle, 'none')
    sS = '';
    sW = sprintf('-W%s/%s/%s/%s ', num2str(fLineWidth), num2str(vColor(1)*255), num2str(vColor(2)*255), num2str(vColor(3)*255));
    sM = '-M ';
    fprintf(hFile, 'psxy %s $default %s %s %s -O -K >> $output\n', [sPrefix '_xy.dat'], sW, sM, sS);
  end
  if ~strcmp(sMarker, 'none')
    switch sMarker
    case '^',         % triangle upwards
      sSymbol = 't';
    case '.',         % point
      sSymbol = 'p';
    case 'hexagram',  % hexagram
      sSymbol = 'a';
    case 'o',         % circle
      sSymbol = 'c';
    case 'square',    % square
      sSymbol = 's';
    case 'x',         % cross
      sSymbol = 'x';
    case '*',         % asterix
      sSymbol = 'a';
    otherwise         % point as default
      sSymbol = 'p';
    end
    vEdgeColor = get(hLine, 'MarkerEdgeColor');
    if strcmp(vEdgeColor, 'auto')
      vEdgeColor = get(hLine, 'Color');
    end
    vFaceColor = get(hLine, 'MarkerFaceColor');
    if ~strcmp(vFaceColor, 'none')
      sG = sprintf('-G%s/%s/%s ', num2str(vFaceColor(1)*255), num2str(vFaceColor(2)*255), num2str(vFaceColor(3)*255));
    else
      sG = '';
    end
    sS = ['-S' sSymbol num2str(fMarkerSize) 'p ']

    sW = sprintf('-W%s/%s/%s/%s ', num2str(fLineWidth), num2str(vEdgeColor(1)*255), num2str(vEdgeColor(2)*255), num2str(vEdgeColor(3)*255));
    fprintf(hFile, 'psxy %s $default %s %s %s -O -K >> $output\n', [sPrefix '_xy.dat'], sW, sG, sS);
  end
  fclose(hFile);
end

author: CelsoReyes


text:
function dy = sparseGalerkin(t,y,ahat,polyorder,usesine)
% Copyright 2015, All Rights Reserved
% Code by Steven L. Brunton
% For Paper, "Discovering Governing Equations from Data:
%        Sparse Identification of Nonlinear Dynamical Systems"
% by S. L. Brunton, J. L. Proctor, and J. N. Kutz

yPool = poolData(y',length(y),polyorder,usesine);
dy = (yPool*ahat)';

author: eurika-kaiser

repo: SINDY-MPC


text:
% Author: Sergiy Iglin
% e-mail: iglin@kpi.kharkov.ua
% or: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru
clear all % Rubik Cube Game Main File
scrsz = get(0,'ScreenSize');
pl = (scrsz(3)-460)/2; % left bound
pu = (scrsz(4)-460)/2; % upper bound
Cube = InitCube; % initialize the cube
hFig = figure('Position',[pl pu 460 460],...
  'Color',get(0,'DefaultUicontrolBackgroundColor'),...
  'Resize','off','MenuBar','none',...
  'NumberTitle','off','Name','Rubik Cube Game',...
  'WindowButtonMotionFcn','ChangePointer(hFig,Cube);');
hAxes = axes('Units','pixels','Position',[20 20 420 420],...
  'CameraPosition',[4.5 4 3.5],...
  'XLim',[-3.7 3.7],'YLim',[-3.7 3.7],'ZLim',[-3.7 3.7],...
  'XTickLabel','','YTickLabel','','ZTickLabel','',...
  'XTick',[],'YTick',[],'ZTick',[]);
hMenuGame = uimenu(hFig,'Label','Game');
hMenuNewGame = uimenu(hMenuGame,'Label','New Game',...
  'Callback','Cube=CreateNewGame(hAxes,Cube);');
hMenuExit = uimenu(hMenuGame,'Label','Exit',...
  'Callback','delete(hFig);','Separator','on');
hMenuHelp = uimenu(hFig,'Label','Help');
hMenuHow = uimenu(hMenuHelp,'Label','How to build the Rubik cube?',...
  'Callback','web jjorg.chem.unc.edu/personal/monroe/cube/rubik.html -browser');
hMenuAuthor = uimenu(hMenuHelp,'Label','Author','Separator','on',...
  'Callback','web iglin.exponenta.ru/All/Rubik.html -browser');
Cube = CreateNewGame(hAxes,Cube); % create new game and draw

author: Sable

repo: mcbench-benchmarks