text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function run_TLD

dir_names = {'06_car'};

for i = 1:numel(dir_names)
    close all;
    run_TLD_one(dir_names{i});
end

function run_TLD_one(seq_name)

addpath(genpath('.'));
init_workspace;

src_dir = sprintf('/home/yuxiang/Projects/Multiview_Tracking/dataset/TLD/%s/img/', seq_name);
dst_dir = sprintf('/home/yuxiang/Projects/Multiview_Tracking/result/OpenTLD/%s/', seq_name);

opt.source          = struct('camera',0,'input',src_dir,'bb0', []); % camera/directory swith, directory_name, initial_bounding_box (if empty, it will be selected by the user)
opt.output          = dst_dir; mkdir(opt.output); % output directory that will contain bounding boxes + confidence

min_win             = 24; % minimal size of the object's bounding box in the scanning grid, it may significantly influence speed of TLD, set it to minimal size of the object
patchsize           = [15 15]; % size of normalized patch in the object detector, larger sizes increase discriminability, must be square
fliplr              = 0; % if set to one, the model automatically learns mirrored versions of the object
maxbbox             = 1; % fraction of evaluated bounding boxes in every frame, maxbox = 0 means detector is truned off, if you don't care about speed set it to 1
update_detector     = 1; % online learning on/off, of 0 detector is trained only in the first frame and then remains fixed
opt.plot            = struct('pex',1,'nex',1,'dt',1,'confidence',1,'target',1,'replace',0,'drawoutput',3,'draw',0,'pts',1,'help', 0,'patch_rescale',1,'save',0);

% Do-not-change -----------------------------------------------------------

opt.model           = struct('min_win',min_win,'patchsize',patchsize,'fliplr',fliplr,'ncc_thesame',0.95,'valid',0.5,'num_trees',10,'num_features',13,'thr_fern',0.5,'thr_nn',0.65,'thr_nn_valid',0.7);
opt.p_par_init      = struct('num_closest',10,'num_warps',20,'noise',5,'angle',20,'shift',0.02,'scale',0.02); % synthesis of positive examples during initialization
opt.p_par_update    = struct('num_closest',10,'num_warps',10,'noise',5,'angle',10,'shift',0.02,'scale',0.02); % synthesis of positive examples during update
opt.n_par           = struct('overlap',0.2,'num_patches',100); % negative examples initialization/update
opt.tracker         = struct('occlusion',10);
opt.control         = struct('maxbbox',maxbbox,'update_detector',update_detector,'drop_img',1,'repeat',1);

% Run TLD -----------------------------------------------------------------
%profile on;
[bb,conf] = tldExample(opt);
%profile off;
%profile viewer;

% Save results ------------------------------------------------------------
dlmwrite([opt.output '/tld.txt'],[bb; conf]');
disp(['Results saved to ' dst_dir]);

author: yuxng


text:
% Ambiguity function
%
% Computes the ambiguity function of the input signal. An analytic
% signal generator is called if the input signal is real. This version
% takes advantage of the property K(n,-m)=K*(n,m) to increase
% calculation speed.
%
% Usage:
%
%     af = ambf( signal, lag_res, window_length )
%
%
% Parameters:
%
%     af
%
%	  The computed ambiguity function representation. size(af) will
%	  return [a, b], where b is the next largest power of two above
%	  (2*window_length), and a is equal signal length.
%
%    signal
%
%	  Input one dimensional signal to be analysed. An analytic signal
%	  is required for this function, however, if signal is real, a
%	  default analytic transformer routine will be called from this
%	  function before computing the ambiguity function.
%
%    window_length
%
%	  The number of signal samples used for analysis.
%
%    lag_res
%
%	  The number of lag to skip between successive slices of
%	  the analysis.
%
%
%
%  See Also: analyt
%
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
% Test file for @chebfun/minus.m.

function pass = test_minus(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

% A random number to use as an arbitrary additive constant.
alpha = -0.194758928283640 + 0.075474485412665i;

% Check behavior for empty arguments.
f = chebfun(@(x) sin(x), pref);
g = chebfun();
pass(1) = isempty(f - []);
pass(2) = isempty(f - g);

% Turn on splitting, since we'll need it for the rest of the tests.
pref.splitting = 1;

% Test addition with scalars.
f1_op = @(x) sin(x).*abs(x - 0.1);
f1 = chebfun(f1_op, pref);
%keyboard;
pass(3:4) = test_sub_function_and_scalar(f1, f1_op, alpha, x);

% Test subition of two chebfun objects.
g1_op = @(x) cos(x).*sign(x + 0.2);
g1 = chebfun(g1_op, pref);
pass(5:6) = test_sub_function_and_function(f1, f1_op, g1, g1_op, x);

% Test operation for array-valued chebfuns.
f2_op = @(x) [sin(x).*abs(x - 0.1)  exp(x)];
f2 = chebfun(f2_op, pref);
pass(7:8) = test_sub_function_and_scalar(f2, f2_op, alpha, x);

g2_op = @(x) [cos(x).*sign(x + 0.2) tan(x)];
g2 = chebfun(g2_op, pref);
pass(9:10) = test_sub_function_and_function(f2, f2_op, g2, g2_op, x);

% Test operation for transposed chebfuns.
pass(11:12) = test_sub_function_and_scalar(f1.', @(x) f1_op(x).', alpha, x);
pass(13:14) = test_sub_function_and_function(f1.', @(x) f1_op(x).', ...
    g1.', @(x) g1_op(x).', x);

% Check error conditions.
try
    h = f1 - uint8(128);
    pass(15) = strcmp(ME.identifier, 'CHEBFUN:plus:unknown')
catch ME
    pass(15) = true;
end

try
    h = f1 - g1.'
    pass(16) = strcmp(ME.identifier, 'CHEBFUN:plus:matdim')
catch ME
    pass(16) = true;
end

%% Tests for function defined on unbounded domain:

% Functions on [-inf b]:

% Set the domain:
dom = [-Inf 3*pi];
domCheck = [-1e6 3*pi];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

opf = @(x) x.*exp(x);
opg = @(x) (1-exp(x))./x;
oph = @(x) x.*exp(x) - (1-exp(x))./x;
f = chebfun(opf, dom);
g = chebfun(opg, dom);
h = f - g;

hVals = feval(h, x);
hExact = oph(x);
err = hVals - hExact;
pass(17) = norm(err, inf) < 1e1*eps*get(h,'vscale');

end

% Test the subtraction of a CHEBFUN F, specified by F_OP, to and from a scalar
% ALPHA using a grid of points X in the domain of F for testing samples.
function result = test_sub_function_and_scalar(f, f_op, alpha, x)
    g1 = f - alpha;
    g2 = alpha - f;
    result(1) = isequal(g1, -g2);
    g_exact = @(x) f_op(x) - alpha;
    result(2) = norm(feval(g1, x) - g_exact(x), inf) < 100*vscale(g1)*eps;
end

% Test the subraction of two CHEBFUN objects F and G, specified by F_OP and
% G_OP, using a grid of points X in the domains of F and G for testing samples.
function result = test_sub_function_and_function(f, f_op, g, g_op, x)
    h1 = f - g;
    h2 = g - f;
    result(1) = isequal(h1, -h2);
    h_exact = @(x) f_op(x) - g_op(x);
    norm(feval(h1, x) - h_exact(x), inf);
    result(2) = norm(feval(h1, x) - h_exact(x), inf) < 100*vscale(h1)*eps;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  qdd = call_direct_dynamics_fast(input)
%  Auxiliar function for the simulink model SIMULATE_ROBOT_AND_CONTROLLER.
%  Rearranges the inputs coming from the simulink model and calls the
%  function accel.
%  As a result the instantaneous acceleration at each joint is returned.
%
%  See also ACCEL.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function qdd = call_direct_dynamics_2dofplanar(input)

global robot

torque = input(1:2);   % Input torque at each joint
q   = input(3:4);	   % Joint positions
qd  = input(5:6);	   % Joint speeds
fe=[0 0 0 0 0 0]'; %external forces applied
torque

% Compute acceleration
qdd = forwarddynamics_2dofplanar(robot, q, qd, torque, 9.81, fe)

repo: ARTE


text:
function model = LogisticRegressionTrain(dataPos, dataNeg, opt, model)

featPos         = dataPos.feat;
featNeg         = dataNeg.feat;
numPos          = size(featPos, 2);
numNeg          = size(featNeg, 2);
featDim         = size(featPos, 1);
num             = numPos + numNeg;

feat            = [featPos, featNeg];
data.feat       = feat;

feat            = [feat; ones(1, num)];     % add bias term

label           = zeros(1, num);
label(1:numPos) = 1;

if nargin <= 3
    model.w     = randn(featDim + 1, 1) / (featDim + 1);
    maxIter         = 1000;
else
    maxIter = 20;
end

lambda          = 1e-2;
iter            = 0;
lr              = 1 / num;
alpha           = 0.99; % Parameter for momentum
deltaW          = 0;
while (iter < maxIter)
    iter        = iter + 1;

    loss        = label - LogisticRegressionTest(data, model);
    deltaW       = lr * (feat * loss' - lambda * model.w) + alpha * deltaW;
%     deltaW       = lr * (feat * loss' - lambda * model.w);
    model.w     = model.w + deltaW;
%     disp(norm(loss));
end

end

author: flyers

repo: drone-tracking


text:
function [month, day, year] = gdate (jdate)

% convert Julian date to Gregorian (calendar) date

% input

%  jdate = julian day

% output

%  month = calendar month [1 - 12]
%  day   = calendar day [1 - 31]
%  year  = calendar year [yyyy]

%  note: day may include fractional part

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jd = jdate;

z = fix(jd + .5);
fday = jd + .5 - z;

if (fday < 0)
   fday = fday + 1;
   z = z - 1;
end

if (z < 2299161)
   a = z;
else
   alpha = floor((z - 1867216.25) / 36524.25);
   a = z + 1 + alpha - floor(alpha / 4);
end

b = a + 1524;
c = fix((b - 122.1) / 365.25);
d = fix(365.25 * c);
e = fix((b - d) / 30.6001);
day = b - d - fix(30.6001 * e) + fday;

if (e < 14)
   month = e - 1;
else
   month = e - 13;
end

if (month > 2)
   year = c - 4716;
else
   year = c - 4715;
end

author: Sable

repo: mcbench-benchmarks


text:
% Chapter 8 - Planar Systems.
% Program_8b - Phase Portrait (Fig. 8.8(a)).
% Copyright Birkhauser 2013. Stephen Lynch.

% Phase portrait of a linear system of ODE's.
% IMPORTANT - Program_8a is vectorfield.m.
clear;
% sys=inline('[2*x(1)+x(2);x(1)+2*x(2)]','t', 'x');
sys = @(t,x) [2*x(1)+x(2);x(1)+2*x(2)];
vectorfield(sys,-3:.25:3,-3:.25:3)
     hold on
     for x0=-3:1.5:3
         for y0=-3:1.5:3
            [ts,xs] = ode45(sys,[0 5],[x0 y0]);
            plot(xs(:,1),xs(:,2))
         end
     end
     for x0=-3:1.5:3
         for y0=-3:1.5:3
            [ts,xs] = ode45(sys,[0 -5],[x0 y0]);
            plot(xs(:,1),xs(:,2))
         end
     end
     hold off
axis([-3 3 -3 3])
fsize=15;
set(gca,'XTick',-3:1:3,'FontSize',fsize)
set(gca,'YTick',-3:1:3,'FontSize',fsize)
xlabel('x(t)','FontSize',fsize)
ylabel('y(t)','FontSize',fsize)
hold off

% End of Program_8b.


text:
function varargout=shadedErrorBar(x,y,errBar,varargin)
% generate continuous error bar area around a line plot
%
% function H=shadedErrorBar(x,y,errBar, ...)
%
% Purpose
% Makes a 2-d line plot with a pretty shaded error bar made
% using patch. Error bar color is chosen automatically.
%
%
% Inputs (required)
% x - vector of x values [optional, can be left empty]
% y - vector of y values or a matrix of n observations by m cases
%     where m has length(x);
% errBar - if a vector we draw symmetric errorbars. If it has a size
%          of [2,length(x)] then we draw asymmetric error bars with
%          row 1 being the upper bar and row 2 being the lower bar
%          (with respect to y). ** alternatively ** errBar can be a
%          cellArray of two function handles. The first defines which
%          statistic the line should be and the second defines the
%          error bar.
%
% Inputs (optional, param/value pairs)
% 'lineProps' - ['-k' by default] defines the properties of
%             the data line. e.g.:
%             'or-', or {'-or','markerfacecolor',[1,0.2,0.2]}
% 'transparent' - [true  by default] if true, the shaded error
%               bar is made transparent. However, for a transparent
%               vector image you will need to save as PDF, not EPS,
%               and set the figure renderer to "painters". An EPS
%               will only be transparent if you set the renderer
%               to OpenGL, however this makes a raster image.
% 'patchSaturation'- [0.2 by default] The saturation of the patch color.
%
%
%
% Outputs
% H - a structure of handles to the generated plot objects.
%
%
% Examples:
% y=randn(30,80);
% x=1:size(y,2);
%
% 1)
% shadedErrorBar(x,mean(y,1),std(y),'lineprops','g');
%
% 2)
% shadedErrorBar(x,y,{@median,@std},'lineprops',{'r-o','markerfacecolor','r'});
%
% 3)
% shadedErrorBar([],y,{@median,@(x) std(x)*1.96},'lineprops',{'r-o','markerfacecolor','k'});
%
% 4)
% Overlay two transparent lines:
% clf
% y=randn(30,80)*10;
% x=(1:size(y,2))-40;
% shadedErrorBar(x,y,{@mean,@std},'lineprops','-r','transparent',1);
% hold on
% y=ones(30,1)*x; y=y+0.06*y.^2+randn(size(y))*10;
% shadedErrorBar(x,y,{@mean,@std},'lineprops','-b','transparent',1);
% hold off
%
%
% Rob Campbell - November 2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parse input arguments
narginchk(3,inf)

params = inputParser;
params.CaseSensitive = false;
params.addParameter('lineProps', '-k', @(x) ischar(x) | iscell(x));
params.addParameter('transparent', true, @(x) islogical(x) || x==0 || x==1);
params.addParameter('patchSaturation', 0.2, @(x) isnumeric(x) && x>=0 && x<=1);

params.parse(varargin{:});

%Extract values from the inputParser
lineProps =  params.Results.lineProps;
transparent =  params.Results.transparent;
patchSaturation = params.Results.patchSaturation;

if ~iscell(lineProps), lineProps={lineProps}; end

%Process y using function handles if needed to make the error bar dynamically
if iscell(errBar)
    fun1=errBar{1};
    fun2=errBar{2};
    errBar=fun2(y);
    y=fun1(y);
else
    y=y(:).';
end

if isempty(x)
    x=1:length(y);
else
    x=x(:).';
end

%Make upper and lower error bars if only one was specified
if length(errBar)==length(errBar(:))
    errBar=repmat(errBar(:)',2,1);
else
    s=size(errBar);
    f=find(s==2);
    if isempty(f), error('errBar has the wrong size'), end
    if f==2, errBar=errBar'; end
end

if length(x) ~= length(errBar)
    error('length(x) must equal length(errBar)')
end

%Log the hold status so we don't change
initialHoldStatus=ishold;
if ~initialHoldStatus, hold on,  end

H = makePlot(x,y,errBar,lineProps,transparent,patchSaturation);

if ~initialHoldStatus, hold off, end

if nargout==1
    varargout{1}=H;
end

function H = makePlot(x,y,errBar,lineProps,transparent,patchSaturation)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Plot to get the parameters of the line
    H.mainLine=plot(x,y,lineProps{:});

    % Work out the color of the shaded region and associated lines.
    % Here we have the option of choosing alpha or a de-saturated
    % solid colour for the patch surface.
    mainLineColor=get(H.mainLine,'color');
    edgeColor=mainLineColor+(1-mainLineColor)*0.55;

    if transparent
        faceAlpha=patchSaturation;
        patchColor=mainLineColor;
    else
        faceAlpha=1;
        patchColor=mainLineColor+(1-mainLineColor)*(1-patchSaturation);
    end

    %Calculate the error bars
    uE=y-errBar(1,:);
    lE=y+errBar(2,:);

    %Add the patch error bar

    %Make the patch
    yP=[lE,fliplr(uE)];
    xP=[x,fliplr(x)];

    %remove nans otherwise patch won't work
    xP(isnan(yP))=[];
    yP(isnan(yP))=[];

    H.patch=patch(xP,yP,1,'facecolor',patchColor, ...
                  'edgecolor','none', ...
                  'facealpha',faceAlpha);

    %Make pretty edges around the patch.
    H.edge(1)=plot(x,lE,'-','color',edgeColor);
    H.edge(2)=plot(x,uE,'-','color',edgeColor);

    uistack(H.mainLine,'top') % Bring the main line to the top

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function [timeson, timesoff, fbin]=ontimestat(bm, threshold)
% [timeson, timesoff, fbin]=ontimestat(bm, thershold)
% Computes ON and OFF times of the time serires bm. ON (OFF) time is how long is
% the time-serires bm above (below) the threshold.
% bm=f';
% threshold=500;

fbin=(bm>threshold);
timeson=computetimes(fbin);
timesoff=computetimes(1-fbin);
end
function times=computetimes(fbin)
zv = zeros(1,size(fbin,2));
fbin2=[zv; fbin; zv];
pos=find(diff(fbin2)>0);
neg=find(diff(fbin2)<0);
times=neg-pos;
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function rts_propUpdateFunc(~, ~, hFig, tcpipClientRes, tcpipClientMainEng, tcpipClientRcsEng, hResUITable1, hResSumProgBar1, hResUITable2, hResSumProgBar2, hResourceCombo1, hResourceCombo2, hMainEngineUITable, hRcsEngineUITable)
%rts_propUpdateFunc Summary of this function goes here
%   Detailed explanation goes here

    try
        data = get(tcpipClientRes,'UserData');
        resources = getResourcesStruct(data{1}, 3);

        if(~isempty(resources))
            %UPDATE RESOURSE COMBOBOX 1
            if(~isempty(hResourceCombo1))
                updateResComboKeepSelection(hResourceCombo1, resources);
            end

            %UPDATE RESOURSE COMBOBOX 2
            if(~isempty(hResourceCombo2))
                updateResComboKeepSelection(hResourceCombo2, resources);
            end

            %UPDATE RESOURCE TABLE 1
            resource = getResourceNameFromCombo(hResourceCombo1);
            if(~isempty(hResUITable1) && ~isempty(hResSumProgBar1) && isfield(resources, resource))
                generateResourceStorageTable(hResUITable1, resources.(resource));
                [progBar1Str,percent] = generateResourceSummaryStr(resources, resource);
                hResSumProgBar1.setValue(percent);
                hResSumProgBar1.setString(progBar1Str);
                hResSumProgBar1.setStringPainted(true);
            end

            %UPDATE RESOURCE TABLE 2
            resource = getResourceNameFromCombo(hResourceCombo2);
            if(~isempty(hResUITable2) && ~isempty(hResSumProgBar2) && isfield(resources, resource))
                generateResourceStorageTable(hResUITable2, resources.(resource));
                [progBar2Str,percent] = generateResourceSummaryStr(resources, resource);
                hResSumProgBar2.setValue(percent);
                hResSumProgBar2.setString(progBar2Str);
                hResSumProgBar2.setStringPainted(true);
            end
        end

        %UPDATE MAIN ENGINE TABLE
        data = get(tcpipClientMainEng,'UserData');
        r = readPartStrAndDblsFromDblArr(data{1}, 6);
        if(~isempty(r))
            if(~isempty(hMainEngineUITable))
                generateMainEngineStatusTable(hMainEngineUITable, r);
            end
        end

        %UPDATE RCS ENGINE TABLE
        data = get(tcpipClientRcsEng,'UserData');
        r = readPartStrAndDblsFromDblArr(data{1}, 3);
        if(~isempty(r))
            if(~isempty(hRcsEngineUITable))
                generateRCSEngineStatusTable(hRcsEngineUITable, r);
            end
        end

    catch ME
%         disp('----');
%         disp(ME.message);
%         disp(ME.stack(1));
%         disp(ME.stack(2));
%         disp('----');
    end

end

function generateRCSEngineStatusTable(hRcsEngineUITable, tData)
    tableData = {};
    for(i=1:size(tData,1))
        row = tData(i,:);
        split = strsplit(row{1}, '/_/_/');

        tableData{i,1} = split{1};
        tableData{i,2} = row{2};
        tableData{i,3} = row{3};
        tableData{i,4} = row{4};
    end

    setTableData(hRcsEngineUITable, tableData);
end

function generateMainEngineStatusTable(hMainEngineUITable, tData)
    tableData = {};
    for(i=1:size(tData,1))
        row = tData(i,:);
        split = strsplit(row{1}, '/_/_/');

        tableData{i,1} = split{1};
        tableData{i,2} = row{2}*100;
        tableData{i,3} = row{3};
        tableData{i,4} = row{4};
        tableData{i,5} = row{5};
    end

    setTableData(hMainEngineUITable, tableData);
end

function str = getResourceNameFromCombo(hResourceCombo)
    contents = cellstr(get(hResourceCombo,'String'));
    str = contents{get(hResourceCombo,'Value')};
end

function updateResComboKeepSelection(hResourceCombo, resources)
    contents = cellstr(get(hResourceCombo,'String'));
    if(get(hResourceCombo,'Value') <= length(contents))
        sel = contents{get(hResourceCombo,'Value')};
    else
        sel = '';
    end

    res = sort(fieldnames(resources));
    if(~isempty(res))
        set(hResourceCombo,'String', res);
        ind=find(ismember(res,sel));
    else
        return;
    end

    if(isempty(ind) || ind<=0)
        value = 1;
    else
        value = ind;
    end

    set(hResourceCombo,'Value', value);
end

function [str,percent] = generateResourceSummaryStr(resources, resource)
    [current, max] = getCurrentMaxResource(resources, resource);

    percent = uint32(round(100*(current/max)));
    resStr = cap1stLetter(resource);

    strT = ['Total ', resStr ,' = ',num2str(current),' / ',num2str(max)];
    str = java.lang.String(strT);
end

function generateResourceStorageTable(hResUITable, resArr)
    tableData = {};
    for(i=1:size(resArr,1))
        if(resArr{i,6} == 1)
            connected = 'True';
        else
            connected = 'False';
        end

        tableData{i,1} = resArr{i,2};
        tableData{i,2} = connected;
        tableData{i,3} = resArr{i,4};
        tableData{i,4} = resArr{i,5};
    end

    setTableData(hResUITable, tableData);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text: path(path,'/afs/eos/info/ie/ie589k_info/GAOT');

author: Grootzz

repo: GA-BP

save_path: github-repos/MATLAB/Grootzz-GA-BP


text:
function block_size = decide_max_block_size_RBM(para)

if para.useGPU == 1
    gpu = gpuDevice;
    GPU_max_memory =  gpu.FreeMemory/8; % This is the maximum number of double precision numbers that can be stored in the GPU
    % Compute how much training data can we put into GPU
    inputDim = para.layerSize(1);
    outputDim = para.layerSize(end);
    max_sample = GPU_max_memory / (inputDim + outputDim);
    block_size = max_sample / 5;   % Let's be conservative, don't use all the memory of GPU
else
    block_size = 1e5;
end

author: singaxiong

repo: SignalGraph


text:
function n = planeNormal(plane)
%PLANENORMAL Compute the normal to a plane
%
%   N = planeNormal(PLANE)
%   compute the normal of the given plane
%   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
%   N : [dx dy dz]
%
%   See also
%   geom3d, planes3d, createPlane
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/02/2005.
%

%   HISTORY
%   15/04/2013 Extended to N-dim planes by Sven Holcombe

% plane normal
outSz = size(plane);
outSz(2) = 3;
n = zeros(outSz);
n(:) = crossProduct3d(plane(:,4:6,:), plane(:, 7:9,:));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [vMc, fMc_org, fStd_Mc, fMeanMc, fMedianMc] = calc_BstMcMaxCurvVar(mCatalog,fBinning, nSample)
% function [vMc, fMc_org, fStd_Mc, fMeanMc, fMedianMc] = calc_BstMcMaxCurvVar(mCatalog,fBinning, nSample)
%--------------------------------------------------------------------------------------------------------
% Bootstrap EQ catalog and determine Mc using maximum curvature
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Magnitude binning interval
% nSample    : Number of bootstrap samples
%
% Outgoing variables:
% vMc      : Best Mc estimate according to MLS
% vMls     : Vector of maximum likelihood scores
% fMc_org  : Mc estimate from original data
% fStd_Mc  : Standard deviation (assuming normal distribution)
% fSkew    : Skewness of Mc distribution
% vPerc    : Percentiles at [5 10 90 95] percent levels
% fMedianMc  : Median (50 percentile) of vMc
% fMeanMc  : Mean of Mc
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 14.02.03

% Check input
if nargin == 0, error('No catalog input'); end
if nargin == 1, fBinning = 0.1, nSample = 500, disp('Default Bin size 0.1, Default 500 bootstrap samples');end
if nargin == 2, nSample = 500, disp('Default 500 bootstrap samples'); end
if nargin > 3, error('Too many arguments!'); end

% Initialize
vMls = [];
vMc = [];

% Get magnitudes
vMags = mCatalog(:,6);

% Create bootstrap samples using bootstrap matlab toolbox
mMag_bstsamp = bootrsp(vMags,nSample);
% First step:
fMc_org = calc_Mc(mCatalog, 1);
% Determine Mc uncertainty
for nSamp=1:nSample
    mCatalog(:,6) = mMag_bstsamp(:,nSamp);
    fMc = calc_Mc(mCatalog, 1);
    vMc =  [vMc; fMc];
end

% Check for Nan and create output
vSel = isnan(vMc);
vNoNanMc = vMc(~vSel,:);
if ~isempty(vNoNanMc)
    fStd_Mc = std(vNoNanMc);
    fMeanMc = mean(vNoNanMc);
    fMedianMc = median(vNoNanMc);
else
    fStd_Mc = nan;
    fMeanMc = nan;
    fMedianMc = nan;
end

author: CelsoReyes


text:
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function status=mrInitRetFromFSL_NIC(procScans,SessionName,subjectName)
%   status=mrInitRetFromFSL(foundScans,SessionName,subjectName)
%   creates a MLR session out of a mixture of Dicom and analyze files
%   this mixture is the result of the prepareUCSF script that preprocesses
%   the functional data from the UCSF scanner, but leaves the inplanes
%   as Dicom.
%   2006.04.20 Mark Schira (mark@ski.org) wrote it.
% 061407: ARW Modified it to use the new pipeline for processing data from
% the NIC.
% Found scans will contain a sparse list of processed scans that can be
% found in ./Reconned.
% The functional data will be in _4d.hdr files. The inplanes will be in a
% _3d file. There may be other data sets there (T1s etc taken during the
% same session).
% First we will identify the inplane scans. Then we will do the
% functionals.

mrGlobals
HOMEDIREC = pwd;

doCrop=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the inplane anatomies
% We can read the volume siye from the header. Maybe other info as well.
% The inplane data should be in the dataPath named anat.hdr/img

InplaneDir=0;
% find the InplaneDir
for thisScan=1:length(procScans)
    if (~isempty(procScans{thisScan}))
    if strcmp(procScans{thisScan}.Action,'Inplanes')
        InplaneDir=thisScan; %take the last if several are set
    end
    end

end
%get the anatomy
%anatDir=[HOMEDIREC,filesep,procScans{InplaneDir}.DirName];
if (~InplaneDir)
    % Prompt for an analyze file that serves as an inplane ref
    [fname,pathname]=uigetfile('*.img','Pick Analyze/Nifti format 3d file');
	 if isnumeric(fname)
		 % try again?
	    [fname,pathname]=uigetfile('*.img','Pick Analyze/Nifti format 3d file');
		 if isnumeric(fname)
			 return
% 			 keyboard		% sometimes dialog never shows up??? fname=pathname=0 then GetAnatomyFromAnalyze crashes
		 end
	 end
    fname=fullfile(pathname,fname);
else
    fname=[procScans{InplaneDir}.Filename3d,'.img'];
end

[anat, inplanes] = GetAnatomyFromAnalyze(fname);
% %since cropping is anoying
inplanes.cropSize=inplanes.fullSize;
inplanes.crop=cat(1,[1 1],inplanes.fullSize);

%AnatomySize=size(anat);
% Save anat
anatFile = fullfile(HOMEDIREC, 'Inplane', 'anat.mat');
save(anatFile, 'anat', 'inplanes');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the functionals
functScans=[];
for thisScan=1:length(procScans)
    if (~isempty(procScans{thisScan}))
        if strcmp(procScans{thisScan}.Action,'Functional')
            functScans=cat(2,functScans,thisScan); %take the last if several are set

        end
    end

end

% put information about the functionals here. I tried to make this independent of the DICOM format but in the end I needed to get at least the TR from that header.
% I know that the individual analyze vols should have the TR as the 4th
% dimension. Something to work on...

timesRot90 = 1;
for thisScan=1:length(functScans)
    scanParams(thisScan).hdr=spm_vol([procScans{functScans(thisScan)}.Filename4d,'.img']);
    scanParams(thisScan).imatrix=spm_imatrix(scanParams(thisScan).hdr(1).mat); % imatrix decomposes the analyze (nifti) header into its components.
    scanParams(thisScan).fullSize=scanParams(thisScan).hdr(1).dim; %#ok<AGROW>
	 if rem(timesRot90,2)
		 scanParams(thisScan).fullSize = scanParams(thisScan).fullSize([2 1 3]);
     end
%      if (ieNotDefined('procScans{functScans(thisScan)}.SkipVols'))
%          disp('Warning - we did not find a .SkipVols field in procScans');
%
%          procScans{functScans(thisScan)}.SkipVols=0;
%      else
          fprintf('\nSkipping %d vols\n',procScans{functScans(thisScan)}.SkipVols);
%      end

    scanParams(thisScan).junkFirstFrames=procScans{functScans(thisScan)}.SkipVols;
    scanParams(thisScan).nFrames=procScans{functScans(thisScan)}.Volumes-procScans{functScans(thisScan)}.SkipVols;
% We don't crop these days. Perhaps we could get rid of that crop param
% altogether?
    scanParams(thisScan).cropSize=scanParams(thisScan).fullSize(1:2);
    scanParams(thisScan).crop=cat(1,[1 1],scanParams(1).fullSize(1:2));
    scanParams(thisScan).TR=procScans{functScans(thisScan)}.hdr{1}.RepetitionTime;
    scanParams(thisScan).PfileName=procScans{functScans(thisScan)}.hdr{1}.SeriesNumber;
    scanParams(thisScan).framePeriod=scanParams(thisScan).TR/1000; % this is in seconds
    scanParams(thisScan).slices=1:procScans{functScans(thisScan)}.Slices;
	 scanParams(thisScan).effectiveResolution = abs(scanParams(thisScan).imatrix(7:9));		% voxel dimensions (mm)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create/load mrSESSION and dataTYPES, modify them, and save

% If mrSESSION already exits, load it.
sessionFile = fullfile(HOMEDIREC, 'mrSESSION.mat');
if exist(sessionFile, 'file')
    loadSession;
    % if docrop, make sure that the mrSESSION is up-to-date
    if doCrop
        mrSESSION.inplanes = inplanes;
        mrSESSION = UpdateSessionFunctionals(mrSESSION,scanParams);
        saveSession;
    end
end

% If we don't yet have a session structure, make a new one.
if (~exist('mrSESSION','var'))
    mrSESSION = CreateNewSession(HOMEDIREC, inplanes, mrLoadRetVERSION);
end

if isempty(mrSESSION)
    mrSESSION = CreateNewSession(HOMEDIREC, inplanes, mrLoadRetVERSION);
end

% Update mrSESSION.functionals with scanParams corresponding to any new Pfiles.
% Set mrSESSION.functionals(:).crop & cropSize fields
mrSESSION = UpdateSessionFunctionals(mrSESSION,scanParams,0);

% this lines replace then EditSession command. One gui-window less to bother
mrSESSION.description=SessionName;
mrSESSION.subject=subjectName;

saveSession;

% Create/edit dataTYPES
if ieNotDefined('dataTYPES')
    dataTYPES = CreateNewDataTypes(mrSESSION);
else
    dataTYPES = UpdateDataTypes(dataTYPES,mrSESSION);
end

% Save any changes that may have been made to mrSESSION & dataTYPES
saveSession;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extract time series
mrVista('inplane');
v=getSelectedInplane;

for thisScan=1:length(functScans)
    analyzeNum=functScans(thisScan);
    inFile=procScans{analyzeNum}.Filename4d; % This has a trailing .hdr part
    %make sure there is a .hdr at the end...
    [pathstr, name, ext, versn]=fileparts(inFile);
    inFile=[pathstr,filesep,name,'.hdr'];
    disp(inFile)
    volsToSkip=mrSESSION.functionals(thisScan).junkFirstFrames;
    scaleFact=[1];
    v=read4dAnalyzeToTseries(v,inFile,thisScan,volsToSkip,0,timesRot90,0); % No slice flip, 3 90 deg rotations
end

% When we only recon some of thedata, this will crash because functScans
% does not contain entries for the missing data set but we try to access
% them three lines lower down from here.
for thisScan=1:length(dataTYPES(1).blockedAnalysisParams)
    dataTYPES(1).blockedAnalysisParams(thisScan).temporalNormalization=0;
    dataTYPES(1).blockedAnalysisParams(thisScan).nCycles=procScans{functScans(thisScan)}.Cycles;
    dataTYPES(1).eventAnalysisParams(thisScan).eventAnalysis=0;
end
saveSession;

status=1;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef SymbolManager
    % SYMBOLMANAGER allows you to control appearance of symbols on map
    properties(Constant)
        %fields = {'LineWidth','LineStyle','Color',...
        %    'Marker','MarkerSize','MarkerEdgeColor','MarkerFaceColor'};
         fields={'DisplayName';'Color';'LineStyle';'LineWidth';'Marker';'MarkerSize';'MarkerEdgeColor';'MarkerFaceColor'}
    end
    properties
        markerNames ={'+','o','*','.','x','square','diamond','^','v','<','>','pentagram','hexagram','none'};
        markerValues={'+','o','*','.','x','s',     'd',      '^','v','<','>','p',        'h',       'none'};

        lineValues = {'-','--',':','-.','none'};

        colors = {'yellow','magenta','cyan','red','lime','blue','white','black','choose'};
        colors_rgb = {[1 1 0],[1 0 1],[0 1 1], [1 0 0], [0 1 0], [0 0 1], [1 1 1], [0 0 0]};

    end

    methods
        function obj = SymbolManager()

        end

        function C=fill_values(obj,h)

            switch h.Type
                case 'scatter'
                    [MECidx, MarkerEdgeColor] = Xcolor({'flat';'none'}, h.MarkerEdgeColor);
                    [MFCidx, MarkerFaceColor] = Xcolor({'auto','flat','none'}, h.MarkerFaceColor);
                    if size(h.CData,2) == 3
                        [Cidx, CColor] = Xcolor({}, h.CData(1,:));
                    else
                        CColor={'various','choose'};
                        Cidx=1;
                    end

                    mIndex=find(strcmp(h.Marker,obj.markerNames));

                    C = {... fieldname, fieldwidth, fieldtype,
                        'DisplayName',      120, {'Style', 'text','String',h.DisplayName};...
                        'CData',            130,  {'Style','PopupMenu','String',CColor,'Value',Cidx,...
                                              'Callback', @(s,~)setLineEdgeFaceColor(s,'CData')};...
                        'LineStyle',        60, {'Style','text','String',''};...
                        'LineWidth',        50,  {'Style','edit','String',num2str(h.LineWidth),...
                                              'Callback', @(src,~)set(h,'LineWidth',str2double(src.String))};...
                        'Marker',           90,  {'Style','PopupMenu','String',obj.markerNames,'Value',find(strcmp(h.Marker,obj.markerNames)),...
                                              'Callback', @(src,~)set(h,'Marker',obj.markerNames{src.Value})};...
                        'SizeData',         50,  {'Style','edit','String',num2str(mean(h.SizeData)),...
                                              'Callback', @(src,~)set(h,'SizeData',str2double(src.String))};...
                        'MarkerEdgeColor',  130,  {'Style','PopupMenu','String',MarkerEdgeColor,'Value',MECidx,...
                                              'Callback', @(s,~)setLineEdgeFaceColor(s,'MarkerEdgeColor')};...
                        'MarkerFaceColor',  130,  {'Style','PopupMenu','String',MarkerFaceColor,'Value',MFCidx,...
                                              'Callback', @(s,~)setLineEdgeFaceColor(s,'MarkerFaceColor')};...
                        };

                case 'line'
                    [MECidx, MarkerEdgeColor] = Xcolor({'auto';'none'}, h.MarkerEdgeColor);
                    [MFCidx, MarkerFaceColor] = Xcolor({'auto','none'}, h.MarkerFaceColor);
                    [Cidx, CColor] = Xcolor({}, h.Color(1,:));

                    mIndex=find(strcmp(h.Marker,obj.markerNames));
                    assert(~isempty(mIndex),'marker issue %s not in markerNames',h.Marker);

                    C = {... fieldname, fieldwidth, fieldtype, inLine, inScatter, getfun, setfun
                        'DisplayName',  120, {'Style','text','String',h.DisplayName};...
                        'Color',        130,  {'Style','PopupMenu','String',CColor,'Value',Cidx,...
                                              'Callback',@(s,~)setLineEdgeFaceColor(s,'Color')};...
                        'LineStyle',    60,  {'Style','PopupMenu','String',obj.lineValues,'Value',find(strcmp(h.LineStyle,obj.lineValues)),...
                                              'Callback',@(src,~)set(h,'LineStyle',obj.lineValues{src.Value})};...
                        'LineWidth',    50,  {'Style','edit','String',num2str(h.LineWidth),...
                                              'Callback',@(src,~)set(h,'LineWidth',str2double(src.String))};...
                        'Marker',       90,  {'Style','PopupMenu','String',obj.markerNames,'Value',mIndex,...
                                              'Callback',@(src,~)set(h,'Marker',obj.markerNames{src.Value})};...
                        'MarkerSize',   50,  {'Style','edit','String',num2str(h.MarkerSize),...
                                              'Callback',@(src,~)set(h,'MarkerSize',str2double(src.String))};...
                        'MarkerEdgeColor',130,  {'Style','PopupMenu','String',MarkerEdgeColor,'Value',MECidx,...
                                              'Callback',@(s,~)setLineEdgeFaceColor(s,'MarkerEdgeColor')};...
                        'MarkerFaceColor',130,  {'Style','PopupMenu','String',MarkerFaceColor,'Value',MFCidx,...
                                              'Callback', @(s,~)setLineEdgeFaceColor(s,'MarkerFaceColor')};...
                        };
                otherwise
                    C={};
            end
            for i=1:size(C,1)
                C{i,4} = @()h.(C{i,1}); % get function
                C{i,5} = @(val)set(h,C{i,1},val); % set function
            end
            %disp(C)

            function setLineEdgeFaceColor(src,fld)
                if ismember(src.String{src.Value},{'auto','none','flat'})
                    h.(fld)=src.String{src.Value};
                elseif src.Value==numel(src.String) % last option is always CHOOSE
                    try
                        newc=uisetcolor(h.(fld));
                    catch
                        newc=uisetcolor([0 0 0]);
                    end
                    h.(fld)=newc;
                    % add this "new" color to the list
                    src.String(end+1)=src.String(end);
                    nm = FancyColors.name(newc);
                    src.String(end-1)={FancyColors.colorize(nm,nm)};
                else
                    cs = src.String;
                    fancy=endsWith(cs,'');
                    cs(fancy)=extractBetween(cs(fancy),'">',"");

                    nAdded = find(strcmp(cs,obj.colors{1})) - 1; % zero if nothing addeds

                    %nAdded = find(strcmp(src.String,obj.colors{1})) - 1; % zero if nothing added

                    idx=src.Value - nAdded;
                    idx=idx(1);
                    newc = FancyColors.rgb(cs{src.Value});
                    %{
                    if idx>numel(obj.colors_rgb)
                        %newc = FancyColors.rgb(src.String{src.Value});
                    else
                        newc = obj.colors_rgb{idx};
                    end
                    %}
                    h.(fld)= newc;
                end
            end

            function [idx, displist]=Xcolor(additionalOptions, existingValue)
                % find position in list. if it doesn't exist, add it. return index.

                % add field-unique options to top of list
                %displist = [additionalOptions(:) ; obj.colors(:)];
                cl = obj.colors(:);
                for n=1:numel(cl)
                    if ~startsWith(cl(n),'')
                        cl(n)={ FancyColors.colorize(cl{n},cl{n}) };
                    end
                end
                displist = [additionalOptions(:) ; cl];
                lookuplist = [additionalOptions(:) ; obj.colors_rgb(:)];

                % find existing value's position within the lookup list
                idx = find(cellfun(@(x)isequal(x,existingValue), lookuplist));

                % if our value isn't found, then add it to the top of the list, and set index to it.
                if isempty(idx) || idx==0
                    nm = FancyColors.name(existingValue);
                    displist(end+1)=displist(end);
                    displist(end-1)={FancyColors.colorize(nm,nm)};
                    idx=numel(displist)-1;
                end
            end
        end

        function test(obj, ax)

            % show the legend. has dual purpose. it assigns DisplayName and provides feedback for user
            wasVisible=~isempty(ax.Legend) && ax.Legend.Visible == "on";
            if ~wasVisible
                legend(ax,'show');
                ax.Legend.Visible='on';
            end

            f=findobj('Tag','testfigure');

            ch=ax.Children;
            x=20; xspace=5;
            y = 20; h=20;
            fullht = y+10;

            ts = ax.Title.String;
            if isempty(ts)
                ts=ax.Tag;
            end

            if isempty(f)
                f=figure('Tag','testfigure','Name',sprintf('Adjust Symbols for %s',ts));
                f.Position(3)=760;
                f.Position(4)=numel(ch)*fullht+60;
                f.MenuBar='none';
            else
                set(f,'Name',sprintf('Adjust Symbols for %s',ts));
                clf
            end

            xoff=x;
            for itt=numel(ch): -1 : 1
                if ~isvalid(ch(itt))
                    ch(itt)=[];
                    continue;
                end
                if (~strcmp(ch(itt).Type,'line') && ~strcmp(ch(itt).Type,'scatter'))
                    continue
                end
                C=obj.fill_values(ch(itt));
                for i=1:size(C,1)
                    w=C{i,2};
                    misc_nv = C{i,3};
                    uicontrol(gcf,'Position',[xoff,y,w,20],misc_nv{:});
                    xoff = xoff + w + xspace;
                end

                y=y+30; % up to next row
                if itt==1 && f.Position(3) < xoff+5
                    f.Position(3)= xoff+5;
                end
                xoff=x;
            end

            % label the fields
            flds = {'DisplayName';'Color';'LineStyle';'LineWidth';'Marker';'MarkerSize';'MarkerEdgeColor';'MarkerFaceColor'};
            for n=1:numel(obj.fields)
                fn=obj.fields{n};
                w = C{n,2};
                uicontrol(f,'Style','text','String',fn,'HorizontalAlignment','center',...
                    'Position',[xoff, y, w, h]);
                xoff = xoff + w + xspace;
            end

            % hide it again.
            if ~wasVisible
                ax.Legend.Visible='off';
            end
        end

    end
    methods(Static)
        function cb(src,ev,ax)
            sm=SymbolManager;
            sm.test(ax);
        end
    end
end


text:
%% --------------------------
% MemNet_M6R6 for image denoising (gaussian noise)
% edit by yingtai 12/08/2017
% -------------------------------
function test_MemNet_M6R6_GD()
setenv('LC_ALL','C')
addpath /data2/taiying/MSU_Code/119-caffe-matlab/matlab; % change to your caffe path
setenv('GLOG_minloglevel','2')
addpath('../');
addpath('../evaluation_func/');
addpath('../evaluation_func/matlabPyrTools-master/');

%% parameters
gpu_id = 5;
Noise_level = 30;
data_set_id = 1;
thresh_hei = 120; % threshold patch size for inference, since too big image may cost too much memory
thresh_wid = 120;
rf = 16;

pathfolder = ['../../data/GaussianDenoising/'];
if data_set_id == 1
    % S14
    setTestCur = 'S14';
    path = [pathfolder setTestCur '/'];
    d = dir([path '*.bmp']);
    filenum = 14;
end
if data_set_id == 2
    % BSD200
    setTestCur = 'BSD200';
    path = [pathfolder setTestCur '/'];
    d = dir([path '*.jpg']);
    filenum = 200;
end

savepath = ['./results/'];
folderResultCur = fullfile(savepath, [setTestCur,'_Noise',num2str(Noise_level)]);
%%% folder to store results
if ~exist(folderResultCur,'file')
    mkdir(folderResultCur);
end

mean_noised = [];
mean_memnet = [];
% caffe.set_mode_cpu(); % for CPU
caffe.set_mode_gpu(); % for GPU
caffe.set_device(gpu_id);

weights = ['../../model/MemNet_M6R6_80C64_GD.caffemodel'];
model_path = './MemNet_M6R6_80C64_deploy';

noise_set =[];
memnet_set = [];
im_b_set = cell(filenum,1);
im_h_set = cell(filenum,1);
im_gnd_set = cell(filenum,1);

for iii = 1:1:length(d)
    disp(['id: ' num2str(iii)]);
    imageName = d(iii).name;
    imageName = imageName(1:end-4);
    im  = imread([path d(iii).name]);
    randn('seed',54);

    %% test: resize
    if data_set_id == 3
        im = imresize(im,1/2);
    end

    %% rgb -> gray
    im_gray = im;
    if size(im,3)>1
        im_gray = rgb2gray(im);
    end
    im_gnd = im2double(im_gray);

    [hei,wid,channels] = size(im_gnd);
    % add noise
    im_b = single(im_gnd + Noise_level/255*randn(size(im_gray)));

    %% adaptively spilt
    % decide patch numbers
    hei_patch = ceil(hei/(thresh_hei+rf));
    wid_patch = ceil(wid/(thresh_wid+rf));
    hei_stride = ceil(hei/hei_patch);
    wid_stride = ceil(wid/wid_patch);
    use_start_x = 0;
    use_start_y = 0;
    use_end_x = 0;
    use_end_y = 0;

    ext_start_x = 0;
    ext_end_x = 0;
    ext_start_y = 0;
    ext_end_y = 0;

    posext_start_x = 0;
    posext_start_y = 0;
    posext_end_x = 0;
    posext_end_y = 0;

    % extract each patch for inference
    im_h = [];
    for x = 1 : hei_stride : hei
        for y = 1 : wid_stride : wid
            % decide the length of hei and wid for each patch
            use_start_x = x;
            use_start_y = y;
            if x - rf > 1 % add border
                ext_start_x = x-rf;
                posext_start_x = rf+1;
            else
                ext_start_x = x;
                posext_start_x = 1;
            end
            if y-rf > 1
                ext_start_y = y-rf;
                posext_start_y = rf+1;
            else
                ext_start_y = y;
                posext_start_y = 1;
            end

            use_end_x = use_start_x+hei_stride-1;
            use_end_y = use_start_y+wid_stride-1;

            if use_start_x+hei_stride+rf-1 <= hei
                hei_length = hei_stride+rf;
                ext_end_x = use_start_x+hei_length-1;
                posext_end_x = hei_length-rf+posext_start_x-1;

            else
                hei_length = hei-ext_start_x+1;
                ext_end_x = ext_start_x+hei_length-1;
                posext_end_x = hei_length;
                use_end_x = ext_start_x+hei_length-1;
            end
            if use_start_y+wid_stride+rf-1 <= wid
                wid_length = wid_stride+rf;
                ext_end_y = use_start_y+wid_length-1;
                posext_end_y = wid_length-rf+posext_start_y-1;

            else
                wid_length = wid-ext_start_y+1;
                ext_end_y = ext_start_y+wid_length-1;
                posext_end_y = wid_length;
                use_end_y = ext_start_y+wid_length-1;
            end

            subim_input = im_b(ext_start_x : ext_end_x, ext_start_y : ext_end_y);  % input
            data = permute(subim_input,[2, 1, 3]);
            model = [model_path '.prototxt'];
            subim_output = do_cnn(model,weights,data);
            subim_output = subim_output';
            subim_output = subim_output(posext_start_x:posext_end_x,posext_start_y:posext_end_y);
            % fill im_h with sub_output
            im_h(use_start_x:use_end_x,use_start_y:use_end_y) = subim_output;
        end
    end

    im_h1 = single(im_h) * 255;
    im_gnd1 = single(im_gnd) * 255;
    im_b1 = single(im_b) * 255;

    im_b_set{iii} = im_b1;
    im_h_set{iii} = im_h1;
    im_gnd_set{iii} = im_gnd1;

    %% compute PSNR and SSIM and IFC
    noised(1) = compute_psnr(im_gnd1,im_b1);
    memnet(1) = compute_psnr(im_gnd1,im_h1);
    noised(2) = ssim_index(im_gnd1,im_b1);
    memnet(2) = ssim_index(im_gnd1,im_h1);

    noise_set = [noise_set; noised];
    memnet_set = [memnet_set; memnet];
    %% save images
    imwrite(uint8(im_h1),fullfile(folderResultCur,[imageName,'_Noise',num2str(Noise_level),'.png']));
end
mean_noised = [mean_noised; [mean(noise_set(:,1)) mean(noise_set(:,2))]];
mean_memnet = [mean_memnet; [mean(memnet_set(:,1)) mean(memnet_set(:,2))]];

%%% save PSNR and SSIM metrics
PSNR_set = memnet_set(:,1);
SSIM_set = memnet_set(:,2);
save(fullfile(folderResultCur,['PSNR_',setTestCur,'_Noise',num2str(Noise_level),'.mat']),['PSNR_set'])
save(fullfile(folderResultCur,['SSIM_',setTestCur,'_Noise',num2str(Noise_level),'.mat']),['SSIM_set'])

disp(['noise = ' num2str(mean_noised(1,:)) '---- MemNet = ' num2str(mean_memnet(1,:))]);

end

author: tyshiwo

repo: MemNet

save_path: github-repos/MATLAB/tyshiwo-MemNet


text:
% run_simMegaExTEShaped.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% This script is run simply by editing the input parameters and then
% clicking "Run".
%
% DESCRIPTION:
% This script simulates an ExTE-MEGA-PRESS experiment with fully shaped editing
% and refocusing pulses.  Phase cycling of both the editing and refocusing
% pulses is performed.  Simulations are run at various
% locations in space to account for the within-voxel spatial variation of
% the GABA signal.  Summation across phase cycles and spatial positions is
% performed.  As a result of the phase cycling and spatially resolved simulations,
% this code takes a long time to run.  Therefore, the MATLAB parallel computing
% toolbox (parfor loop) was used to accelerate the siumulations.  Accelration
% is currently performed in the direction of the slice selective pulse along
% the x-direction, but this can be changed.  Up to a factor of 12 acceleration
% can be achieved using this approach.  To enable the use of the MATLAB
% parallel computing toolbox, initialize the multiple worked nodes using
% "matlabpool size X" where "X" is the number of available processing
% nodes.  If the parallel processing toolbox is not available, then replace
% the "parfor" loop with a "for" loop.
%
% INPUTS:
% To run this script, edit the parameters below as desired and then click
% "run":
% refocWaveform     = name of refocusing pulse waveform.
% editWaveform      = name of editing pulse waveform.
% editOnFreq        = freqeucny of edit on pulse[ppm]
% editOffFreq       = frequency of edit off pulse[ppm]
% refTp             = duration of refocusing pulses[ms]
% editTp            = duration of editing pulses[ms]
% Bfield            = Magnetic field strength in [T]
% Npts              = number of spectral points
% sw                = spectral width [Hz]
% Bfield            = magnetic field strength [Tesla]
% lw                = linewidth of the output spectrum [Hz]
% thkX              = slice thickness of x refocusing pulse [cm]
% thkY              = slice thickness of y refocusing pulse [cm]
% x                 = vector of X positions to simulate [cm]
% y                 = vector of y positions to simulate [cm]
% taus              = vector of pulse sequence timings  [ms]
% spinSys           = spin system to simulate
% editPhCyc1        = vector of phase cycling steps for 1st editing pulse [degrees]
% editPhCyc2        = vector of phase cycling steps for 2nd editing pulse [degrees]
% refPhCyc1         = vector of phase cycling steps for 1st refocusing pulse [degrees]
% refPhCyc2         = vector of phase cycling steps for 2nd refocusing pulse [degrees]
%
% OUTPUTS:
% outON_posxy       = Simulated ExTE-MEGA-PRESS edit-ON spectrum, spatially resolved.
% outOFF_posxy      = Simulated ExTE-MEGA-PRESS edit-OFF spectrum, spatially resolved.
% outDIFF_posxy     = Simulated ExTE-MEGA-PRESS difference spectrum, spatially resolved.
% outON             = Simulated ExTE-MEGA-PRESS edit-ON spectrum, summed over
%                     all positions.
% outOFF            = Simulated ExTE-MEGA-PRESS edit-OFF spectrum, summed over
%                     all positions.
% outDIFF           = Simulated ExTE-MEGA-PRESS difference spectrum, summed over
%                     all positions.

% ************INPUT PARAMETERS**********************************
refocWaveform='sampleRefocPulse.pta'; %name of refocusing pulse waveform.
editWaveform='sampleEditPulse.pta'; %name of editing pulse waveform.
editOnFreq=1.88; %freqeucny of edit on pulse[ppm]
refTp=5.2; %duration of refocusing pulses[ms]
editTp=14; %duration of editing pulses[ms]
Npts=2048; %number of spectral points
sw=2000; %spectral width [Hz]
lw=2; %linewidth of the output spectrum [Hz]
Bfield=2.89; %Magnetic field strength in [T]
thkX=3.5; %slice thickness of x refocusing pulse [cm]
thkY=3.5; %slice thickness of y refocusing pulse [cm]
x=linspace(-2.0125,2.0125,12); %X positions to simulate [cm]
y=linspace(-2.0125,2.0125,12); %y positions to simulate [cm]
taus_inv=...  %Timing for J-Inverted scan
        [4.9,...    %time from excitation to 1st refoc pulse [ms]
        67.4885,... %time from 1st refoc pulse to 1st editing pulse [ms]
        33.5115,... %time from 1st editing pulse to 2nd refoc pulse [ms]
        33.4885,... %time from 2nd refoc pulse to 2nd editing pulse [ms]
        62.6115];   %time from 2nd editing pulse to ADC onset [ms]
taus_ref=... %Timing for J-refocused scan
        [4.9,...    %time from excitation to 1st refoc pulse [ms]
        50.4885,... %time from 1st refoc pulse to 1st editing pulse [ms]
        50.5115,... %time from 1st editing pulse to 2nd refoc pulse [ms]
        50.4885,... %time from 2nd refoc pulse to 2nd editing pulse [ms]
        45.6115];   %time from 2nd editing pulse to ADC onset [ms]
spinSys='GABA'; %spin system to simulate
centreFreq=3.0; %Centre frequency of MR spectrum [ppm]
editPhCyc1=[0 90]; %phase cycling steps for 1st editing pulse [degrees]
editPhCyc2=[0 90]; %phase cycling steps for 2nd editing pulse [degrees]
refPhCyc1=[0,90]; %phase cycling steps for 1st refocusing pulse [degrees]
refPhCyc2=[0,90]; %phase cycling steps for 2nd refocusing pulse [degrees]
% ************END OF INPUT PARAMETERS**********************************

%Load RF waveforms
refRF=io_loadRFwaveform(refocWaveform,'ref',0);
editRF=io_loadRFwaveform(editWaveform,'inv',0);

gamma=42577000; %gyromagnetic ratio

%Load spin systems
load spinSystems
if strcmp(spinSys,'MM');
    sys=sysGABA;
    sys.shifts(3)=1.7;
    sys.shifts(4)=1.7;
else
    sys=eval(['sys' spinSys]);
end

%Resample refocusing RF pulse from 400 pts to 100 pts to reduce
%computational workload
refRF=rf_resample(refRF,100);

%This is the step where the editing pulse waveform (initially a pulse with
%zero-frequency) is frequency shifted to produce and edit-on and an
%edit-off pulse;
editRFon=rf_freqshift(editRF,editTp,(centreFreq-editOnFreq)*Bfield*gamma/1e6);

Gx=(refRF.tbw/(refTp/1000))/(gamma*thkX/10000); %[G/cm]
Gy=(refRF.tbw/(refTp/1000))/(gamma*thkY/10000); %[G/cm]

[DX,DY]=meshgrid(x,y);

%n=1;
%totalIters=length(x)*length(y)*length(editPhCyc1)*length(editPhCyc2)*length(refPhCyc1)*length(refPhCyc2);

%Initialize structures:
outON_posxy_epc_rpc=cell(length(x),length(y),length(editPhCyc1),length(editPhCyc2),length(refPhCyc1),length(refPhCyc2));
outOFF_posxy_epc_rpc=cell(length(x),length(y),length(editPhCyc1),length(editPhCyc2),length(refPhCyc1),length(refPhCyc2));
outON_posxy_epc=cell(length(x),length(y),length(editPhCyc1),length(editPhCyc2));
outOFF_posxy_epc=cell(length(x),length(y),length(editPhCyc1),length(editPhCyc2));
outON_posxy=cell(length(x),length(y));
outOFF_posxy=cell(length(x),length(y));
outDIFF_posxy=cell(length(x),length(y));
outON=struct([]);
outOFF=struct([]);

%loop through space: Don't forget to initialize the parallel processing
%toolbox workers using 'matlabpool open N' (for N workers, 12 max).

for X=1:length(x);
%parfor X=1:length(x);
    for Y=1:length(y);
        for EP1=1:length(editPhCyc1)
            for EP2=1:length(editPhCyc2)
                for RP1=1:length(refPhCyc1)
                    for RP2=1:length(refPhCyc2)
                        disp(['Executing X-position ' num2str(X) ' of ' num2str(length(x)) ', '...
                            'Y-position ' num2str(Y) ' of ' num2str(length(y)) ', '...
                            'First Edit phase cycle ' num2str(EP1) ' of ' num2str(length(editPhCyc1)) ', '...
                            'Second Edit phase cycle ' num2str(EP2) ' of ' num2str(length(editPhCyc2)) ', '...
                            'First Refoc phase cycle ' num2str(RP1) ' of ' num2str(length(refPhCyc1)) ', '...
                            'Second Refoc phase cycle ' num2str(RP2) ' of ' num2str(length(refPhCyc2)) '!!!']);
                        outON_posxy_epc_rpc{X}{Y}{EP1}{EP2}{RP1}{RP2}=sim_megapress_shaped(Npts,sw,Bfield,lw,taus_ref,sys,...
                            editRFon,editTp,editPhCyc1(EP1),editPhCyc2(EP2),...
                            refRF,refTp,Gx,Gy,x(X),y(Y),refPhCyc1(RP1),refPhCyc2(RP2));
                        outOFF_posxy_epc_rpc{X}{Y}{EP1}{EP2}{RP1}{RP2}=sim_megapress_shaped(Npts,sw,Bfield,lw,taus_inv,sys,...
                            editRFon,editTp,editPhCyc1(EP1),editPhCyc2(EP2),...
                            refRF,refTp,Gx,Gy,x(X),y(Y),refPhCyc1(RP1),refPhCyc2(RP2));

                        if RP1==1 && RP2==1
                            outON_posxy_epc{X}{Y}{EP1}{EP2}=outON_posxy_epc_rpc{X}{Y}{EP1}{EP2}{RP1}{RP2};
                            outOFF_posxy_epc{X}{Y}{EP1}{EP2}=outOFF_posxy_epc_rpc{X}{Y}{EP1}{EP2}{RP1}{RP2};
                        else
                            outON_posxy_epc{X}{Y}{EP1}{EP2}=op_addScans(outON_posxy_epc{X}{Y}{EP1}{EP2},outON_posxy_epc_rpc{X}{Y}{EP1}{EP2}{RP1}{RP2},xor(RP1==length(refPhCyc1),RP2==length(refPhCyc2)));
                            outOFF_posxy_epc{X}{Y}{EP1}{EP2}=op_addScans(outOFF_posxy_epc{X}{Y}{EP1}{EP2},outOFF_posxy_epc_rpc{X}{Y}{EP1}{EP2}{RP1}{RP2},xor(RP1==length(refPhCyc1),RP2==length(refPhCyc2)));
                        end
                    end %end of 1st refocusing phase cycle loop
                end %end of 2nd refocusing phase cycle loop.

                if EP1==1 && EP2==1
                    outON_posxy{X}{Y}=outON_posxy_epc{X}{Y}{EP1}{EP2};
                    outOFF_posxy{X}{Y}=outOFF_posxy_epc{X}{Y}{EP1}{EP2};
                else
                    outON_posxy{X}{Y}=op_addScans(outON_posxy{X}{Y},outON_posxy_epc{X}{Y}{EP1}{EP2});
                    outOFF_posxy{X}{Y}=op_addScans(outOFF_posxy{X}{Y},outOFF_posxy_epc{X}{Y}{EP1}{EP2});
                end
                outDIFF_posxy{X}{Y}=op_subtractScans(outON_posxy{X}{Y},outOFF_posxy{X}{Y});
            end %end of 1st editing phase cycle loop.
        end %end of 2nd editing phase cycle loop.

        outON=op_addScans(outON,outON_posxy{X}{Y});
        outOFF=op_addScans(outOFF,outOFF_posxy{X}{Y});

    end %end of spatial loop (parfor) in y direction.
end %end of spatial loop (parfor) in x direction.

outDIFF=op_subtractScans(outON,outOFF);

figure
sim_make2DSimPlot(outON_posxy,2.75,3.25);
figure
sim_make2DSimPlot(outOFF_posxy,2.75,3.25);
figure
sim_make2DSimPlot(outDIFF_posxy,2.75,3.25);

author: CIC-methods

repo: FID-A

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


text:
function h = SideAxes(a,location,s,varargin)

%SideAxes - Add side axes to existing axes.
%
%  USAGE
%
%    h = SideAxes(a,location,s,<options>)
%
%    Using cell arrays will overlay variable pairs.
%
%    a              optional target axes (default = gca)
%    s              size, as proportion of target axes
%    location       'top', 'bottom', 'left' or 'right'
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'gap'         gap between target and side axes, as proportion of total
%                   width or height (default = 0.1)
%    =========================================================================
%
%  OUTPUT
%
%    h              handle to the new axes

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Defaults
gap = 0.1;

% Check number of parameters
if nargin < 2,
	error('Incorrect number of parameters (type ''help SideAxes'' for details).');
end

% Optional axes
if nargin == 2,
	s = location;
	location = a;
	a = gca;
elseif isstring_FMAT(a),
	varargin = {s,varargin{:}};
	s = location;
	location = a;
	a = gca;
end

% Check number of options
if mod(length(varargin),2) ~= 0,
	error('Incorrect number of parameters (type ''help SideAxes'' for details).');
end

% Check parameters
if ~ishandle(a),
	error('Incorrect axes (type ''help SideAxes'' for details).');
end
if ~isdscalar(s,'>0','<1'),
	error('Incorrect size (type ''help SideAxes'' for details).');
end
location = lower(location);
if ~isstring_FMAT(location,'top','bottom','left','right'),
	error('Incorrect location (type ''help SideAxes'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help SideAxes'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'gap',
			gap = varargin{i+1};
			if ~isdscalar(gap,'>=0','<1'),
				error('Incorrect value for property ''gap'' (type ''help SideAxes'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SideAxes'' for details).']);
	end
end

% Resize target axes to make room for side axes

main = a;
position = get(main,'position');
x = position(1);
y = position(2);
currentWidth = position(3); % current width
currentHeight = position(4); % current height

% The following variables will be used only when adding left or right side axes
horizontalGap = currentWidth * gap;
newWidth = currentWidth - horizontalGap;
newMainWidth = newWidth * (1-s);
newSideWidth = newWidth * s;

% The following variables will be used only when adding left or right side axes
verticalGap = currentHeight * gap;
newHeight = currentHeight - verticalGap;
newMainHeight = newHeight * (1-s);
newSideHeight = newHeight * s;

xLims = get(main,'xlim');
yLims = get(main,'ylim');

switch(location),
	case 'top',
		set(main,'position',[x y currentWidth newMainHeight]);
		h = axes('position',[x y+newMainHeight+verticalGap currentWidth newSideHeight]);
		xlim(h,xLims);
	case 'bottom',
		set(main,'position',[x y+newSideHeight+verticalGap currentWidth newMainHeight]);
		h = axes('position',[x y currentWidth newSideHeight]);
		xlim(h,xLims);
	case 'right',
		set(main,'position',[x y newMainWidth currentHeight]);
		h = axes('position',[x+newMainWidth+horizontalGap y newSideWidth currentHeight]);
		ylim(h,yLims);
	case 'left',
		set(main,'position',[x+newSideWidth+horizontalGap y newMainWidth currentHeight]);
		h = axes('position',[x y newSideWidth currentHeight]);
		ylim(h,yLims);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef OptimalControler < handle
    % Let s :=[x(0), .., x(n), u(0)....u(n-1)] be decision variables, then
    % optimal control problem is composed by
    % (i) objective function V(s) = s'*H*S in quadratic form
    % (ii) equality constraints specified by C_eq1 and C_eq2 where C_eq1*s =C_eq2
    % (iii) inequality constraints specified by C_ineq1 and C_ineq2 where C_ineq1 * s <= C_ineq2
    %
    % For potential future extensions, all constraints will be managed by ConstraintManager class.
    % Any time you add a new constraint, those constraint is pushed into the manager with a key (name of constraint).

    properties (SetAccess = private)
        sys; %system
        Xc; Uc; % constraints set for statespace and input space
        x_min; x_max; % lower and upper bound of Xc
        N; % prediction horizon
        Ak; % S.T.M of closed-roop system with LQR feedback
        n_opt; % dim. of optimization parameter s :=[x(0), .., x(n), u(0)....u(n-1)]
        H; % positive definite matrix for objective function V(s) = s'*H*s
        constraint_manager;
    end

    %% Public Methods
    methods (Access = public)

        function obj = OptimalControler(sys, Xc, Uc, N)
            obj.sys = sys;
            obj.Xc = Xc;
            obj.x_min = min(Xc.V, [], 1)';
            obj.x_max = max(Xc.V, [], 2)';
            obj.Uc = Uc;
            obj.N = N;

            obj.n_opt = obj.sys.nx*(obj.N+1)+obj.sys.nu*obj.N;
            obj.constraint_manager = ConstraintManager();

            obj.H = obj.construct_costfunction();
            [C_eq1, C_eq2] = obj.construct_dynamics_constraint();
            [C_ineq1, C_ineq2] = obj.construct_ineq_constraint(Xc, Uc);

            %% Let's change initial and dynamics constraints!!
            obj.constraint_manager.add_eq_constraint('dynamics', C_eq1, C_eq2);
            obj.constraint_manager.add_ineq_constraint('feasible', C_ineq1, C_ineq2);
        end

        function add_initial_eq_constraint(obj, x_init)
            % E * x0 = x_init
            idx_x0_start = 1;
            idx_x0_end = obj.sys.nx;

            C_eq1_init = zeros(obj.sys.nx, obj.n_opt);
            C_eq1_init(:, idx_x0_start:idx_x0_end) = eye(obj.sys.nx);
            C_eq2_init = x_init;
            obj.constraint_manager.add_eq_constraint('initial', C_eq1_init, C_eq2_init);
        end

        function add_terminal_constraint(obj, Xadd)
            [C_ineq1_add, C_ineq2_add] = add_ineq_constraint(obj, Xadd, obj.N+1);
            obj.constraint_manager.add_ineq_constraint('terinal', C_ineq1_add, C_ineq2_add);
        end

        function add_initial_constraint(obj, Xadd)
            [C_ineq1_add, C_ineq2_add] = add_ineq_constraint(obj, Xadd, 1);
            obj.constraint_manager.add_ineq_constraint('initial', C_ineq1_add, C_ineq2_add);
        end

        function [x_seq, u_seq] = solve(obj)
            quadprog_solved = 0;
            [C_eq1, C_eq2] = obj.constraint_manager.combine_all_eq_constraints();
            [C_ineq1, C_ineq2] = obj.constraint_manager.combine_all_ineq_constraints();

            options = optimoptions('quadprog', 'Display', 'none');
            [var_optim, ~, exitflag] = quadprog(obj.H, [], C_ineq1, C_ineq2, C_eq1, C_eq2, [], [], [], options);
            x_seq = reshape(var_optim(1:obj.sys.nx*(obj.N+1)), obj.sys.nx, obj.N+1);
            u_seq = reshape(var_optim(obj.sys.nx*(obj.N+1)+1:obj.n_opt), obj.sys.nu, obj.N);

        end

    end

    %% Methods Used in Constoructor
    methods (Access = private)

        function H = construct_costfunction(obj)
            % compute H
            Q_block = [];
            R_block = [];
            for itr=1:obj.N
                Q_block = blkdiag(Q_block, obj.sys.Q);
                R_block = blkdiag(R_block, obj.sys.R);
            end
            H = blkdiag(Q_block, obj.sys.P, R_block);
        end

        function [C_eq1, C_eq2] = construct_dynamics_constraint(obj)
            % compute C_eq1 and C_eq2
            function C_ss_eq1 = single_step_dynamics_eq1(k)
                % A x(k) - E x(k+1) + Bu(k) = 0
                idx_xk_start = obj.sys.nx * k + 1;
                idx_xk_end = obj.sys.nx * (k + 1);

                idx_xkp1_start = obj.sys.nx * (k + 1) + 1;
                idx_xkp1_end = obj.sys.nx * (k + 2);

                idx_uk_start  = obj.sys.nx * (obj.N+1) + obj.sys.nu * k + 1;
                idx_uk_end  = obj.sys.nx * (obj.N+1) + obj.sys.nu * (k + 1);

                C_ss_eq1 = zeros(obj.sys.nx, obj.n_opt);

                C_ss_eq1(:, idx_xk_start:idx_xk_end) = obj.sys.A;
                C_ss_eq1(:, idx_xkp1_start:idx_xkp1_end) = - eye(obj.sys.nx);
                C_ss_eq1(:, idx_uk_start:idx_uk_end) = obj.sys.B;
            end

            C_eq1 = [];
            for k = 0:obj.N-1
                C_ss_eq1 = single_step_dynamics_eq1(k);
                C_eq1 = [C_eq1; C_ss_eq1];
            end
            C_eq2 = zeros(size(C_eq1, 1), 1);
        end

        function [C_ineq1, C_ineq2] = construct_ineq_constraint(obj, Xc, Uc)
            % compute C_ineq
            [F, G, nc] = convert_Poly2Mat(Xc, Uc);

            F_block = [];
            G_block = [];
            for itr = 1:obj.N
                G_block = blkdiag(G_block, G);
            end
            for itr = 1:obj.N+1
                F_block = blkdiag(F_block, F);
            end
            C_ineq1 = [F_block, [G_block; zeros(nc, obj.sys.nu*obj.N)]];
            nc_total = size(C_ineq1, 1);
            C_ineq2 = ones(nc_total, 1);
        end

    end

    methods (Access = private)

        function [C_ineq1_add, C_ineq2_add] = add_ineq_constraint(obj, Xadd, k_add)
            % add a new constraint at time step k
            if Xadd.contains(zeros(2, 1)) % If Xadd contains the origin, the contraint can be expressed as C1*x<=1
                [F_add, ~, nc_add] = convert_Poly2Mat(Xadd, Polyhedron());
                C_ineq2_add = ones(nc_add, 1);

            else % in other cases, expressed in a general affine form C1*x<=C2
                F_add = Xadd.A;
                nc_add = size(F_add, 1);
                C_ineq2_add = Xadd.b;
            end

            C_ineq1_add = zeros(nc_add, obj.n_opt);
            C_ineq1_add(:, (k_add-1)*obj.sys.nx+1:k_add*obj.sys.nx) = F_add;
        end

    end
end

author: HiroIshida

repo: robust-tube-mpc


text:
function [ a_lu, pivot, info ] = r8ge_fa ( a, n )

%*****************************************************************************80
%
%% R8GE_FA factors a general matrix.
%
%  Discussion:
%
%    R8GE_FA is a simplified version of the LINPACK routine DGEFA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, real A(LDA,N), the matrix to be factored.
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Output, real A_LU(LDA,N), contains an upper
%    triangular matrix and the multipliers which were used to obtain
%    it.  The factorization can be written A = L * U, where L is a
%    product of permutation and unit lower triangular matrices and
%    U is upper triangular.
%
%    Output, integer PIVOT(N), a vector of pivot indices.
%
%    Output, integer INFO, singularity flag.
%    0, no singularity detected.
%    nonzero, the factorization failed on the INFO-th step.
%
  info = 0;
  a_lu(1:n,1:n) = a(1:n,1:n);

  for k = 1 : n-1
%
%  Find L, the index of the pivot row.
%
    l = k;
    for i = k+1 : n
      if ( abs ( a_lu(l,k) ) < abs ( a_lu(i,k) ) )
        l = i;
      end
    end

    pivot(k) = l;
%
%  If the pivot index is zero, the algorithm has failed.
%
    if ( a_lu(l,k) == 0.0E+00 )
      info = k;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8GE_FA - Warning!\n' );
      fprintf ( 1, '  Zero pivot on step %d\n', info );
      return
    end
%
%  Interchange rows L and K if necessary.
%
    if ( l ~= k )
      [ a_lu(l,k), a_lu(k,k) ] = r8_swap ( a_lu(l,k), a_lu(k,k) );
    end
%
%  Normalize the values that lie below the pivot entry A(K,K).
%
    a_lu(k+1:n,k) = -a_lu(k+1:n,k) / a_lu(k,k);
%
%  Row elimination with column indexing.
%
    for j = k+1 : n

      if ( l ~= k )
        [ a_lu(l,j), a_lu(k,j) ] = r8_swap ( a_lu(l,j), a_lu(k,j) );
      end

      a_lu(k+1:n,j) = a_lu(k+1:n,j) + a_lu(k+1:n,k) * a_lu(k,j);

    end

  end

  pivot(n) = n;

  if ( a_lu(n,n) == 0.0E+00 )
    info = n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8GE_FA - Warning!\n' );
    fprintf ( 1, '  Zero pivot on step %d\n', info );
    return
  end

  return
end


text:
function setMatlabPool(val, prop)
% change status of matlab pool
% val: open|close
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

if(nargin == 0)
    val = 'close';
end

if(prop.flagParallel)
    % parallel computation toolbox exists
    if(strcmp(val,'open'))
        if(~prop.openPool)
            if(matlabpool('size') > 0)
                warning('Closing an unopened MatLab pool');
                matlabpool close
            end
            matlabpool open
            prop.openPool = true;
        else
            warning('A currently used MatLab pool is already open');
        end

    elseif(strcmp(val,'close'))
        if(prop.openPool)
            matlabpool close
        else
            if(matlabpool('size') > 0)
                warning('Closing an unopened MatLab pool');
                matlabpool close
            end
        end
        prop.openPool = false;
    else
        error('setMatlabPool(): Undetermined value');
    end
end

end

author: thomaskuestner


text:
% RES = pyrBandIndices(INDICES, BAND_NUM)
%
% Return indices for accessing a subband from a pyramid
% (gaussian, laplacian, QMF/wavelet, steerable).

% Eero Simoncelli, 6/96.

function indices =  pyrBandIndices(pind,band)

if ((band > size(pind,1)) | (band < 1))
  error(sprintf('BAND_NUM must be between 1 and number of pyramid bands (%d).', ...
      size(pind,1)));
end

if (size(pind,2) ~= 2)
  error('INDICES must be an Nx2 matrix indicating the size of the pyramid subbands');
end

ind = 1;
for l=1:band-1
  ind = ind + prod(pind(l,:));
end

indices = ind:ind+prod(pind(band,:))-1;

repo: SelfExSR