text:
%sim_evolve.m
%Jamie Near, 2014.
%
% USAGE:
% d_out = sim_evolve(d_in,H,t)
%
% DESCRIPTION:
% This function simulates free evolution of the spin system under the
% effects of chemical shift and scalar coupling.
%
% INPUTS:
% d_in      = input density matrix structure.
% H         = Hamiltonian operator structure.
% t         = duration of evolution (s)
%
% OUTPUTS:
% d_out     = output density matrix following free evolution.

function d_out = sim_evolve(d_in,H,t)

for n=1:length(H) %JN - loop through the different parts of the spin-system
    p=expm(1i*H(n).HAB*t);
    d_out{n} = p' * d_in{n} * p;  %Faster than doing expm twice.  Suggested by Martin Froeling.
end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
clear all; close all;
figure(1)
% underdetermined
n=20; m=100
A=rand(n,m); b=rand(n,1);

cvx_begin;
variable x2(m)
minimize( norm(x2,2) );
subject to
A*x2 == b;
cvx_end;

cvx_begin;
variable x1(m)
minimize( norm(x1,1) );
subject to
A*x1 == b;
cvx_end;

subplot(3,1,1), bar(x2)
subplot(3,1,2), bar(x1)
subplot(3,2,5), hist(x2,40)
subplot(3,2,6), hist(x1,40)

subplot(3,1,1), axis([0 100 -.2 .2]), subplot(3,1,2), axis([0 100 -.4 .4])
subplot(3,2,5), axis([-.15 .15 0 82]), set(gca,'Ytick',[0 40 80],'Xtick',[-0.1 0 0.1])
h = findobj(gca,'Type','patch'); h.FaceColor = [0.6 0.6 0.6]; h.EdgeColor = 'k';
subplot(3,2,6), axis([-.15 .15 0 82]), set(gca,'Ytick',[0 40 80],'Xtick',[-0.1 0 0.1])
h = findobj(gca,'Type','patch'); h.FaceColor = [0.6 0.6 0.6]; h.EdgeColor = 'k';

%%
clear A
clear b
clear x
figure(2)

% overdetermined
n=500; m=100;
A=rand(n,m);
b=rand(n,1);
xdag=pinv(A)*b;

lam=[0 0.1 0.5];
for j=1:3

    cvx_begin;
    variable x(m)
    minimize( norm(A*x-b,2) + lam(j)*norm(x,1) );
    cvx_end;

    subplot(4,1,j),bar(x)
    subplot(4,3,9+j), hist(x,20)
end

for j=1:3

    subplot(4,1,j), axis([0 100 -.15 .15])
    subplot(4,3,9+j), axis([-.15 .15 0 80]), set(gca,'Ytick',[0 40 80])
    h = findobj(gca,'Type','patch');
    h.FaceColor = [0.6 0.6 0.6];
    h.EdgeColor = 'k';
end

%%
% matrix overdetermined system

clear A
clear b
clear x
figure(3)

% overdetermined
n=300; m=60; p=20;
A=rand(n,m); b=rand(n,p);

lam=[0 0.1];
for j=1:2
    cvx_begin;
    variable x(m,p)
    minimize(norm(A*x-b,2) + lam(j)*norm(x,1));
    cvx_end;
    subplot(2,1,j), pcolor(x.'), colormap(hot), colorbar
end

%%
subplot(2,1,1), set(gca,'Fontsize',[15],'Xtick',[1 20 40 60],'Ytick',[1 10 20])
subplot(2,1,2), set(gca,'Fontsize',[15],'Xtick',[1 20 40 60],'Ytick',[1 10 20])

author: dynamicslab


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

% Last Modified by GUIDE v2.5 10-Jun-2012 21:47:30

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @PlayJAABAResults_OpeningFcn, ...
                   'gui_OutputFcn',  @PlayJAABAResults_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1}) && exist(varargin{1}) %#ok<EXIST>
    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 PlayJAABAResults is made visible.
function PlayJAABAResults_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 PlayJAABAResults (see VARARGIN)

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

persistent expdir;
persistent classifierparamsfile;

[handles.expdir,...
  handles.classifierparamsfile] = ...
  myparse(varargin,'expdir','',...
  'classifierparamsfile','');

if isempty(handles.expdir),
  handles.expdir = uigetdir(expdir,'Experiment directory');
  if ~ischar(handles.expdir),
    delete(handles.figure_PJR);
    return;
  end
  expdir = handles.expdir;
end
if isempty(handles.classifierparamsfile),
  if isempty(classifierparamsfile), classifierparamsfile = ''; end
  if iscell(classifierparamsfile),
    classifierparamsfile = classifierparamsfile{1};
  end
  [filename,path] = uigetfile(classifierparamsfile,'Classifier parameters file','MultiSelect','on');
  if ~ischar(filename) && ~iscell(filename),
    delete(handles.figure_PJR);
    return;
  end
  if iscell(filename),
    handles.classifierparamsfile = cellfun(@(x) fullfile(path,x),filename,'UniformOutput',false);
  else
    handles.classifierparamsfile = fullfile(path,filename);
  end
  classifierparamsfile = handles.classifierparamsfile;
end

handles = InitializeState(handles);

handles = InitializePlots(handles);

handles = SetCurrentFly(handles);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes PlayJAABAResults wait for user response (see UIRESUME)
% uiwait(handles.figure_PJR);

function handles = InitializeState(handles)

score_norm_prctile = 80;

%% read config parameters
handles.classifierparams = ReadClassifierParamsFile(handles.classifierparamsfile);
handles.nbehaviors = numel(handles.classifierparams);
handles.behaviors = cell(1,handles.nbehaviors);
handles.scorefns = cell(1,handles.nbehaviors);
for i = 1:handles.nbehaviors,
  if iscell(handles.classifierparams(i).behaviors.names),
    handles.behaviors{i} = sprintf('%s_',handles.classifierparams(i).behaviors.names{:});
    handles.behaviors{i} = handles.behaviors{i}(1:end-1);
  else
    handles.behaviors{i} = handles.classifierparams(i).behaviors.names;
  end
  handles.behaviors{i} = regexprep(handles.behaviors{i},'(_|^)([a-z])','${upper($2)}');
  scorefn = handles.classifierparams(i).file.scorefilename;
  handles.scorefns{i} = regexprep(scorefn,'\.mat$','');
end

%% create trx structure

%global STORED_TRX;
STORED_TRX = [];

if isempty(STORED_TRX) || ~strcmp(STORED_TRX.expdirs,handles.expdir),
handles.trx = Trx('trxfilestr',handles.classifierparams(1).file.trxfilename,...
  'moviefilestr',handles.classifierparams(1).file.moviefilename,...
  'perframedir',handles.classifierparams(1).file.perframedir);
if isfield(handles.classifierparams,'perframe') && isfield(handles.classifierparams(1).perframe,'params'),
  handles.trx.SetPerFrameParams(handles.classifierparams(1).perframe.params);
end
if isfield(handles.classifierparams,'perframe') && isfield(handles.classifierparams(1).perframe,'landmark_params'),
  handles.trx.SetLandmarkParams(handles.classifierparams(1).perframe.landmark_params);
end

handles.trx.AddExpDir(handles.expdir,'openmovie',false);
STORED_TRX = handles.trx;
else

  handles.trx = STORED_TRX;
end

%% for reading video

moviename = fullfile(handles.expdir,handles.trx.moviefilestr);
[handles.readframe,handles.nframes,handles.fid,handles.headerinfo] = get_readframe_fcn(moviename);

%% current fly, frame, etc.

handles.target = 1;
handles.t = handles.trx(handles.target).firstframe;
handles.nflies = handles.trx.nflies;

%% playback fps

handles.fps = 3;
handles.dozoom = true;

%% choose weights for each behavior

handles.behaviorweight = nan(handles.nbehaviors,1);
for behaviori = 1:handles.nbehaviors,
  tmp = handles.trx.(handles.scorefns{behaviori});
  if iscell(tmp),
    tmp = [tmp{:}];
  end
  handles.behaviorweight(behaviori) = nnz(tmp>0);
end
handles.behaviorweight = handles.behaviorweight / sum(handles.behaviorweight);
handles.behaviorweight = handles.behaviorweight(:);

%% choose colors for each behavior

tmp = jet(256)*.7;
handles.behavior_colors = tmp(round(linspace(1,256,handles.nbehaviors)),:);

%% initial zoom

a = handles.trx.a;
if iscell(a),
  a = [a{:}];
end
meana = nanmedian(a)*4;
awidthradius = 3;
aheightradius = 3;
aborder = 1;
aback = .25;
handles.pxwidthradius = awidthradius*meana;
handles.pxheightradius = aheightradius*meana;
handles.pxborder = aborder*meana;
handles.pxback = meana*aback;
handles.meana = meana;

%% get the score normalizations

handles.score_norms = nan(handles.nbehaviors,1);
for i = 1:handles.nbehaviors,
  tmp = handles.trx.(handles.scorefns{i});
  if iscell(tmp),
    tmp = [tmp{:}];
  end
  handles.score_norms(i) = prctile(abs(tmp),score_norm_prctile);
end

%% not playing

handles.isplaying = false;

function handles = InitializePlots(handles)

trx_linewidth = 1;
target_linewidth = 2;
trxcolor = 'k';
textcolor = 'k';
fontsize = 30;
bar_width = .8;
barbackcolor = [.5,.5,.5];

%% main axis
% frame
im = handles.readframe(handles.t);
ncolors = size(im,3);
if ncolors == 1,
  handles.him = imagesc(im2double(im),'Parent',handles.haxmain,[0,1]);
  colormap(handles.haxmain,'gray');
else
  handles.him = image(im2double(im),'Parent',handles.haxmain);
end
hold(handles.haxmain,'on');

% plot current fly
handles.htrx = plot(handles.haxmain,nan,nan,'.-','LineWidth',trx_linewidth,'color',trxcolor,'HitTest','off');
handles.htarget = plot(handles.haxmain,nan,nan,'-','LineWidth',target_linewidth,'color','k','HitTest','off');
set(handles.haxmain,'XTick',[],'YTick',[]);

% plot other flies
handles.hothers = nan(1,handles.nflies);
for i = 1:handles.nflies,
  handles.hothers(i) = plot(handles.haxmain,nan,nan,'mo','markerfacecolor','m');
  set(handles.hothers(i),'ButtonDownFcn',@(hObject,eventdata) PlayJAABAResults('fly_ButtonDownFcn',hObject,eventdata,guidata(hObject),i));
end

% zoom
set(handles.haxmain,'Units','pixels');
handles.mainpos = get(handles.haxmain,'Position');
set(handles.haxmain,'Units','normalized');
pxheightradius1 = min((handles.pxwidthradius)*handles.mainpos(4) / handles.mainpos(3),(handles.headerinfo.nr-1)/2);
pxwidthradius1 = min((handles.pxheightradius)*handles.mainpos(3) / handles.mainpos(4),(handles.headerinfo.nc-1)/2);
handles.pxheightradius = max(pxheightradius1,handles.pxheightradius);
handles.pxwidthradius = max(pxwidthradius1,handles.pxwidthradius);
axis(handles.haxmain,'image');
%axis(handles.haxmain,handles.ax);

% target info
handles.hinfo = text(0,0,'','Parent',handles.haxmain,'Color',textcolor,...
  'FontUnits','pixels','FontSize',fontsize,...
  'HorizontalAlignment','left',...
  'VerticalAlignment','bottom',...
  'HitTest','off');

% current behaviors
handles.htext = text(nan,nan,'','Parent',handles.haxmain,'Color',textcolor,...
  'FontUnits','pixels','FontSize',fontsize);

%% bar plots
handles.hbar = nan(1,handles.nbehaviors);
handles.hbarback = nan(1,handles.nbehaviors);
for i = 1:handles.nbehaviors,
  handles.hbarback(i) = patch(i+bar_width/2*[-1,1,1,-1,-1],zeros(1,5)+.1,barbackcolor,'EdgeColor','none',...
    'Parent',handles.haxbottom);
  if i == 1,
    hold(handles.haxbottom,'on');
  end
  handles.hbar(i) = patch(i+bar_width/2*[-1,1,1,-1,-1],zeros(1,5)+.1,handles.behavior_colors(i,:),'EdgeColor','none','Parent',handles.haxbottom);
end

% zero line
plot(handles.haxbottom,[0,handles.nbehaviors+1],[0,0],'--','color',[.5,.5,.5]);
set(handles.haxbottom,'XTick',1:handles.nbehaviors,'XTickLabel',handles.behaviors,'YTick',[],...
  'XLim',[0,handles.nbehaviors+1],'YLim',[-1,1],'XColor','w','Color','k','TickLength',[0,0]);

% slider
set(handles.slider_frame,'Min',1,'Max',handles.nframes,...
  'Value',1,...
  'SliderStep',[1/(handles.nframes-1),100/(handles.nframes-1)]);

RecursiveSetKeyPressFcn(handles.figure_PJR);

%% set current fly

function handles = SetCurrentFly(handles)

handles.x = handles.trx(handles.target).x;
handles.y = handles.trx(handles.target).y;
handles.a = handles.trx(handles.target).a;
handles.b = handles.trx(handles.target).b;
handles.theta = handles.trx(handles.target).theta;
fil = [0.000263865082737   0.106450771973592   0.786570725887342   0.106450771973592   0.000263865082737];
fil = fil / sum(fil(:));
handles.smooththeta = imfilter(handles.theta,fil,'same','replicate');
handles.nframes_curr = handles.trx(handles.target).nframes;
handles.firstframe = handles.trx(handles.target).firstframe;
handles.endframe = handles.trx(handles.target).endframe;
handles.scores = nan(handles.nbehaviors,handles.nframes_curr);
handles.sex = handles.trx(handles.target).sex;
for behaviori = 1:handles.nbehaviors,
  handles.scores(behaviori,:) = min(1,max(-1,handles.trx(handles.target).(handles.scorefns{behaviori})/handles.score_norms(behaviori)));
end

if handles.t < handles.firstframe,
  handles.t = handles.firstframe;
end
if handles.t > handles.endframe,
  handles.t = handles.endframe;
end

SetCurrentFrame(handles);

%% set current frame

function SetCurrentFrame(handles,hObject)

if nargin < 2,
  hObject = nan;
end
if hObject ~= handles.slider_frame,
  set(handles.slider_frame,'Value',handles.t);
end
if hObject ~= handles.edit_frame,
  set(handles.edit_frame,'String',handles.t);
end
UpdatePlots(handles);

%% update plots

function UpdatePlots(handles,resetax)

persistent ax;

if nargin > 1 && resetax,
  ax = [];
end

trx_rad = 50;
printnone = true;

im = handles.readframe(handles.t);
set(handles.him,'CData',im2double(im));
i = handles.t+handles.trx(handles.target).off;

i0 = max(1,min(handles.nframes_curr,i-trx_rad(1)));
i1 = max(1,min(handles.nframes_curr,i+trx_rad(1)));

set(handles.htrx,'XData',handles.x(i0:i1),...
  'YData',handles.y(i0:i1));

updatefly(handles.htarget,handles.x(i),handles.y(i),...
  handles.theta(i),handles.a(i),handles.b(i));

for behaviori = 1:handles.nbehaviors,
  set(handles.hbar(behaviori),'YData',[0,0,1,1,0]*handles.scores(behaviori,i));
  set(handles.hbarback(behaviori),'YData',[0,0,1,1,0]*sign(handles.scores(behaviori,i)));
end
w = max(0,handles.scores(:,i)).*handles.behaviorweight;
z = sum(w);
if z == 0,
  colortmp = [0,0,0];
else
  colortmp = sum(bsxfun(@times,w,handles.behavior_colors),1)/z*max(handles.scores(:,i));
end
set(handles.htarget,'Color',colortmp);

tailpos = [handles.x(i)-(handles.pxback+2*handles.a(i))*cos(handles.theta(i)),...
  handles.y(i)-(handles.pxback+2*handles.a(i))*sin(handles.theta(i))];
sectori = floor(modrange((handles.smooththeta(i)-pi/8),0,2*pi)/pi*4)+1;

if ismember(sectori,[1,7,8]),
  halign = 'right';
elseif ismember(sectori,[2,6]),
  halign = 'center';
elseif ismember(sectori,[3,4,5]),
  halign = 'left';
else
  error('Sanity check: sectori > 8 or < 1');
end
if ismember(sectori,[1,2,3]),
  valign = 'bottom';
elseif ismember(sectori,[4,8]),
  valign = 'middle';
elseif ismember(sectori,[5,6,7]),
  valign = 'top';
else
  error('Sanity check: sectori > 8 or < 1');
end

labelis = find(handles.scores(:,i) > 0);
s = cell(1,numel(labelis));
for iii = 1:numel(labelis),
  s{iii} = sprintf('\\color[rgb]{%f %f %f}%s',handles.behavior_colors(labelis(iii),:),handles.behaviors{labelis(iii)});
end
if printnone && isempty(s),
  s = {'None'};
end

set(handles.htext,'String',s,...
  'Position',tailpos,...
  'HorizontalAlignment',halign,...
  'VerticalAlignment',valign);

set(handles.hinfo,'String',sprintf('Target %d, sex = %s',handles.target,handles.sex{i}));

for fly = 1:handles.nflies,
  if fly == handles.target || ...
      handles.t > handles.trx(fly).endframe || ...
      handles.t < handles.trx(fly).firstframe,
    set(handles.hothers(fly),'XData',nan,'YData',nan);
  else
    i = handles.t + handles.trx(fly).off;
    set(handles.hothers(fly),'XData',handles.trx(fly).x(i),...
      'yData',handles.trx(fly).y(i));
  end
end

outofbounds = handles.dozoom && (isempty(ax) || ...
  handles.x(i)-handles.pxborder < ax(1) || ...
  handles.x(i)+handles.pxborder > ax(2) || ...
  handles.y(i)-handles.pxborder < ax(3) || ...
  handles.y(i)+handles.pxborder > ax(4));
if outofbounds,
  ax = ResetAxis(handles.x,handles.y,i,handles.headerinfo.nr,handles.headerinfo.nc,handles.pxwidthradius,handles.pxheightradius,handles.pxborder);
  axis(handles.haxmain,ax);
  set(handles.hinfo,'Position',ax([1,4]));
end

%%

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

% --- Executes on slider movement.
function slider_frame_Callback(hObject, eventdata, handles)
% hObject    handle to slider_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,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
v = round(get(hObject,'Value'));
if v > handles.endframe,
  v = handles.endframe;
  set(hObject,'Value',v);
elseif v < handles.firstframe,
  v = handles.firstframe;
  set(hObject,'Value',v);
end
handles.t = v;
SetCurrentFrame(handles,hObject);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function slider_frame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider_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: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
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
v = str2double(get(hObject,'String'));
if isnan(v),
  set(hObject,'Value',num2str(handles.t));
  return;
end
v = max(handles.firstframe,v);
v = min(handles.endframe,v);
v = round(v);
handles.t = v;
SetCurrentFrame(handles,hObject);
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 handles = stop(handles)

if handles.isplaying,
  pushbutton_play_Callback(handles.pushbutton_play,[],handles);
  handles = guidata(handles.pushbutton_play);
end

function playstop(handles)

pushbutton_play_Callback(handles.pushbutton_play,[],handles);

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

if handles.isplaying,
  handles.isplaying = false;
  guidata(hObject,handles);
else
  handles.isplaying = true;
  guidata(hObject,handles);
  play(handles);
end

function play(handles)

set(handles.pushbutton_play,'BackgroundColor',[.5,0,0],'String','Stop');

tic;
dt_sec = 1/handles.fps;

while true,

  if ~handles.isplaying,
    break;
  end

  if handles.t == handles.endframe,
    break;
  end

  handles = guidata(handles.pushbutton_play);
  dt = max(1,round(handles.fps*dt_sec));
  %fprintf('dt_sec = %f, dt = %d\n',dt_sec,dt);
  handles.t = handles.t+dt;
  guidata(handles.pushbutton_play,handles);
  SetCurrentFrame(handles);
  drawnow;
  dt_sec = toc;
  tic;
  if dt == 1,
    pause_time = 1/handles.fps - dt_sec;
    if pause_time > 0,
      %fprintf('pausing for %f s\n',pause_time);
      pause(pause_time);
    end
  end

end

set(handles.pushbutton_play,'BackgroundColor',[0,.5,0],'String','Play');
handles = guidata(handles.pushbutton_play);
handles.isplaying = false;
guidata(handles.pushbutton_play,handles);

function ax = ResetAxis(x,y,i,nr,nc,pxwidthradius,pxheightradius,pxborder)

n = numel(x);

minx = x(i);
maxx = x(i);
miny = y(i);
maxy = y(i);
w = pxwidthradius-2*pxborder;
h = pxheightradius-2*pxborder;
for j = i+1:n,
  minx = min(minx,x(j));
  maxx = max(maxx,x(j));
  dx = ceil(maxx-minx)+1;
  if dx > w,
    break;
  end
  miny = min(miny,y(j));
  maxy = max(maxy,y(j));
  dy = ceil(maxy-miny)+1;
  if dy > h,
    break;
  end
end
mux = (minx+maxx)/2;
muy = (miny+maxy)/2;
ax = [mux+pxwidthradius*[-1,1],muy+pxheightradius*[-1,1]];

if ax(2)>nc && ax(1)<1,
  ax(1:2) = (nc+1)/2+pxwidthradius*[-1,1];
elseif ax(2)>nc,
  ax(2) = nc;
  ax(1) = nc-(2*pxwidthradius+1)+1;
elseif ax(1)<1,
  ax(1) = 1;
  ax(2) = 2*pxwidthradius+1;
end
if ax(2)>nc || ax(1)<1,
  ax(1:2) = (nc+1)/2+pxwidthradius*[-1,1];
end
if ax(4)>nr && ax(3)<1,
  ax(3:4) = (nr+1)/2+pxwidthradius*[-1,1];
elseif ax(4)>nr,
  ax(4) = nr;
  ax(3) = nr-(2*pxheightradius+1)+1;
elseif ax(3)<1,
  ax(3) = 1;
  ax(4) = 2*pxheightradius+1;
end
if ax(4)>nr || ax(3)<1,
  ax(3:4) = (nr+1)/2+pxheightradius*[-1,1];
end

% --- Executes when figure_PJR is resized.
function figure_PJR_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to figure_PJR (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isfield(handles,'pxwidthradius'),
  return;
end

awidthradius = 3;
aheightradius = 3;
aborder = 1;
aback = .25;
meana = handles.meana;
handles.pxwidthradius = awidthradius*meana;
handles.pxheightradius = aheightradius*meana;
handles.pxborder = aborder*meana;
handles.pxback = meana*aback;

set(handles.haxmain,'Units','pixels');
handles.mainpos = get(handles.haxmain,'Position');
set(handles.haxmain,'Units','normalized');
pxheightradius1 = min((handles.pxwidthradius)*handles.mainpos(4) / handles.mainpos(3),(handles.headerinfo.nr-1)/2);
pxwidthradius1 = min((handles.pxheightradius)*handles.mainpos(3) / handles.mainpos(4),(handles.headerinfo.nc-1)/2);
handles.pxheightradius = max(pxheightradius1,handles.pxheightradius);
handles.pxwidthradius = max(pxwidthradius1,handles.pxwidthradius);
UpdatePlots(handles,true);
guidata(hObject,handles);

function fly_ButtonDownFcn(hObject,eventdata,handles,fly)

if fly == handles.target,
  return;
end

if handles.isplaying,
  handles = stop(handles);
end
handles.target = fly;
handles = SetCurrentFly(handles);
guidata(hObject,handles);

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

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

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

handles.dozoom = strcmpi(get(hObject,'Checked'),'on');
guidata(hObject,handles);
UpdateZoom(handles);

function UpdateZoom(handles)

if handles.dozoom,
  set(handles.menu_zoom_showwholevideo,'Checked','off');
  set(handles.menu_zoom_keepflyinaxes,'Checked','on');
  UpdatePlots(handles,true);
else
  set(handles.menu_zoom_showwholevideo,'Checked','on');
  set(handles.menu_zoom_keepflyinaxes,'Checked','off');
  axis(handles.haxmain,'image');
  set(handles.hinfo,'Position',[1,handles.headerinfo.nr]);
end

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

handles.dozoom = strcmpi(get(hObject,'Checked'),'off');
guidata(hObject,handles);
UpdateZoom(handles);

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

CleanUp(handles);

function CleanUp(handles)

if isfield(handles,'fid'),
  fclose(handles.fid);
end

% --- Executes when user attempts to close figure_PJR.
function figure_PJR_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure_PJR (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
delete(hObject);

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

while true,

  v = inputdlg('Playback speed (fps):','Playback speed',1,{num2str(handles.fps)});
  v = v{1};
  if ~ischar(v),
    return;
  end
  v = str2double(v);
  if isnan(v) || v <= 0,
    continue;
  end
  break;
end
handles.fps = v;
guidata(hObject,handles);

% --- Executes on key press with focus on figure_PJR and none of its controls.
function figure_PJR_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to figure_PJR (see GCBO)
% eventdata  structure with the following fields (see FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)

switch eventdata.Key,

  case 'leftarrow',
    if ~handles.isplaying,
      handles.t = max(handles.t-1,handles.firstframe);
      guidata(hObject,handles);
      SetCurrentFrame(handles);
    end
  case 'rightarrow',
    if ~handles.isplaying,
      handles.t = min(handles.t+1,handles.endframe);
      guidata(hObject,handles);
      SetCurrentFrame(handles);
    end
  case 'uparrow',
    if ~handles.isplaying,
      handles.t = max(handles.t-30,handles.firstframe);
      guidata(hObject,handles);
      SetCurrentFrame(handles);
    end
  case 'downarrow',
    if ~handles.isplaying,
      handles.t = min(handles.t+30,handles.endframe);
      guidata(hObject,handles);
      SetCurrentFrame(handles);
    end
  case {'space','p'},
    playstop(handles);
end

function RecursiveSetKeyPressFcn(hfig)

hchil = findall(hfig,'-property','KeyPressFcn');
goodidx = true(1,numel(hchil));
for i = 1:numel(hchil),
  if strcmpi(get(hchil(i),'Type'),'uicontrol') && strcmpi(get(hchil(i),'Style'),'edit'),
    goodidx(i) = false;
  end
end
set(hchil(goodidx),'KeyPressFcn',get(hfig,'KeyPressFcn'));

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

s = cellstr(num2str((1:handles.nflies)'));
[sel,ok] = listdlg('ListString',s,'SelectionMode','single',...
  'InitialValue',handles.target,'Name','Switch to target',...
  'PromptString','Switch to target: ');
if ~ok,
  return;
end
handles.target = sel;
handles = SetCurrentFly(handles);
guidata(hObject,handles);

author: kristinbranson

repo: JAABA


text:
function k = findrows(A, b)
%FINDROWS Find indices of a given row within a matrix.
%
%   FINDROWS(A, B) returns a column vector with the indices of the rows
%   in the matrix A that are identical to the row vector B.  If no rows
%   in A are identical to B, an empty vector is returned.
%
%   The methods uses a for-loop, but it uses less memory and is in many
%   cases a lot faster than the vectorized methods
%
%      find( all( A == repmat(b, size(A, 1), 1), 2 ) )
%      find( all( A == b(ones(size(A, 1), 1),:), 2 ) )
%
%   See also FIND, FINDCOLS.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:51:15 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   k = find( A(:,1) == b(1) );
   for j = 2:size(A, 2)
      k = k( A(k,j) == b(j) );
      if isempty(k)
         return
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function r4lib_test023 ( )

%*****************************************************************************80
%
%% R4LIB_TEST023 tests R4_SIGN and R4_SIGN3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2014
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;

  r4_test = [ -1.25, -0.25, 0.0, +0.5, +9.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4LIB_TEST023\n' );
  fprintf ( 1, '  R4_SIGN returns the sign of a number.\n' );
  fprintf ( 1, '  R4_SIGN3 returns the three-way sign of a number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    R4  R4_SIGN(R4)  R4_SIGN3(R4)\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    r4 = r4_test(test);
    s1 = r4_sign ( r4 );
    s2 = r4_sign3 ( r4 );
    fprintf ( 1, '  %8f  %8f  %8f\n', r4, s1, s2 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bndinfo2, err] = transferSuperpixelLabels(bndinfo, wseg2)

wseg1 = bndinfo.wseg;

labim = zeros(size(wseg1));
labim(wseg1>0) = bndinfo.labels(wseg1(wseg1>0));

stats = regionprops(wseg2, 'PixelIdxList');
spind = {stats(:).PixelIdxList};

lab2 = zeros(max(wseg2(:)), 1);

for k = 1:numel(lab2)
    pixlab = labim(spind{k});
    pixlab = pixlab(pixlab>0);
    if numel(pixlab)>numel(spind{k})*0.01
        lab2(k) = mode(pixlab);
    end
end

bndinfo2 = bndinfo;
bndinfo2.wseg = uint16(wseg2);
bndinfo2.labels = lab2;

labim2 = zeros(size(wseg2));
labim2(wseg2>0) = bndinfo2.labels(wseg2(wseg2>0));

ind = (labim2 > 0) & (labim > 0);

if nargout > 1
    err = mean(labim(ind)~=labim2(ind));
    disp(num2str([max(wseg1(:)) max(wseg2(:))]))
    disp(['Error: ' num2str(err)]);
end

author: Cloud-CV

repo: object-proposals


text:
function edgeBoxesSweeps()
% Parameter sweeps for Edges Boxes object proposals.
%
% Running the parameter sweeps requires altering internal flags.
% The sweeps are not well documented, use at your own discretion.
%
% Structured Edge Detection Toolbox      Version 3.0
% Copyright 2014 P. Dollar and L. Zitnick.  [pdollar-at-microsoft.com]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the MSR-LA Full Rights License [see license.txt]

% define parameter sweeps
rt = 'D:\code\research\edges\';
expNms = {'alpha','beta','minScore','edgeMinMag','edgeMergeThr',...
  'clusterMinMag','maxAspectRatio','minBoxArea','gamma','kappa'};
expNms=expNms(1:end); opts=createExp(rt,expNms); maxn=inf;

% run training and testing jobs
jobs = createJobs(rt,opts,maxn);
tic, for i=1:length(jobs), edgeBoxes(jobs{i}{:}); end; toc

% plot all results
for e=1:length(expNms)
  eval={'data',boxesData('split','val'),'names',{opts{e}.name},...
    'resDir',[rt 'boxes/sweeps/'],'maxn',maxn,'fName',expNms{e}};
  boxesEval(eval{:},'thrs',.7);
  boxesEval(eval{:},'thrs',.5:.05:1,'cnts',1000);
end

end

function jobs = createJobs( rt, opts, maxn )
% create jobs
M='models/forest/modelBsds'; M=load(M); M=M.model; M.opts.nThreads=1;
data=boxesData('split','val'); fs=data.imgs(1:min(end,maxn));
opts=[opts{:}]; N=length(opts); jobs=cell(1,N); D=zeros(1,N);
for e=1:N, opts(e).name=[rt 'boxes/sweeps/' opts(e).name '-val.mat']; end
for e=1:N, D(e)=exist(opts(e).name,'file')==2; jobs{e}={fs,M,opts(e)}; end
[~,K]=unique({opts.name},'stable'); D=D(K); jobs=jobs(K); jobs=jobs(~D);
fprintf('nJobs = %i\n',length(jobs));
end

function opts = createExp( rt, expNm )

% if expNm is a cell, call recursively and return
if( iscell(expNm) )
  N=length(expNm); opts=cell(1,N);
  for e=1:N, opts{e}=createExp(rt,expNm{e}); end; return;
end

% setup opts
opts=edgeBoxes(); opts.minScore=0;
N=100; optsDefault=opts; opts=opts(ones(1,N));
switch expNm
  case 'alpha'
    vs=45:5:75; N=length(vs);
    for e=1:N, opts(e).alpha=vs(e)/100; end
  case 'beta'
    vs=60:5:90; N=length(vs);
    for e=1:N, opts(e).beta=vs(e)/100; end
  case 'minScore'
    vs=[0 5 10 25 50 100]; N=length(vs);
    for e=1:N, opts(e).minScore=vs(e)/1000; end
  case 'edgeMinMag'
    vs=[0 50 100 200 400]; N=length(vs);
    for e=1:N, opts(e).edgeMinMag=vs(e)/1000; end
  case 'edgeMergeThr'
    vs=[25 50 100 200 400]; N=length(vs);
    for e=1:N, opts(e).edgeMergeThr=vs(e)/100; end
  case 'clusterMinMag'
    vs=[0 50 100 200 400]; N=length(vs);
    for e=1:N, opts(e).clusterMinMag=vs(e)/100; end
  case 'maxAspectRatio'
    vs=1:5; N=length(vs);
    for e=1:N, opts(e).maxAspectRatio=vs(e); end
  case 'minBoxArea'
    vs=[100 250 500 1000 2500 5000]; N=length(vs);
    for e=1:N, opts(e).minBoxArea=vs(e); end
  case 'gamma'
    vs=[25 50 100 200 400 10000]; N=length(vs);
    for e=1:N, opts(e).gamma=vs(e)/100; end
  case 'kappa'
    vs=50:25:200; N=length(vs);
    for e=1:N, opts(e).kappa=vs(e)/100; end
  otherwise, error('invalid exp: %s',expNm);
end

% produce final set of opts and find default opts
O=1:N; opts=opts(O); d=0;
for e=1:N, if(isequal(optsDefault,opts(e))), d=e; break; end; end
if(d==0), disp(expNm); assert(false); end; opts(d).name='Default';
for e=1:N, if(e~=d), opts(e).name=[expNm int2str2(vs(e),5)]; end; end

end


text:
function [f, df, d2f] = nlp_costfcn(om, x)
%NLP_COSTFCN  Evaluates objective function, gradient and Hessian.
%   [F, DF, D2F] = NLP_COSTFCN(OM, X)
%
%   Objective function evaluation routine, suitable for use with MIPS,
%   FMINCON, etc. Computes objective function value, gradient and Hessian.
%
%   Inputs:
%     OM : Opt-Model object
%     X : optimization vector
%
%   Outputs:
%     F   : value of objective function
%     DF  : (optional) gradient of objective function (column vector)
%     D2F : (optional) Hessian of objective function (sparse matrix)
%
%   Examples:
%       f = nlp_costfcn(om, x);
%       [f, df] = nlp_costfcn(om, x);
%       [f, df, d2f] = nlp_costfcn(om, x);
%
%   See also NLP_CONSFCN, NLP_HESSFCN.

%   MP-Opt-Model
%   Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%%----- evaluate objective function -----
%% general nonlinear costs
if nargout == 3
    [f, df, d2f]    = om.eval_nln_cost(x);
    if om.qdc.NS
        [fq, dfq, d2fq] = om.eval_quad_cost(x);
        f = f + sum(fq);
        df = df + dfq;
        d2f = d2f + d2fq;
    end
elseif nargout == 2
    [f, df]   = om.eval_nln_cost(x);
    if om.qdc.NS
        [fq, dfq] = om.eval_quad_cost(x);
        f = f + sum(fq);
        df = df + dfq;
    end
else
    f  = om.eval_nln_cost(x);
    if om.qdc.NS
        fq = om.eval_quad_cost(x);
        f = f + sum(fq);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% OPENMAS HIT-BOX ENUMERATION CLASSES (OMAS_hitBoxType.m) %%%%%%%%%%%%%%%%
% This file contains the enumeration definition of the simulation
% recognised object types.

% Author: James A. Douthwaite 25/01/19

classdef OMAS_hitBoxType < uint8
   % Assign uint8 type to enumeration values
   enumeration
      none(0);           % Generic
      spherical(1);      % Spherical bounding box assumption
      AABB(2);           % Axis aligned bounding box
      OBB(3);            % Object orientated bounding box
      capsule(4);        % Capsule enclosed volume
      mesh(5);           % Mesh enclosed volume
   end
end

repo: OpenMAS


text:
%% Reading and writing midi files
% This short demo will show you how to read in a midi file, do some simple
% manipulation, and then write out new modified files.
%
% Before you try these function, you will need to make sure that
% KaraokeMidiJava.jar can be found by matlab.  The easiest way to do this
% is to edit Matlab's |classpath.txt| and add the path to the midi jar.
% You will, I believe, need to restart Matlab once you have updated the
% classpath file.
%
% |readmidi_java.m| and |writemidi_java.m| were inspired by the |readmidi.m| and
% |writemidi.m| of the
% <https://www.jyu.fi/hum/laitokset/musiikki/en/research/coe/materials/miditoolbox/ Midi Toolbox>.
% My versions can take the same inputs, but also
% have a few extensions, which you can find by reading the help information
% for each function.
%

%% Reading in a midi file
% First read in the file, which is "The sun, whose rays are all ablaze"
% from _The Mikado_ by Gilbert and Sullivan.  (Thanks to the
% <http://math.boisestate.edu/gas/mikado/html/index.html Gilbert & Sullivan archives> for this file!)

% Read a midi file into a note matrix, getting the optional
% track column.  (Note that .kar files are midi karaoke files.  They use
% exactly the same file format as .mid files.  The .kar basically just
% indicates that this is a midi file which contains midi lyric messages.)
%
% The columns are:
%       (1) - note start in beats
%       (2) - note duration in beats
%       (3) - channel
%       (4) - midi pitch (60 --> C4 = middle C)
%       (5) - velocity
%       (6) - note start in seconds
%       (7) - note duration in seconds
note_matrix = readmidi_java('202.kar');

% see how many notes we have
fprintf('Number of notes: %d\n', size(note_matrix,1));
% plot the note onsets
figure(1);
clf;
plot(note_matrix(:,1));
ylabel('Note onset (beats)');
xlabel('Note index');
title('Note onsets in beats');
figure(2);
clf;
plot(note_matrix(:,6));
ylabel('Note onset (seconds)');
xlabel('Note index');
title('Note onsets in seconds');

%% Writing out a midi file
% Now write out a couple of midi files that just contain the notes from the
% vocals melody.

% The vocal melody is in channel 3.
vocal_matrix = note_matrix(note_matrix(:,3) == 3,:);
% Get the tempo information from 202.kar
tempos = get_tempos('202.kar');
fprintf('Number of tempos: %d\n', size(tempos,1));

% Write out a midi file that preserves the beat information.  Use the last
% tempo from 202.kar
writemidi_java(vocal_matrix,'just_vocals_beats.mid',120,tempos(end,1));
% Write out another version that preserves the timing.  This may be helpful
% if you want to synthesize the file.  202.kar has
% multiple tempo changes, so specifying a single tempo in writemidi_java
% won't help.  This function writes a midi file that preserves the times in
% the note matrix (columns 6 & 7), but not the beats.
writemidi_seconds(vocal_matrix,'just_vocals_seconds.mid');

%%  Checking written midi file
% Read the vocal file in and make sure that it matches what we wrote.
% The |just_vocals_beats.mid| file should have preserved the beat information
% in |note_matrix|, but the timing information will be wrong because of the
% tempo changes in |202.kar|.  The |just_vocals_seconds.mid| file should do the
% opposite.

vm_beats = readmidi_java('just_vocals_beats.mid');
vm_seconds = readmidi_java('just_vocals_seconds.mid');

figure(1);
clf
plot(vocal_matrix(:,1), 'b.-');
hold on;
plot(vm_beats(:,1),'r');
plot(vm_seconds(:,1),'g');
legend('manipulated note matrix', 'correct beats matrix',...
    'correct times matrix','Location','NorthWest');
legend('boxoff');
ylabel('Note onset (beats)');
xlabel('Note index');
title('Vocal line notes in beats');

figure(2);
clf
plot(vocal_matrix(:,6), 'b.-');
hold on;
plot(vm_beats(:,6),'r');
plot(vm_seconds(:,6),'g');
legend('manipulated note matrix', 'correct beats matrix',...
    'correct times matrix','Location','NorthWest');
legend('boxoff');
ylabel('Note onset (seconds)');
xlabel('Note index');
title('Vocal line notes in seconds');

% 2010-05-03 Christine Smit csmit@ee.columbia.edu

author: Sable

repo: mcbench-benchmarks


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtFfmMaxwellCFIEpec.m                           |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Francois Alouges            |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Solve PEC scatering problem with CFIE         |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N    = 1e3
tol  = 1e-3
typ  = 'RWG'
gss  = 3
beta = 0.5;

% Spherical mesh
sphere = mshSphere(N,1);
sigma  = dom(sphere,gss);
figure
plot(sphere)
axis equal

% Frequency adjusted to maximum esge size
stp = sphere.stp;
k   = 1/stp(2);
c   = 299792458;
f   = (k*c)/(2*pi);
disp(['Frequency : ',num2str(f/1e6),' MHz']);

% Incident direction and field
X0 = [0 0 -1];
E  = [0 1  0]; % Polarization (+x for Theta-Theta and +y for Phi-Phi)
H  = cross(X0,E);

% Incident Plane wave (electromagnetic field)
PWE{1} = @(X) exp(1i*k*X*X0') * E(1);
PWE{2} = @(X) exp(1i*k*X*X0') * E(2);
PWE{3} = @(X) exp(1i*k*X*X0') * E(3);

PWH{1} = @(X) exp(1i*k*X*X0') * H(1);
PWH{2} = @(X) exp(1i*k*X*X0') * H(2);
PWH{3} = @(X) exp(1i*k*X*X0') * H(3);

% Incident wave representation
plot(sphere,real(PWH{1}(sphere.vtx)))
title('Incident wave')
xlabel('X');   ylabel('Y');   zlabel('Z');
hold off
view(0,10)

%%% PREPARE OPERATOR
disp('~~~~~~~~~~~~~ PREPARE OPERATOR ~~~~~~~~~~~~~')

% Green kernel function --> G(x,y) = exp(ik|x-y|)/|x-y|
Gxy = '[exp(ikr)/r]';
Hxy = {'grady[exp(ikr)/r]1','grady[exp(ikr)/r]2','grady[exp(ikr)/r]3'};

% Finite elements
u = fem(sphere,'RWG');
v = fem(sphere,'RWG');

% Finite element mass matrix --> \int_Sx psi(x)' psi(x) dx
Id = integral(sigma,u,v);

% Finite element boundary operator
tic
T = 1i*k/(4*pi)*integral(sigma, sigma, v, Gxy, k, u, tol) ...
    - 1i/(4*pi*k)*integral(sigma, sigma, div(v), Gxy, k, div(u), tol) ;
toc

% Regularization
tic
Tr = 1i*k/(4*pi)*regularize(sigma, sigma, v, '[1/r]', u) ...
      - 1i/(4*pi*k)*regularize(sigma, sigma, div(v), '[1/r]', div(u));
T  = T + Tr;
toc

% Finite element boundary operator
tic
nxK = 1/(4*pi) * integral(sigma, sigma, nx(v), Hxy, k, u, tol);
toc

% Regularization
tic
nxKr = 1/(4*pi) * regularize(sigma, sigma, nx(v), 'grady[1/r]', u);
nxK  = nxK + nxKr;
toc

% Left hand side
LHS  = - beta * T  + (1-beta) * (0.5*Id - nxK);
LHSr = - beta * Tr + (1-beta) * (0.5*Id - nxKr);

% Right hand side
RHS = beta*integral(sigma,v,PWE) - (1-beta)*integral(sigma,nx(v),PWH);

%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% ILU preconditionner
tic
[L,U] = ilu(LHSr);
toc

% Iterative solver
tic
J = mgcr(@(V) LHS*V,RHS,[],tol,1000,L,U);
toc

%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
Ninf  = 1e3;
theta = 2*pi/1e3 .* (1:Ninf)';
nu    = [sin(theta),zeros(size(theta)),cos(theta)];

% Green kernel function
Ginf = '[exp(-ikxy)]';

% Finite element infinite operator --> \int_Sy exp(ik*nu.y) * psi(y) dx
Tinf = integral(nu,sigma,Ginf,k,v,tol);
sol  = 1i*k/(4*pi)*cross(nu, cross([Tinf{1}*J, Tinf{2}*J, Tinf{3}*J], nu));

% Radiation infinie de reference, convention e^(+ikr)/r
nMax = 100; refInf = zeros(Ninf,1);
if E(1) == 1
    for jj = 1:Ninf
        refInf(jj,:) = sphereMaxwell(1, -f, theta(jj), 0.0, nMax);
    end
else
    for jj = 1:Ninf
        [~,refInf(jj)] = sphereMaxwell(1, -f, theta(jj), pi/2, nMax);
    end
end
refInf = refInf ./ sqrt(4*pi);

% Radiations infinies en X
if E(1) == 1
    sol = sin(theta)'.*sol(:,3) - cos(theta)'.*sol(:,1);
else
    sol = sol(:,2);
end

% Erreur
eL2   = norm(refInf-sol,2)/norm(refInf,2)
eLINF = norm(refInf-sol,'inf')/norm(refInf,'inf')

% Representation graphique
figure
subplot(1,2,1)
plot(theta,20*log10(abs(sol)),'b',theta,20*log10(abs(refInf)),'r--')

subplot(1,2,2)
plot(theta,real(sol),'--b', theta,imag(sol),'--r', theta, real(refInf),':b', theta,imag(refInf),':r');
drawnow

%%% SURFACIC RADIATION
disp('~~~~~~~~~~~~~ SURFACIC RADIATION ~~~~~~~~~~~~~')

% Mesh Interpolation
Jmsh = feval(u,J,sphere);
V    = sqrt(sum(real(cell2mat(Jmsh)).^2,2));

% Graphical representation
figure
plot(sphere)
hold on
plot(sphere,V)
axis equal
title('|J| surfacic')
xlabel('X');   ylabel('Y');   zlabel('Z');
colorbar

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y = PhaseMod(x, bits)
% function y = PhaseMod(x, bits)
%
% This function performs phase modulation using l bits per symbol.  Thus the
% modulation is (2^bits)-PSK.  The return vector y is complex baseband.

l=bits;

M = 2^l;

if size(x,1)==1
   x = x.';
end

a = size(x,1);

x = num2str(x);

s = reshape(x,a/l,l);

ss = bin2dec(s);

y = exp(j*2*pi*ss/M);


text:
function varargout = calcTensor(odf,T,varargin)
% compute the average tensor for an ODF
%
% Syntax
%   [TVoigt, TReuss, THill] = calcTensor(odf,T)
%   THill = calcTensor(odf,T,'Hill')
%   TGeo = calcTensor(odf,T,'geometric')
%
% Input
%  odf - @SO3Fun
%  T   - @tensor
%
% Output
%  T    - @tensor
%
% Options
%  Voigt     - Boigt mean
%  Reuss     - Reuss mean
%  Hill      - Hill mean
%  geometric - geometric mean
%
% See also
% tensor/mean EBSD/calcTensor

% decide between the quadrature based method and the harmonic method
if ~check_option(varargin,'quadrature')

  % the harmonic route is directly implemented into tensor/mean
  [varargout{1:nargout}] = mean(T,odf,varargin{:});

else % quadrature based method

  % define a grid
  res = get_option(varargin,'resolution',2.5*degree);
  S3G = equispacedSO3Grid(odf.CS,odf.SS,'resolution',res);

  % evaluate the ODF
  f = eval(odf,S3G,varargin{:});
  f = f ./ sum(f(:));

  % compute the means
  [varargout{1:nargout}] = mean(S3G*T,'weights',f(:),varargin{:});

end

end

author: mtex-toolbox

repo: mtex

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


text:
%% Testing SCIP Option Setting Methods
% 17/5/15
clc
clear
% Default Problem
a = 100;
fun = @(x) a*(x(2) - x(1)^2)^2 + (1 - x(1))^2;
lb = [-1.5; -1.5]; ub = [2;2]; x0 = [-2; 1];
%Convert to Ins List
zs = zeros(size(x0));
x = scipvar(size(x0));
fs = fun(x); nl.obj_instr = fs.ins; nl.obj_val = fun(x0); nl.x0 = x0;
opts = optiset('solver','scip');

%% No scipopts Field
sopts = [];
sopts.maxiter = 1;
sopts.display = 3;
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Not a cell array
sopts.scipopts = 'a';
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Cell array of one column
sopts.scipopts = {'a';1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Empty Name
sopts.scipopts = {'',1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Empty Val (will skip it)
sopts.scipopts = {'a',[]};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Unknown Option
sopts.scipopts = {'limits/times',1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Wrong Datatype for Option
sopts.scipopts = {'limits/time','a'};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Outside Option Limits
sopts.scipopts = {'limits/time',-1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Bool Option logical
sopts.scipopts = {'display/lpinfo',false};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Bool Option double
sopts.scipopts = {'display/lpinfo',1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Char Option not char
clc
sopts.scipopts = {'lp/initalgorithm',1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Char Option not valid char
clc
sopts.scipopts = {'lp/initalgorithm','a'};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% Char Option
clc
sopts.scipopts = {'lp/initalgorithm','p'};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% String Option not string
clc
sopts.scipopts = {'nlp/solver',1};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% String Option unknown
clc
sopts.scipopts = {'nlp/solver','a'};
scip([],zs,[],[],[],lb,ub,[],[],[],nl,sopts)

%% opti_scipnl test
clc
fun = @(x) -x(1)*x(2)*x(3);
nlcon = @(x) -x(1)^2 - 2*x(2)^2 - 4*x(3)^2 + 48;
cl = 0;
cu = inf;
x0 = [1;1;1];
opts = optiset('display','iter','solverOpts',scipset('scipopts',{'nlp/disable',true}));
opti_scipnl(fun,[],[],[],[],[],nlcon,cl,cu,[],x0,opts);

%% opti test
clc
fun = @(x) (x(1)-x(2))^2 + (x(3)-1)^2 + (x(4)-1)^4 + (x(5)-1)^6;
nlcon = @(x) [x(1) + x(2) + x(3) + 4*x(4) - 7;
              x(3) + 5*x(5) - 6];
cl = [0;0];
cu = [0;0];
x0 = [10;7;2;-3;0.8];

sopts = scipset('scipopts',{'limits/solutions',1}); %find first feasible solution
opts = optiset('solver','scip','display','iter','solverOpts',sopts);
Opt = opti('fun',fun,'nl',nlcon,cl,cu,'x0',x0,'opts',opts)

[x,f,e,i] = solve(Opt)

%% MIQCQP1 [-2.5429] [non-convex]
clc
% Objective
H = eye(2); f = [-2;-2];
% Linear Constraints
A = [-1 1; 1 3]; b = [2;5];
% Quadratic Constraint
Q = [1 0;0 1]; l = [0;-2];
qrl = 3.5; qru = 5;
% Bounds
lb = [0;0]; ub = [40;inf];
% Integrality
xtype = 'IC';
% Options
sopts = scipset('scipopts',{'propagating/obbt/freq',1;...
                            'propagating/obbt/itlimitfactor',0;
                            'constraints/quadratic/maxproprounds',10;
                            'separating/maxrounds',20;
                            'limits/solutions',1});
opts = optiset('solver','scip','display','iter','solverOpts',sopts);
%Build & Solve
Opt = opti('H',H,'f',f,'ineq',A,b,'qcrow',Q,l,qrl,qru,'bounds',lb,ub,'xtype',xtype,'opts',opts)
[x,fval,exitflag,info] = solve(Opt)
% Plot
plot(Opt)

%% Global Emphasis
clc
fun = @(x) -x(1)*x(2)*x(3);
nlcon = @(x) -x(1)^2 - 2*x(2)^2 - 4*x(3)^2 + 48;
cl = 0;
cu = inf;
x0 = [1;1;1];

sopts = scipset('globalEmphasis','hardlp');
opts = optiset('solver','scip','display','iter','solverOpts',sopts);
Opt = opti('fun',fun,'nl',nlcon,cl,cu,'x0',x0,'opts',opts)

[x,f,e,i] = solve(Opt)

%% Presolving Emphasis
clc
fun = @(x) -x(1)*x(2)*x(3);
nlcon = @(x) -x(1)^2 - 2*x(2)^2 - 4*x(3)^2 + 48;
cl = 0;
cu = inf;
x0 = [1;1;1];

sopts = scipset('presolvingEmphasis','fast');
opts = optiset('solver','scip','display','iter','solverOpts',sopts);
Opt = opti('fun',fun,'nl',nlcon,cl,cu,'x0',x0,'opts',opts)

[x,f,e,i] = solve(Opt)

%% Heuristics Emphasis
clc
fun = @(x) -x(1)*x(2)*x(3);
nlcon = @(x) -x(1)^2 - 2*x(2)^2 - 4*x(3)^2 + 48;
cl = 0;
cu = inf;
x0 = [1;1;1];

sopts = scipset('heuristicsEmphasis','aggressive');
opts = optiset('solver','scip','display','iter','solverOpts',sopts);
Opt = opti('fun',fun,'nl',nlcon,cl,cu,'x0',x0,'opts',opts)

[x,f,e,i] = solve(Opt)

%% Separating Emphasis
clc
fun = @(x) -x(1)*x(2)*x(3);
nlcon = @(x) -x(1)^2 - 2*x(2)^2 - 4*x(3)^2 + 48;
cl = 0;
cu = inf;
x0 = [1;1;1];

sopts = scipset('separatingEmphasis','aggressive');
opts = optiset('solver','scip','display','iter','solverOpts',sopts);
Opt = opti('fun',fun,'nl',nlcon,cl,cu,'x0',x0,'opts',opts)

[x,f,e,i] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef SettingsLevelSetHorizontalFibers < SettingsLevelSetCreator

    properties (Access = public)
        levelOfFibers
        volume
    end

    methods (Access = public)

        function obj = SettingsLevelSetHorizontalFibers(varargin)
            obj.loadParams('paramsLevelSetCreator_HorizontalFibers.json')
            if nargin == 1
                obj.loadParams(varargin{1})
            end
        end

    end

end


text:
function pass = test_linearSystem1(pref)
% A linear CHEBOP test. This test tests a system of coupled ODEs, both with and
% without breakpoints in the domain of the CHEBOP. Problems are solved using
% chebcolloc1, chebcolloc2 and ultraS discretizations.

% Test 2x2 system (sin/cos)
% Toby Driscoll (extended by AB)

% Note: This test is taken from chebop_systemsolve1 in V4.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

%% Smooth domain, chebcolloc1:
d = [-pi pi];
A = chebop(@(x,u,v) [u - diff(v) ; diff(u) + v],d);
A.lbc = @(u,v) u + 1;
A.rbc = @(u,v) v;
x = chebfun('x',d);

pref.discretization = @chebcolloc1;
u = mldivide(A, 0, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = A.rbc(u1,u2);

pass(1) = norm( u1 - cos(x),inf) < 100*tol;
pass(2) = norm( u2 - sin(x),inf) < 100*tol;
pass(3) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;

%% Smooth domain, chebcolloc2
pref.discretization = @chebcolloc2;
u = mldivide(A, 0, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = A.rbc(u1,u2);

pass(4) = norm( u1 - cos(x),inf) < 100*tol;
pass(5) = norm( u2 - sin(x),inf) < 100*tol;
pass(6) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;

%% Smooth domain, ultraS
pref.discretization = @ultraS;
u = mldivide(A, 0, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = A.rbc(u1,u2);

pass(7) = norm( u1 - cos(x),inf) < 100*tol;
pass(8) = norm( u2 - sin(x),inf) < 100*tol;
pass(9) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;
%% Piecewise, chebcolloc1:
A.domain = [-pi 0 pi];
pref.discretization = @chebcolloc1;
u = mldivide(A, 0, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = A.rbc(u1,u2);
bcValLeft = bcFunLeft(d(1));
bcValRight = bcFunRight(d(2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(10) = norm( u1 - cos(x),inf) < 2000*tol;
pass(11) = norm( u2 - sin(x),inf) < 2000*tol;
pass(12) = norm(bcValLeft) < tol && norm(bcValRight) < tol;
pass(13) = norm(u1jump) < tol && norm(u2jump) < tol;

%% Piecewise, chebcolloc2:
pref.discretization = @chebcolloc2;
u = mldivide(A, 0, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = A.rbc(u1,u2);
bcValLeft = bcFunLeft(d(1));
bcValRight = bcFunRight(d(2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(14) = norm( u1 - cos(x),inf) < 2000*tol;
pass(15) = norm( u2 - sin(x),inf) < 2000*tol;
pass(16) = norm(bcValLeft) < tol && norm(bcValRight) < tol;
pass(17) = norm(u1jump) < tol && norm(u2jump) < tol;

%% Piecewise, ultraS:
pref.discretization = @ultraS;
u = mldivide(A, 0, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = A.rbc(u1,u2);
bcValLeft = bcFunLeft(d(1));
bcValRight = bcFunRight(d(2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(18) = norm( u1 - cos(x),inf) < 2000*tol;
pass(19) = norm( u2 - sin(x),inf) < 2000*tol;
pass(20) = norm(bcValLeft) < tol && norm(bcValRight) < tol;
pass(21) = norm(u1jump) < tol && norm(u2jump) < tol;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function midrad(c)
%MIDRAD       Display of interval hessians in midrad notation
%
%   midrad(c)
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 02/11/06     S.M. Rump  SparseInfNanFlag removed
% modified 06/04/09     S.M. Rump  Comment
% modified 08/26/12     S.M. Rump  global variables removed
%

  if nargin<2
    name = inputname(1);
    if isempty(name)                    % happens for display(hessianinit(random))
      name = 'ans';
    end
  end

  if isa(c.x,'intval')
    display_gen(c,name,'midrad')
  else
    display_gen(c,name,'display_')
  end


benchmarkDir = 'data/MOT16/train/';
allMets = evaluateTracking('c5-train.txt', 'results/', benchmarkDir);


text:
function e = sarimax_errors(parameters,p,q,constant,seasonal,y,x,sigma)
% PURPOSE:
%   Likelihood function for armaxfilter
%
% USAGE:
%   [LLF, LIKELIHOODS, ERRORS] = armaxfilter_likelihood(PARAMETERS,P,Q,CONSTANT,Y,X,M)
%   [LLF, LIKELIHOODS, ERRORS] = armaxfilter_likelihood(PARAMETERS,P,Q,CONSTANT,Y,X,M)
%
% INPUTS:
%   PARAMETERS  - A vector of GARCH process aprams of the form [constant, arch, garch]
%   P           - Vector containing lag indices of the AR component
%   Q           - Vector containing lag indices of the MA component
%   CONSTANT    - Value indicating whether the model contains a constant (1) or not (0)
%   Y           -
%   X           -
%   M           - Index to first element to use in the recursive residual calculation
%   SIGMA       - T by 1 vector
%
% OUTPUTS:
%   LLF         - Minus 1 times the log likelihood
%   LIKELIHOODS - Time series of likelihoods
%   ERRORS      - Time series of model errors
%
% COMMENTS:
%
% See also armaerrors

% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 4/1/2004

m = size(x,2);
[armaParameters, p, q] = sarma2arma(parameters(m+1:length(parameters)), p, q, seasonal);
parameters = [parameters(1:m);armaParameters];
e = armaxerrors(parameters,p,q,constant,y,x,m,sigma);

author: bashtage

repo: mfe-toolbox


text:
function [ iou ] = IOU( xi, X )
%IOU Intersection over union similarity measure for two bounding boxes.
%   Each sample must be composed of x = [ULx ULy BRx BRy]

    % Check area on xi
    xi_height = (xi(4) - xi(2) + 1);
    xi_width = (xi(3) - xi(1) + 1);
    xi_area = xi_height * xi_width;

    nSamples = size(X,1);
    iou = zeros(nSamples, 1);
    for j = 1:nSamples

        xj = X(j,:);

        % Check area on xj
        xj_height = (xj(4) - xj(2) + 1);
        xj_width = (xj(3) - xj(1) + 1);
        xj_area = xj_height * xj_width;

        % Check intersection
        count_intersect = rectint([xi(2), xi(1), xi_height, xi_width], [xj(2), xj(1), xj_height, xj_width]);

        % Calculate overlap score (intersection over union)
        iou(j) = count_intersect / (xi_area + xj_area - count_intersect);
    end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function tfi_hackOffSkull(normImg,fslStripImg,stripImg);
%
% tfi_hackOffSkull(normImg,fslStripImg,stripImg);
%
% Problem: there are 2 attempts to produce a skull-stripped
% image in the tfi path (one produced by fsl bet, and one
% produced by SkullStrip.tcl), and neither of them
% reliably do all the things needed.
%
% Solution: combine info from the two. This loads
% both the existing volumes, sets the norm img
% to zero wherever the fsl stripped img is zero,
% and saves the result as stripImg.
%
% This is a very crude hack designed to test the
% remaining code until Jonas updates SkullStrip.
%
%
% ras 02/05.
[img1 mmpervox1] = loadAnalyze(normImg);
[img2 mmpervox2] = loadAnalyze(fslStripImg);
if ~isequal(mmpervox1,mmpervox2)
    error('Volumes not equal size.')
end

img1(img2==0) = 0;
saveAnalyze(img1,stripImg(1:end-4),mmpervox1);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function C = xor(A,B)
%XOR Logical XOR for sptensors.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Observations for sparse matrix case.
% The result of xor(a,5) is dense!
% The result of xor(a,0) is dense!
% The result of xor(a,full(a)) is dense!
% The result of xor(a,b) is sparse.

%% Case 1: One argument is a scalar
if isscalar(B) || isa(B,'tensor')
    C = xor(full(A),B);
    return;
end
if isscalar(A)
    C = xor(A,full(B));
    return;
end

%% Case 2: Both x and y are tensors of some sort
if ~isequal(size(A),size(B))
    error('Must be tensors of the same size');
end

if isa(A,'sptensor') && isa(B,'sptensor')
    C = sptensor([A.subs; B.subs], 1, size(A), @(x) length(x) == 1);
    return;
end

if isa(B,'tensor')
    Bsubs = find(B ~= 0);
    C = sptensor([A.subs; Bsubs], 1, size(A), @(x) length(x) == 1);
    return;
end

%% Otherwise
error('The arguments must be two sptensors or an sptensor and a scalar.');

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function Script_CECLM_general()

addpath(genpath('../'));

% Replace this with the location of the 300W data location
if(exist([getenv('USERPROFILE') '/Dropbox/AAM/test data/'], 'file'))
    root_test_data = [getenv('USERPROFILE') '/Dropbox/AAM/test data/'];
else
    root_test_data = 'F:\Dropbox\AAM\test data/';
end
[images, detections, labels] = Collect_wild_imgs(root_test_data);

%% loading the CE-CLM model and parameters
[patches, pdm, clmParams, early_term_params] = Load_CECLM_general();
% Use the multi-hypothesis model, as bounding box tells nothing about
% orientation
views = [0,0,0; 0,-30,0; 0,30,0; 0,0,30; 0,0,-30;];
views = views * pi/180;

%% Setup recording

num_points = numel(pdm.M)/3;

shapes_all = zeros(size(labels,2),size(labels,3), size(labels,1));
labels_all = zeros(size(labels,2),size(labels,3), size(labels,1));
lhoods = zeros(numel(images),1);
all_lmark_lhoods = zeros(num_points, numel(images));
all_views_used = zeros(numel(images),1);

% Change if you want to visualize the outputs
verbose = false;
output_img = false;

if(output_img)
    output_root = './ceclm_gen_out/';
    if(~exist(output_root, 'dir'))
        mkdir(output_root);
    end
end
if(verbose)
    f = figure;
end

%% Fitting the model to the provided images

tic
for i=1:numel(images)

    image = imread(images(i).img);
    image_orig = image;

    if(size(image,3) == 3)
        image = rgb2gray(image);
    end

    bbox = detections(i,:);

    [shape,~,~,lhood,lmark_lhood,view_used] =...
        Fitting_from_bb_multi_hyp(image, [], bbox, pdm, patches, clmParams, views, early_term_params);

    all_lmark_lhoods(:,i) = lmark_lhood;
    all_views_used(i) = view_used;

    shapes_all(:,:,i) = shape;
    labels_all(:,:,i) = labels(i,:,:);

    if(mod(i, 200)==0)
        fprintf('%d done\n', i );
    end

    lhoods(i) = lhood;

    if(output_img)
        v_points = sum(squeeze(labels(i,:,:)),2) > 0;
        DrawFaceOnImg(image_orig, shape, sprintf('%s/%s%d.jpg', output_root, 'fit', i), bbox, v_points);
    end

    if(verbose)
        v_points = sum(squeeze(labels(i,:,:)),2) > 0;
        DrawFaceOnFig(image_orig, shape, bbox, v_points);
    end
end
toc

experiment.errors_normed = compute_error(labels_all, shapes_all + 1.0);
experiment.lhoods = lhoods;
experiment.shapes = shapes_all;
experiment.labels = labels_all;
experiment.all_lmark_lhoods = all_lmark_lhoods;
experiment.all_views_used = all_views_used;

fprintf('Done: mean normed error %.3f median normed error %.4f\n', ...
    mean(experiment.errors_normed), median(experiment.errors_normed));

%%
output_results = 'results/results_ceclm_general.mat';
save(output_results, 'experiment');

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function res = indfrequency(this, f)
% Method for getting the index closest to given frequency
% FORMAT  res = indfrequency(this, f)
% this       - MEEG object
% f          - vector of frequencies (in Hz)
%
% res        - vector of sample indices matching indices
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: indfrequency.m 3254 2009-07-07 15:18:54Z vladimir $

if ~strncmpi(transformtype(this),'TF',2)
    error('Only TF datasets are supported');
end

res = NaN(1,length(f));
fdiff = mean(diff(frequencies(this)));
if nsamples(this) > 0
    F = frequencies(this);
    for i = 1:length(f)
        [m, res(i)] = min(abs(F-f(i)));
        if m > fdiff
            warning('Could not find an index matching the requested frequency %d Hz', f(i));
            res(i) = NaN;
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clc; clear all; close all;

% prepare the input mesh
[V,F] = readOBJ('../data/spot.obj');
[V,F,S] = loop(V,F);

% specify handles and handle displacement
handles = [1837; 2274; 1144; 1454];
handles_disp = [0.2,0,0;0.2,0,0;0.2,0,0;-0.8,0,0];

% solve for deformation
A = cotmatrix(V,F);
B = zeros(size(V,1),3);
tic;
preF = [];
[d1, preF] = ...
    min_quad_with_fixed(...
    A,B,handles,handles_disp,[],[],preF);
toc;

figure(1)
subplot(1,2,1)
tsurf(F,V)
hold on
sct(V(handles,:),'filled', 'r')
qvr(V(handles,:), handles_disp, 'b')
axis equal
title('initial')
subplot(1,2,2)
tsurf(F,V+d1)
axis equal
title('deformation')

% solve for new deformation
handles_disp_new = [0.2,0.2,0;0.2,0.2,0;0.2,-0.2,0;-1,0.3,0];

tic;
[d2, preF] = ...
    min_quad_with_fixed(...
    A,B,handles,handles_disp_new,[],[],preF);
toc;

figure(2)
subplot(1,2,1)
tsurf(F,V)
hold on
sct(V(handles,:),'filled', 'r')
qvr(V(handles,:), handles_disp_new, 'g')
axis equal
title('initial')
subplot(1,2,2)
tsurf(F,V+d2)
axis equal
title('deformation')

author: odedstein

repo: sgi-introduction-course


text:
% OP_SU_EV_TP: assemble the matrix A = [a(i,j)], a(i,j) = 1/2 (sigma (u_j), epsilon (v_i)), exploiting the tensor product structure.
%
%   mat = op_su_ev_tp (spu, spv, msh, lambda, mu);
%   [rows, cols, values] = op_su_ev_tp (spu, spv, msh, lambda, mu);
%
% INPUT:
%
%   spu:     object representing the space of trial functions (see sp_vector)
%   spv:     object representing the space of test functions (see sp_vector)
%   msh:     object that defines the domain partition and the quadrature rule (see msh_cartesian)
%   lambda, mu: function handles to compute the Lame' coefficients
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = op_su_ev_tp (space1, space2, msh, lambda, mu)

  for icomp = 1:space1.ncomp_param
    for idim = 1:msh.ndim
      size1 = size (space1.scalar_spaces{icomp}.sp_univ(idim).connectivity);
      size2 = size (space2.scalar_spaces{icomp}.sp_univ(idim).connectivity);
      if (size1(2) ~= size2(2) || size1(2) ~= msh.nel_dir(idim))
        error ('One of the discrete spaces is not associated to the mesh')
      end
    end
  end

  A = spalloc (space2.ndof, space1.ndof, 5*space1.ndof);

  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    sp1_col = sp_evaluate_col (space1, msh_col, 'value', false, ...
                               'gradient', true, 'divergence', true);
    sp2_col = sp_evaluate_col (space2, msh_col, 'value', false, ...
                               'gradient', true, 'divergence', true);

    for idim = 1:msh.rdim
      x{idim} = reshape (msh_col.geo_map(idim,:,:), msh_col.nqn, msh_col.nel);
    end

    A = A + op_su_ev (sp1_col, sp2_col, msh_col, lambda (x{:}), mu (x{:}));
  end

  if (nargout == 1)
    varargout{1} = A;
  elseif (nargout == 3)
    [rows, cols, vals] = find (A);
    varargout{1} = rows;
    varargout{2} = cols;
    varargout{3} = vals;
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function failed_bug2396

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_read_headshape

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2396'));

%% define the files

mrifile = 'nobias_KR.nii';
labelfile = 'segmentation/brainsuite/nobias_KR.label.nii';

dfsfile = {
  'segmentation/brainsuite/nobias_KR.brain.dfs'
  'segmentation/brainsuite/nobias_KR.inner_skull.dfs'
  'segmentation/brainsuite/nobias_KR.outer_skull.dfs'
  'segmentation/brainsuite/nobias_KR.scalp.dfs'
  };

meshfile = {
  'segmentation/brainvisa/nobias_KR_BEM_brain.mesh'
  'segmentation/brainvisa/nobias_KR_BEM_head.mesh'
  'segmentation/brainvisa/nobias_KR_BEM_skull.mesh'
  'segmentation/brainvisa/nobias_KR_head.mesh'
  };

% the following contain metadata that goes with the mesh files
% segmentation/brainvisa/nobias_KR_BEM_brain.mesh.minf
% segmentation/brainvisa/nobias_KR_BEM_head.mesh.minf
% segmentation/brainvisa/nobias_KR_BEM_skull.mesh.minf
% segmentation/brainvisa/nobias_KR_head.mesh.minf

surffile = {
  'segmentation/freesurfer/lh.BEM_brain_surface.surf'
  'segmentation/freesurfer/lh.BEM_inner_skull_surface.surf'
  'segmentation/freesurfer/lh.BEM_outer_skin_surface.surf'
  'segmentation/freesurfer/lh.BEM_outer_skull_surface.surf'
  };

%% read the content from the files

mri = ft_read_mri(mrifile);

% it would be more consistent to read it with ft_read_atlas, as if it were a ft_datatype_segmentation
seg = ft_read_mri(labelfile);

% dress up the segmentation
seg.anatomylabel = {};
seg.anatomylabel{19} = 'brain';
seg.anatomylabel{18} = 'csf';
seg.anatomylabel{17} = 'skull';
seg.anatomylabel{16} = 'scalp';
seg.coordsys = 'ctf'; % it seems to be in CTF coordinates, i.e. positive y axis through nasion

cfg = [];
cfg.atlas = seg; % use itself as the atlas to look up the number-name mappings
ft_sourceplot(cfg, seg);

%% read the content from the files

dfs = {};
figure
for i=1:4
  dfs{i} = ft_read_headshape(dfsfile{i});
  ft_plot_mesh(dfs{i}, 'edgecolor', 'none', 'facecolor', [1 1 1]/2, 'facealpha', 0.3);
end

mesh = {};
figure
for i=1:4
  mesh{i} = ft_read_headshape(meshfile{i});
  ft_plot_mesh(mesh{i}, 'edgecolor', 'none', 'facecolor', [1 1 1]/2, 'facealpha', 0.3);
  % note that mesh number 4 is ugly and has internal structures
end

surf = {};
figure
for i=1:4
  surf{i} = ft_read_headshape(surffile{i});
  ft_plot_mesh(surf{i}, 'edgecolor', 'none', 'facecolor', [1 1 1]/2, 'facealpha', 0.3);
end

% plot them together, this should form a consistent set
ft_plot_mesh(dfs{1}, 'edgecolor', 'none', 'facecolor', [1 1 1]/2, 'facealpha', 0.3);
ft_plot_mesh(mesh{3}, 'edgecolor', 'none', 'facecolor', [1 1 1]/2, 'facealpha', 0.3);
ft_plot_mesh(surf{3}, 'edgecolor', 'none', 'facecolor', [1 1 1]/2, 'facealpha', 0.3);

dfsall = cellfun(@ft_read_headshape, dfsfile);


text:
function img_new = dtiResliceAnalyzeSn3d(img, sn3dParams, nonLinearFlag, mmPerVoxOut)
% img_new = dtiResliceAnalyzeSn3d(img, sn3dParams, [nonLinearFlag], [mmPerVoxOut])
%
% Applies spm-style sn3d parameters to a volume loaded using loadAnalyze.
% These sn3d parameters include a linear (affine) component and a non-linear component.
%
% INPUT:
%   img                 XxYxZ image array.
%   sn3dParams          The sn3d parameter structure from SPM.
%   nonLinearFlag       If 1 (default), non-linear component is included.
%                       If 0, then non-linear component is not included.
%   mmPerVoxOut         mm per voxel in resliced array.
%
% OUTPUT:
%   img_new             XxYxZ image array.
%
% REQUIRES:
% spm99 or spm2 in the path.
%
% HISTORY:
% 2004.09.13 ASH & RFD Wrote it.
%
% Example:
% sn3dParams = load('dti_analyses/may021126_B0_sn3d.mat');
% B0 = loadAnalyze('/snarp/u1/dti/adultData/may021126/dti_analyses/may021126_B0.hdr');
% B0_new = dtiResliceAnalyzeSn3d(B0, sn3dParams, mmPerVoxOut);
% figure; imagesc(B0_new(:,:,23)); axis equal; colormap gray; axis xy;

if(~exist('nonLinearFlag') | isempty(nonLinearFlag))
    nonLinearFlag = 1;
end
if(~exist('mmPerVoxOut') | isempty(mmPerVoxOut))
    mmPerVoxOut = [2 2 2];
end

% Create a bounding box for the resliced data.
% The following is the spm default bounding box. Note that the box is
% defined in Talairach space (units = mm).
bb = [-78 -112 -50;
       78  76   85];

x = (bb(1,1):mmPerVoxOut(1):bb(2,1));
y = (bb(1,2):mmPerVoxOut(2):bb(2,2));
z = (bb(1,3):mmPerVoxOut(3):bb(2,3));

mm = sn3dParams.Dims(6,:);

[X,Y,Z] = meshgrid(x, y, z);
clear x y z;

talCoords = [X(:) Y(:) Z(:)];
newSize = size(X);
clear X Y Z;

imgCoords = mrAnatGetImageCoordsFromSn(sn3dParams, talCoords, nonLinearFlag);
imgCoords = imgCoords./repmat(mm',1,size(imgCoords,2));
imgOrigin = sn3dParams.Dims(5,:)'/2;
imgCoords = imgCoords+repmat(imgOrigin,1,size(imgCoords,2));
% We now have image coords for the space defined by the bounding box (in
% Talairach space). These image coords are in the space of the original
% image that was used to compute the sn3d params. But that may not be the
% same space as the tensor image.

% Following lines correspond to canonical transform.
% This is unnecessary because that has already been applied in Analyze format
% imgCoords = xformToTensor*[imgCoords;ones(1,size(imgCoords,2))];
% imgCoords = imgCoords(1:3,:);
% imgCoords(1,:) = imgCoords(1,:).*mm(1);
% imgCoords(2,:) = imgCoords(2,:).*mm(2);
% imgCoords(3,:) = imgCoords(3,:).*mm(3);

% Permute dimensions because of Analyze format
img = permute(img,[2,1,3]);

% Interpolate coordinate grid
img_new = myCinterp3(img, [size(img,1) size(img,2)], size(img,3), imgCoords', 0.0);
img_new = reshape(img_new, [newSize, 1]);

return;


text:
classdef IN1KMOP7 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'mnv3';
            config.args.objs = 'err&params';
            config.args.normalized_objectives = false;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef ConvergenceVariables < handle

    properties (GetAccess = public, SetAccess = private)
        values
        nVar
    end

    methods (Access = public)

        function obj = ConvergenceVariables(n)
            obj.nVar = n;
        end

        function reset(obj)
            obj.values = [];
            obj.nVar = 0;
        end

        function append(obj,value)
            obj.values(obj.nVar+1) = value;
            obj.nVar = obj.nVar + 1;
        end

        function set(obj,value)
            obj.values(end) = value;
        end

    end

end


text:
function determ = rectangle_adj_determinant ( row_num, col_num )

%*****************************************************************************80
%
%% RECTANGLE_ADJ_DETERMINANT returns the determinant of the RECTANGLE_ADJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ROW_NUM, COL_NUM, the number of rows and
%    columns in the rectangle.
%
%    Output, real DETERM, the determinant.
%

%
%  If ROW_NUM == 1 or COL_NUM == 1 we have a case of the LINE_ADJ matrix.
%
  if ( row_num == 1 )

         if ( mod ( row_num, 4 ) == 1 )
      determ =   0.0;
    elseif ( mod ( row_num, 4 ) == 2 )
      determ = - 1.0;
    elseif ( mod ( row_num, 4 ) == 3 )
      determ =   0.0;
    elseif ( mod ( row_num, 4 ) == 0 )
      determ = + 1.0;
    end

  elseif ( col_num == 1 )

         if ( mod ( col_num, 4 ) == 1 )
      determ =   0.0;
    elseif ( mod ( col_num, 4 ) == 2 )
      determ = - 1.0;
    elseif ( mod ( col_num, 4 ) == 3 )
      determ =   0.0;
    elseif ( mod ( col_num, 4 ) == 0 )
      determ = + 1.0;
    end
%
%  Otherwise, we can form at least one square, hence a null vector,
%  hence the matrix is singular.
%
  else

    determ = 0.0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function filteredM = tsmooth(m,fSize,fSd)
%
% filteredM = tsmooth(m,fSize,fSd)
%
%AUTHOR:  Engel, Wandell
%DATE:    1993-4
%PURPOSE:
%   Take a set of time series values in the columns of m and smooth
% each one by convolution with a Gaussian kernel of support fSize and
% standard deviation of fSd.
%

kernel = mkGaussKernel([1 fSize],[1 fSd]);
s = size(m);

% 7/09/97 Lea updated to 5.0
filteredM = [];

for i = (1:s(2))
 tmp = evenConv(m(:,i),kernel);
% tmp = mrDeTrend(tmp);
 filteredM = [filteredM;tmp];
end


text:
function out = saveobj(obj)
%SAVEOBJ (overloaded)

% We have to save the persistent variables in the SDPVAR class
obj.savedata.internalsdpvarstate = yalmip('getinternalsdpvarstate');
obj.savedata.version = yalmip('version');
out = obj;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function out = not( in )
%NOT for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

[out, idx] = flattenCellMatrix(in.data);
out = cellfun(@(x) ~x, out, 'UniformOutput', false);
out = reconFlatCellMatrix(out,idx);
out = TRAFO(out,in.meta);

end

author: thomaskuestner


text:
function channel_loc2tvs

% start at the root and convert the channel_loc.txt files for BIDS
% SPM12 dependency

folders = dir('sub*')
for f=7:length(folders)
    cd(folders(f).name)
    file = dir('*.txt')
    fileID = fopen(file.name);
    C = textscan(fileID,'%d %f32 %f32 %f32 %s %s');
    fclose(fileID);

    % do the electrodes.tsv
    clear name
    for i=1:74
        name(i,:) = cell2mat(C{5}(i));
    end
    x = C{2}(1:74);
    y = C{3}(1:74);
    z = C{4}(1:74);
    T = table(name,x,y,z);
    writetable(T,[pwd filesep 'eeg' filesep folders(f).name '_task-facerecognition_electrodes.txt'],'Delimiter','\t')
    movefile([pwd filesep 'eeg' filesep folders(f).name '_task-facerecognition_electrodes.txt'],[pwd filesep 'eeg' filesep folders(f).name '_task-facerecognition_electrodes.tsv'])

    % do the electrodes.json
    json = struct('EEGCoordinateSystem','T1w', ...
        'EEGCoordinateUnits','mm','IntendedFor',['anat' filesep folders(f).name '_T1w.nii.gz']);
    spm_jsonwrite(['eeg' filesep folders(f).name '_task-facerecognition_electrodes.json'],json,struct('indent','  '))

    % do the fid.json
    json = struct(...
        'EEGCoordinateSystem','T1w', ...
        'EEGCoordinateUnits', 'mm', ...
        'AnatomicalMRICoordinateSystem', 'T1w', ...
        'AnatomicalMRICoordinateUnits', 'mm', ...
        'LandmarkCoordinates', ' ', ...
        'LandmarkCoordinateSystem','T1w', ...
        'LandmarkCoordinateUnits','mm');

    clear name
    for i=75:77
        tmp = cell2mat(C{5}(i));
        name(i-74,:) = tmp(1:3);
    end
        json.LandmarkCoordinates = ['{"' name(1,:) '":[' num2str(double([C{2}(75),C{3}(75),C{4}(75)])) '], "' ...
       name(2,:) '":[' num2str(double([C{2}(76),C{3}(76),C{4}(76)])) '], "' ...
        name(3,:) '":[' num2str(double([C{2}(77),C{3}(77),C{4}(77)])) ']}'];

    spm_jsonwrite(['eeg' filesep folders(f).name '_task-facerecognition_fid.json'],json,struct('indent','  '))

    % clean up
    delete(file.name)
    cd ..
end

author: bids-standard

repo: bids-examples


text:
function [out] = saturation_3(S,Smax,p1,In)
%saturation_3

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Saturation excess from a store with different degrees of
%               saturation (exponential variant)
% Constraints:  -
% @(Inputs):    S    - current storage [mm]
%               Smax - maximum contributing storage [mm]
%               p1   - linear scaling parameter [-]
%               In   - incoming flux [mm/d]

out = (1-(1/(1+exp((S/Smax + 0.5)/p1)))).*In;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
% test for low rank approximation of an image

name = 'lena';
n = 128;
M = load_image(name);
M = crop(M,n);

[U,S,V] = svd(M);

%% display reconstruction
m_list = [5 10 50 100];
for i = 1:length(m_list)
    m = m_list(i);
    s = diag(S);
    s(m+1:end) = 0;
    M1 = U*diag(s)*V';
    subplot(2,2,i);
    imagesc(M1); axis image; axis off;
    title([num2str(m) ' eigv']);
end
colormap gray(256);

return;

% test for pseudo inverse
n = 300; m = 4*n;
D = rand(n,m);
tic;
B = D' * (D*D')^(-1);
t = toc;
tic;
B1 = pinv(D);
t1 = toc;
e = norm(B-B1, 'fro')/norm(B, 'fro');
disp(['speed-up=' num2str((t1-t)/t*100) '%, err=' num2str(e*100)]);

author: gpeyre

repo: matlab-toolboxes