text:
function [data, xAxis, yAxis, misc] = impload(filename)
% Reads intensity map data from PerkinElmer block structured files.
% This version is compatible with '1994' standard IMP files.
% IMP files have 2 constant data intervals.
%
% [data, xAxis, yAxis, misc] = impload(filename):
%   data: length(x) * length(y) 2D array
%   xAxis: vector for horizontal axis (e.g. micrometers)
%   yAxis: vector for vertical axis (e.g. micrometers)
%   misc: miscellanous information in name,value pairs

% Copyright (C)2007 PerkinElmer Life and Analytical Sciences
% Stephen Westlake, Seer Green
%
% History
% 2007-04-24 SW     Initial version

% Block IDs
DS3C2IData            = 5108;   % 3D DS: header info
DS3C2IPts             = 5109;   % 3D DS: point coordinates (CvCoOrdArray)
DS3C2IPtsCont         = 5110;   % 3D DS: point coordinates (continuator for large CvCoOrdArray)
DS3C2IHistory         = 5111;   % 3D DS: history record
DSet3DC2DIBlock       = 5112;   % 3D DS: the block containing all other 3d blocks.

fid = fopen(filename,'r');
if fid == -1
    error('Cannot open the file.');
    return
end

% Fixed file header of signature and description
signature = setstr(fread(fid, 4, 'uchar')');
if ~strcmp(signature, 'PEPE')
    error('This is not a PerkinElmer block structured file.');
    return
end
description = setstr(fread(fid, 40, 'uchar')');

% Initialize a variable so we can tell if we have read it.
xLen = int16(0);

% The rest of the file is a list of blocks
while ~feof(fid)
    blockID = fread(fid,1,'int16');
    blockSize = fread(fid,1,'int32');

    % feof does not go true until after the read has failed.
    if feof(fid)
        break
    end

    switch blockID
        case DS3C2IData % Dataset header
            len = fread(fid,1,'int16');
            alias = setstr(fread(fid, len, 'uchar')');

            xDelta = fread(fid, 1, 'double');
            yDelta = fread(fid, 1, 'double');
            x0 = fread(fid, 1, 'double');
            y0 = fread(fid, 1, 'double');
            xLen = fread(fid, 1, 'int32');
            yLen = fread(fid, 1, 'int32');

            len = fread(fid, 1, 'int16');
            xLabel = setstr(fread(fid, len, 'uchar')');
            len = fread(fid, 1, 'int16');
            yLabel = setstr(fread(fid, len, 'uchar')');
            len = fread(fid, 1, 'int16');
            zLabel = setstr(fread(fid, len, 'uchar')');

        case DS3C2IPts  % This contains the first block
            count = blockSize / 8;
            data = fread(fid, count, 'double');

        case DS3C2IPtsCont % Continuation blocks
            countCont = blockSize / 8;
            data(count + 1 : count + countCont) = fread(fid, countCont, 'double');
            count = count + countCont;

        otherwise               % unknown block, just seek past it
            fseek(fid, blockSize, 'cof');
    end
end
fclose(fid);

if xLen == 0
    error('The file does not contain spectral image data.');
    return
end

% Reshape the linear array we have read in and flip so the
% data is in the correct aspect for image(data) with 1,1 at the top.
data = flipud(transpose(reshape(data, xLen, yLen)));

% Expand the axes specifications into vectors
% The yAxis is reversed to match the flip.
xEnd = x0 + (xLen - 1) * xDelta;
yEnd = y0 + (yLen - 1) * yDelta;
xAxis = x0 : xDelta : xEnd;
yAxis = yEnd : -yDelta : y0;

% Return the other details as name,value pairs
misc(1,:) = {'xLabel', xLabel};
misc(2,:) = {'yLabel', yLabel};
misc(3,:) = {'zLabel', zLabel};
misc(4,:) = {'alias', alias};

author: Sable


text:
function getSGLogs(tg, folder)

% Author: Alexander Wischnewski     Date: 21-12-2018
%
% Description:
%   copies the log data from the speedgoat to the host computer and
%   converts them to the .mat format
%
% Input:
%   tg:     Speedgoat target object
%   folder: If given, the files are copied to the speciefied folder.
%           Otherwise the current folder is used.

% standard folder for logs
std_folder = 'C:/RoboSim/Logs/';
% check if standard folder exits
if(~exist(std_folder, 'dir'))
    mkdir(std_folder);
end

% configure time formats for automatic folder creation
format_day = 'yy_mm_dd';
format_time = 'HH_MM_SS';

% if no target folder was specified, use standard folder
if(nargin<2)
    % check if a folder with the current date exists
    if(~exist([std_folder datestr(now, format_day)], 'dir'))
       % create a new folder if this does not exits
       mkdir([std_folder datestr(now, format_day)]);
    end
    % create a new folder with the current time stamp
    folder = [std_folder datestr(now, format_day) '/'...
        datestr(now, format_time)];
    mkdir(folder);
end

% get ftp object on target
ftp = SimulinkRealTime.openFTP(tg);

% switch to logs folder
cd(ftp,'C:\Logs\');

% get directory of log files
d = dir(ftp);

% retrieve all log files into current folder
disp('Copy files from RMC');
for i = 3:numel(d)
    fileName = [d(i).name];
    disp(['Copy ' fileName]);
    mget(ftp, fileName, folder);
end

disp('Convert data to .mat');
old_folder = pwd;
cd(folder);
convertSGLogs(pwd);
cd(old_folder);

disp('Finished');

author: TUMFTM


text:
% Demonstrate how negative values arise from sampling variability when
% there is no effect

% trials per class
Ntrl_class = 100;
% sampling repetitions
Nrep = 1000;

Inobc = zeros(1,Nrep);
Ibc = zeros(1,Nrep);
for ri=1:Nrep
    % generate data
    % the class has no effect on the distribution the data is drawn from
    x = randn(2*Ntrl_class,1);
    cx = copnorm(x);
    y = [zeros(Ntrl_class,1); ones(Ntrl_class,1)];

    % without bias correction
    Inobc(ri) = mi_model_gd(cx,y,2,false,true);
    % with bias correction
    Ibc(ri) = mi_model_gd(cx,y,2,true,true);
end

%%
figure
subplot(1,2,1)
hist(Inobc,20);
vline(mean(Inobc),'r')
title('mi\_model\_gd, no bias correction')
set(gca,'FontSize',12)

subplot(1,2,2)
hist(Ibc,20);
vline(mean(Ibc),'r')
title('mi\_model\_gd, with bias correction')
set(gca,'FontSize',12)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [id,dist] = find(S3G,ori,varargin)
% return indece and distance of all nodes within a eps neighborhood
%
% Syntax
%   % find the closes point
%   [ind,dist] = find(SO3G,ori)
%
%   % find points with a radius
%   [ind,dist] = find(SO3G,ori,radius)
%
%   % find cube corners
%   cubeInd = find(SO3G,ori,'cube')
%
% Input
%  SO3G   - @homochoricSO3Grid
%  ori    - @orientation
%  radius - double
%
% Output
%  ind  - index of the closes grid point
%  cubeInd - Nx8 list of indeces of the cube corners containing ori
%  dist - misorientation angle
%

% project oris to fundamental Region
ori = project2FundamentalRegion(ori, S3G.CS, S3G.SS);

% translate input (ori) into cubochoric coordinates
xyz = quat2cube(ori);

% N intervals of length hres along each edge of the cube
N = round(2 * pi / S3G.res);
hres = pi^(2/3) / N;

if nargin == 2 % closest point

    % calculate grid index along each axis of the cube
    % let the grid have N points along each axis
    % then xyz/hres takes values in [-N/2,N/2]
    % 0 is included in the grid iff N is odd
    % so xyz/hres+N/2-1/2 takes values in -1/2,...,N-1/2
    % after rounding this should yield values in 0,...,N-1
    sub  = mod(round(xyz/hres + N/2 - 0.5),N) + 1;    % [ix, iy, iz] each from 1 to N
    id = sub2ind(S3G,sub(:,1),sub(:,2),sub(:,3));

    if nargout == 2
        dist = zeros(size(id));     % initialize distances
        dist(id>0) = angle(ori.subSet(id>0), S3G.subSet(id(id>0)), 'noSymmetry');
    end

    % project those not beeing in fR again --> Quaternion Projected
    if any(id==0)

        % take outside grid points
        subxyz = hres * (sub(id==0,:) - (N+1)/2);

        % and project them back into the fundamental region
        q2 = project2FundamentalRegion(cube2quat(subxyz),S3G.CS,S3G.SS);

        % find closest grid point
        xyz2 = quat2cube(q2);
        sub2 = mod(round(xyz2/hres + N/2 - 0.5),N) + 1;

        % insert eventually new found neighbors on grid into id
        id2 = sub2ind(S3G,sub2(:,1),sub2(:,2),sub2(:,3));

        % get indice of zeros in id
        iszero = find(~id);

        % compute distances where id is zero and id2 isnt zero
        % (new found valid neighbors)
        if nargout == 2
            dist(iszero(id2>0)) = angle(ori.subSet(iszero(id2>0)), S3G.subSet(id2(id2>0)));
        end

        % plug in new found indice
        id(~id) = id2;

    end

    % if still not in fR, then search for best neighbor (of 8 surrounding)
    % notice that at least one of those 8 neighbors will lie in S3G
    if any(id==0)

        % get indice of surrounding vertices of the whole grid for those points
        % who dont have a valid nearest neighbor yet
        % for more info see docu of getS3Gvertices
        idx = getS3GVertices(S3G, xyz(id==0,:));

        % initialize set of nearest neighbors all as invalid
        S3Gvertices = rotation.nan(sum(id==0), 8);

        % and only overwrite this ones where idx > 0 (valid neighbor found)
        S3Gvertices(idx>0) = S3G.subSet(idx(idx>0));

        % calculate distances ...
        d = angle(repmat(ori(id==0), 1, 8), S3Gvertices, 'noSymmetry');

        % ... to find the nearest point
        [dmin, idmin] = min(d, [], 2);

        % overwrite the distance, if 2 outputs
        if nargout == 2
            dist(id==0) = dmin;
        end

        % overwrite new found indice
        % all indice now have a valid value ~= 0 (>= 1 of 8 neighbors in fR)
        id(id==0) = idx(sub2ind(size(idx), [1:sum(id==0)]', idmin));

    end

elseif ischar(varargin{1}) % cube c

    % get indice of surrounding vertices of the whole grid for those points
    % who dont have a valid nearest neighbor yet
    % for more info see docu of getS3Gvertices
    id = getS3GVertices(S3G, xyz);

    % calculate coordinate-wise distances to the floor corner
    % (attained by rounding coordinates down to the grid)
    % care that it mathers if the grid is odd or even since we then have to
    % round towards multiples of hres or the same shifted by hres/2
    dist = mod((N+1)/2*hres+xyz,hres);

else % neighborhood

end

end

author: mtex-toolbox

repo: mtex

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


text:
function test_bug2471

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_timelockgrandaverage

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2471.mat'));

cfg = [];
cfg.channel = 'all';
cfg.keepindividual = 'no';
cfg.method = 'within';
out_data = ft_timelockgrandaverage(cfg,block{:});


text:
function E = stToEdges( S, nms, suppress )
% Convert sketch tokens to edges.
%
% USAGE
%  E = stToEdges( S, [nms], [suppress] )
%
% INPUTS
%  S          - [h x w x (nTokens+1)] sketch token probability maps
%  nms        - [1] if true apply non-maximum suppression to edges
%  suppress   - [1] if true suppress boundary confidence in order to
%                   compensate padding effect
%
% OUTPUTS
%  E          - [h x w] edge probability map
%
% EXAMPLE
%
% See also stDetect
%
% Sketch Token Toolbox     V0.95
% Copyright 2013 Joseph Lim [lim@csail.mit.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see bsd.txt]

    % get default arguments
    if(nargin<2 || isempty(nms)),
        nms=1;
    end
    if(nargin<3 || isempty(suppress)),
        suppress=1;
    end

    % extract edge probabilities
    E = 1-S(:,:,end);

    % apply nms
    if nms
        O=edgeOrient(E,4);
        E=edgeNms(E,O,1);
    end

    % apply boundary suppression
    if suppress
        S=5;
        for s=1:S,
            E([s end-s+1],:,:)=E([s end-s+1],:,:)*(s-1)/S;
        end
        for s=1:S,
            E(:,[s end-s+1],:)=E(:,[s end-s+1],:)*(s-1)/S;
        end
    end

end

function E = edgeNms( E, O, r )
    % suppress locations where edge is stronger in orthogonal direction
    E1=imPad(E,r+1,'replicate');
    Dx=cos(O);
    Dy=sin(O);

    [ht,wd]=size(E1);
    [cs,rs]=meshgrid(r+2:wd-r-1,r+2:ht-r-1);
    for i=-r:r,
        if i==0
            continue;
        end
        cs0=i*Dx+cs;
        dcs=cs0-floor(cs0);
        cs0=floor(cs0);

        rs0=i*Dy+rs;
        drs=rs0-floor(rs0);
        rs0=floor(rs0);

        E2 = (1-dcs).*(1-drs) .* E1(rs0+0+(cs0-1)*ht);
        E2 = E2 + dcs.*(1-drs) .* E1(rs0+0+(cs0-0)*ht);
        E2 = E2 + (1-dcs).*drs .* E1(rs0+1+(cs0-1)*ht);
        E2 = E2 + dcs.*drs .* E1(rs0+1+(cs0-0)*ht);
        E(E*1.01<E2) = 0;
    end
end

function O = edgeOrient( E, r )
    % compute very approximate orientation map from edge map
    E2=convTri(E,r);
    f=[-1 2 -1];

    Dx=conv2(E2,f,'same');
    Dy=conv2(E2,f','same');
    F=conv2(E2,[1 0 -1; 0 0 0; -1 0 1],'same')>0;
    Dy(F)=-Dy(F);
    O=mod(atan2(Dy,Dx),pi);
end

author: zouchuhang

repo: LayoutNet


text:
function r8col_to_r8vec_test ( )

%*****************************************************************************80
%
%% R8COL_TO_R8VEC_TEST tests R8COL_TO_R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8COL_TO_R8VEC_TEST\n' );
  fprintf ( 1, '  R8COL_TO_R8VEC converts an array of columns to a vector.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : m
    for j = 1 : n
      a(i,j) = 10 * i + j;
    end
  end

  r8mat_print ( m, n, a, '  The array of columns:' );

  x = r8col_to_r8vec ( m, n, a );

  r8vec_print ( m*n, x, '  The resulting vector of columns:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [c, s] = foopsi(y, g, lam, solver)
%% Infer the most likely discretized spike train underlying an AR(1) fluorescence trace
% Solves the sparse non-negative deconvolution problem
%  min 1/2|c-y|^2 + lam |s|_1 subject to s_t = c_t-g c_{t-1} >= 0

%% inputs:
%   y:  T*1 vector, One dimensional array containing the fluorescence intensities
%withone entry per time-bin.
%   g:  scalar, Parameter of the AR(1) process that models the fluorescence ...
%impulse response.
%   lam:  scalar, sparsity penalty parameter lambda.
%   solver: string, optimization solver

%% outputs
%   c: T*1 vector, the inferred denoised fluorescence signal at each time-bin.
%   s: T*1 vector, discetized deconvolved neural activity (spikes)

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging

%% initialization
y = reshape(y, [], 1);
T = length(y);
if ~exist('g', 'var') || isempty(g);   g = estimate_time_constant(y, 2); end

if ~exist('lam', 'var') || isempty(lam);   lam = 0; end
if ~exist('solver', 'var') || isempty(smin);   solver = 'SDPT3'; end

% construct the deconvolution matrix (s=G*c)
len_g = length(g);
g = reshape(g, 1, []);
indc = bsxfun(@minus, (1:T)', 0:len_g);
indr = repmat((1:T)', [1, len_g+1]);
v = ones(T,1)*[1,-g];
ind = (indc>0);
G = sparse(indr(ind), indc(ind), v(ind), T, T);

%% run optimization
% cvx_solver(solver);

cvx_begin quiet
variable c(T)
minimize(0.5*(c-y)'*(c-y) + lam*(1-sum(g))*norm(c,1))
subject to
G*c >=0;
cvx_end

s = reshape(G*c, 1, []);
s(1) = 0;
c = reshape(c,1, []);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [out1,out2,out3] = ginput(arg1)
    %GINPUT Graphical input from mouse.
    %   [X,Y] = GINPUT(N) gets N points from the current axes and returns
    %   the X- and Y-coordinates in length N vectors X and Y.  The cursor
    %   can be positioned using a mouse (or by using the Arrow Keys on some
    %   systems).  Data points are entered by pressing a mouse button
    %   or any key on the keyboard except carriage return, which terminates
    %   the input before N points are entered.
    %
    %   [X,Y] = GINPUT gathers an unlimited number of points until the
    %   return key is pressed.
    %
    %   [X,Y,BUTTON] = GINPUT(N) returns a third result, BUTTON, that
    %   contains a vector of integers specifying which mouse button was
    %   used (1,2,3 from left) or ASCII numbers if a key on the keyboard
    %   was used.

    %   Copyright 1984-2001 The MathWorks, Inc.
    %   $Revision: 1399 $  $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    out1 = []; out2 = []; out3 = []; y = [];
    c = computer;
    if ~strcmp(c(1:2),'PC')  &&  ~strcmp(c(1:2),'MA')
        tp = get(groot,'TerminalProtocol');
    else
        tp = 'micro';
    end

    if ~strcmp(tp,'none') & ~strcmp(tp,'x')  &&  ~strcmp(tp,'micro')
        if nargout == 1
            if nargin == 1
                out1 = trmginput(arg1);
            else
                out1 = trmginput;
            end
        elseif nargout == 2  ||  nargout == 0
            if nargin == 1
                [out1,out2] = trmginput(arg1);
            else
                [out1,out2] = trmginput;
            end
            if  nargout == 0
                out1 = [out1 out2 ];
            end
        elseif nargout == 3
            if nargin == 1
                [out1,out2,out3] = trmginput(arg1);
            else
                [out1,out2,out3] = trmginput;
            end
        end
    else

        fig = gcf;
        figure_w_normalized_uicontrolunits(gcf);

        if nargin == 0
            how_many = -1;
            b = [];
        else
            how_many = arg1;
            b = [];
            if  ischar(how_many) ...
                    | size(how_many,1) ~= 1 | size(how_many,2) ~= 1 ...
                    | ~(fix(how_many) == how_many) ...
                    | how_many < 0
                error('Requires a positive integer.')
            end
            if how_many == 0
                ptr_fig = 0;
                while(ptr_fig ~= fig)
                    ptr_fig = get(groot,'PointerWindow');
                end
                scrn_pt = get(groot,'PointerLocation');
                loc = get(fig,'Position');
                pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
                out1 = pt(1); y = pt(2);
            elseif how_many < 0
                error('Argument must be a positive integer.')
            end
        end

        % Remove figure button functions
        state = uisuspend(fig);
        pointer = get(gcf,'pointer');
        set(gcf,'pointer','crosshair');
        fig_units = get(fig,'units');
        char = 0;

        while how_many ~= 0
            % Use no-side effect WAITFORBUTTONPRESS
            waserr = 0;
            try
                keydown = wfbp;
            catch
                waserr = 1;
            end
            if(waserr == 1)
                if(ishandle(fig))
                    set(fig,'units',fig_units);
                    uirestore(state);
                    error('Interrupted');
                else
                    error('Interrupted by figure deletion');
                end
            end

            ptr_fig = get(groot,'CurrentFigure');
            if(ptr_fig == fig)
                if keydown
                    char = get(fig, 'CurrentCharacter');
                    button = abs(get(fig, 'CurrentCharacter'));
                    scrn_pt = get(groot, 'PointerLocation');
                    set(fig,'units','pixels')
                    loc = get(fig, 'Position');
                    pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
                    set(fig,'CurrentPoint',pt);
                else
                    button = get(fig, 'SelectionType');
                    if strcmp(button,'open')
                        button = b(length(b));
                    elseif strcmp(button,'normal')
                        button = 1;
                    elseif strcmp(button,'extend')
                        button = 2;
                    elseif strcmp(button,'alt')
                        button = 3;
                    else
                        error('Invalid mouse selection.')
                    end
                end
                pt = get(gca, 'CurrentPoint');

                how_many = how_many - 1;

                if(char == 13) % & how_many ~= 0)
                    % if the return key was pressed, char will == 13
                    % and that's our signal to break out of here whether
                    % or not we have collected all the requested data
                    % points.
                    % If this was an early breakout, don't include
                    % the <Return> key info in the return arrays.
                    % We will no longer count it if it's the last input.
                    break;
                end

                out1 = [out1;pt(1,1)];
                y = [y;pt(1,2)];
                b = [b;button];
            end
        end

        uirestore(state);
        set(fig,'units',fig_units);

        if nargout > 1
            out2 = y;
            if nargout > 2
                out3 = b;
            end
        else
            out1 = [out1 y];
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key = wfbp
    %WFBP   Replacement for WAITFORBUTTONPRESS that has no side effects.

    fig = gcf;
    current_char = [];

    % Now wait for that buttonpress, and check for error conditions
    waserr = 0;
    try
        h=findall(fig,'type','uimenu','accel','C');   % Disabling ^C for edit menu so the only ^C is for
        set(h,'accel','');   % interrupting the function.
        pause(0.1)
        keydown = waitforbuttonpress;
        current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
        if~isempty(current_char) & (keydown == 1)           % If the character was generated by the
            if(current_char == 3)                       % current keypress AND is ^C, set 'waserr'to 1
                waserr = 1;                             % so that it errors out.
            end
        end

        set(h,'accel','C');                                 % Set back the accelerator for edit menu.
    catch
        waserr = 1;
    end
    drawnow;
    if(waserr == 1)
        set(h,'accel','C');                                % Set back the accelerator if it errored out.
        error('Interrupted');
    end

    if nargout>0, key = keydown; end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: CelsoReyes


text:
function noise = noisevector(vlength,xcfunction,noisevariance)
% function noise = noisevector(vlength,xcfunction,noisevariance)
%
% returns a noise vector of 'vlength' samples
% with characteristics of the autocorrelation function that you specify
%   (must be same sampling rate as timeseries!)
%
% the vector is normalized such that the variance is what you specify
% and the mean is zero
%
% vlength: vector length
% xcfunction: any cross-correlation (i.e., noise autocorrelation) function
%   - also called a periodogram
%   - a 1/f function is generally good for fMRI data
% noisevariance: desired variance of noise
%
% 2/11/01 Tor Wager

if isempty(xcfunction), xcfunction = [1 0];,end

% define series of random "shocks".  Last one is time 1, paramount one is time 2, etc.
%       pad with zeros, because 1st shock has no history to influence it.
xclength = size(xcfunction,2);
a = randn(1,vlength); a(end + 1:end + xclength) = 0;
% the noise is the series of shocks multiplied by the autocorrelation function
% so the current value of the system = weighted sum of past shocks, up to size of autocorr function
for i = 1:vlength
    noise(i) = dot(a(end-i-xclength+1:end-i),xcfunction');
end

% make sure that the variance is one and the mean is zero

noise = noise - mean(noise);

noise = noise / sqrt(var(noise));

% now make the variance match your estimate of the variance

noise = noise * sqrt(noisevariance);

%disp(['Noise variance is ' num2str(var(noise))])
%figure;plot(noise)

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function trans = createTranslation(varargin)
%CREATETRANSLATION Create the 3*3 matrix of a translation.
%
%   TRANS = createTranslation(DX, DY);
%   Returns the translation corresponding to DX and DY.
%   The returned matrix has the form :
%   [1 0 TX]
%   [0 1 TY]
%   [0 0  1]
%
%   TRANS = createTranslation(VECTOR);
%   Returns the matrix corresponding to a translation by the vector [x y].
%
%
%   See also
%   transforms2d, transformPoint, createRotation, createScaling

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-04-06
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% process input arguments
if isempty(varargin)
    tx = 0;
    ty = 0;
elseif length(varargin)==1
    var = varargin{1};
    tx = var(1);
    ty = var(2);
else
    tx = varargin{1};
    ty = varargin{2};
end

% create the matrix representing the translation
trans = [1 0 tx ; 0 1 ty ; 0 0 1];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function[] = hindawi()
%
% Simulation code for [1],
%
% [1] Paul Rodriguez "Review: Total Variation Regularization Algorithms for
%  		      Images Corrupted With Different Noise Models"
%
%
%
% Legal:
%
% Authors
%   Paul Rodriguez    prodrig@pucp.edu.pe
%

l2TV         = 1;
l2TV_nqp     = 2;
l1TV         = 3;
l1TV_adapt   = 4;
poiTV_denoi  = 5;
poiTV_deconv = 6;
gammaTV      = 7;
mixTV        = 8;
nTV          = mixTV;

% Test Images (for each case)

l2Imgs = {'goldhl', 'Clena'};           % To be corrupted with Gaussian noise
l1Imgs = {'Cbarb', 'boats'};            % To be corrupted with Salt & Pepper noise
poiImgs = {'cman512', 'cpeppers'};     % To be corrupted with Poisson noise
gammaImgs = {'tank', 'Clena'};          % To be corrupted with Gamma noise
mixImgs = {'cman512'};                  % To be corrupted with mixed Gaussian + Impulse noise

Tests = cell(nTV-1, 1);

Tests{l2TV} = l2Imgs;
Tests{l2TV_nqp} = l2Imgs;
Tests{l1TV} = l1Imgs;
Tests{l1TV_adapt} = l1Imgs;
Tests{poiTV_denoi} = poiImgs;
Tests{poiTV_deconv} = poiImgs;

%  Tests{gammaTV} = gammaImgs;
%  Tests{mixTV} = mixImgs;

% kernel

  kernel = fspecial('disk',3.2);

  K = @(x) imfilter(x, kernel, 'symmetric','conv');

  KT = @(x) K(x);
  KC = {K, KT};

% Noise level (define the noise level to corrupt images)

l2Noise = [0.1 0.2];
l1Noise = [0.3 0.5 0.8];
randomNoise = [0.1 0.2 0.3];
poiNoise = [5 30 100 255];

Noise = cell(nTV-1, 1);

Noise{l2TV} = l2Noise;
Noise{l2TV_nqp} = l2Noise;
Noise{l1TV} = l1Noise;
Noise{l1TV_adapt} = l1Noise;
Noise{poiTV_denoi} = poiNoise;
Noise{poiTV_deconv} = poiNoise;

% TV Noise Model (to corrupt imput images)

NoiseModel = cell(nTV-1, 1);

NoiseModel{l2TV} = @(Img, sigma) imnoise(K(Img),'gaussian', 0, (sigma^2) );
NoiseModel{l2TV_nqp} = @(Img, sigma) imnoise(K(Img),'gaussian', 0, (sigma^2) );
NoiseModel{l1TV} = @(Img, spnoise) imnoise(Img, 'salt & pepper', spnoise);
NoiseModel{l1TV_adapt} = @(Img, spnoise) imnoise(Img, 'salt & pepper', spnoise);

NoiseModel{poiTV_denoi} = @(Img, M) poissrnd( M*NormalizeMax(Img) );
NoiseModel{poiTV_deconv} = @(Img, M) poissrnd( M*NormalizeMax(K(Img)) );

% ----------------------
% Define TV parameters
% ----------------------

% Loops

TVLoops = zeros(nTV-1, 1);

TVLoops(l2TV)           = 5;
TVLoops(l2TV_nqp)       = 5;
TVLoops(l1TV)           = 4;
TVLoops(l1TV_adapt)     = 4;

TVLoops(poiTV_denoi)    = 6;
TVLoops(poiTV_deconv)   = 6;

% regularization parameters

TVLambda = zeros(nTV-1, 5);

TVLambda(l2TV,1) = 0.05; TVLambda(l2TV,2) = 0.1;
TVLambda(l2TV_nqp,1) = 0.05; TVLambda(l2TV_nqp,2) = 0.1;

TVLambda(l1TV,1) = 1.2; TVLambda(l1TV,2) = 1.4; TVLambda(l1TV,3) = 1.5;
TVLambda(l1TV_adapt,1) = 1.0; TVLambda(l1TV_adapt,2) = 1.0; TVLambda(l1TV_adapt,3) = 1.0;

TVLambda(poiTV_denoi,1:4) = 2./[0.35, 0.15, 0.075, 0.025 ];

%  Wrappers to TV functions)

RestoreTV = cell(nTV-1, 1);

RestoreTV{l2TV} = @(b, lambda, loops, KC) runl2TV(b, lambda, loops, KC);
RestoreTV{l2TV_nqp} = @(b, lambda, loops, KC) runl2TV_nqp(b, lambda, loops, KC);
RestoreTV{l1TV} = @(b, lambda, loops, dummy) runl1TV(b, lambda, loops, []);
RestoreTV{l1TV_adapt} = @(b, lambda, loops, dummy) runl1TV_Adapt(b, lambda, loops);

RestoreTV{poiTV_denoi} = @(b, lambda, loops, dummy) run_poiTV(b, lambda, loops, {});
RestoreTV{poiTV_deconv} = @(b, lambda, loops, KC) run_poiTV(b, lambda, loops, KC);

% additional vars

t1 = zeros(nTV-1, 2, 2);
t2 = zeros(nTV-1, 2, 2);

nmpdef;

%%%%%%%%%%%%%%%%%%%%%%%%
%        Setup         %
%%%%%%%%%%%%%%%%%%%%%%%%

SSIM_CODE = exist('ssim_index');
if( SSIM_CODE == 0 )
  disp('NOTE:');
  disp('  The function ssim_index (code for [SSIM]) is not in your path...');
  disp('  You may download it from:');
  disp('  http://www.ece.uwaterloo.ca/~z70wang/research/ssim/');
  disp('  ');
  disp('  Disabling SSIM reports');
  disp('  ');
  disp('  [SSIM]  Z. Wang, A. Bovik, H. Sheikh and E. Simoncelli, ');
  disp('         "Image quality assessment: From error visibility to');
  disp('         structural similarity " ');
  disp('         IEEE Transactions on Image Processing, 2004, 13:4(600-612).');

end

t=5;

str_all = sprintf(' \n');

while( iscell(Tests{t}) )

  str = sprintf(' \n'); str_all = [ str_all str ];
  str = sprintf(' Img \t\t Noise \t\t SNR (db) \t\t Time (s) \t\t SSIM \n');
  str_all = [ str_all str ];

  NImgs = length( Tests{t} );

  for k = 1: NImgs,

  switch lower( Tests{t}{k} )

    case{'lena'}
      I = double( imread('gray_imgs/lena_gray_512.png') ) / 255;
      SSIM_COLOR = 1;

    case{'peppers'}
      I = double( imread('gray_imgs/peppers_gray.png') ) / 255;
      SSIM_COLOR = 1;

    case{'boats'}
      I = double( imread('gray_imgs/boats_gray.png') ) / 255;
      SSIM_COLOR = 1;

    case{'goldhl'}
      I = double( imread('gray_imgs/goldhill_gray.png') ) / 255;
      SSIM_COLOR = 1;

    case{'cman'}
      I = double( imread('gray_imgs/cameraman_256x256.tiff') ) / 255;
      SSIM_COLOR = 1;

    case{'cman512'}
      I = double( imread('gray_imgs/cameraman_512x512.tiff') ) / 255;
      SSIM_COLOR = 1;

    case{'clena'}
      I = double( imread('color_imgs/lena_color_512.png') ) / 255;
      SSIM_COLOR = 0;

    case{'cbarb'}
      I = double( imread('color_imgs/barbara_color.png') ) / 255;
      SSIM_COLOR = 0;

    case{'cpeppers'}
      I = double( imread('color_imgs/peppers_color.png') ) / 255;
      SSIM_COLOR = 0;

    case{'cmandrill'}
      I = double( imread('color_imgs/mandrill_color.png') ) / 255;
      SSIM_COLOR = 0;

  end % --- END (switch) ---

  for n = 1:length(Noise{t})

    % add noise
    b = NoiseModel{t}( I, Noise{t}(n) );
    figure; imagesc( Normalize(b) );  colormap gray; axis image; axis off;

    % Restore via TV
    [u t1(t,k,n) t2(t,k,n)] = RestoreTV{t}(b, TVLambda(t,n), TVLoops(t), KC );

    % Restore performance
    [snrRec(t,k,n), ssimRec(t,k,n)] = computePerf(I, u, SSIM_CODE*SSIM_COLOR, 0);

    if(n==1)
      str = sprintf('%s \t\t %1.2f \t\t %2.2f \t\t %2.2f (%2.2f) \t\t %1.3f \n', ...
                    lower( Tests{t}{k} ), Noise{t}(n), snrRec(t,k,n), t1(t,k,n), t2(t,k,n), ssimRec(t,k,n) );
    else
      str = sprintf('     \t\t %1.2f \t\t %2.2f \t\t %2.2f (%2.2f) \t\t %1.3f \n', ...
                     Noise{t}(n), snrRec(t,k,n), t1(t,k,n), t2(t,k,n), ssimRec(t,k,n) );
    end
    str_all = [ str_all str ];

  end % --- END (FOR(n)) ---

  end % --- END (FOR(k)) ---

  disp(str_all);

 t = t + 1;

 if(t==6) break; end

end

end

% ============================================

%%%%%%%%%%%%%%%%%%%%%%%%
%     Normalize        %
%%%%%%%%%%%%%%%%%%%%%%%%

function[u] = Normalize(x)

u = (x - min(x(:)))/(max(x(:)) - min(x(:)));

end

% ============================================

function[u] = NormalizeMax(x)

u = x/max(abs(x(:)));

end

% ============================================

function[u, timePerf, t2] = runl2TV(b, lambda, loops, KC)

    pars = irntvInputPars('l2tv');

    pars.pcgtol_ini = 1e-4;
    pars.loops      = loops;
%      pars.U0         = b;

    pars.adapt_epsR   = 1;
    pars.epsR_cutoff  = 0.01;
    pars.adapt_epsF   = 1;
    pars.epsF_cutoff  = 0.05;

    tic;
    u = irntv(b, KC, lambda, pars);
    timePerf = toc;

    t2 = 0;

end

% ============================================

function[u, timePerf, t2] = runl2TV_nqp(b, lambda, loops, KC)

    pars = irntvInputPars('l2tv_nqp');

    pars.pcgtol_ini = 1e-4;
    pars.loops      = loops;
%      pars.U0         = b;

    pars.adapt_epsR   = 1;
    pars.epsR_cutoff  = 0.01;
    pars.adapt_epsF   = 1;
    pars.epsF_cutoff  = 0.1;

    pars.loops_NQP    = 15;
    pars.gamma_NQP    = 0.5e-3;
    pars.alpha_NQP    = 5e-1;
    pars.vmax_NQP     = 1+1;            % maximum value of output image

    tic;
    u = irntv(b+1, KC, lambda, pars)-1;
    timePerf = toc;

    t2 = 0;

end

% ============================================

function[u, timePerf, t2] = runl1TV(b, lambda, loops, KC)

    if(nargin < 6)
      flagTitle = 0;
    end

    pars = irntvInputPars('l1tv');

    pars.pcgtol_ini = 1e-4;
    pars.epsF       = 1e-2;
    pars.epsR       = 1e-4;
    pars.loops      = loops;

    tic;
    u = irntv(b, KC, lambda, pars);
    timePerf = toc;

    t2 = 0;

end

% ============================================

function[u, timePerf, t2] = runl1TV_Adapt(b, lambda, loops)

tic;

S0 = adaptMedian(b);

t2 = toc;

InputDims = size(S0);

if(length(InputDims) == 2)

    lambda_adapt  = 1e-6*(S0 == 0) + lambda*5*(S0 ~= 0);

else

    for d=1:3,
      lambda_adapt(:,:,d)  = 1e-6*(S0(:,:,d) == 0) + lambda*5.0*(S0(:,:,d) ~= 0);
    end

end

% ---------- Setup for IRN_adapt  ------------------

pars_irn_adapt = irntvInputPars('l1tv_adapt');

pars_irn_adapt.adapt_epsR   = 1;
pars_irn_adapt.epsR_cutoff  = 0.01;
pars_irn_adapt.adapt_epsF   = 1;
pars_irn_adapt.epsF_cutoff  = 0.1;

pars_irn_adapt.pcgtol_ini = 1e-4;

pars_irn_adapt.loops      = 1;

pars_irn_adapt.U0      = b;   % necessary the for adapt case

  tic;
  u = irntv(b, [], lambda_adapt, pars_irn_adapt);

  for m = 2:loops

    pars_irn_adapt.U0 = u;
    lambda_adapt = lambda_adapt*1.2;

    u = irntv(b, [], lambda_adapt, pars_irn_adapt);

  end

  timePerf = toc;

end

% ============================================

% ============================================

function[u, timePerf, t2] = run_poiTV(b, lambda, loops, KC)

    pars = irntvInputPars('tv_poisson_adapt');

    pars.pcgtol_ini = 1e-4;
    pars.loops      = 1;
%      pars.U0         = b;

    pars.adapt_epsR   = 1;
    pars.epsR_cutoff  = 0.1;
    pars.adapt_epsF   = 1;
    pars.epsF_cutoff  = 0.15;

    pars.loops_NQP    = 35;
    pars.gamma_NQP    = 5e-3;
    pars.alpha_NQP    = 5e-1;
    pars.vmax_NQP     = 1+1;            % maximum value of output image

    T = lambda;

    tic;
    u = irntv(b+0.01, KC, lambda, pars)-0.01;

    for ln=2:loops,

      [sig1, cv1, mu1, siginv1, cvinv1, S1, lmu1] = estNoise(u, 2*3+1, 100, 1);
      [sig2, cv2, mu2, siginv2, cvinv2, S2, lmu2] = estNoise(lmu1, 2*3+1, 100, 1);

      pars.U0 = u;
      vmin = min( pars.U0(:) );
      if(vmin < 0) pars.U0 = pars.U0 - vmin + 0.01; end

      alpha = 1.0;
      [Nrows Ncols Ndims] = size(u);

      mask = zeros(size(S2));

      for d = 1:Ndims,
        maskTmp = zeros(Nrows*Ncols,1);

        vS1 = S2(:,:,d); vS1 = vS1(:);

        tau = unimodal(S1(:,:,d));
        p2 = find( vS1 < alpha*tau );
        maskTmp(p2) = (1 - Normalize( vS1(p2) ) )*(1-0.85) + 0.85;

        p3 = find(vS1 >= alpha*tau);
        maskTmp(p3) = 1.0;

        mask(:,:,d) = reshape(maskTmp, [Nrows Ncols]);

      end % _END_ FOR(d)

      T = T.*mask;

      u = irntv(b+0.1, KC, T, pars)-0.1;

    end

    % ====================================
    timePerf = toc;

    t2 = 0;

end

% ============================================

function[snrRec, ssimRec] = computePerf(Img, u, flagSSIM, flagTitle)

    snrRec = snr(Img, u);
    if(flagSSIM)
      ssimRec = ssim_index(255*Normalize(Img), 255*Normalize(u));
    else
      ssimRec = NaN;
    end

    figure; imagesc( Normalize(u) );  colormap gray; axis image; axis off;
    if(flagTitle)
      title(sprintf('Restored Image\n SNR: %4.1fdB. SSIM: %1.2f ', ...
               snrRec, ssimRec));
    end

end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
% MPC applied to LOTKA-VOLTERRA system using
% a SINDYc model for different prediction horizon lengths.

clear all, close all, clc
figpath = '../FIGURES/';
datapath = '../DATA/';
addpath('../utils');

%% Generate Data
% Parameters: SINDy
polyorder = 3;
usesine = 0;

% True Parameters of Lotka-Volterra model
a = .5;
b = .025;
d = .005;
g = .5;
n = 2;
x0=[60; 50];
dt = .01;

% Choose forcing function to excite system for model identification
% forcing = @(x,t) [(0.33*(sin(1*t)+sin(.1*t)))];
forcing = @(x,t) [(2*(sin(1*t)+sin(.1*t))).^2];

% Integrate excited system
tspan=[0:dt:100];
u = forcing(0,tspan);
N = length(tspan);
options = odeset('RelTol',1e-10,'AbsTol',1e-10*ones(1,n));
[t,x]=ode45(@(t,x) lotkacontrol(t,x,forcing(x,t),a,b,d,g),tspan,x0,options);
plot(t,x,'LineWidth',1.5)
xlabel('Time')
ylabel('Population size')
legend('Prey','Predator')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', [figpath,'EX_LOTKA_Dynamics.eps']);

%% SINDYc Model Identification
% Compute Derivative
xold = x;
clear dx
x = xold;
eps = 0.1;
for i=1:length(x)
    dx(i,:) = lotkacontrol(0,x(i,:),u(i),a,b,d,g);
end
x = [xold u'];
dx(:,3) = 0*dx(:,2);
dx = dx + eps*randn(size(dx));
n = 3;

% Sparse regression
clear Theta Xi
Theta = poolData(x,n,polyorder,usesine);
m = size(Theta,2);

lambda = 0.001;      % lambda is our sparsification knob.
Xi = sparsifyDynamics(Theta,dx,lambda,n);
poolDataLIST({'x','y','u'},Xi,n,polyorder,usesine);

%% FIGURE 1:  Lotka-Volterra // Validation
x0 = x(end,1:2);
tspan = [100 200];
[tA,xA]=ode45(@(t,x)lotkacontrol(t,x,forcing(x,t),a,b,d,g),tspan,x0,options);   % true model
[tB,xB]=ode45(@(t,x)sparseGalerkinControl(t,x,forcing(x,t),Xi(:,1:2),polyorder,usesine),tspan,x0,options);  % approximate

h = figure;
subplot(2,1,1), box on
plot(t,x(:,1),'Color',[.4 .4 .4],'LineWidth',1.5), hold on
plot(tA,xA(:,1),'k','LineWidth',1.5), hold on
plot(tB,xB(:,1),'r--','LineWidth',1.5)
grid on
ylim([0 110])
% xlabel('Time','FontSize',13)
ylabel('Prey, x_1','FontSize',13)
set(gca,'FontSize',13, 'LineWidth',1)
subplot(2,1,2), box on
plot(t,x(:,2),'Color',[.4 .4 .4],'LineWidth',1.5), hold on
plot(tA,xA(:,2),'k','LineWidth',1.5), hold on
plot(tB,xB(:,2),'r--','LineWidth',1.5)
l1=legend('Training','Validation','SINDYc');
set(l1,'Location','NorthWest')
grid on
ylim([0 60])
ylabel('Predator, x_2','FontSize',13)
set(gca,'FontSize',13, 'LineWidth',1)
xlabel('Time','FontSize',13)
set(gca,'FontSize',13)

set(h,'Units','Inches');
set(gcf,'Position',[1 1 6. 5.5])
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperSize',[pos(3), pos(4)])
print(h,'-dpdf', [figpath,'EX_LOTKA_ControlValidation.pdf'],'-r0');
print(h,'-depsc2', [figpath,'EX_LOTKA_ControlValidation.eps'],'-r0');

clear ph
figure;hold on, box on,
ccolors = get(gca,'colororder');
plot([100,100],[0 260],':','Color',[0.4,0.4,0.4],'LineWidth',1)
text(5,120,'Training', 'FontSize',12)
text(105,120,'Prediction', 'FontSize',12)
plot(t,x(:,1),'Color',[.4 .4 .4],'LineWidth',1.5); hold on
plot(t,x(:,2),'Color',[.4 .4 .4],'LineWidth',1.5)
ph(1) = plot(tA,xA(:,1),'k','LineWidth',1.5);
plot(tA,xA(:,2),'k','LineWidth',1.5)
ph(2) = plot(tB,xB(:,1),'--','Color',ccolors(1,:),'LineWidth',1.5);
ph(3) = plot(tB,xB(:,2),'--','Color',ccolors(2,:),'LineWidth',1.5);
xlabel('Time')
ylabel('Population size')
l1=legend(ph,'Validation','SINDYc','SINDYc');
axis tight
ylim([-15 260])
set(gca,'xtick',[50,100,150,200])
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-dpdf', [figpath,'EX_LOTKA_ControlValidation2']);
print('-depsc2', [figpath,'EX_LOTKA_ControlValidation2']);

%% Apply Model predictive controlto system using SINDYc model
pest.ahat = Xi(:,1:2);
pest.polyorder = polyorder;
pest.usesine = usesine;
p.a = a; % True model parameters
p.b = b;
p.d = d;
p.g = g;

% Choose prediction horizon over which the optimization is performed
% Nvec = [1,3,5,7,10,15,20,25,30,35,40,45,50];
Nvec = [2,4,6,8,9,11,12,13];

for i = 1:length(Nvec)
    Ts          = 0.1;              % Sampling time
    N           = Nvec(i);          % Control / prediction horizon (number of iterations)
    Duration    = 100;              % Run control for 100 time units
    Nvar        = 2;
    Q           = [1 0];            % State weights
    R           = 0.5;%0.5;         % Control variation du weights
    Ru = 0.5;%0.01;                 % Control weights
    B = [0; 1];                     % Control vector (which state is controlled)
    C = eye(Nvar);                  % Measurement matrix
    D = 0;                          % Feedforward (none)
    x0n=x0';%[100; 50];             % Initial condition
    uopt0 = 0;                      % Set initial control input to zero

    % Constraints on control optimization
    LB = [];%-100*ones(N,1);        % Lower bound of control input
    UB = [];%100*ones(N,1);         % Upper bound of control input

    % Reference state, which shall be achieved
    xref1 = [g/d;a/b]; % critical point
    % xref1 = [50;0];      % Reference values
    % xref2 = [50;0];
    % xref_vec = [xref1(1)*ones(size(0:Ts:10)),xref2(1)*ones(size(10+Ts:Ts:Duration));
    %             xref1(2)*ones(size(0:Ts:10)),xref2(2)*ones(size(10+Ts:Ts:Duration))];

    % Options for optimization routine
    options = optimoptions('fmincon','Algorithm','sqp','Display','none');

    % Start simulation
    fprintf('Simulation started.  It might take a while...\n')
    x        = x0n;
    Ton      = 30;      % Time when control starts
    uopt     = uopt0.*ones(N,1);
    xHistory = x;       % Stores state history
    uHistory = uopt(1); % Stores control history
    tHistory = 0;       % Stores time history
    rHistory = xref1;   % Stores reference (could be trajectory and vary with time)
    tic
    for ct = 1:(Duration/Ts)   % For each iteration: take measurements & optimize control input & apply control input
        if ct*Ts>30            % Turn control on
            if ct*Ts==Ton+Ts
                disp('Start control.')
            end

            % Set references
            xref = xref1;

            % NMPC with full-state feedback
            COSTFUN = @(u) lotkaObjectiveFCN(u,x,Ts,N,xref,uopt(1),pest,diag(Q),R,Ru);
            CONSFUN = @(u) lotkaConstraintFCN(u,x,Ts,N,pest);
            uopt = fmincon(COSTFUN,uopt,[],[],[],[],LB,UB,CONSFUN,options);
            %  uopt = fmincon(COSTFUN,uopt,[],[],[],[],LB,UB,[],options);
            %  %use this without constraint functions CONSFUN

        else                    % If control is off
            uopt = uopt0.*ones(N,1);
            xref = [nan; nan];
        end

        % Integrate system: Apply control & Step one timestep forward
        x = rk4u(@lotkacontrol_discrete,x,uopt(1),Ts/1,1,[],p); %10, 2
        xHistory = [xHistory x];
        uHistory = [uHistory uopt(1)];
        tHistory = [tHistory tHistory(end)+Ts/1];
        rHistory = [rHistory xref];

    end
    fprintf('Simulation finished!\n')
    toc

    %% Show results
    clear ph

    figure;hold on, box on,
    ccolors = get(gca,'colororder');
    plot([Ton+tspan(1),Ton+tspan(1)],[-15 260],':','Color',[0.4,0.4,0.4],'LineWidth',1)
    text(31+tspan(1),210,'Control', 'FontSize',12)
    text(31+tspan(1),190,'turned on', 'FontSize',12)
    plot(tHistory+tspan(1),zeros(length(tHistory),1),'-k','LineWidth',0.5)
    plot(tHistory+tspan(1),xref1(1)*ones(length(tHistory),1),'--','Color',ccolors(1,:),'LineWidth',1)
    plot(tHistory+tspan(1),xref1(2)*ones(length(tHistory),1),'--','Color',ccolors(2,:),'LineWidth',1)
    ph(1) = plot(tHistory+tspan(1),xHistory(1,:),'-','Color',ccolors(1,:),'LineWidth',1.5);
    ph(2) = plot(tHistory+tspan(1),xHistory(2,:),'-','Color',ccolors(2,:),'LineWidth',1.5);
    ph(3) = plot(tHistory+tspan(1),uHistory,'-k','LineWidth',1.5);
    xlabel('Time')
    ylabel('Population size')
    legend(ph,'Prey','Predator','Control')
    axis tight
    ylim([-15 260])
    xlim([100,200.0001])
    %ylim([min(uHistory)-5 260])
    set(gca,'xtick',[50,100,150,200])
    set(gca,'LineWidth',1, 'FontSize',14)
    set(gcf,'Position',[100 100 300 200])
    set(gcf,'PaperPositionMode','auto')
    print('-depsc2', [figpath,'EX_LOTKA_DynamicsControlled_cnstrnd_N',num2str(N),'.eps']);

    %% Save Results
    Results.t = tHistory;
    Results.x = xHistory;
    Results.u = uHistory;
    Results.J = evalObjectiveFCN(uHistory,xHistory,rHistory,diag(Q),R,Ru);

    save(fullfile(datapath,['EX_LOTKA_MPC_SINDYc_N',num2str(N),'.mat']),'Results')

end
return

author: eurika-kaiser

repo: SINDY-MPC


text:
function segmentation = PTKGetSurfaceFromSegmentation(segmentation, exclude_direction)
    % PTKGetSurfaceFromSegmentation. Finds the surface of a segmented 3D binary volume.
    %
    % The input and output images are raw image matrices containing a mask of 0s
    % and 1s where 1 is a segmented voxel.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if nargin < 2
        exclude_direction = [];
    end

    segmentation = int8(segmentation > 0);

    % Find voxels that are nonzero and do not have 6 nonzero neighbours
    filter = zeros(3, 3, 3);
    if isempty(exclude_direction)
        filter(:) = [0 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 0 0 0 0];
    elseif exclude_direction == PTKImageOrientation.Coronal
        filter(:) = [0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0];
    else
        error('PTKGetSurfaceFromSegmentation:UnsupportedExcludeDirection', 'Unsupported exclude direction');
    end
    max_neighbours = sum(filter(:));
    image_conv = convn(segmentation, filter, 'same');
    segmentation = (segmentation & (image_conv < max_neighbours));
end

author: tomdoel

repo: pulmonarytoolkit


text:
function x = Advanced_SC_decoder(llr, frozen_bits, node_type_matrix, lambda_offset_)
%Theoretically, this decoder has much lower latency compared with
%traditional SC decoder because 9 kinds of specific constituent nodes are
%identified and decoded with alternative methods instead of successive
%cancellation. However, due to complex logical that is employed to
%distingguish above nodes, the decoding speed is even slightly lower than traditional SC in terms of MATLAB implementation.
N = length(llr);
m = log2(N);
C = zeros(2 * N - 1, 2);
P = zeros(2 * N - 1, 1);
P(end - N + 1 : end) = llr;
node_type_matrix_cnt = 1;
phi = 0;
while(phi < N)
    if (phi + 1) == node_type_matrix(node_type_matrix_cnt, 1)
        switch node_type_matrix(node_type_matrix_cnt, 3)
            case -1%RATE 0
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;%This psi leads you to a known (higher level with already obtained LLRs) layer.
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end

                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = zeros(M, 1);
                else
                    C(M : 2 * M - 1, 2) = zeros(M, 1);
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %
                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = zeros(M, 1);
                %                 else
                %                     C(M : 2 * M - 1, 2) = zeros(M, 1);
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 1%RATE 1
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_r1 = P(M : 2 * M - 1) < 0;
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_r1;
                else
                    C(M : 2 * M - 1, 2) = x_r1;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_r1 = P(M : 2 * M - 1) < 0;
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_r1;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_r1;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 2%REP
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_rep = (sum(P(M : 2 * M - 1)) < 0) * ones(M , 1);
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_rep;
                else
                    C(M : 2 * M - 1, 2) = x_rep;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_rep = (sum(P(M : 2 * M - 1)) < 0) * ones(M , 1);
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_rep;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_rep;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 3%spc
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_spc = SPC(P(M : 2 * M - 1));
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_spc;
                else
                    C(M : 2 * M - 1, 2) = x_spc;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_spc = SPC(P(M : 2 * M - 1));
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_spc;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_spc;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 4% I type
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_typeI = typeI(P(M : 2 * M - 1));
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_typeI;
                else
                    C(M : 2 * M - 1, 2) = x_typeI;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_typeI = typeI(P(M : 2 * M - 1));
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_typeI;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_typeI;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 5% II type

                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_typeII = typeII(P(M : 2 * M - 1));
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_typeII;
                else
                    C(M : 2 * M - 1, 2) = x_typeII;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_typeII = typeII(P(M : 2 * M - 1));
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_typeII;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_typeII;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 6% III type
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_typeIII = typeIII(P(M : 2 * M - 1));
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_typeIII;
                else
                    C(M : 2 * M - 1, 2) = x_typeIII;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_typeIII = typeIII(P(M : 2 * M - 1));
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_typeIII;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_typeIII;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 7% IV type
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_typeIV = typeIV(P(M : 2 * M - 1));
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_typeIV;
                else
                    C(M : 2 * M - 1, 2) = x_typeIV;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;

                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_typeIV = typeIV(P(M : 2 * M - 1));
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_typeIV;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_typeIV;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
            case 8% V type
                M = node_type_matrix(node_type_matrix_cnt, 2);
                reduced_layer = log2(M); %This reduced layer denotes where to stop.
                psi = phi;
                psi_ = phi;%This psi_ is used for bits recursion
                for i = 1 : reduced_layer
                    psi_ = floor(psi_/2);
                end
                layer = 0;%This layer denotes where to start.

                if phi ~= 0
                    while(mod(psi, 2) == 0)
                        psi = floor(psi/2);
                        layer = layer + 1;
                    end
                end

                if phi == 0
                    for i_layer = m - 1 : -1 : reduced_layer%NOT FULL LAYER = 0
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        for beta = 0 : index_1 - 1
                            P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                                min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                        end
                    end
                else
                    for i_layer = layer: -1 : reduced_layer
                        index_1 = lambda_offset_(i_layer + 1);
                        index_2 = lambda_offset_(i_layer + 2);
                        if i_layer == layer
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                            end
                        else
                            for beta = 0 : index_1 - 1
                                P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                                    min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                            end
                        end
                    end
                end
                x_typev = typeV(P(M : 2 * M - 1));
                if mod(psi_, 2) == 0
                    C(M : 2 * M - 1, 1) = x_typev;
                else
                    C(M : 2 * M - 1, 2) = x_typev;
                    while(mod(psi_, 2) == 1)
                        psi_ = floor(psi_/2);
                        index_1 = lambda_offset_(reduced_layer + 1);
                        index_2 = lambda_offset_(reduced_layer + 2);
                        for beta = 0 : index_1 - 1
                            C(2 * beta + index_2, 1 + mod(psi_, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                            C(2 * beta + 1 + index_2, 1 + mod(psi_, 2)) = C(beta + index_1, 2);
                        end
                        reduced_layer = reduced_layer + 1;
                    end
                end
                phi = phi + M;
                %                 M = node_type_matrix(node_type_matrix_cnt, 2);
                %                 reduced_layer = log2(M);
                %                 psi = phi;
                %                 for i = 1 : reduced_layer
                %                     psi = floor(psi/2);
                %                 end
                %                 P = recursive_calc_llr(m - reduced_layer, psi, P, C(: , 1), m, lambda_offset);
                %                 x_typeV = typeV(P(M : 2 * M - 1));
                %                 if mod(psi, 2) == 0
                %                     C(M : 2 * M - 1, 1) = x_typeV;
                %                 else
                %                     C(M : 2 * M - 1, 2) = x_typeV;
                %                     [C(:, 1), C(:, 2)] = recursive_calc_bits(m - reduced_layer, psi, C(:, 1), C(:, 2), m, lambda_offset);
                %                 end
                %                 phi = phi + M;
        end
        node_type_matrix_cnt = node_type_matrix_cnt + 1;
        if node_type_matrix_cnt == size(node_type_matrix, 1)
            node_type_matrix_cnt = 1;
        end
    else
        layer = 0;
        if phi == 0
            layer = m - 1;
        else
            psi = phi;
            while(mod(psi, 2) == 0)
                psi = floor(psi/2);
                layer = layer + 1;
            end
        end

        if phi == 0
            for i_layer = m - 1 : -1 : 0
                index_1 = lambda_offset_(i_layer + 1);
                index_2 = lambda_offset_(i_layer + 2);
                for beta = 0 : index_1 - 1
                    P(beta + index_1) = sign(P(2*beta + index_2)) *  sign(P(2*beta + 1 + index_2)) * ...
                        min(abs(P(2*beta + index_2)), abs(P(2*beta + 1 + index_2)));
                end
            end
        else
            for i_layer = layer: -1 : 0
                index_1 = lambda_offset_(i_layer + 1);
                index_2 = lambda_offset_(i_layer + 2);
                if i_layer == layer
                    for beta = 0 : index_1 - 1
                        P(beta + index_1) = (1 - 2*C(beta + index_1, 1)) * P(2 * beta + index_2) + P(2 * beta + 1 + index_2);
                    end
                else
                    for beta = 0 : index_1 - 1
                        P(beta + index_1) = sign(P(2 * beta + index_2)) *  sign(P(2 * beta + 1 + index_2)) * ...
                            min(abs(P(2 * beta + index_2)), abs(P(2 * beta + 1 + index_2)));
                    end
                end
            end
        end

        if frozen_bits(phi + 1) == 1
            C(1, 1 + mod(phi, 2)) = 0;
        else
            C(1, 1 + mod(phi, 2)) = P(1) < 0;
        end

        if mod(phi, 2) == 1
            layer = 0;
            psi = phi;
            while(mod(psi, 2) == 1)
                psi = floor(psi/2);
                index_1 = lambda_offset_(layer + 1);
                index_2 = lambda_offset_(layer + 2);
                for beta = 0 : index_1 - 1
                    C(2 * beta + index_2, 1 + mod(psi, 2)) = mod(C(beta + index_1, 1) + C(beta + index_1, 2), 2);
                    C(2 * beta + 1 + index_2, 1 + mod(psi, 2)) = C(beta + index_1, 2);
                end
                layer = layer + 1;
            end
        end

        %         P = recursive_calc_llr(m, phi, P, C(: , 1), m, lambda_offset);
        %         if frozen_bits(phi + 1) == 1
        %             C(1, mod(phi, 2) + 1) = 0;
        %         else
        %             C(1, mod(phi, 2) + 1) = P(1) < 0;
        %         end
        %
        %         if mod(phi, 2) == 1
        %             [C(:, 1), C(:, 2)] = recursive_calc_bits(m, phi, C(:, 1), C(:, 2), m, lambda_offset);
        %         end

        phi = phi + 1;
    end
end
x = C(end - N + 1 : end, 1);
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
%DEMO_AUDITORYFILTERBANK  Construct an auditory filterbank
%
%   In this file we construct a uniform filterbank using a the impulse
%   response of a 4th order gammatone for each channel. The center frequencies
%   are equidistantly spaced on an ERB-scale, and the width of the filter are
%   choosen to match the auditory filter bandwidth as determined by Moore.
%
%   Each channel is subsampled by a factor of 8 (a=8), and to generate a
%   nice plot, 4 channels per Erb has been used.
%
%   The filterbank covers only the positive frequencies, so we must use
%   |filterbankrealdual| and |filterbankrealbounds|.
%
%   .. figure::
%
%      Classic spectrogram
%
%      A classic spectrogram of the spoken sentense. The dynamic range has
%      been set to 50 dB, to highlight the most important features.
%
%   .. figure::
%
%      Auditory filterbank representation
%
%      Auditory filterbank representation of the spoken sentense using
%      gammatone filters on an Erb scale.  The dynamic range has been set to
%      50 dB, to highlight the most important features.
%
%   See also: freqtoaud, audfiltbw, gammatonefir, ufilterbank, filterbankrealdual
%
%   References: glasberg1990daf

% Use part of the 'cocktailparty' spoken sentense.
f=cocktailparty;
f=f(20001:80000,:);
fs=44100;
a=8;
channels_per_erb=2;
filterlength=5000;
dynrange_for_plotting=50;

% Determine minimal transform length
Ls=length(f);
L=ceil(filterlength/a)*a;

% Number of channels, slightly less than 1 ERB(Cambridge) per channel.
M=ceil(freqtoerb(fs/2)*channels_per_erb);

% Compute center frequencies.
fc=erbspace(0,fs/2,M);

%% --------------- display classic spectrogram -------------------
figure(1);
sgram(f,fs,dynrange_for_plotting);

%% ---------------  Gammatone filters ----------------------------

if 1

g_gam=gammatonefir(fc,fs,filterlength,'peakphase');

% In production code, it is not necessary to call 'filterbankrealbounds',
% this is just for veryfying the setup.
disp('Frame bound ratio for gammatone filterbank, should be close to 1 if the filters are choosen correctly.');
filterbankrealbounds(g_gam,a,L)

% Create reconstruction filters
gd_gam=filterbankrealdual(g_gam,a,L);

% Analysis transform
coef_gam=ufilterbank(f,g_gam,a);

% Synthesis transform
r_gam=2*real(ifilterbank(coef_gam,gd_gam,a,Ls));

disp('Relative error in reconstruction, should be close to zero.');
norm(f-r_gam)/norm(f)

figure(2);
plotfilterbank(coef_gam,a,fc,fs,dynrange_for_plotting,'audtick');

F  = frame('ufilterbankreal',g_gam,a,M);
c2 = frana(F,f);
Ls=length(f);

[r_iter,relres,iter] = frsyniter(F,c2,Ls);
disp('Relative error in interative reconstruction, should be close to zero.');
norm(f-r_iter)/norm(f)

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function mrGrayXformRoi(roiFile, roiXform, newRoiFile)
%
%  mrGrayXformRoi(roiFile, roiXform, newRoiFile)
%
% Convert a mrVista ROI from one vAnatomy space to another. This is
% usually used to convert an ROI defined in an old vAnatomy file to one
% compatible with a new NIFTI anatomy file.
%
% See also:
%   mrGrayConvertClassToNifti
%
% HISTORY:
% 2008.09.05 RFD wrote it.

%% Get input file names
opts = {'*.mat','mrVista ROI'; '*.*','All Files (*.*)'};
if(~exist('roiFile','var'))
    [f, p]=uigetfile(opts, 'Pick a mrVista ROI file');
    if(isequal(f,0)|| isequal(p,0)) , disp('user canceled.'); return; end
    if(isempty(p)), p = pwd; end
    roiFile = fullFile(p,f);
end

if(~exist('roiXform','var')||isempty(roiXform))
    default = fullfile(fileparts(roiFile),'t1_xformToVanat.mat');
    [f, p]=uigetfile(opts, 'Pick an ROI xform file',default);
    if(isequal(f,0)|| isequal(p,0)), disp('user canceled.'); return; end
    roiXform = fullFile(p,f);
end

if(~exist('newRoiFile','var')||isempty(newRoiFile))
    [p,f,e] = fileparts(roiFile);
    default = fullfile(p,[f '_NEW' e]);
    [f,p] = uiputfile('*.mat','Save new ROI file as...',default);
    if(isequal(f,0)|| isequal(p,0)), disp('user canceled.'); return; end
    newRoiFile = fullfile(p,f);
end
disp(['New ROI will be saved in ' newRoiFile]);

oldRoi = load(roiFile);
roiXform = load(roiXform);
oldImg = zeros(roiXform.vAnatSize);
oldInds = sub2ind(roiXform.vAnatSize, oldRoi.ROI.coords(1,:), oldRoi.ROI.coords(2,:), oldRoi.ROI.coords(3,:));
oldImg(oldInds) = 1;
newImg = mrAnatResliceSpm(double(oldImg), roiXform.xform, roiXform.bb, roiXform.mm, [1 1 1 0 0 0], false);
newImg = mrAnatRotateAnalyze(newImg);
newInds = find(newImg>=0.5);
[x,y,z] = ind2sub(roiXform.sz,newInds);

ROI = oldRoi.ROI;
ROI.coords = [x'; y'; z'];

save(newRoiFile,'ROI');

return;

% % For debugging:
% vAnat = readVolAnat('OLD/vAnatomy.dat');
% r = vAnat; g=r; b=r;
% r(oldInds) = 255;
% showMontage(cat(4,r,g,b)./255);
%
% vAnatNew = readVolAnat('t1.nii.gz');
% r = vAnatNew; g=r; b=r;
% r(newInds) = 255;
% showMontage(cat(4,r,g,b)./255);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function x = blend_rst_2dn ( r, s, t, n, bound_rst )

%*****************************************************************************80
%
%% BLEND_RST_2DN extends vector data on faces into a cube.
%
%  Diagram:
%
%    010-----r10-----110        011-----r11-----111
%      |       .       |          |       .       |
%      |       .       |          |       .       |
%    0s0.....rs0.....1s0        0s1.....rs1.....1s1     S
%      |       .       |          |       .       |     |
%      |       .       |          |       .       |     |
%    000-----r00-----100        001-----r01-----101     +----R
%           BOTTOM                      TOP
%
%    011-----0s1-----001        111-----1s1-----101
%      |       .       |          |       .       |
%      |       .       |          |       .       |
%    01t.....0st.....00t        11t.....1st.....10t          T
%      |       .       |          |       .       |          |
%      |       .       |          |       .       |          |
%    010-----0s0-----000        110-----1s0-----100     S----+
%           LEFT                       RIGHT
%
%    001-----r01-----101        011-----r11-----111
%      |       .       |          |       .       |
%      |       .       |          |       .       |
%    00t.....r0t.....100        01t.....r1t.....11t     T
%      |       .       |          |       .       |     |
%      |       .       |          |       .       |     |
%    000-----r00-----100        010-----r10-----110     +----R
%           FRONT                       BACK
%
%  Discussion:
%
%    BLEND_RST_2DN is NOT equivalent to a trilinear finite element
%    method, since the data is sampled everywhere along the corners,
%    edges, and faces, rather than at a finite number of nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Gordon,
%    Blending-Function Methods of Bivariate and Multivariate Interpolation
%    and Approximation,
%    SIAM Journal on Numerical Analysis,
%    Volume 8, Number 1, March 1971, pages 158-177.
%
%    William Gordon and Charles Hall,
%    Transfinite Element Methods: Blending-Function Interpolation over
%    Arbitrary Curved Element Domains,
%    Numerische Mathematik,
%    Volume 21, Number 1, 1973, pages 109-129.
%
%    William Gordon and Charles Hall,
%    Construction of Curvilinear Coordinate Systems and Application to
%    Mesh Generation,
%    International Journal of Numerical Methods in Engineering,
%    Volume 7, 1973, pages 461-477.
%
%    Joe Thompson, Bharat Soni, Nigel Weatherill,
%    Handbook of Grid Generation,
%    CRC Press, 1999.
%
%  Parameters:
%
%    Input, real R, S, T, the (R,S,T) coordinates of the point
%    to be evaluated.
%
%    Input, integer N, the dimension of the vector space.
%
%    External, BOUND_RST, is a function which is given (R,S,T)
%    coordinates and an component value I, and returns XI, the value
%    of the I-th component of the N-vector at that point.  BOUND_RST
%    will only be called for "faces", that is, for values (R,S,T) where
%    at least one of R, S and T is either 0.0 or 1.0.  BOUND_RST has
%    the form:
%      function xi = bound_rst ( r, s, t, i )
%
%    Output, real X(N), the interpolated value at the
%    point (R,S,T).
%
  for i = 1 : n
%
%  Get the I-th coordinate component at the corners.
%
    x000 = bound_rst ( 0.0, 0.0, 0.0, i );
    x001 = bound_rst ( 0.0, 0.0, 1.0, i );
    x010 = bound_rst ( 0.0, 1.0, 0.0, i );
    x011 = bound_rst ( 0.0, 1.0, 1.0, i );
    x100 = bound_rst ( 1.0, 0.0, 0.0, i );
    x101 = bound_rst ( 1.0, 0.0, 1.0, i );
    x110 = bound_rst ( 1.0, 1.0, 0.0, i );
    x111 = bound_rst ( 1.0, 1.0, 1.0, i );
%
%  Get the I-th coordinate component at the edges.
%
    xr00 = bound_rst ( r, 0.0, 0.0, i );
    xr01 = bound_rst ( r, 0.0, 1.0, i );
    xr10 = bound_rst ( r, 1.0, 0.0, i );
    xr11 = bound_rst ( r, 1.0, 1.0, i );

    x0s0 = bound_rst ( 0.0, s, 0.0, i );
    x0s1 = bound_rst ( 0.0, s, 1.0, i );
    x1s0 = bound_rst ( 1.0, s, 0.0, i );
    x1s1 = bound_rst ( 1.0, s, 1.0, i );

    x00t = bound_rst ( 0.0, 0.0, t, i );
    x01t = bound_rst ( 0.0, 1.0, t, i );
    x10t = bound_rst ( 1.0, 0.0, t, i );
    x11t = bound_rst ( 1.0, 1.0, t, i );
%
%  Get the I-th component on the faces.
%
    x0st = bound_rst ( 0.0, s, t, i );
    x1st = bound_rst ( 1.0, s, t, i );
    xr0t = bound_rst ( r, 0.0, t, i );
    xr1t = bound_rst ( r, 1.0, t, i );
    xrs0 = bound_rst ( r, s, 0.0, i );
    xrs1 = bound_rst ( r, s, 1.0, i );
%
%  Interpolate the I-th coordinate component of the interior point.
%
    x(i) = blend_123 ( r, s, t, x000, x001, x010, x011, x100, x101, x110, x111, ...
      xr00, xr01, xr10, xr11, x0s0, x0s1, x1s0, x1s1, x00t, x01t, x10t, x11t, ...
      x0st, x1st, xr0t, xr1t, xrs0, xrs1 );

  end

  return
end


text:
function [C, d] = spot_ellipsoid(A, b)

n = size(A, 2);

% tic
pr = spotsosprog;
[pr, t] = pr.newFree(1);
[pr, Z] = pr.newFree(n, n);
[pr, C] = pr.newSym(n);
[pr, d] = pr.newFree(n);
z = diag(Z);
pr = pr.withPSD([C, Z'; Z, diag(z)]);
lor = [b - A * d, A * C']';
pr = pr.withLor(lor);

l = ceil(log2(n));
m = 2^l - n;
if m > 0
  x = [z; repmat(t, m, 1)];
else
  x = z;
end
[pr, g] = geo_mean_recursive(pr, x);
pr = pr.withEqs(2^(l/2) * t - g);

% fprintf(1, 'setup: %f\n', toc);

% tic
% disp('mosek')
% solver = @spot_mosek;
% sol = pr.minimize(-t, solver);
% sol.eval(lor)
% double(sol.eval(C))
% double(sol.eval(d))
% fprintf(1, 'mosek: %f\n', toc);

% tic
solver = @spot_sedumi;
sol = pr.minimize(-t, solver);
% sol.eval(lor)
% fprintf(1, 'sedumi: %f\n', toc);

% tic
C = double(sol.eval(C));
d = double(sol.eval(d));
% fprintf(1, 'extract: %f\n', toc);

end

function [pr, g] = geo_mean_recursive(pr, x)
  n = length(x);
  if n > 1
    [pr, y] = pr.newFree(n/2);
    pr = pr.withRLor([reshape(x, 2, n/2); reshape(y, 1, n/2)]);
    [pr, g] = geo_mean_recursive(pr, y);
  else
    g = x;
  end
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function z = round2(x,y)
%ROUND2 rounds number to nearest multiple of arbitrary precision.
%   Z = ROUND2(X,Y) rounds X to nearest multiple of Y.
%
%Example 1: round PI to 2 decimal places
%   >> round2(pi,0.01)
%   ans =
%         3.14
%
%Example 2: round PI to 4 decimal places
%   >> round2(pi,1e-4)
%   ans =
%         3.1416
%
%Example 3: round PI to 8-bit fraction
%   >> round2(pi,2^-8)
%   ans =
%         3.1406
%
%Examples 4-6: round PI to other multiples
%   >> round2(pi,0.05)
%   ans =
%         3.15
%   >> round2(pi,2)
%   ans =
%         4
%   >> round2(pi,5)
%   ans =
%         5
%
% See also ROUND.

%% defensive programming
error(nargchk(2,2,nargin))
error(nargoutchk(0,1,nargout))
if numel(y)>1
  error('Y must be scalar')
end

%%
z = round(x/y)*y;

author: Sable

repo: mcbench-benchmarks


text:
%DEMO_LOCSTAT - performing the image enhancement based on local statistics.
% Using NLFILTER function to save memory.
%Copyright by Alex Zuo (zuoxinian@yahoo.com.cn).

clear,clc;
%% Read original SEM images
Isem=imread('testimages/sem.jpg');
%% Caculate the global mean and variance.
f=double(Isem);
M=mean2(f);D=std2(f);
%% Perform the enhancing task.
Bsize=[3 3];
k=[0.4 0.02 0.4];
E=4.0;
tic
fsem_enh=nlfilter(f,Bsize,@mylocstat,M,D,E,k);
t_nlfilter=toc %display time-consuming
Isem_enh=im2uint8(mat2gray(fsem_enh));
%% Visualize the result.
subplot(1,2,1),imshow(Isem),title('SEM Image')
subplot(1,2,2),imshow(Isem_enh),title('Enhanced SEM Image with NLFILTER')
%Save the result for publish.
%imwrite(Isem_enh,'outfiles/sem_enhanced.jpg')
%print(gcf,'-deps','outfiles/sem_locstatenh.eps')
%figure,imshow(Isem)
%print(gcf,'-deps','outfiles/sem.eps')