text:
function varargout = Bedienpult(varargin)
%BEDIENPULT M-file for Bedienpult.fig
%      BEDIENPULT, by itself, creates a new BEDIENPULT or raises the existing
%      singleton*.
%
%      H = BEDIENPULT returns the handle to a new BEDIENPULT or the handle to
%      the existing singleton*.
%
%      BEDIENPULT('Property','Value',...) creates a new BEDIENPULT using the
%      given property value pairs. Unrecognized properties are passed via
%      varargin to Bedienpult_OpeningFcn.  This calling syntax produces a
%      warning when there is an existing singleton*.
%
%      BEDIENPULT('CALLBACK') and BEDIENPULT('CALLBACK',hObject,...) call the
%      local function named CALLBACK in BEDIENPULT.M with the given input
%      arguments.
%
%      *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 Bedienpult

% Last Modified by GUIDE v2.5 18-Apr-2011 17:03:05

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Bedienpult_OpeningFcn, ...
                   'gui_OutputFcn',  @Bedienpult_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 Bedienpult is made visible.
function Bedienpult_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   unrecognized PropertyName/PropertyValue pairs from the
%            command line (see VARARGIN)

  % ------------------------------------------------------------
  % If Cruise.mdl is not open: open it and bring GUI to front
  % Also set model parameters to match GUI settings
  % ------------------------------------------------------------

  model_open(handles)

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

% Set the sliders value to correspond the initial value
set(handles.AltitudeSlider,'Value',...
    str2double(get(handles.AltitudeCurrentValue,'String')))
set(handles.BearingSlider,'Value',...
    str2double(get(handles.BearingCurrentValue,'String')))
set(handles.CameraSlider,'Value',...
    str2double(get(handles.CameraCurrentValue,'String')))
set(handles.VelocitySlider,'Value',...
    str2double(get(handles.VelocityCurrentValue,'String')))

% Update handles structure
guidata(hObject, handles);

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

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

%endfunction Bedienpult_OutputFcn

% ------------------------------------------------------------
% Ensure that the Simulink model is open
% ------------------------------------------------------------

function model_open(handles)
% Make sure the diagram is still open
if isempty(find_system('Name','Cruise')),
  % Open the model file 'EEG_BERT_FACE.mdl'
  open_system('Cruise');

  set_param('Cruise/Altitude','Value',...
      get(handles.AltitudeCurrentValue,'String'))
  set_param('Cruise/Azimuth', 'Value',...
      get(handles.CameraCurrentValue,'String'))
   set_param('Cruise/Bearing', 'Value',...
      get(handles.BearingCurrentValue,'String'))
   set_param('Cruise/Velocity', 'Value',...
      get(handles.VelocityCurrentValue,'String'))

end
%endfunction model_open

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

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

% Ensure model is open
  model_open(handles)

% Get the new value for the Altitude value
  NewStrVal = get(hObject,'String');
  NewVal = str2double(NewStrVal);

% Check that the entered value falls within the allowable range
  if  isempty(NewVal) || (NewVal< 0),
    % Revert to last value, as indicated by AltitudeSlider
    OldVal = get(handles.AltitudeSlider,'Value');
    set(hObject,'String',OldVal)
  else
    % Set the value of the AltitudeSlider to the new value
    set(handles.AltitudeSlider,'Value',NewVal)

    % Set the Value parameter of the Altitude Constant Block to the
    % new value
    set_param('Cruise/Altitude','Value',NewStrVal)
  end

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function AltitudeCurrentValue_CreateFcn(hObject, eventdata, handles)
% hObject    handle to AltitudeCurrentValue (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 slider movement.
function AltitudeSlider_Callback(hObject, eventdata, handles)
% hObject    handle to AltitudeSlider (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
% Ensure model is open
model_open(handles)

% Get the new value for the Altitude value from the slider
NewVal = get(hObject,'Value');

% Set the value of the AltitudeCurrentValue to the new value set by the
% slider
set(handles.AltitudeCurrentValue,'String',NewVal)

% Set the value of the Altitude Constant Block to the new value
set_param('Cruise/Altitude','Value',num2str(NewVal))

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

% Start Simulation
set_param('Cruise','SimulationCommand','start');

% --- Executes on button press in StopButton.
function StopButton_Callback(hObject, eventdata, handles)
% hObject    handle to StopButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set_param('Cruise','SimulationCommand','stop');

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

% Hints: get(hObject,'String') returns contents of CameraCurrentValue as text
%        str2double(get(hObject,'String')) returns contents of CameraCurrentValue as a double
% Ensure model is open
  model_open(handles)

% Get the new value for the Altitude value
  NewStrVal = get(hObject,'String');
  NewVal = str2double(NewStrVal);

% Check that the entered value falls within the allowable range
  if  isempty(NewVal) || (NewVal< 0),
    % Revert to last value, as indicated by AltitudeSlider
    OldVal = get(handles.CameraSlider,'Value');
    set(hObject,'String',OldVal)
  else
    % Set the value of the AltitudeSlider to the new value
    set(handles.CameraSlider,'Value',NewVal)

    % Set the Value parameter of the Altitude Constant Block to the
    % new value
    set_param('Cruise/Azimuth','Value',NewStrVal)
  end

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function CameraCurrentValue_CreateFcn(hObject, eventdata, handles)
% hObject    handle to CameraCurrentValue (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 slider movement.
function CameraSlider_Callback(hObject, eventdata, handles)
% hObject    handle to CameraSlider (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

% Ensure model is open
model_open(handles)

% Get the new value for the Altitude value from the slider
NewVal = get(hObject,'Value');

% Set the value of the AltitudeCurrentValue to the new value set by the
% slider
set(handles.CameraCurrentValue,'String',NewVal)

% Set the value of the Altitude Constant Block to the new value
set_param('Cruise/Azimuth','Value',num2str(NewVal))

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

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

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

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

% Ensure model is open
  model_open(handles)

% Get the new value for the Altitude value
  NewStrVal = get(hObject,'String');
  NewVal = str2double(NewStrVal);

% Check that the entered value falls within the allowable range
  if  isempty(NewVal) || (NewVal< 0) || (NewVal>360),
    % Revert to last value, as indicated by AltitudeSlider
    OldVal = get(handles.BearingSlider,'Value');
    set(hObject,'String',OldVal)
  else
    % Set the value of the AltitudeSlider to the new value
    set(handles.BearingSlider,'Value',NewVal)

    % Set the Value parameter of the Altitude Constant Block to the
    % new value
    set_param('Cruise/Bearing','Value',NewStrVal)

    % Set the value for Camera to the same value as Bearing
    set_param('Cruise/Azimuth','Value',NewStrVal)
    set(handles.CameraCurrentValue,'String',NewStrVal)
  end

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function BearingCurrentValue_CreateFcn(hObject, eventdata, handles)
% hObject    handle to BearingCurrentValue (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 slider movement.
function BearingSlider_Callback(hObject, eventdata, handles)
% hObject    handle to BearingSlider (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
model_open(handles)

% Get the new value for the Altitude value from the slider
NewVal = get(hObject,'Value');

% Set the value of the AltitudeCurrentValue to the new value set by the
% slider
set(handles.BearingCurrentValue,'String',NewVal)

% Set the value of the Altitude Constant Block to the new value
set_param('Cruise/Bearing','Value',num2str(NewVal))
set_param('Cruise/Azimuth','Value',num2str(NewVal))
set(handles.CameraCurrentValue,'String',num2str(NewVal))
set(handles.CameraSlider,'Value',NewVal)

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

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

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

% Hints: get(hObject,'String') returns contents of VelocityCurrentValue as text
%        str2double(get(hObject,'String')) returns contents of VelocityCurrentValue as a double
% Ensure model is open
  model_open(handles)

% Get the new value for the Altitude value
  NewStrVal = get(hObject,'String');
  NewVal = str2double(NewStrVal);

% Check that the entered value falls within the allowable range
  if  isempty(NewVal) || (NewVal< 0) || (NewVal>360),
    % Revert to last value, as indicated by AltitudeSlider
    OldVal = get(handles.VelocitySlider,'Value');
    set(hObject,'String',OldVal)
  else
    % Set the value of the AltitudeSlider to the new value
    set(handles.VelocitySlider,'Value',NewVal)

    % Set the Value parameter of the Altitude Constant Block to the
    % new value
    set_param('Cruise/Velocity','Value',NewStrVal)

    % Set the value for Camera to the same value as Bearing
    set_param('Cruise/Velocity','Value',NewStrVal)
    set(handles.VelocityCurrentValue,'String',NewStrVal)
  end

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function VelocityCurrentValue_CreateFcn(hObject, eventdata, handles)
% hObject    handle to VelocityCurrentValue (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 slider movement.
function VelocitySlider_Callback(hObject, eventdata, handles)
% hObject    handle to VelocitySlider (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
model_open(handles)

% Get the new value for the Altitude value from the slider
NewVal = get(hObject,'Value');

% Set the value of the AltitudeCurrentValue to the new value set by the
% slider
set(handles.VelocityCurrentValue,'String',NewVal)

% Set the value of the Altitude Constant Block to the new value
set_param('Cruise/Velocity','Value',num2str(NewVal))
set(handles.VelocitySlider,'Value',NewVal)

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

% Get the value from Bearing
    NewStrVal = get(handles.BearingCurrentValue,'String');

 % Set the value for Camera to the same value as Bearing
    set_param('Cruise/Azimuth','Value',NewStrVal)
    set(handles.CameraCurrentValue,'String',NewStrVal)

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

% Hints: get(hObject,'String') returns contents of TiltCurrentValue as text
%        str2double(get(hObject,'String')) returns contents of TiltCurrentValue as a double
  model_open(handles)

% Get the new value for the Altitude value
  NewStrVal = get(hObject,'String');
  NewVal = str2double(NewStrVal);

% Check that the entered value falls within the allowable range
  if  isempty(NewVal) || (NewVal< 0) || (NewVal>360),
    % Revert to last value, as indicated by AltitudeSlider
    OldVal = get(handles.TiltSlider,'Value');
    set(hObject,'String',OldVal)
  else
    % Set the value of the AltitudeSlider to the new value
    set(handles.TiltSlider,'Value',NewVal)

    % Set the Value parameter of the Altitude Constant Block to the
    % new value
    set_param('Cruise/Tilt','Value',NewStrVal)

    set(handles.TiltCurrentValue,'String',NewStrVal)
  end

guidata(hObject, handles);

% --- Executes during object creation, after setting all properties.
function TiltCurrentValue_CreateFcn(hObject, eventdata, handles)
% hObject    handle to TiltCurrentValue (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 slider movement.
function TiltSlider_Callback(hObject, eventdata, handles)
% hObject    handle to TiltSlider (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
model_open(handles)

% Get the new value for the Altitude value from the slider
NewVal = get(hObject,'Value');

% Set the value of the AltitudeCurrentValue to the new value set by the
% slider
set(handles.TiltCurrentValue,'String',NewVal)

% Set the value of the Altitude Constant Block to the new value
set_param('Cruise/Tilt','Value',num2str(NewVal))
set(handles.TiltSlider,'Value',NewVal)

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


text:
function [Psi2 P] = biasVardistPsi2Compute(biaskern, vardist, Z)

% BIASVARDISTPSI2COMPUTE one line description
% FORMAT
% DESC description
% RETURN Psi2 : description
% RETURN P : description
% ARG biasKern : the kernel structure associated with the bias kernel.
% ARG vardist : description
% ARG Z : description
%
%
% SEEALSO : others
%
%
% COPYRIGHT : Michalis K. Titsias, 2009
%

% VARGPLVM

Psi2 = repmat(vardist.numData*(biaskern.variance^2),size(Z,1),size(Z,1));
P = [];

author: SheffieldML

repo: vargplvm


text:
function H = makeBaseLdpccc(N, numOfOne)
% Create LDPC-CC base parity check matrix
%
%  N        : Matrix dimension (N X N)
%  numOfOne : Number of ones per column or row
%
%  H        : Base LDPC-CC matrix
%
% Bagawan S. Nugroho 2007

% Number of ones per column or row
onePerCol = numOfOne;
onePerRow = numOfOne;

% Evenboth: For each column and row, place 1s uniformly at random

% Distribute 1s uniformly at random within column
for i = 1:N
   onesInCol(:, i) = randperm(N)';
end

% Create non zero elements (1s) index
r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
tmp = repmat([1:N], onePerCol, 1);
c = reshape(tmp, N*onePerCol, 1);

% Make the number of 1s between rows as uniform as possible

% Order row index
[r, ix] = sort(r);

% Order column index based on row index
for i = 1:N*onePerCol
   cSort(i, :) = c(ix(i));
end

% Create new row index with uniform weight
tmp = repmat([1:N], onePerRow, 1);
r = reshape(tmp, N*onePerCol, 1);

% Create sparse matrix H
% Remove any duplicate non zero elements index using logical AND
S = and(sparse(r, cSort, 1, N, N), ones(N, N));
H = full(S);

% Check rows that have no 1 or only have one 1
for i = 1:N

   n = randperm(N);
   % Add two 1s if row has no 1
   if length(find(r == i)) == 0
      H(i, n(1)) = 1;
      H(i, n(2)) = 1;
   % Add one 1 if row has only one 1
   elseif length(find(r == i)) == 1
      H(i, n(1)) = 1;
   end

end % for i

% Add 0/1 bits (see text)
H = [ones(N, 1) H zeros(N, 1)];

fprintf('LDPC-CC base matrix is created.\n');

author: Sable

repo: mcbench-benchmarks


text:
function vscale = getvscl( f )
%GETVSCL(F)   Estimate the vertical scale of a function.
%   VSCALE = GETVSCL(F) estimates the vertical scale (also known as the
%   dynamical range) of a function. This is required because a CHEBTECH does not
%   store its interpolation data. If F is an array-valued CHEBTECH with K
%   columns, then VSCALE is a row vector of length K.

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

% Check isempty:
if ( isempty( f.coeffs ) )
    vscale = 0;
else
    % compute values
    values = f.coeffs2vals( f.coeffs );
    vscale = max(abs(values), [], 1);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = yalmiptest(prefered_solver,auto)
%YALMIPTEST Runs a number of test problems.
%
%   YALMIPTEST is recommended when a new solver or a new version
%   of YALMIP installed.
%
%   EXAMPLES
%    YALMIPTEST               % Without argument, default solver used
%    YALMIPTEST('solver tag') % Test with specified solver
%    YALMIPTEST(options)      % Test with specific options structure from
%
%   See also SDPSETTINGS

if ~exist('sedumi2pen.m')
    disp('Add /yalmip/extras etc to your path first...')
    disp('Read the Installation notes.')
    return
end

if ~exist('callsedumi.m')
    disp('Still missing paths...Just do an addpath(genpath(''yalmiprootdirectory''));')
    disp('Read the Installation notes.')
    return
end

% SDPT3 has a function called constraint.m which causes issues
detected = which('constraint.m');
if isa(detected,'cell')
    if length(detected)>0
        if isempty(strfind(detected{1},'extras\@constraint'))
            clc
            disp('You seem to have some other toolbox with a function called constraint.m');
            disp('Delete that toolbox, or delete the function/class, or change path so that YALMIP is on top.');
            disp(detected{1})
        end
        return
    end
end

detected = which('yalmip.m','-all');
% Will not work in Octave as Octave only reports first item found?
if isa(detected,'cell')
    if length(detected)>1
        clc
        disp('You seem to have multiple installations of YALMIP in your path.')
        disp('Please correct this...');
        disp(detected)
        return
    end
end

% Pagination really doesn't work well with solvers
more off

if exist('OCTAVE_VERSION', 'builtin')
    OctaveRunning = 1;
else
    OctaveRunning = 0;
end

donttest = 0;
if (nargin==1) && isa(prefered_solver,'char') && strcmp(prefered_solver,'test')
    donttest = 0;
    prefered_solver = '';
else
    donttest = 1;
end

if nargin==0
    prefered_solver = '';
else
    if ~(isa(prefered_solver,'struct') | isa(prefered_solver,'char'))
        error('Argument should be a solver tag, or a sdpsettings structure');
    end
    if isa(prefered_solver,'char')
        donttest = 1;
    end
end

if ~(exist('callsedumi')==2)
    clc
    disp('The directory yalmip/solvers is not in your path.')
    disp('These must be in path:')
    disp('                      yalmip/');
    disp('                      yalmip/extras');
    disp('                      yalmip/operators');
    disp('                      yalmip/modules');
    disp('                      yalmip/solvers');
    disp('See installation guide')
    return
end

foundstring = {'not found','found','internal'};
teststring = {'-failed','+passed'};
if ~donttest
    header = {'Solver','Version','Status','Unit test'};
else
    header = {'Solver','Version','Status'};
end

[solvers,found] = getavailablesolvers(0);
j = 1;
status = ones(length(solvers),1);
s = {solvers.tag};
for i = 1:length(solvers)
    if solvers(i).show
        same = find(strcmpi(solvers(i).tag,s));
        % Find all instances of same solver different versions
        for k = setdiff(same,i)
            % No reason to show for all versions
            solvers(k).show = 0;
        end

        data{j,1} = upper(solvers(i).tag);
        found_versions = same(find(found(same)));
        if ~isempty(found_versions)
        	idx = min(found_versions);
            version = solvers(idx).version;
            if ~any(strcmpi(version,{'geometric','standard'}))
                data{j,2} = [solvers(idx).version ' ' solvers(idx).subversion];
            end
        else
        	idx = i;
        end
        status(j) = found(idx)+1+solvers(idx).builtin;
        data{j,3} = foundstring{found(idx)+1+solvers(idx).builtin};
        j = j+1;
    end
end

if isa(prefered_solver,'char')
    ops = sdpsettings('Solver',prefered_solver);
else
    ops = prefered_solver;
end
ops.saveyalmipmodel = 1;

if ~((nargin==2) & (ops.verbose==0))
    [sortedName,loc] = sort({data{:,1}});
    loc = [loc(find(status(loc)==3)) loc(find(status(loc)==2)) loc(find(status(loc)==1))];
    dataSorted = reshape({data{loc,:}},[],3);
    yalmiptable({'Searching for installed solvers'},header,dataSorted);
    disp(' ')
end
if nargin<2
    disp('Press any key to continue test')
    pause
end

i=1;
test{i}.fcn  = 'test_core';
test{i}.desc = 'Core functionalities';
i = i+1;

test{i}.fcn  = 'test_linear_programming';
test{i}.desc = 'Linear programming (LP)';
i = i+1;

test{i}.fcn  = 'test_quadratic_programming';
test{i}.desc = 'Quadratic programming (QP)';
i = i+1;

test{i}.fcn  = 'test_socp_programming';
test{i}.desc = 'Second-order cone programming (SOCP)';
i = i+1;

test{i}.fcn  = 'test_semidefinite_programming';
test{i}.desc = 'Semidefinite programming (SDP)';
i = i+1;

test{i}.fcn  = 'test_geometric_programming';
test{i}.desc = 'Geometric programming (GP)';
i = i+1;

test{i}.fcn  = 'test_nonlinear_programming';
test{i}.desc = 'Nonlinear programming (NLP)';
i = i+1;

test{i}.fcn  = 'test_nonlinear_semidefinite_programming';
test{i}.desc = 'Nonlinear SDP (NLSDP)';
i = i+1;

test{i}.fcn  = 'test_exponential_cone_programming';
test{i}.desc = 'Exponential cone programming (ECP)';
i = i+1;

test{i}.fcn  = 'test_milinear_programming';
test{i}.desc = 'Mixed-integer LP (MIQP)';
i = i+1;

test{i}.fcn  = 'test_miquadratic_programming';
test{i}.desc = 'Mixed-integer QP (MIQP)';
i = i+1;

test{i}.fcn  = 'test_misocp_programming';
test{i}.desc = 'Mixed-integer SOCP (MISOCP)';
i = i+1;

test{i}.fcn  = 'test_nonconvex_quadratic_programming';
test{i}.desc = 'Global nonconvex quadratic programming';
i = i+1;

test{i}.fcn  = 'test_nonconvex_global_programming';
test{i}.desc = 'Global nonconvex programming';
i = i+1;

pass_strings = {'Error','Passed','Solver not available'};

% Run test-problems
for i = 1:length(test)
    try
        if ops.verbose
            disp(' ');
            disp(['Testing function ' test{i}.fcn]);
            disp(' ');
        end
        % First make call to figure out solver
        info = eval([test{i}.fcn '(ops)']);
        if ~OctaveRunning
            sols{i} = addLink(upper(cleanversion(info.yalmipmodel.solver.tag)));
        else
            sols{i} = cleanversion(info.yalmipmodel.solver.tag);
        end
        % And now run all the tests
        info = eval([test{i}.fcn '(ops)']);
        pass(i) = 1;
        results{i}='Success';
    catch
        pass(i) = 0;
        results{i} = 'Failed';
        sols{i} = '';
    end
end

clear data;
header = {'Test','Status', 'Solver'};
for i = 1:length(pass)
    data{i,1} = test{i}.desc;
    data{i,2} = results{i};
    data{i,3} = sols{i};
end
if ops.verbose
    disp(' ');
end
formats{1}.data.just = 'right';
formats{2}.data.just = 'right';
formats{3}.data.just = 'right';

formats{1}.header.just = 'right';
formats{2}.header.just = 'right';
formats{3}.header.just = 'right';

clc
yalmiptable([],header,data,formats)

% Test if any LMI solver is installed.
x = sdpvar(2);[p,aux1,aux2,m] = export(x>=0,[],[],[],[],0);
if ~isempty(m)
    only_lmilab = strcmpi(m.solver.tag,'lmilab');
    only_fmincon = strcmpi(m.solver.tag,'fmincon-standard');
else
    only_lmilab = 0;
    only_fmincon = 0;
end
if isempty(m)
    disp('You do not have any LMI solver installed')
    disp(' If you intend to solve LMIs you must install a solver.')
else
    if only_lmilab
        disp('You do not have any good LMI solver installed')
        disp(' (only found LMILAB which should be avoided in YALMIP).')
        disp('If you intend to solve LMIs, please install a better solver.')
    elseif only_fmincon
        disp('You do not have any LMI solver installed')
        disp(' (YALMIP will use a nonlinear solver which cannot be expected to work)')
        disp(' If you intend to solve LMIs you must install a solver.')
    end
end
x = binvar(1);[p,aux1,aux2,m] = export(x>=0,x,[],[],[],0);
if isempty(m)
    disp('You do not have any LP/MILP solver installed')
    disp(' If you intend to solve LPs/MILPs, you have to install one.')
else
    only_bnb = strcmpi(m.solver.tag,'bnb');
    if only_bnb
        disp('You do not have any MILP solver installed')
        disp(' (only found internal BNB).')
        disp(' If you intend to solve MILP, please install a better solver.')
    end
end
x = binvar(1);[p,aux1,aux2,m] = export(x>=0,x^2,[],[],[],0);
if isempty(m)
    disp('You do not have any QP/SOCP/MIQP/MISOCP solver installed')
    disp(' If you intend to solve QP/SOCP/MIQP/MISOCP you must install solver.')
else
    only_bnb = strcmpi(m.solver.tag,'bnb');
    if only_bnb
        disp('You do not have any MIQP/MISOCP solver installed (only found internal BNB)')
        disp(' If you intend to solve MIQP/MISOCP, please install a better solver.')
    end
end
disp('See guide on interfaced solvers')

function info = test_core(ops)
if nargout > 0
    % Fake
    info.yalmipmodel.solver.tag = '';
    return
end
% Test the sdpvar implementation
pass = 1;
sol.info = yalmiperror(0,'YALMIP');
try
    x = sdpvar(2,2);
    x = sdpvar(2,2,'symmetric');
    x = sdpvar(2,2,'full');
    x = sdpvar(2,2,'toeplitz');
    x = sdpvar(2,2,'hankel');
    x = sdpvar(2,2,'skew');
    if ~ishermitian(sdpvar(2,2,'hermitian','complex'))
        error('bug')
    end
    if ~issymmetric(sdpvar(2,2,'symmetric','complex'))
        error('bug')
    end
    if ~isreal(real(sdpvar(2,2,'symmetric','complex')))
        error('bug')
    end
    if isreal(sqrt(-1)*real(sdpvar(2,2,'symmetric','complex')))
        error('bug')
    end
    x = sdpvar(2,1,'','co');
    if ~isreal(x'*x)
        error('bug')
    end
    x = sdpvar(2,2,'','co');
    if ~isreal(diag(x'*x))
        error('bug')
    end
    x = sdpvar(1,1);
    y = sdpvar(2,2);
    x*eye(2);
    eye(2)*x;
    y*3;
    3*y;
    x = sdpvar(2,3);
    y = sdpvar(2,3);
    assign(x,randn(2,3));
    z = replace(x,x(1,1:2),[8 9]);
    z = x+y;
    z = x-y;
    z = x+1;
    z = x-1;
    z = x+ones(2,3);
    z = x-ones(2,3);
    z = ones(2,3)-x;
    z = ones(2,3)-x;
    z = eye(2)*x;
    z = x*eye(3);
    z = diag(x);
    z = trace(x(1:2,1:2));
    z = diff(x);
    z = fliplr(x);
    z = flipud(x);
    z = kron(x,eye(3));
    z = kron(eye(3),x);
    z = rot90(x);
    z = sum(x);
    z = diff(x);
    z = x';
    z = x.';
    z = tril(x);
    z = triu(x);
    z = [x y];
    z = [x;y];
    sdpvar x y
    diag([x y])*[x^-1;y^-1];
    assert(isequal([x x;x x]*x-[x x;x x].*x,zeros(2)))
    assert(isequal(trace([x x;x x]*[x y;y x])-(x*x+x*y+y*x+x*x),0))

    % Regression ??
    yalmip('clear')
    sdpvar x

    (1+x+x^4)*(1-x^2);

    % Regression complex multiplcation
    A = randn(10,5)+sqrt(-1)*randn(10,5);
    b = randn(10,1)+sqrt(-1)*randn(10,1);
    x = sdpvar(5,1);
    res = A*x-b;
    assert(nnz(clean([res res]'*[res res]-res'*res,1e-8))==0)
    assert(isreal(clean(res'*res,1e-8)))

    assert(isreal(x*x'))

    result = 'N/A';
catch
    sol.info = 'Problems';
    result = 'N/A';
    pass = 0;
end

function sol = test_semidefinite_programming(ops)
%%
t = sdpvar(1,1);
Y = sdpvar(2,2);
F = [Y<=t*eye(2), Y>=[1 0.2;0.2 1]];
sol = optimize(F,t,ops);
if nargout > 0
    return
end
assert(ismember(sol.problem,[0 3 4 5]) && abs(value(t)-1.2)<=1e-4);
%%
A = [1 0;0.4 1];
B = [0.4;0.08];
L = [1.9034 1.1501];
Y = sdpvar(2,2);
F = [Y Y*(A-B*L)';(A-B*L)*Y Y]>=0;
F = F+[L*Y*L'<=1];
sol = optimize(F,-logdet(Y),ops);
Y = double(Y);
assert(ismember(sol.problem,[0 3 4 5]) && norm(value(Y)-[2.9957 -4.1514;-4.1514 6.2918])<=1e-3)
%
x = sdpvar(1,1);
y = sdpvar(1,1);
z = sdpvar(1,1);
X = [[x 1 2];[1 y 3];[2 3 100]];
F = [X>=0,x>=10,y>=0,z>=0, x<=1000, y<=1000,z<=1000];
sol = optimize(F,x+y+z,ops);
x   = double(x);
y   = double(y);
z   = double(z);
assert(ismember(sol.problem,[0 3 4 5]) && norm(value([x;y;z])-[10;0.1787;0])<=1e-3)
%%
x = sdpvar(1,1);
z = sdpvar(1,1);
X = [[x 2];[2 z]];
F = [X>=0, x>=0,z>=0,x<=10,z<=10];
sol = optimize(F,x+z,ops);
assert(~sol.problem)
%%
% Upper bound on maxcut of a n-cycle
n = 15;
Q = zeros(n);
for i = 1:n-1
    Q(i,i+1) = 1;Q(i+1,i)  = 1;
end
Q(n,1) = 1;Q(1,n) = 1;
Q = 0.25*(diag(Q*ones(n,1))-Q);
t = sdpvar(1,1);
tau = sdpvar(n,1);
F = t>=0;
M = [[-Q zeros(n,1)];[zeros(1,n) t]];
for i = 1:n
    ei = zeros(n,1);ei(i,1) = 1;
    M = M+tau(i)*[ei*ei' zeros(n,1);zeros(1,n) -1];
end
F = F+[M>=0];
sol = optimize(F,t,ops);
t   = double(t);
tau = double(t);
pass = ismember(sol.problem,[0 3 4 5]);
assert(sol.problem == 0 && abs(value(t)-14.8361)<=1e-3)
%%
n = 5;
P = magic(n);
Z = sdpvar(n,n,'toeplitz');
t = sdpvar(n,n,'full');
F = (P(:)-Z(:)<=t(:))+(P(:)-Z(:)>=-t(:));
sol = optimize(F,sum(sum(t)),ops);
assert(~sol.problem)
assert(abs(value(sum(sum(t)))-156)<=1e-3)
%%
n = 5;
P = magic(n);
Z = sdpvar(n,n,'toeplitz');
t = sdpvar(n,n,'full');
resid = P-Z;resid = resid(:);
sol = optimize([],resid'*resid,ops);
assert(~sol.problem)
assert(abs(value(resid'*resid)-1300)<=1e-2)
%%
A = magic(6);
A = A*A';
P = sdpvar(6,6);
sol = optimize((A'*P+P*A <= -P) + (P>=eye(6)),trace(P),ops);
assert(sol.problem==1)

function sol = test_linear_programming(ops)
%%
N = 5;
A = [2 -1;1 0];
B = [1;0];
C = [0.5 0.5];
[H,S] = create_CHS(A,B,C,N);
x = [2;0];
t = sdpvar(2*N,1);
U = sdpvar(N,1);
Y = H*x+S*U;
F = (U<=1)+(U>=-1);
F = F+(Y(N)>=-1);
F = F+(Y(N)<=1);
F = F+([Y;U]<=t)+([Y;U]>=-t);
sol = optimize(F,sum(t),ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(sum(t))-12.66666)<=1e-3)
%%
t = sdpvar(1,1);
F = [t>=0, t<=-10];
sol = optimize(F,t,ops);
assert(sol.problem~=0)

function sol = test_socp_programming(ops)
x = sdpvar(2,1);
a = [0;1];
b = [1;1];
F = norm(x-a)<=1;
F = F+[norm(x-b) <= 1];
sol = optimize(F,sum(x),ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(sum(x))-0.58578)<=1e-3)
z = sdpvar(3,1);
x = sdpvar(3,1);
y = sdpvar(3,1);
a = [0;1;0];
b = [1;1;0];
F = norm(x-a)<=1;
F = F+[norm(x-b)<=1];
F = F+[x(1)==0.35];
F = F+[z(2:3)==[5;6]];
sol = optimize(F,sum(x),ops);
assert(~sol.problem)
x = double(x);
y = double(y);
z = double(z);
assert(abs(sum(x)-0.27592)<=1e-3);
z = sdpvar(2,1);
x = sdpvar(2,1);
y = sdpvar(3,1);
a = [0;1];
b = [1;1];
F = norm(x-a)<=1;
F = F+[norm(x-b)<=1];
F = F+[x(1)==0.35];
F = F+[z(1,end)>=5];
F = F+[z(2,end)<=100];
F = F+[z(2)==5];
sol = optimize(F,sum(x),ops);
assert(~sol.problem)
x = double(x);
assert(abs(sum(x)-0.59)<=1e-3)

function sol = test_misdp_programming(ops)
x = intvar(4,1);
e = magic(4)*x-1;
sdpvar t
obj = t;
sol = optimize([t e';e eye(4)]>=0,obj,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-4)<=1e-3)

function sol = test_misocp_programming(ops)
x = intvar(4,1);
obj = norm(magic(4)*x-1,2);
sol = optimize([-5 <= x <= 5],obj,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-2)<=1e-3)

function sol = test_miquadratic_programming(ops)
x = intvar(4,1);
obj = norm(magic(4)*x-1,2)^2
sol = optimize([-5 <= x <= 5],obj,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-4)<=1e-3)

function sol = test_milinear_programming(ops)
x = intvar(4,1);
obj = norm(magic(4)*x-1,1)
sol = optimize([-5 <= x <= 5],obj,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-2)<=1e-3)

function sol = test_quadratic_programming(ops)
x = sdpvar(10,1);
sol = optimize([sum(x)==2, -1 <= x <= 1],x'*x,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-.4)<=1e-3)

function sol = test_nonconvex_quadratic_programming(ops)
x = sdpvar(10,1);
ops.forceglobal = 1;
sol = optimize([sum(x)==2, -1 <= x <= 1],-x'*x,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-10)<=1e-3)

function sol = test_nonconvex_global_programming(ops)
x = sdpvar(3,1);
ops.forceglobal = 1;
sol = optimize([sum(x.^3)==2, -1 <= x <= 1],-x'*x,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(obj)-2.28)<=1e-1)

function sol = test_nonlinear_semidefinite_programming(ops)
%%
A = [-1 2;-3 -4];
P = sdpvar(2,2);
alpha = sdpvar(1,1);
F = (P>=eye(2))+(A'*P+P*A <= -2*alpha*P)+(alpha >= 0);
sol = optimize([F,P(:) <= 100],-alpha,ops);
if nargout > 0
    return
end
assert(~sol.problem);
assert(abs(value(alpha)-2.5)<= 1e-3)

function sol = test_geometric_programming(ops)
%%
t1 = sdpvar(1,1);
t2 = sdpvar(1,1);
t3 = sdpvar(1,1);
t = [t1 t2 t3];
obj = (40*t1^-1*t2^-0.5*t3^-1)+(20*t1*t3)+(40*t1*t2*t3);
F = ((1/3)*t1^-2*t2^-2+(4/3)*t2^0.5*t3^-1 <= 1);
F = [F, t>=0];
sol = optimize(F,obj,ops);
if nargout > 0
    return
end
assert(~sol.problem)
assert(abs(value(D)-133.1371)<=1e-3)

function sol = test_nonlinear_programming(ops)
sdpvar x y
sol = optimize(x^2 + x^4 + exp(x) <= 1, x^2+y^2,ops)
if nargout > 0
    return
end
assert(sol.problem == 0);

function sol = test_exponential_cone_programming(ops)
sdpvar x y z
sol = optimize([expcone([x;2;z]),x==1],z,ops)
if nargout > 0
    return
end
assert(sol.problem == 0);

function html = addLink(x)

if length(x)>0
    html = ['' upper(x) ''];
else
    html = '';
end

function x = cleanversion(x)

s = strfind(x,'-');
if ~isempty(s)
    x=x(1:s-1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function rootLayer = optimizeGraph(rootLayer)
%OPTIMIZEGRAPH Do graph optimizations, such as merging redundant vl_nnwsum

  % find all weighted sum layers (vl_nnwsum)
  wsums = rootLayer.find(@vl_nnwsum) ;

  % decide whether each wsum can be merged or not (unless explicitly set by
  % user, with layer.optimize = true/false). if a wsum is used as input to
  % a *single* other wsum, then it can be merged; if it is used as input to
  % any other layer, or more than 1 wsum, then it cannot be merged.
  allLayers = rootLayer.find() ;

  for i = 1:numel(allLayers)
    inputs = allLayers{i}.inputs ;

    if isequal(allLayers{i}.func, @vl_nnwsum)
      % a wsum, flag its input wsums; if it was flagged before, then it
      % cannot be merged (because it's used as input to 2 or more wsums).
      for j = 1:numel(inputs)
        in = inputs{j} ;
        if isa(in, 'Layer') && isequal(in.func, @vl_nnwsum) && ~islogical(in.optimize)
          if isempty(in.optimize)
            in.optimize = 'flag' ;
          elseif isequal(in.optimize, 'flag')
            in.optimize = false ;
          end
        end
      end
    else
      % not a wsum, its input wsums cannot be merged unless set explicitly
      for j = 1:numel(inputs)
        in = inputs{j} ;
        if isa(in, 'Layer') && isequal(in.func, @vl_nnwsum) && ~islogical(in.optimize)
          in.optimize = false ;
        end
      end
    end
  end

  % remove flags
  for i = 1:numel(wsums)
    if isequal(wsums{i}.optimize, 'flag')
      wsums{i}.optimize = true ;
    end
  end

  % now optimize each wsum, by merging any of its valid wsum inputs
  for i = 1:numel(wsums)
    inputs = wsums{i}.inputs ;

    % make sure there's a 'weights' property at the end, with correct size
    assert(strcmp(inputs{end-1}, 'weights') && ...
      numel(inputs{end}) == numel(inputs) - 2) ;

    % separate inputs to the sum, and weights
    origWeights = inputs{end} ;
    inputs = inputs(1:end-2) ;
    weights = cell(size(inputs)) ;

    for k = 1:numel(inputs)
      in = inputs{k} ;
      if isa(in, 'Layer') && isequal(in.func, @vl_nnwsum) && in.optimize
        % merge weights and store results
        inputs{k} = in.inputs(1:end-2) ;
        weights{k} = origWeights(k) * in.inputs{end} ;
      else
        % any other input (Layer or constant), wrap it in a single cell
        inputs{k} = {in} ;
        weights{k} = origWeights(k) ;
      end
    end

    % merge the results in order
    inputs = [inputs{:}] ;
    weights = [weights{:}] ;

    % store the merged inputs list, and the weights as a name-value pair
    wsums{i}.enableCycleChecks = false ;  % faster
    wsums{i}.inputs = [inputs, {'weights', weights}] ;
    wsums{i}.enableCycleChecks = true ;
  end

end

repo: MFT


text:
% pop_rejkurt() - rejection of artifact in a dataset using kurtosis
%                 of activity (i.e. to detect peaky distribution of
%                 activity).
%
% Usage:
%   >> pop_rejkurt( INEEG, typerej) % pop-up interative window mode
%   >> [OUTEEG, locthresh, globthresh, nrej] = ...
%		= pop_rejkurt( INEEG, typerej, elec_comp, ...
%                   locthresh, globthresh, superpose, reject, vistype);
%
% Graphical interface:
%   "Electrode|Component" - [edit box] electrodes or components indices to take
%                 into consideration for rejection. Same as the 'elec_comp'
%                 parameter from the command line.
%   "Single-channel limit|Single-component limit" - [edit box] kurtosis limit
%                  in terms of standard-dev. Same as 'locthresh' command line
%                 parameter.
%   "All-channel limit|All-component limit" - [edit box] kurtosis limit
%                  in terms of standard-dev (all channel regrouped). Same as
%                 'globthresh' command line parameter.
%   "visualization type" - [popup menu] can be either 'REJECTRIALS'|'EEGPLOT'.
%                 This correspond to the command line input option 'vistype'.
%   "Display with previous rejection(s)" - [checkbox] This checkbox set the
%                 command line input option 'superpose'.
%   "Reject marked trial(s)" - [checkbox] This checkbox set the command
%                  line input option 'reject'.
%
% Inputs:
%   INEEG      - Input dataset
%   typerej    - Type of rejection (0 = independent components; 1 = eeg
%              data). Default is 1. For independent components, before
%              thresholding, the activity is normalized for each
%              component.
%   elec_comp  - [e1 e2 ...] electrodes (number) to take into
%              consideration for rejection.
%   locthresh  - Activity kurtosis limit in terms of standard-dev.
%   globthresh - Global limit (for all channel). Same unit as above.
%   superpose  - [0] do not superpose pre-labelling with previous
%              pre-labelling (stored in the dataset). [1] consider
%              both pre-labelling (using different colors). Default is [0].
%   reject     - [0] do not reject labelled trials (but still
%              store the labels. [1] reject labelled trials.
%              Default is [1].
%   vistype    - Visualization type. [0] calls rejstatepoch() and [1] calls
%              eegplot() default is [0].When added to the command line
%              call it will not display the plots if the option 'plotflag'
%              is not set.
%   topcommand - [] Deprecated argument , keep to ensure backward compatibility
%   plotflag   - [1,0] [1]Turns plots 'on' from command line, [0] off.
%              (Note for developers: When called from command line
%              it will make 'calldisp = plotflag') {Default: 0}
%
% Outputs:
%   OUTEEG     - output dataset with updated kurtosis array
%   locthresh  - electrodes probability of activity thresholds in terms
%              of standard-dev.
%   globthresh - global threshold (where all electrode activity are
%              regrouped).
%   nrej       - number of rejected sweeps
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: rejkurt(), rejstatepoch(), pop_rejepoch(), eegplot(), eeglab()

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

% 01-25-02 reformated help & license -ad
% 03-07-02 added srate argument to eegplot call -ad
% 03-08-02 add eeglab options -ad

function [EEG, locthresh, globthresh, nrej, com] = pop_rejkurt( EEG, icacomp, elecrange, ...
                       		locthresh, globthresh, superpose, reject, vistype, topcommand, plotflag);
nrej = []; com = '';
if nargin < 1
   help pop_rejkurt;
   return;
end;
if nargin < 2
   icacomp = 1;
end;
if exist('reject') ~= 1
    reject = 1;
end;
if icacomp == 0
	if isempty( EEG.icasphere )
	    ButtonName=questdlg( 'Do you want to run ICA now ?', ...
                         'Confirmation', 'NO', 'YES', 'YES');
    	switch ButtonName,
        	case 'NO', disp('Operation cancelled'); return;
        	case 'YES', [ EEG com ] = pop_runica(EEG);
    	end % switch
	end;
end;

if nargin < 3

    % which set to save
    % -----------------
    promptstr = { [ fastif(icacomp, 'Electrode', 'Component') ' (indices; Ex: 2 6:8 10):' ], ...
                  [ fastif(icacomp, 'Single-channel', 'Single-component') ' limit(s) (std. dev(s): Ex: 2 2 2 2.5):'], ...
                  [ fastif(icacomp, 'All-channel', 'All-component') ' limit(s) (std. dev(s): Ex: 2.1 2 2 2):'], ...
                  'Visualization mode',...
                  'Display previous rejection marks', ...
                  'Reject marked trial(s)'};
    inistr = { fastif(icacomp, ['1:' int2str(EEG.nbchan)], ['1:' int2str(size(EEG.icaweights,1))])...
               fastif(icacomp, '3', '5'),  ...
               fastif(icacomp, '3', '5'), ...
               ''...
               '1', ...
               '0',};
    vismodelist= {'REJECTTRIALS','EEGPLOT'};
    g1  = [1 0.1 0.75];
    g2 = [1 0.26 0.9];
    g3 = [1 0.22 0.85];
    geometry = {g1 g1 g1 g2 [1] g3 g3};

    uilist = {...
        { 'Style', 'text', 'string', promptstr{1}} {} { 'Style','edit'      , 'string' ,inistr{1} 'tag' 'cpnum'}...
        { 'Style', 'text', 'string', promptstr{2}} {} { 'Style','edit'      , 'string' ,inistr{2} 'tag' 'singlelimit'}...
        { 'Style', 'text', 'string', promptstr{3}} {} { 'Style','edit'      , 'string' ,inistr{3} 'tag' 'alllimit'}...
        { 'Style', 'text', 'string', promptstr{4}} {} { 'Style','popupmenu' , 'string' , vismodelist 'tag' 'specmethod' }...
        {}...
        { 'Style', 'text', 'string', promptstr{5}} {} { 'Style','checkbox'  ,'string'  ,' ' 'value'  str2double(inistr{5})  'tag' 'rejmarks' }...
        { 'Style', 'text', 'string', promptstr{6}} {} { 'Style','checkbox'  ,'string'  ,' ' 'value'  str2double(inistr{6})  'tag' 'rejtrials'} ...
        };
    figname = fastif(~icacomp, 'Trial rejection using comp. kurtosis -- pop_rejkurt()', 'Trial rejection using data kurtosis -- pop_rejkurt()');
    result = inputgui( geometry,uilist,'pophelp(''pop_rejkurt'');', figname);

    size_result  = size( result );
    if size_result(1) == 0, locthresh = []; globthresh = []; return; end;
    elecrange    = result{1};
    locthresh    = result{2};
    globthresh   = result{3};
    switch result{4}, case 1, vistype=0; otherwise, vistype=1; end;
    superpose    = result{5};
    reject       = result{6};
end;

if ~exist('vistype','var'),   vistype   = 0;   end;
if ~exist('reject','var'),    reject    = 0;    end;
if ~exist('superpose','var'), superpose = 1; end;

if isstr(elecrange) % convert arguments if they are in text format
    calldisp   = 1;
	elecrange  = eval( [ '[' elecrange ']'  ]  );
	locthresh  = eval( [ '[' locthresh ']'  ]  );
	globthresh = eval( [ '[' globthresh ']' ]  );
else
    calldisp = 0;
end;

if exist('plotflag','var') && ismember(plotflag,[1,0])
    calldisp = plotflag;
else
    plotflag = 0;
end

if isempty(elecrange)
	error('No electrode selectionned');
end;

% compute the joint probability
% -----------------------------
if icacomp == 1
	fprintf('Computing kurtosis for channels...\n');
    tmpdata = eeg_getdatact(EEG);
    if isempty(EEG.stats.kurtE )
		[ EEG.stats.kurtE rejE ] = rejkurt( tmpdata, locthresh, EEG.stats.kurtE, 1);
	end;
	[ tmp rejEtmp ] = rejkurt( tmpdata(elecrange, :,:), locthresh, EEG.stats.kurtE(elecrange, :), 1);
    rejE    = zeros(EEG.nbchan, size(rejEtmp,2));
	rejE(elecrange,:) = rejEtmp;

	fprintf('Computing all-channel kurtosis...\n');
	tmpdata2 = permute(tmpdata, [3 1 2]);
	tmpdata2 = reshape(tmpdata2, size(tmpdata2,1), size(tmpdata2,2)*size(tmpdata2,3));
	[ EEG.stats.kurt rej ] = rejkurt( tmpdata2, globthresh, EEG.stats.kurt, 1);
else
	fprintf('Computing joint probability for components...\n');
    % test if ICA was computed
    % ------------------------
    icaacttmp = eeg_getica(EEG);
    if isempty(EEG.stats.icakurtE )
		[ EEG.stats.icakurtE rejE ] = rejkurt( icaacttmp, locthresh, EEG.stats.icakurtE, 1);
	end;
	[ tmp rejEtmp ] = rejkurt( icaacttmp(elecrange, :,:), locthresh, EEG.stats.icakurtE(elecrange, :), 1);
	rejE    = zeros(size(icaacttmp,1), size(rejEtmp,2));
	rejE(elecrange,:) = rejEtmp;

	fprintf('Computing global joint probability...\n');
	tmpdata = permute(icaacttmp, [3 1 2]);
	tmpdata = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3));
	[ EEG.stats.icakurt rej] = rejkurt( tmpdata, globthresh, EEG.stats.icakurt, 1);
end;
rej = rej' | max(rejE, [], 1);
fprintf('%d/%d trials marked for rejection\n', sum(rej), EEG.trials);

if calldisp
	if vistype == 1 % EEGPLOT -------------------------
	    if icacomp == 1 macrorej  = 'EEG.reject.rejkurt';
	        			macrorejE = 'EEG.reject.rejkurtE';
	    else			macrorej  = 'EEG.reject.icarejkurt';
	        			macrorejE = 'EEG.reject.icarejkurtE';
	    end;

		colrej = EEG.reject.rejkurtcol;
		eeg_rejmacro; % script macro for generating command and old rejection arrays

	    if icacomp == 1
	        eegplot( tmpdata(elecrange,:,:), 'srate', ...
		      EEG.srate, 'limits', [EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:});
	    else
	        eegplot( icaacttmp(elecrange,:,:), 'srate', ...
		      EEG.srate, 'limits', [EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:});
	    end;
    else % REJECTRIALS -------------------------
	  	if icacomp	== 1
			[ rej, rejE, n, locthresh, globthresh] = ...
				rejstatepoch( tmpdata, EEG.stats.kurtE(elecrange,:), 'global', 'on', 'rejglob', EEG.stats.kurt, ...
						'threshold', locthresh, 'thresholdg', globthresh, 'normalize', 'off'  );
		else
			[ rej, rejE, n, locthresh, globthresh] = ...
				rejstatepoch( icaacttmp, EEG.stats.icakurtE(elecrange,:), 'global', 'on', 'rejglob', EEG.stats.icakurt, ...
						'threshold', locthresh, 'thresholdg', globthresh, 'normalize', 'off' );
		end;
		nrej = n;
	end;
else
	% compute rejection locally
	rejtmp = max(rejE(elecrange,:),[],1);
	rej = rejtmp | rej;
	nrej =  sum(rej);
	fprintf('%d trials marked for rejection\n', nrej);
end;
if ~isempty(rej)
	if icacomp	== 1
		EEG.reject.rejkurt = rej;
		EEG.reject.rejkurtE = rejE;
	else
		EEG.reject.icarejkurt = rej;
		EEG.reject.icarejkurtE = rejE;
	end;
    if reject
        EEG = pop_rejepoch(EEG, rej, 0);
    end;
end;
nrej = sum(rej);

com = [ com sprintf('%s = pop_rejkurt(%s,%s);', inputname(1), ...
		inputname(1), vararg2str({icacomp,elecrange,locthresh,globthresh,superpose,reject,vistype, [], plotflag})) ];
if nargin < 3 & nargout == 2
	locthresh = com;
end;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function ribbons = riverplot_ribbon_matrix(layer1, layer2, sim_matrix, varargin)
%
% This is a subfunction of the riverplot toolbox.
%  - Draws a set of ribbon FROM layer1 TO layer2
%  - layers are structures created by riverplot_draw_layer.m
%  - Ribbon covers a specified percentage of each rect structure
%  (rectangle) within each layer, can be determined by a symmetrical
%  similarity matrix or an asymmetric one (sim_matrix).
%
% :Usage:
% ::
%
%     myribbon = riverplot_ribbon_matrix(layer1, layer2, sim_matrix, [optional inputs])
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016  Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **layer1:**
%        Cell array of structures with handles and decription of imrect class objects,
%        created with riverplot_draw_layer.m.  Each object is a rectangular
%        graphic that denotes a spatial brain pattern or mask created with riverplot_rect.m.
%        layer1 is the LEFTMOST layer in the plot.
%
%   **layer2:**
%        Same as layer1, but this is the RIGHTMOST layer of the pair
%
%
%   **sim_matrix:**
%       A [j x k] rectangular matrix of similarities to plot, where
%       j is num. of elements in layer2, and k is num. elements in layer1
%       Any similarity metric of choice is OK, but all values should be between
%       0 and 1], so negative values cannot be used.
%
%
% :Optional Inputs:
%   **coveragetype:**
%        'coveragetype' followed by 'relative' 'absolute' or 'normalized'
%
%       Relative: cover 100% of each rectangle, in proportion to relative associations.
%       Thickness of ribbon is proportional to input similarity only within
%       each element, not across them
%
%       absolute coverage: Ribbon width is proportional to similarity on 0-1
%       scale, with no normalization
%
%       normalized coverage [default]: Ribbon width for layer 1 is relative to layer 1
%       so that each layer is fully covered, but on layer2 is proportional to
%       similarity across all elements, normalized by max across all
%
%       sum-normalized coverage: Ribbon width for layer 1 is relative to layer 1
%       so that each layer is fully covered, but on layer2 is proportional to
%       summed similarity across all elements, that is normalized by max
%       of the summ across all similarities for each layer2 item.
%
%   **color, colors:**
%        followed by cell array of {[r g b] [r g b]...} color spec for
%        ribbons. Colors are specified for each element of layer1.
%
%   **steepness**
%        followed by value for sigmoid steepness, on [0 1] inteval
%
%   **ax:**
%        followed by axis handle to plot
%
%   **from_bottom:**
%        align ribbon with bottom edges of rects
%
% :Outputs:
%
%   **myribbon:**
%        structure with handles and info for ribbon
%
%
% :Examples:
% ::
%
% rect1 = riverplot_rect('color', [.7 .3 .3], 'position', [1 2 1 3]);
% rect2 = riverplot_rect('color', [.3 .3 .6], 'position', [5 4 1 3]);
% myribbon = riverplot_ribbon(rect1, rect2);
% myribbon = riverplot_ribbon(rect1, rect2, 'coverage', [1 .5], 'color', [.2 .5 .2]);
% myribbon = riverplot_ribbon(rect1, rect2, 'coverage', [.5 .2], 'color', [.2 .5 .2], 'from_bottom');
%
% :References:
%   CITATION(s) HERE
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    Created by Tor Wager, July 2, 2016
%    8/21/2017 Stephan Geuter - added sum-normalization
% ..

% ..
%    DEFAULTS AND INPUTS
% ..

ax = gca;
[n2, n1] = size(sim_matrix);
mycolor = scn_standard_colors(n1);
coveragetype = 'normalized';
steepness = .05;

% arguments passed on to subfunctions:
% from_bottom

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'ax', ax = varargin{i+1}; varargin{i+1} = [];
            %case 'position', myposition = varargin{i+1}; varargin{i+1} = [];
            case {'color', 'colors'}, mycolor = varargin{i+1}; varargin{i+1} = [];

            case 'coveragetype', coveragetype = varargin{i+1}; varargin{i+1} = [];

            % arguments passed on to subfunctions: Do nothing, but these
            % are used.
            case 'from_bottom'

            case 'coverage', coverage = varargin{i+1}; varargin{i+1} = [];

            case 'steepness', steepness = varargin{i+1}; varargin{i+1} = [];

            case ' '
                % do nothing, avoid warning from calling function

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

axes(ax)

ribbons = cell(n2, n1);

if length(layer2) ~= size(sim_matrix, 1), error('layer2 size does not match sim matrix column size'); end
if length(layer1) ~= size(sim_matrix, 2), error('layer1 size does not match sim matrix row size'); end

% For determining offset and position:
my_total_layer1 = sum(abs(sim_matrix));
my_total_layer2 = sum(abs(sim_matrix'))';

my_total_layer1(my_total_layer1 == 0) = 1;  % prevent dividing by zero
my_total_layer2(my_total_layer2 == 0) = 1;  % prevent dividing by zero

% for both offsets and coverages. offsets are starting points, coverages
% are the thickness of each ribbon
% coverages should always be positive - take absolute values (SG 8/21/2017)
fullcoverages_layer1 = abs(bsxfun(@(x,y) x ./ y, sim_matrix, my_total_layer1));
fullcoverages_layer2 = abs(bsxfun(@(x,y) x ./ y, sim_matrix, my_total_layer2));

switch coveragetype

    case 'relative'
        % Relative: cover 100% of each rectangle, in proportion to relative associations.
        % Thickness of ribbon is proportional to input similarity only within
        % each element, not across them

        coverages_layer1 = fullcoverages_layer1;
        coverages_layer2 = fullcoverages_layer2;

    case 'absolute'
        % absolute coverage: Ribbon width is proportional to similarity on 0-1
        % scale, with no normalization
        coverages_layer1 = abs(sim_matrix);
        coverages_layer2 = abs(sim_matrix);

    case 'normalized'
        % normalized coverage [default]: Ribbon width for layer 1 is relative to layer 1
        % so that each layer is fully covered, but on layer2 is proportional to
        % similarity across all elements, normalized by max across all

        mymax = max(abs(sim_matrix(:)));
        if mymax <=0, mymax = 1; end % in case empty

        %mymax = max(my_total_layer2);

%         coverages_layer1 = sim_matrix ./ mymax;
        coverages_layer1 = fullcoverages_layer1;

        coverages_layer2 = abs(sim_matrix) ./ mymax;

    case 'sumnormalized'
        % sum-normalized coverage: Ribbon width for layer 1 is relative to layer 1
        % so that each layer is fully covered, but on layer2 is proportional to
        % summed similarity across all elements, that is normalized by max
        % of the summ across all similarities for each layer2 item.

        mymax = max(sum(abs(sim_matrix)));
        if mymax <=0, mymax = 1; end % in case empty

        coverages_layer1 = fullcoverages_layer1;
        coverages_layer2 = abs(sim_matrix) ./ mymax;

    otherwise
        error('Illegal string entered for coveragetype');

end

% offsets are in one direction, take abs(coverages). SG 8/21/2017
offsets_layer1 = [zeros(1, n1); cumsum(abs(fullcoverages_layer1))];
offsets_layer2 = [zeros(1, n2); cumsum(abs(fullcoverages_layer2'))]';

% where there is only one target, center offsets
wh = sum(coverages_layer2 ~= 0, 2) == 1;
mycov = sum(coverages_layer2, 2) ./ 2; % half the coverage
mycov = 0.5 - mycov;  % midpoint 0.5 - half the coverage interval
offsets_layer2(wh, :) = repmat(mycov(wh), 1, size(offsets_layer2, 2));

for i = 1:n1

    %myoffset = i .* (1 ./ (1.5*n1));  % adjust based on number of sources to reduce overlap in ribbon endings

    for j = 1:n2

        myoffset = [offsets_layer1(j, i) offsets_layer2(j, i)];

        mycoverage = [coverages_layer1(j, i) coverages_layer2(j, i)];

        for k=1:2 %cludgey fix for riverplots that extend beyond layer 2 rectangle - this should be fixed better above
           if myoffset(k)+mycoverage(k)>1 %we dont want ribbons extending beyond range.. scoot down a bit
              myoffset(k)=max(0,myoffset(k)- (myoffset(k)+mycoverage(k)-1) ); % 8/30/17 SG scoot only down to match the max of box (=1)
           end
        end

        % 8/30/17 Stephan Geuter
        % reverse offsets in order to reduce the number of ribbon
        % crossings.
        myoffset = 1-myoffset-mycoverage;

        % this would be the relative coverage way
        %mycoverage = [sim_matrix(j, i) ./ my_total_layer1(i) sim_matrix(j, i) ./ my_total_layer2(j)];

%         myoffset = j .* (1 ./ (2*n2));  % adjust based on number of targets to reduce overlap in ribbon endings

        if sim_matrix(j, i) ~= 0

            ribbons{j, i} = riverplot_ribbon(layer1{i}, layer2{j}, 'coverage', mycoverage, ...
                'y_offset', myoffset, 'color', mycolor{i}, 'from_bottom', 'steepness', steepness);

            % riverplot_remove_ribbons(ribbons)  % for debugging.
        end

    end

end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [xh] = lin_scale( xh, mNorm )

hNorm = sqrt(sum(xh.^2));

if hNorm,
    s = mNorm*1.2/hNorm;
    xh = xh.*s;
end

author: thomas-koehler

repo: SupER


text:
function ccfi_0_tabulate ( )

%*****************************************************************************80
%
%% CCFI_0_TABULATE tabulates CCFI_0 quadrature rules for the Laguerre integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CCFI_0_TABULATE\n' );
  fprintf ( 1, '  Tabulate CCFI_0 quadrature rules for the Laguerre integral.\n' );
  fprintf ( 1, '  Density function rho(x) = exp(-x).\n' );
  fprintf ( 1, '  Region: 0 <= x < +oo.\n' );
  fprintf ( 1, '  Exactness: NONE.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Print the first 5 rules:\n' );
  fprintf ( 1, '   I          X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  ell = 1.0;

  for n = 1 : 5

    [ x, w ] = ccfi_0 ( n, ell );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6f  %14.6f\n', i, x(i), w(i) );
    end
    fprintf ( 1, ' Sum                  %14.6f\n', sum ( w(1:n) ) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Print the first 5 nested rules:\n' );
  fprintf ( 1, '   I          X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  ell = 1.0;

  for j = 1 : 5

    n = 2^j - 1;

    [ x, w ] = ccfi_0 ( n, ell );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6f  %14.6f  %14.6f\n', i, x(i), w(i) );
    end
    fprintf ( 1, ' Sum                  %14.6f\n', sum ( w(1:n) ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [f] = mci_nmm_fx_delay (x,u,P,M)
% State equations for a neural mass model of erps with first order delays
% FORMAT [f] = mci_nmm_fx_delay (x,u,P,M)
%
% x      - state vector
%   x(:,1) - voltage (spiny stellate cells)
%   x(:,2) - voltage (pyramidal cells) +ve
%   x(:,3) - voltage (pyramidal cells) -ve
%   x(:,4) - current (spiny stellate cells)    depolarizing
%   x(:,5) - current (pyramidal cells)         depolarizing
%   x(:,6) - current (pyramidal cells)         hyperpolarizing
%   x(:,7) - voltage (inhibitory interneurons)
%   x(:,8) - current (inhibitory interneurons) depolarizing
%   x(:,9) - voltage (pyramidal cells)
%
% f        - dx(t)/dt  = f(x(t))
%
% Prior fixed parameter scaling [Defaults]
%
%  M.pF.E = [32 16 4];           % extrinsic rates (forward, backward, lateral)
%  M.pF.H = [1 4/5 1/4 1/4]*128; % intrinsic rates (g1, g2 g3, g4)
%  M.pF.D = [2 16];              % propogation delays (intrinsic, extrinsic)
%  M.pF.G = [4 32];              % receptor densities (excitatory, inhibitory)
%  M.pF.T = [8 16];              % synaptic constants (excitatory, inhibitory)
%  M.pF.R = [1 1/2];             % parameter of static nonlinearity
%
%__________________________________________________________________________
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: mci_nmm_fx_delay.m 6548 2015-09-11 12:39:47Z will $

if ~isstruct(P)
%     if isfield(M,'V')
%         % Project params out of subspace
%         P=M.V*P(:);
%     end
    P=spm_unvec(P,M.pE);
end

% get dimensions and configure state variables
%--------------------------------------------------------------------------
n = length(P.A{1});         % number of sources
x = spm_unvec(x,M.x);       % neuronal states

% [default] fixed parameters
%--------------------------------------------------------------------------
E = [1 1/2 1/8]*32;         % extrinsic rates (forward, backward, lateral)
G = [1 4/5 1/4 1/4]*128;    % intrinsic rates (g1 g2 g3 g4)
D = [2 16];                 % propagation delays (intrinsic, extrinsic)
H = [4 32];                 % receptor densities (excitatory, inhibitory)
T = [8 16];                 % synaptic constants (excitatory, inhibitory)
R = [2 1]/3;                % parameters of static nonlinearity

% [specified] fixed parameters
%--------------------------------------------------------------------------
if isfield(M,'pF')
    try, E = M.pF.E; end
    try, G = M.pF.H; end
    try, D = M.pF.D; end
    try, H = M.pF.G; end
    try, T = M.pF.T; end
    try, R = M.pF.R; end
end

% test for free parameters on intrinsic connections
%--------------------------------------------------------------------------
try
    G = G.*exp(P.H);
end
G     = ones(n,1)*G;

% exponential transform to ensure positivity constraints
%--------------------------------------------------------------------------
A{1}  = exp(P.A{1})*E(1);
A{2}  = exp(P.A{2})*E(2);
A{3}  = exp(P.A{3})*E(3);
C     = exp(P.C);

% intrinsic connectivity and parameters
%--------------------------------------------------------------------------
Te    = T(1)/1000*exp(P.T(:,1));         % excitatory time constants
Ti    = T(2)/1000*exp(P.T(:,2));         % inhibitory time constants
He    = H(1)*exp(P.G(:,1));              % excitatory receptor density
Hi    = H(2)*exp(P.G(:,2));              % inhibitory receptor density

% delays
%--------------------------------------------------------------------------
De = D(2).*exp(P.D)/1000;                % extrinsic
Di = D(1)/1000;                          % intrinsic

% presynaptic i/p from delayed voltages using first order Taylor approx
% v(t-tau) = v(t) - tau dv/dt
%--------------------------------------------------------------------------
pyr_int = presynaptic (x(:,9)-Di*(x(:,5)-x(:,6)),P,R);  % Intrinsic delays
stl_int = presynaptic (x(:,1)-Di*x(:,4),P,R);
inh_int = presynaptic (x(:,7)-Di*x(:,8),P,R);
pyr_ext = presynaptic (x(:,9)-De*(x(:,5)-x(:,6)),P,R);  % Extrinsic delays

% input
%==========================================================================
if isfield(M,'u')

    % endogenous input
    %----------------------------------------------------------------------
    U = u(:)*64;

else
    % exogenous input
    %----------------------------------------------------------------------
    U = C*u(:)*2;
end

% State: f(x)
%==========================================================================

%f=zeros(size(M.x));

% Supragranular layer (inhibitory interneurons): Voltage & depolarizing current
%--------------------------------------------------------------------------
f(:,7) = x(:,8);
f(:,8) = (He.*((A{2} + A{3})*pyr_ext + G(:,3).*pyr_int) - 2*x(:,8) - x(:,7)./Te)./Te;

% Granular layer (spiny stellate cells): Voltage & depolarizing current
%--------------------------------------------------------------------------
f(:,1) = x(:,4);
f(:,4) = (He.*((A{1} + A{3})*pyr_ext + G(:,1).*pyr_int + U) - 2*x(:,4) - x(:,1)./Te)./Te;

% Infra-granular layer (pyramidal cells): depolarizing current
%--------------------------------------------------------------------------
f(:,2) = x(:,5);
f(:,5) = (He.*((A{2} + A{3})*pyr_ext + G(:,2).*stl_int) - 2*x(:,5) - x(:,2)./Te)./Te;

% Infra-granular layer (pyramidal cells): hyperpolarizing current
%--------------------------------------------------------------------------
f(:,3) = x(:,6);
f(:,6) = (Hi.*G(:,4).*inh_int - 2*x(:,6) - x(:,3)./Ti)./Ti;

% Infra-granular layer (pyramidal cells): Voltage
%--------------------------------------------------------------------------
f(:,9) = x(:,5) - x(:,6);

%f      = spm_vec(f);
f=f(:);

end

function [S] = presynaptic (x,P,R)

% pre-synaptic inputs: s(V)
%--------------------------------------------------------------------------
R     = R.*exp(P.S);
S     = 1./(1 + exp(-R(1)*(x - R(2)))) - 1./(1 + exp(R(1)*R(2)));

end

author: spm


text:
function [MethodList, pathmodels] = qMRLab_static_Models
% Static listing of all Models in the qMRLab Model folder
pathmodels = {'Diffusion/' 'Diffusion/' 'Diffusion/' 'FieldMaps/' 'FieldMaps/' 'FieldMaps/' 'Magnetization_transfer/' 'Magnetization_transfer/' 'Magnetization_transfer/' 'Magnetization_transfer/' 'Magnetization_transfer/' 'Noise/' 'Noise/' 'Processing/' 'QSM/' 'T1_relaxometry/' 'T1_relaxometry/' 'T1_relaxometry/' 'T2_relaxometry/' 'T2_relaxometry/' 'UnderDevelopment/'};
MethodList = {'charmed' 'dti' 'noddi' 'b0_dem' 'b1_dam' 'b1_afi' 'mt_ratio' 'mt_sat' 'qmt_bssfp' 'qmt_sirfse' 'qmt_spgr' 'denoising_mppca' 'noise_level' 'filter_map' 'qsm_sb' 'inversion_recovery' 'vfa_t1' 'mp2rage' 'mtv' 'mwf' 'mono_t2' 'CustomExample'};

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


function [ Graph ] = getMatlabGraph( pointByPose, poseLength )
%GETMATLABGRAPH Summary of this function goes here
%   Detailed explanation goes here

    symMat = zeros(poseLength, poseLength);

    for i=1:length(pointByPose)

        combinations = combntns(pointByPose(i,:), 2);

        for c = 1:length(combinations)
           [c1,c2] = combinations(c,:);
           symMat(c1, c2) = symMat(c1, c2) + 1;
        end

    end

    % symming
    symMat = symMat + symMat';

    Graph = graph(symMat);

end


text:
function ft_save_workspace(dae1Ot8o_dir)

% FT_SAVE_WORKSPACE saves every variable in the base workspace to a .mat file with
% the same name as the variable in the workspace itself. For example, the variable
% "ans" would be saved to the file "ans.mat". Prior to calling this function, you
% might want to clean up your workspace using CLEAR or KEEP.
%
% Use as
%   ft_save_workspace(dirname)
%
% If the directory does not yet exist, this function will create it for you. If you
% leave it empty, the files will be saved to the present working directory.
%
% For example, the following will save all variables to a time-stamped
% sub-directory that is created inside the present working directory:
%
%   ft_save_workspace(datestr(now))
%
% See also SAVE, LOAD, SAVEFIG, CLEAR, KEEP

% Copyright (C) 2018, 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$

% use long and random variable names to avoid clashes with variables in the base workspace
if nargin<1
  dae1Ot8o_dir = '.';
end

if ~isdir(dae1Ot8o_dir)
  ft_info('creating directory ''%s''\n', dae1Ot8o_dir);
  mkdir(dae1Ot8o_dir);
end

% get a list of all the variables in the base workspace
dae1Ot8o_w = evalin('base', 'whos');

% copy each of the variables over to the local workspace
for dae1Ot8o_i=1:numel(dae1Ot8o_w)
  dae1Ot8o_var = dae1Ot8o_w(dae1Ot8o_i).name;
  dae1Ot8o_val = evalin('base', dae1Ot8o_var);
  eval(sprintf('%s = dae1Ot8o_val;', dae1Ot8o_var));
  dae1Ot8o_file = fullfile(dae1Ot8o_dir, [dae1Ot8o_var '.mat']);
  ft_info('saving variable ''%s'' to ''%s''\n', dae1Ot8o_var, dae1Ot8o_file);
  save(dae1Ot8o_file, dae1Ot8o_var);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ValidationConfiguratorUIElements < handle
    properties
        validationTypeRadioButtonGroup;
        holdoutRadio;
        crossValidationRadio;
        moveRightButton;
        moveLeftButton;
        trainFilesList;
        testFilesList;
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
tic;
CENTERS = {'igs_gps'; 'emr'; 'esoc'; 'code'; 'code1d'; 'code5s'; 'gfz'; 'cnes'; 'wuhan'};

state = Prj_Settings('last_settings.ini');
sta_co = GNSS_Station(); sta_co(1) = [];
sta_u2 = GNSS_Station(); sta_u2(1) = [];
for i = 1 : numel(CENTERS)
    state.selected_center = CENTERS{i};
    goGPS(state, 0);
    [res_zwd{i}, time{i}] = rec(1).getZwd();
    sta_u2(i) = rec(1);
    sta_co(i) = rec(2);
end
toc;

sta = sta_co;
fh = sta.showZtd;
ax = fh.Children(end);
axes(ax);
ylim([212.5 220]);
title(sprintf('%s of %s station \nusing different GPS orbits providers - combined engine\\fontsize{5} \n', ax.Title.String, rec(1).getMarkerName4Ch));
legend(CENTERS, 'Interpreter', 'none', 'Location', 'NorthWest');
fh = gcf; Core_UI.addBeautifyMenu(fh); Core_UI.beautifyFig(fh, 'light');

sta = sta_u2;
fh = sta.showZtd;
ax = fh.Children(end);
axes(ax);
ylim([212.5 220]);
title(sprintf('%s of %s station \nusing different GPS orbits providers - uncombined engine\\fontsize{5} \n', ax.Title.String, rec(1).getMarkerName4Ch));
legend(CENTERS, 'Interpreter', 'none', 'Location', 'NorthWest');
fh = gcf; Core_UI.addBeautifyMenu(fh); Core_UI.beautifyFig(fh, 'light');

%%
%sta_co_ca = sta_co;
%sta_u2_ca = sta_u2;

author: goGPS-Project


text:
%ChromosomeAnimation
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 7/20/2011
classdef ChromosomeAnimation < edu.stanford.covert.cell.sim.analysis.Animation
    properties (SetAccess = protected)
        title = 'Whole Cell DNA Animation'
        description = 'Animation of dynamics of DNA in a single Mycoplasma genitalium whole-cell simulation.'
        author = 'Jonathan Karr, Covert Lab, Department of Bioengineering, Stanford University'
    end

    methods
        function this = ChromosomeAnimation(varargin)
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            this = this@edu.stanford.covert.cell.sim.analysis.Animation(varargin{:});

            this.movieWidth = round((600/90 * this.movieDPI) / 2) * 2;
            this.movieHeight = round((300/90 * this.movieDPI) / 2) * 2;
        end

        function [times, states] = loadSimulationData(this)
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

            stateNames = {
                'Time'       'values'
                'Chromosome' 'polymerizedRegions'
                'Chromosome' 'monomerBoundSites'
                'Chromosome' 'complexBoundSites'
                };
            simDir = SimulationDiskUtil.getSimulation([this.simGroupId filesep num2str(this.simId)]);
            states = DiskLogger.load(simDir, stateNames, [], [], [], 'extract');
            times = states.Time.values;
        end

        function svg = drawFrame(this, timeIdx)
            chr = this.simulation.state('Chromosome');

            %open svg
            svg = [];
            svg = [svg sprintf('<?xml version="1.0" standalone="no"?>\n')];
            svg = [svg sprintf('<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">\n')];
            svg = [svg sprintf('<svg width="%d" height="%d" version="1.1" xmlns="http://www.w3.org/2000/svg" viewBox="%d %d %d %d">\n', 600, 300, 0, 0, 600, 300)];

            %white background
            svg = [svg sprintf('<rect x="%d" y="%d" width="%d" height="%d" style="fill:white; stroke:none"/>\n', ...
                0, 0, 600, 300)];

            %polymerized Regions
            [posStrnds, lens] = find(this.states.Chromosome.polymerizedRegions);
            lens = lens(posStrnds(:, 3) == timeIdx, :);
            posStrnds = posStrnds(posStrnds(:, 3) == timeIdx, 1:2);
            svg = [svg sprintf(this.drawArc(posStrnds, lens, [95 105], 'black', 2))];

            %bound proteins
            [monPosStrnds, monIdxs] = find(this.states.Chromosome.monomerBoundSites);
            [cpxPosStrnds, cpxIdxs] = find(this.states.Chromosome.complexBoundSites);
            posStrnds = [
                monPosStrnds(monPosStrnds(:, 3) == timeIdx, 1:2)
                cpxPosStrnds(cpxPosStrnds(:, 3) == timeIdx, 1:2)
                ];
            lens = [
                chr.monomerDNAFootprints(monIdxs(monPosStrnds(:, 3) == timeIdx, :))
                chr.complexDNAFootprints(cpxIdxs(cpxPosStrnds(:, 3) == timeIdx, :))
                ];
            svg = [svg sprintf(this.drawArc(posStrnds, lens, [90 110], 'red', 1))];

            %close svg
            svg = [svg sprintf('</svg>\n')];
        end

        function svg = drawArc(this, posStrnds, lens, r0, strokeColor, strokeWidth)
            chr = this.simulation.state('Chromosome');

            svg = [];
            for j = 1:size(posStrnds, 1)
                theta = posStrnds(j, 1) / chr.sequenceLen * 2 * pi - pi/2;
                dTheta = lens(j) / chr.sequenceLen * 2 * pi;
                largeArc = dTheta >= pi;
                r = r0(iseven(posStrnds(j, 2)) + 1);
                if posStrnds(j, 2) <= 2
                    cx = 50 + 100;
                else
                    cx = 3*50 + 3*100;
                end
                cy = 50 + 100;
                if lens(j) == chr.sequenceLen
                    svg = [svg sprintf('<circle cx="%d" cy="%d" r="%d" stroke="%s" stroke-width="%f" fill="none"/>\n', ...
                        cx, cy, r, ...
                        strokeColor, strokeWidth)]; %#ok<*AGROW>
                else
                    svg = [svg sprintf('<path d="M%f,%f A%f,%f %d %d,%d %f %f" stroke="%s" stroke-width="%d" fill="none"/>\n', ...
                        cx + r*cos(theta), cy + r*sin(theta), ...
                        r, r, 0, largeArc, 1, cx + r*cos(theta+dTheta), cy + r*sin(theta+dTheta), ...
                        strokeColor, strokeWidth)];
                end
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function pgma_write_test ( file_out_name )

%*****************************************************************************80
%
%% PGMA_WRITE_TEST tests the ASCII PGM write routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_OUT_NAME, the name of the output file.
%
  nrow = 300;
  ncol = 300;
%
%  Set the data.
%
  g = pgma_example ( nrow, ncol );
%
%  Write the data to the file.
%
  pgma_write ( file_out_name, g );

  return
end


text:
% msTrain

LOAD = 0;

nsegments = [5 7 10 15 20 25 35 50 60 70 80 90 100 Inf]; % number of segments per segmentation
ncv = 1; % number of cross-validation sets (typically 1)
labeltol = 0.9; % required percentage of single-label pixels for segment to be good (default 0.9)
nclasses = 23;  % number of categories

imdir = '~/data/msrc/MSRC_ObjCategImageDatabase_v2/Images/';  % where the training/test images are stored
datadir = '~/data/eccv08/msrc21/';  % where existing data is stored
outdir = datadir;  % where to save data, results, etc.

if LOAD
    load(fullfile(datadir, 'trainTestFn.mat'));
    load(fullfile(datadir, 'msrc_imsegs.mat'));
    load(fullfile(datadir, 'trainLabels.mat'));
    load(fullfile(outdir, 'msSuperpixelLabels.mat'));
    load(fullfile(outdir, 'msSuperpixelData.mat'));
    load(fullfile(outdir, 'msEdgeData.mat'));
    load(fullfile(outdir, 'msEdgeClassifier.mat'));
end

% assume that imsegs structures have been created
% See msCreateSuperpixels

% assume that "train" and "test" have been loaded, which include
% indices into imsegs for training and testing data
allind = union(train, test);
trainind1 = train(1:50);  % use first 50 training images to train segmentations
trainind = train;
trainind2 = setdiff(trainind, trainind1);
testind = test;

nimages = numel(imsegs);

% labels is a cell array of pixel maps, with values equal to the category (0 for unlabeled)
disp('Converting labels to superpixels')
if ~exist('splabels', 'var')
    splabels = msLabelMap2Sp(labels, {imsegs.segimage});
    save(fullfile(outdir, 'msSuperpixelLabels.mat'), 'splabels');
end

disp('Getting superpixel features')
if ~exist('spfeatures', 'var')
    spfeatures = mcmcGetAllSuperpixelData(imdir, imsegs);
    save(fullfile(outdir, 'msSuperpixelData.mat'), 'spfeatures');
end

disp('Getting edge features')
if ~exist('efeatures', 'var')
    [efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imsegs);
    save(fullfile(outdir, 'msEdgeData.mat'), 'efeatures', 'adjlist');
end

disp('Training edge classifier')
if ~exist('eclassifier', 'var')
    eclassifier = msTrainEdgeClassifier(efeatures(trainind1), ...
        adjlist(trainind1), imsegs(trainind1), splabels(trainind1));
    ecal = msCalibrateEdgeClassifier(efeatures(trainind2), adjlist(trainind2), ...
        imsegs(trainind2), eclassifier, splabels(trainind2), ncv);
    ecal = ecal{1};
    save(fullfile(outdir, 'msEdgeClassifier.mat'), 'eclassifier', 'ecal');
end

disp('Getting multiple segmentations')
if ~exist('smaps', 'var')
    for f = allind
        pE{f} = test_boosted_dt_mc(eclassifier, efeatures{f});
        pE{f} = 1 ./ (1+exp(ecal(1)*pE{f}+ecal(2)));
        smaps{f} = msCreateMultipleSegmentations(pE{f}, adjlist{f}, ...
            imsegs(f).nseg, nsegments);
    end
    save(fullfile(outdir, 'msMultipleSegmentations.mat'), 'pE', 'smaps');
end

disp('Getting segment features')
if ~exist('segfeatures', 'var')
    for f = allind
        if mod(f, 50)==0, disp(num2str(f)), end
        im = im2double(imread([imdir '/' imsegs(f).imname]));
        imdata = mcmcComputeImageData(im, imsegs(f));

        for k = 1:numel(nsegments)
            segfeatures{f, k} = mcmcGetSegmentFeatures(imsegs(f), ...
                spfeatures{f}, imdata, smaps{f}(:, k), (1:max(smaps{f}(:, k))));
            [mclab{f, k}, mcprc{f, k}, allprc{f, k}, trainw{f, k}] = ...
                msSegs2labels(imsegs(f), smaps{f}(:, k), splabels{f}, nclasses);
            seglabel{f, k} = mclab{f, k}.*(mcprc{f, k}>labeltol);
            seggood{f,k} =  1*(mcprc{f, k}>labeltol) + (-1)*(mcprc{f, k}<labeltol);
        end
    end
    save(fullfile(outdir, 'msSegmentFeatures.mat'), 'segfeatures', 'seglabel', 'seggood', 'trainw');
end

disp('Training segmentation classifier')
if ~exist('segclassifier', 'var')
    traindata = segfeatures(trainind, :);
    trainlabels = seggood(trainind, :);
    trainweights = trainw(trainind, :);
    segclassifier = mcmcTrainSegmentationClassifier2(traindata, trainlabels, trainweights);
    save(fullfile(outdir, 'msSegmentationClassifier.mat'), 'segclassifier');
end

disp('Training label classifier')
if ~exist('labelclassifier', 'var')
    traindata = segfeatures(trainind, :);
    trainlabels = seglabel(trainind, :);
    trainweights = trainw(trainind, :);
    labelclassifier = msTrainLabelClassifier(traindata, trainlabels, trainweights, classnames, Inf);
    %labelclassifier = mcmcTrainSegmentClassifier2(traindata, trainlabels, trainweights);
    save(fullfile(outdir, 'msLabelClassifier.mat'), 'labelclassifier');
end

save(fullfile(outdir, 'msClassifiers.mat'), 'eclassifier', 'segclassifier', 'labelclassifier', 'ecal');

author: Cloud-CV

repo: object-proposals


text:
% Figure 8.11: Approximate linear discrimination via support vector classifier
% Section 8.6.1, Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 10/16/05
% (a figure is generated)
%
% The goal is to find a function f(x) = a'*x - b that classifies the non-
% separable points {x_1,...,x_N} and {y_1,...,y_M} by doing a trade-off
% between the number of misclassifications and the width of the separating
% slab. a and b can be obtained by solving the following problem:
%           minimize    ||a||_2 + gamma*(1'*u + 1'*v)
%               s.t.    a'*x_i - b >= 1 - u_i        for i = 1,...,N
%                       a'*y_i - b <= -(1 - v_i)     for i = 1,...,M
%                       u >= 0 and v >= 0
% where gamma gives the relative weight of the number of misclassified
% points compared to the width of the slab.

% data generation
n = 2;
randn('state',2);
N = 50; M = 50;
Y = [1.5+0.9*randn(1,0.6*N), 1.5+0.7*randn(1,0.4*N);
     2*(randn(1,0.6*N)+1), 2*(randn(1,0.4*N)-1)];
X = [-1.5+0.9*randn(1,0.6*M),  -1.5+0.7*randn(1,0.4*M);
      2*(randn(1,0.6*M)-1), 2*(randn(1,0.4*M)+1)];
T = [-1 1; 1 1];
Y = T*Y;  X = T*X;
g = 0.1;            % gamma

% Solution via CVX
cvx_begin
    variables a(n) b(1) u(N) v(M)
    minimize (norm(a) + g*(ones(1,N)*u + ones(1,M)*v))
    X'*a - b >= 1 - u;
    Y'*a - b <= -(1 - v);
    u >= 0;
    v >= 0;
cvx_end

% Displaying results
linewidth = 0.5;  % for the squares and circles
t_min = min([X(1,:),Y(1,:)]);
t_max = max([X(1,:),Y(1,:)]);
tt = linspace(t_min-1,t_max+1,100);
p = -a(1)*tt/a(2) + b/a(2);
p1 = -a(1)*tt/a(2) + (b+1)/a(2);
p2 = -a(1)*tt/a(2) + (b-1)/a(2);

graph = plot(X(1,:),X(2,:), 'o', Y(1,:), Y(2,:), 'o');
set(graph(1),'LineWidth',linewidth);
set(graph(2),'LineWidth',linewidth);
set(graph(2),'MarkerFaceColor',[0 0.5 0]);
hold on;
plot(tt,p, '-r', tt,p1, '--r', tt,p2, '--r');
axis equal
title('Approximate linear discrimination via support vector classifier');
% print -deps svc-discr2.eps

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function faure_test ( )

%*****************************************************************************80
%
%% FAURE_TEST tests the FAURE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FAURE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FAURE library.\n' );

  faure_test005 ( );
  faure_test006 ( );
  faure_test01 ( );
  faure_test02 ( );
  faure_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FAURE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [slice, sliceInd, subX, subY, subZ] = extractSlice(volume,centerX,centerY,centerZ,normX,normY,normZ,radius)

%%extractSlice extracts an arbitray slice from a volume.
% [slice, sliceInd, subX, subY, subZ] = extractSlice extracts an arbitrary
% slice given the center and normal of the slice.  The function outputs the
% intensities, indices and the subscripts of the slice base on the input volume.
% If you are familiar with IDL, this is the equivalent function to EXTRACT_SLICE.
%
% Note that output array 'sliceInd' contains integers and NaNs if the particular
% locations are outside the volume, while output arrays 'subX, subY, subZ'
% contain floating points and does not contain NaNs when locations are
% outside the volume.
%
% Example:
% %Extracts slices from a mri volume by varying the slice orientation from
% sagittal to transverse while shifting the center of the slice from left to right.
% %(pardon me for demonstrating with an image without orientation labels and out-of-scale pixels.)
%
% load mri;
% D = squeeze(D);
% for i = 0:30;
% pt = [64 i*2 14];
% vec = [0 30-i i];
% [slice, sliceInd,subX,subY,subZ] = extractSlice(D,pt(1),pt(2),pt(3),vec(1),vec(2),vec(3),64);
% surf(subX,subY,subZ,slice,'FaceColor','texturemap','EdgeColor','none');
% axis([1 130 1 130 1 40]);
% drawnow;end;
%
%
% $Revision: 1.0 $ $Date: 2011/07/02 18:00$ $Author: Pangyu Teng $
% $Updated $ $Date: 2011/08/012 07:00$ $Author: Pangyu Teng $

if nargin < 7
   display('requires at least 7 parameters');
   return;
end

if nargin < 8,
    % sets the size for output slice radius*2+1.
    radius = 50;
end

isDebug = 0;
if isDebug,
    clear all;close all;
    isDebug = 1;
    load mri; D = squeeze(D);
    volume = D;
    pt = [63 63 24];
    vec = [0 0 1];
    radius = 30;
end

pt = [centerX,centerY,centerZ];
vec = [normX,normY,normZ];

%initialize needed parameters
%size of volume.
volSz=size(volume);
%a very small value.
epsilon = 1e-12;

%assume the slice is initially at [0,0,0] with a vector [0,0,1] and a silceSize.
hsp = surf(linspace(-radius,radius,2*radius+1),linspace(-radius,radius,2*radius+1),zeros(2*radius+1));
hspInitialVector = [0,0,1];

%normalize vectors;
hspVec = hspInitialVector/norm(hspInitialVector);
hspVec(hspVec ==0) = epsilon;
vec = vec/norm(vec);
vec(vec == 0)=epsilon;

%this does not rotate the surface, but initializes the subscript z in hsp.
rotate(hsp,[0,0,1],0);

if isDebug,
    %get the coordinates
    xdO = get(hsp,'XData');
    ydO = get(hsp,'YData');
    zdO = get(hsp,'ZData');
end

%find how much rotation is needed, below are the references.
%http://www.siggraph.org/education/materials/HyperGraph/modeling/mod_tran/3drota.htm
%http://www.mathworks.com/help/toolbox/sl3d/vrrotvec.html
hspVecXvec = cross(hspVec, vec)/norm(cross(hspVec, vec));
acosineVal = acos(dot(hspVec, vec));

%help prevents errors (i.e. if vec is same as hspVec),
hspVecXvec(isnan(hspVecXvec)) = epsilon;
acosineVal(isnan(acosineVal)) = epsilon;

%rotate to the requested orientation
rotate(hsp,hspVecXvec(:)',180*acosineVal/pi);

%get the coordinates
xd = get(hsp,'XData');
yd = get(hsp,'YData');
zd = get(hsp,'ZData');

%translate;
subX = xd + pt(1);
subY = yd + pt(2);
subZ = zd + pt(3);

%round the subscripts to obtain its corresponding values and indices in the volume.
xd = round(subX);
yd = round(subY);
zd = round(subZ);

%delete the surf
delete(hsp);

%obtain the requested slice intensitis and indices from the input volume.
xdSz=size(xd);
sliceInd=ones(xdSz)*NaN;
slice=ones(xdSz)*NaN;
for i = 1:xdSz(1)
    for j = 1:xdSz(2)
        if xd(i,j) > 0 && xd(i,j)<= volSz(1) &&...
             yd(i,j) > 0 && yd(i,j)<= volSz(2) &&...
             zd(i,j) > 0 && zd(i,j)<= volSz(3),
         sliceInd(i,j) = sub2ind(volSz,xd(i,j),yd(i,j),zd(i,j));
         slice(i,j) = volume(xd(i,j),yd(i,j),zd(i,j));
        end
    end
end

if isDebug,
    plot3(xdO,ydO,zdO,'b+'); hold on;
    plot3(xd,yd,zd,'m.');
    figure(2);
    imagesc(slice);axis tight;axis equal;
end


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrrspwt
%TFRRSPWT Unit test for the function TFRRSPWV.

%	F. Auger - December 1995, O. Lemoine - April 1996.

N=128;

% Reality of the TFR
sig=noisecg(N);
[tfr,rtfr]=tfrrspwv(sig);
if sum(any(abs(imag(rtfr))>sqrt(eps)))~=0,
 error('tfrrspwv test 1 failed');
end

% Energy conservation
sig=fmlin(N);
[tfr,rtfr]=tfrrspwv(sig);
Es=norm(sig)^2;
Etfr=sum(mean(rtfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrrspwv test 2 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
[tfr,rtfr]=tfrrspwv(sig);
[ik,jk]=find(abs(rtfr)>sqrt(eps));
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrrspwv test 3 failed');
end;

% frequency localization
f0=30;
sig=fmconst(N+6,f0/N);
[tfr rtfr]=tfrrspwv(sig,N/2+2,N,tftb_window(N+1,'rect'));
if any(find(rtfr>max(rtfr)/N)~=2*f0+1)|(abs(mean(rtfr)-1.0)>sqrt(eps)),
 error('tfrrspwv test 4 failed');
end;

% A RSPWVD with a Dirac time-smoothing window is a RPWVD
sig=noisecg(N);
[tfr1 rtfr1]=tfrrspwv(sig,1:N,N/2,1);
[tfr2 rtfr2]=tfrrpwv(sig,1:N,N/2);
if any(any(abs(rtfr1-rtfr2)>sqrt(eps))),
 error('tfrrspwv test 5 failed');
end;

N=111;

% Reality of the TFR
sig=noisecg(N);
[tfr,rtfr]=tfrrspwv(sig);
if sum(any(abs(imag(rtfr))>sqrt(eps)))~=0,
 error('tfrrspwv test 6 failed');
end

% Energy conservation
sig=fmlin(N);
[tfr,rtfr]=tfrrspwv(sig);
Es=norm(sig)^2;
Etfr=sum(mean(rtfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrrspwv test 7 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
[tfr,rtfr]=tfrrspwv(sig);
[ik,jk]=find(abs(rtfr)>sqrt(eps));
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrrspwv test 8 failed');
end;

% frequency localization
f0=30;
sig=fmconst(N+6,f0/N);
[tfr rtfr]=tfrrspwv(sig,round(N/2)+2,N,tftb_window(N,'rect'));
if any(find(rtfr>max(rtfr)/N)~=2*f0+1)|(abs(mean(rtfr)-1.0)>sqrt(eps)),
 error('tfrrspwv test 9 failed');
end;

% A RSPWVD with a Dirac time-smoothing window is a RPWVD
sig=noisecg(N);
[tfr1 rtfr1]=tfrrspwv(sig,1:N,round(N/2),1);
[tfr2 rtfr2]=tfrrpwv(sig,1:N,round(N/2));
if any(any(abs(rtfr1-rtfr2)>sqrt(eps))),
 error('tfrrspwv test 10 failed');
end;

author: HeLiangHIT


text:
function  [ DAG, K2Score ] = k2( LGObj,Order,u )
%Input: a set of variables x1,...xn; a given order among them, an upper limit u on the number of parents for a node; a database on x1 ... xn.
%Output: a DAG with oriented arcs, K2Score saves the maximum score upon each node.
%For i := 1 to n do
%1.  Pai(xi) = null; Ok := true;
%2. Pold := g(xi; Pai(xi));
%3. While Ok and |pai(xi)| <= u do
%3.1. Let z be the node in the set of predecessors of xi that does not belong to Pai(xi) which maximizes g(xi; Pai(xi)) + {z} );
%3.2. Pnew := g(xi; Pai(xi) + {z});
%3.3. If Pnew > Pold. Then Pold := Pnew;  Pai(xi) := Pai(xi)+{z} Else Ok := false.

% This code is written by Lowell Guangdi, Email: lowellli121@gmail.com

LG = struct( LGObj );
Dim = LG.VarNumber;
DAG = zeros( Dim );
K2Score = zeros( 1,Dim );

  for p = 2:Dim  % the first node have no parents
      Parent = zeros( Dim,1 ); OKToProcced = 1;
      Pold = -Inf;  %  the initiate stage
      while OKToProcced == 1 && sum( Parent )<= u
          %3.1.maximizes g(xi; Pai(xi)) + {z} )
          LocalMax = -Inf; LocalNode = 0;
          for q = ( p-1):( -1 ):1
              if Parent( Order( q ) )==0
                 Parent( Order( q ) ) = 1;
                 LocalScore = GClosedFun( LGObj, Order( p ), find( Parent' == 1 ) );
                 if LocalScore > LocalMax
                     LocalMax = LocalScore ;
                     LocalNode = Order( q );
                 end
                 Parent( Order( q ) ) = 0;
              end
          end
         %3.2. Pnew := g(xi; Pai(xi)+{z});
         Pnew = LocalMax;
         %3.3. If Pnew > Pold. Then Pold := Pnew;  Pai(xi) := Pai(xi)+{z}
         %Else Ok := false.
         if Pnew > Pold
             Pold = Pnew;
             Parent( LocalNode ) = 1;
         else
             OKToProcced = 0;
         end
      end
      % Save the local resutl upon node Order( p )
      K2Score( Order( p ) ) = Pold;
      DAG( :,Order( p ) ) = Parent;
  end
end