text:
function mpc = case16am
%CASE16AM  Power flow data for 15 bus distribution system from Das, et al
%   Please see CASEFORMAT for details on the case file format.
%
%   Data from ...
%       Das D, Kothari DP, Kalam A (1995) Simple and efficient method for load
%       flow solution of radial distribution networks. Int J Electr Power
%       Energy Syst 17:335-346. doi: 10.1016/0142-0615(95)00050-0
%       URL: https://doi.org/10.1016/0142-0615(95)00050-0

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 10;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [ %% (Pd and Qd are specified in kW & kVAr here, converted to MW & MVAr below)
	1	3	0	0	0	0	1	1	0	12.66	1	1	1;
	2	1	0	0	0	0	1	1	0	12.66	1	1.1	0.9;
	3	1	2000	1600	0	0	1	1	0	12.66	1	1.1	0.9;
	4	1	3000	400	0	0	1	1	0	12.66	1	1.1	0.9;
	5	1	2000	-400	0	0	1	1	0	12.66	1	1.1	0.9;
	6	1	1500	1200	0	0	1	1	0	12.66	1	1.1	0.9;
	7	1	4000	2700	0	0	1	1	0	12.66	1	1.1	0.9;
	8	1	5000	1800	0	0	1	1	0	12.66	1	1.1	0.9;
	9	1	1000	900	0	0	1	1	0	12.66	1	1.1	0.9;
	10	1	600	-500	0	0	1	1	0	12.66	1	1.1	0.9;
	11	1	4500	-1700	0	0	1	1	0	12.66	1	1.1	0.9;
	12	1	1000	900	0	0	1	1	0	12.66	1	1.1	0.9;
	13	1	1000	-1100	0	0	1	1	0	12.66	1	1.1	0.9;
	14	1	1000	900	0	0	1	1	0	12.66	1	1.1	0.9;
	15	1	2100	-800	0	0	1	1	0	12.66	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	0	0	10	-10	1	100	1	10	0	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [  %% (r and x specified in ohms here, converted to p.u. below)
	1	2	0	1e-8	0	0	0	0	0	0	1	-360	360;    %% original reactance of 0 set to 1e-8
	2	3	0.1202	0.1603	0	0	0	0	0	0	1	-360	360;
	3	4	0.1282	0.1763	0	0	0	0	0	0	1	-360	360;
	3	5	0.1442	0.2885	0	0	0	0	0	0	1	-360	360;
	5	6	0.0641	0.0641	0	0	0	0	0	0	1	-360	360;
	2	7	0.1763	0.1763	0	0	0	0	0	0	1	-360	360;
	7	8	0.1282	0.1763	0	0	0	0	0	0	1	-360	360;
	7	9	0.1763	0.1763	0	0	0	0	0	0	1	-360	360;
	8	10	0.1763	0.1763	0	0	0	0	0	0	1	-360	360;
	8	11	0.1282	0.1763	0	0	0	0	0	0	1	-360	360;
	2	12	0.1763	0.1763	0	0	0	0	0	0	1	-360	360;
	12	13	0.1442	0.1923	0	0	0	0	0	0	1	-360	360;
	12	14	0.1282	0.1763	0	0	0	0	0	0	1	-360	360;
	14	15	0.0641	0.0641	0	0	0	0	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0	20	0;
];

%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3;      %% in Volts
Sbase = mpc.baseMVA * 1e6;              %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);

%% convert loads from kW to MW
mpc.bus(:, [PD, QD]) = mpc.bus(:, [PD, QD]) / 1e3;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function LinearVariables = depends(varargin)

if nargin > 1
    LinearVariables = [];
    for i = 1:nargin
        LinearVariables_i = depends(varargin{i});
        LinearVariables = [LinearVariables;LinearVariables_i(:)];
    end
    LinearVariables = unique(LinearVariables);
    return
else
    LinearVariables = depends(lmi(varargin{1}));
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [scoresAll, labelsAll, nGTall] = assignGTmulti(pred,annolist_gt,thresh)

nJoints = 14;
% LSP to MPII format map
jidxMap = [0:5 10:15 8:9];
assert(length(jidxMap) == nJoints);

% part detections scores
scoresAll = cell(nJoints,1);
% positive / negative labels
labelsAll = cell(nJoints,1);
% number of annotated GT joints per image
nGTall = zeros(nJoints,length(annolist_gt));

for i = 1:nJoints
    scoresAll{i} = cell(length(annolist_gt),1);
    labelsAll{i} = cell(length(annolist_gt),1);
end

for imgidx = 1:length(annolist_gt)
    % distance between predicted and GT joints
    dist = inf(length(pred(imgidx).annorect),length(annolist_gt(imgidx).annorect),nJoints);
    % score of the predicted joint
    score = nan(length(pred(imgidx).annorect),nJoints);
    % body joint prediction exist
    hasPred = false(length(pred(imgidx).annorect),nJoints);
    % body joint is annotated
    hasGT = false(length(annolist_gt(imgidx).annorect),nJoints);

    % iterate over predicted poses
    for ridxPred = 1:length(pred(imgidx).annorect)
        % predicted pose
        rectPred = pred(imgidx).annorect(ridxPred);
        pointsPred = rectPred.annopoints.point;
        % iterate over GT poses
        for ridxGT = 1:length(annolist_gt(imgidx).annorect) % GT
            % GT pose
            rectGT = annolist_gt(imgidx).annorect(ridxGT);
            % compute reference distance as head size
            refDist = util_get_head_size(rectGT);
            pointsGT = [];
            if (~isempty(rectGT.annopoints))
                pointsGT = rectGT.annopoints.point;
            end
            % iterate over all possible body joints
            for i = 1:nJoints
                % predicted joint in LSP format
                ppPred = util_get_annopoint_by_id(pointsPred, jidxMap(i));
                if (~isempty(ppPred))
                    score(ridxPred,i) = ppPred.score;
                    hasPred(ridxPred,i) = true;
                end
                % GT joint in LSP format
                ppGT = util_get_annopoint_by_id(pointsGT, jidxMap(i));
                if (~isempty(ppGT))
                    hasGT(ridxGT,i) = true;
                end
                % compute distance between predicted and GT joint locations
                if (hasPred(ridxPred,i) && hasGT(ridxGT,i))
                    dist(ridxPred,ridxGT,i) = norm([ppGT.x ppGT.y] - [ppPred.x ppPred.y])/refDist;
                end
            end % joints
        end % GT poses
    end % predicted poses

    % number of annotated joints
    nGT = repmat(sum(hasGT,2)',size(hasPred,1),1);
    % compute PCKh
    match = dist <= thresh;
    pck = sum(match,3)./nGT;
    % preserve best GT match only
    [val,idx] = max(pck,[],2);
    for ridxPred = 1:length(idx)
        pck(ridxPred,setdiff(1:size(pck,2),idx(ridxPred))) = 0;
    end
    [val,predToGT] = max(pck,[],1);
    predToGT(val == 0) = 0;

    % assign predicted poses to GT poses
    for ridxPred = 1:length(pred(imgidx).annorect)
        if (ismember(ridxPred,predToGT)) % pose matches to GT
            % GT pose that matches the predicted pose
            ridxGT = find(predToGT == ridxPred);
            s = score(ridxPred,:);
            m = squeeze(match(ridxPred,ridxGT,:));
            hp = hasPred(ridxPred,:);
            idxs = find(hp);
            for i = 1:length(idxs)
                scoresAll{idxs(i)}{imgidx} = [scoresAll{idxs(i)}{imgidx};s(idxs(i))];
                labelsAll{idxs(i)}{imgidx} = [labelsAll{idxs(i)}{imgidx};m(idxs(i))];
            end
        else % no matching to GT
            s = score(ridxPred,:);
            m = false(size(match,3),1);
            hp = hasPred(ridxPred,:);
            idxs = find(hp);
            for i = 1:length(idxs)
                scoresAll{idxs(i)}{imgidx} = [scoresAll{idxs(i)}{imgidx};s(idxs(i))];
                labelsAll{idxs(i)}{imgidx} = [labelsAll{idxs(i)}{imgidx};m(idxs(i))];
            end
        end
    end % prediction assignment

    % save number of GT joints
    for ridxGT = 1:length(annolist_gt(imgidx).annorect)
        hg = hasGT(ridxGT,:);
        nGTall(:,imgidx) = nGTall(:,imgidx) + hg';
    end
end
end

author: Guanghan

repo: GNet-pose

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


text:
  function num = qsublisten(callback, varargin)

% QSUBLISTEN checks whether jobs, submitted by qsubfeval, have been
% completed. Whenever a job returns, it executes the provided callback function
% (should be a function handle), with the job ID as an input argument. Results
% can then be retrieved by calling QSUBGET. If a cell-array is provided as
% a the 'filter' option (see below), the second input argument passed to the
% callback function will be an index into this cell-array (to facilitate
% checking which job returned in the callback function).
%
% Note that this function is blocking; i.e., it only returns after a
% certain criterion has been met.
%
% Arguments can be supplied with key-value pairs:
%     maxnum      = maximum number of jobs to collect, function will return
%                   after this is reached. Default = Inf; so it is highly
%                   recommended you provide something here, since with
%                   maxnum=Inf the function will never return.
%     filter      = regular expression filter for job IDs to respond to.
%                   The default tests for jobs generated from the current
%                   MATLAB process. A cell-array of strings can be
%                   provided; in that case, exact match is required.
%     sleep       = number of seconds to sleep between checks (default=0)
%
% This function returns the number of jobs that were collected and for
% which the callback function was called.

% Copyright (C) 2012, Eelke Spaak
%
% 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$

maxnum = ft_getopt(varargin, 'maxnum', Inf);
filter = ft_getopt(varargin, 'filter', [generatesessionid '.*']);
sleep  = ft_getopt(varargin, 'sleep', 0);

if ischar(filter)
  regexpFilt = 1;
elseif iscellstr(filter)
  regexpFilt = 0;
else
  error('filter should either be a regexp string or cell-array of exact-match strings');
end

% keep track of which job IDs we have already recognized and fired the callback for
foundJobs = [];

curPwd = getcustompwd();

num = 0;
while (num < maxnum)

  files = dir();
  for k = 1:numel(files)

    % preliminary filter to get just the qsub-specific output files
    jobid = regexp(files(k).name, '^(.*)\.o.*$', 'tokens');

    if ~isempty(jobid) && ~contains(foundJobs, jobid{1}{1})
      jobid = jobid{1}{1};

      % wait until not only the stdout file exists, but also the stderr and
      % _output.mat. If we fire the callback before all three files are
      % present, a subsequent call to qsubget will fail
      outputfile   = fullfile(curPwd, sprintf('%s_output.mat', jobid));
      logerr       = fullfile(curPwd, sprintf('%s.e*', jobid));
      while ~exist(outputfile,'file') || ~isfile(logerr)
        pausejava(0.01);
      end

      if (regexpFilt && ~isempty(regexp(jobid, filter, 'once'))) || (~regexpFilt && ~isempty(find(strcmp(jobid, filter))))

        if (~regexpFilt && nargin(callback)>1)
          % also provide an index into the filter array
          callback(jobid, find(strcmp(jobid, filter)));
        else
          callback(jobid);
        end

        num = num+1;
        foundJobs = [foundJobs '|' jobid];
      end

    end
  end

  if (sleep > 0)
    pausejava(sleep);
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper function that detects a file, even with a wildcard in the filename
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function status = isfile(name)
tmp = dir(name);
status = length(tmp)==1 && ~tmp.isdir;
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function h = plotMesh(mesh, style, az, el)

if nargin < 2
    style = 'solid';
end

if strcmpi(style, 'mesh')
    h = trimesh(mesh.F', mesh.V(1,:)', mesh.V(2,:)', mesh.V(3,:)', 'FaceColor', 'none', 'EdgeColor', 'w', ...
        'AmbientStrength', 0.4, 'EdgeLighting', 'flat');
    set(gcf, 'Color', 'k', 'Renderer', 'OpenGL');
    set(gca, 'Projection', 'perspective');
    axis equal;
    axis off;
    view(az,el);
elseif strcmpi(style, 'solid')
    h = trimesh(mesh.F', mesh.V(1,:)', mesh.V(2,:)' ,mesh.V(3,:)', 'FaceColor', 'w', 'EdgeColor', 'none', ...
        'AmbientStrength', 0.3, 'DiffuseStrength', 0.6, 'SpecularStrength', 0.0, 'FaceLighting', 'flat');
    set(gcf, 'Color', 'w', 'Renderer', 'OpenGL');
    set(gca, 'Projection', 'perspective');
    axis equal;
    axis off;
    view(az,el);
    camlight('HEADLIGHT');
elseif strcmpi(style, 'solidphong')
    mesh = normals(mesh);
    h = trimesh(mesh.F', mesh.V(1,:)', mesh.V(2,:)' ,mesh.V(3,:)', 'FaceColor', 'w', 'EdgeColor', 'none', ...
        'AmbientStrength', 0.3, 'DiffuseStrength', 0.6, 'SpecularStrength', 0.0, 'FaceLighting', 'gouraud', ...
        'VertexNormals', -mesh.Nv(1:3,:)', 'BackFaceLighting', 'reverselit');
    set(gcf, 'Color', 'w', 'Renderer', 'OpenGL');
    set(gca, 'Projection', 'perspective');
    axis equal;
    axis off;
    view(az,el);
    camlight('HEADLIGHT');
elseif strcmpi(style, 'soliddoublesided')
    mesh = normals(mesh);
    lx = cos(az) * cos(el);
    ly = cos(az) * sin(el);
    lz = sin(az);
    lightdir = [lx ly lz]';
    mesh.C = zeros( size(mesh.V, 2), 3 );
    for i=1:size(mesh.V, 2)
        mesh.C(i, :) = .3 + .6 * max( sum( lightdir .* mesh.Nv(:, i) ), sum( -lightdir .* mesh.Nv(:, i) ) );
    end
    h = trimesh(mesh.F', mesh.V(1,:)', mesh.V(2,:)' ,mesh.V(3,:)', 'EdgeColor', 'none', ...
        'FaceVertexCData', mesh.C, 'FaceColor', 'interp' );
    set(gcf, 'Color', 'w', 'Renderer', 'OpenGL');
    set(gca, 'Projection', 'perspective');
    axis equal;
    axis off;
    view(az,el);
    camlight;
end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


function varargout = test_popup_menu(varargin)
% TEST_POPUP_MENU M-file for test_popup_menu.fig
%      TEST_POPUP_MENU, by itself, creates a new TEST_POPUP_MENU or raises the existing
%      singleton*.
%
%      H = TEST_POPUP_MENU returns the handle to a new TEST_POPUP_MENU or the handle to
%      the existing singleton*.
%
%      TEST_POPUP_MENU('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in TEST_POPUP_MENU.M with the given input arguments.
%
%      TEST_POPUP_MENU('Property','Value',...) creates a new TEST_POPUP_MENU or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before test_popup_menu_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to test_popup_menu_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 13-Mar-2004 16:54:29

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

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

% Find the value of the popup menu
value = get(handles.Popup1,'Value');

% Place the value into the text field
str = ['Option ' num2str(value)];
set (handles.Label1,'String',str);

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

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


text:
function y = SSST (x, p)
% SSST  Spherical Surface Spline in Tension
%
%   y = SSST (x, p)
%
%   x is cos(theta) in range -1 <= x <= 1
%   p is tension (p >= 0)
%
% Returns the Green's' function for a spherical surface spline
% in tension, following Wessel and Becker [2008].
% If p == 0 or not given then we return Parker's [1994] minimum
% curvature solution instead.

% $Id: SSST.m,v 1.1.1.1 2008/05/09 21:34:52 myself Exp $
% P. Wessel, SOEST, U of Hawaii, April 2008 (pwessel@hawaii.edu)

if nargin == 1
	p = 0;
end
if (p == 0)     % Just do Parker's dilog solution
    y = dilog (0.5 - 0.5 * x);
    return
end

% Here we have tension
v = tension2nu (p);

y = zeros(size(x));     % Initialize output array
A = pi / sin(v*pi);     % Constant scale
k = find (abs(x) < (1-eps));  % Where Pv solution works
if (~isempty(k))
    y(k) = A*Pv(-x(k),v) - log(1-x(k));
end
% Deal with special case x == -1
k = find ((x-eps) <= -1);
if (~isempty(k))
    y(k) = A - log(2);
end
% Deal with special case x == +1
k = find ((x+eps) >= +1);
if (~isempty(k))
    y(k) = pi*cot(v*pi) + 2*(0.577215664901 + psi(1+v)) - log(2);
end
y = real(y);    % Knock off any insignificant imaginary noise

% Sub-functions used by SSST
function P = Pv (x, v)

P = zeros(size(x));
for i = 1:length(x)
    if (x(i) == -1)
        p = inf;
    else
        [p q k] = PvQv (x(i), v);
    end
    P(i) = p;
end

function [Pv Qv iter] = PvQv (x, v)
% Based on recipe in "An Atlas of Functions" by
% Spanier and Oldham, 1987
iter = 0;
if (x == -1)
    Pv = -inf;
    Qv = -inf;
    return
end
if (x == +1)
    Pv = 1;
    Qv = inf;
    return
end
a = 1;
R = 1;
K = 4 * sqrt (abs(v - v^2));
if (abs(1+v) + floor (1+v)) == 0
	a = 1.0e99;
	v = -1 - v;
end
s = sin (0.5*pi*v);
c = cos (0.5*pi*v);
w = (0.5 + v)^2;
while v <= 6.0
	v = v + 2;
	R = R * (v - 1)/v;
end
X = 1.0 / (4 + 4*v);
g = 1 + 5*X*(1 - 3*X*(0.35+6.1*X));
R = R*(1 - X*(1 - g*X/2))/sqrt (8*X);
g = 2*x;
u = g;
f = 1;
t = 1;
k = 0.5;
X = 1 + (1e8/(1 - x.^2));

t = t .* x.^2 * (k^2 - w) / ((k + 1)^2 - 0.25);
k = k + 1;
f = f + t;
u = u .* x.^2 * (k^2 - w) / ((k + 1)^2 - 0.25);
k = k + 1;
g = g + u;
while (k < K || abs (X*t) > abs(f))
        iter = iter + 1;
	t = t .* x.^2 * (k^2 - w) / ((k + 1)^2 - 0.25);
	k = k + 1;
	f = f + t;
	u = u .* x.^2 * (k^2 - w) / ((k + 1)^2 - 0.25);
	k = k + 1;
	g = g + u;
end
f = f + (x.^2.*t ./ (1 - x.^2));
g = g + (x.^2.*u ./ (1 - x.^2));
Pv = ((s*g*R) + (c*f/R))/sqrt(pi);
Qv = a*sqrt(pi)*((c*g*R) - (s*f/R))/2;

function [f] = psi(z)
%Psi     Psi (or Digamma) function valid in the entire complex plane.
%
%                 d
%        Psi(z) = --log(Gamma(z))
%                 dz
%
%usage: [f] = psi(z)
%
%        Z may be complex and of any size.
%
%        This program uses the analytical derivative of the
%        Log of an excellent Lanczos series approximation
%        for the Gamma function.
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%
siz = size(z);
z=z(:);
zz=z;

f = zeros(size(z));

%reflection point
p=find(real(z)<0.5);
if ~isempty(p)
   z(p)=1-z(p);
end

%Lanczos approximation for the complex plane

g=607/128; % best results when 4<=g<=5

c = [  0.99999999999999709182;
      57.156235665862923517;
     -59.597960355475491248;
      14.136097974741747174;
      -0.49191381609762019978;
        .33994649984811888699e-4;
        .46523628927048575665e-4;
       -.98374475304879564677e-4;
        .15808870322491248884e-3;
       -.21026444172410488319e-3;
        .21743961811521264320e-3;
       -.16431810653676389022e-3;
        .84418223983852743293e-4;
       -.26190838401581408670e-4;
        .36899182659531622704e-5];

n=0;
d=0;
for k=size(c,1):-1:2
    dz=1./(z+k-2);
    dd=c(k).*dz;
    d=d+dd;
    n=n-dd.*dz;
end
d=d+c(1);
gg=z+g-0.5;
%log is accurate to about 13 digits...

f = log(gg) + (n./d - g./gg) ;

if ~isempty(p)
   f(p) = f(p)-pi*cot(pi*zz(p));
end

p=find(round(zz)==zz & real(zz)<=0 & imag(zz)==0);
if ~isempty(p)
   f(p) = Inf;
end

f=reshape(f,siz);

function y = dilog (x)
% DILOG   The dilogarithm
%
%   y = dilog (x)
%
% Compute dilog(x) (defined for x >= 0) by the method of Parker's
% Appendix A of his Geophysical Inverse Theory.  The function
% is needed for x in the range 0 <= x <= 1 when solving the
% spherical spline interpolation in section 2.07 of Parker.

y = zeros (size (x));

pisqon6 = pi * pi / 6.0;
k = find (x <= 0.0);
if (~isempty(k))
    y(k) = pisqon6;
end
k = find (x > 0.0 & x < 0.5);
if (~isempty(k))
	y(k) = -log (1.0 - x(k));
	ysq = y(k) .* y(k);
	z = y(k) .* (1.0 + y(k) .* (-0.25 + y(k) .* (0.027777777777213 + ...
        ysq .* (-2.7777776990e-04 + ysq .* (4.724071696e-06 + ...
        ysq .* (-9.1764954e-08 + 1.798670e-09 .* ysq))))));
	y(k) = pisqon6 - z + y(k) .* log (x(k));
end
k = find (x >= 0.5 & x < 2.0);
if (~isempty(k))
    y(k) = -log (x(k));
	ysq = y(k) .* y(k);
	z = y(k) .* (1.0 + y(k) .* (-0.25 + y(k) .* (0.027777777777213 + ...
        ysq .* (-2.7777776990e-04 + ysq .* (4.724071696e-06 + ...
        ysq .* (-9.1764954e-08 + 1.798670e-09 .* ysq))))));
    y(k) = z;
end
k = find (x >= 2.0);
if (~isempty(k))
    y(k) = log (x(k));
	ysq = y(k) .* y(k);
	z = y(k) .* (1.0 + y(k) .* (-0.25 + y(k) .* (0.027777777777213 + ...
        ysq .* (-2.7777776990e-04 + ysq .* (4.724071696e-06 + ...
        ysq .* (-9.1764954e-08 + 1.798670e-09 .* ysq))))));
        y(k) = -z - 0.5 * ysq;
end

function nu = tension2nu (p)
nu = (-1 + sqrt (1 - 4*p^2))/2;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [] = ReqName_callback(block)

% Callback processing function for mask parameter
% - number of Preconditions used in the Requirment Name Set
set_param(block, 'LinkStatus', 'none');
load_system('Simulink'); % Load Simulink library (background)
load_system('MIL_Test');

No = str2double(get_param(block,'No')); % Get # of Reqs
OldNo = length(find_system(block, 'regexp', 'on','LookUnderMasks','all','Name','<Precondition name>'));

% Add or remove model items
for i = 1:max(No,OldNo)
    if OldNo > No && i > No && i <= OldNo % Remove
        % Remove lines
        set_param([block '/switch <signal1>'], 'Inputs', num2str(No+1));
        PreLines = get_param([block '/<Precondition name>' num2str(i)],'LineHandles');

        delete_line(PreLines.Inport(1));
        delete_line(PreLines.Outport(2));

        % Remove blocks
        delete_block([block '/<Precondition name>' num2str(i)]);

    elseif No > OldNo && i > OldNo && i <= No % Add
        % Add blocks
        add_block('MIL_Test/Test Data/Test Data Architecture/<Test data generator>/<Requirement name>/<Precondition name>',[block '/<Precondition name>' num2str(i)]);
        set_param([block '/<Precondition name>' num2str(i)], 'LinkStatus', 'none');
        if i == 2
             pos = get_param([block '/<Precondition name>'], 'Position');
        else pos = get_param([block '/<Precondition name>' num2str(i-1)], 'Position');
        end
        pos1 = pos(1);
        pos3 = pos(3);
        pos2 = pos(2)+100;
        pos4 = pos(4)+100;
        set_param([block '/<Precondition name>' num2str(i)], 'Position', [pos1 pos2 pos3 pos4]);
        set_param([block '/switch <signal1>'], 'Inputs', num2str(i+1));
        % Add lines

        add_line(block,['Out Bus/1'],['<Precondition name>' num2str(i) '/1'],'autorouting','on');
        add_line(block,['<Precondition name>' num2str(i) '/2'],['switch <signal1>/' num2str(i+2)],'autorouting','on');

    end

end
if No == 1
        InPos = get_param([block '/<Precondition name>'], 'Position');
else
        InPos = get_param([block '/<Precondition name>' num2str(No)], 'Position');
end
       InP1 = InPos(1)+130;
       InP3 = InPos(3)+120;
       InP2 = InPos(2)+80;
       InP4 = InPos(4)+70;
       set_param([block '/Initialisation & Stabilisation'], 'Position',[InP1 InP2 InP3 InP4]);
parentOfReq = get_param(block, 'Parent');
system_name1(parentOfReq);

author: Sable

repo: mcbench-benchmarks


text:
%%train network config
param.caffe_path=fullfile(fileparts(pwd),'matcaffe');
param.solver_netfile=fullfile(fileparts(pwd),'model','solver_bn.prototxt');
param.fintune_model=fullfile(fileparts(pwd),'pretrain_model','googlenet_bn.caffemodel');
param.test_net_file=fullfile(fileparts(pwd),'model','test_bn.prototxt');
param.test_batch_size=32;
param.result_save_file='LPW.txt';
param.save_model_file='LPW';
param.save_model_name='LPW_iter';
param.test_interval=500;
param.train_maxiter=50000;
param.data_path=fullfile(fileparts(pwd),'generate_data');
param.traindata_filename='train_data_LPW';
param.testdata_filename='test_data_LPW';
param.use_data_split_index=1;
param.data_split_num=1;
param.test_person_num=1975;
param.use_gpu=1;
param.gpu_id=0;
%%
for split_index=param.use_data_split_index:param.data_split_num
    if ~exist(strcat(param.save_model_file,num2str(split_index)),'file')
        mkdir(strcat(param.save_model_file,num2str(split_index)));
    end
    load(strcat(param.data_path,param.traindata_filename,num2str(split_index),'/train_data.mat'));
    %% find caffe
    cur_path=pwd;
    addpath(genpath(param.caffe_path))
    cd(param.caffe_path);
    caffe.reset_all;
    caffe.init_log(fullfile(cur_path,'log'));
    addpath(param.caffe_path);
    if param.use_gpu
        caffe.set_mode_gpu;
    else
        caffe.set_mode_cpu;
    end
    cd(cur_path);
    caffe_solver=caffe.get_solver(param.solver_netfile,param.gpu_id);
	net=caffe.get_net(param.test_net_file,'test');
    caffe_solver.use_caffemodel(param.fintune_model);
    input_data_shape=caffe_solver.nets{1}.blobs('data').shape;
    batch_size=input_data_shape(4);
    batch_data=zeros(input_data_shape,'single');
    batch_label=zeros(1,1,1,batch_size,'single');
    train_x_axis=[];
    train_y_axis=[];
    label_train=[label_train_cam2_view1;label_train_cam2_view2;label_train_cam2_view3;label_train_cam2_view4; ...
        label_train_cam3_view1;label_train_cam3_view2;label_train_cam3_view3;label_train_cam3_view4];
    train_image_name=[train_image_name_cam2_view1,train_image_name_cam2_view2,train_image_name_cam2_view3,train_image_name_cam2_view4, ...
        train_image_name_cam3_view1,train_image_name_cam3_view2,train_image_name_cam3_view3,train_image_name_cam3_view4];
    iter=0;
    while iter<param.train_maxiter
        [batch_data,batch_label]=get_train_minibatch_for_LPW(batch_data,batch_label,batch_size, ...
            label_train,train_image_name);
        caffe_solver.nets{1}.blobs('data').set_data(batch_data);
        caffe_solver.nets{1}.blobs('label').set_data(batch_label);
        caffe_solver.step(1);
        iter = caffe_solver.iter;
        if mod(iter,10)==0
            soft_loss=caffe_solver.nets{1}.blobs('loss3/loss').get_data;
            acc=caffe_solver.nets{1}.blobs('loss3/top-1').get_data;
            loss1=caffe_solver.nets{1}.blobs('loss1/loss').get_data;
            loss2=caffe_solver.nets{1}.blobs('loss2/loss').get_data;
            train_x_axis=[train_x_axis,iter];
            train_y_axis=[train_y_axis,soft_loss];
            plot(train_x_axis,train_y_axis);
            drawnow;
            fprintf('iter= %d ,softmaxloss= %f ,loss1=%f  loss2=%f\n',iter,soft_loss,loss1,loss2);
        end
        if mod(iter,param.test_interval)==0
            caffe_solver.nets{1}.save(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
            net.copy_from(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
            net.blobs('data').reshape([224 224 3 param.test_batch_size]);
            load(strcat(param.data_path,param.testdata_filename,num2str(split_index),'/test_data.mat'));
            test_script;
        end
        if iter>10000 && mod(iter,200)==0
            caffe_solver.nets{1}.save(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
        end
    end
    %test_network_for_LPW;
end

author: liuyuisanai

repo: Quality-Aware-Network


text:
% EEG_REJMACRO - Internal EEGLAB macro for all pop_ functions that
%                  perform data rejection.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGLAB

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

% 01-25-02 reformated help & license -ad
% 03-08-02 include past rejections in eegplot -ad

% this macro is used by pop_eegplot, pop_eegthresh, pop_const, pop_jointprob, pop_eegthreshfreq, pop_kurtosis

% inputs:
%   rej - input rejection in this workspace (string)
%   rejE - electrode rejection in this workspace (string)
%   elecrange - range of selected electrodes
%   reject - superpose with previous rejections
% optional : nbpnts = number of points per trial

% outputs:
%   rejection arrays ('oldrej', 'oldrejE'), 'command' for eeglot

if ~exist('elecrange')
	help eeg_rejmacro;
	error('Error: eeg_rejmacro cannot be called from the command line');
end;

if ~exist('nbpnts') nbpnts = EEG.pnts; end

% mix all type of rejections
% --------------------------
if icacomp
    nChan = EEG.nbchan;
else
    nChan = size(EEG.icaweights,1);
end
if superpose == 2
    com2 = ['if ~isempty(TMPREJ), ' ...
            '  icaprefix = ' fastif(icacomp, '''''', '''ica''') ';' ...
            '  for indextmp = 1:length(EEG.reject.disprej),' ...
            '     eval([ ''colortmp = EEG.reject.rej'' EEG.reject.disprej{indextmp} ''col;'' ]);' ...
            '     [tmprej tmprejE] = eegplot2trial(TMPREJ,' int2str(nbpnts) ', EEG.trials, colortmp, []);' ...
            '     if ~isempty(tmprejE),' ...
            '          tmprejE2 = zeros(' int2str(nChan) ', length(tmprej));' ...
            '          tmprejE2([' int2str(elecrange) '],:) = tmprejE;' ...
            '     else,' ...
            '          tmprejE2 = [];' ...
            '     end;' ...
            '     eval([ ''EEG.reject.'' icaprefix ''rej'' EEG.reject.disprej{indextmp} ''= tmprej;'' ]);' ...
            '     eval([ ''EEG.reject.'' icaprefix ''rej'' EEG.reject.disprej{indextmp} ''E = tmprejE2;'' ]);' ...
            '  end;' ];
%           '     colortmp, tmpstr = [ ''EEG.reject.'' icaprefix ''rej'' EEG.reject.disprej{indextmp} ], sum(tmprej),' ...
else
    com2 = [ 'if ~isempty(TMPREJ),' ...
             '  icaprefix = ' fastif(icacomp, '''''', '''ica''') ';' ...
             '  [tmprej tmprejE] = eegplot2trial(TMPREJ,' int2str(nbpnts) ', EEG.trials, [' num2str(colrej) '], []);' ...
             '  if ~isempty(tmprejE),' ...
             '     tmprejE2 = zeros(' int2str(nChan) ', length(tmprej));' ...
             '     tmprejE2([' int2str(elecrange) '],:) = tmprejE;' ...
             '  else,' ...
             '     tmprejE2 = [];' ...
             '  end;' ...
             macrorej '= tmprej;' macrorejE '= tmprejE2;' ...
             ... % below are manual rejections
             '  tmpstr = [ ''EEG.reject.'' icaprefix ''rejmanual'' ];' ...
             '  if ~isempty(tmprej) eval([ ''if ~isempty('' tmpstr ''),'' tmpstr ''='' tmpstr ''| tmprej; else, '' tmpstr ''=tmprej; end;'' ]); end;' ...
             '  if ~isempty(tmprejE2) eval([ ''if ~isempty('' tmpstr ''E),'' tmpstr ''E='' tmpstr ''E| tmprejE2; else, '' tmpstr ''E=tmprejE2; end;'' ]); end;' ];
%             '  size(tmprejE2), eval([''disp(size('' tmpstr ''E))'']),' ...
end
% text commented below to fix BUG 478
%             '  [tmprej tmprejE] = eegplot2trial(TMPREJ,' int2str(nbpnts) ', EEG.trials, [' num2str(EEG.reject.rejmanualcol) '], []);' ...
%             '  if ~isempty(tmprejE),' ...
%             '     tmprejE2 = zeros(EEG.nbchan, length(tmprej));' ...
%             '     tmprejE2([' int2str(elecrange) '],:) = tmprejE;' ...
%             '  else,' ...
%             '     tmprejE2 = [];' ...
%             '  end;' ...
if reject
    com2 = [ com2 ...
           '[EEGTMP LASTCOM] = pop_rejepoch(EEG, tmprej, 1);' ...
		   'if ~isempty(LASTCOM),'...
           '   EEG = eegh(LASTCOM, EEG);' ...
           '   [ALLEEG EEG CURRENTSET tmpcom] = pop_newset(ALLEEG, EEGTMP, CURRENTSET);' ...
           '   if ~isempty(tmpcom),' ...
           '     eegh(tmpcom);' ...
           '     eeglab(''redraw'');' ...
            '  end;' ...
	       'end; clear EEGTMP tmpcom; end;' ] ;
else
	com2 = [com2 ...
			' warndlg2(strvcat(''Epochs (=trials) marked for rejection have been noted.'',' ...
			'''To actually reject these epochs, use the same menu item to'',' ...
            '''inspect/reject data and select "Reject marked trials" checkbox''), ''Warning'');' ...
			'[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); eeglab(''redraw''); end;' ];
end;
if ~exist('topcommand')
	topcommand = [];
end

% the first part is used to convert the eegplot output
command = [  com2 topcommand 'clear indextmp colortmp icaprefix tmpcom tmprej tmprejE tmprejE2 TMPREJ;' ];

if all(colrej == EEG.reject.rejmanualcol)
	oldrej = [];  % for manual rejection, old rejection are
	oldrejE = []; % the current rejection
else
	oldrej  = eval(macrorej);
	oldrejE = eval(macrorejE);
end

switch superpose
 case 0, rejeegplot = trial2eegplot(  rej, rejE, nbpnts, colrej);
 case 1, rejeegplottmp = trial2eegplot(  oldrej, oldrejE, nbpnts, min(colrej+0.15, [1 1 1]));
         if ~isempty(rejeegplottmp), rejeegplot = [ rejeegplottmp ];
		 else rejeegplot = []; end
         rejeegplottmp = trial2eegplot(  rej, rejE, nbpnts, colrej);
         if ~isempty(rejeegplottmp), rejeegplot = [ rejeegplot; rejeegplottmp ]; end
 case 2,
  rejeegplot = [];
  for index = 1:length(EEG.reject.disprej)
	  if ~isempty(EEG.reject.disprej{index})
		  eval([ 'colortmp = EEG.reject.rej' EEG.reject.disprej{index} 'col;']);
		  if any(colortmp ~= colrej) % test if current rejection (if color different)
			  if icacomp == 0 % ica
				  currentname = [ 'EEG.reject.icarej' EEG.reject.disprej{index} ];
			  else
				  currentname = [ 'EEG.reject.rej' EEG.reject.disprej{index} ];
			  end
			  currentcolor =  [ 'EEG.reject.rej' EEG.reject.disprej{index} 'col' ];
			  %if strcmp(EEG.reject.disprej{index}, 'manual')
			  %	  currentcolor = [ 'min(' currentcolor '+0.15, [1 1 1])' ];
			  %end; % using this test, manual rejections won't be added to current rej
			  eval( [ 'rejeegplottmp = trial2eegplot( ' currentname ',' currentname ...
					  'E, nbpnts,' currentcolor ');' ]);
			  if ~isempty(rejeegplottmp), rejeegplot = [ rejeegplot; rejeegplottmp ]; end
		  end
	  end
  end
  rejeegplottmp = trial2eegplot(  rej, rejE, nbpnts, colrej);
  if ~isempty(rejeegplottmp), rejeegplot = [ rejeegplot; rejeegplottmp ]; end
end
if ~isempty(rejeegplot)
	rejeegplot = rejeegplot(:,[1:5,elecrange+5]);
else
	rejeegplot = [];
end
eegplotoptions = { 'events', EEG.event, 'winlength', 5, 'winrej', ...
				   rejeegplot, 'xgrid', 'off', 'wincolor', EEG.reject.rejmanualcol, ...
				   'colmodif', { { EEG.reject.rejmanualcol EEG.reject.rejthreshcol EEG.reject.rejconstcol ...
                                   EEG.reject.rejjpcol     EEG.reject.rejkurtcol   EEG.reject.rejfreqcol } } };

if ~isempty(EEG.chanlocs) && icacomp == 1
    if exist('elecrange')
        eegplotoptions = { eegplotoptions{:}  'eloc_file', EEG.chanlocs(elecrange) };
    else
        eegplotoptions = { eegplotoptions{:}  'eloc_file', EEG.chanlocs };
    end
else
    if exist('elecrange')
        for index = 1:length(elecrange)
            tmpstruct(index).labels = int2str(elecrange(index));
        end
        eegplotoptions = { eegplotoptions{:}  'eloc_file' tmpstruct };
    end
end
if ~reject
	eegplotoptions = { eegplotoptions{:}  'butlabel', 'UPDATE MARKS' };
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function res = qdist(a,b)
%QDIST        Implements  q(a,b)  metrical distance
%  Name  qdist  is used to avoid ambiguities with variable  q
%  This functions for non-interval input only for completeness
%
%     res = qdist(a,b)
%
% for real input          abs(a-b)
% for complex input       qdist(real(a),real(b)) + qdist(imag(a),imag(b))
%

% written  03/23/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if isreal(a) & isreal(b)
    res = abs(a-b);
  else
    res = abs(real(a)-real(b)) + abs(imag(a)-imag(b));
  end

repo: OpenMAS


text:
function spline_test08 ( )

%*****************************************************************************80
%
%% TEST08 tests BASIS_MATRIX_OVERHAUSER_NONUNI, BASIS_MATRIX_TMP.
%
%  Discussion:
%
%    YDATA(1:NDATA) = ( TDATA(1:NDATA) - 2 )**2 + 3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 4;
  ndata = 4;
  nsample = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  BASIS_MATRIX_OVERHAUSER_NONUNI sets up the\n' );
  fprintf ( 1, '    basis matrix for the nonuniform Overhauser\n' );
  fprintf ( 1, '    spline.\n' );

  tdata = [ 0.0E+00, 1.0E+00, 2.0E+00, 3.0E+00 ];

  alpha = ( tdata(3) - tdata(2) ) / ( tdata(3) - tdata(1) );
  beta =  ( tdata(3) - tdata(2) ) / ( tdata(4) - tdata(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  ALPHA = %14f\n', alpha );
  fprintf ( 1, '  BETA  = %14f\n', beta  );

  mbasis = basis_matrix_overhauser_nonuni ( alpha, beta );

  for i = 1 : ndata
    ydata(i) = ( tdata(i) - 2.0E+00 )^2 + 3.0E+00;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%14f  %14f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %14f  %14f\n', mark, tval, yval );

    end

  end

  tdata(1:4) = [ 0.0E+00, 1.0E+00, 2.0E+00, 5.0E+00 ];

  alpha = ( tdata(3) - tdata(2) ) / ( tdata(3) - tdata(1) );
  beta =  ( tdata(3) - tdata(2) ) / ( tdata(4) - tdata(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  ALPHA = %14f\n', alpha );
  fprintf ( 1, '  BETA  = %14f\n', beta  );

  mbasis = basis_matrix_overhauser_nonuni ( alpha, beta );

  for i = 1 : ndata
    ydata(i) = ( tdata(i) - 2.0E+00 )^2 + 3.0E+00;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%14f  %14f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %14f  %14f\n', mark, tval, yval );

    end

  end

  tdata(1:4) = [ 0.0E+00, 3.0E+00, 4.0E+00, 5.0E+00 ];

  alpha = ( tdata(3) - tdata(2) ) / ( tdata(3) - tdata(1) );
  beta =  ( tdata(3) - tdata(2) ) / ( tdata(4) - tdata(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  ALPHA = %14f\n', alpha );
  fprintf ( 1, '  BETA  = %14f\n', beta  );

  mbasis = basis_matrix_overhauser_nonuni ( alpha, beta );

  for i = 1 : ndata
    ydata(i) = ( tdata(i) - 2.0E+00 )^2 + 3.0E+00;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '    TDATA, YDATA\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ndata
    fprintf ( 1, '%14f  %14f\n', tdata(i), ydata(i) );
  end

  left = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '    T, Spline(T)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ndata

    if ( i == 0 )
      tlo = tdata(1) - 0.5E+00 * ( tdata(2) - tdata(1) );
      thi = tdata(1);
    elseif ( i < ndata )
      tlo = tdata(i);
      thi = tdata(i+1);
    elseif ( ndata <= i )
      tlo = tdata(ndata);
      thi = tdata(ndata) + 0.5E+00 * ( tdata(ndata) - tdata(ndata-1) );
    end

    if ( i < ndata )
      jhi = nsample - 1;
    else
      jhi = nsample;
    end

    for j = 0 : jhi

      tval = ( ( nsample - j ) * tlo   ...
             + (           j ) * thi ) ...
             /   nsample;

      yval = basis_matrix_tmp ( left, n, mbasis, ndata, tdata, ydata, tval );

      if ( 0 < i & j == 0 )
        mark = '*';
      else
        mark = ' ';
      end

      fprintf ( 1, '%c  %14f  %14f\n', mark, tval, yval );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [varargout] = subsref(obj,s)

n = numel(s);
if n == 0,
  varargout{1} = obj;
  return;
end

% which flies are referenced?
i = 1;
switch s(i).type,
  case '()',

    switch numel(s(i).subs),

      % one argument means fly index
      case 1,
        if ischar(s(i).subs{1}) && strcmp(s(i).subs{1},':'),
          idx = 1:obj.nflies;
        elseif isnumeric(s(i).subs{1}),
          idx = s(i).subs{1};
          if any(idx < 1 | idx > obj.nflies | idx ~= round(idx)),
            error('Indices must be whole numbers between 1 and %d',obj.nflies);
          end
        else
          error('Bad indexing argument');
        end

      % two arguments means experiment + fly
      case 2,

        % experiment indices
        j = 1;
        if ischar(s(i).subs{j}) && strcmp(s(i).subs{j},':'),
          exps = 1:obj.nexpdirs;
        elseif ischar(s(i).subs{j}) && ismember(s(i).subs{j},obj.expdirs),
          exps = find(strcmp(s(i).subs{j},obj.expdirs),1);
        elseif iscell(s(i).subs{j}),
          [didfind,exps] = ismember(s(i).subs{j},obj.expdirs);
          if any(~didfind),
            error(['Bad experiment indexing arguments: ',sprintf('%s ',s(i).subs{j}{~didfind})]);
          end
        elseif isnumeric(s(i).subs{j}),
          exps = s(i).subs{j};
          if any(exps < 1 | exps > obj.nexpdirs | exps ~= round(exps)),
            error('Experiment indices must be whole numbers between 1 and nexpdirs = %d',obj.nexpdirs);
          end
        else
          error('Bad experiment indexing argument');
        end

        % fly indices
        j = 2;
        if ischar(s(i).subs{j}) && strcmp(s(i).subs{j},':'),
          idx = [obj.exp2flies{exps}];
        elseif isnumeric(s(i).subs{j}),
          idx = cell2mat(cellfun(@(x)x(s(i).subs{j}),obj.exp2flies(exps),'UniformOutput',false));
        else
          error('Bad fly indexing argument');
        end

      otherwise
        error('Number of indexing arguments must be <= 2');
    end

    % go on to next argument
    i = 2;

  case '.',

    if ismember(s(i).subs,methods(obj)) || ismember(s(i).subs,properties(obj)),
      [varargout{1:nargout}] = builtin('subsref',obj,s);
      return;
    end

    idx = 1:obj.nflies;

    % use this argument
    i = 1;

  otherwise
    error('{} referencing not supported');
end

if n < i,
  % TODO: maybe implement this?
  error('Field name must be specified');
end

% get field
if ~strcmpi(s(i).type,'.'),
  error('Field name must be specified');
end
fn = s(i).subs;
i = i + 1;

% get indices
elidx = [];
if n < i,
  t = ':';
else
  if ~strcmp(s(i).type,'()'),
    error('Frame indices must be indexed by ()');
  end
  if numel(s(i).subs) > 2,
    error('Can only index as (i,t) or (t)');
  elseif numel(s(i).subs) == 2,
    elidx = s(i).subs{1};
    t = s(i).subs{2};
  else
    t = s(i).subs{1};
  end
end

if strcmp(fn,'nframes'),
  res = {obj.nframes(idx)};
elseif strcmp(fn,'firstframe'),
  res = {obj.firstframes(idx)};
elseif strcmp(fn,'endframe'),
  res = {obj.endframes(idx)};
elseif strcmp(fn,'off'),
  res = {1-obj.firstframes(idx)};
elseif strcmp(fn,'trkInfo')
  res = {obj.trkInfo{idx}};
elseif (numel(fn)>=4) && strcmp(fn(1:4),'kpts')
  res = {obj.(fn){idx}};
else
  res = cell(1,numel(idx));
  for j = 1:numel(idx),
    x = obj.GetPerFrameData(fn,idx(j));
    if ~isempty(elidx),
      if ~ischar(elidx),
        x = x(elidx,:);
      end
      if ischar(t),
        res{j} = x;
      else
        res{j} = x(:,t);
      end
    else
      if ischar(t),
        res{j} = x;
      else
        res{j} = x(t);
      end
    end
  end
end
if nargout>0 && nargout < numel(res),
  if numel(res) == 1,
    varargout{1:nargout} = res{1};
  else
    varargout{1:nargout} = res;
  end
else
  varargout = res;
end

author: kristinbranson

repo: JAABA


text:
function name = dtiGetCurAnatName(handles)
% OBSOLETE
%
% name = dtiGetCurAnatName(handles)
%
% HISTORY:
% 2004.02.09 RFD (bob@white.stanford.edu) wrote it.
%

warning('Use dtiGet(handles,''currentAnatomyName'')');

opts = get(handles.popupBackground,'String');
name = opts{get(handles.popupBackground,'Value')};

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function finetuned_model_path = spp_finetune_voc(opts)
% [] = spp_finetune_voc(opts)
%   finetune fc layers
%
% this version read conv feature maps from disk each iteration for compatibility and small
% memory usage. load all conv feature maps into memory will accelerate finetuning
% greatly
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

% -------------------- CONFIG --------------------
if ~exist('opts', 'var')
    opts.net_file               = fullfile(pwd, 'data\cnn_model\Zeiler_conv5\Zeiler_conv5');
    opts.finetune_net_def_file  = fullfile(pwd, 'model-defs\pascal_finetune_fc_spm_solver.prototxt');
    opts.spp_params_def         = fullfile(pwd, 'data\cnn_model\Zeiler_conv5\spp_config');

    opts.flip = true;
    opts.finetune_cache_name = 'Zeiler_conv5_ft(5s_flip)_ss';

    opts.with_selective_search = true;
    opts.with_edge_box = false;

    opts.devkit = './datasets/VOCdevkit2007';
    opts.imdb_train   = {   imdb_from_voc(opts.devkit, 'trainval', '2007', opts.flip) };
    opts.feat_cache_train   = { 'Zeiler_conv5' };

    opts.imdb_test    = { imdb_from_voc(opts.devkit, 'test', '2007') };
    opts.feat_cache_test = { 'Zeiler_conv5' };

    opts.with_hard_samples = false;  % only for train
else
    if ~iscell(opts.imdb_test)
        opts.imdb_test = {opts.imdb_test};
        opts.feat_cache_test = {opts.feat_cache_test};
    end
end

opts.random_scale = 1;

data_param.img_num_per_iter = 128; % should be same with the prototxt file
data_param.random_scale = opts.random_scale;
data_param.iter_per_batch = 500; % for load data effectively
data_param.fg_fraction = 0.25;
data_param.fg_threshold = 0.5;
data_param.bg_threshold = [0.1 0.5];

solver_param.test_iter = 500; % better times of data_param.iter_per_batch
solver_param.test_interval = 2000; % better times of data_param.iter_per_batch
if isfield(opts, 'max_iter')
    solver_param.max_iter = opts.max_iter;
else
    solver_param.max_iter = 1000000;
end

opts.imdb_train = opts.imdb_train(:)';
opts.feat_cache_train = opts.feat_cache_train(:)';
opts.imdb_test = opts.imdb_test(:)';
opts.feat_cache_test = opts.feat_cache_test(:)';
roidb_train = cellfun(@(x) x.roidb_func(x, opts.with_hard_samples, opts.with_selective_search, opts.with_edge_box), opts.imdb_train, 'UniformOutput', false);
roidb_test = cellfun(@(x) x.roidb_func(x, true, opts.with_selective_search, opts.with_edge_box), opts.imdb_test, 'UniformOutput', false);
% ------------------------------------------------

work_dir = fullfile(pwd, 'finetuning', opts.finetune_cache_name);
timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
log_file = fullfile(work_dir, 'output', [opts.finetune_cache_name, '_' timestamp, '.txt']);
caffe_log_file = fullfile(work_dir, opts.finetune_cache_name);

mkdir_if_missing(work_dir);
mkdir_if_missing(fileparts(log_file));
diary(log_file);

fprintf('data_param:\n');
disp(data_param);
fprintf('solver_param:\n');
disp(solver_param);
fprintf('opts_param:\n');
disp(opts);

%% get all data info
[fg_windows_train, bg_windows_train] = cellfun(@(x) SetupData(x, data_param), roidb_train, 'UniformOutput', false);
[fg_windows_test, bg_windows_test] = cellfun(@(x) SetupData(x, data_param), roidb_test, 'UniformOutput', false);
fg_windows_num_train = cellfun(@(y) cellfun(@(x) size(x, 1), y, 'UniformOutput', true), ...
    fg_windows_train, 'UniformOutput', false);
fg_windows_num_test = cellfun(@(y) cellfun(@(x) size(x, 1), y, 'UniformOutput', true), ...
    fg_windows_test, 'UniformOutput', false);
bg_windows_num_train = cellfun(@(y) cellfun(@(x) size(x, 1), y, 'UniformOutput', true), ...
    bg_windows_train, 'UniformOutput', false);
bg_windows_num_test = cellfun(@(y) cellfun(@(x) size(x, 1), y, 'UniformOutput', true), ...
    bg_windows_test, 'UniformOutput', false);

fg_windows_train = CombineImdbWindows(fg_windows_train);
bg_windows_train = CombineImdbWindows(bg_windows_train);
fg_windows_test = CombineImdbWindows(fg_windows_test);
bg_windows_test = CombineImdbWindows(bg_windows_test);

fg_windows_num_train = cell2mat(fg_windows_num_train');
bg_windows_num_train = cell2mat(bg_windows_num_train');
fg_windows_num_test = cell2mat(fg_windows_num_test');
bg_windows_num_test = cell2mat(bg_windows_num_test');

fprintf('total fg_windows_num_train : %d\n', sum(fg_windows_num_train));
fprintf('total bg_windows_num_train : %d\n', sum(bg_windows_num_train));
fprintf('total fg_windows_num_test : %d\n', sum(fg_windows_num_test));
fprintf('total bg_windows_num_test : %d\n', sum(bg_windows_num_test));
%% init caffe solver
caffe('set_random_seed', 6);
spp_model.spp_pooler = spp_load_pooling_params(opts.spp_params_def);
InitSolver(opts.finetune_net_def_file, opts.net_file, caffe_log_file);

% finetune
caffe('set_mode_gpu');
caffe('set_phase_train');

%% get test sample
% fix the random seed for repeatability
prev_rng = seed_rand();
data_param_test = data_param;
data_param_test.iter_per_batch = solver_param.test_iter;
data_param_test.random_scale = 0; % for test use the best scale
[data_test, label_test] = PrepareBatchData(spp_model.spp_pooler, opts.imdb_test, fg_windows_test, ...
                    fg_windows_num_test, bg_windows_test, bg_windows_num_test, data_param_test, data_param.img_num_per_iter, opts.feat_cache_test);
rng(prev_rng);

%% finetune
prev_rng = seed_rand();
data_train = {}; label_train = {};
train_error_sum = 0;
train_error_count = 0;
min_test_error = 1;
iter_ = 1;
while(iter_ < solver_param.max_iter)
    if isempty(data_train)
        [data_train, label_train] = ...
        PrepareBatchData(spp_model.spp_pooler, opts.imdb_train, fg_windows_train, ...
        fg_windows_num_train, bg_windows_train, bg_windows_num_train, data_param, data_param.img_num_per_iter, opts.feat_cache_train);
    end
    rst = caffe('train', {data_train{1}, label_train{1}});
    data_train(1) = []; label_train(1) = [];
    train_error_sum = train_error_sum + rst(1);
    train_error_count = train_error_count + 1;

    % test per test_interval iterations
    if ~mod(iter_, solver_param.test_interval)
        test_error = nan(solver_param.test_iter, 1);
        for it = 1:solver_param.test_iter
            rst = caffe('test', {data_test{it}, label_test{it}});
            test_error(it) = rst(1);
        end
        ShowState(iter_, train_error_sum / train_error_count, mean(test_error));
        if (mean(test_error) < min_test_error)
            min_test_error = mean(test_error);
            finetuned_model_path = fullfile(opts.finetune_rst_dir, sprintf('FT_iter_%d', iter_));
            caffe('snapshot', finetuned_model_path);
            fprintf('best! save as %s\n', finetuned_model_path);
        end
        train_error_sum = 0;
        train_error_count = 0;

        diary; diary; % flush diary
    end

    iter_ = caffe('get_solver_iter');
end

rng(prev_rng);
caffe('release_solver');

end

function [data, label] = PrepareBatchData(spp_pooler, imdb, fg_windows, fg_windows_num, bg_windows, bg_windows_num, data_param, img_num_per_iter, feat_cache)
    % windows -- (image_id_in_total, image_imdb_idx, image_id_in_imdb,
    % label, ov, l, t, r, b)

    nTimesMoreData = 10;

    fg_num_each = int32(data_param.fg_fraction * img_num_per_iter);
    bg_num_each = img_num_per_iter - fg_num_each;
    fg_num_total = fg_num_each * data_param.iter_per_batch;
    bg_num_total = bg_num_each * data_param.iter_per_batch;
    total_img = size(fg_windows, 1);

    % random perm image and get top n image with total more than patch_num
    % windows
    img_idx = randperm(total_img);
    fg_windows_num = fg_windows_num(img_idx);
    fg_windows_cumsum = cumsum(fg_windows_num);
    bg_windows_num = bg_windows_num(img_idx);
    bg_windows_cumsum = cumsum(bg_windows_num);
    img_idx_end = max(find(fg_windows_cumsum > fg_num_total * nTimesMoreData, 1, 'first'), ...
        find(bg_windows_cumsum > bg_num_total * nTimesMoreData, 1, 'first'));
    if isempty( img_idx_end )
        img_idx_end = length(img_idx);
    end
    select_img_idx = img_idx(1:img_idx_end);

    fg_window = cell2mat(fg_windows(select_img_idx));
    bg_window = cell2mat(bg_windows(select_img_idx));

    % random perm all windows, and drop redundant data
    window_idx = randperm(size(fg_window, 1));
    fg_window = fg_window(window_idx(1:fg_num_total), :);
    window_idx = randperm(size(bg_window, 1));
    bg_window = bg_window(window_idx(1:bg_num_total), :);

    fg_data = cell(1, length(select_img_idx)); bg_data = cell(1, length(select_img_idx));
    fg_label = cell(length(select_img_idx), 1); bg_label = cell(length(select_img_idx), 1);
    random_scale = data_param.random_scale;
    fg_window_this_img = cell(length(select_img_idx), 1);
    bg_window_this_img = cell(length(select_img_idx), 1);
    for i = 1:length(select_img_idx)
        fg_window_this_img{i} = fg_window(fg_window(:, 1) == select_img_idx(i), :);
        bg_window_this_img{i} = bg_window(bg_window(:, 1) == select_img_idx(i), :);
    end
    assert(sum(cellfun(@(x) size(x, 1), fg_window_this_img, 'UniformOutput', true)) == fg_num_total);
    assert(sum(cellfun(@(x) size(x, 1), bg_window_this_img, 'UniformOutput', true)) == bg_num_total);
    random_seed = randi(10^6, length(select_img_idx));
    parfor i = 1:length(select_img_idx)
        rng(random_seed(i), 'twister')
        window_this_img = [fg_window_this_img{i}; bg_window_this_img{i}];
        if isempty(window_this_img)
            continue;
        end
        image_id_in_total = select_img_idx(i);
        assert(all(window_this_img(:, 1) == image_id_in_total));
        image_imdb_idx = window_this_img(1, 2);
        assert(all(window_this_img(:, 2) == image_imdb_idx));
        image_id_in_imdb = window_this_img(1, 3);
        assert(all(window_this_img(:, 3) == image_id_in_imdb));

        [fg_data{i}, bg_data{i}] = LoadCachedPoolXFeature(spp_pooler, imdb{image_imdb_idx}, feat_cache{image_imdb_idx}, ...
            image_id_in_imdb, fg_window_this_img{i}(:, 6:9), bg_window_this_img{i}(:, 6:9), random_scale);
        fg_label{i} = fg_window_this_img{i}(:, 4);
        bg_label{i} = bg_window_this_img{i}(:, 4);
    end
    valid_fg_data = cellfun(@(x) ~isempty(x), fg_label, 'UniformOutput', true);
    valid_bg_data = cellfun(@(x) ~isempty(x), bg_label, 'UniformOutput', true);
    fg_data = cell2mat(fg_data(valid_fg_data));
    bg_data = cell2mat(bg_data(valid_bg_data));
    fg_label = cell2mat(fg_label(valid_fg_data));
    bg_label = cell2mat(bg_label(valid_bg_data));

    % random perm all windows
    fg_rnd_idx = randperm(size(fg_label, 1));
    bg_rnd_idx = randperm(size(bg_label, 1));
    fg_data = fg_data(:, fg_rnd_idx);
    fg_label = fg_label(fg_rnd_idx, :);
    bg_data = bg_data(:, bg_rnd_idx);
    bg_label = bg_label(bg_rnd_idx, :);

    % split to iter_per_batch
    fg_split_div = ones(1, data_param.iter_per_batch) * double(fg_num_each);
    bg_split_div = ones(1, data_param.iter_per_batch) * double(bg_num_each);
    fg_data = mat2cell(fg_data, size(fg_data, 1), fg_split_div);
    fg_label = mat2cell(fg_label', size(fg_label, 2), fg_split_div);
    bg_data = mat2cell(bg_data, size(bg_data, 1), bg_split_div);
    bg_label = mat2cell(bg_label', size(bg_label, 2), bg_split_div);

    % merge fg and bg_data
    data = cellfun(@(x, y) [x, y], fg_data, bg_data, 'UniformOutput', false);
    label = cellfun(@(x, y) [x, y], fg_label, bg_label, 'UniformOutput', false);
end

function [fg_windows, bg_windows] = SetupData(roidb, data_param)
    fg_windows = cell(length(roidb.rois), 1);
    bg_windows = cell(length(roidb.rois), 1);
    for i = 1:length(roidb.rois)
        roi = roidb.rois(i);
        [ov, label] = max(roi.overlap, [], 2);
        ov = full(ov);
        fg_mask = ov >= data_param.fg_threshold;
        fg_ov = ov(fg_mask);
        fg_label = label(fg_mask);
        fg_img_idx = fg_label * 0 + i;
        fg_window = roi.boxes(fg_mask, :);
        fg_windows{i} = [fg_img_idx, fg_label, fg_ov, fg_window];

        bg_mask = ~fg_mask & ov >= data_param.bg_threshold(1) & ov < data_param.bg_threshold(2);
        bg_ov = ov(bg_mask) * 0;
        bg_label = label(bg_mask) * 0;
        bg_img_idx = bg_label * 0 + i;
        bg_window = roi.boxes(bg_mask, :);
        bg_windows{i} = [bg_img_idx, bg_label, bg_ov, bg_window];
    end
end

function [fg_feat, bg_feat] = LoadCachedPoolXFeature(spp_pooler, imdb, feat_cache, img_id, fg_windows, bg_windows, random_scale)
    feat = spp_load_cached_poolX_features(spp_pooler, feat_cache, imdb.name, imdb.image_ids{img_id}, [fg_windows; bg_windows], random_scale);
    fg_feat = feat(:, 1:size(fg_windows, 1));
    bg_feat = feat(:, (size(fg_windows, 1)+1):end);
end

function ShowState(iter, train_error, test_error)
    fprintf('\n------------------------- Iteration %d -------------------------\n', iter);
    fprintf('Error - Training : %.4f - Testing : %.4f\n', train_error, test_error);
end

function windows = CombineImdbWindows(windows_in_cell)
    % windows -- (image_id_in_total, image_imdb_idx, image_id_in_imdb,
    % label, ov, l, t, r, b)

    image_imdb_idx = cellfun(@(x, y) ones(size(y, 1), 1) * x, num2cell(1:length(windows_in_cell)), windows_in_cell, 'UniformOutput', false);
    image_imdb_idx = cat(1, image_imdb_idx{:});
    windows = cat(1, windows_in_cell{:});
    windows = cellfun(@(x, y, z) [ones(size(z, 1), 1)*x, ones(size(z, 1), 1)*y, z], ...
        num2cell(1:size(windows, 1))', num2cell(image_imdb_idx), windows, 'UniformOutput', false);

end

function InitSolver(net_def_file, net_file, log_file)
    cur_dir = pwd;
    cd(fileparts(net_def_file));
    caffe('init_solver', net_def_file, net_file, log_file);
    cd(cur_dir);
end

author: ShaoqingRen


text:
%-----------------------------------------------------------------------
% Job saved on 01-Dec-2015 13:23:09 by cfg_util (rev $Rev: 6627 $)
% spm SPM - SPM12 (12.1)
% cfg_basicio BasicIO - Unknown
% dtijobs DTI tools - Unknown
% impexp_NiftiMrStruct NiftiMrStruct - Unknown
% cfg_logextract LogExt - Unknown
% tbx_dpc dPC Tools - Unknown
% prt PRoNTo - Unknown
% perfusion Perfusion - Unknown
% menu_cfg ConfGUI - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.menu_cfg.menu_entry.conf_files.type = 'cfg_files';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.name = 'File Set';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.tag = 'files';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.filter = '\.gz$';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.ufilter = '.*';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.dir = '';
matlabbatch{1}.menu_cfg.menu_entry.conf_files.num = [0 Inf];
matlabbatch{1}.menu_cfg.menu_entry.conf_files.check = [];
matlabbatch{1}.menu_cfg.menu_entry.conf_files.rewrite_job = [];
matlabbatch{1}.menu_cfg.menu_entry.conf_files.help = {'Select a set of files.'};
matlabbatch{1}.menu_cfg.menu_entry.conf_files.def = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_files.type = 'cfg_files';
matlabbatch{2}.menu_cfg.menu_entry.conf_files.name = 'Output directory';
matlabbatch{2}.menu_cfg.menu_entry.conf_files.tag = 'outdir';
matlabbatch{2}.menu_cfg.menu_entry.conf_files.filter = 'dir';
matlabbatch{2}.menu_cfg.menu_entry.conf_files.ufilter = '.*';
matlabbatch{2}.menu_cfg.menu_entry.conf_files.dir = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_files.num = [0 1];
matlabbatch{2}.menu_cfg.menu_entry.conf_files.check = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_files.rewrite_job = [];
matlabbatch{2}.menu_cfg.menu_entry.conf_files.help = {'Output files will be placed in this folder. Leave empty to put them into the same folder as the original files.'};
matlabbatch{2}.menu_cfg.menu_entry.conf_files.def = [];
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.type = 'cfg_menu';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.name = 'Keep original files';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.tag = 'keep';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.labels = {
                                                       'Yes'
                                                       'No'
                                                       }';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.values = {
                                                       true
                                                       false
                                                       }';
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.check = [];
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.rewrite_job = [];
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.help = {};
matlabbatch{3}.menu_cfg.menu_entry.conf_menu.def = [];
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.type = 'cfg_exbranch';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.name = 'Gunzip Files';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.tag = 'cfg_gunzip_files';
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{1}(1) = cfg_dep('Files: File Set (cfg_files)', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('()',{1}));
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{2}(1) = cfg_dep('Files: Output directory (cfg_files)', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('()',{1}));
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.val{3}(1) = cfg_dep('Menu: Keep original files (cfg_menu)', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('()',{1}));
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.prog = @cfg_run_gunzip_files;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.vout = @cfg_vout_gunzip_files;
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.check = [];
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.rewrite_job = [];
matlabbatch{4}.menu_cfg.menu_struct.conf_exbranch.help = {'Gunzip each file in a set of files.'};

author: spm


text:
function t_pnes_master(quiet)
%T_PNES_MASTER  Tests of PNE solvers via PNES_MASTER().

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

if nargin < 1
    quiet = 0;
end

%%  alg         name        check       opts
cfg = {
    {'DEFAULT', 'default',  []          []  },
};

n = 188;

t_begin(n*length(cfg), quiet);

for k = 1:length(cfg)
    alg   = cfg{k}{1};
    name  = cfg{k}{2};
    check = cfg{k}{3};
    opts  = cfg{k}{4};
    if ~isempty(check) && ~have_feature(check)
        t_skip(n, sprintf('%s not installed', name));
    else
        %% default start point
        x0 = [-1;0;0];
        opt = struct( ...
            'verbose', 0, ...
            'alg', alg, ...
            'stop_at', 0.7, ...
            'parameterization', 1, ...
            'nleqs_opt', struct('tol', 1e-9), ...
            'adapt_step', 0, ...
            'step_max', 10, ...
            'target_lam_tol', 1e-6, ...
            'nose_tol', 1e-6, ...
            'adapt_step_tol', 1e-2 );
%         opt.plot = struct('level', 2, 'idx', 2);

        t = sprintf('%s - TARGET_LAM = 0.7 (natural) : ', name);
        [x, f, e, out, jac] = pnes_master(@f1p, x0, opt);
        it = 14;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-1.931782106; -1.268217894; 0.7], 6, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 0.7, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'TARGET_LAM', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Reached desired lambda 0.7 in %d continuation steps', it), [t 'out.done_msg']);
        t_ok(isfield(out, 'corrector'), [t 'out.corrector exists']);
        t_is(out.corrector.iterations, 3, 12, [t 'out.corrector.iterations']);
        t_str_match(out.corrector.message, sprintf('Newton''s method converged in %d iterations.', 3), [t 'out.corrector.message']);

        t = sprintf('%s - TARGET_LAM = 0.7 (arc len) : ', name);
        opt.adapt_step = 1;
        opt.parameterization = 2;
        [x, f, e, out, jac] = pnes_master(@f1p, x0, opt);
        it = 9;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-1.931782106; -1.268217894; 0.7], 6, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 0.7, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'TARGET_LAM', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Reached desired lambda 0.7 in %d continuation steps', it), [t 'out.done_msg']);

        t = sprintf('%s - TARGET_LAM = 0.7 (pseudo arc len) : ', name);
        opt.parameterization = 3;
        [x, f, e, out, jac] = pnes_master(@f1p, x0, opt);
        it = 9;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-1.931782106; -1.268217894; 0.7], 6, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 0.7, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'TARGET_LAM', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Reached desired lambda 0.7 in %d continuation steps', it), [t 'out.done_msg']);

        t = sprintf('%s - FULL : ', name);
        opt.stop_at = 'FULL';
        p = struct('fcn', @f1p, 'x0', x0, 'opt', opt);
        [x, f, e, out, jac] = pnes_master(p);
        it = 34;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [2;-1;0], 10, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 1.04127275, 8, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'TARGET_LAM', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Traced full continuation curve in %d continuation steps', it), [t 'out.done_msg']);

        t = sprintf('%s - FULL (max_it) : ', name);
        opt.max_it = 25;
        p = struct('fcn', @f1p, 'x0', x0, 'opt', opt);
        [x, f, e, out, jac] = pnes_master(p);
        opt.max_it = 2000;
        it = 25;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [0.64152370; -4.588447338; 0.82448727], 8, [t 'x - final']);
        t_is(f, [0;0], 8, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 1.04127275, 8, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 0, 12, [t 'length(out.events)']);
        t_str_match(out.done_msg, sprintf('Reached maximun number of continuation steps (opt.max_it = %d)', it), [t 'out.done_msg']);

        t = sprintf('%s - NOSE (arc len): ', name);
        p.opt.stop_at = 'NOSE';
        p.opt.parameterization = 2;
        [x, f, e, out, jac] = pnes_master(p);
        it = 18;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-0.5; -4.75; 1.04166667], 6, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 1.04166666667, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'NOSE', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Reached limit in %d continuation steps, lambda = 1.042.', it), [t 'out.done_msg']);

        t = sprintf('%s - NOSE (pseudo arc len) : ', name);
        p.opt.parameterization = 3;
        [x, f, e, out, jac] = pnes_master(p);
        it = 18;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-0.5; -4.75; 1.04166667], 6, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [-3;4;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 1.04166666667, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'NOSE', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Reached limit in %d continuation steps, lambda = 1.042.', it), [t 'out.done_msg']);

        t = sprintf('%s - NOSE (opp dir) : ', name);
        p.x0 = [1;-1;0];
        [x, f, e, out, jac] = pnes_master(p);
        it = 20;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-0.5; -4.75; 1.04166667], 5, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_is(out.x(:,1), [2;-1;0], 8, [t 'out.x(:,1)']);
        t_is(out.max_lam, 1.04166666667, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.x), [3,it+1], 12, [t 'size(out.x)']);
        t_is(size(out.x_hat), [3,it+1], 12, [t 'size(out.x_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 1, 12, [t 'length(out.events)']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'NOSE', [t 'out.events.name']);
        t_str_match(out.done_msg, sprintf('Reached limit in %d continuation steps, lambda = 1.042.', it), [t 'out.done_msg']);

        t = sprintf('%s - FULL warmstart, before SWITCH : ', name);
        opt.parameterization = 3;
        opt.adapt_step = 1;
        opt.adapt_step_ws = 0.25;
        opt.stop_at = 'FULL';
        opt.adapt_step_tol = 5e-3;
        opt.callbacks = {@pne_callback_test1};
        opt.events = {{'SWITCH!', @pne_event_test1, 1e-6}};
        opt.output_fcn = @pne_output_fcn_test1;
        % opt.verbose = 4;
        % opt.plot = struct('level', 2, 'yname', 'y', 'idx', [1;2]);
        x0 = [-1;0;0];
        [x, f, e, out, jac] = pnes_master(@f1p, x0, opt);
        it = 10;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [-2; -1; 2/3], 6, [t 'x - final']);
        t_is(f, [0;0], 10, [t 'f']);
        t_ok(isfield(out, 'warmstart') && ~isempty(out.warmstart), [t 'out.warmstart exists']);
        t_is(out.max_lam, 2/3, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.y), [2,it+1], 12, [t 'size(out.y)']);
        t_is(size(out.y_hat), [2,it+1], 12, [t 'size(out.y_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_ok(isstruct(out.events), [t 'out.events is struct']);
        t_is(out.events.k, it, 12, [t 'out.events.k']);
        t_is(out.events.idx, 1, 12, [t 'out.events.idx']);
        t_str_match(out.events.name, 'SWITCH!', [t 'out.events.name']);
        t_str_match(out.events.msg, 'ZERO detected for SWITCH! event', [t 'out.events.msg']);
        t_is(out.warmstart.cont_steps, it, 12, [t 'out.warmstart.cont_steps']);
        t_is(out.warmstart.default_step, 0.64432407, 8, [t 'out.warmstart.default_step']);
        t_ok(isa(out.warmstart.parm, 'function_handle'), [t 'out.warmstart.parm is function']);
        t_ok(isa(out.warmstart.default_parm, 'function_handle'), [t 'out.warmstart.default_parm is function']);
        t_ok(isstruct(out.warmstart.cbs), [t 'out.warmstart.cbs is struct']);
        t_ok(isstruct(out.warmstart.cbs.default), [t 'out.warmstart.cbs.default is struct']);
        t_is(out.warmstart.cbs.default.iterations, it, 12, [t 'out.warmstart.cbs.default.iterations']);
        t_ok(isstruct(out.warmstart.events), [t 'out.warmstart.events is struct']);
        t_is(out.warmstart.events.k, it, 12, [t 'out.warmstart.events.k']);
        t_is(out.warmstart.events.idx, 1, 12, [t 'out.warmstart.events.idx']);
        t_str_match(out.warmstart.events.name, 'SWITCH!', [t 'out.warmstart.events.name']);
        t_str_match(out.warmstart.events.msg, 'ZERO detected for SWITCH! event', [t 'out.warmstart.events.msg']);
        t_str_match(out.done_msg, sprintf('Reached switching point in %d continuation steps', it), [t 'out.done_msg']);
        t_ok(isfield(out, 'corrector'), [t 'out.corrector exists']);
        t_is(out.corrector.iterations, 3, 12, [t 'out.corrector.iterations']);
        t_str_match(out.corrector.message, sprintf('Newton''s method converged in %d iterations.', 3), [t 'out.corrector.message']);

        t = sprintf('%s - FULL warmstart, after SWITCH : ', name);
        ws = out.warmstart;
        ws.x  = [ws.x(end)/2;  ws.x];
        ws.xp = [ws.xp(end)/2; ws.xp];
        ws.z = [0; ws.z];
        ws.zp = [0; ws.zp];
        x0 = ws.x;      %% ignored for warm start
        opt.warmstart = ws;
        opt.events = {{'SNOUT!', @pne_event_test2, 1e-6}};
        opt.callbacks = {};
        [x, f, e, out, jac] = pnes_master(@f2p, x0, opt);
        it = 49;
        t_is(e, 1, 12, [t 'exitflag']);
        t_is(x, [0; 6; -5; 0], 6, [t 'x - final']);
        t_is(f, [0;0;0], 10, [t 'f']);
        t_is(out.y(:,1), [-3;4], 8, [t 'out.y(:,1)']);
        t_is(out.max_lam, 1.04166666667, 10, [t 'out.max_lam']);
        t_is(out.iterations, it, 12, [t 'out.iterations']);
        t_is(size(out.lam), [1,it+1], 12, [t 'size(out.lam)']);
        t_is(size(out.lam_hat), [1,it+1], 12, [t 'size(out.lam_hat)']);
        t_is(size(out.y), [2,it+1], 12, [t 'size(out.y)']);
        t_is(size(out.y_hat), [2,it+1], 12, [t 'size(out.y_hat)']);
        t_is(size(out.steps), [1,it+1], 12, [t 'size(out.steps)']);
        t_is(length(out.events), 3, 12, [t 'length(out.events)']);
        t_is(out.events(1).k, 10, 12, [t 'out.events(1).k']);
        t_is(out.events(1).idx, 1, 12, [t 'out.events(1).idx']);
        t_str_match(out.events(1).name, 'SWITCH!', [t 'out.events(1).name']);
        t_str_match(out.events(1).msg, 'ZERO detected for SWITCH! event', [t 'out.events(1).msg']);
        t_is(out.events(2).k, 28, 12, [t 'out.events(2).k']);
        t_is(out.events(2).idx, 1, 12, [t 'out.events(2).idx']);
        t_str_match(out.events(2).name, 'SNOUT!', [t 'out.events(2).name']);
        t_str_match(out.events(2).msg, 'ZERO detected for SNOUT! event', [t 'out.events(2).msg']);
        t_is(out.events(3).k, it, 12, [t 'out.events(3).k']);
        t_is(out.events(3).idx, 1, 12, [t 'out.events(3).idx']);
        t_str_match(out.events(3).name, 'TARGET_LAM', [t 'out.events(3).name']);
        t_str_match(out.events(3).msg, 'ZERO detected for TARGET_LAM event', [t 'out.events(3).msg']);
        t_str_match(out.done_msg, sprintf('Traced full continuation curve in %d continuation steps', it), [t 'out.done_msg']);
        t_ok(isfield(out, 'corrector'), [t 'out.corrector exists']);
        t_is(out.corrector.iterations, 3, 12, [t 'out.corrector.iterations']);
        t_str_match(out.corrector.message, sprintf('Newton''s method converged in %d iterations.', 3), [t 'out.corrector.message']);
    end
end

t_end;

% lam = 1;
% m = 6;
% xs = [-m:0.2:m];
% ys = xs;
% [xx, yy] = meshgrid(xs, ys);
% zz0 = zeros([size(xx), 2]);
% zz = zz0;
% for i = 1:length(xs)
%     for j = 1:length(ys)
%         zz(i, j, :) = f1p([xx(i, j); yy(i, j); lam]);
% %         zz(i, j, :) = f2([xx(i, j); yy(i, j)]);
%     end
% end
%
% figure
% ax = gca;
% surf(ax, xx, yy, squeeze(zz(:, :, 1)))
% hold on
% surf(ax, xx, yy, squeeze(zz(:, :, 2)))
% surf(ax, m*[-1 -1; 1 1], m*[-1 1; -1 1], [0 0; 0 0])
% hold off

%% 2-d problem with 2 solutions
%% from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/
function [f, J] = f1(x)
f = [  x(1)   + x(2) - 1;
      -x(1)^2 + x(2) + 5    ];
if nargout > 1
    J = [1 1; -2*x(1) 1];
end

%% parameterized 2-d problem
%% based on https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/
function [f, J] = f1p(x)
if nargout < 2
    f = f1(x(1:2));
else
    [f, J] = f1(x(1:2));
    J = [J [6;0]];
end
f = f + [6*x(3); 0];

%% parameterized 3-d problem, based on f1p
function [f, J] = f2p(x)
x1 = [x(3)+2; x(2)-2; x(4)];
if nargout < 2
    f = [2*x(1)-x(4); f1p(x1)];
else
    [f, JJ] = f1p(x1);
    J = [  2      0        0       -1;
          [0; 0] JJ(:, 2) JJ(:, 1) JJ(:, 3) ];
end
f = [ 2*x(1)-x(4); f ];

%% example custom event function 1 (target lambda == 2/3)
function efv = pne_event_test1(cx, opt)
tlam = 2/3;
efv = cx.x(end) - tlam;

%% example custom event function 2 (nose point)
function efv = pne_event_test2(cx, opt)
efv = cx.z(end);

%% example custom callback function (exit for warmstart on SWITCH! event)
function [nx, cx, s] = pne_callback_test1(k, nx, cx, px, s, opt)
if k <= 0 || s.done, return; end    %% skip if initialize, finalize or done
tlam = 2/3;
ev = pne_detected_event(s.events, 'SWITCH!');
if ~isempty(ev)
    if ev.zero              %% prepare to terminate
        s.done = 1;
        s.done_msg = sprintf('Reached switching point in %d continuation steps', k);
        s.warmstart = struct(); %% signal that we want to exit, then resume
    else                    %% set step-size & parameterization to terminate next time
        cx.this_parm = @pne_pfcn_natural;   %% change to natural parameterization
        cx.this_step = tlam - cx.x(end);
        ev.msg = sprintf('%s\n  step %d to overshoot lambda = %g, reduce step size and set natural param', ev.msg, k, tlam);
    end
end

function [names, vals] = pne_output_fcn_test1(x, x_hat)
%% [names, vals] = pne_output_fcn_default(x, x_hat)
%% names = pne_output_fcn_default()
names = {'y_hat', 'y'};
if nargin
    if length(x) == 3
        k = [1;2];
    else %% == 3
        k = [2;3];
    end
    vals = {x_hat(k), x(k)};
end

% %% another 2-d problem
% %% from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg
% function [f, J] = f2(x)
% m = 10;
% m1 = 0; %2;
% m2 = 0; %5;
% f = [  x(1)^2 +   x(1)*x(2)   - 10 + m1;
%        (x(2)   + 3*x(1)*x(2)^2 - 57) / m  ] + m2;
% if nargout > 1
%     J = sparse([    2*x(1)+x(2)    x(1);
%                     (3*x(2)^2) / m (6*x(1)*x(2)+1) / m    ]);
% end


text:
classdef CachedCalculationResults
    %CACHEDCALCULATIONRESULTS Summary of this class goes here
    %   Detailed explanation goes here
    % Usage:
    %     cached_params = [];
    %     cached_params.Y = Y;
    %     cached_params.k = k;
    %     cached_file = '../../../../CachedResults/cached_nmf_adjacency_graph.mat';
    %     ccr = CachedCalculationResults(cached_file);
    %     A = ccr.GetMatchedResult(cached_params);
    %     if isempty(A)
    %         [A,~,~] = create_adjacency_graph(Y,'nn',k,1);
    %         if N*B > 1e6
    %             ccr.AddResult(cached_params,A);
    %         end
    %     end

    properties
        FileName
        MaxSignatureSize
        MaxEntries
    end

    methods
        function obj = CachedCalculationResults(file_name,max_param_size)
            if nargin < 2
                max_param_size = 1e5;
            end
            obj.FileName = file_name;
            obj.MaxSignatureSize = max_param_size;
            obj.MaxEntries = 50;
        end

        function AddResult(obj, params, result)
            [list_params,list_results] = LoadCachedResults(obj);
            if length(list_params) >= obj.MaxEntries
                list_params(1) = [];
                list_results(1) = [];
            end
            list_params = [list_params,{obj.GetSignature(params)}];
            list_results = [list_results,{result}];
            save(obj.FileName,'list_params','list_results','-v7.3');
        end

        function matched = GetMatchedResult(obj,params)
            [list_params,list_results] = LoadCachedResults(obj);
            signature = obj.GetSignature(params);

            matched = [];
            for i = 1:length(list_params)
                if isequal(list_params{i},signature)
                    matched = list_results{i};
                    break;
                end
            end

        end

        function [list_params,list_results] = LoadCachedResults(obj)
            if ~exist(obj.FileName,'file')
                list_params = {};
                list_results = {};
                save(obj.FileName,'list_params','list_results');
            end
            S = load(obj.FileName);
            list_params = S.list_params;
            list_results = S.list_results;
        end

        function new_params = GetSignature(obj,params)
            new_params = params;
            names = fieldnames(new_params);
            for i = 1:length(names)
                X = new_params.(names{i});
                x = X(:);
                if length(x) > obj.MaxSignatureSize
                    step = ceil(length(x) / obj.MaxSignatureSize);
                    x = x(1:step:end);
                end
                new_params.(names{i}) = x;
            end
        end

    end

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function c = acorr(y, maxlag)

if isvector(y)
  if nargin < 2
    maxlag = length(y);
  end
  y = y - mean(y);
  z = xcorr(y,maxlag,'coeff');
  l0 = (length(z)+1)/2;
  c = z(l0:end);
else
  if nargin < 2
    maxlag = size(y,1);
  end
  y = bsxfun(@minus, y, mean(y,1));
  c = zeros(maxlag+1,size(y,2));
  for n=1:size(y,2)
    z = xcorr(y(:,n),maxlag,'coeff');
    l0 = (size(z,1)+1)/2;
    c(:,n) = z(l0:end);
  end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%show_dose_temp.m

global planC
indexS = planC{end};

planC{indexS.IM}(1).IMDosimetry.beams = [IMwQIB.beams];

planC{indexS.IM}.IMDosimetry.goals = IMwQIB(1).goals;

planC{indexS.IM}.IMDosimetry.params = IMwQIB(1).params;

planC{indexS.IM}.IMDosimetry.isFresh = 1;

planC{indexS.IM}.IMDosimetry.name = 'Eclipse Leaf Seq';

planC{indexS.IM}.IMDosimetry.solutions = [];

planC{indexS.IM}.IMDosimetry.assocScanUID = planC{indexS.scan}(1).scanUID;

planC{indexS.IM}.IMDosimetry.solutions = [beamletWeightC{:}];

[dose3D] = getIMDose(planC{indexS.IM}(1).IMDosimetry, [beamletWeightC{:}], structNumsV);

showIMDose(dose3D,'QIB reCalc')

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
clear;
clc;
close all;

addpath('../../library/nav_lib');

%% INPUTS: INITIAL CONDITION
pos_user_true = [1000, 100]';
pos_user_inital = [0, 0]';
pos_user = pos_user_inital;
anchor_pos  = [0,1000; 0 -1000; 2000, 100]';
pr = vecnorm(anchor_pos - pos_user_true);

 pos_user = ch_multilateration(anchor_pos, pos_user , pr, 2);

pos_user

author: yandld


text:
function [x,I,exitFlag]=getMinNrFluxes(model, toMinimize, params,scores)
% getMinNrFluxes
%   Returns the minimal set of fluxes that satisfy the model using
%   mixed integer linear programming.
%
%	model         a model structure
%   toMinimize    either a cell array of reaction IDs, a logical vector
%                 with the same number of elements as reactions in the model,
%                 of a vector of indexes for the reactions that should be
%                 minimized (opt, default model.rxns)
%   params        parameter structure as used by getMILPParams (opt)
%   scores        vector of weights for the reactions. Negative scores
%                 should not have flux. Positive scores are not possible in this
%                 implementation, and they are changed to max(scores(scores<0)).
%                 Must have the same dimension as toMinimize (find(toMinimize)
%                 if it is a logical vector) (opt, default -1 for all reactions)
%
%   x             the corresponding fluxes for the full model
%   I             the indexes of the reactions in toMinimize that were used
%                 in the solution
%   exitFlag      1: optimal solution found
%                -1: no feasible solution found
%                -2: optimization time out
%
%   NOTE: Uses 1000 mmol/gDW/h as an arbitary large flux. Could possibly
%   cause problems if the fluxes in the model are larger than that.
%
%   Usage: [x,I,exitFlag]=getMinNrFluxes(model, toMinimize, params, scores)

exitFlag=1;

if nargin<2
    toMinimize=model.rxns;
elseif ~islogical(toMinimize) && ~isnumeric(toMinimize)
    toMinimize=convertCharArray(toMinimize);
else
    toMinimize=model.rxns(toMinimize);
end

%For passing parameters to the solver
if nargin<3
    params=struct();
end

if nargin<4
    %It says that the default is -1, but that is to fit with other code
    scores=ones(numel(toMinimize),1)*1;
else
    if numel(scores)~=numel(toMinimize)
        EM='The number of scores must be the same as the number of reactions to minimize';
        dispEM(EM);
    end

    %Change positive scores to have a small negative weight. This is a
    %temporary solution.
    scores(scores>=0)=max(scores(scores<0));

    %It says that the default is -1, but that is to fit with other code
    scores=scores*-1;
end

%Check if the model is in irreversible format
if any(model.rev)
    %Convert the model to irreversible format
    irrevModel=convertToIrrev(model);

    %Find the indexes for the reactions in toMinimize
    [indexes, I]=ismember(strrep(irrevModel.rxns,'_REV',''),toMinimize);
else
    irrevModel=model;

    %Find the indexes for the reactions in toMinimize
    [indexes, I]=ismember(irrevModel.rxns,toMinimize);
end

indexes=find(indexes);
%Adjust scores to fit with reversible
scores=scores(I(indexes));

%Add binary constraints in the following manner: -  Add one unique
%"metabolite" for each integer reaction as a substrate.
%   These metabolites can have net production
%-  Add reactions for the production of each of those metabolites. The
%   amount produced in one reaction unit must be larger than the largest
%   possible flux in the model (but not too large to avoid bad scaling)

%Calculate a solution to the problem without any constraints. This is to
%get an estimate about the magnitude of fluxes in the model and to get a
%feasible start solution.
sol=solveLP(irrevModel,1);

%Return an empty solution if the non-constrained problem couldn't be solved
if isempty(sol.x)
    x=[];
    I=[];
    exitFlag=-1;
    return;
end

%Take the maximal times 5 to have a safe margin. If it's smaller than 1000,
%then use 1000 instead.
maxFlux=max(max(sol.x)*5,1000);

intArray=speye(numel(irrevModel.rxns))*-1;
intArray=intArray(indexes,:);
prob.a=[irrevModel.S;intArray];
a=[sparse(numel(irrevModel.mets),numel(indexes));speye(numel(indexes))*maxFlux];
prob.a=[prob.a a];
prob.ints.sub=numel(irrevModel.rxns)+1:numel(irrevModel.rxns)+numel(indexes);

prob.c=[zeros(numel(irrevModel.rxns),1);scores(:);zeros(size(prob.a,1),1)]; %Minimize the number of fluxes
prob.A=[prob.a -speye(size(prob.a,1))];
prob.blc=[irrevModel.b(:,1);zeros(numel(indexes),1)];
if size(irrevModel.b,2)==2
    prob.buc=[irrevModel.b(:,2);inf(numel(indexes),1)];
else
    prob.buc=[irrevModel.b(:,1);inf(numel(indexes),1)];
end
prob.blx=[irrevModel.lb;zeros(numel(indexes),1)];
prob.bux=[irrevModel.ub;ones(numel(indexes),1)];
prob.lb = [prob.blx; prob.blc];
prob.ub = [prob.bux; prob.buc];
prob.osense=1;
prob.csense=repmat('E', 1, size(prob.a,1),1);
prob.b=zeros(size(prob.a,1), 1);

%Use the output from the linear solution as starting point. Only the values
%for the integer variables will be used, but all are supplied.
prob.sol.int.xx=zeros(numel(prob.c),1);
prob.sol.int.xx(prob.ints.sub(sol.x(indexes)>10^-12))=1;
prob.x0=[];
prob.vartype=repmat('C', size(prob.A,2), 1);
prob.vartype(prob.ints.sub) = 'I'; % with .lb = 0 and .ub = 1, they are binary
% integers (glpk in octave only allows 'continuous' or '', not 'binary')
prob=rmfield(prob,{'blx','bux','blc','buc'});

% Optimize the problem
res = optimizeProb(prob,params);
isFeasible=checkSolution(res);

if ~isFeasible
    x=[];
    I=[];
    exitFlag=-1;
    return;
end

xx=res.full(1:numel(irrevModel.rxns));
I=res.full(numel(xx)+1:end);

%Check if Mosek aborted because it reached the time limit
%TODO: modify for cobra/gurobi
% if strcmp('MSK_RES_TRM_MAX_TIME',res.rcode)
%     exitFlag=-2;
% end

%Map back to original model from irrevModel
x=xx(1:numel(model.rxns));
if numel(irrevModel.rxns)>numel(model.rxns)
    x(model.rev~=0)=x(model.rev~=0)-xx(numel(model.rxns)+1:end);
end

I=ismember(toMinimize,strrep(irrevModel.rxns(indexes(I>10^-12)),'_REV',''));
end

author: SysBioChalmers

repo: RAVEN