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

% Last Modified by GUIDE v2.5 13-Jan-2012 17:19:05

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @equalizer_OpeningFcn, ...
                   'gui_OutputFcn',  @equalizer_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 equalizer is made visible.
function equalizer_OpeningFcn(hObject, eventdata, handles, varargin)
global s Fs
global g1 g2 n hs
global po
global pla pau
global hdls
global vol
global Q f filts_type
global hp hpp

% g1 g2 limits og gains
% po buffer size in seconds

% s - signal
% Fs - sampling frequency
% g1 g2 limits og gains
% po buffer size in seconds

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

% Update handles structure
guidata(hObject, handles);

hdls=handles; % save all handles to global variable

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

g1=-10;
g2=10;
n=10; % mber of sliders
f=[60  170 310  600 1e3 3e3   6e3  12e3  14e3  16e3]; % frequencyes
%f=0.7*exp(1:10);
%f=24*exp(1:6);
% f=[43 118  321 873 2374 6454 17540];
% f=[60  310  600 1e3   6e3  12e3   16e3]; % frequencyes
Q=1*ones(1,n); % Q-factors
filts_type=[1 2*ones(1,n-2) 3]; % filt types: 1-low shelving, 2-peak, 3-high shelving

po=0.2; % buffers size, seconds
pla=false; % true if play
pau=false; % true if pause

vol=0.7; % volume
ps=get(handles.szg,'position');
pu=get(handles.szg,'units');

sz=ps(4);

set(handles.vol,'value',vol);

hp=semilogx(NaN,NaN,'r-','parent',handles.axes1);
set(handles.axes1,'NextPlot','add');
hpp=semilogx(NaN,NaN,'ok','parent',handles.axes1);
xlabel(handles.axes1,'frequrency, Hz');
ylabel(handles.axes1,'gain, dB');

pg1=get(handles.g1,'position');
pg2=get(handles.g2,'position');
szy=pg2(2)-pg1(2);

x=ps(1)+sz;
y=ps(2)-1*sz-szy;

pa=get(handles.axes1,'position');
sxa=pa(3);

hs=zeros(1,n);

for nc=1:n
    hs(nc)=uicontrol('Style', 'slider',...
            'Min',g1,'Max',g2,'Value',0,...
            'units',pu,...
            'Position', [x y 1.5*sz szy],...
            'Callback', 'recalculate_filter;');

    x=x+sxa/n;
end

Fs=44100;
recalculate_filter;

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

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

function file_Callback(hObject, eventdata, handles)
fln=get(handles.file,'string');
set(handles.reading,'visible','on');
drawnow;
read_audio(fln);
set(handles.reading,'visible','off');
drawnow;

% --- Executes during object creation, after setting all properties.
function file_CreateFcn(hObject, eventdata, handles)
% hObject    handle to file (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 browse.
function browse_Callback(hObject, eventdata, handles)
[FileName,PathName,FilterIndex] = uigetfile({'*.mp3;*.wav;*.wma;*.avi;*.mpg','media files (*.mp3,*.wav,*.wma,*.avi,*.mpg)';'*.*',  'All Files (*.*)'},'File Selector') ;
isok=true;
if length(FileName)==1
     if FileName==0
         isok=false;
     end
end
if isok
    fln=[PathName FileName];
    set(handles.file,'string',fln);
    set(handles.reading,'visible','on');
    drawnow;
    read_audio(fln);
    set(handles.reading,'visible','off');
    drawnow;

end

% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (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
% disp('slider releized');

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

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

% --- Executes on slider movement.
function vol_Callback(hObject, eventdata, handles)
% hObject    handle to vol (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

global vol
vol=get(hObject,'Value');

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

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

% --- Executes on button press in szg.
function szg_Callback(hObject, eventdata, handles)
global n hs
for nc=1:n
    set(hs(nc),'value',0);
end
drawnow;
recalculate_filter;

% --- Executes on button press in play.
function play_Callback(hObject, eventdata, handles)
open_device;

% --- Executes on button press in pause.
function pause_Callback(hObject, eventdata, handles)
global pau
pau=~pau;
if pau
    set(handles.pause,'string','Unpause');
    drawnow;
else
    set(handles.pause,'string','Pause');
    drawnow;
end

% --- Executes on button press in stop.
function stop_Callback(hObject, eventdata, handles)
global pla
pla=false;

% --- Executes on button press in sset.
function sset_Callback(hObject, eventdata, handles)
global po pso pco s
L=size(s,1);
v=get(handles.slider1,'value');
pco=round(v*L/pso);
if pco<1
    pco=1;
end


text:
function F = OCP_F(u,x,p,discretizationMethod,isMEnabled,parIdx) %#codegen

    if isMEnabled
        fdt = OCP_GEN_fdt(u,x,p,parIdx);
        % forced to 'Euler'
        M = OCP_GEN_M(u,x,p);
        invM = inv(M);
        invMf = invM*fdt;
        % F
        F = invMf - x;
    else % M disabled
        switch discretizationMethod
            case 'RK2'
                fdt   = OCP_GEN_fdt(u,x,p,parIdx);
                f12dt = OCP_GEN_fdt(u,x-fdt/2,p,parIdx);

                F  =  f12dt  - x;
            case 'RK4'
                fdt   = OCP_GEN_fdt(u,x,p,parIdx);   % k1
                f12dt = OCP_GEN_fdt(u,x-fdt/2,p,parIdx); % k2
                f22dt = OCP_GEN_fdt(u,x-f12dt/2,p,parIdx); % k3
                f3dt  = OCP_GEN_fdt(u,x-f22dt,p,parIdx); % k4

                F = (fdt+2*f12dt+2*f22dt+f3dt)/6 - x;

            otherwise % 'Euler'
                fdt = OCP_GEN_fdt(u,x,p,parIdx);
                F  = fdt - x;
        end
    end
end

author: deng-haoyang

repo: ParNMPC


text:
function value = daub8_condition ( n )

%*****************************************************************************80
%
%% DAUB8_DETERMINANT returns the L1 condition of the DAUB8 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real VALUE, the L1 condition.
%
  c = [ ...
    0.2303778133088964, ...
    0.7148465705529154, ...
    0.6308807679298587, ...
   -0.0279837694168599, ...
   -0.1870348117190931, ...
    0.0308413818355607, ...
    0.0328830116668852, ...
   -0.0105974017850690 ];

  a_norm = sum ( abs ( c(1:8) ) );
  b_norm = a_norm;
  value = a_norm * b_norm;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% CONTRAST2_FSU uses the fsuClusterMatlab command to run the CONTRAST2 computation.
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    [] allows us to specify an output directory;
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%    'm" indicates this is an "MPI-like" job;
%    'w' indicates that the MATLAB session should wait for completion;
%    4 is the number of processors requested;
%    @contrast2_fun indicates the function we want to execute;
%    {} would contain input arguments if there were any.
%
%    The function "contrast2_fun" must correspond to a MATLAB M-file "contrast2_fun.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2010
%
%  Author:
%
%    John Burkardt
%

%
%  Read a black and white image.
%
  x = imread ( 'surfsup.tif' );

  fsuClusterMatlab ( [], [], 'm', 'w', 4, @contrast2_fun, { x } )


text:
function [fids, positions] = op_CSIExtractFids(MRSIStruct, kSpaceFile)
    table = readTable(kSpaceFile);
    TR = table.TR;
    Kx = table.Kx;
    Ky = table.Ky;
    t = table.t;
    [MRSIStruct, prevPermute, prevShape] = reshapeDimensions(MRSIStruct, {'t', 'ky', 'kx'});

end

author: CIC-methods

repo: FID-A

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


text:
% resample the set of particles.
% A particle has a probability proportional to its weight to get
% selected. A good option for such a resampling method is the so-called low
% variance sampling, Probabilistic Robotics pg. 109
function newParticles = resample(particles)

numParticles = length(particles);

w = [particles.weight];

% normalize the weight
w = w / sum(w);

% consider number of effective particles, to decide whether to resample or not
useNeff = false;
%useNeff = true;
if useNeff
  neff = 1. / sum(w.^2);
  neff
  if neff > 0.5*numParticles
    newParticles = particles;
    for i = 1:numParticles
      newParticles(i).weight = w(i);
    end
    return;
  end
end

newParticles = struct;

% TODO: implement the low variance re-sampling

% the cummulative sum
cs = cumsum(w);
weightSum = cs(length(cs));

% initialize the step and the current position on the roulette wheel
step = weightSum / numParticles;
position = unifrnd(0, weightSum);
idx = 1;

% walk along the wheel to select the particles
for i = 1:numParticles
  position += step;
  if (position > weightSum)
    position -= weightSum;
    idx = 1;
  end
  while (position > cs(idx))
    idx++;
  end
  newParticles(i) = particles(idx);
  newParticles(i).weight = 1/numParticles;
end

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function setglobal(var)

% SETGLOBAL assigns global variables from a structure
%
% Use as
%   var = getglobal;
%   setglobal(var);
%
% See also GETGLOBAL

% Copyright (C) 2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isempty(var)
  return
end

list = fieldnames(var);

for i=1:length(list)
  eval(sprintf('global %s', list{i}));
  eval(sprintf('%s = var.%s;', list{i}, list{i}));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [template] = tapas_load_template(fname)
%% Loads a file that work as a template for filling.
%
% Input
%   fname       -- File name.
%
% Output
%   string      -- An array containing the string.
%

% aponteeduardo@gmail.com
% copyright (C) 2017
%

fid = fopen(fname);
try
    template = {};
    while feof(fid) == 0
       template{end + 1} = fgets(fid);
    end
catch err
    % Make sure that the file is closed properly.
    try
        fclose(fid);
    end
    rethrow(err);
end
fclose(fid);
template = [template{:}];

end

author: translationalneuromodeling

repo: tapas


text:
% Script comes from Convolutional Pose Machine
function genJSON(dataset)

    addpath('util/jsonlab/');

    if(strcmp(dataset, 'MPI'))
        mat = load('../dataset/MPI/mpii_human_pose_v1_u12_1/mpii_human_pose_v1_u12_1.mat');
        RELEASE = mat.RELEASE;
        trainIdx = find(RELEASE.img_train);

        % in MPI: (0 - r ankle, 1 - r knee, 2 - r hip, 3 - l hip, 4 - l knee,
        %          5 - l ankle, 6 - pelvis, 7 - thorax, 8 - upper neck, 9 - head top,
        %          10 - r wrist, 11 - r elbow, 12 - r shoulder, 13 - l shoulder, 14 - l elbow, 15 - l wrist)"

        tompson = load('../dataset/MPI/Tompson_valid/detections');
        tompson_i_p = [tompson.RELEASE_img_index; tompson.RELEASE_person_index];

        count = 1;
        makeFigure = 0;
        validationCount = 0;

        for i = trainIdx
            numPeople = length(RELEASE.annolist(i).annorect);
            fprintf('image: %d (numPeople: %d) last: %d\n', i, numPeople, trainIdx(end));

            for p = 1:numPeople
                loc = find(sum(~bsxfun(@minus, tompson_i_p, [i;p]))==2, 1);
                if(~isempty(loc))
                    validationCount = validationCount + 1;
                    %fprintf('Tomspon''s validation! %d\n', validationCount);
                    isValidation = 1;
                else
                    isValidation = 0;
                end
                joint_all(count).dataset = 'MPI';
                joint_all(count).isValidation = isValidation;

                try % sometimes no annotation at all....
                    anno = RELEASE.annolist(i).annorect(p).annopoints.point;
                catch
                    continue;
                end

                % set image path
                %joint_all(count).img_paths = sprintf('dataset/MPI/images/%s', RELEASE.annolist(i).image.name);
                joint_all(count).img_paths = RELEASE.annolist(i).image.name;

                img = imread(['../dataset/MPI/images/', joint_all(count).img_paths]);
                [h,w,~] = size(img);
                %[h,w,~] = size(imread(['../dataset/MPI/images/', joint_all(count).img_paths]));

                joint_all(count).img_width = w;
                joint_all(count).img_height = h;
                joint_all(count).objpos = [RELEASE.annolist(i).annorect(p).objpos.x, RELEASE.annolist(i).annorect(p).objpos.y];
                % set part label: joint_all is (np-3-nTrain)

                % for this very center person
                for part = 1:length(anno)
                    joint_all(count).joint_self(anno(part).id+1, 1) = anno(part).x;
                    joint_all(count).joint_self(anno(part).id+1, 2) = anno(part).y;
                    try % sometimes no is_visible...
                        if(anno(part).is_visible == 0 || anno(part).is_visible == '0')
                            joint_all(count).joint_self(anno(part).id+1, 3) = 0;
                        else
                            joint_all(count).joint_self(anno(part).id+1, 3) = 1;
                        end
                    catch
                        joint_all(count).joint_self(anno(part).id+1, 3) = 1;
                    end
                end
                % pad it into 16x3
                dim_1 = size(joint_all(count).joint_self, 1);
                dim_3 = size(joint_all(count).joint_self, 3);
                pad_dim = 16 - dim_1;
                joint_all(count).joint_self = [joint_all(count).joint_self; zeros(pad_dim, 3, dim_3)];
                % set scale
                joint_all(count).scale_provided = RELEASE.annolist(i).annorect(p).scale;

                % for other person on the same image
                count_other = 1;
                joint_others = cell(0,0);
                for op = 1:numPeople
                    if(op == p), continue; end
                    try % sometimes no annotation at all....
                        anno = RELEASE.annolist(i).annorect(op).annopoints.point;
                    catch
                        continue;
                    end
                    joint_others{count_other} = zeros(16,3);
                    for part = 1:length(anno)
                        joint_all(count).joint_others{count_other}(anno(part).id+1, 1) = anno(part).x;
                        joint_all(count).joint_others{count_other}(anno(part).id+1, 2) = anno(part).y;
                        try % sometimes no is_visible...
                            if(anno(part).is_visible == 0 || anno(part).is_visible == '0')
                                joint_all(count).joint_others{count_other}(anno(part).id+1, 3) = 0;
                            else
                                joint_all(count).joint_others{count_other}(anno(part).id+1, 3) = 1;
                            end
                        catch
                            joint_all(count).joint_others{count_other}(anno(part).id+1, 3) = 1;
                        end
                    end
                    joint_all(count).scale_provided_other(count_other) = RELEASE.annolist(i).annorect(op).scale;
                    joint_all(count).objpos_other{count_other} = [RELEASE.annolist(i).annorect(op).objpos.x RELEASE.annolist(i).annorect(op).objpos.y];

                    count_other = count_other + 1;
                end

                if(makeFigure) % visualizing to debug
                    %imshow(joint_all(count).img_paths);
                    imshow(img);
                    hold on;
                    visiblePart = joint_all(count).joint_self(:,3) == 1;
                    invisiblePart = joint_all(count).joint_self(:,3) == 0;
                    plot(joint_all(count).joint_self(visiblePart, 1), joint_all(count).joint_self(visiblePart,2), 'gx');
                    plot(joint_all(count).joint_self(invisiblePart,1), joint_all(count).joint_self(invisiblePart,2), 'rx');
                    plot(joint_all(count).objpos(1), joint_all(count).objpos(2), 'cs');
                    if(~isempty(joint_all(count).joint_others))
                        for op = 1:size(joint_all(count).joint_others, 3)
                            visiblePart = joint_all(count).joint_others(:,3,op) == 1;
                            invisiblePart = joint_all(count).joint_others(:,3,op) == 0;
                            plot(joint_all(count).joint_others(visiblePart,1,op), joint_all(count).joint_others(visiblePart,2,op), 'mx');
                            plot(joint_all(count).joint_others(invisiblePart,1,op), joint_all(count).joint_others(invisiblePart,2,op), 'cx');
                        end
                    end
                    close all;
                end
                joint_all(count).annolist_index = i;
                joint_all(count).people_index = p;
                joint_all(count).numOtherPeople = length(joint_all(count).joint_others);
                count = count + 1;

                %if(count==10), break; endscale_provided
            end
            %if(count==10), break; end
        end
        opt.FileName = 'json/MPI_annotations.json';
        opt.FloatFormat = '%.3f';
        savejson('root', joint_all, opt);

    elseif(strcmp(dataset, 'LEEDS'))
        % in cpp: real scale = param_.target_dist()/meta.scale_self = (41/35)/scale_input
        targetDist = 41/35; % in caffe cpp file 41/35
        oriTrTe = load('../dataset/LEEDS/lsp_dataset/joints.mat');
        extTrain = load('../dataset/LEEDS/lspet_dataset/joints.mat');

        flag_showJoints = 0;
        flag_showScaledImg = 0;

        % in LEEDS:
        % 1  Right ankle
        % 2  Right knee
        % 3  Right hip
        % 4  Left hip
        % 5  Left knee
        % 6  Left ankle
        % 7  Right wrist
        % 8  Right elbow
        % 9  Right shoulder
        % 10 Left shoulder
        % 11 Left elbow
        % 12 Left wrist
        % 13 Neck
        % 14 Head top
        % 15,16 DUMMY
        % We want to comply to MPII: (1 - r ankle, 2 - r knee, 3 - r hip, 4 - l hip, 5 - l knee, 6 - l ankle, ..
        %                             7 - pelvis, 8 - thorax, 9 - upper neck, 10 - head top,
        %                             11 - r wrist, 12 - r elbow, 13 - r shoulder, 14 - l shoulder, 15 - l elbow, 16 - l wrist)
        ordering = [1 2 3, 4 5 6, 15 16, 13 14, 7 8 9, 10 11 12]; % should follow MPI 16 parts..?
        oriTrTe.joints(:,[15 16],:) = 0;
        oriTrTe.joints = oriTrTe.joints(:,ordering,:);
        oriTrTe.joints(3,:,:) = 1 - oriTrTe.joints(3,:,:);
        oriTrTe.joints = permute(oriTrTe.joints, [2 1 3]);

        extTrain.joints([15 16],:,:) = 0;
        extTrain.joints = extTrain.joints(ordering,:,:);

        count = 1;

        %path = {'lspet_dataset/images/im%05d.jpg', 'lsp_dataset/images/im%04d.jpg'};
        path = {'lspet_dataset/images/im%05d.png', 'lsp_dataset/images/im%04d.jpg'};
        %local_path = {'../dataset/LEEDS/lspet_dataset/images/im%05d.jpg', '../dataset/LEEDS/lsp_dataset/images/im%04d.jpg'};
        local_path = {'../dataset/LEEDS/lspet_dataset/images/im%05d.png', '../dataset/LEEDS/lsp_dataset/images/im%04d.jpg'};

        num_image = [10000, 1000];

        for dataset = 1:2
            for im = 1:num_image(dataset)
                path_this = sprintf(local_path{dataset}, im);
                if exist(path_this, 'file') == 0
                    continue
                end
                disp(im);
                disp(count);

                % trivial stuff for LEEDS
                joint_all(count).dataset = 'LEEDS';
                joint_all(count).isValidation = 0;
                joint_all(count).img_paths = sprintf(path{dataset}, im);
                joint_all(count).numOtherPeople = 0;
                joint_all(count).annolist_index = count;
                joint_all(count).people_index = 1;
                % joints and w, h
                if(dataset == 1)
                    %joint_this = extTrain.joints(:,:,im);
                    joint_this = extTrain.joints(:,:,count);
                else
                    joint_this = oriTrTe.joints(:,:,im);
                end

                [h,w,~] = size(imread(path_this));

                joint_all(count).img_width = w;
                joint_all(count).img_height = h;
                joint_all(count).joint_self = joint_this;
                % infer objpos
                invisible = (joint_all(count).joint_self(:,3) == 0);
                if(dataset == 1) %lspet is not tightly cropped
                    joint_all(count).objpos(1) = (min(joint_all(count).joint_self(~invisible, 1)) + max(joint_all(count).joint_self(~invisible, 1))) / 2;
                    joint_all(count).objpos(2) = (min(joint_all(count).joint_self(~invisible, 2)) + max(joint_all(count).joint_self(~invisible, 2))) / 2;
                else
                    joint_all(count).objpos(1) = w/2;
                    joint_all(count).objpos(2) = h/2;
                end
                % visualize
                if (flag_showJoints == 1)
                 figure(1); clf; imshow(path_this);
                 hold on;
                 plot(joint_all(count).joint_self([1 2 3], 1), joint_all(count).joint_self([1 2 3], 2), 'wx', 'Linewidth', 3);
                 plot(joint_all(count).joint_self([4 5 6], 1), joint_all(count).joint_self([4 5 6], 2), 'bx', 'Linewidth', 3);
                 plot(joint_all(count).joint_self([9 10], 1), joint_all(count).joint_self([9 10], 2), 'gx', 'Linewidth', 3);
                 plot(joint_all(count).joint_self([11 12 13], 1), joint_all(count).joint_self([11 12 13], 2), 'mx', 'Linewidth', 3);
                 plot(joint_all(count).joint_self([14 15 16], 1), joint_all(count).joint_self([14 15 16], 2), 'yx', 'Linewidth', 3);
                 plot(joint_all(count).joint_self(invisible, 1), joint_all(count).joint_self(invisible, 2), 'rx');
                 plot(joint_all(count).objpos(1), joint_all(count).objpos(2), 'cs');
                 pause;
                end
                % increase counter and display info
                count = count + 1;
                fprintf('processing %s\n', path_this);
            end
        end

        joint_all = insertMPILikeScale(joint_all, targetDist);

        if (flag_showScaledImg == 1)
         for i = 1:length(joint_all)
             path = ['../dataset/LEEDS/', joint_all(i).img_paths];
             figure(1); clf;
             img = imread(path);
             scale_in_cpp = targetDist/joint_all(i).scale_provided;
             img = imresize(img, scale_in_cpp);
             img = [zeros(400,size(img,2),3); img];
             img = [zeros(size(img,1),400,3), img];
             img = [img; zeros(400,size(img,2),3)];
             img = [img, zeros(size(img,1),400,3)];
             imshow(img);
             hold on;
             fprintf('scale: %f\n', scale_in_cpp);
             objpos = joint_all(i).objpos * scale_in_cpp + 400;
             plot(objpos(1), objpos(2), 'cs');
             line([objpos(1)-368/2, objpos(1)-368/2], [objpos(2)+368/2, objpos(2)-368/2], 'lineWidth', 3);
             line([objpos(1)+368/2, objpos(1)-368/2], [objpos(2)+368/2, objpos(2)+368/2], 'lineWidth', 3);
             line([objpos(1)+368/2, objpos(1)+368/2], [objpos(2)-368/2, objpos(2)+368/2], 'lineWidth', 3);
             line([objpos(1)-368/2, objpos(1)+368/2], [objpos(2)-368/2, objpos(2)-368/2], 'lineWidth', 3);
             pause;
         end
        end

        opt.FileName = 'json/LEEDS_annotations.json';
        opt.FloatFormat = '%.3f';
        savejson('root', joint_all, opt);

    elseif(strcmp(dataset, 'FLIC'))
        % note FLIC is OC
        targetDist = 41/35;
        constant = 0.0110;
        annotation = load('../dataset/FLIC/examples.mat');
        count_flip = 0;
        anno = {annotation.examples.filepath};
        [~,~,IC] = unique(anno);

        count = 1;
        for i = 1:length(annotation.examples)
            if(annotation.examples(i).istest)
                continue;
            end
            fprintf('FLIC images %d, conut = %d\n', i, count);
            joint_all(count).dataset = 'FLIC';
            joint_all(count).isValidation = 0;

            img_path = sprintf('images/%s', annotation.examples(i).filepath);
            joint_all(count).img_paths = img_path;

            current_IC = IC(i);
            fellow_idx = find(IC == current_IC);
            fellow_idx = setdiff(fellow_idx, i);

            joint_all(count).numOtherPeople = length(fellow_idx);
            joint_all(count).annolist_index = i;
            joint_all(count).people_index = 1;

            coords = annotation.examples(i).coords;
            coords(isnan(coords)) = [];
            coords = reshape(coords, [2 11]);
            coords(:,[9 10]) = [];
            coords = coords';
            coords = [coords, ones(9,1)];

            if(coords(4,1) > coords(1,1))
                coords([1 2 3 4 5 6 7 8],:) = coords([4 5 6 1 2 3 8 7],:);
                count_flip = count_flip + 1;
            end

            joint_all(count).joint_self = coords;

            img = imread(sprintf('../dataset/FLIC/%s', img_path));

            for k = 1:length(fellow_idx)
                coords_fellow = annotation.examples(fellow_idx(k)).coords;
                coords_fellow(isnan(coords_fellow)) = [];
                coords_fellow = reshape(coords_fellow, [2 11]);
                coords_fellow(:,[9 10]) = [];
                coords_fellow = coords_fellow';
                coords_fellow = [coords_fellow, ones(9,1)];
                if(coords_fellow(4,1) > coords_fellow(1,1))
                    coords_fellow([1 2 3 4 5 6 7 8],:) = coords_fellow([4 5 6 1 2 3 8 7],:);
                end
                box = annotation.examples(fellow_idx(k)).torsobox;
                torso_height = box(4) - box(2);
                scale = constant*torso_height;

                joint_all(count).joint_others{k} = coords_fellow;
                joint_all(count).scale_provided_other(k) = scale;
                joint_all(count).objpos_other{k} = [(box(1)+box(3))/2, (box(2)*0.7 + box(4)*0.3)];
            end

            joint_all(count).img_width = size(img, 2);
            joint_all(count).img_height = size(img, 1);

%             imshow(img); hold on;
%             plot(coords(:,1), coords(:,2), 'bx');
%             plot(coords(1:3,1), coords(1:3,2), 'rx');
%             plot(coords(end,1), coords(end,2), 'yx');
            box = annotation.examples(i).torsobox;
%             plot([box(1) box(1) box(3) box(3) box(1)], [box(2) box(4) box(4) box(2) box(2)]);
            torso_height = box(4) - box(2);

            joint_all(count).objpos = [(box(1)+box(3))/2, (box(2)*0.7 + box(4)*0.3)];

            % in cpp: real scale = param_.target_dist()/meta.scale_self = (41/35)/scale_input

            scale = constant * torso_height;
%             real_scale = targetDist/scale;
%             img = imresize(img, real_scale);
%
%             x = joint_all(i).objpos(1)*real_scale;
%             y = joint_all(i).objpos(2)*real_scale;
%             if(length(fellow_idx)~=0)
%                 clf; imshow(img); hold on;
%                 plot(joint_all(i).objpos(1)*real_scale, joint_all(i).objpos(2)*real_scale, 'cs');
%                 plot([x+184, x+184, x-184, x-184, x+184], [y+184 y-184 y-184 y+184 y+184], 'b-', 'LineWidth', 3);
%                 plot(coords(:,1)*real_scale, coords(:,2)*real_scale, 'rx');
%                 for k = 1:length(fellow_idx)
%                     plot(joint_all(i).objpos_other(k,1)*real_scale, joint_all(i).objpos_other(k,2)*real_scale, 'cs');
%                     coords = joint_all(i).joint_others(:,:,k);
%                     plot(coords(:,1)*real_scale, coords(:,2)*real_scale, 'yx');
%                 end
%             end

%             lm(i) = min(coords(:,1)*real_scale) - (x-184);
%             rm(i) = (x+184) - max(coords(:,1)*real_scale);
%             tm(i) = min(coords(:,2)*real_scale) - (y-184);
%             bm(i) = (y+184) - max(coords(:,2)*real_scale);

            joint_all(count).scale_provided = scale;
            count = count + 1;
        end
        fprintf('Flipped %d samples\n', count_flip);
        opt.FileName = 'json/FLIC_annotations.json';
        opt.FloatFormat = '%.3f';
        savejson('root', joint_all, opt);
    end

function out = parseFile(in)
    % out is like /media/posenas1//141215/141215_Pose1/Kinect//KINECTNODE1-December-15-2014-13-22-33/color/color-2561794.394.png
    % in should be /media/posenas1/Captures/poseMachineKinect/141215/141215_Pose1/Kinect//KINECTNODE1-December-15-2014-13-22-33/color/color-2561794.394.png
    out = strrep(in, '/media/posenas1/', '/media/posenas1/Captures/poseMachineKinect');

function joint_all = insertMPILikeScale(joint_all, targetDist)
    % calculate scales for each image first
    joints = cat(3, joint_all.joint_self);
    joints([7 8],:,:) = [];
    pa = [2 3 7, 5 4 7, 8 0, 10 11 7, 13 12 7];
    x = permute(joints(:,1,:), [3 1 2]);
    y = permute(joints(:,2,:), [3 1 2]);
    vis = permute(joints(:,3,:), [3 1 2]);
    validLimb = 1:14-1;

    x_diff = x(:, [1:7,9:14]) - x(:, pa([1:7,9:14]));
    y_diff = y(:, [1:7,9:14]) - y(:, pa([1:7,9:14]));
    limb_vis = vis(:, [1:7,9:14]) .* vis(:, pa([1:7,9:14]));
    l = sqrt(x_diff.^2 + y_diff.^2);

    for p = 1:14-1 % for each limb. reference: 7th limb, which is 7 to pa(7) (neck to head)
        valid_compare = limb_vis(:,7) .* limb_vis(:,p);
        ratio = l(valid_compare==1, p) ./ l(valid_compare==1, 7);
        r(p) = median(ratio(~isnan(ratio), 1));
    end

    numFiles = size(x_diff, 1);
    all_scales = zeros(numFiles, 1);

    boxSize = 368;
    psize = 64;
    nSqueezed = 0;

    for file = 1:numFiles %numFiles
        l_update = l(file, validLimb) ./ r(validLimb);
        l_update = l_update(limb_vis(file,:)==1);
        distToObserve = quantile(l_update, 0.75);
        scale_in_lmdb = distToObserve/35; % can't get too small. 35 is a magic number to balance to MPI
        scale_in_cpp = targetDist/scale_in_lmdb; % can't get too large to be cropped

        visibleParts = joints(:, 3, file);
        visibleParts = joints(visibleParts==1, 1:2, file);
        x_range = max(visibleParts(:,1)) - min(visibleParts(:,1));
        y_range = max(visibleParts(:,2)) - min(visibleParts(:,2));
        scale_x_ub = (boxSize - psize)/x_range;
        scale_y_ub = (boxSize - psize)/y_range;

        scale_shrink = min(min(scale_x_ub, scale_y_ub), scale_in_cpp);

        if scale_shrink ~= scale_in_cpp
            nSqueezed = nSqueezed + 1;
            fprintf('img %d: scale = %f %f %f shrink %d\n', file, scale_in_cpp, scale_shrink, min(scale_x_ub, scale_y_ub), nSqueezed);
        else
            fprintf('img %d: scale = %f %f %f\n', file, scale_in_cpp, scale_shrink, min(scale_x_ub, scale_y_ub));
        end

        joint_all(file).scale_provided = targetDist/scale_shrink; % back to lmdb unit
    end

    fprintf('total %d squeezed!\n', nSqueezed);

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function [ bic ] = ml_gmm_bic(Data,Priors,Mu,Sigma,cov_type)
%ML_GMM_BIC Bayesian Information criterion
%
%  input ------------------------------------------------------------------
%
%   o Data:     D x N array representing N datapoints of D dimensions.
%
%   o Priors:   1 x K array representing the prior probabilities of the
%               K GMM components.
%   o Mu:       D x K array representing the centers of the K GMM components.
%
%   o Sigma:    D x D x K array representing the covariance matrices of the
%               K GMM components.
%
%   o cov_type: string, covariance type = 'full','diag' or 'iso'
%
%  output -----------------------------------------------------------------
%
%   o bic   : (1 x 1)
%
%

[D,N]       = size(Data);
K           = length(Priors);
num_param   = K-1 + K * D;

if strcmp(cov_type,'full') == true
    num_param = num_param + K * ( D * ( D - 1)/2 );
elseif strcmp(cov_type,'diag') == true
    num_param = num_param + K * D;
elseif strcmp(cov_type,'iso') == true
    num_param = num_param + K * 1;
else
   error(['no such covariance type: ' cov_type '  only full | diag | isot ']);
end

% compute the loglikelihood of the data given the model
% loglik: (N x 1)

loglik = ml_LogLikelihood_gmm(Data,Priors,Mu,Sigma);
bic    = - 2 * loglik + num_param * log(N);

end

author: epfl-lasa


text:
function func = fcond(c, varargin)
% FUNC = FCOND(COND-FUNC, FUNCTION_1, ...)
%   Make a function FUNC that takes any number of ARGS, and invokes
%   FUNCTO_n with them, where n = COND-FUNC(ARGS).
%

  functions = varargin;
  function varargout = fcondFunc(varargin)
    idx = c(varargin{:});
    f = functions{idx};
    [varargout{1:nargout}] = f(varargin{:});
  end

  func = @fcondFunc;
end

author: Sable

repo: mcbench-benchmarks


text:
function [ y, next ] = p24_equil ( neqn, next )

%*****************************************************************************80
%
%% P24_EQUIL returns equilibrium solutions of problem p24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
function r = XYRectangle(x1,x2,y1,y2,z)

% XYRECTANGLE  Rectangle in the XY plane.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

s1 = xSegment(x1,x2,y1,z);
s2 = ySegment(x2,y1,y2,z);
s3 = xSegment(x2,x1,y2,z);
s4 = ySegment(x1,y2,y1,z);

r = [s1 s2 s3 s4];

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [b,c,d,x,y,z,qfac] = vox2rasToQform(vox2ras)
% [a,b,c,x,y,z,qfac] = vox2rasToQform(vox2ras)
%
% Converts a vox2ras matrix to NIFTI qform parameters.
% Note: the vox2ras should be 6 DOF. This code mostly just
% follows CH's mriToNiftiQform() in mriio.c.
%
%  hdr.pixdim(1) = qfac;
%  hdr.quatern_b = b;
%  hdr.quatern_c = c;
%  hdr.quatern_d = d;
%  hdr.qoffset_x = x;
%  hdr.qoffset_y = y;
%  hdr.qoffset_z = z;
%  hdr.qform_code = NIFTI_XFORM_SCANNER_ANAT=1;
%

%
% vox2rasToQform.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

a = [];
if(nargin ~= 1)
  fprintf('[a,b,c,x,y,z,qfac] = vox2rasToQform(vox2ras)\n');
  return;
end

x = vox2ras(1,4);
y = vox2ras(2,4);
z = vox2ras(3,4);

d = sqrt(sum(vox2ras(:,1:3).^2));
Mdc = vox2ras(1:3,1:3) ./ repmat(d,[3 1]);
if(det(Mdc) == 0.0)
  fprintf('ERROR: vox2ras determinant is 0\n');
  return;
end

r11 = Mdc(1,1);
r21 = Mdc(2,1);
r31 = Mdc(3,1);
r12 = Mdc(1,2);
r22 = Mdc(2,2);
r32 = Mdc(3,2);
r13 = Mdc(1,3);
r23 = Mdc(2,3);
r33 = Mdc(3,3);

if(det(Mdc) > 0.0) qfac = 1.0;
else
  r13 = -r13;
  r23 = -r23;
  r33 = -r33;
  qfac = -1.0;
end

%  /* following mat44_to_quatern() */

a = r11 + r22 + r33 + 1.0;
if(a > 0.5)
  a = 0.5 * sqrt(a);
  b = 0.25 * (r32-r23) / a;
  c = 0.25 * (r13-r31) / a;
  d = 0.25 * (r21-r12) / a;
else
  xd = 1.0 + r11 - (r22+r33);
  yd = 1.0 + r22 - (r11+r33);
  zd = 1.0 + r33 - (r11+r22);
  if(xd > 1.0)
    b = 0.5 * sqrt(xd);
    c = 0.25 * (r12+r21) / b;
    d = 0.25 * (r13+r31) / b;
    a = 0.25 * (r32-r23) / b;
  elseif( yd > 1.0 )
    c = 0.5 * sqrt(yd);
    b = 0.25 * (r12+r21) / c;
    d = 0.25 * (r23+r32) / c;
    a = 0.25 * (r13-r31) / c;
  else
    d = 0.5 * sqrt(zd);
    b = 0.25 * (r13+r31) / d;
    c = 0.25 * (r23+r32) / d;
    a = 0.25 * (r21-r12) / d;
  end
  if(a < 0.0)
    a = -a;
    b = -b;
    c = -c;
    d = -d;
  end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [] = fig2idtf(filename,...
                        surface_vertices, faces, face_vertex_data,...
                        line_vertices, line_edges, line_colors,...
                        pointset_points, pointset_colors)
%FIG2IDTF   Save figure in IDTF format.
%
% usage
%   count = FIG2IDTF(fid, surface_vertices, faces, face_vertex_data,...
%                         line_vertices, line_edges, line_colors,...
%                         pointset_points, pointset_colors)
%
% input
%   filename = string of filename (including extension)
%
%   faces = {1 x #surfaces}
%   surface_vertices = {1 x #surfaceplots}
%   face_vertex_data = {1 x #surfaceplots}
%
%   line_vertices = {1 x #lineseries}
%   line_lines = {1 x #lineseries}
%   line_colors = {1 x #lineseries}
%
%   pointset_points = {1 x #pointsets}
%   pointset_colors = {1 x #pointsets}
%
% output
%   count = number of lines written to file with handle fid
%
% reference
%   IDTF (Intermediate Data Text File) Format Description, Version 100,
%   Intel Corporation, 2005, available at:
%       http://u3d.svn.sourceforge.net/viewvc/u3d/releases/Gold12Update/Docs/IntermediateFormat/IDTF%20Format%20Description.pdf
%
% See also FIG2U3D, IDTF2U3D, FIG2PDF3D, U3D_PRE_SURFACE, U3D_PRE_LINE,
%          U3D_PRE_QUIVERGROUP, U3D_PRE_CONTOURGROUP.
%
% File:      fig2idtf.m
% Author:    Ioannis Filippidis, jfilippidis@gmail.com
% Date:      2012.06.10 - 2012.06.14
% Language:  MATLAB R2012a
% Purpose:   export to u3d multiple surfaces, lines, quivergroups
% Copyright: Ioannis Filippidis, 2012-
%
% acknowledgment
%   Based on save_idtf by Alexandre Gramfort.
%   This can be found on the MATLAB Central File Exchange:
%       http://www.mathworks.com/matlabcentral/fileexchange/25383-matlab-mesh-to-pdf-with-3d-interactive-object
%   and is covered by the BSD License.

% depends
%   idtf_model_nodes, populate_mesh_resource_str,
%   populate_line_resource_str, populate_point_resource_str,
%   shaders_materials_modifiers, check_file_extension, verbatim

n_meshes = size(faces, 2);
n_lines = size(line_vertices, 2);
n_pointsets = size(pointset_points, 2);

%% nodes
nodes = idtf_model_nodes(n_meshes, n_lines, n_pointsets);

%% resources
mesh_resources = populate_mesh_resource_str(faces, surface_vertices, face_vertex_data);
n_resources = n_meshes;
line_resources = populate_line_resource_str(line_vertices, line_edges, line_colors, n_resources);
n_resources = n_resources +n_lines;
point_resources = populate_point_resource_str(pointset_points, pointset_colors, n_resources);

mesh_line_resources = [mesh_resources, line_resources, point_resources];

resource_list = resource_list_model_str;
total_resource_number = n_meshes +n_lines +n_pointsets;
resources = sprintf(resource_list, total_resource_number,...
                             mesh_line_resources);

%% final
file_info = file_info_str;
[shaders, materials, modifiers] = shaders_materials_modifiers(surface_vertices, faces, face_vertex_data);

str = [file_info, nodes, resources, shaders, materials, modifiers];

idtffile = check_file_extension(filename, '.idtf');
fid = fopen(idtffile, 'w');
count = fprintf(fid, str);
fclose(fid);

disp(['Number of lines written to IDTF file: ', num2str(count) ] )

function [str] = file_info_str
str = verbatim;
%{
FILE_FORMAT "IDTF"
FORMAT_VERSION 100

%}

function [str] = resource_list_model_str
% total_node_number, aggregate_string_of_resources
str = verbatim;
%{

RESOURCE_LIST "MODEL" {
     RESOURCE_COUNT %d
     %s
}

%}


text:
% PAC_CONT - compute phase-amplitude coupling (power of first input
%         correlation with phase of second). There is no graphical output
%         to this function.
%
% Usage:
%   >> pac_cont(x,y,srate);
%   >> [pac timesout pvals] = pac_cont(x,y,srate,'key1', 'val1', 'key2', val2' ...);
%
% Inputs:
%    x       = [float array] 1-D data vector of size (1xtimes)
%    y       = [float array] 1-D data vector of size (1xtimes)
%    srate   = data sampling rate (Hz)
%
% Optional time information inputs:
%   'winsize'   = If cycles==0: data subwindow length (fastest, 2^n<frames);
%                 If cycles >0: *longest* window length to use. This
%                 determines the lowest output frequency. Note that this
%                 parameter is overwritten if the minimum frequency has been set
%                 manually and requires a longer time window {~frames/8}
%   'ntimesout' = Number of output times (int<frames-winframes). Enter a
%                 negative value [-S] to subsample original time by S.
%   'timesout'  = Enter an array to obtain spectral decomposition at
%                 specific time values (note: algorithm find closest time
%                 point in data and this might result in an unevenly spaced
%                 time array). Overwrite 'ntimesout'. {def: automatic}
%   'tlimits'   = [min max] time limits in ms.
%
% Optional PAC inputs:
%   'method'    = ['modulation'|'plv'|'corr'|'glm'] (see reference).
%   'freqphase' = [min max] frequency limits. Default [minfreq 50],
%                 minfreq being determined by the number of data points,
%                 cycles and sampling frequency. Use 0 for minimum frequency
%                 to compute default minfreq. You may also enter an
%                 array of frequencies for the spectral decomposition
%                 (for FFT, closest computed frequency will be returned; use
%                 'padratio' to change FFT freq. resolution).
%   'freqamp'   = [float array] array of frequencies for the second
%                 argument. 'freqs' is used for the first argument.
%                 By default it is the same as 'freqs'.
%   'filterfunc' = ['eegfilt'|'iirfilt'|'eegfftfilt'] filtering function.
%                 Default is iirfilt. Warning, filtering may dramatically
%                 affect the result. With the 'corr' method, make sure you
%                 have a large window size because each window is filtered
%                 independently.
%   'filterphase' = @f_handle. Function handle to filter the data for the
%                 phase information. For example, @(x)iirfilt(x, 1000, 2,
%                 20). Note that 'freqphase' is ignore in this case.
%   'filteramp' = @f_handle. Function handle to filter the data for the
%                 amplitude information. Note that 'freqamp' is ignore in
%                 this case.
%
% Inputs for statistics:
%   'alpha'     = [float] p-value threshold. Default is none (no statistics).
%   'mcorrect'  = ['none'|'fdr'] method to correct for multiple comparison.
%                 Default is 'none'.
%   'baseline'  = [min max] baseline period for the Null distribution. Default
%                 is the whole data range. Note that this option is ignored
%                 for instanstaneous statistics.
%   'instantstat' = ['on'|'off'] performs statistics for each time window
%                 independently. Default is 'off'.
%   'naccu'     = [integer] number of accumulations for surrogate
%                 statistics.
%   'statlim'   = ['parametric'|'surrogate'] use a parametric methods to
%                 asseess the limit of the surrogate distribution or use
%                 the tail of the distribution ('surrogate' method)
%
% Other inputs:
%   'title'     = [string] figure title. Default is none.
%   'vert'      = [float array] array of time value for which to plot
%                 vertical lines. Default is none.
%
% Outputs:
%  pac      = Phase-amplitude coupling values.
%  timesout = vector of time indices
%  pvals    = Associated p-values
%
% Author: Arnaud Delorme and Makoto Miyakoshi, SCCN/INC, UCSD 2012-
%
% References:
% Methods used here are introduced and compared in:
% Penny, Duzel, Miller, Ojemann. (20089). Testing for Nested Oscilations.
% J Neuro Methods. 174:50-61
%
% Modulation Index is defined in:
% Canolty, Edwards, Dalal, Soltani, Nagarajan, Kirsch, et al. (2006). Modulation index is defined in High Gamma Power Is Phase-Locked to Theta
% Oscillations in Human Neocortex. Science. 313:1626-8.
%
% PLV (Phase locking value) is defined in:
% Lachaux, Rodriguez, Martiniere, Varela. (1999). Measuring phase synchrony
% in brain signal. Hum Brain Mapp. 8:194-208.
%
% corr (correlation) method is defined in:
% Brunce, Eckhorn. (2004). Task-related coupling from high- to
% low-frequency signals among visual cortical areas in human subdural
% recordings. Int J Psychophysiol. 51:97-116.
%
% glm (general linear model) is defined in
% Penny, Duzel, Miller, Ojemann. (20089). Testing for Nested Oscilations.
% J Neuro Methods. 174:50-61

% Copyright (C) 2012 Arnaud Delorme, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [m_raw pvals indexout] = pac_cont(X, Y, srate, varargin);

if nargin < 1
    help pac_cont;
    return;
end

% deal with 3-D inputs
% --------------------
if ndims(X) == 3 || ndims(Y) == 3, error('Cannot process 3-D input'); end
if size(X,1) > 1, X = X'; end
if size(Y,1) > 1, Y = Y'; end
if size(X,1) ~= 1 || size(Y,1) ~= 1, error('Cannot only process vector input'); end
frame = size(X,2);
pvals = [];

g = finputcheck(varargin, ...
    { ...
    'alpha'         'real'     [0 0.2]                [];
    'baseline'      'float'    []                       [];
    'freqphase'     'real'     [0 Inf]                  [0 srate/2];
    'freqamp'       'real'     [0 Inf]                  [];
    'mcorrect'      'string'   { 'none' 'fdr' }         'none';
    'method'        'string'   { 'plv' 'modulation' 'glm' 'corr' } 'modulation';
    'naccu'         'integer'  [1 Inf]                   250;
    'instantstat'   'string'   {'on','off'}              'off';
    'newfig'        'string'   {'on','off'}              'on';
    'nofig'         'string'   {'on','off'}              'off';
    'statlim'       'string'   {'surrogate','parametric'}  'parametric';
    'timesout'      'real'     []                        []; ...
    'filterfunc'    'string'   { 'eegfilt' 'iirfilt' 'eegfiltfft' }   'eegfiltfft'; ...
    'filterphase'   ''         {}                        [];
    'filteramp'     ''         {}                        [];
    'ntimesout'     'integer'  []                        200; ...
    'tlimits'       'real'     []                        [0 frame/srate];
    'title'         'string'   []                        '';
    'vert'          'real'     []                    [];
    'winsize'       'integer'  [0 Inf]                   max(pow2(nextpow2(frame)-3),4) }, 'pac');

if ischar(g), error(g); end

if ~isempty(g.filterphase)
     x_freqphase = feval(g.filterphase, X(:)');
else x_freqphase = feval(g.filterfunc, X(:)', srate, g.freqphase(1), g.freqphase(end));
end
if ~isempty(g.filteramp)
     x_freqamp   = feval(g.filteramp, Y(:)');
else x_freqamp   = feval(g.filterfunc, Y(:)', srate, g.freqamp(  1), g.freqamp( end));
end
z_phasedata = hilbert(x_freqphase);
z_ampdata   = hilbert(x_freqamp);
phase       = angle(z_phasedata);
amplitude   = abs(  z_ampdata);
z           = amplitude.*exp(i*phase); % this is the pac measure

% get time windows
% ----------------
g.verbose = 'on';
g.causal  = 'off';
[ timesout1 indexout ] = gettimes(frame, g.tlimits, g.timesout, g.winsize, g.ntimesout, g.causal, g.verbose);

% scan time windows
% -----------------
if ~isempty(g.alpha)
    m_raw = zeros(1,length(indexout));
    pvals = zeros(1,length(indexout));
end
fprintf('Computing PAC:\n');
for iWin = 1:length(indexout)
    x_phaseEpoch = x_freqphase(indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
    x_ampEpoch   = x_freqamp(  indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
    z_phaseEpoch = z_phasedata(indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
    z_ampEpoch   = z_ampdata(  indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
    z_epoch      = z(          indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);

    numpoints=length(x_phaseEpoch);
    if rem(iWin,10) == 0,  verboseprintf(g.verbose, ' %d',iWin); end
    if rem(iWin,120) == 0, verboseprintf(g.verbose, '\n'); end

    % Choose method
    % -------------
    if strcmpi(g.method, 'modulation')

        % Modulation index
        m_raw(iWin) = abs(sum(z_epoch))/numpoints;

    elseif strcmpi(g.method, 'plv')

        if iWin == 145
            %dsfsd;
        end

        %amplitude_filt = sgolayfilt(amplitude, 3, 101);
        if ~isempty(g.filterphase)
             amplitude_filt = feval(g.filterphase, z_ampEpoch);
        else amplitude_filt = feval(g.filterfunc , z_ampEpoch, srate, g.freqphase(1), g.freqphase(end));
        end
        z_amplitude_filt = hilbert(amplitude_filt);

        phase_amp_modulation = angle(z_amplitude_filt);
        m_raw(iWin) = abs(sum(exp(i*(x_phaseEpoch - phase_amp_modulation)))/numpoints);

    elseif strcmpi(g.method, 'corr')

        if iWin == inf %145
            figure; plot(abs(z_ampdata))
            hold on; plot(x_phasedata/10000000000, 'r')
            x = X(indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
            hold on; plot(x, 'g');
            dsfsd;
        end
        [r_ESC pval_corr] = corrcoef(x_phaseEpoch, abs(z_ampEpoch));
        m_raw(iWin)   = r_ESC(1,2);
        pvals(iWin)   = pval_corr(1,2);

    elseif strcmpi(g.method, 'glm')

        [b dev stats] = glmfit(x_phaseEpoch', abs(z_ampEpoch)', 'normal');
        GLM_beta      = stats.beta(2,1);
        pvals(iWin)   = stats.p(2,1);
        m_raw(iWin)   = b(1);

    end

    %% compute statistics (instantaneous)
    % -----------------------------------
    if ~isempty(g.alpha) && strcmpi(g.instantstat, 'on') && ~strcmpi(g.method, 'corr') && ~strcmpi(g.method, 'glm')

        % compute surrogate values
        numsurrogate=g.naccu;
        minskip=srate;
        maxskip=numpoints-srate; % max variation half a second
        if maxskip < 1
            error('Window size shorter than 1 second; too short for computing surrogate data');
        end
        skip=ceil(numpoints.*rand(numsurrogate*4,1));
        skip(skip>maxskip)=[];
        skip(skip<minskip)=[];
        skip=skip(1:numsurrogate,1);
        surrogate_m=zeros(numsurrogate,1);
        for s=1:numsurrogate
            surrogate_amplitude=[amplitude(skip(s):end) amplitude(1:skip(s)-1)]; % consider circular shifts
            surrogate_m(s)=abs(mean(surrogate_amplitude.*exp(i*phase)));
            %disp(numsurrogate-s)
        end

        if strcmpi(g.statlim, 'surrogate')
            pvals(iWin) = stat_surrogate_pvals(surrogate_m, m_raw(iWin), 'upper');
            %fprintf('Raw PAC is %3.2f (p-value=%1.3f)\n', m_raw(iWin), pvals(iWin));

        else
            % Canolty method below

            %% fit gaussian to surrogate data, uses normfit.m from MATLAB Statistics toolbox
            [surrogate_mean,surrogate_std]=normfit(surrogate_m);

            %% normalize length using surrogate data (z-score)
            m_norm_length=(abs(m_raw(iWin))-surrogate_mean)/surrogate_std;
            pvals(iWin)  = normcdf(0, m_norm_length, 1);
            m_norm_phase=angle(m_raw(iWin));
            m_norm=m_norm_length*exp(i*m_norm_phase);

            % compare parametric and non-parametric methods (return similar
            % results)
            if iWin == length(indexout)
                figure;
                plot(-log10(pvals));
                hold on; plot(-log10(pvals2), 'r');
            end
        end
    end

end
fprintf('\n');

% Computes alpha
% --------------
if ~isempty(g.alpha) && strcmpi(g.instantstat, 'off')
    if isempty(g.baseline)
        g.baseline = [ timesout1(1) timesout1(end) ];
    end
    baselineInd = find(timesout1 >= g.baseline(1) & timesout1 <= g.baseline(end));

    m_raw_base = abs(m_raw(baselineInd));

    if strcmpi(g.statlim, 'surrogate')
        for index = 1:length(m_raw)
            pvals(index) = stat_surrogate_pvals(m_raw_base, m_raw(index), 'upper');
        end
    else
        [surrogate_mean,surrogate_std]=normfit(m_raw_base);
        m_norm_length=(abs(m_raw)-surrogate_mean)/surrogate_std;
        pvals = normcdf(0, m_norm_length, 1);
    end
    if strcmpi(g.mcorrect, 'fdr')
        pvals = fdr(pvals);
    end
end

%% plot results
% -------------
if strcmpi(g.nofig, 'on')
    return
end
if strcmpi(g.newfig, 'on')
    figure;
end
if ~isempty(g.alpha)
    plotcurve(timesout1, m_raw, 'maskarray', pvals < g.alpha);
else plotcurve(timesout1, m_raw);
end
xlabel('Time (ms)');
ylabel('PAC (0 to 1)');
title(g.title);

% plot vertical lines
% -------------------
if ~isempty(g.vert)
    hold on;
    yl = ylim;
    for index = 1:length(g.vert)
        plot([g.vert(index) g.vert(index)], yl, 'g');
    end
end

% -------------
% gettime function identical to timefreq function
% DO NOT MODIFY
% -------------
function [ timevals, timeindices ] = gettimes(frames, tlimits, timevar, winsize, ntimevar, causal, verbose);
timevect = linspace(tlimits(1), tlimits(2), frames);
srate    = 1000*(frames-1)/(tlimits(2)-tlimits(1));

if isempty(timevar) % no pre-defined time points
    if ntimevar(1) > 0
        % generate linearly space vector
        % ------------------------------
        if (ntimevar > frames-winsize)
            ntimevar = frames-winsize;
            if ntimevar < 0
                error('Not enough data points, reduce the window size or lowest frequency');
            end
            verboseprintf(verbose, ['Value of ''timesout'' must be <= frame-winsize, ''timesout'' adjusted to ' int2str(ntimevar) '\n']);
        end
        npoints = ntimevar(1);
        wintime = 500*winsize/srate;
        if strcmpi(causal, 'on')
             timevals = linspace(tlimits(1)+2*wintime, tlimits(2), npoints);
        else timevals = linspace(tlimits(1)+wintime, tlimits(2)-wintime, npoints);
        end
        verboseprintf(verbose, 'Generating %d time points (%1.1f to %1.1f ms)\n', npoints, min(timevals), max(timevals));
    else
        % subsample data
        % --------------
        nsub     = -ntimevar(1);
        if strcmpi(causal, 'on')
            timeindices = [ceil(winsize+nsub):nsub:length(timevect)];
        else timeindices = [ceil(winsize/2+nsub/2):nsub:length(timevect)-ceil(winsize/2)-1];
        end
        timevals    = timevect( timeindices ); % the conversion at line 741 leaves timeindices unchanged
        verboseprintf(verbose, 'Subsampling by %d (%1.1f to %1.1f ms)\n', nsub, min(timevals), max(timevals));
    end
else
    timevals = timevar;
    % check boundaries
    % ----------------
    wintime = 500*winsize/srate;
    if strcmpi(causal, 'on')
         tmpind  = find( (timevals >= tlimits(1)+2*wintime-0.0001) & (timevals <= tlimits(2)) );
    else tmpind  = find( (timevals >= tlimits(1)+wintime-0.0001) & (timevals <= tlimits(2)-wintime+0.0001) );
    end
    % 0.0001 account for numerical innacuracies on opteron computers
    if isempty(tmpind)
        error('No time points. Reduce time window or minimum frequency.');
    end
    if  length(timevals) ~= length(tmpind)
        verboseprintf(verbose, 'Warning: %d out of %d time values were removed (now %3.2f to %3.2f ms) so the lowest\n', ...
            length(timevals)-length(tmpind), length(timevals), timevals(tmpind(1)), timevals(tmpind(end)));
        verboseprintf(verbose, '         frequency could be computed with the requested accuracy\n');
    end
    timevals = timevals(tmpind);
end

% find closet points in data
% --------------------------
timeindices = round(eeg_lat2point(timevals, 1, srate, tlimits, 1E-3));
if length(timeindices) < length(unique(timeindices))
    timeindices = unique_bc(timeindices)
    verboseprintf(verbose, 'Warning: duplicate times, reduce the number of output times\n');
end
if length(unique(timeindices(2:end)-timeindices(1:end-1))) > 1
    verboseprintf(verbose, 'Finding closest points for time variable\n');
    verboseprintf(verbose, 'Time values for time/freq decomposition is not perfectly uniformly distributed\n');
else
    verboseprintf(verbose, 'Distribution of data point for time/freq decomposition is perfectly uniform\n');
end
timevals    = timevect(timeindices);

function verboseprintf(verbose, varargin)
if strcmpi(verbose, 'on')
    fprintf(varargin{:});
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%% change gaussian width
tmp_gsiz = str2double(get(edit_gsiz, 'string'));
neuron_raw.options.gSiz = tmp_gsiz;
if ~exist('neuron', 'var')
    neuron.options.gsiz = tmp_gsiz / neuron_raw.options.ssub;
end

author: zhoupc


text:
function pdf = point_distance_3d_pdf ( x, a, b )

%*****************************************************************************80
%
%% POINT_DISTANCE_3D_PDF evaluates the point distance PDF in the 3D.
%
%  Discussion:
%
%    It is assumed that a set of points has been generated in 3D
%    according to a Poisson process.  The number of points in a region
%    of size VOLUME is a Poisson variate with mean value B * VOLUME.
%
%    For a point chosen at random, we may now find the nearest
%    Poisson point, the second nearest and so on.  We are interested
%    in the PDF that governs the expected behavior of the distances
%    of rank A = 1, 2, 3, ... with Poisson density B.
%
%  Formula:
%
%    PDF(X)(A,B) = 3 * ( (4/3) * B * PI )**A * X**( 3 * A - 1 )
%      * EXP ( - (4/3) * B * PI * X * X * X ) / ( A - 1 )!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996, pages 580.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    0.0 <= X.
%
%    Input, integer A, indicates the degree of nearness of the point.
%    A = 1 means the nearest point, A = 2 the second nearest, and so on.
%    0 < A.
%
%    Input, real B, the Poisson point density.  0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( a < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POINT_DISTANCE_3D_PDF - Fatal error!\n' );
    fprintf ( 1, '  Input parameter A < 1.\n' );
    error ( 'POINT_DISTANCE_3D_PDF - Fatal error!' );
  end

  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POINT_DISTANCE_3D_PDF - Fatal error!\n' );
    fprintf ( 1, '  Input parameter B <= 0.0.\n' );
    error ( 'POINT_DISTANCE_3D_PDF - Fatal error!' );
  end

  if ( x < 0.0 )
    pdf = 0.0;
  else
    pdf = 3.0 * ( ( 4.0 / 3.0 ) * b * pi )^a ...
      * x^( 3 * a - 1 ) * exp ( - ( 4.0 / 3.0 ) * b * pi * x^3 ) ...
      / i4_factorial ( a - 1 );
  end

  return
end


text:
function pass = test_times(pref)

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

dom = [-1, 1];
diffOp = operatorBlock.diff(dom);
V = chebpoly(1:6);
A = linop(diffOp);
AV = A*V;

B = linop(2*diffOp);
BV = B*V;

err(1) = norm(2*AV - BV);
tol = 1e-14;
pass(1) = err < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = processingLLCs(process, varargin)
% This function contains two different processes for manipulating the
% localized loop constraints (LLCs) to formulate a model-specific and
% objective-function-specific MILP for finding loopless flux distributions
% with a minimal number of binary variables
%
% USAGE:
%    1. Preprocess loop information after calling `addLoopLawConstraints`:
%       [solveLP, MILPproblem, loopInfo] = processingLLCs('preprocess', loopInfo, LPproblem, model, nRxns, osenseStr, MILPproblem)
%
%    2. Update the loop constraints for a specific objective vector
%       [solveLP, MILPproblem] = processingLLCs('update', loopInfo, osenseStr, MILPproblem, objVector)
%
% INPUTS:
%    loopInfo:      structure containing info about the loops, initially outputed
%                   from `addLoopLawConstraints` and updated by the current function
%    LPproblem:     original COBRA LP problem structure
%    model:         the COBRA model from which LPproblem is constructed
%    nRxns:         number of reactions in the model (default size(model.S, 2))
%    osenseStr:     optimization sense of the current problem (e.g., FBA max v_biomass, or FVA min v_PYK)
%                   'max' (defaulted) or 'min'
%    MILPproblem:   COBRA MILP problem generated from LPproblem using `addLoopLawConstraints`
%                   (default generated from LPproblem if not given in the preprocessing call)
%    objVector:     nRxns-by-1 objective vector for the current optimization problem
%                   for updating LLCs (e.g., FVA min v_PYK given fixed v_biomass)
%
% OUTPUTS:
%    solveLP:       true if solving LP is sufficient to gaurantee the objective function value
%                   is the same as solving MILP with loop constraints
%    MILPproblem:   updated MILP problem with LLCs
%    loopInfo:      updated loopInfo with info about the loops (outputed only for the preprocess call)

solveLP = true;
switch process
    case 'preprocess'
        loopInfo = preprocessLLCs(varargin{1:5});
        MILPproblem = varargin{6};
        % initial update for MILP. No specific objective reactions
        if loopInfo.alwaysLLC
            % need to solve MILP if the problem constraints or original objective function necessitate the need
            solveLP = false;
            % update bounds and rhs
            MILPproblem = updateLLCs(MILPproblem, loopInfo, []);
        end
        varargout = {solveLP, MILPproblem, loopInfo};
    case 'update'
        loopInfo = varargin{1};
        osenseStr = varargin{2};
        if isempty(osenseStr)
            osenseStr = 'max';
        end
        MILPproblem = varargin{3};
        if numel(varargin) < 4
            objVector = [];
        else
            objVector = varargin{4};
        end

        if isempty(objVector)
            % if updating but no object vector is given, just restore the original bounds
            MILPproblem = restoreOriginalBounds(MILPproblem, loopInfo.rhs0, loopInfo.var, loopInfo.BDg);
        else
            % update with a specific objective function
            [rxnID, ~, cCoeff] = find(objVector(:));
            osense = strcmp(osenseStr, 'min') - strcmp(osenseStr, 'max');
            % need to solve MILP if the problem constraints or original objective function necessitate the need or
            % the reactions being minimized has their reverse direction in loops or
            % the reactions being maximized has their forward direction in loops
            if loopInfo.alwaysLLC || any(loopInfo.rxnInLoops(rxnID(cCoeff * osense > 0), 1)) || any(loopInfo.rxnInLoops(rxnID(cCoeff * osense < 0), 2))
                solveLP = false;
                % restore the original bounds
                MILPproblem = restoreOriginalBounds(MILPproblem, loopInfo.rhs0, loopInfo.var, loopInfo.BDg);
                % find the reactions in the objective function which are being minimized and has their reverse
                % direction in loops or the reactions being maximized has their forward direction in loops
                rxnIdLLC = false(numel(rxnID));
                for j = 1:numel(rxnID)
                    if (loopInfo.rxnInLoops(rxnID(j), 1) && cCoeff(j) * osense > 0) ...
                            || (loopInfo.rxnInLoops(rxnID(j), 2) && cCoeff(j) * osense < 0)
                        rxnIdLLC(j) = true;
                    end
                end
                rxnIdLLC = rxnID(rxnIdLLC);
                % update bounds and rhs
                MILPproblem = updateLLCs(MILPproblem, loopInfo, rxnIdLLC);
            end

        end
        varargout = {solveLP, MILPproblem};
end

end

function loopInfo = preprocessLLCs(loopInfo, LPproblem, model, nRxns, osenseStr)
if nargin < 4 || isempty(nRxns)
    nRxns = size(model.S, 2);
end
if nargin < 5 || isempty(osenseStr)
    osenseStr = 'max';
end
if strcmp(osenseStr, 'min')
    model.c = -model.c;
end
rxnInLoops = loopInfo.rxnInLoops;
conComp  = loopInfo.conComp;

% determine the set of reactions for which LLCs are always required
% condition I in Prop. 2 in Chan et al., 2017
cond1 = rxnInLoops(:, 2) & model.c > 0;
% condition II in the paper in Prop. 2 in Chan et al., 2017
cond2 = rxnInLoops(:, 1) & model.c < 0;
% condition III in the paper in Prop. 2 in Chan et al., 2017
[cond3A1, cond3A2, cond3B] = deal(false(nRxns, 1));
for i = (size(model.S, 1) + 1):size(LPproblem.A, 1)
    % for constraint p with sum(a_pj * v_j) <= b_p
    if ~strcmp(LPproblem.csense(i), 'G')  % '<=' or '=' constraint
        % if reaction j has its forward direction in cycles and a_pj < 0
        cond3A1 = cond3A1 | (rxnInLoops(:, 2) & LPproblem.A(i, 1:nRxns)' < 0);
        % if reaction j has its reverse direction in cycles and a_pj > 0
        cond3A2 = cond3A2 | (rxnInLoops(:, 1) & LPproblem.A(i, 1:nRxns)' > 0);
        % if the constraint involves 2 or more reactinos or RHS < 0
        cond3B = cond3B | (nnz(LPproblem.A(i, 1:nRxns)) > 1 | LPproblem.b(i) < 0);
    end
    if ~strcmp(LPproblem.csense(i), 'L')  % '>=' or '=' constraint
        cond3A1 = cond3A1 | (rxnInLoops(:, 2) & LPproblem.A(i, 1:nRxns)' > 0);
        cond3A2 = con3A2 | (rxnInLoops(:, 1) & LPproblem.A(i, 1:nRxns)' < 0);
        cond3B = cond3B | (nnz(LPproblem.A(i, 1:nRxns)) > 1 | LPproblem.b(i) > 0);
    end
end
% reactions satisfying (3A1 or 3A2) and 3B
cond3 = (cond3A1 | cond3A2) & cond3B;
% condition III for bound constraints can be simplified as follows:
cond3 = cond3 | (model.lb > 0 & rxnInLoops(:, 2)) | (model.ub < 0 & rxnInLoops(:, 1));
% reactions that are required to be constrained by loopless constraints all the time
rxnInLoopsAlwaysOn = cond1 | cond2 | cond3;
% LLCs are always required if the set is non-empty
alwaysLLC = any(rxnInLoopsAlwaysOn);
% the corresponding set of reactions in the same connected components as
% the always-on reactions
conCompAlwaysOn = false(max(conComp), 1);
conCompAlwaysOn(conComp(rxnInLoopsAlwaysOn)) = true;
if loopInfo.printLevel
    fprintf('Reactions in internal nullspace can be divided into %d connected components.\n', max(conComp))
end

% get an initial feasible and loopless solution in case MipStart is needed
model2 = model;
model2.lb = model2.lb(1:size(model2.S, 2));
model2.ub = model2.ub(1:size(model2.S, 2));
model2.c = zeros(size(model2.S, 2), 1);
model2.b = zeros(size(model2.S, 1), 1);
sFeas = optimizeCbModel(model2, 'max', 'one');
x0 = sFeas.x;
[loopInfo.alwaysLLC, loopInfo.rxnInLoopsAlwaysOn, loopInfo.conCompAlwaysOn, loopInfo.x0] ...
    = deal(alwaysLLC, rxnInLoopsAlwaysOn, conCompAlwaysOn, x0);
end

function MILPproblemLLC = updateLLCs(MILPproblemLLC, loopInfo, rxnID)
% apply LLCs by relaxing constraints and pre-assign values to variables
if nargin < 3
    rxnID = [];
end
conCompOn = loopInfo.conCompAlwaysOn;
conCompOn(loopInfo.conComp(rxnID)) = true;

bigM = inf;
if ~loopInfo.useRxnLink
    % use connections from nullspace
    for jCon = 1:numel(loopInfo.conCompAlwaysOn)
        if ~conCompOn(jCon)
            % relax constraints not affecting optimality and feasibility
            MILPproblemLLC.b(loopInfo.con.vU(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = bigM;
            MILPproblemLLC.b(loopInfo.con.gU(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = bigM;
            MILPproblemLLC.b(loopInfo.con.vL(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = -bigM;
            MILPproblemLLC.b(loopInfo.con.gL(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = -bigM;
            % fix variables not affecting optimality and feasibility
            MILPproblemLLC.lb(loopInfo.var.g(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = 0;
            MILPproblemLLC.ub(loopInfo.var.g(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = 0;
            MILPproblemLLC.ub(loopInfo.var.z(loopInfo.rxnInLoopIds(loopInfo.conComp == jCon))) = 0;
        end
    end
else
    % use connections from EFMs
    rxnOn = loopInfo.rxnInLoopsAlwaysOn;
    rxnOn(rxnID) = true;

    % reactions in cycles not sharing EFMs with the current rxns and
    % not being one of the reactions required to have no flux through cycles
    id = ~any(loopInfo.rxnLink(rxnOn, :), 1)' & any(loopInfo.rxnInLoops, 2);
    % the loop constraints on them are relaxed
    MILPproblemLLC.b(loopInfo.con.vU(loopInfo.rxnInLoopIds(id))) = bigM;
    MILPproblemLLC.b(loopInfo.con.gU(loopInfo.rxnInLoopIds(id))) = bigM;
    MILPproblemLLC.b(loopInfo.con.vL(loopInfo.rxnInLoopIds(id))) = -bigM;
    MILPproblemLLC.b(loopInfo.con.gL(loopInfo.rxnInLoopIds(id))) = -bigM;

    % pre-determine variables not connected to the reaction for FVA
    % except reactions required to be always constrained
    rxnKeep = loopInfo.conComp == 0;
    for jCon = 1:numel(conCompOn)
        if conCompOn(jCon)
            rxnKeep(loopInfo.conComp == jCon) = true;
        end
    end
    MILPproblemLLC.lb(loopInfo.var.g(loopInfo.rxnInLoopIds(~rxnKeep))) = 0;
    MILPproblemLLC.ub(loopInfo.var.g(loopInfo.rxnInLoopIds(~rxnKeep))) = 0;
    MILPproblemLLC.ub(loopInfo.var.z(loopInfo.rxnInLoopIds(~rxnKeep))) = 0;
end
end

function MILPproblemLLC = restoreOriginalBounds(MILPproblemLLC, rhs0, varInd, BDg)
    MILPproblemLLC.b = rhs0;
    MILPproblemLLC.ub(varInd.z) = 1;
    MILPproblemLLC.ub(varInd.g) = BDg;
    MILPproblemLLC.lb(varInd.g) = -BDg;
end

author: opencobra

repo: cobratoolbox


text:
function hRD = plugin_rd(board, design)
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Construct reference design object
hRD = hdlcoder.ReferenceDesign('SynthesisTool', 'Xilinx Vivado');

% Create the reference design for the SOM-only
% This is the base reference design that other RDs can build upon
hRD.ReferenceDesignName = sprintf('%s (%s)', upper(board), design);

% Determine the board name based on the design
hRD.BoardName = sprintf('AnalogDevices ADRV9361-Z7035 AGC');

% Tool information
hRD.SupportedToolVersion = {'2018.2'};

% Add AGC control
hRD.addParameter( ...
    'ParameterID',   'en_agc', ...
    'DisplayName',   'Control AGC', ...
    'DefaultValue',  'Rx', ...
    'ParameterType',  hdlcoder.ParameterType.Dropdown, ...
    'Choice',       {'Rx', 'Tx'});

% Get the root directory
rootDir = fileparts(strtok(mfilename('fullpath'), '+'));
tmp = strsplit(rootDir,'/');
rootDir = fullfile('/',tmp{1:end-2});
rootDirBSP = fullfile('hdl_wa_bsp','vendor','AnalogDevices','vivado');

% Design files are shared
hRD.SharedRD = true;
hRD.SharedRDFolder = rootDir;

board = 'ccfmc_lvds_agc';

hRD.addCustomVivadoDesign( ...
    'CustomBlockDesignTcl', fullfile('targeting_models', 'tuneAGC-ad9361',  'ccfmc_lvds_agc', 'system_project_rx.tcl'), ...
    'CustomTopLevelHDL',    fullfile('targeting_models', 'tuneAGC-ad9361',  'ccfmc_lvds_agc', 'system_top.v'));

hRD.BlockDesignName = 'system';

% custom constraint files
board_type = strsplit(board,'_');
hRD.CustomConstraints = {...
    fullfile(rootDirBSP, 'projects', 'adrv9361z7035', 'common', strcat(board_type{1}, '_constr.xdc')), ...
	fullfile(rootDirBSP, 'projects', 'adrv9361z7035', 'common', 'adrv9361z7035_constr.xdc'), ...
	fullfile(rootDirBSP, 'projects', 'adrv9361z7035', 'common', strcat('adrv9361z7035_constr_', board_type{2}, '.xdc')), ...
    };

% custom source files
hRD.CustomFiles = {...
	fullfile(rootDirBSP, 'library')...,
	fullfile(rootDirBSP, 'library','xilinx')...,
	fullfile(rootDirBSP, 'projects','common')...,
	fullfile(rootDirBSP, 'projects','scripts')...,
	fullfile(rootDirBSP, 'projects','fmcomms2')...,
	fullfile(rootDirBSP, 'projects','adrv9361z7035', 'common')...,
    fullfile('targeting_models', 'tuneAGC-ad9361', 'ccfmc_lvds_agc')...,
    };

%% Add interfaces
% add clock interface
hRD.addClockInterface( ...
    'ClockConnection',   'util_ad9361_divclk/clk_out', ...
    'ResetConnection',   'util_ad9361_divclk_reset/peripheral_aresetn');

author: analogdevicesinc


text:
%{
load('dataset/trafficdb/traffic_patches.mat');
[M,m,n,p] = convert_video3d_to_2d(im2double(imgdb{100}));
out = run_algorithm('MC', 'LMaFit', M, [])
show_results(M.*out.Omega,out.L,out.S,out.O,p,m,n);
%}

MIdx = M(Idx);
rank = 10;
[X,Y] = lmafit_mc_adp(size(M,1),size(M,2),rank,Idx,MIdx,[]);
L = X*Y; % low-rank
S = (M - L); % sparse

author: andrewssobral

repo: lrslibrary


text:
function strings = GetParameterStrings(modelname)
%
% function strings = GetParameterStrings(modelname)
%
% Given an input modelname, this function returns the names of the model
% parameters
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

if (strcmp(modelname, 'StickIsoV_B0'))
    strings = {'ficvf', 'di', 'dh', 'fiso', 'diso', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'StickTortIsoV_B0'))
    strings = {'ficvf', 'di', 'fiso', 'diso', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'WatsonSHStick'))
    strings = {'ficvf', 'di', 'dh', 'kappa', 'theta', 'phi'};
elseif (strcmp(modelname, 'WatsonSHStickIsoV_B0'))
    strings = {'ficvf', 'di', 'dh', 'kappa', 'fiso', 'diso', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'WatsonSHStickIsoVIsoDot_B0'))
    strings = {'ficvf', 'di', 'dh', 'kappa', 'fiso', 'diso', 'irfrac', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'WatsonSHStickTortIsoV_B0'))
    strings = {'ficvf', 'di', 'kappa', 'fiso', 'diso', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'WatsonSHStickTortIsoVIsoDot_B0'))
    strings = {'ficvf', 'di', 'kappa', 'fiso', 'diso', 'irfrac', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'BinghamStickTortIsoV_B0'))
    strings = {'ficvf', 'di', 'kappa', 'beta', 'psi', 'fiso', 'diso', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'WatsonSHCylSingleRadTortIsoV_GPD_B0'))
    strings = {'ficvf', 'di', 'rad', 'kappa', 'fiso', 'diso', 'b0', 'theta', 'phi'};
elseif (strcmp(modelname, 'CylSingleRadIsoDotTortIsoV_GPD_B0'))
    strings = {'ficvf', 'di', 'rad', 'irfrac', 'fiso', 'diso', 'b0', 'theta', 'phi'};
else
    error(['Parameter strings yet to be defined for this model:', modelname]);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function D = EuDist2(fea_a,fea_b,bSqrt)
%EUDIST2 Efficiently Compute the Euclidean Distance Matrix by Exploring the
%Matlab matrix operations.
%
%   D = EuDist(fea_a,fea_b)
%   fea_a:    nSample_a * nFeature
%   fea_b:    nSample_b * nFeature
%   D:      nSample_a * nSample_a
%       or  nSample_a * nSample_b
%
%    Examples:
%
%       a = rand(500,10);
%       b = rand(1000,10);
%
%       A = EuDist2(a); % A: 500*500
%       D = EuDist2(a,b); % D: 500*1000
%
%   version 2.1 --November/2011
%   version 2.0 --May/2009
%   version 1.0 --November/2005
%
%   Written by Deng Cai (dengcai AT gmail.com)

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

if (~exist('fea_b','var')) || isempty(fea_b)
    aa = sum(fea_a.*fea_a,2);
    ab = fea_a*fea_a';

    if issparse(aa)
        aa = full(aa);
    end

    D = bsxfun(@plus,aa,aa') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
    D = max(D,D');
else
    aa = sum(fea_a.*fea_a,2);
    bb = sum(fea_b.*fea_b,2);
    ab = fea_a*fea_b';

    if issparse(aa)
        aa = full(aa);
        bb = full(bb);
    end

    D = bsxfun(@plus,aa,bb') - 2*ab;
    D(D<0) = 0;
    if bSqrt
        D = sqrt(D);
    end
end

author: zzstefan

repo: BrainNetClass


text:
%% Example
% An example script, it is used to show how to use the different
% functions of the KUKA iiwa matlab toolbox
% First start the server on the KUKA iiwa controller
% Then run the following script in Matlab

% Mohammad SAFEEA, 3rd of May 2017

close all,clear all,clc

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else

    %% move to initial position
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % initial confuguration
relVel=0.15; % relative velocity
movePTPJointSpace( t_Kuka , pinit, relVel); % point to point motion in joint space

      %% PTP motion
      jPos={}
      homePos={}
     [ jPos ] = getJointsPos( t_Kuka ) % get current joints position

      for ttt=1:7  % home position
          homePos{ttt}=0;
      end

      pause(0.1)
      relVel=0.15;
      movePTPJointSpace( t_Kuka , homePos, relVel); % go to home position
      pause(0.1)
      movePTPJointSpace( t_Kuka , jPos, relVel); % return back to original position

      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function phi = cf_merton(u,lnS,T,r,d,sigma,a,b,lambda)
% characteristic function for the Merton model
    phi = exp(cf_black(u,lnS,T,r,d,sigma) ...
        + cf_lognormjump(u,a,b,lambda,T));
end


text:
function Cont=Contingency(Mem1,Mem2)
%CONTINGENCY Form contigency matrix for two vectors
% C=Contingency(Mem1,Mem2) returns contingency matrix for two
% column vectors Mem1, Mem2. These define which cluster each entity
% has been assigned to.
%
% See also RANDINDEX.
%

%(C) David Corney (2000)   		D.Corney@cs.ucl.ac.uk

if nargin < 2 | min(size(Mem1)) > 1 | min(size(Mem2)) > 1
   error('Contingency: Requires two vector arguments')
   return
end

Cont=zeros(max(Mem1),max(Mem2));

for i = 1:length(Mem1);
   Cont(Mem1(i),Mem2(i))=Cont(Mem1(i),Mem2(i))+1;
end


text:
function x = poisson_solution ( nrow, ncol, rhs_num )

%*****************************************************************************80
%
%% POISSON_SOLUTION returns the solution of a Poisson linear system.
%
%  Discussion:
%
%    The Poisson matrix is associated with an NROW by NCOL rectangular
%    grid of points.
%
%    Assume that the points are numbered from left to right, bottom to top.
%
%    If the K-th point is in row I and column J, set X = I + J.
%
%    This will be the solution to the linear system.
%
%  Example:
%
%    NROW = 3, NCOL = 3
%
%    ^
%    |  7  8  9
%    J  4  5  6
%    |  1  2  3
%    |
%    +-----I---->
%
%    Solution vector X = ( 2, 3, 4, 3, 4, 5, 4, 5, 6 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, Charles Van Loan,
%    Matrix Computations, second edition,
%    Johns Hopkins University Press, Baltimore, Maryland, 1989
%    (Section 4.5.4).
%
%  Parameters:
%
%    Input, integer NROW, NCOL, the number of rows and columns
%    in the grid.
%
%    Input, integer RHS_NUM, the number of right hand sides.
%
%    Output, real X(NROW*NCOL,RHS_NUM), the solution.
%
  n = nrow * ncol;

  x = zeros ( n, rhs_num );

  k = 0;
  for j = 1 : nrow
    for i = 1 : ncol
      k = k + 1;
      x(k,1) = i + j;
    end
  end

  return
end

author: johannesgerer


text:
function [whitestats PCstats stabilitystats residualstats] = est_validateMVAR(varargin)
%
% Validate a fitted VAR model. With two inputs, this function generates a
% GUI where the validation scheme can be specified. Validation consists of
% statistical tests for "whiteness" of fitted VAR model residuals [1,2],
% consistency of the fitted model [1,3] and stability of fitted model [1-3]
%
% Input:
%
%   ALLEEG:     Array of EEGLAB data structures containing fitted MODEL
%   typeproc:   reserved for future use. Use 0
%
% Optional:
%
%     'whitenessCriteria':    Cell array containing names of residual whiteness
%                             test criteria to evaluate. See [1, 2] for details.
%                             Possible Values: {'Ljung-Box','ACF','Box-Pierce','Li-McLeod'}
%                             Default Value  : all
%                             Data Input Type: cell array
%
%
%     'checkWhiteness':       Whether or not to check whiteness.
%                             Default Value  : true
%                             Data Input Type: boolean
%
%     'checkConsistency':     Whether or not to check consistency. See [1,3]
%                             for details.
%                             Default Value  : true
%                             Data Input Type: boolean
%
%     'checkStability':       Whether or not to check stability. See
%                             [1-3] for details.
%                             Default Value  : true
%                             Data Input Type: boolean
%
%     'alpha':                significance level for determining whiteness
%                             Data Input Range: [0 1]
%                             Default Value   : 0.05
%                             Data Input Type : real number (double)
%
%     'prctWinToSample':      percent of time windows to randomly select
%                             Data Input Range: [0 100]
%                             Default Value   : 100
%                             Data Input Type : real number (double)
%
%     'verb':                 verbosity level (0=no output, 1=text, 2=gui)
%
%
% Output:
%
%     whitestats:             Structure containing whiteness statistics.
%                             See est_checkMVARWhiteness() for details on
%                             structure format
%
%     PC:                     Vector of percent consistency estimates for
%                             each window.
%
%     stability:              Vector of stability estimates for each window
%
% See Also: est_checkMVARWhiteness(), est_checkMVARStability(),
%           est_checkMVARConsistency, pop_est_fitMVAR()
%
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 3.6 and 6.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% [2] Lutkepohl, H. (2007) New Introduction to Time Series Analysis.
%   Springer.
%
% [3] Ding M, Bressler SL, Yang W, Liang H (2000) Short-window spectral
%   analysis of cortical event-related potentials by adaptive multivariate
%   autoregressive modeling: data preprocessing, model validation, and
%   variability assessment. Biol. Cybern. 83:35-45
%
% Author: Tim Mullen, 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

verb = arg_extract(varargin,{'verb','VerbosityLevel'},[],0);
prctWinToSample = arg_extract(varargin,{'prctWinToSample','WindowSamplePercent'},[],100);

g = arg_define([0 1],varargin, ...
    arg_norep({'EEG','ALLEEG'},mandatory,[],'EEGLAB dataset'), ...
    arg_subtoggle({'checkWhiteness','CheckResidualWhiteness'},{'verb',verb,'prctWinToSample',prctWinToSample},@est_checkMVARWhiteness,sprintf('Check residual whiteness. \nResiduals are "white" if they are uncorrelated.'),'cat','Validation Methods','suppress',{'prctWinToSample','VerbosityLevel'}), ...
    arg_subtoggle({'checkResidualVariance','CheckResidualVariance'},{'verb',verb,'prctWinToSample',prctWinToSample,'whitenessCriteria',{}},@est_checkMVARWhiteness,'Compute residual autocorrelation and variance','cat','Validation Methods','suppress',{'WhitenessCriteria','MultipleComparisonsCorrection','SignificanceLevel','prctWinToSample','VerbosityLevel'}), ...
    arg_subtoggle({'checkConsistency','CheckConsistency'},{'verb',verb,'prctWinToSample',prctWinToSample},@est_checkMVARConsistency,sprintf('Check model consistency. \nA "consistent" model is able to generate data with similar covariance structure as the original data'),'cat','Validation Methods','suppress',{'prctWinToSample','VerbosityLevel'}), ...
    arg_subtoggle({'checkStability','CheckStability'},{'verb',verb},@est_checkMVARStability,sprintf('Check model stability. \nA stable model is one that cannot "blow up" (e.g. generate data that increases to infinity) -- this is a requirement for causal modeling.'),'cat','Validation Methods','suppress',{'prctWinToSample','VerbosityLevel'}), ...
    arg({'prctWinToSample','WindowSamplePercent'},100,[1 100],'Percent of windows to sample','cat','Data Reduction'), ...
    arg_nogui({'winStartIdx','WindowStartIndices'},[],[],'Starting indices for windows. This is a vector of sample points (start of windows) at which to estimate windowed VAR model. Default is empty (use all windows)','cat','Data Reduction'), ...
    arg({'verb','VerbosityLevel'},2,{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical'), ...
    arg({'plot','PlotResults'},true,[],'Plot results') ...
    );

% commit EEG variable to workspace
% [data g] = hlp_splitstruct(g,{'g.EEG'});
% arg_toworkspace(data);
% clear data;

% initialize default output
[whitestats PCstats stabilitystats residualstats] = deal([]);

if ~isfield(g.EEG.CAT,'MODEL')
    error('est_validateMVAR:NoModel','g.EEG.CAT.MODEL must be present in the dataset');
else
    MODEL = g.EEG.CAT.MODEL;
end

% determine which windows to use
if isempty(g.winStartIdx)
    % starting point of each window (points)
    g.winStartIdx  = round(MODEL.winStartTimes*g.EEG.srate)+1;
end

if g.prctWinToSample<100
    % randomly select percentage of windows to work with
    randwin = randperm(length(g.winStartIdx));
    randwin = sort(randwin(1:ceil(length(g.winStartIdx)*g.prctWinToSample/100)));
    g.winStartIdx = g.winStartIdx(randwin);
    g.prctWinToSample = 100;
end

if g.verb==2
    % create waitbar
    waitbarTitle = sprintf('Validating Model %s...', ...
        fastif(isempty(g.EEG.condition),'',['for ' g.EEG.condition]));

    multiWaitbar(waitbarTitle,'Reset');
    multiWaitbar(waitbarTitle,'ResetCancel',true);
    multiWaitbar(waitbarTitle, ...
        'Color', hlp_getNextUniqueColor('reset'), ...
        'CanCancel','on', ...
        'CancelFcn',@(a,b) disp('[Cancel requested. Please wait...]'));
    %     wb_cleanup = onCleanup(@() multiWaitbar(waitbarTitle,'Close'));
end

% determine number of checks we'll do (for progbar)
numchecks = sum([g.checkConsistency.arg_selection  ...
    g.checkWhiteness.arg_selection    ...
    g.checkStability.arg_selection    ...
    g.checkResidualVariance.arg_selection]);
curcheck  = 0;

% check residual whiteness
% -------------------------------------------------------------------------
if g.checkWhiteness.arg_selection
    [whitestats residualstats] = est_checkMVARWhiteness(...
                                    'EEG',g.EEG,                        ...
                                    'MODEL',MODEL,                      ...
                                     g.checkWhiteness,                  ...
                                    'winStartIdx',g.winStartIdx,        ...
                                    'prctWinToSample',g.prctWinToSample,...
                                    'verb',g.verb);
    if isempty(whitestats)
        % whiteness checks cancelled
        g.checkWhiteness.arg_selection = false;
        numchecks = numchecks - 1;
        curcheck  = curcheck  - 1;
    else
        winTimes = whitestats.winStartTimes;
    end
    if g.verb==2 && ~updateWaitbar(), return; end
end

% check residual variance
% -------------------------------------------------------------------------
if g.checkResidualVariance.arg_selection
    if ~g.checkWhiteness.arg_selection
        [~, residualstats] = est_checkMVARWhiteness(...
                                'EEG',g.EEG,                            ...
                                'MODEL',MODEL,                          ...
                                 g.checkResidualVariance,               ...
                                'whitenessCriteria',{},                 ...
                                'winStartIdx',g.winStartIdx,            ...
                                'prctWinToSample',g.prctWinToSample,    ...
                                'verb',g.verb);
    end
    % bookeeping
    if isempty(residualstats)
        % whiteness checks cancelled
        g.checkResidualVariance.arg_selection = false;
        numchecks = numchecks - 1;
        curcheck  = curcheck  - 1;
    else
        winTimes = residualstats.winStartTimes;
    end
    if g.verb==2 && ~updateWaitbar(), return; end
end

% check model consistency
% -------------------------------------------------------------------------
if g.checkConsistency.arg_selection
    PCstats = est_checkMVARConsistency(...
                    'EEG',g.EEG,                        ...
                    'MODEL',MODEL,                      ...
                     g.checkConsistency,                ...
                    'winStartIdx',g.winStartIdx,        ...
                    'prctWinToSample',g.prctWinToSample,...
                    'verb',g.verb);
    % bookeeping
    if isempty(PCstats)
        % consistency checks cancelled
        g.checkConsistency.arg_selection = false;
        numchecks = numchecks - 1;
        curcheck  = curcheck  - 1;
    else
        winTimes = PCstats.winStartTimes;
    end
    if g.verb==2 && ~updateWaitbar(), return; end
end

% check model stability
% -------------------------------------------------------------------------
if g.checkStability.arg_selection
    stabilitystats = est_checkMVARStability(...
                        'EEG',g.EEG,                        ...
                        'MODEL',MODEL,                      ...
                         g.checkStability,                  ...
                        'winStartIdx',g.winStartIdx,        ...
                        'prctWinToSample',g.prctWinToSample,...
                        'verb',g.verb);
    % bookeeping
    if isempty(stabilitystats)
        % consistency checks cancelled
        g.checkStability.arg_selection = false;
        numchecks = numchecks - 1;
        curcheck  = curcheck  - 1;
    else
        winTimes = stabilitystats.winStartTimes;
    end
    if g.verb==2 && ~updateWaitbar(), return; end
end

% plot results
% -------------------------------------------------------------------------
if g.plot && numchecks > 0
    if g.checkWhiteness.arg_selection
        whitenessCriteria = g.checkWhiteness.whitenessCriteria;
    else
        whitenessCriteria = {};
    end
    vis_plotModelValidation({whitestats},{PCstats},{stabilitystats}, ...
        'whitenessCriteria',whitenessCriteria,                       ...
        'checkWhiteness',g.checkWhiteness.arg_selection,             ...
        'checkConsistency',g.checkConsistency.arg_selection,         ...
        'checkStability',g.checkStability.arg_selection,             ...
        'conditions',{g.EEG.condition},                              ...
        'windowTimes',winTimes);
end

% % clean up
if g.verb==2
    multiWaitbar(waitbarTitle,'Close');
end

    % subfunction updateWaitbar
    % ---------------------------------------------------------------------
    function nocancel = updateWaitbar()
        nocancel = true;
        curcheck = curcheck + 1;
        drawnow;
        cancel = multiWaitbar(waitbarTitle,curcheck/numchecks);
        if cancel && hlp_confirmWaitbarCancel(waitbarTitle, ...
                     'Are you sure you want to cancel all validation?')
            nocancel = false;
            return;
        end
    end

end

repo: REST