text:
clear; clf;
		% Introducerea punctelor definite
t=[0,0.3,0.8,1.1,1.6,2.3]';
y=[0.5,0.82,1.14,1.25,1.35,1.4]';
		% Se reprezinta cu O-uri negre puntele date
plot(t,y,'Ko');
hold on;
        %Generarea matricei coeficientilor sistemului
A=[ones(size(t)),exp(-t),t.*exp(-t)];
        % Rezolvarea sistemului de ecuatii
a=A\y
        % Extragerea coeficientilor functiei de aproximare
a0=a(1); a1=a(2); a2=a(3);
        % Generarea punctelor in care se calculeaza
        % valorile functiei de aproximare
tt=linspace(0,2.5,200);
		% Calcularea valorilor functiei de aproximare
yy=a0+a1.*exp(-tt)+a2.*tt.*exp(-tt);
        % Reprezentarea grafica a functiei de aproximare
plot(tt,yy,'b','LineWidth', 1.5);
        % Personalizarea reprezentarii grafice
grid on;
xlabel('t');
ylabel('y');


text:
function [ landmarks ] = readLandmarks( filename )
%READLANDMARKS Summary of this function goes here
%   Detailed explanation goes here

fileID = fopen( filename );
%C = textscan(fileID,'%f %f','headerLines', 1);
C = textscan(fileID,'%f','CommentStyle','%');
fclose(fileID);
landmarks = cell2mat(C);

end

author: anilbas


text:
% ANFIS Feed-Forward operation.
function y = anfis_scatter_forward(x,mean,sigma,b,ThetaL4)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        								 							     				                                       %
%   			                 	NETWORK FUNCTIONALITY SECTION					   %
%      																					                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[NumInVars NumInTerms] = size(mean);

NumRules = NumInTerms;

% LAYER 1 - INPUT TERM NODES

 In2 = x*ones(1,NumInTerms);
 Out1 = 1./(1 + (abs((In2-mean)./sigma)).^(2*b));

% LAYER 2 - PRODUCT NODES
 Out2 = prod(Out1.',2);
 S_2 = sum(Out2);

 % LAYER 3 - NORMALIZATION NODES
  Out3 = Out2./S_2;

% LAYERSS 4 - 5: CONSEQUENT NODES - SUMMING NODE
 Aux1 = [x; 1]*Out3';

 a = reshape(Aux1,(NumInVars+1)*NumRules,1);  % New Input Training Data shaped as a column vector.

 y = a'*ThetaL4;

author: Sable

repo: mcbench-benchmarks


text:
Dataflow {
    TemporalMap(16,16) K;
    SpatialMap(Sz(R),1) Y;
    TemporalMap(Sz(S),1) X;
    TemporalMap(1,1) C;
    Cluster(16, P);
    SpatialMap(1,1) K;
    TemporalMap(Sz(R),1) Y;
    TemporalMap(Sz(S),1) X;
    TemporalMap(Sz(R),7) R;
    TemporalMap(Sz(S),7) S;
}

author: maestro-project

repo: maestro


text:
% undoradial    remove radial distortion
%
% [xl] = undoradial(x,K,kc)
%
% x ... 3xN coordinates of the distorted pixel points
% K ... 3x3 camera calibration matrix
% kc ... 4x1 vector of distortion parameters
%
% xl ... linearized pixel coordinates
%        these coordinates should obey the linear pinhole model
%
% It calls comp_distortion_oulu: undistort pixel coordinates.
% function taken from the CalTech camera calibration toolbox

function [xl] = undoradial(x_kk,K,kc)

cc(1) = K(1,3);
cc(2) = K(2,3);
fc(1) = K(1,1);
fc(2) = K(2,2);

% First: Subtract principal point, and divide by the focal length:
x_distort = [(x_kk(1,:) - cc(1))/fc(1);(x_kk(2,:) - cc(2))/fc(2)];

if norm(kc) ~= 0,
	% Third: Compensate for lens distortion:
	xn = comp_distortion_oulu(x_distort,kc);
else
   xn = x_distort;
end;

% back to the linear pixel coordinates
xl = K*[xn;ones(size(xn(1,:)))];

author: strawlab

repo: MultiCamSelfCal


text:
clear all, close all, clc

n = 1000;
q = n/4;
X = zeros(n,n);
X(q:(n/2)+q,q:(n/2)+q) = 1;

subplot(2,2,1), imshow(X); colormap gray, axis off

Y = imrotate(X,10,'bicubic');  % rotate 10 degrees
Y = Y - Y(1,1);
nY = size(Y,1);
startind = floor((nY-n)/2);
Xrot = Y(startind:startind+n-1, startind:startind+n-1);
subplot(2,2,2), imshow(Xrot); colormap gray, axis off

[U,S,V] = svd(X);  % svd well-aligned square
[U,S,V] = svd(Xrot);  % svd rotated square

subplot(2,2,3), semilogy(diag(S),'-ko')
ylim([1.e-16 1.e4]), grid on
set(gca,'YTick',[1.e-16 1.e-12 1.e-8 1.e-4 1. 1.e4])
set(gca,'XTick',[0 250 500 750 1000])

subplot(2,2,4), semilogy(diag(S),'-ko')
ylim([1.e-16 1.e4]), grid on
set(gca,'YTick',[1.e-16 1.e-12 1.e-8 1.e-4 1. 1.e4])
set(gca,'XTick',[0 250 500 750 1000])

author: dynamicslab


text:
function DSO3 = refine(DSO3)

% step 1: compute center of any two orientation connected by an edge

% all combinations of vertices of all tetrahegons
v1 = DSO3.tetra(:,[1 1 1 2 2 3]).';
v2 = DSO3.tetra(:,[2 3 4 3 4 4]).';

% unique edges
vv = sort([v1(:),v2(:)],2);
[v,~,iCenterEdges] = unique(vv,'rows');
iCenterEdges = length(DSO3) + reshape(iCenterEdges,6,[]).';

% step 2: set up corner tetrahegons
% bases
base = iCenterEdges(:,[1 2 3  1 4 5  2 4 6  3 5 6]);
base = reshape(base.',3,[]);
corner = reshape(DSO3.tetra.',1,[]);

% split center octaeder into four tetrahegons
octaeder = iCenterEdges(:,[1 3 4 5  3 4 5 6  1 2 3 4  2 3 4 6]);
octaeder = reshape(octaeder.',4,[]);

% set up tetraeder
newTetra = [[base;corner],octaeder].';

% new vertices as mean of all edges
centerEdges = mean2(DSO3.subSet(v(:,1)),DSO3.subSet(v(:,2)));

% set up refined grid
DSO3.a = [DSO3.a(:);centerEdges.a];
DSO3.b = [DSO3.b(:);centerEdges.b];
DSO3.c = [DSO3.c(:);centerEdges.c];
DSO3.d = [DSO3.d(:);centerEdges.d];
DSO3.i = zeros(size(DSO3.d));

% set up new tetrahegons
DSO3.tetra = sort(newTetra,2);

% compute neighbouring list
DSO3.tetraNeighbour = calcNeighbour(DSO3.tetra);

% the neighbours of the four corner tetrahegons
%1 2 3 4
%cornerNeighbour =
%DSO3.tetraNeighbour

% set up new lookup
%DSO3.lookup = DSO3.lookup*8;
%DSO3.lookup = calcLookUp(DSO3,5*degree);
DSO3.lookup = [];
res = 40*degree;
for i = 1:3
  res = res / 2;
  DSO3.lookup = calcLookUp(DSO3,res);
end

end

% -----------------------------------------------------------------
function v = mean2(v1,v2)
% compute the mean of two orientations

dv = inv(v1) .* v2;
dv.SS = specimenSymmetry;

[dv,omega] = dv.project2FundamentalRegion;

% half angle
dv.a = dv.a .* cos(omega./4) ./ cos(omega./2);
dv.b = dv.b .* sin(omega./4) ./ sin(omega./2);
dv.c = dv.c .* sin(omega./4) ./ sin(omega./2);
dv.d = dv.d .* sin(omega./4) ./ sin(omega./2);

% new vertices
v = v1 .* dv;

end

author: mtex-toolbox

repo: mtex

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


text:
function x = matset_withcrop(x,y,varargin)

nd = length(varargin);
szx = size(x);

if nd == 1 && nnz(szx~=1) <= 1,
  % special case for vectors
  inbounds = varargin{1} >= 1 & varargin{1} <= max(szx);
  x(varargin{1}(inbounds)) = y(inbounds);
else
  idxx = cell(1,nd);
  idxy = cell(1,nd);
  for d = 1:nd,
    inbounds = varargin{d} >= 1 & varargin{d} <= szx(d);
    idxx{d} = varargin{d}(inbounds);
    idxy{d} = inbounds;
  end
  x(idxx{:}) = y(idxy{:});
end

author: kristinbranson

repo: JAABA


text:
function []= prepare_h5_data(data_folder, h5_save_folder,h5_save_list)
% Generate the hdf5 file and a file list used to train the network
% --------------------------------------------------------
% Copyright (c) 2017, Lequan Yu
% Licensed under The MIT License
% --------------------------------------------------------

    % re-sample the data into same resolution, default is False.
    use_isotropic=0;
    patchSize = 64;
    if nargin<3
        data_folder = '../data';
        h5_save_folder = '../h5_data';
        h5_save_list = '../train.list';

    if ~exist(h5_save_folder,'dir')
        mkdir(h5_save_folder);
    end

    addpath('./util');
    fid = fopen(h5_save_list, 'w');

    %% generate the hdf5 file
    for id = 0:0
        tic;
        img_path = [data_folder, '/', 'training_axial_crop_pat', num2str(id), '.nii.gz'];
        seg_path = [data_folder, '/', 'training_axial_crop_pat', num2str(id), '-label.nii.gz'];

        img_nii = load_nii(img_path);
        seg_nii = load_nii(seg_path);

        %pre-process the images (intensity and resize)
        [img,seg] = pre_process_isotropic(img_nii,seg_nii,use_isotropic,id);

        % crop the heart patches (ROI) from whole image and randomly crop patches
        [img,seg,rimgs,rsegs] = Create_ROIs(img, seg,patchSize);

        %% Do data augmentation (permute & rotate & flip)
        % We only do augmentation at axial plane
        % You can do agumentation in all three planes by setting it_p = 1:3
        pp= [1 2 3; 2 3 1; 3 1 2];
        for it_p = 1:1
            permute_img = permute(img,pp(it_p,:));
            permute_seg = permute(seg,pp(it_p,:));
            for it_r = 1:4
                % rotate
                rotate_img = rot90(permute_img, it_r - 1);
                rotate_seg  = rot90(permute_seg, it_r - 1);
                % flip
                for it_f = 0:1
                    if it_f == 0
                        flip_img = rotate_img;
                        flip_seg = rotate_seg;
                    else
                        flip_img = flip(rotate_img, it_f);
                        flip_seg = flip(rotate_seg, it_f);
                    end

                    % h5 file path
                    h5_path    = [h5_save_folder, '/', num2str(id), '_p', num2str(it_p),...
                        '_r', num2str(it_r), 'f', num2str(it_f), '.h5'];
                    if exist(h5_path,'file')
                        delete(h5_path);
                    end

                    [d1, d2, d3, d4, d5] = size(flip_img);
                    data_dims = [d1, d2, d3, d4, d5];
                    [d1, d2, d3, d4, d5] = size(flip_seg);
                    seg_dims = [d1, d2, d3, d4, d5];

                    % store as h5 format
                    h5create(h5_path,'/data',data_dims,'Datatype','single','Deflate',0,'ChunkSize',data_dims);
                    h5create(h5_path,'/label',seg_dims,'Datatype','uint8','Deflate',0,'ChunkSize',seg_dims);
                    h5write(h5_path,'/data', single(flip_img));
                    h5write(h5_path,'/label',uint8(flip_seg));
                    fprintf(fid,'%s\n', h5_path);
                end
            end
        end

        %% save randomly cropped patches.
        % Default, we do not use them to train the model
        for i= 1: 0 %size(random_img, 5)
            h5_path  = [h5_save_folder, '/',num2str(id), '_random_', num2str(i), '.h5'];
            if exist(h5_path,'file')
                delete(h5_path);
            end

            write_img = rimgs{i};
            write_seg  = rsegs{i};
            [d1, d2, d3, d4, d5] = size(write_img);
            data_dims = [d1, d2, d3, d4, d5];
            [d1, d2, d3, d4, d5] = size(write_seg);
            seg_dims = [d1, d2, d3, d4, d5];

            % store as h5 format
            h5create(h5_path,'/data',data_dims,'Datatype','single','Deflate',0,'ChunkSize',data_dims);
            h5create(h5_path,'/label',seg_dims,'Datatype','uint8','Deflate',0,'ChunkSize',seg_dims);
            h5write(h5_path,'/data', single(write_img));
            h5write(h5_path,'/label',uint8(write_seg));
            fprintf(fid,'%s\n', h5_path);
        end
        toc;
    end
    fclose(fid);
end

author: yulequan

repo: HeartSeg

save_path: github-repos/MATLAB/yulequan-HeartSeg


text:
function q=axisAng2Quat(u,theta)
%%AXISANG2QUAT Get a unit quaternion representative of a rotation of a 3D
%              vector an angle of theta counterlockwise (right-handed) or
%              clockwise (left-handed) about an axis given by u. The
%              handedness of theta matches the handedness of the
%              quaternion produced.
%
%INPUTS: u A unit vector representing the axis about which a vector should
%          be rotated using the rotation matrix.
%    theta The angle in radians by which a vector is to be rotated about
%          the axis u when multiplied by the rotation matrix. The
%          handedness of the rotation angle matches the handedness of the
%          quaternion produced.
%
%OUTPUTS: q A 4X1 unit quaternion corresponding to the supplied axis and
%           angle. The handedness of the quaternion (which is important to
%           know when using the quaternion with other functions, such as
%           quatMult) is the same as that of the angle. The quaternion is
%           ordered in terms of hypercomplex numbers as
%           q(1)+i*q(2)+j*q(3)+k*q(4).
%
%Quaternions and rotations are discussed in [1].
%
%REFERENCES:
%[1] M. D. Shuster, "A survey of attitude representations," The Journal of
%    Astronautical Sciences, vol. 41, no. 4, pp. 439-517, Oct. -Dec. 1993.
%    wherein the convention for the quaternions is left-handed.
%
%August 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

q=[cos(theta/2);sin(theta/2)*u];

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function s1=copyStructFieldByField(s1,s2)

if ~strcmp(class(s1),class(s2))
    if ~isobject(s1) && ~isstruct(s2)
        if ~isstruct(s1) && ~isobject(s2)
            return;
        end
    end
end

if isstruct(s1) || isobject(s1)

    if ~isstruct(s2) && ~isobject(s2)
        return;
    end

    for jj=1:length(s2)
        fields = fieldnames(s2(jj));
        for ii=1:length(fields)
            field_class = eval(sprintf('class(s2(jj).%s)',fields{ii}));
            if ~isproperty(s1,fields{ii}) || length(s1)<jj
                if isobject(s1)
                    continue;
                end
                if strcmp(field_class,'cell')
                    field_arg='{}';
                else
                    field_arg='[]';
                end
                eval(sprintf('s1(jj).%s = %s(%s);',fields{ii},field_class,field_arg));
            end
            eval(sprintf('s1(jj).%s = copyStructFieldByField(s1(jj).%s, s2(jj).%s);', fields{ii}, fields{ii}, fields{ii}));
        end
    end

else

    s1 = s2;

end

author: BUNPC


text:
function A = catstruct(varargin)
% CATSTRUCT   Concatenate or merge structures with different fieldnames
%   X = CATSTRUCT(S1,S2,S3,...) merges the structures S1, S2, S3 ...
%   into one new structure X. X contains all fields present in the various
%   structures. An example:
%
%     A.name = 'Me' ;
%     B.income = 99999 ;
%     X = catstruct(A,B)
%     % -> X.name = 'Me' ;
%     %    X.income = 99999 ;
%
%   If a fieldname is not unique among structures (i.e., a fieldname is
%   present in more than one structure), only the value from the last
%   structure with this field is used. In this case, the fields are
%   alphabetically sorted. A warning is issued as well. An axample:
%
%     S1.name = 'Me' ;
%     S2.age  = 20 ; S3.age  = 30 ; S4.age  = 40 ;
%     S5.honest = false ;
%     Y = catstruct(S1,S2,S3,S4,S5) % use value from S4
%
%   The inputs can be array of structures. All structures should have the
%   same size. An example:
%
%     C(1).bb = 1 ; C(2).bb = 2 ;
%     D(1).aa = 3 ; D(2).aa = 4 ;
%     CD = catstruct(C,D) % CD is a 1x2 structure array with fields bb and aa
%
%   The last input can be the string 'sorted'. In this case,
%   CATSTRUCT(S1,S2, ..., 'sorted') will sort the fieldnames alphabetically.
%   To sort the fieldnames of a structure A, you could use
%   CATSTRUCT(A,'sorted') but I recommend ORDERFIELDS for doing that.
%
%   When there is nothing to concatenate, the result will be an empty
%   struct (0x0 struct array with no fields).
%
%   NOTE: To concatenate similar arrays of structs, you can use simple
%   concatenation:
%     A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ;

%   NOTE: This function relies on unique. Matlab changed the behavior of
%   its set functions since 2013a, so this might cause some backward
%   compatibility issues when dulpicated fieldnames are found.
%
%   See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL, ORDERFIELDS

% version 4.1 (feb 2015), tested in R2014a
% (c) Jos van der Geest
% email: jos@jasen.nl

% History
% Created in 2005
% Revisions
%   2.0 (sep 2007) removed bug when dealing with fields containing cell
%                  arrays (Thanks to Rene Willemink)
%   2.1 (sep 2008) added warning and error identifiers
%   2.2 (oct 2008) fixed error when dealing with empty structs (thanks to
%                  Lars Barring)
%   3.0 (mar 2013) fixed problem when the inputs were array of structures
%                  (thanks to Tor Inge Birkenes).
%                  Rephrased the help section as well.
%   4.0 (dec 2013) fixed problem with unique due to version differences in
%                  ML. Unique(...,'last') is no longer the deafult.
%                  (thanks to Isabel P)
%   4.1 (feb 2015) fixed warning with narginchk

narginchk(1,Inf) ;
N = nargin ;

if ~isstruct(varargin{end}),
    if isequal(varargin{end},'sorted'),
        narginchk(2,Inf) ;
        sorted = 1 ;
        N = N-1 ;
    else
        error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ;
    end
else
    sorted = 0 ;
end

sz0 = [] ; % used to check that all inputs have the same size

% used to check for a few trivial cases
NonEmptyInputs = false(N,1) ;
NonEmptyInputsN = 0 ;

% used to collect the fieldnames and the inputs
FN = cell(N,1) ;
VAL = cell(N,1) ;

% parse the inputs
for ii=1:N,
    X = varargin{ii} ;
    if ~isstruct(X),
        error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ;
    end

    if ~isempty(X),
        % empty structs are ignored
        if ii > 1 && ~isempty(sz0)
            if ~isequal(size(X), sz0)
                error('catstruct:UnequalSizes','All structures should have the same size.') ;
            end
        else
            sz0 = size(X) ;
        end
        NonEmptyInputsN = NonEmptyInputsN + 1 ;
        NonEmptyInputs(ii) = true ;
        FN{ii} = fieldnames(X) ;
        VAL{ii} = struct2cell(X) ;
    end
end

if NonEmptyInputsN == 0
    % all structures were empty
    A = struct([]) ;
elseif NonEmptyInputsN == 1,
    % there was only one non-empty structure
    A = varargin{NonEmptyInputs} ;
    if sorted,
        A = orderfields(A) ;
    end
else
    % there is actually something to concatenate
    FN = cat(1,FN{:}) ;
    VAL = cat(1,VAL{:}) ;
    FN = squeeze(FN) ;
    VAL = squeeze(VAL) ;

    [UFN,ind] = unique(FN, 'last') ;
    % If this line errors, due to your matlab version not having UNIQUE
    % accept the 'last' input, use the following line instead
    % [UFN,ind] = unique(FN) ; % earlier ML versions, like 6.5

    if numel(UFN) ~= numel(FN),
        warning('catstruct:DuplicatesFound','Fieldnames are not unique between structures.') ;
        sorted = 1 ;
    end

    if sorted,
        VAL = VAL(ind,:) ;
        FN = FN(ind,:) ;
    end

    A = cell2struct(VAL, FN);
    A = reshape(A, sz0) ; % reshape into original format
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [aTS, bTS] = cRecon(scanParams, fName, slice, compFlag, rotFlag, shifts)

% [aTS, bTS] = cRecon(scanParams, fName, slice, compFlag, shifts);
%
% Load and crop complex recon data.
% Returns pair of single-slice time-series arrays of size:
% ny x nx x nFrames -- aTS is real and bTS is imaginary;
% the latter is returned only if the optional compFlag
% input is set. Time-duration and spatial cropping are performed,
% but no trend removal. Inputs are pFile name [fName], and
% slice index [slice]. Optional input rotFlag performs a 90
% degree rotation to make functional images match inplanes;
% use MakeMovie.m or similar to check if necessary. Optional
% input shifts allows correction of fallback errors or similar.
% Assumes pwd points to valid session
% directory, and uses global mrSESSION.
%
% DBR 12/00

global mrSESSION
if ~exist('compFlag'), compFlag = 0; end
if ~exist('rotFlag'), rotFlag = 0; end

% Do recon, extract complex time series:
[aTS, bTS] = PaulyRecon(scanParams, fName, slice, compFlag, rotFlag);

% Remove junk frames:
f0 = mrSESSION.junkFirstFrames+1;
totalFrames = size(aTS, 1);
nFrames = totalFrames - f0 + 1;

aTS = aTS(f0:end, :, :);
x0 = mrSESSION.tseriesCrop(1, 1) + shifts(2);
xN = mrSESSION.tseriesCrop(2, 1) + shifts(2);
y0 = mrSESSION.tseriesCrop(1, 2) + shifts(1);
yN = mrSESSION.tseriesCrop(2, 2) + shifts(1);
% Crop in time/space and shuffle to standard t-series shape:
aTS = aTS(:, y0:yN, x0:xN);
aTS = reshape(aTS, nFrames, (yN-y0+1)*(xN-x0+1));
if compFlag
  bTS = bTS(f0:end, :, :);
  bTS = bTS(:, y0:yN, x0:xN);
  bTS = reshape(bTS, nFrames, (yN-y0+1)*(xN-x0+1));
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function acq = acqfsn2_vbmc(Xs,vp,gp,optimState,fmu,fs2,fbar,vtot)
%ACQFSN2_VBMC Acquisition fcn. for noisy prospective uncertainty search.

% Xs is in *transformed* coordinates

% Probability density of variational posterior at test points
p = max(vbmc_pdf(vp,Xs,0),realmin);

% Estimate observation noise at test points from nearest neighbor
[~,pos] = min(sq_dist(bsxfun(@rdivide,Xs,optimState.gplengthscale),gp.X_rescaled),[],2);
sn2 = gp.sn2new(pos);

z = optimState.ymax;

% Prospective uncertainty search corrected for noisy observations
acq = -vtot.*(1 - sn2./(vtot+sn2)) .* exp(fbar-z) .* p;

end

%SQ_DIST Compute matrix of all pairwise squared distances between two sets
% of vectors, stored in the columns of the two matrices, a (of size n-by-D)
% and b (of size m-by-D).
function C = sq_dist(a,b)

n = size(a,1);
m = size(b,1);
mu = (m/(n+m))*mean(b,1) + (n/(n+m))*mean(a,1);
a = bsxfun(@minus,a,mu); b = bsxfun(@minus,b,mu);
C = bsxfun(@plus,sum(a.*a,2),bsxfun(@minus,sum(b.*b,2)',2*a*b'));
C = max(C,0);

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function grains = plus(grains,xy)
% shift grains in x/y direction
%
% Syntax
%
%   % shift in x direction
%   grains = grains + [100,0]
%
% Input
%  grains- @grain2d
%  xy - x and y coordinates of the shift
%
% Output
%  grains - @grain2d

if isa(xy,'grain2d'), [xy,grains] = deal(grains,xy); end

if isa(xy,'vector3d'), xy = [xy.x,xy.y]; end

grains.V = grains.V + repmat(xy,size(grains.V,1),1);


text:
seed = 1235;
rng(seed);

[Batch,~,~] = generate_data(5000,[10,10]',[10,10]',0.0,[0,0]',[1,1]',0,1,0,0,0);
% generate_data(T,Nneurons,Dt,NeuronNoise,SeqNoiseTime,SeqNoiseNeuron,shared,diff,stretch,bin,seed)

L = 250; K = 4;
figure(1)
[W,H,temp,Loadings,Power] = seqNMF(Batch, ...
    'K', K, 'L', L, 'lambda',.0005, ...
    'showPlot', 1, 'maxiter', 100, 'shift', 1,...
    'lambdaOrthoW', 0,'lambdaOrthoH',0);
SimpleWHPlot(W,H,[],0)
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'SharedNeurons';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_noshared.pdf']));

[W,H,temp,Loadings,Power] = seqNMF(Batch, ...
    'K', K, 'L', L, 'lambda',.00000, ...
    'showPlot', 1, 'maxiter', 100, 'shift', 1,...
    'lambdaOrthoW', 0,'lambdaOrthoH',0);
SimpleWHPlot(W,H,[],0)
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'SharedNeurons';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_cnmf.pdf']));

[Batch,~,~] = generate_data(5000,[10,10]',[10,10]',0.0,[0,0]',[1,1]',1,1,0,0,0);
[W,H,temp,Loadings,Power] = seqNMF(Batch, ...
    'K', K, 'L', L, 'lambda',.0005, ...
    'showPlot', 1, 'maxiter', 100, 'shift', 1,...
    'lambdaOrthoW', 0,'lambdaOrthoH',0);
SimpleWHPlot(W,H,[],0)
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'SharedNeurons';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_sharedDiff.pdf']));

[Batch,~,~] = generate_data(5000,[10,10]',[10,10]',0.0,[0,0]',[1,1]',1,0,0,0,0);
[W,H,temp,Loadings,Power] = seqNMF(Batch, ...
    'K', K, 'L', L, 'lambda',.0005, ...
    'showPlot', 1, 'maxiter', 100, 'shift', 1,...
    'lambdaOrthoW', 0,'lambdaOrthoH',0);
SimpleWHPlot(W,H,[],0)
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'SharedNeurons';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_sharedSame.pdf']));
%% Example parts-based vs events-based factorizations (W or H soft orthogonality)

% change these parameters to switch between parts-based and events-based
lambdaOrthoH = 0; % favor events-based (these can take any value, don't need to be zero and one)
lambdaOrthoW = 1; % favor parts-based

% rng(236); % fixed rng seed for reproduceability
X = trainNEURAL;
K = 3;
L = 2/3; % units of seconds
Lneural = ceil(L*VIDEOfs);
Lsong = ceil(L*SONGfs);
shg
display('For parts-based, set lambdaOrthoW=1 and lambdaOrthoH=0; for events-based, vice versa')
[W, H, ~,loadings,power]= seqNMF(X,'K',K,'L',Lneural,...
            'lambda', .00, 'maxiter', 100, 'showPlot', 1,...
            'lambdaOrthoH', lambdaOrthoH, 'lambdaOrthoW', lambdaOrthoW);

p = .05; % desired p value for factors

display('Testing significance of factors on held-out data')
[pvals,is_significant] = test_significance(testNEURAL,W,p);

W = W(:,is_significant,:);
H = H(is_significant,:);

% plot, sorting neurons by latency within each factor
[max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
indSort = hybrid(:,3);
tstart = 180; % plot data starting at this timebin
figure; WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end),0)
title(['lambdaOrthoW=' num2str(lambdaOrthoW) ', lambdaOrthoH=' num2str(lambdaOrthoH)])

% savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
% saveTitle = 'PartsBasedOrthW1';
% papersize = [8 6];
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_raw.pdf']));

% figure; WHPlot(W(indSort,:,:),H(:,tstart:end), ...
%     helper.reconstruct(W(indSort,:,:),H(:,tstart:end)),...
%     0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end),0)
% title('SeqNMF reconstruction')

% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_reconstruction.pdf']));

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function [dat, mrk] = proc_resample(dat, target_fs, varargin)
%PROC_RESAMPLE - resample the EEG data to a target sampling frequency
%
%Synposis:
% [dat, mrk] = proc_resample(dat, target_fs, <mrk, N>)
%
%Arguments:
%     dat - Structure with fields x and fs. This fields
%     will be overwritten
%     target_fs - the target sampling frequency
%
%     optional (given with keywords or as opt.xxx):
%     mrk - Structure with field time. It will be updated.
%     N - remove the first and last N samples from the resampled data to
%     avoid edge effects (default is 0)
%
%Returns:
%     dat: updated data structure
%     mrk: updated marker structure
%
%Description:
% Resamples the field dat.x such that is has the desired sampling frequency.
% dat.t and mrk (if given) will be updated as well.

% Sven Daehne, 06-2011

props= {'N'     0     '!INT[1]'
        'mrk'   []    'STRUCT'};

if nargin==0,
  dat= props; return
end

misc_checkType(dat, 'STRUCT(x fs)');
dat = misc_history(dat);

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

N = opt.N;                  % # samples
Nt = N/target_fs*1000;      % corresponding time in ms
mrk = opt.mrk;

p = round(target_fs*1000);
q = round(dat.fs*1000);

if ndims(dat.x)==3
  nTrials = size(dat.x,3);
  for n= 1:nTrials,
    X(:,:,n) = resample(dat.x(:,:,n), p, q);
  end
else
  X = resample(dat.x, p, q);
end
dat.x = X;
dat.fs = target_fs;
n_samples = size(dat.x,1);
if isfield(dat, 't'),
  t = linspace(dat.t(1), dat.t(end), n_samples);
  dat.t = t; % time in ms
end

% remove the first and the last N samples to avoid edge effects
dat.x = dat.x((N+1):end-N, :, :);
if isfield(dat, 't')
    dat.t = dat.t((N+1):end-N);
end

% Adjust marker if necessary
if N>0 && ~isempty(mrk)
    mrk.time = mrk.time-Nt;
    % Remove markers that refer to data that was cut-off the beginning or
    % end
    rem_idx = [find(mrk.time<=0) find(mrk.time>length(dat.x)/dat.fs*1000)];
    mrk = mrk_selectEvents(mrk,'not',rem_idx);
end

author: bbci


text:
%test_VAcPlant(t, ID) is an example, which shows how to use the VAModel in MATLAB
%In this example, Luyben's multi-loop control structure is implemented and tested
%Eight process disturbances are available and transients are generated at the end of simulation
%On a PIII-1GHz PC, the VA plant model runs approaximately 80 times faster than the real time
%To use this example:
%       t sets the simulation time (in munite)
%       ID is an integer selected between 0 and 8
%           0: no disturbance
%           1: setpoint of the reactor outlet temperature decreases 8 degC (from 159 to 151)
%           2: setpoint of the reactor outlet temperature increases 6 degC (from 159 to 165)
%           3: setpoint of the H2O composition in the column bottom increases 9% (from 9% to 18%)
%           4: the vaporizer liqiud inlet flowrate increases 0.44 kmol/min (from 2.2 to 2.64)
%           5: HAc fresh feed stream lost for 5 minutes
%           6: O2 fresh feed stream lost for 5 minutes
%           7: C2H6 composition changes from 0.001 to 0.003 in the C2H4 fresh feed stream
%           8: column feed lost for 5 minutes
%       Note: in this example, all the disturbances occur at 10 minute after the simulation starts
%-----------------------------------------------------------------------------
%Copyright: Rong Chen and Kedar David, June 2002
%-----------------------------------------------------------------------------
function test_VAcPlant(minute,selected_ID)

%-----------------------------------------------------------------------------
close all
format short g
%-----------------------------------------------------------------------------

%-----------------------------------------------------------------------------
%The recommended simulation sampling time is 1/3 second, which corresponds to a frequency of 180 samples in one minute
model_sampling_frequency=180;
%The recommended historical data sampling time is 1 minute, which corresponds to a frequency of 1 sample in one minute
storage_sampling_frequency=1;
%The recommended display sampling time is 1 minute, which corresponds to a frequency of 1 sample in one minute
display_sampling_frequency=1;
%-----------------------------------------------------------------------------

%-----------------------------------------------------------------------------
% Inintialze process states and MV's
% 246 state variables
states=zeros(246,1);
% 26 manipulated variables
MVs=zeros(26,1);
% set the current process time to 0 (in minute)
time=0;
% set initialization flag 1
is_initial=1;
% set disturbance_ID 0
disturbance_ID=0;
% get the base operation
[dstatedt,states,MVs,y_ss]=VAModel(states,MVs,time,is_initial,disturbance_ID);
%-----------------------------------------------------------------------------

%-----------------------------------------------------------------------------
% Implement Luyben and Tyreus's control structure for the VA plant
% Note that: in this example,
%   1. for each transmitter, a 3 second time lag is used
%   2. for two analyzer transmitters (on the gas recycle and column bottom), a 10 minute deadtime is used
%   3. for other transmitters, no deadtime is used
%   4. for three controllers using analyzers on the gas recycle and column bottom, a 10 minute sampling time is used
%   5. for other controllers, a 1 second sampling time is used
%   6. for each transmitter, a 1 second sampling time is used
%-----------------------------------------------------------------------------
transmitter_lag=0.05; %in minute
transmitter_deadtime=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 10 0 10 0 0 0 0 0 0]; %in minute
transmitter_sampling_frequency=[60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60];
controller_sampling_frequency=[60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 60 0.1 0.1 60 0.1 60 60 60 60 60 60];
%-----------------------------------------------------------------------------
%transmitter and controller initialization
%SP: setpoint (not scaled)
%K: controller gain which is dimentionless
%Ti: reset time (in minute)
%act: 1 if positive process gain, -1 if negative process gain
%mode: 1 for automatic, 2 for manual
%Ponly: 1 for proportional only control, 0 for PI control
%hc: controller sampling time (in minute), should be 1/"controller_sampling_frequency"
%ht: transmitter sampling time (in minute), should be 1/"transmitter_sampling_frequency"
%uLO: low limit of MV
%uHI: high limit of MV
%ded: should be "transmitter_deadtime" defined above
%tau: should be "transmitter_lag" defined above
%yLO: low limit of measurement
%yHI: high limit of measurement
%Nded: measurement deadtime storage, which should be 1+transmitter_deadtime*transmitter_sampling_frequency
%Ntau: measurement time constant storage, which should be 1 all the time
%-----------------------------------------------------------------------------
%MV1: F_O2 - O2 composition
SP_O2=0.075;
K_O2=5;
Ti_O2=10;
act_O2=1;
mode_O2=1;
Ponly_O2=0;
ht_O2=1/transmitter_sampling_frequency(1);
hc_O2=1/controller_sampling_frequency(1);
uLO_O2=0;
uHI_O2=2.268;
ded_O2=transmitter_deadtime(1);
tau_O2=transmitter_lag;
yLO_O2=0;
yHI_O2=0.2;
Nded_O2=1+transmitter_deadtime(1)*transmitter_sampling_frequency(1);
Ntau_O2=1;
xxx_O2(1)= (y_ss(37) -yLO_O2) / (yHI_O2 -yLO_O2);
for i=2:Nded_O2
    xxx_O2(i)=xxx_O2(1);
end
for i=1:Ntau_O2
    yyy_O2(i)= xxx_O2(Nded_O2);
end
xi_O2= (MVs(1) - uLO_O2) /(uHI_O2 -uLO_O2)-K_O2*act_O2*((SP_O2-yLO_O2)/(yHI_O2-yLO_O2)-yyy_O2(Ntau_O2));
%-----------------------------------------------------------------------------
%controller initialization
%MV2: F_C2H4 - Recycle Pressure
SP_C2H4=128;
K_C2H4=0.3;
Ti_C2H4=20;
act_C2H4=1;
mode_C2H4=1;
Ponly_C2H4=0;
ht_C2H4=1/transmitter_sampling_frequency(2);
hc_C2H4=1/controller_sampling_frequency(2);
uLO_C2H4=0;
uHI_C2H4=7.56;
ded_C2H4=transmitter_deadtime(2);
tau_C2H4=transmitter_lag;
yLO_C2H4=0;
yHI_C2H4=200;
Nded_C2H4=1+transmitter_deadtime(2)*transmitter_sampling_frequency(2);
Ntau_C2H4=1;
xxx_C2H4(1)= (y_ss(12) -yLO_C2H4) / (yHI_C2H4 -yLO_C2H4);
for i=2:Nded_C2H4
    xxx_C2H4(i)= xxx_C2H4(1);
end
for i=1:Ntau_C2H4
    yyy_C2H4(i)= xxx_C2H4(Nded_C2H4);
end
xi_C2H4= (MVs(2) - uLO_C2H4) /(uHI_C2H4 -uLO_C2H4)-K_C2H4*act_C2H4*((SP_C2H4-yLO_C2H4)/(yHI_C2H4-yLO_C2H4)-yyy_C2H4(Ntau_C2H4));
%-----------------------------------------------------------------------------
%controller initialization
%MV3: HAc feed stream - HAc tank level
SP_HAc=0.5;
K_HAc=2;
Ti_HAc=100;
act_HAc=1;
mode_HAc=1;
Ponly_HAc=1;
ht_HAc=1/transmitter_sampling_frequency(3);
hc_HAc=1/controller_sampling_frequency(3);
uLO_HAc=0;
uHI_HAc=4.536;
ded_HAc=transmitter_deadtime(3);
tau_HAc=transmitter_lag;
yLO_HAc=0;
yHI_HAc=1;
Nded_HAc=1+transmitter_deadtime(3)*transmitter_sampling_frequency(3);
Ntau_HAc=1;
xxx_HAc(1)= (y_ss(23) -yLO_HAc) / (yHI_HAc -yLO_HAc);
for i=2:Nded_HAc
    xxx_HAc(i)= xxx_HAc(1);
end
for i=1:Ntau_HAc
    yyy_HAc(i)= xxx_HAc(Nded_HAc);
end
xi_HAc= (MVs(3) - uLO_HAc) /(uHI_HAc -uLO_HAc)-K_HAc*act_HAc*((SP_HAc-yLO_HAc)/(yHI_HAc-yLO_HAc)-yyy_HAc(Ntau_HAc));
%-----------------------------------------------------------------------------
%controller initialization
%MV4: Qv - Level
SP_LevelVap=0.7;
K_LevelVap=0.1;
Ti_LevelVap=30;
act_LevelVap=-1;
mode_LevelVap=1;
Ponly_LevelVap=0;
ht_LevelVap=1/transmitter_sampling_frequency(4);
hc_LevelVap=1/controller_sampling_frequency(4);
uLO_LevelVap=0;
uHI_LevelVap=1433400;
ded_LevelVap=transmitter_deadtime(4);
tau_LevelVap=transmitter_lag;
yLO_LevelVap=0;
yHI_LevelVap=1;
Nded_LevelVap=1+transmitter_deadtime(4)*transmitter_sampling_frequency(4);
Ntau_LevelVap=1;
xxx_LevelVap(1)= (y_ss(2) -yLO_LevelVap) / (yHI_LevelVap -yLO_LevelVap);
for i=2:Nded_LevelVap
    xxx_LevelVap(i)= xxx_LevelVap(1);
end
for i=1:Ntau_LevelVap
    yyy_LevelVap(i)= xxx_LevelVap(Nded_LevelVap);
end
xi_LevelVap= (MVs(4) - uLO_LevelVap) /(uHI_LevelVap -uLO_LevelVap)-K_LevelVap*act_LevelVap*((SP_LevelVap-yLO_LevelVap)/(yHI_LevelVap-yLO_LevelVap)-yyy_LevelVap(Ntau_LevelVap));
%-----------------------------------------------------------------------------
%controller initialization
%MV5: F_Vaporizer - Vaporizer Pressure
SP_PresVap=128;
K_PresVap=5;
Ti_PresVap=10;
act_PresVap=-1;
mode_PresVap=1;
Ponly_PresVap=0;
ht_PresVap=1/transmitter_sampling_frequency(5);
hc_PresVap=1/controller_sampling_frequency(5);
uLO_PresVap=0;
uHI_PresVap=50;
ded_PresVap=transmitter_deadtime(5);
tau_PresVap=transmitter_lag;
yLO_PresVap=0;
yHI_PresVap=200;
Nded_PresVap=1+transmitter_deadtime(5)*transmitter_sampling_frequency(5);
Ntau_PresVap=1;
xxx_PresVap(1)= (y_ss(1) -yLO_PresVap) / (yHI_PresVap -yLO_PresVap);
for i=2:Nded_PresVap
    xxx_PresVap(i)= xxx_PresVap(1);
end
for i=1:Ntau_PresVap
    yyy_PresVap(i)= xxx_PresVap(Nded_PresVap);
end
xi_PresVap= (MVs(5) - uLO_PresVap) /(uHI_PresVap -uLO_PresVap)-K_PresVap*act_PresVap*((SP_PresVap-yLO_PresVap)/(yHI_PresVap-yLO_PresVap)-yyy_PresVap(Ntau_PresVap));
%-----------------------------------------------------------------------------
%controller initialization
%MV6: Q_heater - Vaporizer Outlet Temperature
SP_ReactorInletTemp=150;
K_ReactorInletTemp=1;
Ti_ReactorInletTemp=5;
act_ReactorInletTemp=1;
mode_ReactorInletTemp=1;
Ponly_ReactorInletTemp=0;
ht_ReactorInletTemp=1/transmitter_sampling_frequency(6);
hc_ReactorInletTemp=1/controller_sampling_frequency(6);
uLO_ReactorInletTemp=0;
uHI_ReactorInletTemp=15000;
ded_ReactorInletTemp=transmitter_deadtime(6);
tau_ReactorInletTemp=transmitter_lag;
yLO_ReactorInletTemp=120;
yHI_ReactorInletTemp=170;
Nded_ReactorInletTemp=1+transmitter_deadtime(6)*transmitter_sampling_frequency(6);
Ntau_ReactorInletTemp=1;
xxx_ReactorInletTemp(1)= (y_ss(4) -yLO_ReactorInletTemp) / (yHI_ReactorInletTemp -yLO_ReactorInletTemp);
for i=2:Nded_ReactorInletTemp
    xxx_ReactorInletTemp(i)= xxx_ReactorInletTemp(1);
end
for i=1:Ntau_ReactorInletTemp
    yyy_ReactorInletTemp(i)= xxx_ReactorInletTemp(Nded_ReactorInletTemp);
end
xi_ReactorInletTemp= (MVs(6) - uLO_ReactorInletTemp) /(uHI_ReactorInletTemp -uLO_ReactorInletTemp)-K_ReactorInletTemp*act_ReactorInletTemp*((SP_ReactorInletTemp-yLO_ReactorInletTemp)/(yHI_ReactorInletTemp-yLO_ReactorInletTemp)-yyy_ReactorInletTemp(Ntau_ReactorInletTemp));
%-----------------------------------------------------------------------------
%controller initialization
%MV7: T_Shell - T_RCTOUT
SP_TRCT=159.17;
K_TRCT=3;
Ti_TRCT=10;
act_TRCT=1;
mode_TRCT=1;
Ponly_TRCT=0;
ht_TRCT=1/transmitter_sampling_frequency(7);
hc_TRCT=1/controller_sampling_frequency(7);
uLO_TRCT=110;
uHI_TRCT=150;
ded_TRCT=transmitter_deadtime(7);
tau_TRCT=transmitter_lag;
yLO_TRCT=0;
yHI_TRCT=200;
Nded_TRCT=1+transmitter_deadtime(7)*transmitter_sampling_frequency(7);
Ntau_TRCT=1;
xxx_TRCT(1)= (y_ss(5)-yLO_TRCT) / (yHI_TRCT -yLO_TRCT);
for i=2:Nded_TRCT
    xxx_TRCT(i)= xxx_TRCT(1);
end
for i=1:Ntau_TRCT
    yyy_TRCT(i)= xxx_TRCT(Nded_TRCT);
end
xi_TRCT= (MVs(7) - uLO_TRCT) /(uHI_TRCT -uLO_TRCT)-K_TRCT*act_TRCT*((SP_TRCT-yLO_TRCT)/(yHI_TRCT-yLO_TRCT)-yyy_TRCT(Ntau_TRCT));
%-----------------------------------------------------------------------------
%controller initialization
%MV8: Separator Liquid Exit - Separator Level
SP_LevelSep=0.5;
K_LevelSep=5;
Ti_LevelSep=100;
act_LevelSep=-1;
mode_LevelSep=1;
Ponly_LevelSep=1;
ht_LevelSep=1/transmitter_sampling_frequency(8);
hc_LevelSep=1/controller_sampling_frequency(8);
uLO_LevelSep=0;
uHI_LevelSep=4.536;
ded_LevelSep=transmitter_deadtime(8);
tau_LevelSep=transmitter_lag;
yLO_LevelSep=0;
yHI_LevelSep=1;
Nded_LevelSep=1+transmitter_deadtime(8)*transmitter_sampling_frequency(8);
Ntau_LevelSep=1;
xxx_LevelSep(1)= (y_ss(9)-yLO_LevelSep) / (yHI_LevelSep -yLO_LevelSep);
for i=2:Nded_LevelSep
    xxx_LevelSep(i)= xxx_LevelSep(1);
end
for i=1:Ntau_LevelSep
    yyy_LevelSep(i)= xxx_LevelSep(Nded_LevelSep);
end
xi_LevelSep= (MVs(8) - uLO_LevelSep) /(uHI_LevelSep -uLO_LevelSep)-K_LevelSep*act_LevelSep*((SP_LevelSep-yLO_LevelSep)/(yHI_LevelSep-yLO_LevelSep)-yyy_LevelSep(Ntau_LevelSep));
%-----------------------------------------------------------------------------
%controller initialization
%MV9: T_Shell_Separator - Separator Liquid Temperature
SP_TempSep=40;
K_TempSep=5;
Ti_TempSep=20;
act_TempSep=-1;
mode_TempSep=1;
Ponly_TempSep=0;
ht_TempSep=1/transmitter_sampling_frequency(9);
hc_TempSep=1/controller_sampling_frequency(9);
uLO_TempSep=0;
uHI_TempSep=80;
ded_TempSep=transmitter_deadtime(9);
tau_TempSep=transmitter_lag;
yLO_TempSep=0;
yHI_TempSep=80;
Nded_TempSep=1+transmitter_deadtime(9)*transmitter_sampling_frequency(9);
Ntau_TempSep=1;
xxx_TempSep(1)= (y_ss(10)-yLO_TempSep) / (yHI_TempSep -yLO_TempSep);
for i=2:Nded_TempSep
    xxx_TempSep(i)= xxx_TempSep(1);
end
for i=1:Ntau_TempSep
    yyy_TempSep(i)= xxx_TempSep(Nded_TempSep);
end
xi_TempSep= (MVs(9) - uLO_TempSep) /(uHI_TempSep -uLO_TempSep)-K_TempSep*act_TempSep*((SP_TempSep-yLO_TempSep)/(yHI_TempSep-yLO_TempSep)-yyy_TempSep(Ntau_TempSep));
%-----------------------------------------------------------------------------
%controller initialization
%MV10: Separator Vapor Exit is fixed at 16.1026 kmol/min
%-----------------------------------------------------------------------------
%controller initialization
%MV11: Q_Compressor_Cooler - Compressor Outlet Temperature (80degC)
SP_CompOutletTemp=80;
K_CompOutletTemp=1;
Ti_CompOutletTemp=5;
act_CompOutletTemp=-1;
mode_CompOutletTemp=1;
Ponly_CompOutletTemp=0;
ht_CompOutletTemp=1/transmitter_sampling_frequency(11);
hc_CompOutletTemp=1/controller_sampling_frequency(11);
uLO_CompOutletTemp=0;
uHI_CompOutletTemp=50000;
ded_CompOutletTemp=transmitter_deadtime(11);
tau_CompOutletTemp=transmitter_lag;
yLO_CompOutletTemp=70;
yHI_CompOutletTemp=90;
Nded_CompOutletTemp=1+transmitter_deadtime(11)*transmitter_sampling_frequency(11);
Ntau_CompOutletTemp=1;
xxx_CompOutletTemp(1)= (y_ss(11) -yLO_CompOutletTemp) / (yHI_CompOutletTemp -yLO_CompOutletTemp);
for i=2:Nded_CompOutletTemp
    xxx_CompOutletTemp(i)= xxx_CompOutletTemp(1);
end
for i=1:Ntau_CompOutletTemp
    yyy_CompOutletTemp(i)= xxx_CompOutletTemp(Nded_CompOutletTemp);
end
xi_CompOutletTemp= (MVs(11) - uLO_CompOutletTemp) /(uHI_CompOutletTemp -uLO_CompOutletTemp)-K_CompOutletTemp*act_CompOutletTemp*((SP_CompOutletTemp-yLO_CompOutletTemp)/(yHI_CompOutletTemp-yLO_CompOutletTemp)-yyy_CompOutletTemp(Ntau_CompOutletTemp));
%-----------------------------------------------------------------------------
%controller initialization
%MV12: Absorber Liquid Exit - Absorber Level
SP_LevelAbs=0.5;
K_LevelAbs=5;
Ti_LevelAbs=100;
act_LevelAbs=-1;
mode_LevelAbs=1;
Ponly_LevelAbs=1;
ht_LevelAbs=1/transmitter_sampling_frequency(12);
hc_LevelAbs=1/controller_sampling_frequency(12);
uLO_LevelAbs=0;
uHI_LevelAbs=4.536;
ded_LevelAbs=transmitter_deadtime(12);
tau_LevelAbs=transmitter_lag;
yLO_LevelAbs=0;
yHI_LevelAbs=1;
Nded_LevelAbs=1+transmitter_deadtime(12)*transmitter_sampling_frequency(12);
Ntau_LevelAbs=1;
xxx_LevelAbs(1)= (y_ss(13)-yLO_LevelAbs) / (yHI_LevelAbs -yLO_LevelAbs);
for i=2:Nded_LevelAbs
    xxx_LevelAbs(i)= xxx_LevelAbs(1);
end
for i=1:Ntau_LevelAbs
    yyy_LevelAbs(i)= xxx_LevelAbs(Nded_LevelAbs);
end
xi_LevelAbs= (MVs(12) - uLO_LevelAbs) /(uHI_LevelAbs -uLO_LevelAbs)-K_LevelAbs*act_LevelAbs*((SP_LevelAbs-yLO_LevelAbs)/(yHI_LevelAbs-yLO_LevelAbs)-yyy_LevelAbs(Ntau_LevelAbs));
%-----------------------------------------------------------------------------
%controller initialization
%MV13: Circulation Flowrate is fixed to 15.1198 kmol/min
%-----------------------------------------------------------------------------
%controller initialization
%MV14: Q_Cooler_Circulation - Circulation Temperature
SP_CircOutletTemp=25;
K_CircOutletTemp=1;
Ti_CircOutletTemp=5;
act_CircOutletTemp=-1;
mode_CircOutletTemp=1;
Ponly_CircOutletTemp=0;
ht_CircOutletTemp=1/transmitter_sampling_frequency(14);
hc_CircOutletTemp=1/controller_sampling_frequency(14);
uLO_CircOutletTemp=0;
uHI_CircOutletTemp=30000;
ded_CircOutletTemp=transmitter_deadtime(14);
tau_CircOutletTemp=transmitter_lag;
yLO_CircOutletTemp=10;
yHI_CircOutletTemp=40;
Nded_CircOutletTemp=1+transmitter_deadtime(14)*transmitter_sampling_frequency(14);
Ntau_CircOutletTemp=1;
xxx_CircOutletTemp(1)= (y_ss(14) -yLO_CircOutletTemp) / (yHI_CircOutletTemp -yLO_CircOutletTemp);
for i=2:Nded_CircOutletTemp
    xxx_CircOutletTemp(i)= xxx_CircOutletTemp(1);
end
for i=1:Ntau_CircOutletTemp
    yyy_CircOutletTemp(i)= xxx_CircOutletTemp(Nded_CircOutletTemp);
end
xi_CircOutletTemp= (MVs(14) - uLO_CircOutletTemp) /(uHI_CircOutletTemp -uLO_CircOutletTemp)-K_CircOutletTemp*act_CircOutletTemp*((SP_CircOutletTemp-yLO_CircOutletTemp)/(yHI_CircOutletTemp-yLO_CircOutletTemp)-yyy_CircOutletTemp(Ntau_CircOutletTemp));
%-----------------------------------------------------------------------------
%controller initialization
%MV15: Scrub Flowrate is fixed to 0.756 kmol/min
%-----------------------------------------------------------------------------
%controller initialization
%MV16: Q_Cooler_Scrub - Scrub Temperature
SP_ScrubOutletTemp=25;
K_ScrubOutletTemp=1;
Ti_ScrubOutletTemp=5;
act_ScrubOutletTemp=-1;
mode_ScrubOutletTemp=1;
Ponly_ScrubOutletTemp=0;
ht_ScrubOutletTemp=1/transmitter_sampling_frequency(16);
hc_ScrubOutletTemp=1/controller_sampling_frequency(16);
uLO_ScrubOutletTemp=0;
uHI_ScrubOutletTemp=5000;
ded_ScrubOutletTemp=transmitter_deadtime(16);
tau_ScrubOutletTemp=transmitter_lag;
yLO_ScrubOutletTemp=10;
yHI_ScrubOutletTemp=40;
Nded_ScrubOutletTemp=1+transmitter_deadtime(16)*transmitter_sampling_frequency(16);
Ntau_ScrubOutletTemp=1;
xxx_ScrubOutletTemp(1)= (y_ss(15) -yLO_ScrubOutletTemp) / (yHI_ScrubOutletTemp -yLO_ScrubOutletTemp);
for i=2:Nded_ScrubOutletTemp
    xxx_ScrubOutletTemp(i)= xxx_ScrubOutletTemp(1);
end
for i=1:Ntau_ScrubOutletTemp
    yyy_ScrubOutletTemp(i)= xxx_ScrubOutletTemp(Nded_ScrubOutletTemp);
end
xi_ScrubOutletTemp= (MVs(16) - uLO_ScrubOutletTemp) /(uHI_ScrubOutletTemp -uLO_ScrubOutletTemp)-K_ScrubOutletTemp*act_ScrubOutletTemp*((SP_ScrubOutletTemp-yLO_ScrubOutletTemp)/(yHI_ScrubOutletTemp-yLO_ScrubOutletTemp)-yyy_ScrubOutletTemp(Ntau_ScrubOutletTemp));
%-----------------------------------------------------------------------------
%controller initialization
%MV17: CO2 - CO2 composition
SP_CO2=0.0076393;
K_CO2=1;
Ti_CO2=100;
act_CO2=-1;
mode_CO2=1;
Ponly_CO2=1;
ht_CO2=1/transmitter_sampling_frequency(17);
hc_CO2=1/controller_sampling_frequency(17);
uLO_CO2=0;
uHI_CO2=22.68;
ded_CO2=transmitter_deadtime(17);
tau_CO2=transmitter_lag;
yLO_CO2=0;
yHI_CO2=0.5;
Nded_CO2=1+transmitter_deadtime(17)*transmitter_sampling_frequency(17);
Ntau_CO2=1;
xxx_CO2(1)= (y_ss(31) -yLO_CO2) / (yHI_CO2 -yLO_CO2);
for i=2:Nded_CO2
    xxx_CO2(i)= xxx_CO2(1);
end
for i=1:Ntau_CO2
    yyy_CO2(i)= xxx_CO2(Nded_CO2);
end
xi_CO2= (MVs(17) - uLO_CO2) /(uHI_CO2 -uLO_CO2)-K_CO2*act_CO2*((SP_CO2-yLO_CO2)/(yHI_CO2-yLO_CO2)-yyy_CO2(Ntau_CO2));
%-----------------------------------------------------------------------------
%controller initialization
%MV18: Purge - C2H6 composition
SP_C2H6=0.25;
K_C2H6=1;
Ti_C2H6=100;
act_C2H6=-1;
mode_C2H6=1;
Ponly_C2H6=1;
ht_C2H6=1/transmitter_sampling_frequency(18);
hc_C2H6=1/controller_sampling_frequency(18);
uLO_C2H6=0;
uHI_C2H6=0.02268;
ded_C2H6=transmitter_deadtime(18);
tau_C2H6=transmitter_lag;
yLO_C2H6=0;
yHI_C2H6=1;
Nded_C2H6=1+transmitter_deadtime(18)*transmitter_sampling_frequency(18);
Ntau_C2H6=1;
xxx_C2H6(1)= (y_ss(33) -yLO_C2H6) / (yHI_C2H6 -yLO_C2H6);
for i=2:Nded_C2H6
    xxx_C2H6(i)= xxx_C2H6(1);
end
for i=1:Ntau_C2H6
    yyy_C2H6(i)= xxx_C2H6(Nded_C2H6);
end
xi_C2H6= (MVs(18) - uLO_C2H6) /(uHI_C2H6 -uLO_C2H6)-K_C2H6*act_C2H6*((SP_C2H6-yLO_C2H6)/(yHI_C2H6-yLO_C2H6)-yyy_C2H6(Ntau_C2H6));
%-----------------------------------------------------------------------------
%controller initialization
%MV19: bypass - FEHE Cold Exit Temp
SP_FEHE=134;
K_FEHE=5;
Ti_FEHE=10;
act_FEHE=1;
mode_FEHE=1;
Ponly_FEHE=0;
ht_FEHE=1/transmitter_sampling_frequency(19);
hc_FEHE=1/controller_sampling_frequency(19);
uLO_FEHE=0;
uHI_FEHE=1;
ded_FEHE=transmitter_deadtime(19);
tau_FEHE=transmitter_lag;
yLO_FEHE=0;
yHI_FEHE=200;
Nded_FEHE=1+transmitter_deadtime(19)*transmitter_sampling_frequency(19);
Ntau_FEHE=1;
xxx_FEHE(1)= (y_ss(8) -yLO_FEHE) / (yHI_FEHE -yLO_FEHE);
for i=2:Nded_FEHE
    xxx_FEHE(i)= xxx_FEHE(1);
end
for i=1:Ntau_FEHE
    yyy_FEHE(i)= xxx_FEHE(Nded_FEHE);
end
xi_FEHE= (MVs(19) - uLO_FEHE) /(uHI_FEHE -uLO_FEHE)-K_FEHE*act_FEHE*((SP_FEHE-yLO_FEHE)/(yHI_FEHE-yLO_FEHE)-yyy_FEHE(Ntau_FEHE));
%-----------------------------------------------------------------------------
%controller initialization
%MV20: LR - %H2O in bottom
SP_H2OCol=0.09344;
K_H2OCol=0.5;
Ti_H2OCol=60;
act_H2OCol=-1;
mode_H2OCol=1;
Ponly_H2OCol=0;
ht_H2OCol=1/transmitter_sampling_frequency(20);
hc_H2OCol=1/controller_sampling_frequency(20);
uLO_H2OCol=0;
uHI_H2OCol=7.56;
ded_H2OCol=transmitter_deadtime(20);
tau_H2OCol=transmitter_lag;
yLO_H2OCol=0;
yHI_H2OCol=0.2;
Nded_H2OCol=1+transmitter_deadtime(20)*transmitter_sampling_frequency(20);
Ntau_H2OCol=1;
xxx_H2OCol(1)= (y_ss(28) -yLO_H2OCol) / (yHI_H2OCol -yLO_H2OCol);
for i=2:Nded_H2OCol
    xxx_H2OCol(i)= xxx_H2OCol(1);
end
for i=1:Ntau_H2OCol
    yyy_H2OCol(i)= xxx_H2OCol(Nded_H2OCol);
end
xi_H2OCol= (MVs(20) - uLO_H2OCol) /(uHI_H2OCol -uLO_H2OCol)-K_H2OCol*act_H2OCol*((SP_H2OCol-yLO_H2OCol)/(yHI_H2OCol-yLO_H2OCol)-yyy_H2OCol(Ntau_H2OCol));
%-----------------------------------------------------------------------------
%controller initialization
%MV21: Qr - Temp
SP_TempCol=110;
K_TempCol=20;
Ti_TempCol=30;
act_TempCol=1;
mode_TempCol=1;
Ponly_TempCol=0;
ht_TempCol=1/transmitter_sampling_frequency(21);
hc_TempCol=1/controller_sampling_frequency(21);
uLO_TempCol=0;
uHI_TempCol=100000;
ded_TempCol=transmitter_deadtime(21);
tau_TempCol=transmitter_lag;
yLO_TempCol=0;
yHI_TempCol=120;
Nded_TempCol=1+transmitter_deadtime(21)*transmitter_sampling_frequency(21);
Ntau_TempCol=1;
xxx_TempCol(1)= (y_ss(22) -yLO_TempCol) / (yHI_TempCol -yLO_TempCol);
for i=2:Nded_TempCol
    xxx_TempCol(i)= xxx_TempCol(1);
end
for i=1:Ntau_TempCol
    yyy_TempCol(i)= xxx_TempCol(Nded_TempCol);
end
xi_TempCol= (MVs(21) - uLO_TempCol) /(uHI_TempCol -uLO_TempCol)-K_TempCol*act_TempCol*((SP_TempCol-yLO_TempCol)/(yHI_TempCol-yLO_TempCol)-yyy_TempCol(Ntau_TempCol));
%-----------------------------------------------------------------------------
%controller initialization
%MV22: Q_Condenser - Decanter Temperature (45.845 degC) is perfetly controlled in the code
SP_DecanterTemp=45.845;
K_DecanterTemp=1;
Ti_DecanterTemp=5;
act_DecanterTemp=-1;
mode_DecanterTemp=1;
Ponly_DecanterTemp=0;
ht_DecanterTemp=1/transmitter_sampling_frequency(22);
hc_DecanterTemp=1/controller_sampling_frequency(22);
uLO_DecanterTemp=0;
uHI_DecanterTemp=150000;
ded_DecanterTemp=transmitter_deadtime(22);
tau_DecanterTemp=transmitter_lag;
yLO_DecanterTemp=40;
yHI_DecanterTemp=50;
Nded_DecanterTemp=1+transmitter_deadtime(22)*transmitter_sampling_frequency(25);
Ntau_DecanterTemp=1;
xxx_DecanterTemp(1)= (y_ss(20) -yLO_DecanterTemp) / (yHI_DecanterTemp -yLO_DecanterTemp);
for i=2:Nded_DecanterTemp
    xxx_DecanterTemp(i)= xxx_DecanterTemp(1);
end
for i=1:Ntau_DecanterTemp
    yyy_DecanterTemp(i)= xxx_DecanterTemp(Nded_DecanterTemp);
end
xi_DecanterTemp= (MVs(22) - uLO_DecanterTemp) /(uHI_DecanterTemp -uLO_DecanterTemp)-K_DecanterTemp*act_DecanterTemp*((SP_DecanterTemp-yLO_DecanterTemp)/(yHI_DecanterTemp-yLO_DecanterTemp)-yyy_DecanterTemp(Ntau_DecanterTemp));
%-----------------------------------------------------------------------------
%controller initialization
%MV22: Organic Level is controlled by Organic Product Flow at 0.5
SP_Organic=0.5;
K_Organic=1;
Ti_Organic=100;
act_Organic=-1;
mode_Organic=1;
Ponly_Organic=1;
ht_Organic=1/transmitter_sampling_frequency(23);
hc_Organic=1/controller_sampling_frequency(23);
uLO_Organic=0;
uHI_Organic=2.4;
ded_Organic=transmitter_deadtime(23);
tau_Organic=transmitter_lag;
yLO_Organic=0;
yHI_Organic=1;
Nded_Organic=1+transmitter_deadtime(23)*transmitter_sampling_frequency(23);
Ntau_Organic=1;
xxx_Organic(1)= (y_ss(18)-yLO_Organic) / (yHI_Organic -yLO_Organic);
for i=2:Nded_Organic
    xxx_Organic(i)= xxx_Organic(1);
end
for i=1:Ntau_Organic
    yyy_Organic(i)= xxx_Organic(Nded_Organic);
end
xi_Organic= (MVs(23) - uLO_Organic) /(uHI_Organic -uLO_Organic)-K_Organic*act_Organic*((SP_Organic-yLO_Organic)/(yHI_Organic-yLO_Organic)-yyy_Organic(Ntau_Organic));
%-----------------------------------------------------------------------------
%controller initialization
%MV23: Aqueous Level is controlled by Aqueous Product Flow at 0.5
SP_Aqueous=0.5;
K_Aqueous=1;
Ti_Aqueous=100;
act_Aqueous=-1;
mode_Aqueous=1;
Ponly_Aqueous=1;
ht_Aqueous=1/transmitter_sampling_frequency(24);
hc_Aqueous=1/controller_sampling_frequency(24);
uLO_Aqueous=0;
uHI_Aqueous=2.4;
ded_Aqueous=transmitter_deadtime(24);
tau_Aqueous=transmitter_lag;
yLO_Aqueous=0;
yHI_Aqueous=1;
Nded_Aqueous=1+transmitter_deadtime(24)*transmitter_sampling_frequency(24);
Ntau_Aqueous=1;
xxx_Aqueous(1)= (y_ss(19)-yLO_Aqueous) / (yHI_Aqueous -yLO_Aqueous);
for i=2:Nded_Aqueous
    xxx_Aqueous(i)= xxx_Aqueous(1);
end
for i=1:Ntau_Aqueous
    yyy_Aqueous(i)= xxx_Aqueous(Nded_Aqueous);
end
xi_Aqueous= (MVs(24) - uLO_Aqueous) /(uHI_Aqueous -uLO_Aqueous)-K_Aqueous*act_Aqueous*((SP_Aqueous-yLO_Aqueous)/(yHI_Aqueous-yLO_Aqueous)-yyy_Aqueous(Ntau_Aqueous));
%-----------------------------------------------------------------------------
%MV24: Column Bottom Exit is used to control the Column Bottom Level
SP_ColButtom=0.5;
K_ColButtom=1;
Ti_ColButtom=100;
act_ColButtom=-1;
mode_ColButtom=1;
Ponly_ColButtom=1;
ht_ColButtom=1/transmitter_sampling_frequency(25);
hc_ColButtom=1/controller_sampling_frequency(25);
uLO_ColButtom=0;
uHI_ColButtom=4.536;
ded_ColButtom=transmitter_deadtime(25);
tau_ColButtom=transmitter_lag;
yLO_ColButtom=0;
yHI_ColButtom=1;
Nded_ColButtom=1+transmitter_deadtime(25)*transmitter_sampling_frequency(25);
Ntau_ColButtom=1;
xxx_ColButtom(1)= (y_ss(21)-yLO_ColButtom) / (yHI_ColButtom -yLO_ColButtom);
for i=2:Nded_ColButtom
    xxx_ColButtom(i)= xxx_ColButtom(1);
end
for i=1:Ntau_ColButtom
    yyy_ColButtom(i)= xxx_ColButtom(Nded_ColButtom);
end
xi_ColButtom= (MVs(25) - uLO_ColButtom) /(uHI_ColButtom -uLO_ColButtom)-K_ColButtom*act_ColButtom*((SP_ColButtom-yLO_ColButtom)/(yHI_ColButtom-yLO_ColButtom)-yyy_ColButtom(Ntau_ColButtom));
%-----------------------------------------------------------------------------
%controller initialization
%MV26: Vaporizer Liquid Inlet is fixed at 2.1924 kmol/min
%-----------------------------------------------------------------------------

%-----------------------------------------------------------------------------
%initialize historical data vectors
x_history=[];
y_history=[];
u_history=[];
dx_history=[];
%initialize simulation parameters
tic
timer=0;
is_initial=0;
disturbance_ID=0;
%start simulation
for k=0:(model_sampling_frequency*minute)
    time=k/model_sampling_frequency; %in minute
    %initialize disturbances
    if (k/model_sampling_frequency)>=10
        switch selected_ID
        case 0
            ;
        case 1
            SP_TRCT=151;
        case 2
            SP_TRCT=165;
        case 3
            SP_H2OCol=0.18;
        case 4
            MVs(26)=2.64;
        case 5
            if (k/model_sampling_frequency)<15
                MVs(3)=0;
            end
        case 6
            if (k/model_sampling_frequency)<15
                MVs(1)=0;
            end
        case 7
             disturbance_ID=1;
        case 8
            if (k/model_sampling_frequency)<15
                disturbance_ID=2;
            else
                disturbance_ID=0;
            end
        end
    end
    %calculate state derivatives and measurements
    [dstatedt,states,MVs,y_ss]=VAModel(states,MVs,time,is_initial,disturbance_ID);
    %in this example, we assume the setpoint of the vaporizer pressure control is given by the current gas recycle pressure
    SP_PresVap=y_ss(12);
    %command window display refreshing rate
    if rem(k,display_sampling_frequency*model_sampling_frequency)==0
        [error, indx]=max(abs(dstatedt./states));
%        info=['minute: ' num2str(time) '    max_dxdt_error: ' num2str(error) '    location: ' num2str(indx)];
 %       disp(info);
    end
    %save information
    if rem(k,storage_sampling_frequency*model_sampling_frequency)==0
        x_history=[x_history;states'];
        y_history=[y_history;[y_ss(37);y_ss(12);y_ss(23);y_ss(2);y_ss(1);y_ss(4);y_ss(5);y_ss(9);y_ss(10);MVs(10);y_ss(11);y_ss(13);MVs(13);y_ss(14);MVs(15);y_ss(15);y_ss(31);y_ss(33);y_ss(8);y_ss(28);y_ss(22);y_ss(20);y_ss(18);y_ss(19);y_ss(21);MVs(26);y_ss(27)]'];
        u_history=[u_history;MVs'];
        dx_history=[dx_history;dstatedt'];
    end
    %update states
    states=states+(1/model_sampling_frequency)*dstatedt;
    %controller action
    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(1))==0 & k~=0
        [xxx_O2, yyy_O2]=Transmit(y_ss(37), xxx_O2, yyy_O2, ded_O2, tau_O2, Nded_O2, Ntau_O2, ht_O2, yLO_O2, yHI_O2);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(1))==0 & k~=0
        [MVs(1),xi_O2,flag_O2]=controller(SP_O2, yyy_O2(Ntau_O2), MVs(1), xi_O2, mode_O2, K_O2, Ti_O2, hc_O2,...
            yLO_O2, yHI_O2, uLO_O2, uHI_O2, act_O2, Ponly_O2, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(2))==0 & k~=0
        [xxx_C2H4, yyy_C2H4]=Transmit(y_ss(12), xxx_C2H4, yyy_C2H4, ded_C2H4, tau_C2H4, Nded_C2H4, Ntau_C2H4, ht_C2H4, yLO_C2H4, yHI_C2H4);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(2))==0 & k~=0
        [MVs(2),xi_C2H4,flag_C2H4]=controller(SP_C2H4, yyy_C2H4(Ntau_C2H4), MVs(2), xi_C2H4, mode_C2H4, K_C2H4, Ti_C2H4, hc_C2H4,...
            yLO_C2H4, yHI_C2H4, uLO_C2H4, uHI_C2H4, act_C2H4, Ponly_C2H4, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(3))==0 & k~=0
        [xxx_HAc, yyy_HAc]=Transmit(y_ss(23), xxx_HAc, yyy_HAc, ded_HAc, tau_HAc, Nded_HAc, Ntau_HAc, ht_HAc, yLO_HAc, yHI_HAc);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(3))==0 & k~=0
        [MVs(3),xi_HAc,flag_HAc]=controller(SP_HAc, yyy_HAc(Ntau_HAc), MVs(3), xi_HAc, mode_HAc, K_HAc, Ti_HAc, hc_HAc,...
            yLO_HAc, yHI_HAc, uLO_HAc, uHI_HAc, act_HAc, Ponly_HAc, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(4))==0 & k~=0
        [xxx_LevelVap, yyy_LevelVap]=Transmit(y_ss(2), xxx_LevelVap, yyy_LevelVap, ded_LevelVap, tau_LevelVap, Nded_LevelVap, Ntau_LevelVap, ht_LevelVap, yLO_LevelVap, yHI_LevelVap);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(4))==0 & k~=0
        [MVs(4),xi_LevelVap,flag_LevelVap]=controller(SP_LevelVap, yyy_LevelVap(Ntau_LevelVap), MVs(4), xi_LevelVap, mode_LevelVap, K_LevelVap, Ti_LevelVap, hc_LevelVap,...
            yLO_LevelVap, yHI_LevelVap, uLO_LevelVap, uHI_LevelVap, act_LevelVap, Ponly_LevelVap, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(5))==0 & k~=0
        [xxx_PresVap, yyy_PresVap]=Transmit(y_ss(1), xxx_PresVap, yyy_PresVap, ded_PresVap, tau_PresVap, Nded_PresVap, Ntau_PresVap, ht_PresVap, yLO_PresVap, yHI_PresVap);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(5))==0 & k~=0
        [MVs(5),xi_PresVap,flag_PresVap]=controller(SP_PresVap, yyy_PresVap(Ntau_PresVap), MVs(5), xi_PresVap, mode_PresVap, K_PresVap, Ti_PresVap, hc_PresVap,...
            yLO_PresVap, yHI_PresVap, uLO_PresVap, uHI_PresVap, act_PresVap, Ponly_PresVap, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(6))==0 & k~=0
        [xxx_ReactorInletTemp, yyy_ReactorInletTemp]=Transmit(y_ss(4), xxx_ReactorInletTemp, yyy_ReactorInletTemp, ded_ReactorInletTemp, tau_ReactorInletTemp, Nded_ReactorInletTemp, Ntau_ReactorInletTemp, ht_ReactorInletTemp, yLO_ReactorInletTemp, yHI_ReactorInletTemp);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(6))==0 & k~=0
        [MVs(6),xi_ReactorInletTemp,flag_ReactorInletTemp]=controller(SP_ReactorInletTemp, yyy_ReactorInletTemp(Ntau_ReactorInletTemp), MVs(6), xi_ReactorInletTemp, mode_ReactorInletTemp, K_ReactorInletTemp, Ti_ReactorInletTemp, hc_ReactorInletTemp,...
            yLO_ReactorInletTemp, yHI_ReactorInletTemp, uLO_ReactorInletTemp, uHI_ReactorInletTemp, act_ReactorInletTemp, Ponly_ReactorInletTemp, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(7))==0 & k~=0
        [xxx_TRCT, yyy_TRCT]=Transmit(y_ss(5), xxx_TRCT, yyy_TRCT, ded_TRCT, tau_TRCT, Nded_TRCT, Ntau_TRCT, ht_TRCT, yLO_TRCT, yHI_TRCT);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(7))==0 & k~=0
        [MVs(7),xi_TRCT,flag_TRCT]=controller(SP_TRCT, yyy_TRCT(Ntau_TRCT), MVs(7), xi_TRCT, mode_TRCT, K_TRCT, Ti_TRCT, hc_TRCT,...
            yLO_TRCT, yHI_TRCT, uLO_TRCT, uHI_TRCT, act_TRCT, Ponly_TRCT, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(8))==0 & k~=0
        [xxx_LevelSep, yyy_LevelSep]=Transmit(y_ss(9), xxx_LevelSep, yyy_LevelSep, ded_LevelSep, tau_LevelSep, Nded_LevelSep, Ntau_LevelSep, ht_LevelSep, yLO_LevelSep, yHI_LevelSep);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(8))==0 & k~=0
        [MVs(8),xi_LevelSep,flag_LevelSep]=controller(SP_LevelSep, yyy_LevelSep(Ntau_LevelSep), MVs(8), xi_LevelSep, mode_LevelSep, K_LevelSep, Ti_LevelSep, hc_LevelSep,...
            yLO_LevelSep, yHI_LevelSep, uLO_LevelSep, uHI_LevelSep, act_LevelSep, Ponly_LevelSep, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(9))==0 & k~=0
        [xxx_TempSep, yyy_TempSep]=Transmit(y_ss(10), xxx_TempSep, yyy_TempSep, ded_TempSep, tau_TempSep, Nded_TempSep, Ntau_TempSep, ht_TempSep, yLO_TempSep, yHI_TempSep);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(9))==0 & k~=0
        [MVs(9),xi_TempSep,flag_TempSep]=controller(SP_TempSep, yyy_TempSep(Ntau_TempSep), MVs(9), xi_TempSep, mode_TempSep, K_TempSep, Ti_TempSep, hc_TempSep,...
            yLO_TempSep, yHI_TempSep, uLO_TempSep, uHI_TempSep, act_TempSep, Ponly_TempSep, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(11))==0 & k~=0
        [xxx_CompOutletTemp, yyy_CompOutletTemp]=Transmit(y_ss(11), xxx_CompOutletTemp, yyy_CompOutletTemp, ded_CompOutletTemp, tau_CompOutletTemp, Nded_CompOutletTemp, Ntau_CompOutletTemp, ht_CompOutletTemp, yLO_CompOutletTemp, yHI_CompOutletTemp);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(11))==0 & k~=0
        [MVs(11),xi_CompOutletTemp,flag_CompOutletTemp]=controller(SP_CompOutletTemp, yyy_CompOutletTemp(Ntau_CompOutletTemp), MVs(11), xi_CompOutletTemp, mode_CompOutletTemp, K_CompOutletTemp, Ti_CompOutletTemp, hc_CompOutletTemp,...
            yLO_CompOutletTemp, yHI_CompOutletTemp, uLO_CompOutletTemp, uHI_CompOutletTemp, act_CompOutletTemp, Ponly_CompOutletTemp, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(12))==0 & k~=0
        [xxx_LevelAbs, yyy_LevelAbs]=Transmit(y_ss(13), xxx_LevelAbs, yyy_LevelAbs, ded_LevelAbs, tau_LevelAbs, Nded_LevelAbs, Ntau_LevelAbs, ht_LevelAbs, yLO_LevelAbs, yHI_LevelAbs);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(12))==0 & k~=0
        [MVs(12),xi_LevelAbs,flag_LevelAbs]=controller(SP_LevelAbs, yyy_LevelAbs(Ntau_LevelAbs), MVs(12), xi_LevelAbs, mode_LevelAbs, K_LevelAbs, Ti_LevelAbs, hc_LevelAbs,...
            yLO_LevelAbs, yHI_LevelAbs, uLO_LevelAbs, uHI_LevelAbs, act_LevelAbs, Ponly_LevelAbs, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(14))==0 & k~=0
        [xxx_CircOutletTemp, yyy_CircOutletTemp]=Transmit(y_ss(14), xxx_CircOutletTemp, yyy_CircOutletTemp, ded_CircOutletTemp, tau_CircOutletTemp, Nded_CircOutletTemp, Ntau_CircOutletTemp, ht_CircOutletTemp, yLO_CircOutletTemp, yHI_CircOutletTemp);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(14))==0 & k~=0
        [MVs(14),xi_CircOutletTemp,flag_CircOutletTemp]=controller(SP_CircOutletTemp, yyy_CircOutletTemp(Ntau_CircOutletTemp), MVs(14), xi_CircOutletTemp, mode_CircOutletTemp, K_CircOutletTemp, Ti_CircOutletTemp, hc_CircOutletTemp,...
            yLO_CircOutletTemp, yHI_CircOutletTemp, uLO_CircOutletTemp, uHI_CircOutletTemp, act_CircOutletTemp, Ponly_CircOutletTemp, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(16))==0 & k~=0
        [xxx_ScrubOutletTemp, yyy_ScrubOutletTemp]=Transmit(y_ss(15), xxx_ScrubOutletTemp, yyy_ScrubOutletTemp, ded_ScrubOutletTemp, tau_ScrubOutletTemp, Nded_ScrubOutletTemp, Ntau_ScrubOutletTemp, ht_ScrubOutletTemp, yLO_ScrubOutletTemp, yHI_ScrubOutletTemp);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(16))==0 & k~=0
        [MVs(16),xi_ScrubOutletTemp,flag_ScrubOutletTemp]=controller(SP_ScrubOutletTemp, yyy_ScrubOutletTemp(Ntau_ScrubOutletTemp), MVs(16), xi_ScrubOutletTemp, mode_ScrubOutletTemp, K_ScrubOutletTemp, Ti_ScrubOutletTemp, hc_ScrubOutletTemp,...
            yLO_ScrubOutletTemp, yHI_ScrubOutletTemp, uLO_ScrubOutletTemp, uHI_ScrubOutletTemp, act_ScrubOutletTemp, Ponly_ScrubOutletTemp, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(17))==0 & k~=0
        [xxx_CO2, yyy_CO2]=Transmit(y_ss(31), xxx_CO2, yyy_CO2, ded_CO2, tau_CO2, Nded_CO2, Ntau_CO2, ht_CO2, yLO_CO2, yHI_CO2);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(17))==0 & k~=0
        [MVs(17),xi_CO2,flag_CO2]=controller(SP_CO2, yyy_CO2(Ntau_CO2), MVs(17), xi_CO2, mode_CO2, K_CO2, Ti_CO2, hc_CO2,...
            yLO_CO2, yHI_CO2, uLO_CO2, uHI_CO2, act_CO2, Ponly_CO2, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(18))==0 & k~=0
        [xxx_C2H6, yyy_C2H6]=Transmit(y_ss(33), xxx_C2H6, yyy_C2H6, ded_C2H6, tau_C2H6, Nded_C2H6, Ntau_C2H6, ht_C2H6, yLO_C2H6, yHI_C2H6);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(18))==0 & k~=0
        [MVs(18),xi_C2H6,flag_C2H6]=controller(SP_C2H6, yyy_C2H6(Ntau_C2H6), MVs(18), xi_C2H6, mode_C2H6, K_C2H6, Ti_C2H6, hc_C2H6,...
            yLO_C2H6, yHI_C2H6, uLO_C2H6, uHI_C2H6, act_C2H6, Ponly_C2H6, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(19))==0 & k~=0
        [xxx_FEHE, yyy_FEHE]=Transmit(y_ss(8), xxx_FEHE, yyy_FEHE, ded_FEHE, tau_FEHE, Nded_FEHE, Ntau_FEHE, ht_FEHE, yLO_FEHE, yHI_FEHE);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(19))==0 & k~=0
        [MVs(19),xi_FEHE,flag_FEHE]=controller(SP_FEHE, yyy_FEHE(Ntau_FEHE), MVs(19), xi_FEHE, mode_FEHE, K_FEHE, Ti_FEHE, hc_FEHE,...
            yLO_FEHE, yHI_FEHE, uLO_FEHE, uHI_FEHE, act_FEHE, Ponly_FEHE, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(20))==0 & k~=0
        [xxx_H2OCol, yyy_H2OCol]=Transmit(y_ss(28), xxx_H2OCol, yyy_H2OCol, ded_H2OCol, tau_H2OCol, Nded_H2OCol, Ntau_H2OCol, ht_H2OCol, yLO_H2OCol, yHI_H2OCol);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(20))==0 & k~=0
        [MVs(20),xi_H2OCol,flag_H2OCol]=controller(SP_H2OCol, yyy_H2OCol(Ntau_H2OCol), MVs(20), xi_H2OCol, mode_H2OCol, K_H2OCol, Ti_H2OCol, hc_H2OCol,...
            yLO_H2OCol, yHI_H2OCol, uLO_H2OCol, uHI_H2OCol, act_H2OCol, Ponly_H2OCol, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(21))==0 & k~=0
        [xxx_TempCol, yyy_TempCol]=Transmit(y_ss(22), xxx_TempCol, yyy_TempCol, ded_TempCol, tau_TempCol, Nded_TempCol, Ntau_TempCol, ht_TempCol, yLO_TempCol, yHI_TempCol);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(21))==0 & k~=0
        [MVs(21),xi_TempCol,flag_TempCol]=controller(SP_TempCol, yyy_TempCol(Ntau_TempCol), MVs(21), xi_TempCol, mode_TempCol, K_TempCol, Ti_TempCol, hc_TempCol,...
            yLO_TempCol, yHI_TempCol, uLO_TempCol, uHI_TempCol, act_TempCol, Ponly_TempCol, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(22))==0 & k~=0
        [xxx_DecanterTemp, yyy_DecanterTemp]=Transmit(y_ss(20), xxx_DecanterTemp, yyy_DecanterTemp, ded_DecanterTemp, tau_DecanterTemp, Nded_DecanterTemp, Ntau_DecanterTemp, ht_DecanterTemp, yLO_DecanterTemp, yHI_DecanterTemp);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(22))==0 & k~=0
        [MVs(22),xi_DecanterTemp,flag_DecanterTemp]=controller(SP_DecanterTemp, yyy_DecanterTemp(Ntau_DecanterTemp), MVs(22), xi_DecanterTemp, mode_DecanterTemp, K_DecanterTemp, Ti_DecanterTemp, hc_DecanterTemp,...
            yLO_DecanterTemp, yHI_DecanterTemp, uLO_DecanterTemp, uHI_DecanterTemp, act_DecanterTemp, Ponly_DecanterTemp, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(23))==0 & k~=0
        [xxx_Organic, yyy_Organic]=Transmit(y_ss(18), xxx_Organic, yyy_Organic, ded_Organic, tau_Organic, Nded_Organic, Ntau_Organic, ht_Organic, yLO_Organic, yHI_Organic);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(23))==0 & k~=0
        [MVs(23),xi_Organic,flag_Organic]=controller(SP_Organic, yyy_Organic(Ntau_Organic), MVs(23), xi_Organic, mode_Organic, K_Organic, Ti_Organic, hc_Organic,...
            yLO_Organic, yHI_Organic, uLO_Organic, uHI_Organic, act_Organic, Ponly_Organic, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(24))==0 & k~=0
        [xxx_Aqueous, yyy_Aqueous]=Transmit(y_ss(19), xxx_Aqueous, yyy_Aqueous, ded_Aqueous, tau_Aqueous, Nded_Aqueous, Ntau_Aqueous, ht_Aqueous, yLO_Aqueous, yHI_Aqueous);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(24))==0 & k~=0
        [MVs(24),xi_Aqueous,flag_Aqueous]=controller(SP_Aqueous, yyy_Aqueous(Ntau_Aqueous), MVs(24), xi_Aqueous, mode_Aqueous, K_Aqueous, Ti_Aqueous, hc_Aqueous,...
            yLO_Aqueous, yHI_Aqueous, uLO_Aqueous, uHI_Aqueous, act_Aqueous, Ponly_Aqueous, 0, 0, 0, 0);
    end

    if rem(k,model_sampling_frequency/transmitter_sampling_frequency(25))==0 & k~=0
        [xxx_ColButtom, yyy_ColButtom]=Transmit(y_ss(21), xxx_ColButtom, yyy_ColButtom, ded_ColButtom, tau_ColButtom, Nded_ColButtom, Ntau_ColButtom, ht_ColButtom, yLO_ColButtom, yHI_ColButtom);
    end
    if rem(k,model_sampling_frequency/controller_sampling_frequency(25))==0 & k~=0
        [MVs(25),xi_ColButtom,flag_ColButtom]=controller(SP_ColButtom, yyy_ColButtom(Ntau_ColButtom), MVs(25), xi_ColButtom, mode_ColButtom, K_ColButtom, Ti_ColButtom, hc_ColButtom,...
            yLO_ColButtom, yHI_ColButtom, uLO_ColButtom, uHI_ColButtom, act_ColButtom, Ponly_ColButtom, 0, 0, 0, 0);
    end

end
toc
%plot graphics
my_label=['   %O2  ';'  Pres  ';'  HAc-L ';'  Vap-L ';'  Vap-P ';'  Pre-T ';'  RCT-T ';'  Sep-L ';'  Sep-T ';'  Sep-V ';'  Com-T ';'  Abs-L ';'  Cir-F ';'  Cir-T ';'  Scr-F ';'  Scr-T ';'  %CO2  ';'  %C2H6 ';' FEHE-T ';'  %H2O  ';'  Col-T ';'  Org-L ';'  Aqu-L ';'  Col-L ';' Vap-In ';' Dect-T ';'%VAc E-3'];
MV_label=[' F-O2 ';'F-C2H4';' F-HAc';' Q-Vap';' F-Vap';'Q-Heat';'ShellT';'F-SepL';' T-Sep';'F-SepV';'Q-Comp';'F-AbsL';'F-Circ';'Q-Circ';'F-Scru';'Q-Scru';' F-CO2';' Purge';'bypass';'Reflux';'Q-Rebo';'F-Orga';'F-Aque';' F-Bot';'Q_Cond';'F-Tank'];
setpoint=[SP_O2;SP_C2H4;SP_HAc;SP_LevelVap;SP_PresVap;150;SP_TRCT;SP_LevelSep;SP_TempSep;16.1026;80;SP_LevelAbs;15.1198;25;0.756;25;SP_CO2;SP_C2H6;SP_FEHE;SP_H2OCol;SP_TempCol;0.5;0.5;SP_ColButtom;45.845;2.1924];
warning off
if selected_ID==1
    McAvoy_Plot_1(y_history,u_history,setpoint,my_label,MV_label,storage_sampling_frequency);
elseif selected_ID==2
    McAvoy_Plot_2(y_history,u_history,setpoint,my_label,MV_label,storage_sampling_frequency);
elseif selected_ID==3
    McAvoy_Plot_3(y_history,u_history,setpoint,my_label,MV_label,storage_sampling_frequency);
elseif selected_ID==4
    McAvoy_Plot_4(y_history,u_history,setpoint,my_label,MV_label,storage_sampling_frequency);
else
    Transient_Plot(y_history,u_history,setpoint,my_label,MV_label,storage_sampling_frequency);
end
return


text:
function basis_11_t6_test ( )

%*****************************************************************************80
%
%% BASIS_11_T6_TEST verifies BASIS_11_T6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  node_num = 6;

  t = [ ...
    2.0, 0.0; ...
    4.0, 3.0; ...
    0.0, 4.0; ...
    3.0, 1.5; ...
    2.0, 3.5; ...
    1.0, 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BASIS_11_T6_TEST:\n' );
  fprintf ( 1, '  Verify basis functions for element T6.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Physical Nodes:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_num
    fprintf ( 1, '  %8d  %7f  %7f\n', j, t(1:2,j) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The basis function values at basis nodes\n' );
  fprintf ( 1, '  should form the identity matrix.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : node_num
    for j = 1 : node_num
      [ phi(i,j), dphidx(i,j), dphidy(i,j) ] = basis_11_t6 ( t, i, t(1:2,j) );
    end
  end

  for i = 1 : node_num
    for ( j = 1 : node_num )
      fprintf ( 1, '  %7f', phi(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The X and Y derivatives should sum to 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      dPhidX sum    dPhidY sum\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_num
    sum_x = sum ( dphidx(1:node_num,j) );
    sum_y = sum ( dphidy(1:node_num,j) );
    fprintf ( 1, '  %14f  %14f\n', sum_x, sum_y );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function B = transform_image_V1_V3(I, corners, interporation)

% B = transform_image3(I, corners, interporation)
% calculate the transformation matrix for each area

B = nan(256,256);

if notDefined('interporation'), interporation = [];  end;

% V3v
input_points = corners{4};
% base_points =  [40 20;40 35; 80 35; 80 20];
base_points =  [30 80; 15 80; 15 40; 30 40];
udata = [1 201];  vdata = [1 201];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',256);
B(round(20/200*256):round(100/200*256), round(10/200*256):round(30/200*256))=tmp(round(20/200*256):round(100/200*256),round(10/200*256):round(30/200*256));

% V2v
input_points = corners{2};
% base_points =  [40 50;40 35; 80 35; 80 50];
base_points =  [30 80; 45 80; 45 40; 30 40];
udata = [1 201];  vdata = [1 201];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',256);
B(round(20/200*256):round(100/200*256), round(30/200*256):round(45/200*256))=tmp(round(20/200*256):round(100/200*256),round(30/200*256):round(45/200*256));

% V1
input_points = corners{1};
% base_points =  [40 50;40 65; 80 65; 80 50];
base_points =  [75 80; 45 80; 45 40; 75 40];
udata = [1 201];  vdata = [1 201];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',256);
B(round(20/200*256):round(100/200*256), round(45/200*256):round(75/200*256))=tmp(round(20/200*256):round(100/200*256),round(45/200*256):round(75/200*256));

% V2d
input_points = corners{3};
% base_points =  [40 80;40 65; 80 65; 80 80];
base_points =  [75 80; 90 80; 90 40; 75 40];
% input_points = [corners{2};mean(corners{2}(2:3,:)); mean(corners{2})];
% % base_points =  [40 80;40 65; 80 65; 80 80];
% base_points =  [90 80; 75 80; 75 40; 90 40; 75 60; 82.5 60];
udata = [1 201];  vdata = [1 201];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',256);
B(round(20/200*256):round(100/200*256), round(75/200*256):round(90/200*256))=tmp(round(20/200*256):round(100/200*256),round(75/200*256):round(90/200*256));

% V3d
input_points = corners{5};
% base_points =  [40 80;40 65; 80 65; 80 80];
base_points =  [105 80; 90 80; 90 40; 105 40];
% input_points = [corners{2};mean(corners{2}(2:3,:)); mean(corners{2})];
% % base_points =  [40 80;40 65; 80 65; 80 80];
% base_points =  [90 80; 75 80; 75 40; 90 40; 75 60; 82.5 60];
udata = [1 201];  vdata = [1 201];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',256);
B(round(20/200*256):round(100/200*256), round(90/200*256):round(110/200*256))=tmp(round(20/200*256):round(100/200*256),round(90/200*256):round(110/200*256));


text:
% std_readtopo() - returns the scalp map of a specified ICA component, assumed
%                  to have been saved in a Matlab file, [dataset_name].icatopo,
%                  in the same directory as the dataset file. If this file does
%                  not exist, use std_topo() to create it, else a pre-clustering
%                  function that calls it: pop_preclust() or eeg_preclust().
% Usage:
%   >> [grid, y, x ] = std_readtopo(ALLEEG, setindx, component);
%   >> [grid, y, x ] = std_readtopo(ALLEEG, setindx, component, transform, mode);
%
% Inputs:
%   ALLEEG     - vector of EEG datasets (can also be one EEG set).
%                must contain the dataset of interest (see 'setindx' below).
%   setindx    - [integer] an index of an EEG dataset in the ALLEEG
%                structure, for which to get the component ERP.
%   component  - [integer] index of the component for which the scalp map
%                grid should be returned.
%   transform  - ['none'!'laplacian'|'gradient'] transform scalp map to
%                laplacian or gradient map. Default is 'none'.
%   mode       - ['2dmap'|'preclust'] return either a 2-D array for direct
%                plotting ('2dmap') or an array formated for preclustering
%                with all the NaN values removed (ncomps x points). Default
%                is '2dmap' for 1 component and 'preclust' for several.
%
% Outputs:
%   grid      - square scalp-map color-value grid for the requested ICA component
%               in the specified dataset, an interpolated Cartesian grid as output
%               by topoplot().
%   y         - y-axis values for the interpolated grid
%   x         - x-axis values of the interpolated grid
%
%  See also  std_topo(), std_preclust()
%
% Authors: Arnaud Delorme, Hilit Serby, SCCN, INC, UCSD, February, 2005

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

function [X, yi, xi ] = std_readtopo(ALLEEG, abset, comps, option, mode)

X = [];
yi = [];
xi = [];
if nargin < 4
    option = 'none';
end;
if nargin < 5
    mode = '2Dmap';
end;
filename = correctfile(fullfile( ALLEEG(abset).filepath,[ ALLEEG(abset).filename(1:end-3) 'icatopo']));
tmpfile  = which(filename);
if ~isempty(tmpfile), filename = tmpfile; end;

while getfield(dir(filename), 'bytes') < 1000
    topo = load( '-mat', filename);
    filename = correctfile(topo.file);
    tmpfile  = which(filename);
    if ~isempty(tmpfile), filename = tmpfile; end;
end;

for k = 1:length(comps)

    if length(comps) < 3
        warning off;
        try
            topo = load( '-mat', filename, ...
                         [ 'comp' int2str(comps(k)) '_grid'], ...
                         [ 'comp' int2str(comps(k)) '_x'], ...
                         [ 'comp' int2str(comps(k)) '_y'] );
        catch
            error( [ 'Cannot read file ''' filename '''' ]);
        end;
        warning on;
    elseif k == 1
        try
            topo = load( '-mat', filename);
        catch
            error( [ 'Cannot read file ''' filename '''' ]);
        end;
    end;

    tmp =  getfield(topo, [ 'comp' int2str(comps(k)) '_grid' ]);

    if strcmpi(option, 'gradient')
        [tmpx, tmpy]  = gradient(tmp); % Gradient
        tmp        = tmpx;
        tmp(:,:,2) = tmpy;
    elseif strcmpi(option, 'laplacian')
        tmp = del2(tmp); % Laplacian
    end;

    if length(comps) > 1 | strcmpi(mode, 'preclust')
        tmp = tmp(find(~isnan(tmp))); % remove NaN for more than 1 component
    end;
    if k == 1
        X = zeros([ length(comps) size(tmp) ]) ;
    end
    X(k,:,:,:) =  tmp;
    if k == 1
        yi   = getfield(topo, [ 'comp' int2str(comps(k)) '_y']);
        xi   = getfield(topo, [ 'comp' int2str(comps(k)) '_x']);
    end;
end
X = squeeze(X);

return;

function filename = correctfile(filename)
    comp = computer;
    if filename(2) == ':' & ~strcmpi(comp(1:2), 'PC')
        filename = [filesep filename(4:end) ];
        filename(find(filename == '\')) = filesep;
    end;

    if ~exist(filename)
        [tmpp tmpf ext] = fileparts(filename);
        if exist([tmpf ext])
            filename = [tmpf ext];
        else
            [tmpp2 tmpp1] = fileparts(tmpp);
            if exist(fullfile(tmpp1, [ tmpf ext ]))
                filename = fullfile(tmpp1, [ tmpf ext ]);
            else
                error([ 'Cannot load file ''' [ tmpf ext ] '''' ]);
            end;
        end;
    end;

author: PatternRecognition

repo: OpenBMI


text:
%Demonstrates compressively sampling and regularization by denoising (RED) recovery of an image.
%Requires  matconvnet and fasta in the path
%Google's RED implementation can be found here: https://github.com/google/RED
% RED reference: Romano, Yaniv, Michael Elad, and Peyman Milanfar. "The little engine that could: Regularization by denoising (RED)." SIAM Journal on Imaging Sciences 10.4 (2017): 1804-1844.

addpath(genpath('..'))
addpath('~/matconvnet/matlab');
addpath(genpath('~/fasta-matlab'));
% addpath(genpath('~/gampmatlab'));

denoiser1='DnCNN';%Available options are NLM, Gauss, Bilateral, BLS-GSM, BM3D, fast-BM3D, BM3D-SAPCA, and DnCNN
denoiser2='DnCNN';
filename='boat.png';
SamplingRate=.1;
RED_iters=200;
AMP_iters=10;
% VAMP_iters=10;
imsize=128;
n_DnCNN_layers=20;%Other option is 17
LoadNetworkWeights(n_DnCNN_layers);

ImIn=double(imread(filename));
x_0=imresize(ImIn,imsize/size(ImIn,1));
[height, width]=size(x_0);
n=length(x_0(:));
m=round(n*SamplingRate);
errfxn = @(x_hat) PSNR(x_0,reshape(x_hat,[height width]));

% %Generate Gaussian Measurement Matrix
% M_matrix=1/sqrt(m)*randn(m,n);
% M=@(x) M_matrix*x(:);
% Mt=@(x) M_matrix'*x(:);
% U=[];
% Ut=[];
% d=[];

% %Generate Coded Diffraction Pattern Measurement Matrix
% signvec = exp(1i*2*pi*rand(n,1));
% inds=[1;randsample(n-1,m-1)+1];
% I=speye(n);
% SubsampM=I(inds,:);
% M=@(x) SubsampM*reshape(fft2(reshape(bsxfun(@times,signvec,x(:)),[height,width])),[n,1])*(1/sqrt(n))*sqrt(n/m);
% Mt=@(x) bsxfun(@times,conj(signvec),reshape(ifft2(reshape(SubsampM'*x(:),[height,width])),[n,1]))*sqrt(n)*sqrt(n/m);
% U=@(x) x(:);
% Ut= @(x) x(:);
% d=ones(m,1)*n/m;

% %Generate Real-valued Measurement Matrix with Fast Transformation and
% %approximately i.i.d. Gaussian distribution
% signvec = 2*round(rand(n,1))-1;
% inds=[1;randsample(n-1,m-1)+1];
% I=speye(n);
% SubsampM=I(inds,:);
% M=@(x) SubsampM*reshape(dct2(reshape(bsxfun(@times,signvec,x(:)),[height,width])),[n,1])*sqrt(n/m);
% Mt=@(x) bsxfun(@times,conj(signvec),reshape(idct2(reshape(SubsampM'*x(:),[height,width])),[n,1]))*sqrt(n/m);
% U=@(x) x(:);
% Ut= @(x) x(:);
% d=ones(m,1)*n/m;

%Generate (something close to) a Fast JL Transform matrix
signvec = 2*round(rand(n,1))-1;
inds=[1;randsample(n-1,m-1)+1];
I=speye(n);
SubsampM=I(inds,:);
M=@(x) SubsampM*reshape(dct(bsxfun(@times,signvec,x(:))),[n,1])*sqrt(n/m);
Mt=@(x) bsxfun(@times,conj(signvec),reshape(idct(SubsampM'*x(:)),[n,1]))*sqrt(n/m);
U=@(x) x(:);
Ut= @(x) x(:);
d=ones(m,1)*n/m;

% %Generate (something close to) a Fast JL Transform matrix
% signvec = 2*round(rand(n,1))-1;
% inds=[1;randsample(n-1,m-1)+1];
% I=speye(n);
% SubsampM=I(inds,:);
% M=@(x) SubsampM*reshape(fwht(bsxfun(@times,signvec,x(:))),[n,1])*n/sqrt(m);
% Mt=@(x) bsxfun(@times,conj(signvec),reshape(ifwht(SubsampM'*x(:)),[n,1]))*sqrt(1/m);
% U=@(x) x(:);
% Ut= @(x) x(:);
% d=ones(m,1)*n/m;

%Compressively sample the image
y=M(x_0(:));

%Recover signal using RED
%Set RED options
prox_opts=[];
prox_opts.width=width;
prox_opts.height=height;
prox_opts.denoiser=denoiser1;
prox_opts.prox_iters=1;
fasta_opts=[];
fasta_opts.maxIters=RED_iters/2;
fasta_opts.tol=1e-7;
fasta_opts.recordObjective=true;
% fasta_ops.function=errfxn;

x_init=100*rand(size(x_0(:)));

%Recover Signal using RED
prox_opts.lambda=1;
t0=tic;
x_hat1=x_init;
prox_opts.sigma_hat=50;
[x_hat1,out1]  = redCS( M,Mt,y,x_hat1(:),fasta_opts,prox_opts);
prox_opts.sigma_hat=15;
[x_hat1,out2]  = redCS( M,Mt,y,x_hat1(:),fasta_opts,prox_opts);
x_hat1=reshape(x_hat1,[height, width]);
t1=toc(t0);

% figure(4);
% subplot(1,2,1);plot(out1.objective);
% subplot(1,2,2);plot(out2.objective);

%Recover Signal using D-(V)AMP algorithms
t0=tic;[x_hat2,psnr2]  = DAMP(y,AMP_iters,height,width,denoiser1,M,Mt,errfxn);t2=toc(t0);
% t0=tic;[x_hat2,psnr2] = DVAMP(y,VAMP_iters,height,width,denoiser2,M,Mt,errfxn, U, Ut, d);t2=toc(t0);

%D(V)AMP Recovery Performance
performance1=PSNR(x_0,x_hat1);
performance2=PSNR(x_0,x_hat2);
display([num2str(SamplingRate*100),'% Sampling ', denoiser1, '-RED: PSNR=',num2str(performance1),', time=',num2str(t1)])
display([num2str(SamplingRate*100),'% Sampling ', denoiser2, '-AMP: PSNR=',num2str(performance2),', time=',num2str(t2)])

%Plot Recovered Signals
figure(1); clf;
subplot(1,3,1);
imshow(uint8(x_0));title('Original Image');
subplot(1,3,2);
imshow(uint8(x_hat1));title([denoiser1, '-RED']);
subplot(1,3,3);
imshow(uint8(x_hat2));title([denoiser2, '-AMP']);

author: ricedsp


text:
function ui = mrViewEditCbar(ui, o);
%
% ui = mrViewEditCbar(ui, o);
%
% Edit the colorbar for a given overlay in a mrViewer UI.
%
% ras, 08/2006.
if ~exist('ui','var') | isempty(ui),    ui = mrViewGet; end
if ishandle(ui), ui = get(ui, 'UserData'); end

if ~exist('o','var') | isempty(o)
    o = length(ui.overlays);
end

ui.overlays(o).cbar = cbarEdit(ui.overlays(o).cbar);

% if isfield(ui, 'fig') & ishandle(ui.fig)
%     set(ui.fig, 'UserData', ui);
% end

% update the cbar
mrViewRefresh(ui);

return


text:
function saveResults(destDir, canonicalImageSize)
% alignment results
load(fullfile(destDir, 'final.mat'),'A','E','O','M') ;

%% display
mkdir([destDir '\Err']);
mkdir([destDir '\Res']);
mkdir([destDir '\O']);
FACT = 255;
for i=1:size(A,2)
    im = reshape(A(:,i), canonicalImageSize);
    imMin = min(im(:));
    imMax = max(im(:));
    im = (im - imMin)./(imMax - imMin);
    im = uint8(255*im);
    imwrite(im,[destDir '\Res\' sprintf('Image_%.4d.bmp',i)],'bmp');

    % E
    im = reshape(E(:,i), canonicalImageSize);
    im = abs(im) ./ max(max(abs(im)));
    im = uint8(FACT.*abs(im));
    imwrite(im,[destDir '\Err\' sprintf('Er_%.4d.bmp',i)],'bmp');

    % O
    im = reshape(O(:,i), canonicalImageSize);
    im = abs(im) ./ max(max(abs(im)));
    im = uint8(FACT.*abs(im));
    imwrite(im,[destDir '\O\' sprintf('O_%.4d.bmp',i)],'bmp');
end;

author: andrewssobral

repo: lrslibrary


text:
%==================================================================
%                        General Data File
% Title: Default_title
% Units: SI
% Dimensions: 2D
% Type of problem: Plane_Stress
% Type of Phisics: Stokes
% Micro/Macro: MACRO
%
%==================================================================

%% Data

Data_prb = {
'TRIANGLE';
'SI';
'2D';
'Plane_Stress';
'Stokes';
'MACRO';
};

%% Coordinates
% Node                X                Y                Z

coord = [
1            0            0            0
2            0         0.25            0
3         0.25            0            0
4         0.25         0.25            0
5            0          0.5            0
6          0.5            0            0
7         0.25          0.5            0
8          0.5         0.25            0
9          0.5          0.5            0
10         0.75            0            0
11            0         0.75            0
12         0.75         0.25            0
13         0.25         0.75            0
14         0.75          0.5            0
15          0.5         0.75            0
16            1            0            0
17            0            1            0
18            1         0.25            0
19         0.25            1            0
20         0.75         0.75            0
21            1          0.5            0
22          0.5            1            0
23         0.75            1            0
24            1         0.75            0
25            1            1            0
];

%% Conectivities
% Element        Node(1)                Node(2)                Node(3)                Material

connec = [
1 6 8 3 0
2 8 9 4 0
3 3 4 1 0
4 8 4 3 0
5 9 7 4 0
6 7 5 2 0
7 4 2 1 0
8 7 2 4 0
9 16 18 10 0
10 18 21 12 0
11 10 12 6 0
12 18 12 10 0
13 21 14 12 0
14 14 9 8 0
15 12 8 6 0
16 14 8 12 0
17 9 15 7 0
18 15 22 13 0
19 7 13 5 0
20 15 13 7 0
21 22 19 13 0
22 19 17 11 0
23 13 11 5 0
24 19 11 13 0
25 21 24 14 0
26 24 25 20 0
27 14 20 9 0
28 24 20 14 0
29 25 23 20 0
30 23 22 15 0
31 20 15 9 0
32 23 15 20 0
];

%% Variable Prescribed
% Node            Dimension                Value

velocity = [
1 1 0
1 2 0
5 1 0
5 2 0
6 1 0
6 2 0
16 1 0
16 2 0
17 1 0
17 2 0
21 1 0
21 2 0
22 1 0
22 2 0
25 1 0
25 2 0
];

pressure = [
25 1 0
];

nu=1;
Vol_force = @(x,y){nu*4*(x^3*((6 - 12*y)) + x^4*((-3 + 6*y)) + y*(1 - 3*y + 2*y^2)-6*x*y*(1 - 3*y + 2*y^2)...
+ 3*x^2*((-1 + 4*y - 6*y^2 + 4*y^3)))+(2*x-1)*(y-1);
-4*nu*(-3*((-1 + y)^2)*y^2 - 3*x^2*((1 - 6*y + 6*y^2))+2*x^3*((1 - 6*y + 6*y^2)) +...
x*(1 - 6*y + 12*y^2 - 12*y^3 + 6*y^4))+x*(x-1)};

% nu=1;
% Vol_force = @(x,y){nu*4*(x^3*((6 - 12*y)) + x^4*((-3 + 6*y)) + y*(1 - 3*y + 2*y^2)-6*x*y*(1 - 3*y + 2*y^2)...
% + 3*x^2*((-1 + 4*y - 6*y^2 + 4*y^3)))+nu*y*(1 - y)*(1 - 2*x);
% -4*nu*(-3*((-1 + y)^2)*y^2 - 3*x^2*((1 - 6*y + 6*y^2))+2*x^3*((1 - 6*y + 6*y^2)) +...
% x*(1 - 6*y + 12*y^2 - 12*y^3 + 6*y^4)+x*(1 - x)*(1 - 2*y))};

%% Group Elements
% Element        Group_num

Group = [
];

%% Initial Holes
% Elements that are considered holes initially
% Element

Initial_holes = [
];

%% Boundary Elements
% Elements that can not be removed
% Element

Boundary_elements = [
];

%% Micro gauss post
%
% Element

Micro_gauss_post = [
];

%% Micro Slave-Master
% Nodes that are Slaves
% Nodes             Value (1-Slave,0-Master)

Micro_slave = [
];

%% Nodes solid
% Nodes that must remain
% Nodes

% nodesolid = unique(pointload_complete(:,1));

%% External border Elements
% Detect the elements that define the edge of the domain
% Element               Node(1)           Node(2)

External_border_elements = [
1 3 6
3 1 3
6 5 2
7 2 1
9 16 18
9 10 16
10 18 21
11 6 10
21 22 19
22 17 11
22 19 17
23 11 5
25 21 24
26 24 25
29 25 23
30 23 22
];

%% External border Nodes
% Detect the nodes that define the edge of the domain
% Node

External_border_nodes = [
];

%% Materials
% Materials that have been used
% Material_Num              Mat_density        Young_Modulus        Poisson

Materials = [
];

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function keypointsAll= det2keypointsImgDPM(expidx,parts,annolist_gt,nDet)

fprintf('det2keypointsImg()\n');

if (nargin < 4)
    nDet = 1;
end

p = rcnn_exp_params(expidx);
conf = rcnn_config('sub_dir', '/cachedir/test', 'exp_dir', [p.expDir '/' p.shortName]);

if (isfield(p,'predDir') && ~bTrain)
    predDir = p.predDir;
else
    predDir = [conf.cache_dir '/pred'];
end

if (nDet == 1)
    keypointsAll = repmat(struct('imgname','','det',nan(16,3)), length(annolist_gt), 1);
else
    for imgidx = 1:length(annolist_gt)
        keypointsAll(imgidx).imgname = '';
        keypointsAll(imgidx).det = cell(16,1);
    end
end

if (isfield(p,'nProposals'))
    nProposals = p.nProposals;
else
    nProposals = inf;
end

if (isfield(p,'min_det_score'))
    min_det_score = p.min_det_score;
else
    min_det_score = -inf;
end

for imgidx = 1:length(annolist_gt)

    fprintf('.');

    fname = [predDir '/imgidx_' padZeros(num2str(imgidx-1),5)];

    if (bNMS)
        load(fname,'aboxes');
        boxes = aboxes;
    else
        load(fname,'aboxes_nonms');
        boxes = aboxes_nonms;

        if (nms_thresh_det > 0)
            assert(expidx == 454 || expidx == 464 || expidx == 466 || expidx == 471);
            for i = 1:length(boxes)
                boxes{i}(:,1:4) = rcnn_scale_bbox(boxes{i}(:,1:4),1.0/p.scale,1000,1000);
            end
            idxsNMSall = [];
            for i = 1:length(boxes)
                I = nms_IoMin([boxes{i}(:,1:4) boxes{i}(:,5)], nms_thresh_det);
                idxsNMSall = sort(unique([idxsNMSall; I]));
            end
            for i = 1:length(boxes)
                boxes{i} = boxes{i}(idxsNMSall,:);
            end
        end
        if (isfield(p,'max_detections') && p.max_detections < inf)
            assert(expidx == 464 || expidx == 466 || expidx == 471);
            unProb = zeros(size(boxes{1},1),length(boxes));
            for i = 1:length(boxes)
                unProb(:,i) = boxes{i}(:,5);
            end
            scores = sort(reshape(unProb, size(unProb,1)*size(unProb,2),1));
            unProbThresh = unProb;
            j = 1;
            while (length(find(sum(unProbThresh,2) > 0)) > p.max_detections && j <= length(scores))
                unProbThresh = unProb;
                for cidx = 1:length(boxes)
                    unProbThresh(unProb(:,cidx) < scores(j),cidx) = 0;
                end
                j = j + 1;
                %                     length(find(sum(unProbThresh,2) > 0))
            end
            min_det_score_found = scores(j);
            idxs = find(sum(unProbThresh,2) > 0);
            for i = 1:length(boxes)
                boxes{i} = boxes{i}(idxs,:);
            end
        end
    end

    keypointsAll(imgidx).imgname = annolist_gt(imgidx).image.name;
    assert(length(p.pidxs) == length(boxes));

    for i = 1:length(p.pidxs)

        pidx = p.pidxs(i);
        % part is a joint
        assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
        jidx = parts(pidx+1).pos(1);

        det = boxes{i};

        if (nDet == 1)

            maxtop = min(size(det,1),nProposals);
            [val,idx] = max(det(1:maxtop,5));
            x = mean(det(idx,[1 3]));
            y = mean(det(idx,[2 4]));

            if (val < min_det_score)
                x = inf;
                y = inf;
            end
%             maxtop = min(size(det,1),100);
%             [~,idxs] = sort(det(:,5),'descend');
%             det = det(idxs(1:maxtop),:);
%             x = mean(det(:,[1 3]),2);
%             y = mean(det(:,[2 4]),2);
%             keep = nms_dist([x,y,det(1:maxtop,end)],0.1,50);
%             det = det(keep,:);
%             [val,idx] = max(det(:,5));
%             x = mean(det(idx,[1 3]));
%             y = mean(det(idx,[2 4]));

            keypointsAll(imgidx).det(jidx+1,:) = [[x y] val];
        else
            score = det(:,5);
            [val,idxs] = sort(score,'descend');
            idxs = idxs(1:min(nDet,length(idxs)));
            % use score threshold
            idxs2 = find(score(idxs) >= min_det_score);
            idxs = idxs(idxs2);
            if (~isempty(idxs))
                x = mean(det(idxs,[1 3]),2);
                y = mean(det(idxs,[2 4]),2);
                keypointsAll(imgidx).det{jidx+1,:} = [[x y] det(idxs,5)];
            else
                keypointsAll(imgidx).det{jidx+1,:} = [[inf inf] min_det_score];
            end
        end

    end
    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(keypointsAll));
    end
end
fprintf(' done\n');

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut