text:
function [d,S] = solve_for_d_S(Y,A,R,options)
[N,B] = size(Y);
M = size(A,2);

sigma0 = 0.1; % initial sigma
sigma_max = 1;
sigma_min = 1e-9;
delta_t0 = 1e-9;

disp('Start estimating D and S');
C = kron((A'*A),eye(B));
YAR = (Y-A*R)';
Z = YAR*A;
E = diag(sum(YAR.^2,2));
F1 = YAR*A;
F = zeros(M*B,B);
for i = 1:M
    F((i-1)*B+1:i*B,:) = diag(F1(:,i));
end

d = (1/N)*sum(YAR.^2,2);
d = 1./sqrt(d);

S = eye(M*B,M*B);
if 0
    u = 1/sqrt(B) * ones(B,1);
    S1 = sigma0^2*u*u' + 1e-4*eye(B);
    S1 = inv(S1);
else
    S1 = sigma0^(-2)*eye(B);
end

for j = 1:M
    inds = (j-1)*B+1:j*B;
    S(inds,inds) = diag(1./d) * S1 * diag(1./d);
end

errors2 = eval_obj_fun_D_S(S,d,C,YAR,Z,M,B,N);
for iter = 1:300
    % solve for S
    z = diag(d)*Z;
    z = z(:);
    if 0
        obj_fun_S = @(S) eval_obj_fun_S(S,C,z,M,B);
        fun_der_S = @(S) calc_der_S(S,C,z,M,B);
        proj_fun_S = @(S) project_to_spd(S,d,sigma_min,sigma_max,M,B);

        S = projected_gradient_descent(S, obj_fun_S, fun_der_S, ...
            proj_fun_S, delta_t0);
    else
        if 1
            proj_fun_S = @(S) project_to_spd(S,d,sigma_min,sigma_max,1,B);
            for j = 1:M
                obj_fun_Sj = @(Sj) eval_obj_fun_Sj(Sj,S,C,z,M,B,j);
                fun_der_Sj = @(Sj) calc_der_Sj(Sj,S,C,z,B,j);
                inds = (j-1)*B+1:j*B;
                Sj = S(inds,inds);
                Sj = projected_gradient_descent(Sj, obj_fun_Sj, fun_der_Sj, ...
                    proj_fun_S, 1e-5);
                S(inds,inds) = Sj;
            end
        else
            der_S = calc_der_S(S,C,z,M,B);
            proj_fun_S = @(S) project_to_spd(S,d,sigma_min,sigma_max,1,B);
            for j = 1:M
                obj_fun_Sj = @(Sj) eval_obj_fun_Sj(Sj,S,C,z,M,B,j);
                fun_der_Sj = @(Sj) calc_der_Sj_approx(Sj,der_S,B,j);
                inds = (j-1)*B+1:j*B;
                Sj = S(inds,inds);
                Sj = projected_gradient_descent(Sj, obj_fun_Sj, fun_der_Sj, ...
                    proj_fun_S, 1e-3);
                S(inds,inds) = Sj;
            end
        end
    end

    % solve for D
    Q = S + C;
    G = (1/N)*(E-F'*(Q\F));
    fcn_d = @(d) G*d - 1./d;
    fsolve_opts = optimset('Display','off');
    d = fsolve(fcn_d, 1./sqrt((1/N)*diag(E)), fsolve_opts);

    % calc error
    err = eval_obj_fun_D_S(S,d,C,YAR,Z,M,B,N);
    errors2 = [errors2;err];
    if test_convergence(errors2, 1e-8), break; end

    if mod(iter,10) == 0
        disp(['Process iteration ',num2str(iter),'. The objective ',...
            'function has value ',num2str(err)]);
    end
end

function S = project_to_spd(S,d,sigma_min,sigma_max,M,B)
touch_bd = 0;
S = (S + S')/2;
for j = 1:M
    inds = (j-1)*B+1:j*B;
    S3 = S(inds,inds);
    SigmaInv = diag(d) * S3 * diag(d);
    [V,D] = eig(SigmaInv);
    d1 = diag(D);
    low_bd = 1/(sigma_max^2);
    up_bd = 1/(sigma_min^2);

    if ~isempty(find(d1<low_bd, 1))
        touch_bd = touch_bd + 1;
%         break;
    end

    d1(d1<low_bd) = low_bd;
    d1(d1>up_bd) = up_bd;
    S(inds,inds) = diag(1./d)*V*diag(d1)*V'*diag(1./d);
end

if touch_bd >= 1
%     disp(['Uncertainty amount upper bound touched ',num2str(touch_bd),' times']);

    % Once uncertainty amount upper bound touched. It means the change of
    % the covariance is too large. Setting S to 0 will make the objective
    % function be Inf, thus revoke the delta_t.

%     if touch_bd >= 1
%         S = 0 * eye(M*B);
%     end
end

function val = eval_obj_fun_Sj(Sj,S,C,z,M,B,j)
inds = (j-1)*B+1:j*B;
S(inds,inds) = Sj;
val = eval_obj_fun_S(S,C,z,M,B);

function val = eval_obj_fun_S(S,C,z,M,B)
Q = S + C;
% R = chol(Q);
val1 = 0;
for j = 1:M
    inds = (j-1)*B+1:j*B;
    val1 = val1 + logdet(S(inds,inds));
end

% val = -z'*(Q\z) + logdet(Q) - val1;
R = chol(Q);
y = R'\z;
val = -sum(y.^2) + 2*sum(log(diag(R))) - val1;

function val = eval_obj_fun_D_S(S,d,C,YAR,Z,M,B,N)
z = diag(d)*Z;
z = z(:);

val = eval_obj_fun_S(S,C,z,M,B);
val = val + sum(sum((diag(d)*YAR).^2)) - 2*N*sum(log(d));

function der = calc_der_S(S,C,z,M,B)
Q = S + C;
invQ = inv(Q);
y = invQ*z;
T = y*y';
der = zeros(M*B,M*B);
for j = 1:M
    inds = (j-1)*B+1:j*B;
    Sj = S(inds,inds);
    der(inds,inds) = T(inds,inds) + invQ(inds,inds) - inv(Sj);
end

% der_norm = norm(der,'fro');
% der = der/der_norm;

function der = calc_der_Sj(Sj,S,C,z,B,j)
inds = (j-1)*B+1:j*B;
S(inds,inds) = Sj;

Q = S + C;
invQ = inv(Q);
y = invQ*z;
T = y*y';
der = zeros(B,B);

der = T(inds,inds) + invQ(inds,inds) - inv(Sj);

der_norm = norm(der,'fro');
der = der/der_norm;

function der = calc_der_Sj_approx(Sj,der_S,B,j)
inds = (j-1)*B+1:j*B;
der = der_S(inds,inds);

der_norm = norm(der,'fro');
der = der/der_norm;

function d = solve_for_d1(Y,A,R,S)
[N,B] = size(Y);
M = size(A,2);

YAR = Y - A*R;
E = diag(sum(YAR.^2,1));

Q = S + kron((A'*A),eye(B));

F1 = A'*YAR;
F = zeros(M*B,B);
for i = 1:M
    F((i-1)*B+1:i*B,:) = diag(F1(i,:));
end
G = (1/N)*(E-F'*(Q\F));

fcn_d = @(d) G*d - 1./d;
d = fsolve(fcn_d, 1./sqrt((1/N)*diag(E)));

author: zhouyuanzxcv

repo: Hyperspectral


text:
function report = edl_readreport(filename)
%EDL_READREPORT Reads in a EDL report
%
% $ Syntax $
%   - report = edl_readreport(filename)
%
% $ Description $
%   - report = edl_readreport(filename) reads in a EDL report from
%     a report xml file. The returned report is a struct with
%     following fields:
%       - attribs: the header attributes
%           - guid:     the GUID string identifying the report
%       - entries:  the report item entries
%         at least have following fields:
%           - internal_index:  the internal index
%
% $ History $
%   - Created by Dahua Lin, on Aug 14, 2006
%

%% Read in file

doctag = 'ExpReport';
nodetag = 'Entry';

S = edl_readprops(filename, nodetag);

%% Post-Processing

% doc tag
if ~strcmp(S.tag, doctag)
    error('edl:parseerror', ...
        'Invalid document tag %s for report', S.tag);
end

% doc attribs
if isempty(S.attribs)
    error('edl:parseerror', ...
        'The document element for report has no attributes');
end

attrnames = {'guid'};
tf = isfield(S.attribs, attrnames);
if ~all(tf)
    error('edl:parserror', ...
        'The required header %s does not exist', ...
        attrnames{find(~tf, 1)});
end

report.attribs = S.attribs;

% entries

report.entries = edl_check_internalindices(S.(nodetag));

author: lmthang


text:
% Root class for composite user controls. More form than function, it seems...

% Copyright 2008-2009 Levente Hunyadi
classdef UIControl < hgsetget
    properties (Abstract, Dependent)
        % The graphics handle that represents the entire user control.
        % This handle is used in setting the parent, position, measurement
        % unit, etc. properties of the control. In order to use MatLab's Handle
        % Graphics semantics in a derived user interface control, you should
        % query the Control property.
        Control;
    end
    properties (Dependent)
        NormalizedPosition;
        % The parent control of this UI control.
        Parent;
        % Units in which the Position property is expressed.
        % Note: this property is deliberately put before the Position property
        % to ensure proper initialization when both Units and Position are set
        % in a constructor argument list. For the constructor processes
        % arguments in the order in which they occur within the class
        % definition rather than the order in which they are specified, we
        % always get the desired behavior that Position is interpreted in terms
        % of the specified Units.
        Units = 'normalized';
        Position;
    end
    methods (Abstract)
        % Forces instantiation of underlying MatLab Handle Graphics objects.
        % Some dependent properties refer to MatLab Handle Graphics object
        % properties in an early phase of object construction. These properties
        % might fail if existence of objects they depend on cannot be
        % guaranteed.
        obj = Instantiate(obj, parent);
    end
    methods
        function obj = UIControl(arg, varargin)
            if nargin > 0 && isscalar(arg) && ishandle(arg)
                obj = obj.Instantiate(arg);
            else
                obj = obj.Instantiate();
            end
            if nargin > 0
                if isscalar(arg) && ishandle(arg)
                    obj.Parent = arg;
                    obj = constructor(obj, varargin{:});
                else
                    if ~strcmp('figure', get(obj.Control, 'Type'))
                        obj.Parent = gcf;
                    end
                    obj = constructor(obj, arg, varargin{:});
                end
            end
        end

        function parent = get.Parent(obj)
            parent = get(obj.Control, 'Parent');
        end

        function obj = set.Parent(obj, parent)
            if ~strcmp(get(obj.Control, 'Type'), 'figure')  % figures do not require a parent
                set(obj.Control, 'Parent', parent);
            end
        end

        function units = get.Units(obj)
            units = get(obj.Control, 'Units');
        end

        function obj = set.Units(obj, units)
            set(obj.Control, 'Units', units);
        end

        function position = get.Position(obj)
            position = get(obj.Control, 'Position');
        end

        function obj = set.Position(obj, position)
            set(obj.Control, 'Position', position);
        end

        function obj = set.NormalizedPosition(obj, position)
            units = get(obj.Control, 'Units');
            if isempty(units)
                units = 'normalized';
            end
            switch units
                case 'normalized'
                    set(obj.Control, 'Position', position);
                otherwise
                    set(obj.Control, 'Units', 'normalized');
                    set(obj.Control, 'Position', position);
                    set(obj.Control, 'Units', units);
            end
        end

        function position = get.NormalizedPosition(obj)
            units = get(obj.Control, 'Units');
            if isempty(units)
                units = 'normalized';
            end
            switch units
                case 'normalized'
                    position = get(obj.Control, 'Position');
                otherwise
                    set(obj.Control, 'Units', 'normalized');
                    position = get(obj.Control, 'Position');
                    set(obj.Control, 'Units', units);
            end
        end

        function position = getpixelposition(obj)
            position = getpixelposition(obj.Control);
        end

        function obj = setpixelposition(obj, position)
            setpixelposition(obj.Control, position);
        end
    end
end

repo: REST


text:
function [param] = tps_compute_param(PP,kernel,U,Pm,Q1,Q2,R,lambda,target)
%%=====================================================================
%% $RCSfile: tps_compute_param.m,v $
%% $Author: bjian $
%% $Date: 2008/11/24 08:59:01 $
%% $Revision: 1.1 $
%%=====================================================================

TB = U*PP;
QQ = Q2*Q2';
A = inv(TB'*QQ*TB + lambda*kernel)*TB'*QQ;
tps = A*target;
affine = inv(R)*Q1'*(target-TB*tps);
param = [affine; tps];

author: Sable

repo: mcbench-benchmarks


text:
function obj_size_print ( input_file_name, node_num, face_num, normal_num, ...
  order_max )

%*****************************************************************************80
%
%% OBJ_SIZE_PRINT prints sizes associated with an OBJ file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Input, integer NODE_NUM, the number of vertices defined.
%
%    Input, integer FACE_NUM, the number of faces defined.
%
%    Input, integer NORMAL_NUM, the number of normal vectors defined.
%
%    Input, integer ORDER_MAX, the maximum number of vertices per face.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Object sizes for OBJ file "%s":\n', input_file_name );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Nodes              = %d\n', node_num );
  fprintf ( 1, '  Faces              = %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', order_max );
  fprintf ( 1, '  Normal vectors     = %d\n', normal_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = quad_fun ( n )

%*****************************************************************************80
%
%% QUAD_FUN demonstrates MATLAB's PARFOR command for parallel programming.
%
%  Discussion:
%
%    This function estimates an integral using the composite trapezoidal rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points to use.
%
%    Output, real VALUE, the estimate for the integral.
%
  a = 0.0;
  b = 1.0;

  value = 0.0;
  h = ( b - a ) / ( n - 1 );

  parfor i = 1 : n

    x = ( ( n - i ) * a + ( i - 1 ) * b ) / ( n - 1 );

    fx = f ( x );

    if ( i == 1 )
      value = value + 0.5 * fx * h
    elseif ( i < n )
      value = value +       fx * h;
    elseif ( i == n )
      value = value + 0.5 * fx * h;
    end

  end

  return
end
function value = f ( x )

%*****************************************************************************80
%
%% F is the function to be integrated.
%
%  Discussion:
%
%    The integral of F(X) from 0 to 1 is exactly PI.
%
%  Modified:
%
%    17 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the values where the integrand is to be evaluated.
%
%    Output, real VALUE, the integrand values.
%
  value = 4.0 ./ ( 1 + x.^2 );

  return
end


text:
function a = poisson ( nrow, ncol )

%*****************************************************************************80
%
%% POISSON returns the POISSON matrix.
%
%  Formula:
%
%    if ( I = J )
%      A(I,J) = 4.0
%    elseif ( I = J+1 or I = J-1 or I = J+NROW or I = J-NROW )
%      A(I,J) = -1.0
%    else
%      A(I,J) = 0.0
%
%  Example:
%
%    NROW = NCOL = 3
%
%     4 -1  0 | -1  0  0 |  0  0  0
%    -1  4 -1 |  0 -1  0 |  0  0  0
%     0 -1  4 |  0  0 -1 |  0  0  0
%     ----------------------------
%    -1  0  0 |  4 -1  0 | -1  0  0
%     0 -1  0 | -1  4 -1 |  0 -1  0
%     0  0 -1 |  0 -1  4 |  0  0 -1
%     ----------------------------
%     0  0  0 | -1  0  0 |  4 -1  0
%     0  0  0 |  0 -1  0 | -1  4 -1
%     0  0  0 |  0  0 -1 |  0 -1  4
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A results from discretizing Poisson's equation with the
%    5 point operator on a square mesh of N points.
%
%    A has eigenvalues
%
%      LAMBDA(I,J) = 4 - 2 * COS(I*PI/(N+1))
%                      - 2 * COS(J*PI/(M+1)), I = 1 to N, J = 1 to M.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub, Charles Van Loan,
%    Matrix Computations, second edition,
%    Johns Hopkins University Press, Baltimore, Maryland, 1989
%    (Section 4.5.4).
%
%  Parameters:
%
%    Input, integer NROW, NCOL, the number of rows and columns
%    in the grid.
%
%    Output, real A(NROW*NCOL,NROW*NCOL), the matrix.
%
  n = nrow * ncol;

  a = zeros ( n, n );

  i = 0;

  for i1 = 1 : nrow
    for j1 = 1 : ncol

      i = i + 1;

      if ( 1 < i1 )
        j = i - ncol;
        a(i,j) = -1.0;
      end

      if ( 1 < j1 )
        j = i - 1;
        a(i,j) = -1.0;
      end

      j = i;
      a(i,j) = 4.0;

      if ( j1 < ncol )
        j = i + 1;
        a(i,j) = -1.0;
      end

      if ( i1 < nrow )
        j = i + ncol;
        a(i,j) = -1.0;
      end

    end
  end

  return
end


text:
function title = p01_title ( )

%*****************************************************************************80
%
%% P01_TITLE returns a title for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'Simple quadratic, (x-2)^2+1.';

  return
end


text:
function forwardTransportQuantReverseFigures(model)
% Figure of the qualitatively forward transport reactions that are quantitatively reversible.
% Creates a vertical errorbar figure of the qualitatively forward transport reactions
% that are quantitatively reversible, whether from group contribution or
% Keq, that then need to be assigned to be forward, to limit the growth rate
% i.e. abc transporters or reactions involving protons.
%
% USAGE:
%
%    forwardTransportQuantReverseFigures(model)
%
% INPUT:
%    model:    structure with fields:
%
%                * .transportRxnBool
%                * .dGt0Min
%                * .dGt0Max
%                * .directions.forwardReversible
%
% .. Author: - Ronan M.T. Fleming

[nMet,nRxn]=size(model.S);

directions=model.directions;

%transport reactions
forwardTransportQuantReverseBool=false(nRxn,1);
for n=1:nRxn
    %if reaction directionality changed from forward to Reversible and als
    %a transport reaction
    if directions.forwardReversible(n) && model.transportRxnBool(n)
        %abc transporters or reactions involving protons
        if ~isempty(strfind(model.rxns{n},'abc')) ||  nnz(strcmp('h[c]',model.mets(model.S(:,n)~=0)))~=0 || nnz(strcmp('h[e]',model.mets(model.S(:,n)~=0)))~=0
            forwardTransportQuantReverseBool(n,1)=1;
        end
    end
end

if 0
    X1=1:nRxn;
    %dGrt0
    Y0=(model.DrGt0Min+model.DrGt0Max)/2;
    L0=Y0-model.DrGt0Min;
    U0=model.DrGt0Max-Y0;
    %dGrt
    Y=(model.DrdGtMin+model.DrGtMax)/2;
    L=Y-model.DrGtMin;
    U=model.DrGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
else
    Y0=model.DrGt0Mean;
    L0=model.DrGt0Mean-model.DrGt0Min;
    U0=model.DrGt0Max-model.DrGt0Mean;

    Y=model.DrGtMean;
    L=model.DrGtMean-model.DrGtMin;
    U=model.DrGtMax-model.DrGtMean;
    P=directions.forwardProbability;
end
%sort by probability that a reaction is forward (puts any NaN first)
[tmp,xip]=sort(P,'descend');

%     only take the indices of the problematic reactions, but be sure to
%     take them in order of descending P
xip2=zeros(nnz(forwardTransportQuantReverseBool),1);
p=1;
for n=1:nRxn
    if forwardTransportQuantReverseBool(xip(n))
        xip2(p)=xip(n);
        p=p+1;
    end
end
xip=xip2;
X1=1:length(xip);

%replace the NaN due to zero st dev
forwardReversible_dGf_dG0Fwd=false(nRxn,1);
forwardReversible_dGf_dG0Rev=false(nRxn,1);
for n=1:nRxn
    %only the transport reactions
    if forwardTransportQuantReverseBool(n)
        if model.DrGt0Min(n)==model.DrGt0Max(n)
            if model.DrGt0Min(n)<0
                forwardReversible_dGf_dG0Fwd(n)=1;
            else
                forwardReversible_dGf_dG0Rev(n)=1;
            end
        end
    end
end
nNaNpLHS=nnz(forwardReversible_dGf_dG0Fwd);
nNaNpRHS=nnz( forwardReversible_dGf_dG0Rev);
if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(forwardTransportQuantReverseBool))))
    warning('ExtraCategory');
end
%nans are first in the ordering of indexes
NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
%sorts indices of the zero std dev met by their mean dG0t
[tmp,xipNaNPInd]=sort(Y0(NaNPInd));
%new ordering
xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];

% close all
figure1 = figure;
% Create axes
axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
hold on;

%upper and lower X
minX=min(model.DrGtMin(forwardTransportQuantReverseBool));
maxX=max(model.DrGtMax(forwardTransportQuantReverseBool));
%baselines
Y1(1:length(xip))=minX;
P(P==1)=NaN;
P(P==0.5)=NaN;
Y1(~isnan(P(xip)) & P(xip)~=1 & P(xip)~=0.5)=maxX;
bar_handle6=barh(X1,Y1,1,'BaseValue',minX,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

%dGrt errorbar
%   PLOTERR(X,Y,{LX,UX},{LY,UY}) plots the graph of vector X vs. vector Y
%   with error bars specified by LX and UX in horizontal direction and
%   with error bars specified by LY and UY in vertical direction.
%   H = PLOTERR(...) returns a vector of line handles in this order:
%      H(1) = handle to datapoints
%      H(2) = handle to errorbar y OR errorbar x if error y not specified
%      H(3) = handle to errorbar x if error y specified
% LineSpec=[''LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r'']
LineSpec='r.';
hE = ploterr(Y(xip),X1,{model.DrGtMin(xip) model.DrGtMax(xip)},[],LineSpec,'hhxy',0);
set(hE(2),'LineWidth',5);
% hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%dGrt0 errorbar on top and inside dGrt
LineSpec='b.';
hE = ploterr(Y0(xip),X1,{model.DrGt0Min(xip) model.DrGt0Max(xip)},[],LineSpec,'hhxy',0);
set(hE(2),'LineWidth',5);
% hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');

%mean dGrt0
plot(Y0(xip),X1,'.','LineStyle','none','Color',[0.3412 0.7961 0.1922]);
%zero line
%cumulative probability that reaction is really forward, assuming a
%normal distribution about the mean dGt0
[AX,H1,H2]=plotyy(zeros(1,length(X1)),X1,P(xip),X1);
set(gca,'XTick',[]);
set(AX(1),'XTickMode','manual','YTick',floor(minX/100)*200:50:maxX);%,'TickDirMode','manual','TickDir','out');
set(AX(2),'YTick',1:length(xip));
set(H1,'LineStyle','none')
set(H2,'LineStyle','-','LineWidth',2,'Color','k')
plot(zeros(1,length(X1)),X1,'w','LineWidth',2,'LineStyle','--');
%axis limits
set(AX(1),'FontSize',10,'YColor','k')
set(AX(1),'YTick',[],'YTickLabelMode','manual');
set(AX(1),'XTickMode','manual','XTick',floor(minX/100)*200:50:maxX);
set(get(AX(1),'Xlabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
set(get(AX(1),'Xlabel'),'FontSize',10)
set(AX(1),'Position',[0.05 0.1 0.4 0.8]);
% axis(AX(1),[minX maxX 0 length(X1)+0.5 ]);
axis(AX(1),[-350 350 0 length(X1)+0.5 ]);

set(AX(2),'FontSize',10,'YColor','k')
set(AX(2),'YTick',1:length(X1),'YTickLabelMode','manual');
set(AX(2),'XAxisLocation','top');
set(get(AX(2),'Xlabel'),'String','P(\Delta_{r}G^{\primem}<0)');
set(get(AX(2),'Xlabel'),'FontSize',10);
axis(AX(2),[0 1 0.5 length(X1)+0.5 ]);
set(AX(2),'Position',[0.05 0.1 0.4 0.8]);

% ylabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0) (blue)');
title({'Qualitatively forward transport reactions, but quantitatively reversible'},'FontSize',16)

for n=1:length(xip)
%     textLen(n)=length(model.rxn(xip(n)).officialName)+length(model.rxn(xip(n)).equation);
%     textLen(n)=length(model.rxn(xip(n)).equation);
   textLen(n)=length(model.rxnNames{xip(n)});
end
textLenMax=max(textLen);
for n=1:length(xip)
    equation=printRxnFormula(model,model.rxns{xip(n)},0);
    YTickLabel{n}=[model.rxnNames{xip(n)} blanks(textLenMax-textLen(n)+5) equation{1}];
    %get rid of the cytoplasmic compartment abbreviations
    YTickLabel{n} = strrep(YTickLabel{n}, '[c]', '');
    YTickLabel{n} = strrep(YTickLabel{n}, '_', '\_');
end
set(AX(2),'YTickLabel',YTickLabel,'FontSize',10);
saveas(figure1 ,'forwardTransportQuantReverseBoolSetToForward','fig');

author: opencobra

repo: cobratoolbox


text:
function disp(obj)
% display MoxUnitFunctionHandleTestCase object
%
% disp(obj)
%
% Inputs:
%   obj             MoxUnitFunctionHandleTestCase object
%
% NNO 2015

    disp(str(obj));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function compareInterations

fCase{1}  = 'LatticeExperimentInputCantileverSymmetricMeshSuperEllipsePDE';
folder{1} = '/media/alex/My Passport/LatticeResults/CantileverSymmetricMeshSuperEllipsePDEGradientEpsilonH';

fCase{2} = 'LatticeExperimentInputCantileverSymmetricMeshSuperEllipsePDE';
folder{2} = '/media/alex/My Passport/LatticeResults/CantileverSymmetricMeshSuperEllipsePDEEpsilonEhGradient2000';

fCase{3} = 'LatticeExperimentInputCantileverSymmetricMeshSuperEllipsePDE';
folder{3} = '/media/alex/My Passport/LatticeResults/CantileverSymmetricMeshSuperEllipsePDE10000iteration';

f = figure();
hold on
it = 1;
for iPlot = 1:numel(fCase)
    fName = fullfile(folder{iPlot},fCase{iPlot});
    [xV,yV,cV] = getCostFunctions(fName);
    for iline = 1:size(yV,2)
        p{it} = plot(xV,yV(:,iline),'Color',cV(:,iline));
        it = it + 1;
    end
end
pr = plotPrinter(f,p);
pr.print(fullfile('/home/alex/Dropbox/GregMeeting30Octubre',['Iterations']));
end

function [xV,yV,cV] = getCostFunctions(fName)
fNameMon = fullfile([fName,'CostVsLinf.fig']);
h = openfig(fNameMon);
handles = findobj(h,'Type','line');
xV = get(handles(1),'Xdata');
for i = 1:numel(handles)
    yV(:,i) = get(handles(i),'Ydata');
    cV(:,i) = get(handles(i),'Color') ;
end
close(h)
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function success =  runAIforDicom(inputDicomPath,outputDicomPath,...
    sessionPath,algorithm,cmdFlag,savePlanc,varargin)
% function success = runSegForDicom(inputDicomPath,outputDicomPath,...
%   sessionPath,algorithm,varargin)
%
% This function serves as a wrapper for different types of segmentations.
%---------------------------------------------------------------------------------------
% INPUT:
% inputDicomPath  - path to input DICOM directory which needs to be segmented.
% outputDicomPath - path to write DICOM RTSTRUCT for resulting segmentation.
% sessionPath     - path to write temporary segmentation metadata.
% algorithm       - string which specifies segmentation algorithm
% cmdFlag         - "condaEnv" or "singContainer"
% savePlanC       - flag to write planC to file ("yes" or "no")
% --- Optional---
% varargin{1} - Path to segmentation container.
% varargin{2} - Scan no. (replaces input scan identifier)
% varargin{3} - Output assoc. scan no. (replaces output scan identifier)
% varargin{4} - Flag (true/false) to skip export of structure masks (default:true)
%              Set to false if model requires segmentation masks as input
%%---------------------------------------------------------------------------------------
% Following directories are created within the session directory:
% --- ctCERR: contains CERR file/s of input DICOM.
% --- segmentedOrigCERR: CERR file with resulting segmentation fused with
% original CERR file.
% --- segResultCERR: CERR file with segmentation. Note that CERR file can
% be cropped based on initial segmentation.
%
% EXAMPLE: to run BABS segmentation
% inputDicomPath = '';
% outputDicomPath = '';
% sessionPath = '';
% algorithm = 'BABS';
% cmdFlag = 'singcontainer';
%savePlanc = 'yes'
% success = runSegForDicom(inputDicomPath,outputDicomPath,sessionPath,...
%           algorithm,babsPath);
% ------------------------------------------------------------------------------------
% APA, 12/14/2018
% RKP, 9/11/19 Updates for compatibility with training pipeline
% AI, 2/7/2020 Added separate DICOM export functions for BABS and DL algorithms
% AI, 3/5/2020 Updates to handle multiple algorithms
% AI, 8/12/22  Call runSegForPlanC

if nargin <=7
    scanNumV = [];
else
    scanNumV = varargin{2};
end
if nargin <=8
    assocScanNumV = [];
else
    assocScanNumV = varargin{3};
end
if nargin <=9
    skipMaskExport = true;
else
    if ischar(varargin{4})
        skipMaskExport = logical(eval(varargin{4}));
    else
        skipMaskExport = varargin{4};
    end
end

%% Create session directory to write segmentation metadata
if inputDicomPath(end) == filesep
    [~,folderNam] = fileparts(inputDicomPath(1:end-1));
else
    [~,folderNam] = fileparts(inputDicomPath);
end
dateTimeV = clock;
randStr = sprintf('%6.3f',rand*1000);
sessionDir = ['session',folderNam,num2str(dateTimeV(4)), num2str(dateTimeV(5)),...
    num2str(dateTimeV(6)), randStr];
fullSessionPath = fullfile(sessionPath,sessionDir);
while exist(fullSessionPath,'dir')
    randStr = sprintf('%6.3f',rand*1000);
    sessionDir = [sessionDir, randStr];
    fullSessionPath = fullfile(sessionPath,sessionDir);
end

% Create sub-directories
%-For  CERR files
mkdir(fullSessionPath)
cerrPath = fullfile(fullSessionPath,'dataCERR');
mkdir(cerrPath)
outputCERRPath = fullfile(fullSessionPath,'outputOrigCERR');
mkdir(outputCERRPath)
AIresultCERRPath = fullfile(fullSessionPath,'AIResultCERR');
mkdir(AIresultCERRPath)
%-For structname-to-label map
AIoutputPath = fullfile(fullSessionPath,'AIoutput');
mkdir(AIoutputPath);

%% Import DICOM to CERR
tic
recursiveFlag = true;
importDICOM(inputDicomPath,cerrPath,recursiveFlag);
toc

% Get container path
containerPath = varargin{1};

%% Run inference
newSessionFlag = false;
savePlancFlag = 0;
if strcmpi(savePlanc,'yes')
    savePlancFlag = 1;
end
if ~any(strcmpi(algorithm,'BABS'))

    % Get segmentations
    [~,origScanNumV,allLabelNamesC,dcmExportOptS] = runAIforPlanC(cerrPath,...
        fullSessionPath,algorithm,cmdFlag,newSessionFlag,[],[],...
        containerPath,scanNumV,assocScanNumV,skipMaskExport);

    % Export segmentations to DICOM RTSTRUCT files
    fprintf('\nExporting to DICOM format...');
    tic
    batchExportAISegToDICOM(cerrPath,origScanNumV,allLabelNamesC,outputCERRPath,...
        outputDicomPath,dcmExportOptS,savePlancFlag)
    toc

else

    babsPath = varargin{1};
    success = babsSegmentation(cerrPath,fullSessionPath,babsPath,AIresultCERRPath);

    % Export the RTSTRUCT file
    savePlancFlag = 0;
    if strcmpi(savePlanc,'yes')
        savePlancFlag = 1;
    end
    exportCERRtoDICOM_forBABS(origCerrPath,AIresultCERRPath,outputCERRPath,...
        outputDicomPath,dcmExportOptS,savePlancFlag)

end

% Remove session directory
rmdir(fullSessionPath, 's')

success = 0;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% op_CSIZeroBaseline
%
% Takes mean of ppm baseline and substracts spectra by this baseline. This makes
% the ppm baseline zero.
%
% INPUT
% MRSIStruct            = MRSI structure used in FID-A
% ppmBaselineRange      = ppm range of the baseline
%
% OUTPUT
% MRSIStruct            = MRSI ouput

function MRSIStruct = op_CSIzeroBaseline(MRSIStruct, ppmBaselineRange)
    [MRSIStruct, prevPermute, prevShape] = reshapeDimensions(MRSIStruct, {'t', 'y', 'x'});
    data = getData(MRSIStruct);
    ppm = getPPM(MRSIStruct);
    ppmBaselineIndex = ppm > ppmBaselineRange(1) & ppm < ppmBaselineRange(2);
    for e = 1:getSizeFromDimensions(MRSIStruct, {'extras'})
        for x = 1:getSizeFromDimensions(MRSIStruct, {'x'})
            for y = 1:getSizeFromDimensions(MRSIStruct, {'y'})
                voxelBaselineData = mean(data(ppmBaselineIndex, y, x, e));
                realBaseline = real(voxelBaselineData);
                imagBaseline = imag(voxelBaselineData);
                data(:, y, x, e) = data(:, y, x, e) - realBaseline - 1i * imagBaseline;
            end
        end
    end
    MRSIStruct = setData(MRSIStruct, data);
    MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevShape);

end

author: CIC-methods

repo: FID-A

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


text:
function [dpixc, dveccr, dpixc_ind, blinkmat, N] = generatedata2(sizevec, psf, separ, maxphotvec, offset, Nt)
% generate 2 points with different shift

nx = sizevec(2)-sizevec(1);
ny = sizevec(4)-sizevec(3);
N = 2;

if size(maxphotvec,2)>size(maxphotvec,1); maxphotvec = maxphotvec'; end
maxphotmat = repmat(maxphotvec, 1,Nt);

blinkmat_equal = rand(N, Nt);
blinkmat = blinkmat_equal .* maxphotmat; %different intensities...
center = round([nx ny]/2);
centerim = pixelize(center, 1, sizevec, nx, ny, [],0);
dpixc_ind = newimar(2);
dpixc_ind{1} = clip(dip_image(conv2(centerim,psf,'same')), 0, Inf);
dpixc_ind{2} = clip(shift(dpixc_ind{1}, separ), 0, Inf);

dpixc_nonoise = array2im(dpixc_ind'*blinkmat + offset);

dpixc_dip = noise(dpixc_nonoise,'poisson');
dpixc = double(dpixc_dip);
dveccr = double(squeeze(reshape(dpixc, nx*ny, 1, Nt))); % vectors of resized images

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function testFliplrMatrix
%testFliplrMatrix Unit test for fliplr with matrix input

in = magic(3);
assertEqual(fliplr(in), in(:, [3 2 1]));

author: ricedsp


text:
function [url, topdir] = UFlocation
%UFLOCATION URL and top-level directory of the UF Sparse Matrix Collection
%
% Example:
%   [url, topdir] = UFlocation
%
% See also UFget.

% Copyright 2006-2007, Timothy A. Davis

params = UFget_defaults ;
t = find (params.dir == filesep) ;
topdir = regexprep (params.dir (1:t(end-1)), '[\/\\]', filesep) ;

t = find (params.url == '/') ;
url = params.url (1:t(end)) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = limiter(a,b,limtype)
% Limiter function as defined in Sect. 4.8 and as used in Sect. 10.8
if  b==0,  y=0; else
  if limtype == 2          	% van Albada
    y = (a^2 + a*b)/(a^2 + b^2);
  elseif limtype == 1		% Minmod
    y = max(0, min(a/b,1));
  else				% No limiting, no MUSCL
    y = 0;
  end
end

repo: Aero-matlab


text:
function varargout = chebguiEdit(varargin)
%CHEBGUIEDIT   CHEBGUI edittor.
%   A CHEBGUIEDIT figure gets created when a user right-clicks the input fields
%   of the CHEBGUI figure. It is not intended for use in any other context.

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

% Suppress irritating MLINT warnings:
%#ok<*INUSL,*DEFNU,*INUSD>

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @chebguiEdit_OpeningFcn, ...
                   'gui_OutputFcn',  @chebguiEdit_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);

if ( nargin && ischar(varargin{1}) )
    gui_State.gui_Callback = str2func(varargin{1});
end

if ( nargout )
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

% --- Executes just before chebguiEdit is made visible.
function chebguiEdit_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for chebguiEdit

handles.output = hObject;
mainGuiInput = find(strcmp(varargin, 'chebguiWindow'));
if ( isempty(mainGuiInput) ...
     || (length(varargin) <= mainGuiInput) ...
     || ~ishandle(varargin{mainGuiInput+1}) )
    disp('-----------------------------------------------------');
    disp('Improper input arguments. ')
    disp('chebguiEdit should only be called from chebgui.')
    disp('-----------------------------------------------------');
end
% TODO: Why not throw an error here?

% Remember the handle, and adjust our position
handles.chebguiwindow = varargin{mainGuiInput+1};

% Obtain handles using GUIDATA with the caller's handle
mainHandles = guidata(handles.chebguiwindow);
% Set the edit text to the String of the main GUI's button
handles.outputTarget = varargin{3};
set(handles.edit1, 'String', ...
    get(mainHandles.(varargin{3}), 'String'));

% Get the default font size.
set(handles.edit1, 'FontSize', get(mainHandles.tempedit, 'FontSize'));

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes chebguiEdit wait for user response (see UIRESUME)
uiwait(hObject);

end

% --- Outputs from this function are returned to the command line.
function varargout = chebguiEdit_OutputFcn(hObject, eventdata, handles)
% Get default command line output from handles structure
varargout{1} =[];
end

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
if ( ispc && isequal(get(hObject, 'BackgroundColor'), ...
        get(0, 'defaultUicontrolBackgroundColor')) )
    set(hObject, 'BackgroundColor', 'white');
end
end

function edit1_Callback(hObject, eventdata, handles)
end

% --- Executes on button press in buttonOK.
function buttonOK_Callback(hObject, eventdata, handles)
% Obtain handles using GUIDATA with the caller's handle
mainHandles = guidata(handles.chebguiwindow);
% Set the edit text to the String of the main GUI's button
set(mainHandles.(handles.outputTarget), 'String', ...
    get(handles.edit1, 'String'));
% Store the used fontsize
updateEditFontSize(hObject, eventdata, handles)
% Resume and close
uiresume(handles.figure1);
delete(handles.figure1)
end

% --- Executes on button press in buttonCancel.
function buttonCancel_Callback(hObject, eventdata, handles)
% Store the used fontsize
updateEditFontSize(hObject, eventdata, handles);
% Resume and close
uiresume(handles.figure1);
delete(handles.figure1)
end

% --- Executes on button press in buttonClear.
function buttonClear_Callback(hObject, eventdata, handles)
set(handles.edit1, 'String', '');
end

function figure_CloseRequestFcn(hObject, eventdata, handles)
% Store the used fontsize
updateEditFontSize(hObject, eventdata, handles);
% Resume and close
uiresume(hObject);
end

function fontplusbutton_Callback(hObject, eventdata, handles)
% Increase font size.
fs = get(handles.edit1, 'FontSize') + 1;
set(handles.edit1, 'FontSize', fs);
end

function fontmbutton_Callback(hObject, eventdata, handles)
% Decrease font size.
fs = get(handles.edit1, 'FontSize') - 1;
set(handles.edit1, 'FontSize', fs);
end

function edit1_ButtonDownFcn(hObject, eventdata, handles)
% Is the same as buttonCancel
buttonCancel_Callback(hObject, eventdata, handles)
end

function updateEditFontSize(hObject, eventdata, handles)
mainHandles = guidata(handles.chebguiwindow);
set(mainHandles.tempedit, 'FontSize', get(handles.edit1, 'FontSize'));
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function bst_set( varargin )
% BST_SET: Set a Brainstorm structure.
%
% DESCRIPTION:  This function is used to abstract the way that these structures are stored.
%
% USAGE:
% ====== DIRECTORIES ==================================================================
%    - bst_set('BrainstormHomeDir', BrainstormHomeDir)
%    - bst_set('BrainstormTmpDir',  BrainstormTmpDir)
%    - bst_set('BrainstormDbDir',   BrainstormDbDir)
%    - bst_set('LastUsedDirs',      sDirectories)
%    - bst_set('BrainSuiteDir',     BrainSuiteDir)
%    - bst_set('PythonExe',         PythonExe)
%    - bst_set('PluginCustomPath',  PluginCustomPath)
%
% ====== PROTOCOLS ====================================================================
%    - bst_set('iProtocol',         iProtocol)
%    - bst_set('ProtocolInfo',      sProtocolInfo)
%    - bst_set('ProtocolSubjects',  ProtocolSubjects)
%    - bst_set('isProtocolLoaded',  isProtocolLoaded)
%    - bst_set('isProtocolModified',isProtocolModified)
%    - bst_set('ProtocolStudies',   ProtocolStudies)
%    - bst_set('Study',   iStudy,   sStudy)    : Set a study in current protocol
%    - bst_set('Subject', iSubject, sSubject)  : Set a subject in current protocol
%
% ====== GUI =================================================================
%    - bst_set('Layout',    sLayout)
%    - bst_set('Layout',    PropName, PropValue)
%    - bst_set('Clipboard', Nodes, isCut)  : Copy operation from the tree
%
% ====== CONFIGURATION =================================================================
%    - bst_set('Version',      Version)
%    - bst_set('ByteOrder',    value)        : 'b' for big endian, 'l' for little endian
%    - bst_set('AutoUpdates',  isAutoUpdates)
%    - bst_set('ExpertMode',   isExpertMode)
%    - bst_set('DisplayGFP',   isDisplayGFP)
%    - bst_set('DownsampleTimeSeries',  isDownsampleTimeSeries)
%    - bst_set('GraphicsSmoothing',     isGraphicsSmoothing)
%    - bst_set('ForceMatCompression',   isForceCompression)
%    - bst_set('IgnoreMemoryWarnings',  isIgnoreMemoryWarnings)
%    - bst_set('SystemCopy',            isSystemCopy)
%    - bst_set('DisableOpenGL',         isDisableOpenGL)
%    - bst_set('InterfaceScaling',      InterfaceScaling)
%    - bst_set('TSDisplayMode',         TSDisplayMode)    : {'butterfly','column'}
%    - bst_set('ElectrodeConfig',       ElectrodeConfig, Modality)
%    - bst_set('ElecInterpDist',        ElecInterpDist, Modality)
%    - bst_set('DefaultFormats'         defaultFormats)
%    - bst_set('BFSProperties',         [scalpCond,skullCond,brainCond,scalpThick,skullThick])
%    - bst_set('ImportEegRawOptions',   ImportEegRawOptions)
%    - bst_set('BugReportOptions',      BugReportOptions)
%    - bst_set('DefaultSurfaceDisplay', displayStruct)
%    - bst_set('MagneticExtrapOptions', extrapStruct)
%    - bst_set('TimefreqOptions_morlet',  Options)
%    - bst_set('TimefreqOptions_fft',     Options)
%    - bst_set('TimefreqOptions_psd',     Options)
%    - bst_set('TimefreqOptions_hilbert', Options)
%    - bst_set('TimefreqOptions_plv',     Options)
%    - bst_set('OpenMEEGOptions',         Options)
%    - bst_set('DuneuroOptions',         Options)
%    - bst_set('GridOptions_headmodel',   Options)
%    - bst_set('GridOptions_dipfit',      Options)
%    - bst_set('UniformizeTimeSeriesScales', isUniform)
%    - bst_set('FlipYAxis',             isFlipY)
%    - bst_set('AutoScaleY',            isAutoScaleY)
%    - bst_set('FixedScaleY',           Modality,  Value)
%    - bst_set('XScale',                XScale)
%    - bst_set('YScale',                YScale)
%    - bst_set('ShowXGrid',             isShowXGrid)
%    - bst_set('ShowYGrid',             isShowYGrid)
%    - bst_set('ShowZeroLines',         isShowZeroLines)
%    - bst_set('ShowEventsMode',        ShowEventsMode)
%    - bst_set('Resolution',            [resX,resY])
%    - bst_set('UseSigProcToolbox',     UseSigProcToolbox)
%    - bst_set('RawViewerOptions',      RawViewerOptions)
%    - bst_set('TopoLayoutOptions',     TopoLayoutOptions)
%    - bst_set('StatThreshOptions',     StatThreshOptions)
%    - bst_set('ContactSheetOptions',   ContactSheetOptions)
%    - bst_set('ProcessOptions',        ProcessOptions)
%    - bst_set('MriOptions',            MriOptions)
%    - bst_set('CustomColormaps',       CustomColormaps)
%    - bst_set('DigitizeOptions',       DigitizeOptions)
%    - bst_set('ReadOnly',              ReadOnly)
%    - bst_set('LastPsdDisplayFunction', LastPsdDisplayFunction)
%    - bst_set('PlotlyCredentials',     Username, ApiKey, Domain)
%    - bst_set('KlustersExecutable',    ExecutablePath)
%    - bst_set('ExportBidsOptions'),    ExportBidsOptions)
%
% SEE ALSO bst_get

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2021
%          Martin Cousineau, 2017

global GlobalData;

%% ==== PARSE INPUTS ====
if ((nargin >= 1) && ischar(varargin{1}))
    contextName  = varargin{1};
    if (nargin >= 2)
        contextValue = varargin{2};
    else
        contextValue = [];
    end
else
    error('Usage : bst_set(contextName, contextValue)');
end

% Get required context structure
switch contextName
%% ==== BRAINSTORM CONFIGURATION ====
    case 'Version'
        GlobalData.Program.Version = contextValue;
    case 'BrainstormHomeDir'
        GlobalData.Program.BrainstormHomeDir = contextValue;
    case 'BrainstormDbDir'
        GlobalData.DataBase.BrainstormDbDir = contextValue;
    case 'BrainstormTmpDir'
        GlobalData.Preferences.BrainstormTmpDir = contextValue;

%% ==== PROTOCOL ====
    case 'iProtocol'
        if isnumeric(contextValue)
            GlobalData.DataBase.iProtocol = contextValue;
        else
            error('iProtocol should be a number.');
        end
    case {'ProtocolSubjects', 'ProtocolStudies'}
        for structField = fieldnames(contextValue)'
            GlobalData.DataBase.(contextName)(GlobalData.DataBase.iProtocol).(structField{1}) = contextValue.(structField{1});
        end
        GlobalData.DataBase.isProtocolModified(GlobalData.DataBase.iProtocol) = 1;
    case 'ProtocolInfo'
        for structField = fieldnames(contextValue)'
            GlobalData.DataBase.(contextName)(GlobalData.DataBase.iProtocol).(structField{1}) = contextValue.(structField{1});
        end
    case 'isProtocolLoaded'
        GlobalData.DataBase.isProtocolLoaded(GlobalData.DataBase.iProtocol) = contextValue;
    case 'isProtocolModified'
        GlobalData.DataBase.isProtocolModified(GlobalData.DataBase.iProtocol) = contextValue;

%% ==== SUBJECT ====
    case 'Subject'
        % Get subjects list
        ProtocolSubjects = bst_get('ProtocolSubjects');
        iSubject = varargin{2};
        sSubject = varargin{3};
        % If default subject
        if (iSubject == 0)
            ProtocolSubjects.DefaultSubject = sSubject;
        else
            ProtocolSubjects.Subject(iSubject) = sSubject;
        end
        % Update DataBase
        bst_set('ProtocolSubjects', ProtocolSubjects);

%% ==== STUDY ====
    case 'Study'
        % Get studies list
        ProtocolStudies = bst_get('ProtocolStudies');
        iStudies = varargin{2};
        sStudies = varargin{3};
        iAnalysisStudy = -2;
        iDefaultStudy  = -3;
        for i = 1:length(iStudies)
            % Normal study
            if (iStudies(i) > 0)
                ProtocolStudies.Study(iStudies(i)) = sStudies(i);
            % Inter-subject analysis study
            elseif (iStudies(i) == iAnalysisStudy)
                ProtocolStudies.AnalysisStudy = sStudies(i);
            % Default study
            elseif (iStudies(i) == iDefaultStudy)
                ProtocolStudies.DefaultStudy = sStudies(i);
            end
        end
        % Update DataBase
        bst_set('ProtocolStudies', ProtocolStudies);

%% ==== GUI ====
    % USAGE: bst_set('Layout', sLayout)
    %        bst_set('Layout', PropName, PropValue)
    case 'Layout'
        if (nargin == 2) && isstruct(contextValue)
            GlobalData.Preferences.Layout = contextValue;
            isUpdateScreens = 0;
        elseif (nargin == 3) && ischar(contextValue) && isfield(GlobalData.Preferences, 'Layout') && isfield(GlobalData.Preferences.Layout, contextValue)
            GlobalData.Preferences.Layout.(contextValue) = varargin{3};
            isUpdateScreens = strcmpi(contextValue, 'DoubleScreen');
        else
            error('Invalid call to bst_set.');
        end
        % Update screen configuration
        GlobalData.Program.ScreenDef = gui_layout('GetScreenClientArea');
        % Update layout right now
        gui_layout('Update');
        % If the number of screen was changed: update the maximum size of the Brainstorm window
        if isUpdateScreens
            gui_layout('UpdateMaxBstSize');
        end

    % USAGE: bst_set('FixedScaleY', [])
    %        bst_set('FixedScaleY', Modality, Value)
    case 'FixedScaleY'
        if (nargin == 3) && ~isempty(contextValue) && ~isempty(varargin{3})
            GlobalData.Preferences.FixedScaleY.(contextValue) = varargin{3};
        elseif (nargin == 2) && isempty(contextValue)
            GlobalData.Preferences.FixedScaleY = struct();
        end

    case 'ByteOrder'
        switch(contextValue)
            case {'b','ieee-le','n'}
                GlobalData.Preferences.ByteOrder = 'b';
            case {'l','ieee-be'}
                GlobalData.Preferences.ByteOrder = 'l';
            otherwise
                error('Invalid byte order.');
        end

    case 'Clipboard'
        if (length(varargin) >= 3)
            isCut = varargin{3};
        else
            isCut = 0;
        end
        GlobalData.Program.Clipboard.Nodes = contextValue;
        GlobalData.Program.Clipboard.isCut = isCut;

    case 'ElectrodeConfig'
        Modality = varargin{2};
        ElectrodeConf = varargin{3};
        if isequal(Modality, 'ECOG+SEEG')
            Modality = 'ECOG_SEEG';
        elseif ~ismember(Modality, {'EEG','SEEG','ECOG','MEG'})
            error(['Invalid modality: ' Modality]);
        end
        GlobalData.Preferences.(contextName).(Modality) = ElectrodeConf;

    case 'ElecInterpDist'
        Modality = varargin{2};
        ElecInterpDist = varargin{3};
        if isequal(Modality, 'ECOG+SEEG')
            Modality = 'ECOG_SEEG';
        elseif ~ismember(Modality, {'EEG','SEEG','ECOG','MEG'})
            error(['Invalid modality: ' Modality]);
        end
        GlobalData.Preferences.(contextName).(Modality) = ElecInterpDist;

    case {'UniformizeTimeSeriesScales', 'XScale', 'YScale', 'FlipYAxis', 'AutoScaleY', 'ShowXGrid', 'ShowYGrid', 'ShowZeroLines', 'ShowEventsMode', ...
          'Resolution', 'AutoUpdates', 'ExpertMode', 'DisplayGFP', 'ForceMatCompression', 'GraphicsSmoothing', 'DownsampleTimeSeries', ...
          'DisableOpenGL', 'InterfaceScaling', 'TSDisplayMode', 'UseSigProcToolbox', 'LastUsedDirs', 'DefaultFormats', ...
          'BFSProperties', 'ImportDataOptions', 'ImportEegRawOptions', 'RawViewerOptions', 'MontageOptions', 'TopoLayoutOptions', ...
          'StatThreshOptions', 'ContactSheetOptions', 'ProcessOptions', 'BugReportOptions', 'DefaultSurfaceDisplay', ...
          'MagneticExtrapOptions', 'MriOptions', 'ConnectGraphOptions', 'NodelistOptions', 'IgnoreMemoryWarnings', 'SystemCopy', ...
          'TimefreqOptions_morlet', 'TimefreqOptions_hilbert', 'TimefreqOptions_fft', 'TimefreqOptions_psd', 'TimefreqOptions_plv', ...
          'OpenMEEGOptions', 'DuneuroOptions', 'DigitizeOptions', 'CustomColormaps', 'PluginCustomPath', 'BrainSuiteDir', 'PythonExe', ...
          'GridOptions_headmodel', 'GridOptions_dipfit', 'LastPsdDisplayFunction', 'KlustersExecutable', 'ExportBidsOptions'}
        GlobalData.Preferences.(contextName) = contextValue;

    case 'ReadOnly'
        GlobalData.DataBase.isReadOnly = contextValue;

    case 'PlotlyCredentials'
        if length(varargin) ~= 4
            error('Invalid call to bst_set.');
        end
        [username, apiKey, domain] = varargin{2:4};
        % Default domain: plot.ly
        if isempty(domain)
            domain = 'https://plot.ly';
        end
        % Plotly needs a URL with HTTP and no trailing slash.
        if strfind(domain, 'https://')
            domain = strrep(domain, 'https://', 'http://');
        elseif isempty(strfind(domain, 'http://'))
            domain = ['http://', domain];
        end
        if domain(end) == '/'
            domain = domain(1:end-1);
        end
        % Save credentials
        saveplotlycredentials(username, apiKey);
        saveplotlyconfig(domain);

%% ==== ERROR ====
    otherwise
        error('Invalid context : ''%s''', contextName);

end

author: brainstorm-tools


text:
function rhs_spiral_test ( )

%*****************************************************************************80
%
%% RHS_SPIRAL_TEST samples the right hand side at the initial time.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/navier_stokes_2d_exact/rhs_spiral_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2015
%
%  Author:
%
%    John Burkardt
%
  nu = 1.0;
  rho = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'RHS_SPIRAL_TEST\n' );
  fprintf ( 1, '  Spiral Flow:\n' );
  fprintf ( 1, '  Sample the Navier-Stokes right hand sides\n' );
  fprintf ( 1, '  at the initial time T = 0, using the unit square.\n' );
  fprintf ( 1, '  Kinematic viscosity NU = %g\n', nu );
  fprintf ( 1, '  Fluid density RHO = %g\n', rho );

  n = 1000;
  xy_lo = 0.0;
  xy_hi = 1.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, xy_lo, xy_hi, seed );
  [ y, seed ] = r8vec_uniform_ab ( n, xy_lo, xy_hi, seed );
  t = 0.0;

  [ f, g, h ] = rhs_spiral ( nu, rho, n, x, y, t );

  fprintf ( 1, '\n' );
  fprintf ( 1, '           Minimum       Maximum\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F:  %14.6g  %14.6g\n', min ( f ), max ( f ) );
  fprintf ( 1, '  G:  %14.6g  %14.6g\n', min ( g ), max ( g ) );
  fprintf ( 1, '  H:  %14.6g  %14.6g\n', min ( h ), max ( h ) );

  return
end


text:
% This script is a demonstration of how to train a toy acoustic model for
% English using filterbank features for the TIMIT task. Note that the demo
% does not contain the TIMIT corpus. You need to provide the TIMIT corpus
% root directory.
%
% Warning: this demo assumes that your system has at least 8GB of free
% memory for Matlab.
%
% Author: Xiong Xiao, NTU, Singapore
% Date Created: 10 Oct 2013
% Last Modified: 16 Jun 2015
%

clear
para.IO.nStream = 2;        % Two input streams, one is the MFCC features, and the other is the frame-level phone label.
para.NET.sequential = 0;    % Use frame-based training, not sentence based training
para.NET.batchSize = 256;   % minibatch size
para.NET.learning_rate = 3e-2;      % global learning rate
para.NET.momentum = [0.5];          % momentun
para.NET.L2weight = 3e-4;           % L2 regularization weight
para.useGPU = 0;                    % whether to use GPU
para.displayInterval = 100;         % display training progress after N minibathces
para.checkGradient = 0;             % whether to perform gradient checking before training
para.stopImprovement = 0.1;         % when to stop the training, Most of time I use Ctr+C :)
para.reduceLearnRate = 0.5;         % when to start reducing learning rate
para.reduceLearnRateSpeed = 0.7;    % how fast do we decay the learning rate in each iteration
para.maxItr = 50;                   % maximum number of iterations allowed
para.minItr = 20;                   % maximum number of iterations allowed

para.local.V = 39;                  % put task specific configurations into para.local. Here we specify that we use only 39 phones.
para.local.fs = 16000;              % sampling rate, we can use 16000Hz or 8000Hz.
para.local.doCMN = 1;               % whether to perform cepstral mean normalization for each utterance

[Data_tr, Data_cv, para] = LoadData_TIMIT(para, 'train');

% generate the preprocessing of streams. We need to generate dynamic
% featuers for stream 1, i.e. the fbank stream.
para.preprocessing{1}{1}.name = 'delta';
para.preprocessing{1}{end}.delta_order = 2;   % use upto second derivative of MFCCs
para.preprocessing{1}{end}.inputDim = 40;
para.preprocessing{1}{end}.outputDim = 120;

feat_tmp = cell2mat(Data_tr(1).data(1:10:end));
feat_tmp_delta = FeaturePipe(feat_tmp, para.preprocessing{1});  % generate the MFCCs with dynamic features
tmpProcessing = genDNNPreprocessing_splice_norm(feat_tmp_delta, para.IO.context(1));    % generate the splicing setting, and global mean and variance normalization processing
para.preprocessing{1} = [para.preprocessing{1} tmpProcessing];
feat_tmp2 = FeaturePipe(feat_tmp, para.preprocessing{1});
plot(mean(feat_tmp2')); hold on         % verify that we now have normal distributed input features
plot(std(feat_tmp2')); hold off;

para.preprocessing{2} = {};     % we don't need any preprocessing for the second stream, i.e. the label

inputDim = para.preprocessing{1}{end}.outputDim;
hiddenLayerSize = [1024 1024 1024];        % you can generate deeper network by using something like: hiddenLayerSize = [512 512 200 512]. Then it will generate 4 hidden layers.
outputDim = length(unique(cell2mat(Data_tr(2).data)));
cost_function = 'cross_entropy';
layer = genNetworkFeedForward_v2(inputDim, hiddenLayerSize, outputDim, cost_function);
para.cost_func.layer_idx = length(layer);
para.cost_func.layer_weight = [1];

para.output = sprintf('nnet/DNN_phone_TIMIT');
para.output = sprintf('%s.U%d.%d', para.output, length(Data_tr(1).data), inputDim);
for i=1:length(hiddenLayerSize)
    para.output = sprintf('%s-%d', para.output, hiddenLayerSize(i));
end
para.output = sprintf('%s-%d.L2_%s.LR_%s/nnet', para.output, outputDim, FormatFloat4Name(para.NET.L2weight),FormatFloat4Name(para.NET.learning_rate));
LOG = [];

trainGraph_SGD(layer, Data_tr, Data_cv, para, LOG);

author: singaxiong

repo: SignalGraph


text:
%--------------------------------------------------------------------------
% glmnetMex.m: Lasso and elastic-net regularized generalized linear models
%--------------------------------------------------------------------------
%     [a0,ca,ia,nin,rsq,alm,nlp,jerr] = ...
%        glmnetMex(parm,x,y,jd,vp,ne,nx,nlam,flmin,ulam,thr,isd,w,ka)
%     [a0,ca,ia,nin,dev,alm,nlp,jerr] = ...
%        glmnetMex(parm,x,y,jd,vp,ne,nx,nlam,flmin,ulam,thr,isd,nc,maxit,kopt)
%
% Extremely efficient procedures for fitting the entire lasso or
% elastic-net regularization path for linear regression, logistic and
% multinomial regression models. The algorithm uses cyclical coordinate
% descent in a pathwise as described in the paper on the maintainer's
% website.
%
% NOTES: This is a MEX-file wrapper of GLMnet.f for MATLAB. Should be called
% only by glmnet.m. For details about input and output arguments, see
% GLMnet.f.
%
% LICENSE: GPL-2
%
% DATE: 13 Jul 2009
%
% AUTHORS:
%    Algorithm was designed by Jerome Friedman, Trevor Hastie and Rob Tibshirani
%    Fortran code was written by Jerome Friedman
%    R wrapper (from which the MATLAB wrapper was adapted) was written by Trevor Hasite
%    MATLAB wrapper was written and maintained by Hui Jiang, jiangh@stanford.edu
%    Department of Statistics, Stanford University, Stanford, California, USA.
%
% REFERENCES:
%    Friedman, J., Hastie, T. and Tibshirani, R. (2009)
%    Regularization Paths for Generalized Linear Models via Coordinate Descent.
%    Journal of Statistical Software, 33(1), 2010
%
% SEE ALSO:
%    glmnet, glmnetSet, glmnetPrint, glmnetPlot, glmnetPredict and glmnetCoef methods.
%
% EXAMPLES:
%    parm = 1.0;
%    x = [1 1; 2 2; 3 3];
%    y = [1 3 2]';
%    jd = 0;
%    vp = [1 1];
%    ne = 3;
%    nx = 2;
%    nlam = 100;
%    flmin = 0.0001;
%    ulam = 0;
%    thr = 1.0e-4;
%    isd = 0;
%    w = [1 1 1]';
%    ka = 1;
%    [a0,ca,ia,nin,rsq,alm,nlp,jerr] = glmnetMex(parm,x,y,jd,vp,ne,nx,nlam,flmin,ulam,thr,isd,w,ka);
%
% DEVELOPMENT: 13 Jul 2009: Original version of glmnetMex.m written.


text:
function [post,lik,lli] = logistK_eval(beta,x,y)
% [post,lik,lli] = logistK_eval(beta,x,y)
%
% Evaluate logistic regression model.
%
% INPUT
% 	beta 	dxk model coefficients (as returned by logistK)
% 	x 	dxn matrix of n input column vectors
% 	[y] 	kxn vector of class assignments
%
% OUTPUT
% 	post 	kxn fitted class posteriors
% 	lik 	1xn vector of sample likelihoods
%	lli	log likelihood
%
% Let p(i,j) = exp(beta(:,j)'*x(:,i)),
% Class j posterior for observation i is:
%	post(j,i) = p(i,j) / (p(i,1) + ... p(i,k))
% The likelihood of observation i given soft class assignments
% y(:,i) is:
%	lik(i) = prod(post(:,i).^y(:,i))
% The log-likelihood of the model given the labeled samples is:
%	lli = sum(log(lik))
%
% See also logistK.
%
% David Martin <dmartin@eecs.berkeley.edu>
% May 7, 2002

% Copyright (C) 2002 David R. Martin <dmartin@eecs.berkeley.edu>
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
% 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html.

error(nargchk(2,3,nargin));

% check sizes
if size(beta,1) ~= size(x,1),
  error('Inputs beta,x not the same height.');
end
if nargin > 3 & size(y,2) ~= size(x,2),
  error('Inputs x,y not the same length.');
end

% get sizes
[d,k] = size(beta);
[d,n] = size(x);

% class posteriors
post = zeros(k,n);
bx = zeros(k,n);
for j = 1:k,
  bx(j,:) = beta(:,j)'*x;
end
for j = 1:k,
  post(j,:) = 1 ./ sum(exp(bx - repmat(bx(j,:),k,1)),1);
end
clear bx;

% likelihood of each sample
if nargout > 1,
  y = y ./ repmat(sum(y,1),k,1); % L1-normalize class assignments
  lik = prod(post.^y,1);
end

% total log likelihood
if nargout > 2,
  lli = sum(log(lik+eps));
end;

% eof


text:
function zz=v_lpccw2zz(cw)
%V_LPCPZ2ZZ LPC: Power spectrum roots to LPC poles ZZ=(CW)
% pz are the roots of the power spectrum polynomial pp(cos(w))

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_lpccw2zz.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zs=sqrt(cw.^2-1);
zz=cw-sign(real(conj(cw).*zs)).*zs;

author: ImperialCollegeLondon

repo: sap-voicebox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: Regularized Parametric Image Registration,  by spline-moment matrix
%
%   - data                 HNSP, Omega=(0,2)x(0,1), level=5, m=[256,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             SSD
%   - pre-registration     splineTransformation2D, regularized by spline-moment matrix
%   - optimization         Gauss-Newton
% ===============================================================================

clear, close all, help(mfilename);
%%
setup2DHNSPData;  level = 5;  omega = ML{level}.omega; m = ML{level}.m;
p = [8,8];
imgModel('reset','imgModel','splineInter','regularizer','none','theta',0);
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
distance('reset','distance','SSD');
trafo('reset','trafo','splineTransformation2D','omega',omega,'m',m,'p',p);
w0 = trafo('w0');                        % get starting guess and stopping

% update the regularizer for parametric image registration
hd = prod((omega(2:2:end)-omega(1:2:end))./m);
Mi = @(i) 1e2*toeplitz([96,-54,0,6,zeros(1,p(i)-4)]);
Qi = @(i) toeplitz([120.8,59.55,6,0.05,zeros(1,p(i)-4)])/7;
M1 = kron(speye(p(2)),sparse(Mi(1)));
M2 = kron(sparse(Mi(2)),speye(p(1)));
M  = hd*sparse(kron(speye(2),...
  kron(Qi(2),Mi(1))+2*kron(Mi(2),Mi(1))+kron(Mi(2),Qi(1))));

%% set up elastic regularization matrix
mu     = 2;
lambda = 1;
B = getElasticMatrixNodal(omega,m,mu,lambda);
A = B'*B;
[y,dy] = splineTransformation2D(w0,getNodalGrid(omega,m),'m',m+1);
M = [dy'*A*dy];

%%
% M = speye(length(w0));
alpha = 10;
beta  = 0;
wRef  = w0;

% initialize objective function
xc = getCellCenteredGrid(omega,m);
Rc = imgModel(R,omega,xc);
fctn  = @(wc) PIRobjFctn(T,Rc,omega,m,beta,alpha*M,wRef,xc,wc);
fctn([]);

% initialize plots
FAIRplots('reset','mode','PIR-regularized','omega',omega,'m',m,'fig',1,'plots',1);
FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

xc = getCellCenteredGrid(omega,m);
Rc = imgModel(R,omega,xc);

% ----- call Gaus-Newton ------------------------------------
GNoptn = {'maxIter',100,'tolY',1e-4,'tolJ',1e-4,'tolG',1,'solver','backslash'};
[wOpt,his] = GaussNewton(fctn,w0,GNoptn{:},'Plots',@FAIRplots);

% plot iteration history
his.str{1} = sprintf('iteration history PIR: distance=%s, y=%s',distance,trafo);
[ph,th] = plotIterationHistory(his,'J',1:4);
%==============================================================================


text:
function loss=cost231(d,nf,kw)
%%
% d is the distance from the AP
% nf is the number of the floors encoutered
% kw is the vector cointaining the number of the type walls met
%  - the first element is for Light wall: plasterboard,particle board or thin
%    (<10 cm), light concrete wall.
%  - the second element is for Heavy wall: thick (>10cm), concrete or brick

% Loss for wall types in dB
Lw=0;
Lwa=[3.4 6.9];
for i=1:length(kw)
    Lw=Lw+kw(i)*Lwa(i);
end

% Loss for the floors in dB
Lfi=18.3;
Lf=20*log10(d)+nf^((nf+2)/(nf+1))*Lfi;

% Free space Path Loss in dB
lambda=0.125;
do=1;
Lo=20*log10((4*pi*do/lambda)^2);

loss=Lw+Lf+Lo;


text:
% hessMultdalgl - function that computes H*x for DAL with grouped
%                 L1 regularization
%
% Copyright(c) 2009 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt

function yy = hessMultdalgl(xx, A, eta, Hinfo)

yy = Hinfo.hloss*xx;
if ~isempty(Hinfo.precomp)
    % general case
    for p=Hinfo.precomp
        % project & reshape xx's for all spans
        xk = reshape(p.AJ'*xx,[],numel(p.jj));
        % dot product between vn's and xk's, scale by ff's and multiply by vn's (also add xk*(1-ff))
        tmp = bsxfun(@times,p.vn,sum(p.vn.*xk).*p.ff) + bsxfun(@times,xk,p.omff);
        % block-diagonalize tmp
        bd = sparse(p.idxu,p.idxv,tmp,numel(xk),numel(p.jj));
        % map through AJ again, sum and multiply by eta(1), add results to yy
        yy = yy + sum(p.AJ*bd,2)*eta(1);
    end
else
    % very sparse case
    blks =Hinfo.blks;
    hloss=Hinfo.hloss;
    I    =Hinfo.I;
    vv   =Hinfo.vv;
    nm   =Hinfo.nm;
    lambda=Hinfo.lambda;
    for kk=1:length(I)
        jj=I(kk);
        J=Hinfo.blkival{jj};
        vn=vv(J)/nm(jj);
        ff=lambda/nm(jj);
        AJ=A.slice(J);
        xk=AJ'*xx;
        yy = yy + eta(1)*(AJ*((1-ff)*xk + ff*(vn'*xk)*vn));
    end
end

B=Hinfo.B;
if ~isempty(B)
  yy = yy + eta(2)*(B*(B'*xx));
end


text:
classdef SOP_F21 < PROBLEM
% <single> <real> <expensive/none>
% Shekel's family

%------------------------------- Reference --------------------------------
% X. Yao, Y. Liu, and G. Lin, Evolutionary programming made faster, IEEE
% Transactions on Evolutionary Computation, 1999, 3(2): 82-102.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 1;
            obj.D = 4;
            obj.lower    = zeros(1,obj.D);
            obj.upper    = zeros(1,obj.D) + 10;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            a = [4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7];
            c = [0.1;0.2;0.2;0.4;0.4];
            PopObj = zeros(size(PopDec,1),1);
            for i = 1 : size(PopDec,1)
                PopObj(i) = -sum(1./(sum((repmat(PopDec(i,:),5,1)-a).^2,2)+c));
            end
        end
        %% Generate the minimum objective value
        function R = GetOptimum(obj,N)
            R = -10.16;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function y_one_hot = one_hot(y)

y_one_hot = zeros(size(y, 1), 30);
for i = 1:30 % TODO: 30 3D models in the graph
    rows = y == i;
    y_one_hot( rows, i ) = 1;
end

author: jhonykaesemodel


text:
% Isoparametric Formulation Implementation
% clear memory
clear all;

figure(5)

% Displacements
subplot(1,2,1); hold on;
%plot(xx,displacements_exact,'-r');

% Stress
subplot(1,2,2); hold on;
%plot(xx,stress_exact,'-r');

%% Physical Parameters
% Initial Area
A0 = 2; % m^2

% Modulus of elasticity
E = 8; % N/m^2

% Length of bar
Ltotal = 4; % m

%% Build Elements
% number of elements,
%numberElements = 4;
test = [1 2 4 8];
for numberElements = test

% number of nodes,
numberNodes = 2*numberElements+1;

% node coordinates,
x = linspace(2,6,numberNodes)';

% Generation of coordinates and connectivities
L = zeros(numberElements,1); elementNodes = zeros(numberElements,3);
L(1) = Ltotal/numberElements; elementNodes(1,:) = [1,2,3];
for i = 2:numberElements
    L(i) = L(1); elementNodes(i,:) = elementNodes(i-1,:) + 2;
end
nodeCoordinates = x;

% element center coordinate
x_mid = 0.5*(x(elementNodes(:,3))-x(elementNodes(:,1)))+x(elementNodes(:,1));

% Evaluate area the center of each element.
% A: area of cross section
A = A0*x_mid;

% for structure:
   % displacements: displacement vector
   % force : force vector
   % stiffness: stiffness matrix
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
% computation of the system stiffness matrix and force vector
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:) ;
  detJacobian=L(e)/2;
  invJacobian=1/detJacobian;
  ngp = 3;
  [w,xi]=gauss1d(ngp);
  xc=0.5*(nodeCoordinates(elementDof(1))+nodeCoordinates(elementDof(3)));
  for ip=1:ngp;
      [shape,naturalDerivatives]=shapeFunctionL3(xi(ip));
      B=naturalDerivatives*invJacobian;
      stiffness(elementDof,elementDof)=...
      stiffness(elementDof,elementDof)+ B'*B*w(ip)*detJacobian*E*A(e);
      force(elementDof) = force(elementDof)+...
          8*(xc+detJacobian*xi(ip))*shape'*detJacobian*w(ip);
  end
end

%% BCs and solution
% prescribed dofs
prescribedDof = find(nodeCoordinates == 2); % fixed at node x = 2 m
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness, ...
    numberNodes,prescribedDof,force)

%% Stress
elementCoords=zeros(numberElements*ngp,1);
stress=zeros(numberElements*ngp,1);
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:);
  detJacobian=L(e)/2;
  invJacobian=1/detJacobian;
  [w,xi]=gauss1d(ngp);
  xc=0.5*(nodeCoordinates(elementDof(1))+nodeCoordinates(elementDof(3)));
  for ip=1:ngp;
      [shape,naturalDerivatives]=shapeFunctionL3(xi(ip));
      B=naturalDerivatives*invJacobian;
      elementCoords(ip+(e-1)*3,1)=xc+xi(ip)*detJacobian;
      stress(ip+(e-1)*3,1)=E*B*displacements(elementDof,1);
  end
end

%% plot figures

switch numberElements
    case 1
        line='r*--';
    case 2
        line='g*--';
    case 4
        line='k*--';
    case 8
        line='m*--';
end

% Displacements
subplot(1,2,1); plot(nodeCoordinates,displacements,line)

% Stress
subplot(1,2,2); axis([2,6,0,26]); plot(elementCoords,stress,line)

end %end for
subplot(1,2,1); legend('Exact','NEL=1','NEL=2','NEL=4','NEL=8',4); hold off
subplot(1,2,2); legend('Exact','NEL=1','NEL=2','NEL=4','NEL=8',1); hold off


text:
classdef KalmanFilterX < FilterX
% KalmanFilterX class
%
% Summary of KalmanFilterX:
% This is a class implementation of a standard Kalman Filter.
%
% KalmanFilterX Properties: (**)
%   + StatePrior - A structure used to store the state prior
%   + StatePrediction - A structure used to store the state prediction
%   + MeasurementPrediction - A structure used to store the measurement prediction
%   + StatePosterior - A structure used to store posterior information
%   + MeasurementList - A (yDim x 1) matrix used to store the received measurement
%   + ControlInput - A (uDim x 1) matrix used to store the last received control input
%   + KalmanGain - A (xDim x yDim) matrix representing the last computed Kalman Gain
%   + Model - An object handle to StateSpaceModelX object
%       + Transition (*)  = Object handle to TransitionModelX SubClass
%       + Measurement (*)  = Object handle to MeasurementModelX SubClass
%       + Control (*)  = Object handle to ControlModelX SubClass
%
%   (*)  Signifies properties necessary to instantiate a class object
%   (**) xDim, yDim and uDim denote the dimentionality of the state, measurement
%        and control vectors respectively.
%
% KalmanFilterX Methods:
%   + KalmanFilterX  - Constructor method
%   + predict        - Performs KF prediction step
%   + update         - Performs KF update step
%
% (+) denotes puplic properties/methods
%
% See also TransitionModelX, MeasurementModelX and ControlModelX template classes

    properties
        StatePrior
        StatePrediction
        MeasurementPrediction
        StatePosterior
        KalmanGain
        ControlInput
    end

    properties (Dependent)
        MeasurementLikelihoods
    end

    properties (Access=protected)
        MeasurementLikelihoods_ = [];
    end

    methods (Access = protected)
        function initialise_(this, config)
            initialise_@FilterX(this,config);
            if (isfield(config,'StatePrior'))
                this.StatePrior = config.StatePrior;
                this.StatePosterior = this.StatePrior;
            end
        end
        % MeasurementList
        function measurementList = setMeasurementList(this, newMeasurementList)
            measurementList = newMeasurementList;
            this.MeasurementLikelihoods_ = [];
        end
        function StatePrior = setStatePrior(this,newStatePrior)
            if(isa(newStatePrior,'GaussianStateX'))
                StatePrior = newStatePrior;
            else
                StatePrior = GaussianStateX(newStatePrior);
            end
        end
        function StatePrediction = setStatePrediction(this,newStatePrediction)
            if(isa(newStatePrediction,'GaussianStateX'))
                StatePrediction = newStatePrediction;
            else
                StatePrediction = GaussianStateX(newStatePrediction.Mean, newStatePrediction.Covar);
            end
        end
        function MeasurementPrediction = setMeasurementPrediction(this,newMeasurementPrediction)
            if(isa(newMeasurementPrediction,'GaussianStateX'))
                MeasurementPrediction = newMeasurementPrediction;
            else
                MeasurementPrediction = GaussianStateX(newMeasurementPrediction.Mean, newMeasurementPrediction.Covar);
            end
        end
        function MeasurementLikelihoods = getMeasurementLikelihoods(this)
            if(isempty(this.MeasurementLikelihoods_))
                this.MeasurementLikelihoods_ =  mvnpdf(this.MeasurementList.Vectors',this.MeasurementPrediction.Mean',this.MeasurementPrediction.Covar)';
            end
            MeasurementLikelihoods = this.MeasurementLikelihoods_;
        end
        function StatePosterior = setStatePosterior(this,newStatePosterior)
            if(isa(newStatePosterior,'GaussianStateX'))
                StatePosterior = newStatePosterior;
            else
                StatePosterior = GaussianStateX(newStatePosterior.Mean, newStatePosterior.Covar);
            end
        end
    end

    methods
        function this = KalmanFilterX(varargin)
        % KalmanFilterX Constructor method
        %
        % Parameters
        % ----------
        % Model: StateSpaceModelX
        %   An object handle to StateSpaceModelX object.
        % StatePrior: struct, optional
        %   A StateX subclass object describing the state prior. If StatePrior
        %   is not a GaussianStateX instance, then it will be converted in
        %   one using the extracted mean and covariance.
        %
        % Usage
        % -----
        % * kf = KalmanFilterX(___,Name,Value) instantiates an object handle,
        %   configured with the options specified by one or more Name,Value
        %   pair arguments.
        %
        %  See also predict, update, smooth.

            % Call SuperClass method
            %this@FilterX(varargin{:});

            if(nargin==0)
                return;
            end

            % First check to see if a structure was received
            if(nargin==1)
                if(isstruct(varargin{1}))
                    config = varargin{1};
                    this.initialise_(config);
                end
                return;
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            config = parser.Unmatched;
            this.initialise_(config);
        end

        function initialise(this,varargin)
        % initialise Initialise the KalmanFilter with a certain set of
        %   parameters.
        %
        % Parameters
        % ----------
        % Model: StateSpaceModelX
        %   An object handle to StateSpaceModelX object.
        % StatePrior: StateX, optional
        %   A StateX subclass object describing the state prior. If StatePrior
        %   is not a GaussianStateX instance, then it will be converted in
        %   one using the extracted mean and covariance.
        %
        % Usage
        % -----
        % * initialise(kf,___,Name,Value) initialises the KalmanFilterX
        %   object kf with the options specified by one or more Name,Value
        %   pair arguments.
        %
        %  See also predict, update, smooth.

            if(nargin==0)
                error("Not enough input arguments.");
            end

            initialise@FilterX(this);

            % First check to see if a structure was received
            if(nargin==2)
                if(isstruct(varargin{1}))
                    config = varargin{1};
                    this.initialise_(config);
                end
                return;
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.parse(varargin{:});
            config = parser.Unmatched;
            this.initialise_(config);
        end

        function [statePrediction, measurementPrediction] = predict(this, varargin)
        % Predict Perform Kalman Filter prediction step
        %
        % Parameters
        % ----------
        % prior: GaussianStateX, optional
        %   The prior state estimate.
        % timestamp: datetime, optional
        %   A timestamp indicating the time at which prediction is
        %   performed.
        %
        % Returns
        % -------
        % GaussianStateX
        %   The generated state prediction
        % GaussianStateX, optional
        %   The generated measurement prediction
        %
        %  See also update, smooth.

            % Predict state and measurement
            statePrediction = this.predictState(varargin{:});
            if nargin>1 && isa(varargin{1},'StateX')
               % Replace a potential prior with the generated prediction
               % before forwarding the arguments to the measurement
               % prediction. Failure to do so will result in errors!!!
               varargin{1} = statePrediction;
            end
            measurementPrediction = this.predictMeasurement(varargin{:});
        end

        function statePrediction = predictState(this,varargin)
        % predictState Perform Kalman Filter state prediction step
        %
        % Usage
        % -----
        % * predictState(this) calculates the predicted system state and covariance.
        %
        % See also update, smooth.

            timestamp = [];
            timestamp_old = [];
            for i = 1:min([2,nargin-1])
                if isa(varargin{i},'StateX')
                    this.StatePosterior = varargin{i};
                    timestamp_old = this.StatePosterior.Timestamp;
                elseif isdatetime(varargin{i})
                    timestamp = varargin{i};
                end
            end

            if isempty(timestamp)
                dt = this.Model.Transition.TimestepDuration;
                timestamp = this.StatePosterior.Timestamp;
            else
                dt = timestamp - timestamp_old;
            end

            % Extract model parameters
            F = this.Model.Transition.matrix(dt);
            Q = this.Model.Transition.covar(dt);
            if(~isempty(this.Model.Control))
                B   = this.Model.Control.feval();
                Qu  = this.Model.Control.covar();
            else
                this.ControlInput   = 0;
                B   = 0;
                Qu  = 0;
            end

            % Perform state prediction
            [statePredictionMean, statePredictionCovar] = ...
                this.predictState_(this.StatePosterior.Mean, this.StatePosterior.Covar, F, Q, this.ControlInput, B, Qu);

            statePrediction = GaussianStateX(statePredictionMean, statePredictionCovar, timestamp);
            this.StatePrediction = statePrediction;
        end

        function measurementPrediction = predictMeasurement(this, varargin)
        % PREDICTOBS Perform Kalman Filter measurement prediction step
        %
        % Usage
        % -----
        % * predict(this) calculates the predicted measurement,
        %   as well as the associated uncertainty covariances.
        %
        % More details
        % ------------
        % * KalmanFilterX uses the Model class property, which should be an
        %   instance of the TrackingX.Models.StateSpaceModel class, in order
        %   to extract information regarding the underlying state-space model.
        % * State prediction is performed using the Model.Transition property,
        %   which must be a subclass of TrackingX.Abstract.TransitionModel and
        %   provide the following interface functions:
        %   - Model.Transition.feval(): Returns the model transition matrix
        %   - Model.Transition.covariance(): Returns the process noise covariance
        % * Measurement prediction and innovation covariance calculation is
        %   performed usinf the Model.Measurement class property, which should be
        %   a subclass of TrackingX.Abstract.TransitionModel and provide the
        %   following interface functions:
        %   - Model.Measurement.heval(): Returns the model measurement matrix
        %   - Model.Measurement.covariance(): Returns the measurement noise covariance
        %
        % See also update, smooth.

            if nargin>1
                this.StatePrediction = varargin{1};
            end

            % Extract model parameters
            H = this.Model.Measurement.feval();
            R = this.Model.Measurement.covar();

            % Perform prediction
            [measurementPredictionMean, measurementPredictionCovar, this.KalmanGain] = ...
                this.predictMeasurement_(this.StatePrediction.Mean, this.StatePrediction.Covar, H, R);

            measurementPrediction = GaussianStateX(measurementPredictionMean,...
                                                   measurementPredictionCovar,...
                                                   this.StatePrediction.Timestamp);
            this.MeasurementPrediction = measurementPrediction;
        end

        function posterior = update(this, varargin)
        % UPDATE Perform Kalman Filter update step
        %
        % Usage
        % -----
        % * update(this) calculates the corrected sytem state and the
        %   associated uncertainty covariance.
        %
        % See also KalmanFilterX, predict, iterate, smooth.

            if nargin>1
                if isa(varargin{1},'MeasurementX')
                    this.MeasurementList = MeasurementListX(varargin{1});
                elseif isa(varargin{1}, 'StateX')
                    this.StatePrediction = varargin{1};
                    this.MeasurementList =  MeasurementListX(varargin{2});
                end
            end
            if(this.MeasurementList.NumMeasurements)
                timestamp = this.MeasurementList.Timestamp;
            else
                timestamp = this.StatePrediction.Timestamp;
            end

            if(isempty(this.MeasurementPrediction.Mean) || isempty(this.MeasurementPrediction.Covar))
                [measurementPredictionMean, measurementPredictionCovar, this.KalmanGain] = ...
                    this.predictMeasurement_(this.StatePrediction.Mean, this.StatePrediction.Covar, H, R);
                measurementPrediction = GaussianStateX(measurementPredictionMean,...
                                                       measurementPredictionCovar,...
                                                       this.StatePrediction.Timestamp);
                this.MeasurementPrediction = measurementPrediction;
            end

            % Perform single measurement update
            [posteriorMean, posteriorCovar] = this.update_(this.StatePrediction.Mean,this.StatePrediction.Covar,...
                                                           this.MeasurementList.Vectors,this.MeasurementPrediction.Mean,...
                                                           this.MeasurementPrediction.Covar,this.KalmanGain);
            posteriorCovar = (posteriorCovar+posteriorCovar')/2;

            posterior = GaussianStateX(posteriorMean, posteriorCovar, timestamp);
            this.StatePosterior = posterior;
        end

        function posterior = updatePDA(this, assocWeights, varargin)
        % UPDATEPDA Performs KF update step, for multiple measurements
        %           Update is performed according to the generic (J)PDAF equations [1]
        %
        % Usage
        % -----
        %  * updatePDA(assocWeights) Performs KF-PDA update step for multiple
        %    measurements based on the provided (1-by-Nm+1) association weights
        %    matrix assocWeights.
        %
        %   [1] Y. Bar-Shalom, F. Daum and J. Huang, "The probabilistic data association filter," in IEEE Control Models, vol. 29, no. 6, pp. 82-100, Dec. 2009.
        %
        %   See also KalmanFilterX, Predict, Iterate, Smooth, resample.

            if(this.MeasurementList.NumMeasurements)
                timestamp = this.MeasurementList.Timestamp;
            else
                timestamp = this.StatePrediction.Timestamp;
            end

            [posteriorMean, posteriorCovar] = ...
                this.updatePDA_(this.StatePrediction.Mean,this.StatePrediction.Covar,this.MeasurementList.Vectors,...
                                assocWeights,this.MeasurementPrediction.Mean,this.MeasurementPrediction.Covar,this.KalmanGain);

            posterior = GaussianStateX(posteriorMean,posteriorCovar,timestamp);
            this.StatePosterior = posterior;
        end

        function measurementLikelihoods = get.MeasurementLikelihoods(this)
            measurementLikelihoods = getMeasurementLikelihoods(this);
        end

        function statePrior = get.StatePrior(this)
            statePrior = this.StatePrior;
        end

        function set.StatePrior(this, newStatePrior)
            this.StatePrior = setStatePrior(this, newStatePrior);
        end

        function statePrediction = get.StatePrediction(this)
            statePrediction = this.StatePrediction;
        end

        function set.StatePrediction(this, newStatePrediction)
            this.StatePrediction = setStatePrediction(this, newStatePrediction);
        end

        function measurementPrediction = get.MeasurementPrediction(this)
            measurementPrediction = this.MeasurementPrediction;
        end

        function set.MeasurementPrediction(this, newMeasurementPrediction)
            this.MeasurementPrediction = setMeasurementPrediction(this, newMeasurementPrediction);
        end

        function statePosterior = get.StatePosterior(this)
            statePosterior = this.StatePosterior;
        end

        function set.StatePosterior(this, newStatePosterior)
            this.StatePosterior = setStatePosterior(this, newStatePosterior);
        end

    end

    methods (Static)

        function [xPred, PPred, yPred, S, K] = predict_(x,P,F,Q,H,R,u,B,O)
        % PREDICT_ Perform the discrete-time KF state and measurement
        % prediction steps, under the assumption of additive process noise.
        %
        % Parameters
        % ----------
        % x: column vector
        %   The (xDim x 1) state estimate at the previous time-step.
        % P: matrix
        %   The (xDim x xDim) state covariance matrix at the previous
        %   time-step.
        % F: matrix
        %   An (xDim x xDim) state transition matrix.
        % Q: matrix
        %   The (xDim x xDim) process noise covariance matrix.
        % H: matrix
        %   A (xDim x yDim) measurement matrix.
        % R: matrix
        %   The (yDim x yDim) measurement noise covariance matrix.
        % u: column vector, optional
        %   A optional (xDim x 1) control input.
        %   If omitted, no control input is used.
        % B: matrix, optional
        %   An optional (xDim x xDim) control gain matrix.
        %   If omitted, B is assumed to be 1.
        % O: matrix, optional
        %   An optional (xDim x xDim) control noise covariance
        %   matrix. If omitted, Q is assumed to be 0.
        %
        % Returns
        % -------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix.
        % yPred: column vector
        %   The (yDim x 1) predicted measurement estimate.
        % Pxy: matrix
        %   The (xDim x yDim) cross-covariance matrix.
        % S: matrix
        %   The (yDim x yDim) innovation covariance matrix.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            switch(nargin)
                case(6)
                    u = 0;
                    B = 0;
                    O = 0;
                case(7)
                    B = 1;
                    O = 0;
                case(8)
                    O = 0;
            end

           [xPred, PPred] = KalmanFilterX.predictState_(x,P,F,Q,u,B,O);
           [yPred, S, K] = KalmanFilterX.predictMeasurement_(xPred,PPred,H,R);
        end

        function [xPred, PPred] = predictState_(x,P,F,Q,u,B,Qu)
        % PREDICTSTATE_ Perform the discrete-time KF state prediction
        % step, under the assumption of additive process noise.
        %
        % Parameters
        % ----------
        % x: column vector
        %   The (xDim x 1) state estimate at the previous time-step.
        % P: matrix
        %   The (xDim x xDim) state covariance matrix at the previous
        %   time-step.
        % F: matrix
        %   An (xDim x xDim) state transition matrix.
        % Q: matrix
        %   The (xDim x xDim) process noise covariance matrix.
        % u: column vector, optional
        %   An optional (xDim x 1) control input.
        %   If omitted, no control input is used.
        % B: matrix, optional
        %   An optional (xDim x xDim) control gain matrix.
        %   If omitted, B is assumed to be 1.
        % O: matrix, optional
        %   An optional (xDim x xDim) control noise covariance
        %   matrix. If omitted, Q is assumed to be 0.
        %
        % Returns
        % -------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            switch(nargin)
                case(4)
                    u  = 0;
                    B  = 0;
                    Qu = 0;
                case(5)
                    B  = 1;
                    Qu = 0;
                case(6)
                    Qu = 0;
            end

            % Compute predicted state mean and covariance
            xPred = F*x + B*u;
            PPred =F*P*F' + Q + B*Qu*B';
        end

        function [yPred, S, K] = predictMeasurement_(xPred,PPred,H,R)
        % PREDICTOBS_ Perform the discrete-time KF observation prediction
        % step, under the assumption of additive process noise.
        %
        % Parameters
        % ----------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate at the current
        %   time-step.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix at
        %   the current time-step.
        % H: matrix
        %   An (xDim x yDim) measurement matrix.
        % R: matrix
        %   The (yDim x yDim) measurement noise covariance matrix.
        %
        % Returns
        % -------
        % yPred: column vector
        %   The (yDim x 1) predicted measurement estimate.
        % Pxy: matrix
        %   The (xDim x yDim) cross-covariance matrix.
        % S: matrix
        %   The (yDim x yDim) innovation covariance matrix.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            % Compute predicted measurement mean and covariance
            yPred   = H*xPred;
            Pxy     = PPred*H';
            S       = H*PPred*H' + R;
            K       = Pxy/(S);
        end

        function [x,P] = update_(xPred,PPred,y,yPred,S,K)
        % UPDATE_ Perform the discrete-time KF update step, under the
        % assumption of additive process noisem for a single measurement.
        %
        % Parameters
        % ----------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix.
        % y: column vector
        %   The (yDim x 1) measurement vector.
        % yPred: column vector
        %   The (yDim x 1) predicted measurement estimate.
        % S: matrix
        %   The (yDim x yDim) innovation covariance matrix.
        % K: matrix
        %   The (xDim x yDim) Kalman gain matrix at the current
        %   time-step.
        %
        % Returns
        % -------
        % x: column vector
        %   The (xDim x 1) state estimate at the current time-step.
        % P: matrix
        %   The (xDim x xDim) state covariance matrix at the current
        %   time-step.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            % Compute the filtered estimates
            x = xPred + K * (y - yPred);
            P = PPred - K*S*K';
        end

        function [x,P] = updatePDA_(xPred,PPred,Y,W,yPred,S,K)
        % UPDATEPDA_ Perform the discrete-time Probabilistic Data
        % Association (PDA) KF update step, under the assumption of additive process
        % noise, for multiple measurements (as a Gaussian Mixture)
        %
        % Parameters
        % ----------
        % xPred: column vector
        %   The (xDim x 1) predicted state estimate.
        % PPred: matrix
        %   The (xDim x xDim) predicted state covariance matrix.
        % Y: matrix
        %   The (yDim x nY) measurement vector.
        % W: row vector
        %   The (1 x nY+1) measurement association/mixture weights
        %   vector. (dummy measurement assumed at index 1)
        % yPred: column vector
        %   The (yDim x 1) predicted measurement estimate.
        % S: matrix
        %   The (yDim x yDim) innovation covariance matrix.
        % K: matrix
        %   The (xDim x yDim) Kalman gain matrix at the current
        %   time-step.
        %
        % Returns
        % -------
        % x: column vector
        %   The (xDim x 1) state estimate at the current time-step.
        % P: matrix
        %   The (xDim x xDim) state covariance matrix at the current
        %   time-step.
        %
        %October 2017 Lyudmil Vladimirov, University of Liverpool.

            % Get size of observation vector
            nY = size(Y,2);
            if(nY==0)
                x = xPred;
                P = PPred;
                return;
            end

            innov_err = Y - yPred;
            xupd = [xPred, xPred + K*innov_err];
            Pplus = PPred - K*S*K';

            try
                x = xupd*W';
            catch
                asd=2;
            end
            v_x = x - xupd;
            P = W(1)*(PPred + v_x(:,1)*v_x(:,1)');
            for j = 2:nY+1
                P = P + W(j)*(Pplus + v_x(:,j)*v_x(:,j)');
            end

%             % Compute innovation mean and (cross) covariance
%             innov_err       = Y - yPred(:,ones(1,nY));
%             tot_innov_err   = innov_err*W(2:end)';
%             Pc              = PPred - K*S*K';
%             Pgag            = K*((innov_err.*W(ones(yDim,1),2:end))*innov_err' - tot_innov_err*tot_innov_err')*K';
%
%             % Compute filtered estimates
%             x    = xPred + K*tot_innov_err;
%             P    = W(1)*PPred + (1-W(1))*Pc + Pgag;
             P    = (P+P')/2;
        end

        function config = getInitConfig()

        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
J = 5;  % J: number of stages

% get filters
[Faf, Fsf] = FSfarras;
[af, sf] = dualfilt1;

x = zeros(1,256);  % zero signal

% Compute dual-tree complex DWT of zero signal
w = dualtree(x, J, Faf, af);
% Set a single (real) coefficient to 1
w{5}{1}(4) = 1;
% Compute the inverse transform
y1 = idualtree(w, J, Fsf, sf);

% Compute dual-tree complex DWT of zero signal
w = dualtree(x, J, Faf, af);
% Set a single (imaginary) coefficient to 1
w{5}{2}(4) = 1;
% Compute the inverse transform
y2 = idualtree(w, J, Fsf, sf);

% Display real and imaginary parts and magnitude
n = [1:256]/256;
plot(n,y1,n,y2,n,sqrt(y1.^2+y2.^2))
title('COMPLEX 1D WAVELET')
xlabel('t');
ylabel('\psi(t)');

author: thomaskuestner