text:
%................................................................

% MATLAB codes for Finite Element Analysis
% problem16vibrations.m
% Timoshenko beam in free vibrations
% AJM Ferreira 2008
% modified by Manuel Diaz, 2013.07.25

% clear memory
clear all

% E; modulus of elasticity
% G; shear modulus
% I: second moments of area
% L: length of beam
% thickness: thickness of beam
<<<<<<< HEAD
E=30e6; poisson = 0.30;L  = 1;thickness=0.01;
=======
E=10e7; poisson = 0.30;L = 3;thickness=0.001;
>>>>>>> 47f4967135b009f9128c106a6744c97a2ffd18ce
I=thickness^3/12;
EI=E*I;
kapa=5/6;
rho=1;
A=1*thickness;
%

P = 0; % uniform pressure
%P = -1; % uniform pressure
% constitutive matrix
G=E/2/(1+poisson);
C=[   EI   0; 0    kapa*A*G];

% mesh
<<<<<<< HEAD
numberElements     = 50;
=======
numberElements  = 40;
>>>>>>> 47f4967135b009f9128c106a6744c97a2ffd18ce
nodeCoordinates=linspace(0,L,numberElements+1);
xx=nodeCoordinates';x=xx';
for i=1:size(nodeCoordinates,2)-1
    elementNodes(i,1)=i;
    elementNodes(i,2)=i+1
end
% generation of coordinates and connectivities
numberNodes=size(xx,1);

% GDof: global number of degrees of freedom
GDof=2*numberNodes;

% computation of the system stiffness, force, mass
[stiffness,force,mass]=...
    formStiffnessMassTimoshenkoBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,C,P,rho,I,thickness);

% boundary conditions (simply-supported at both bords)
%fixedNodeW =[1 ; numberNodes];
%fixedNodeTX=[];
% boundary conditions (clamped at both bords)
%fixedNodeW =[1 ; numberNodes];
%fixedNodeTX=[1 ; numberNodes];
% boundary conditions (cantilever)
%fixedNodeW =[1];
%fixedNodeTX=[1];
<<<<<<< HEAD
% Free-free
%fixedNodeW =[];
%fixedNodeTX=[];
=======
% Free-Free conditions
fixedNodeW =[];
fixedNodeTX=[];
>>>>>>> 47f4967135b009f9128c106a6744c97a2ffd18ce
prescribedDof=[fixedNodeW; fixedNodeTX+numberNodes];

% solution (optional*)
displacements=solution(GDof,prescribedDof,stiffness,force);

% output displacements/reactions (optional*)
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)

% free vibration problem
activeDof=setdiff([1:GDof]',[prescribedDof]);
modeNumber=4;

[V,D]=eig(stiffness(activeDof,activeDof),...
    mass(activeDof,activeDof));
D = diag(sqrt(D)*L*L*sqrt(rho*A/E/I));
[D,ii] = sort(D);

V1=zeros(GDof,1);
V1(activeDof,1:modeNumber)=V(:,1:modeNumber);
!
% drawing eigenmodes
drawEigenmodes1D(modeNumber,numberNodes,V1,xx,x)

repo: Aero-matlab


text:
% In this small MATLAB script, we compute the least squares solution to a
% regression problem subject to L1 regularization, which rewards "sparse"
% models that have regression coefficients of zero. See, for instance,
% the work in the "Lasso" by the statistician Robert Tibshirani.
%
% Copyright (C) 2008 Peter Carbonetto. All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% Author: Peter Carbonetto
%         Dept. of Computer Science
%         University of British Columbia
%         September 18, 2008

% Experiment parameters.
lambda = 1;                      % Level of L1 regularization.
n      = 100;                    % Number of training examples.
e      = 1;                      % Std. dev. in noise of outputs.
beta   = [ 0 0 2 -4 0 0 -1 3 ]'; % "True" regression coefficients.

% Set the random number generator seed.
seed = 7;
rand('state',seed);
randn('state',seed);

% CREATE DATA SET.
% Generate the input vectors from the standard normal, and generate the
% binary responses from the regression with some additional noise, and then
% transform the results using the logistic function. The variable "beta" is
% the set of true regression coefficients.
m     = length(beta);      % Number of features.
A     = randn(n,m);        % The n x m matrix of examples.
noise = e * randn(n,1);    % Noise in outputs.
y     = A * beta + noise;  % The binary outputs.

% COMPUTE SOLUTION WITH IPOPT.
% Compute the L1-regularized maximum likelihood estimator.
w = lasso(A,y,lambda);
fprintf('Solution:\n');
disp(w);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [f, df, d2f] = eval_legacy_cost(om, x, name, idx)
%EVAL_LEGACY_COST  Evaluates individual or full set of legacy user costs.
%   F = OM.EVAL_LEGACY_COST(X ...)
%   [F, DF] = OM.EVAL_LEGACY_COST(X ...)
%   [F, DF, D2F] = OM.EVAL_LEGACY_COST(X ...)
%   [F, DF, D2F] = OM.EVAL_LEGACY_COST(X, NAME)
%   [F, DF, D2F] = OM.EVAL_LEGACY_COST(X, NAME, IDX_LIST)
%   Evaluates an individual named set or the full set of legacy user
%   costs and their derivatives for a given value of the optimization vector
%   X, based on costs added by ADD_LEGACY_COST.
%
%   Example:
%       [f, df, d2f] = om.eval_legacy_cost(x)
%       [f, df, d2f] = om.eval_legacy_cost(x, name)
%       [f, df, d2f] = om.eval_legacy_cost(x, name, idx)
%
%   See also OPT_MODEL, ADD_LEGACY_COST, PARAMS_LEGACY_COST.

%   MATPOWER
%   Copyright (c) 2008-2017, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if om.cost.N
    done = 0;

    %% collect cost parameters
    if nargin < 3                       %% full set
        [cp, vs] = om.params_legacy_cost();
    elseif nargin < 4 || isempty(idx)   %% name, no idx provided
        dims = size(om.cost.idx.i1.(name));
        if prod(dims) == 1              %% simple named set
            [cp, vs] = om.params_legacy_cost(name);
        elseif nargout == 1             %% indexing required, recurse
            f = 0;          %% initialize cumulative cost
            idx = num2cell(ones(size(dims))); %% initialize idx
            while ~done     %% call eval_legacy_cost() recursively
                f = f + om.eval_legacy_cost(x, name, idx);

                %% increment idx
                D = length(dims);
                idx{D} = idx{D} + 1;    %% increment last dimension
                for d = D:-1:2          %% increment next dimension, if necessary
                    if idx{d} > dims(d)
                        idx{d} = 1;
                        idx{d-1} = idx{d-1} + 1;
                    end
                end
                if idx{1} > dims(1)     %% check if done
                    done = 1;
                end
            end
        else
            error('@opt_model/eval_legacy_cost: legacy cost set ''%s'' requires an IDX_LIST arg when requesting DF output', name)
        end
    else                                %% indexed named set
        [cp, vs] = om.params_legacy_cost(name, idx);
    end

    if ~done
        %% assemble appropriately-sized x vector
        xx = om.varsets_x(x, vs, 'vector');

        %% compute function & derivatives
        if nargout == 1
            f = opf_legacy_user_cost_fcn(xx, cp);
        elseif nargout == 2
            [f, df] = opf_legacy_user_cost_fcn(xx, cp);
        else    %% nargout == 3
            [f, df, d2f] = opf_legacy_user_cost_fcn(xx, cp);
        end
    end
else
    f = 0;
    if nargout > 1
        df = [];
        if nargout > 2
            d2f = [];
        end
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef spherefun < separableApprox
%SPHEREFUN class for representing functions on the unit sphere.
%
%   Class for approximating functions defined on the unit sphere. The
%   functions should be smooth.
%
% SPHEREFUN(F) constructs a SPHEREFUN object representing the function F on
% the unit sphere. F can have the following form:
%    1. A function handle in (x,y,z), e.g., @(x,y,z) x.*y.*z + cos(x).
%    2. A function handle in spherical coordinates (lambda,theta), where
%       lambda is the azimuthal variable and satisfies -pi <= lambda <= pi
%       and theta is the polar angle and satisfies 0 <= theta < pi,
%       e.g., @(lambda,theta) cos(cos(lambda).*sin(theta))
%    3. A matrix of numbers.
% If F is a function handle then it should allow for vectorized evaluations.
%
% If F is a matrix, F = (f_ij), the numbers fij are used as function values
% at tensor equally-spaced points in the intrinsic spherical coordinate
% system, i.e., [-pi,pi]x[0,pi].
%
% SPHEREFUN(F, k) returns a rank k approximation to F.
%
% SPHEREFUN(F, [m n]) returns a representation of F using a degree m
% trigonometric approximation of F in the theta (polar) direction and a
% degree n trigonometric approximation in the lambda (azimuthal) direction.
% The result is compressed in low rank form and the rank k is still
% determined adaptively (satisfying k<=min(m,n)+1).
%
% The SPHEREFUN software system is based on:
%
% A. Townsend, H. Wilber, and G. Wright, Computing with functions in
% spherical and polar geometries I: The sphere, SIAM. J. Sci. Comput., 38-4
% (2016), C403-C425.
%
% See also CHEBFUN2, DISKFUN, SPHEREFUNV

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

% TODO: Include documentation of fixed eps construction

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public, Static = false)

        function f = spherefun(varargin)
            % The main spherefun constructor!

            % Return an empty CHEBFUN:
            if ( (nargin == 0) || isempty(varargin{1}) )
                f.domain = [-pi pi 0 pi];
                return
            end

            f = constructor(f, varargin{:});
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        %f = conj(f);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% HIDDEN METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public, Static = false, Hidden = true)
        % Project a spherefun to have exact BMC-I symmetry so that it is
        % a continuous function on the sphere.
        f = projectOntoBMCI(f);

        % Fast evaluation of a spherefun using non-uniform fft.
        vals = fastSphereEval( f, lambda, theta );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% PUBLIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public, Static = false)

        % The main bulk of the SPHEREFUN constructor:
        g = constructor(g, op, dom, varargin);

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public, Static = true)

        % Poisson solver:
        u = poisson(f, const, m, n);

        % Helmholtz solver:
        u = helmholtz(f, K, m, n);

        % Convert matrix of coefficients to a spherefun:
        f = coeffs2spherefun(CFS);

        varargout = coeffs2vals(U, varargin);
        varargout = vals2coeffs(U, varargin);

        % Convert a function in spherical coordinates to one in Cartesian
        % coordinates on the sphere.
        fdf = sphf2cartf(f, lam, th, coord);

        % Degree l Order m spherical harmonic.
        Y = sphharm(l, m, coord);

        % Plot the outline of the landmasses of earth
        h = plotEarth(linespec);

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Private Static methods implemented by SPHEREFUN class.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = private, Static = true)
        % Fast evaluation of a spherefun using non-uniform fft.
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties (Access = public)
        % DOMAIN: default is [-pi,pi] x [0,pi] which corresponds to using
        % colatitude for the elevation angle (second input argument).
        % Doubled-up sphere will have a domain of [-pi,pi] x [-pi,pi].
        idxPlus
        idxMinus
        nonZeroPoles = 0;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Private constant properties
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties (Constant)
        % TODO: Add support for this constant here.
        % alpha = 50;  % Growth factor control.
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function prob = cwp2Fn(paramS,doseBinsC,volHistC)
% prob = cwp2Fn(paramS,doseBinsC,volHistC);
% ----------------------------------------------------------------------
% AI, 1/31/19

%Extract parameters and corresponding weights
[weight,x,parListC] = getParCoeff(paramS,'weight',doseBinsC,volHistC);

%Compute TCP/NTCP
constIdx = strcmpi(parListC,'constant');
gx = sum(weight(~constIdx).*x(~constIdx));
const = weight(constIdx).*x(constIdx);
prob = 1- exp(const*exp(gx));

%%--- Functions to extract model parameters and weights ----
    function [coeff,par,keepParC] = getParCoeff(paramS,fieldName,doseBinsC,volHistC)

        %Extract relevant parameters
        genFieldC = fields(paramS);
        for i = 1:numel(genFieldC)
            if strcmpi(genFieldC{i},'structures')
                structS = paramS.(genFieldC{i});
                structListC = fieldnames(structS);
                for j = 1:length(structListC)
                    strParamS = structS.(structListC{j});
                    strParamListC = fieldnames(strParamS);
                    for k = 1 : numel(strParamListC)
                        if isfield(strParamS.(strParamListC{k}),'cteg') |  isfield(strParamS.(strParamListC{k}),'weight')
                            parName = [structListC{j},strParamListC{k}];
                            keepParS.(parName) = strParamS.(strParamListC{k});
                        end
                    end
                end
            else
                 if isfield(paramS.(genFieldC{i}),'cteg') |  isfield(paramS.(genFieldC{i}),'weight')
                     parName = genFieldC{i};
                     keepParS.(parName) = paramS.(genFieldC{i});
                 end
            end
        end

       %Compute parameters, extract coefficients
        keepParC = fieldnames(keepParS);
        numStr = 0;
        par = zeros(1,length(keepParC));
        coeff = zeros(1,length(keepParC));
        for n = 1:length(keepParC)
            coeff(n) = keepParS.(keepParC{n}).(fieldName);
            if isnumeric(keepParS.(keepParC{n}).val)
               par(n) = keepParS.(keepParC{n}).val;
            else
                if ~iscell(doseBinsC)  %For single-structure models
                    doseBinsV = doseBinsC;
                    volHistV = volHistC;
                else
                numStr = numStr+1;
                doseBinsV = doseBinsC{numStr};
                volHistV = volHistC{numStr};
                end
                if ~isfield(keepParS.(keepParC{n}),'params')
                    par(n) = eval([keepParS.(keepParC{n}).val,...
                        '(doseBinsV, volHistV)']);
                else
                    %Copy number fo fractions, abRatio
                    keepParS.(keepParC{n}).params.numFractions.val = paramS.numFractions;
                    if isfield(paramS,'abRatio')
                        keepParS.(keepParC{n}).params.abRatio.val = paramS.abRatio;
                    end
                    par(n) = eval([keepParS.(keepParC{n}).val,...
                        '(doseBinsV, volHistV,keepParS.(keepParC{n}).params)']);
                end
            end
        end
    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%CIRCGRID Generate 2d quadrilateral grid for a circle.
%
%   [ GRID ] = CIRCGRID( NS, NR, R, XP, TH_OFFSET ) Generates a quadrilateral
%   grid for a circular domain with NS+NR cells in the radial direction. NS
%   specifies the cell resolution of the inner square (default 4), and NR the
%   number of cells in the radial direction of the outer layer (default 3).
%   The optional arguments R and XP = [x0;y0] specify the radius and center
%   coordinates of the circle (default R = 1 and XP = [0;0]). Furthermore,
%   TH_OFFSET specifies a rotation of the whole grid in radians.
%
%   Examples:
%
%      1) A 64 cell grid for a circle with radius 1 centered at [ 0, 0 ]:
%
%      grid = circgrid();
%
%      2) A 1024 cell grid for a circle with radius 0.5 centered at [ 1, 1 ]:
%
%      grid = circgrid( 16, 12, 0.5, [1;1], 0 );
%
%   See also CYLGRID, BLOCKGRID, HOLEGRID, LINEGRID, RECTGRID, RINGGRID, SPHEREGRID

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function hd = DirectionCalulate(X1,Y1,X2,Y2)
    hd=rem(atan2d(Y2-Y1,X2-X1)+360, 360);
end

author: kavli-ntnu


text:
function shepard_interp_1d_test ( )

%*****************************************************************************80
%
%% SHEPARD_INTERP_1D_TEST tests the SHEPARD_INTERP_1D library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 August 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../r8lib' )
  addpath ( '../test_interp' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SHEPARD_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SHEPARD_INTERP_1D library.\n' );
  fprintf ( 1, '  This test needs the TEST_INTERP library.\n' );

  prob_num = p00_prob_num ( );

  for prob = 1 : prob_num

    for p = 1 : [ 0.0, 1.0, 2.0, 4.0, 8.0 ];
      shepard_interp_1d_test01 ( prob, p );
    end

  end

  prob = 7;
  for p = 1 : [ 0.0, 1.0, 2.0, 4.0, 8.0 ];
    shepard_interp_1d_test02 ( prob, p );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SHEPARD_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../r8lib' )
  rmpath ( '../test_interp' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [out_image_m,out_ref_points_m] = rotate_image( degree, in_image_m, in_ref_points_m )
%
% rotate_image - rotates an image given inside a matrix by the amount of "degree" counter-clockwise
%                using linear interpolation of the output grid points from the back-rotated input points
%                in this way, the output image will never have a "blank" point
%
% Format:   [out_image_m,out_ref_points_m] = rotate_image( degree, in_image_m, in_ref_points_m )
%
% Input:    degree          - rotation degree in dergees, counter-clockwise
%           in_image_m      - input image, given inside a matrix (gray level image only)
%           in_ref_points_m - points on the image wich their output coordinates will be given
%                             after the rotation. given format of this matrix is:
%                             [ x1,x2,...,xn;y1,y2,...,yn]
%
% Output:   out_image_m      - the output image
%           out_ref_points_m - the position of the input handle points after the rotation.
%                              this element is given in "in_ref_points_m" exists
%                              format of the matrix is the same as of "in_ref_points_m"
%
% NOTE:     By definition of rotation, in order to perserve all the image inside the
%           rotated image space, the output image will be a matrix with a bigger size.
%

%  NO INPUT ARGs - Launch demo and exit
if (nargin == 0)
    rotate_image_demo;
    out_image_m = [];
    return;
end

% check input
if ~exist('in_ref_points_m')
    in_ref_points_m = [];
end

% check for easy cases
switch (mod(degree,360))
case 0,
    out_image_m      = in_image_m;
    out_ref_points_m = in_ref_points_m;
    return;
case 90,
    out_image_m           = in_image_m(:,end:-1:1)';
    out_ref_points_m      = in_ref_points_m(end:-1:1,:);
    out_ref_points_m(2,:) = size(out_image_m,1) - out_ref_points_m(2,:);
    return;
case 180,   % TBD for rotation of the ref_points
    out_image_m           = in_image_m(end:-1:1,end:-1:1);
    out_ref_points_m      = in_ref_points_m;
    out_ref_points_m(2,:) = size(out_image_m,2) - out_ref_points_m(2,:);
    out_ref_points_m(1,:) = size(out_image_m,1) - out_ref_points_m(1,:);
    return;
case 270,
    out_image_m           = in_image_m(end:-1:1,:)';
    out_ref_points_m      = in_ref_points_m(end:-1:1,:);
    out_ref_points_m(1,:) = size(out_image_m,2) - out_ref_points_m(1,:);
    return;
otherwise,  % enter the routine and do some calculations
end

% wrap input image by zeros from all sides
zeros_row    = zeros(1,size(in_image_m,2)+2);
zeros_column = zeros(size(in_image_m,1),1);
in_image_m   = [zeros_row; zeros_column,in_image_m,zeros_column; zeros_row ];

% build the rotation matrix
degree_rad = degree * pi / 180;
R = [ cos(degree_rad), sin(degree_rad); sin(-degree_rad) cos(degree_rad) ];

% input and output size of matrices (output size is found by rotation of 4 corners)
in_size_x       = size(in_image_m,2);
in_size_y       = size(in_image_m,1);
in_mid_x        = (in_size_x-1) / 2;
in_mid_y        = (in_size_y-1) / 2;
in_corners_m    = [ [0,0,in_size_x-1,in_size_x-1] - in_mid_x;
                    [0,in_size_y-1,in_size_y-1,0] - in_mid_y ];
out_corners_m   = R * in_corners_m;

% the grid (integer grid) of the output image and the output image
[out_x_r,out_y_r]   = rotated_grid( out_corners_m );
out_size_x          = max( out_x_r ) - min( out_x_r ) + 1;
out_size_y          = max( out_y_r ) - min( out_y_r ) + 1;
out_image_m         = zeros( ceil( out_size_y ),ceil( out_size_x ) );
out_points_span     = (out_x_r-min(out_x_r))*ceil(out_size_y) + out_y_r - min(out_y_r) + 1;
if ~isempty( in_ref_points_m )
    out_ref_points_m    = (R * [in_ref_points_m(1,:)-in_mid_x;in_ref_points_m(2,:)-in_mid_y]);
    out_ref_points_m    = [out_ref_points_m(1,:)-min( out_x_r )+1;out_ref_points_m(2,:)-min( out_y_r )+1];
else
    out_ref_points_m    = [];
end

% % for debug
% out_image_m(out_points_span) = 1;
% return;
% % end of for debug

% the position of points of the output grid in terms of the input grid
in_cords_dp_m   = inv(R) * [out_x_r;out_y_r];

x_span_left     = floor(in_cords_dp_m(1,:) + in_mid_x + 10*eps );
y_span_down     = floor(in_cords_dp_m(2,:) + in_mid_y + 10*eps );
x_span_right    = x_span_left + 1;
y_span_up       = y_span_down + 1;
dx_r            = in_cords_dp_m(1,:) - floor( in_cords_dp_m(1,:) + 10*eps );
dy_r            = in_cords_dp_m(2,:) - floor( in_cords_dp_m(2,:) + 10*eps );

point_span_0_0  = x_span_left*ceil(in_size_y)  + y_span_down + 1; % position of combined index in output matrix
point_span_1_0  = x_span_left*ceil(in_size_y)  + y_span_up + 1;
point_span_0_1  = x_span_right*ceil(in_size_y) + y_span_down + 1;
point_span_1_1  = x_span_right*ceil(in_size_y) + y_span_up + 1;

out_image_m(out_points_span) = ...
    in_image_m( point_span_0_0 ).*(1-dx_r).*(1-dy_r) + ...
    in_image_m( point_span_1_0 ).*(1-dx_r).*(  dy_r) + ...
    in_image_m( point_span_0_1 ).*(  dx_r).*(1-dy_r) + ...
    in_image_m( point_span_1_1 ).*(  dx_r).*(  dy_r);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%              Inner function implementation                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x_r,y_r] = rotated_grid( rect_points_m )
%
% rotated_grid - creates a grid of points bounded inside a rotated RECTANGLE
%
% Format:   [x_m,y_m] = rotated_grid( rect_points_m )
%
% Input:    rect_points_m   -   a set of (x;y) points which define a rectangle ordered clock-wise
%                               ( format: [x1,x2,x3,x4;y1,y2,y3,y4] )
%
% Output:   x_r,y_r         -   2 row vectors which hold the x and y positions of
%                               the output grid
%
% NOTE:     THE ASSUMPTION IS THAT THE RECTANGLE IS ORDERED CLOCK-WISE !!!
%           AND THAT THE GIVEN CO-ORDINATES ARE A RECTANGLE !
%

% make sure that the first point of the clock-wise-ordered rectange is of the most left point
[temp,idx] = min( rect_points_m(1,:) );
if ( idx > 1 )
    rect_points_m = [ rect_points_m(:,idx:end) , rect_points_m(:,1:idx-1) ];
end

% put into variables so it is easier to access/read the numbers
x1 = rect_points_m(1,1);
x2 = rect_points_m(1,2);
x3 = rect_points_m(1,3);
x4 = rect_points_m(1,4);
y1 = rect_points_m(2,1);
y2 = rect_points_m(2,2);
y3 = rect_points_m(2,3);
y4 = rect_points_m(2,4);

% initialization for grid creation
clipped_top     = floor( y2 );
clipped_bottom  = ceil( y4 );
fraction_bottom = clipped_bottom - y4;
rows            = ( clipped_top - clipped_bottom );
left_crossover  = y1 - y4;
right_crossover = y3 - y4;

% calculate the position of the edges (left and right) along the y axis
m = [0:rows] + fraction_bottom ;
switch (y1)
case y2, x_left = repmat( ceil( x4 ),size(m) );
case y4, x_left = repmat( ceil( x2 ),size(m) );
otherwise
    x_left = ( m >= left_crossover ).*ceil( x2 - (x1-x2)/(y1-y2)*(rows-m+2*fraction_bottom) ) + ...
        ( m < left_crossover ).*ceil( x4 + (x1-x4)/(y1-y4)*m );
end
switch (y3)
case y2,    x_right = repmat( floor( x4 ),size(m) );
case y4,    x_right = repmat( floor( x2 ),size(m) );
otherwise
    x_right = ( m >= right_crossover ).*floor( x2 - (x3-x2)/(y3-y2)*(rows-m+2*fraction_bottom) ) + ...
        ( m < right_crossover ).*floor( x4 + (x3-x4)/(y3-y4)*m );
end

% build the output vectors (initialize)
vec_length = sum(x_right-x_left+1);
x_r = zeros(1,vec_length );
y_r = zeros(1,vec_length );

% build the grid into the output vectors
cursor = 1;
for n = 1:length(m)
    if ( x_right(n) >= x_left(n) )
        span        = cursor:(x_right(n) - x_left(n) + cursor);
        x_r( span ) = x_left(n):x_right(n);
        y_r( span ) = m(n) + y4;
        cursor      = cursor + x_right(n) - x_left(n) + 1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%              Demo implementation of this routine                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rotate_image_demo

% plot the "child" image, and get it's matrix
close all;
h           = imagesc;
in_image_m  = get( h,'cdata' );
set( get( h,'parent' ),'ydir','reverse' );
title( 'original image' );

% create targets on the image
[sy,sx]     = size( in_image_m );
hold on;
in_ref_points_m = [ [0.05 0.05 0.5 0.95 0.75 0.95]*sx; [0.05 0.7 0.95 0.7 0.3 0.05]*sy ];
plot( in_ref_points_m(1,:),in_ref_points_m(2,:),'k','linewidth',2 );
hold off;

% loop over selected angles and plot the roated image with it's targets
for degree = [0 15 25 30 45 60 75 90]
    [out_image_m,out_ref_points_m] = rotate_image( degree, in_image_m, in_ref_points_m );
    figure;
    imagesc( out_image_m );
    title( sprintf( 'Rotated image by %d degrees',degree ) );
    hold on;
    plot( out_ref_points_m(1,:),out_ref_points_m(2,:),'k','linewidth',2 );
    hold off;
end

author: Sable

repo: mcbench-benchmarks


text:
function [vfick,vfick2] = dquat2fick(dquat)

% DQUAT2FICK  transforms a rotation dual quaternion into Fick coordinates
%
%     [VFICK,VFICK2] = DQUAT2FICK(ROTMATRIX) returns the Fick coordinates
%       VFICK and VFICK2 [deg] corresponding to the rotation dual
%       quaternion DQUAT.
%        - DQ is a rotation dual quaternion. It is a 8*N array (column i
%            represents the rotation i). DQ = cos(theta/2)+n*sin(theta/2)
%            has the unitary n vecor expressed in the following reference
%            frame:
%            The positive rotation around x corresponds to a clockwise
%            torsional rotation. The positive rotation around y corresponds
%            to a downward rotation. The positive rotation about z
%            corresponds to a leftward rotation.
%        - VFICK (resp. VFICK2) is a 3*N array (column i represents
%           rotation i). VFICK(:,i) can be written as VFICK = [H V T] where
%           (from the subject's perspective):
%            - H is the horizontal fick angle (positive when to the left)
%            - V is the vertical fick angle (around the inter-aural axis).
%               It is positive when down and is comprise between -90 and
%               90deg.
%            - T is the torsional fick angle (around the line of signt). It
%               is positive when clockwise.
%         - VFICK2 is the second solution VFICK2 = [H2 V2 T2] and V2 is
%             comprised between 90 and 270deg.
%
% See also DQUAT2ROTMATRIX, ROTMATRIX2FICK, FICK2DQUAT

rotMatrix = dquat2rotMatrix(dquat);
[vfick,vfick2] = rotMatrix2fick(rotMatrix);


text:
function [psnr,ssim_val]=compute_psnr(im1,im2,shave_border)
if size(im1, 3) == 3,
    im1 = rgb2ycbcr(im1);
    im1 = im1(:, :, 1);
end

if size(im2, 3) == 3,
    im2 = rgb2ycbcr(im2);
    im2 = im2(:, :, 1);
end

imdff = double(im1) - double(im2);
if shave_border > 0
    imdff = shave(imdff,[shave_border,shave_border]);
end
imdff = imdff(:);
im1 = shave(im1,[shave_border,shave_border]);
im2 = shave(im2,[shave_border,shave_border]);
ssim_val = ssim_index(im1,im2);
rmse = sqrt(mean(imdff.^2));
psnr = 20*log10(255/rmse);

author: huangzehao

repo: caffe-vdsr


text:
function tapas_hgf_binary_plotTraj(r)
% Plots the estimated or generated trajectories for the binary HGF perceptual model
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_hgf_binary_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Optional plotting of responses (true or false)
ploty = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Time axis
if size(r.u,2) > 1
    t = r.u(:,end)';
else
    t = ones(1,size(r.u,1));
end

ts = cumsum(t);
ts = [0, ts];

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = (length(r.p_prc.p)+1)/5;
end

% Upper levels
for j = 1:l-1

    % Subplots
    subplot(l,1,j);

    if plotsd == true
        upperprior = r.p_prc.mu_0(l-j+1) +sqrt(r.p_prc.sa_0(l-j+1));
        lowerprior = r.p_prc.mu_0(l-j+1) -sqrt(r.p_prc.sa_0(l-j+1));
        upper = [upperprior; r.traj.mu(:,l-j+1)+sqrt(r.traj.sa(:,l-j+1))];
        lower = [lowerprior; r.traj.mu(:,l-j+1)-sqrt(r.traj.sa(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mu_0(l-j+1); r.traj.mu(:,l-j+1)], 'b', 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(l-j+1), 'ob', 'LineWidth', 2); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu_', num2str(l-j+1)]);
end

% Input level
subplot(l,1,l);

plot(ts, [tapas_sgm(r.p_prc.mu_0(2), 1); tapas_sgm(r.traj.mu(:,2), 1)], 'r', 'LineWidth', 2);
hold all;
plot(0, tapas_sgm(r.p_prc.mu_0(2), 1), 'or', 'LineWidth', 2); % prior
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0.6 0]); % inputs
plot(ts(2:end), r.traj.wt(:,1), 'k') % implied learning rate
if (ploty == true) && ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    if ~isempty(find(strcmp(fieldnames(r),'c_sim'))) && strcmp(r.c_sim.obs_model,'tapas_beta_obs')
        y = r.y(:,1);
    else
        y = r.y(:,1) -0.5; y = 1.16 *y; y = y +0.5; % stretch
        if ~isempty(find(strcmp(fieldnames(r),'irr')))
            y(r.irr) = NaN; % weed out irregular responses
            plot(ts(r.irr),  1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
            plot(ts(r.irr), -0.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
        end
    end
    plot(ts(2:end), y, '.', 'Color', [1 0.7 0]); % responses
    title(['Response y (orange), input u (green), learning rate (fine black), and posterior expectation of input s(\mu_2) ', ...
           '(red) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ...
           num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ...
      'FontWeight', 'bold');
    ylabel('y, u, s(\mu_2)');
    axis([0 ts(end) -0.15 1.15]);
else
    title(['Input u (green), learning rate (fine black), and posterior expectation of input s(\mu_2) ', ...
           '(red) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ...
           num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ...
      'FontWeight', 'bold');
    ylabel('u, s(\mu_2)');
    axis([0 ts(end) -0.1 1.1]);
end
plot(ts(2:end), 0.5, 'k');
xlabel('Trial number');
hold off;

author: translationalneuromodeling

repo: tapas


text:
% ANISODIFF - Anisotropic diffusion.
%
%
%  diff = anisodiff(im, niter, kappa, lambda, option)
%
%
%         im     - input image
%         niter  - number of iterations.
%         kappa  - conduction coefficient 20-100 ?
%         lambda - max value of .25 for stability
%         option - 1 Perona Malik diffusion equation No 1
%                  2 Perona Malik diffusion equation No 2
%
% Return
%         diff   - diffused image.
%
% kappa controls conduction as a function of gradient.  If kappa is low
% then mall intensity gradients are able to block conduction and hence diffusion
% across step edges.  A large value reduces the influence of intensity
% gradients on conduction.
%
% lambda controls speed of diffusion (you usually want it at a maximum of
% 0.25)
%
% Diffusion equation 1 preserve high contrast edges over low contrast ones.
% Diffusion equation 2 favours wide regions over smaller ones.

% Reference:
% P. Perona and J. Malik.
% Scale-space and edge detection using anisotropic diffusion.
% IEEE Transactions on Pattern Analysis and Machine Intelligence,
% 12(7):629-639, July 1990.

function diff = anisodiff(im, niter, kappa, lambda, option)

if ndims(im)==3
  error('Anisodiff only operates on 2D grey-scale images');
end

im = double(im);
[rows,cols] = size(im);
diff = im;

%{
var = 2;
x = (-4:4);
g = exp(-x.*x/(2*var)); g  = g/sum(g);

blurred = conv2(im,g,'same');
im_b = conv2(blurred,g','same');

%}

for i = 1:niter
 % fprintf('\rIteration %d',i);

  % Construct diffl which is the same as diff but
  % has an extra padding of zeros around it.
  diffl = zeros(rows+2, cols+2);
  diffl(2:rows+1, 2:cols+1) = diff;

  % North, South, East and West differences
  deltaN = diffl(1:rows,2:cols+1)   - diff;

  deltaS = diffl(3:rows+2,2:cols+1) - diff;

  deltaE = diffl(2:rows+1,3:cols+2) - diff;

  deltaW = diffl(2:rows+1,1:cols)   - diff;
  %deltaN = diff;deltaW;

  % Conduction

  if option == 1
    cN = exp(-(deltaN/kappa).^2);

    cS = exp(-(deltaS/kappa).^2);
    cE = exp(-(deltaE/kappa).^2);
    cW = exp(-(deltaW/kappa).^2);

  elseif option == 2
    cN = 1./(1 + (deltaN/kappa).^2);
    cS = 1./(1 + (deltaS/kappa).^2);
    cE = 1./(1 + (deltaE/kappa).^2);
    cW = 1./(1 + (deltaW/kappa).^2);
  end

  % APPLYING FOUR-POINT-TEMPLETE FOR numerical solution of DIFFUSION P.D.E.

  diff = diff + lambda*(cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW);
figure();
%  Uncomment the following to see a progression of images
 subplot(ceil(sqrt(niter)),ceil(sqrt(niter)), i)
 figure();
 pause(2)
imagesc(diff), colormap(gray), axis image

end
fprintf('\n');


text:
function symb_pvec = amplexpr(pvec)
%AMPLEXPR Converts SDPVAR variable to AMPL string

for pi = 1:size(pvec,1)
    for pj = 1:size(pvec,2)
        p = pvec(pi,pj);

        if isa(p,'double')
            symb_p = num2str(p);
        else
            LinearVariables = depends(p);
            x = recover(LinearVariables);
            exponent_p = full(exponents(p,x));
            names = cell(length(LinearVariables),1);
            for i = 1:length(LinearVariables)
                names{i}=['x[' num2str(LinearVariables(i)) ']'];
            end

            symb_p = '';
            if all(exponent_p(1,:)==0)
                symb_p = num2str(getbasematrix(p,0));
                exponent_p = exponent_p(2:end,:);
            end

            for i = 1:size(exponent_p,1)
                coeff = getbasematrixwithoutcheck(p,i);
                switch full(coeff)
                    case 1
                        coeff='+';
                    case -1
                        coeff = '-';
                    otherwise
                        if coeff >0
                            coeff = ['+' num2str2(coeff)];
                        else
                            coeff=num2str2(coeff);
                        end
                end
                if strcmp(symb_p,'') & (strcmp(coeff,'+') | strcmp(coeff,'-'))
                    symb_p = [symb_p coeff symbmonom(names,exponent_p(i,:))];
                else
                    symb_p = [symb_p coeff '*' symbmonom(names,exponent_p(i,:))];
                end
            end
            if symb_p(1)=='+'
                symb_p = symb_p(2:end);
            end
        end

        symb_p = strrep(symb_p,'+*','+');
        symb_p = strrep(symb_p,'-*','-');
        symb_pvec{pi,pj} = symb_p;
    end
end

function s = symbmonom(names,monom)
s = '';
for j = 1:length(monom)
    if monom(j)>0
        if strcmp(s,'')
            s = [s names{j}];
        else
            s = [s '*' names{j}];
        end
    end
    if monom(j)>1
        s = [s '^' num2str(monom(j))];
    end
end

function s = num2str2(x)
s = num2str(x);
if isequal(s,'1')
    s = '';
end
if isequal(s,'-1')
    s = '-';
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [hFig, iDS, iFig] = view_statcluster( StatFile, DisplayMode, Modality, hFig )
% VIEW_STATCLUSTER: Display the significant clusters in cluster-based stat file.
%
% USAGE: [hFig, iDS, iFig] = view_statcluster(StatFile, DisplayMode=[], Modality=[], hFig=[])
%
% INPUT:
%     - StatFile    : Relative path to stat file to visualize
%     - DisplayMode : {'clustindex_time', 'clustsize_time', 'longest'}
%     - hFig        : If defined, display file in existing figure
%
% OUTPUT :
%     - hFig : Matlab handle to the figure that was created or updated
%     - iDS  : DataSet index in the GlobalData variable
%     - iFig : Indice of returned figure in the GlobalData(iDS).Figure array

% @=============================================================================
% 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, 2015

global GlobalData;

% ===== READ INPUTS =====
if (nargin < 4) || isempty(hFig)
    hFig = [];
end
if (nargin < 3) || isempty(Modality)
    Modality = [];
end
if (nargin < 2) || isempty(DisplayMode)
    DisplayMode = [];
end
iDS = [];
iFig = [];

% ===== LOAD DATA =====
FileType = file_gettype(StatFile);
switch (FileType)
    case 'pdata'
        % Load data file
        iDS = bst_memory('LoadDataFile', StatFile);
        if isempty(iDS)
            return;
        end
        bst_memory('LoadRecordingsMatrix', iDS);
        % Get clusters
        StatClusters = GlobalData.DataSet(iDS).Measures.StatClusters;
        % Row names = channel names
        RowNames = {GlobalData.DataSet(iDS).Channel.Name};
        nRows = length(RowNames);
        % Default display
        if isempty(DisplayMode)
            DisplayMode = 'clustindex_time';
        end
        FreqLabels = [];

    case 'presults'
        % Load results file
        [iDS, iRes] = bst_memory('LoadResultsFileFull', StatFile);
        if isempty(iDS)
            return;
        end
        % Get clusters
        StatClusters = GlobalData.DataSet(iDS).Results(iRes).StatClusters;
        nRows = size(GlobalData.DataSet(iDS).Results(iRes).ImageGridAmp,1);
        RowNames = 1:nRows;
        % Default display
        if isempty(DisplayMode)
            DisplayMode = 'clustsize_time';
        end
        FreqLabels = [];

    case 'ptimefreq'
        % Load timefreq file
        [iDS, iTf] = bst_memory('LoadTimefreqFile', StatFile);
        if isempty(iDS)
            return;
        end
        % Get clusters
        StatClusters = GlobalData.DataSet(iDS).Timefreq(iTf).StatClusters;
        RowNames     = GlobalData.DataSet(iDS).Timefreq(iTf).RowNames;
        nRows = length(RowNames);
        % Default display
        if isempty(DisplayMode)
            DisplayMode = 'clustsize_time';
        end
        FreqLabels = panel_freq('FormatFreqLabels', GlobalData.DataSet(iDS).Timefreq(iTf).Freqs);

    case 'pmatrix'
        % Load matrix file
        [iDS, iMat] = bst_memory('LoadMatrixFile', StatFile);
        if isempty(iDS)
            return;
        end
        % Get clusters
        StatClusters = GlobalData.DataSet(iDS).Matrix(iMat).StatClusters;
        RowNames     = GlobalData.DataSet(iDS).Matrix(iMat).Description;
        nRows = length(RowNames);
        % Default display
        if isempty(DisplayMode)
            DisplayMode = 'clustindex_time';
        end
        FreqLabels = [];

    otherwise
        error('Unsupported file type.');
end
% No clusters available
if isempty(StatClusters)
    error(['No cluster information available from this file:' StatFile]);
end
% Get time vector
Time = bst_memory('GetTimeVector', iDS);

% ===== THRESHOLD CLUSTER MAPS =====
% Get significant clusters
[sClusters, PosClust, NegClust] = panel_stat('GetSignificantClusters', StatClusters);
% Set all the values that are not is a signficant cluster to zero
if isempty(PosClust)
    StatClusters.posclusterslabelmat = zeros(nRows, length(Time));
else
    StatClusters.posclusterslabelmat(~ismember(StatClusters.posclusterslabelmat, [PosClust.ind])) = 0;
end
if isempty(NegClust)
    StatClusters.negclusterslabelmat = zeros(nRows, length(Time));
else
    StatClusters.negclusterslabelmat(~ismember(StatClusters.negclusterslabelmat, [NegClust.ind])) = 0;
end

% ===== BUILD CLUSTER IMAGE =====
% Switch display mode
switch lower(DisplayMode)
    case 'clustindex_time'
        % Combine the two maps into one
        ClusterMap = StatClusters.posclusterslabelmat - StatClusters.negclusterslabelmat;
        % Apply function to map cluster indices with colors in the same way as the panel Stat
        clustPos = ClusterMap(ClusterMap > 0);
        ClusterMap(ClusterMap > 0) = max(clustPos(:)) - clustPos + 1;
        clustNeg = abs(ClusterMap(ClusterMap < 0));
        ClusterMap(ClusterMap < 0) = -(max(clustNeg(:)) - clustNeg + 1);
        % Replace non-stat values with NaN
        ClusterMap(ClusterMap == 0) = NaN;
        % Create the labels
        Labels = cell(1,4);
        Labels{1} = RowNames;
        Labels{2} = [];
        Labels{3} = Time;
        Labels{4} = FreqLabels;
        % Create the image volume: [N1 x N2 x Ntime x Nfreq]
        M = reshape(ClusterMap, size(ClusterMap,1), 1, size(ClusterMap,2), size(ClusterMap,3));
        % Show the image
        [hFig, iDS, iFig] = view_image_reg(M, Labels, [1,3], {'Signals','Time (s)'}, StatFile, hFig, 'cluster', 1, '$freq');
        % Configure the colormap to show positive and negative values
        bst_colormaps('SetColormapAbsolute', 'Image', 0);

    case 'clustsize_time'
        % For each cluster, compute the number of channels/sources involved at each time point
        if ~isempty(sClusters)
            F = zeros(length(sClusters), length(Time));
            LinesColor = cell(length(sClusters), 1);
            LinesLabels = cell(length(sClusters), 1);
            for i = 1:length(sClusters)
                F(i,:) = sum(sum(double(sClusters(i).mask), 1), 3);
                LinesColor{i} = sClusters(i).color;
                LinesLabels{i} = sprintf('p=%1.5f, c=%d, s=%d', sClusters(i).prob, round(sClusters(i).clusterstat), sClusters(i).clustsize);
            end
        else
            F = zeros(1, length(Time));
            LinesColor{1} = [0 0 0];
            LinesLabels{1} = 'Empty';
        end
        % Display signals
        [hFig, iDS, iFig] = view_timeseries_matrix(StatFile, F, [], 'Clusters', 'Cluster size', LinesLabels, LinesColor, hFig, []);

    % Topography or source map with longest values for each cluster
    case 'longest'
        % Plot topography/source map
        switch (FileType)
            case 'pdata'
                % For each cluster, compute the number of significant time points for each channel/source
                if ~isempty(sClusters)
                    F = zeros(nRows, 1);
                    for i = 1:length(sClusters)
                        [lenMax, iRow] = max(sum(sClusters(i).mask,2));
                        F(iRow,:) = lenMax;
                    end
                else
                    F = zeros(nRows, 1);
                end

                % Get channels from selected modality
                iChannels = channel_find(GlobalData.DataSet(iDS).Channel, Modality);
                % Display topography
                [hFig, iDS, iFig] = view_topography(StatFile, Modality, '2DSensorCap', F(iChannels,:), 0, 'NewFigure', []);
                % Set colormap
                bst_colormaps('ConfigureColorbar', hFig, 'timefreq', 'samples', 'samples');
                ColormapInfo = getappdata(hFig, 'Colormap');
                ColormapInfo.Type = 'timefreq';
                ColormapInfo.AllTypes = {'timefreq'};
                ColormapInfo.DisplayUnits = 'samples';
                setappdata(hFig, 'Colormap', ColormapInfo);

            case 'presults'
                % Display surface
                [hFig, iDS, iFig] = view_surface_data(GlobalData.DataSet(iDS).Results(iRes).SurfaceFile, StatFile, [], hFig);
                % Get axes object
                hAxes = findobj(hFig, '-depth', 1, 'Tag', 'Axes3D');
                % Get surface in memory
                % sSurf = bst_memory('GetSurface', GlobalData.DataSet(iDS).Results(iRes).SurfaceFile);
                % Get displayed surface
                hPatch = findobj(hAxes, '-depth', 1, 'Tag', 'AnatSurface');
                Vertices = get(hPatch, 'Vertices');
                % Remove previous markers
                delete(findobj(hAxes, '-depth', 1, 'Tag', 'ptClusterMax'));
                % Search maximum duration for each cluster
                for i = 1:length(sClusters)
                    % Get maximum
                    [lenMax, iRow] = max(sum(sClusters(i).mask,2));
                    if (lenMax == 0)
                        continue;
                    end
                    % Get point coordinates
                    pt = Vertices(iRow,:);
                    % Mark new point
                    line(pt(1)*1.005, pt(2)*1.005, pt(3)*1.005, ...
                         'MarkerFaceColor', sClusters(i).color, ...
                         'MarkerEdgeColor', sClusters(i).color, ...
                         'Marker',          '+',  ...
                         'MarkerSize',      12, ...
                         'LineWidth',       2, ...
                         'Parent',          hAxes, ...
                         'Tag',             'ptClusterMax');
                end

            otherwise
                error('Not supported.');
        end
    otherwise
        error(['Invalid display mode: "' DisplayMode '"']);
end

% Add tag in the figure appdata
StatInfo.StatFile    = StatFile;
StatInfo.DisplayMode = DisplayMode;
setappdata(hFig, 'StatInfo', StatInfo);

% Show "stat" tab
gui_brainstorm('ShowToolTab', 'Stat');
% Update figure selection
bst_figures('SetCurrentFigure', hFig, '2D');

% Save reload call
ReloadCall = {'view_statcluster', StatFile, DisplayMode, Modality, hFig};
setappdata(hFig, 'ReloadCall', ReloadCall);

% Stat panel: reload clusters list
panel_stat('CurrentFigureChanged_Callback', hFig);

% Last update for topography
if strcmpi(DisplayMode, 'longest')
    if strcmpi(FileType, 'pdata')
        figure_topo('ColormapChangedCallback', iDS, iFig);
        figure_topo('CurrentTimeChangedCallback', iDS, iFig);
%     elseif strcmpi(FileType, 'presults')
%         panel_surface('UpdateSurfaceData', hFig, 1);
    end
end

author: brainstorm-tools


text:
function [Z,H,T,perm] = dendrogram_subtreepixels(data, method, p , x_mu, y_mu, savethis, nameappend)
% [Z,H,T,perm] = dendrogram_subtreepixels(data, method, p , x_mu, y_mu, savethis, nameappend)

% dveccr= reshape(dpixc,15^2, 1000);
if ~exist('savethis', 'var')
    savethis = 0;
end
if ~exist('nameappend', 'var')
    nameappend=[];
else
    nameappend=['_' nameappend];
end

sized = size(data);
dveccr= reshape(data,sized(1)*sized(2), sized(3));
ccd = (corrcoef(dveccr'));
e=eye(size(ccd));
ccdflip=e-(ccd-e);
ccds=squareform(1-ccdflip);

Z = linkage(ccds,method);
[H,T,perm] = dendrogram(Z,0);
co= get(gca, 'colororder');
set(H,'color','k')
ylim([0 max(Z(:,3))])

zvalues = input('Z values from dendrogram[val1*10^4 val2*10^4 ...]\n');
clear('cl')
for ii=1:length (zvalues)
    cl(ii).ixZ = find(round(10^4*Z(:,3))==zvalues(ii));
end

imz2 = zeros([size(dveccr,1) 1]);
[H,T,perm] = dendrogram(Z,0);
co= get(gca, 'colororder');
set(H,'color','k')
ylim([0 max(Z(:,3))])

lcl=length(cl);
lco=length(co);
rat = lcl/lco;

if rat > 1 %more clusters than colors
    co = repmat(co, ceil(rat), 1);
end

for ii=1:lcl
    [cl(ii).ixZ_vec, cl(ii).endleaves_vec] = recursivesubtree(Z, cl(ii).ixZ, [], []);
    set(H(cl(ii).ixZ_vec), 'color', co(ii,:))
    imz2(cl(ii).endleaves_vec, ii)=1;
end
set(H,'linewidth',2)

if savethis
    name=['dendrogram-' method nameappend];
    fprintf('saving the dendrogram figure:'' %s ''\n',name)
    SaveImageFULL(name)
    save ([name '_zvalues'], 'zvalues')
end

imzRGB = reshape(imz2*co(1:lcl,:),sized(1), sized(2), 3);
dipshow(joinchannels('RGB', imzRGB));
hold on
if exist ('p', 'var')
    if ~isempty(p)
%         scatter(p.x_vec-0.5, p.y_vec-0.5,200,'xw')
%         scatter(p.x_vec, p.y_vec,'w')
        scatter(p.x_vec, p.y_vec,[],[.9 .9 .9])
    end
end
if and(exist ('x_mu', 'var'), exist ('y_mu', 'var'))
%     scatter(x_mu-1, y_mu-1,[],'xw')
    scatter(x_mu-1, y_mu-1,[],[.9 .9 .9],'x')
end
if savethis
    name=['clustfig-' method nameappend];
    fprintf('saving the clustered figure:'' %s ''\n',name)
    SaveImageFULL(name)
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [lin_fun, pp] = RAWCRF(image_raw, image_jpg, N, threshold_outliers)
%
%       [lin_fun, pp] = RAWCRF(image_raw, image_jpg, N, threshold_outliers)
%
%
%
%        Input:
%           -image_raw:
%           -image_raw:
%           -N:
%           -threshold_outliers:
%
%        Output:
%           -lin_fun:
%           -pp:
%

if(isempty(image_raw))
    error('RAWCRF: a stack cannot be empty!');
end

if(isempty(image_jpg))
    error('RAWCRF: a stack_exposure cannot be empty!');
end

if(~exist('threshold', 'var'))
    threshold_outliers = 0.05;
end

if(~exist('N', 'var'))
    N = -1;
end

threshold_outliers_inv = 1.0 - threshold_outliers;

col = size(image_raw, 3);

lin_fun = zeros(256, col);

err = -1;
N_max = 6;
thr = [threshold_outliers, threshold_outliers_inv];

if(N < 0)
    for N_tmp=1:N_max
        pp = RAWCRFn(image_raw, image_jpg, N_tmp, thr);

        imgOut = RemoveCRF(image_jpg, 'poly', pp);

        tmp_err = abs(imgOut - image_raw).^2;
        tmp_err = mean(tmp_err(:));

        if(err < 0)
            err = tmp_err;
            N = N_tmp;
        else
            if(tmp_err < err)
                err = tmp_err;
                N = N_tmp;
            end
        end

    end
else
    pp = RAWCRFn(image_raw, image_jpg, N, thr);
end

x_val = (0:255) / 255;

for i=1:col
    lin_fun(:, i) = polyval(pp(:, i), x_val);
end

end

function pp = RAWCRFn(image_raw, image_jpg, N, thr)
    col = size(image_raw, 3);

    pp = zeros(N + 1, col);

    bFlag = 0;

    %bDebug = 1;

    for i=1:col
        slice_raw = image_raw(:,:,i);
        slice_jpg = image_jpg(:,:,i);

        indx = find(slice_jpg > thr(1) & slice_jpg < thr(2));

        if(~isempty(indx))
            x = slice_jpg(indx);
            y = slice_raw(indx);

%
%             if(bDebug)
%                 figure(i);
%                 c = zeros(256,256);
%                 for k =1:length(x)
%                     tx = ClampImg(round(x(k) * 255) + 1, 1, 256);
%                     ty = ClampImg(round(y(k) * 255) + 1, 1, 256);
%                     c(256 - ty + 1, tx) = c(256 - ty + 1, tx) + 1;
%                 end
%                 imshow(c);
%             end

            pp(:, i) = polyfit(x, y, N);
        else
            bFlag = 1;
            disp('RAWCRF: no enough data for estimating the CRF');
        end
    end

    if(bFlag)
        pp = zeros(N + 1, col);
    end
end

author: banterle


text:
function var = loadgenericdata(varfile, vartype, fields, varname, args)
%LOADGENERICDATA Loads data of a specified type from file or data structure.
%
%   VAR = LOADGENERICDATA(VARFILE, VARTYPE)
%   VAR = LOADGENERICDATA(VARFILE, VARTYPE, FIELDS)
%   VAR = LOADGENERICDATA(VARFILE, VARTYPE, FIELDS, VARNAME)
%   VAR = LOADGENERICDATA(VARFILE, VARTYPE, FIELDS, VARNAME, ARGS)
%
%   Loads data from a variable or M-file or MAT-file and checks that it
%   matches a specified type.
%
%   Inputs:
%       VARFILE : Variable containing the data structure or a string
%                 containing the name of a function M-file or a MAT-file
%                 on the MATLAB/Octave path. If no file extension is provided
%                 it will attempt to load a MAT-file with that name and,
%                 if not found, will call a function by that name to
%                 get the data. The function M-file should return a
%                 single argument containing the data. A MAT-file should
%                 either contain a single variable with the desired data
%                 or provide the variable name in VARNAME.
%       VARTYPE : String or cell array of strings with, in order of
%                 priority, the data structure type to be returned.
%                 Valid values are 'struct', 'cell' and 'array'.
%       FIELDS :  (optional) String or cell array of strings containing
%                 a list of required fields in case the VARTYPE is struct.
%                 If a required field is missing it will throw an error.
%       VARNAME   (optional) String containing the name of the variable to
%                 extract when loading a MAT-file. If not provided, the
%                 default is to extract the 1st variable, regardless of name.
%       ARGS :    (optional) Scalar or cell array of values that are passed
%                 as input function arguments to VARFILE if it is an M-file.
%
%   Output:
%       VAR : Returned data structure of the first matching VARTYPE.

%   MOST
%   Copyright (c) 2013-2016, Power Systems Engineering Research Center (PSERC)
%   by Daniel Munoz-Alvarez and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

if isempty(varfile)
    var = [];
    return;
end
if nargin < 5
    args = [];
    if nargin < 4
        varname = '';
        if nargin < 3
            fields = {};
        end
    end
end

if ischar(varfile)      %% it's a file name
    %% check for explicit extension
    l = length(varfile);
    if l > 2
        if strcmp(varfile(l-1:l), '.m')
            rootname = varfile(1:l-2);
            extension = '.m';
        elseif l > 4
            if strcmp(varfile(l-3:l), '.mat')
                rootname = varfile(1:l-4);
                extension = '.mat';
            end
        end
    end

    %% set extension if not specified explicitly
    if ~exist('rootname', 'var')
        rootname = varfile;
        if exist([rootname '.mat'], 'file') == 2
            extension = '.mat';
        elseif exist([rootname '.m'], 'file') == 2
            extension = '.m';
        else
            error('loadgenericdata: No file named ''%s.mat'' or ''%s.m'' found in MATLAB/Octave path.', rootname, rootname);
        end
    end

    %% attempt to read file
    if strcmp(extension,'.mat')         %% from MAT file
        s = load(rootname);
        if isempty(varname)
            f = fieldnames(s);
            var = s.(f{1});             %% use first variable if not specified
        elseif isfield(s, varname)
            var = s.(varname);          %% use specified variable
        else
            error('loadgenericdata: No variable named ''%s'' in ''%s.mat''', varname, rootname);
        end
    elseif strcmp(extension,'.m')       %% from M file
        if isempty(args)
            var = feval(rootname);
        elseif iscell(args)
            var = feval(rootname, args{:});
        else
            var = feval(rootname, args);
        end
    end
else                    %% it's a data structure, not a file name
    var = varfile;
end

%% check the provided types
if ~iscell(vartype)
    vartype = {vartype};    %% convert string to single element cell array
end
if ~iscell(fields) && ~isempty(fields)
    fields = {fields};      %% convert string to single element cell array
end
done = false;
for i = length(vartype)
    switch vartype{i}
        case 'struct'
            if isstruct(var)
                for j = 1:length(fields)
                    if ~ischar(fields{j})
                        error('loadgenericdata: FIELDS input must be a string or cell array of strings.')
                    end
                    if ~isfield(var, fields{j})
                        error('loadgenericdata:missingfield', 'loadgenericdata: field ''%s'' missing from struct data.', fields{j});
                    end
                end
                done = true; break;     %% struct OK
            end

        case 'cell'
            if iscell(var)
                done = true; break;     %% cell array OK
            end

        case 'array'
            if isnumeric(var)
                done = true; break;     %% array OK
            end
        otherwise
            error('loadgenericdata: ''%s'' is not a valid value for VARTYPE', vartype{i});
    end
end
if ~done
    error('loadgenericdata: data not found in any of the specified formats');
end


text:
function [Vertices_sm, A] = tess_smooth(Vertices, a, nIterations, VertConn, isKeepSize, Faces)
% TESS_SMOOTH: Smooths a surface.
%
% USAGE:  [Vertices_sm, A] = tess_smooth(Vertices, a, nIterations, VertConn, isKeepSize)
%         [Vertices_sm, A] = tess_smooth(Vertices, a, nIterations, VertConn)
%
% INPUT:
%    - Vertices    : [N,3] vertices of matrix to smooth
%    - a           : scalar smooth weighting parameter (0-1 less-more smoothing)
%    - nIterations : number of times to apply the smoothing
%    - VertConn    : Vertex connectivity sparse matrix
%    - isKeepSize  : If 1, the final surface is scaled so that the convex envelope
%                    has the same bounding box as the initial surface
% OUTPUT:
%    - Vertices_sm : vertices list of smoothed surface
%    - A           : smoothing matrix

% @=============================================================================
% 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, 2012-2013

% Parse inputs
if (nargin < 5) || isempty(isKeepSize)
    isKeepSize = 1;
end
% Check matrices orientation
if (size(Vertices, 2) ~= 3)
    error('Vertices must have 3 columns (X,Y,Z).');
end

% Get initial bounding box
if isKeepSize
    initBounds = [min(Vertices); max(Vertices)];
end
% Calculate smoothing matrix
A = spones(VertConn);
sumA = sum(A);
sumA(sumA==0) = eps;
A = spdiags((a./sumA)', 0, size(A,1), size(A,2)) * A;
A = A + spdiags((1-a) * ones(size(A,1), 1), 0, size(A,1), size(A,2));

% Smooth vertices matrix
Vertices_sm = double(Vertices);
for i = 1:nIterations
    Vertices_sm = A * Vertices_sm;
end

% Scale final surface
if isKeepSize
    % Compute scale
    finalBounds = [min(Vertices_sm); max(Vertices_sm)];
    initBounds = initBounds - repmat(mean(Vertices),2,1);
    scale = diff(initBounds,[],1) ./ diff(finalBounds,[],1);
    % Center and apply scale
    center = mean(Vertices_sm);
    Vertices_sm = bst_bsxfun(@minus, Vertices_sm, center);
    Vertices_sm = bst_bsxfun(@times, Vertices_sm, scale);
    Vertices_sm = bst_bsxfun(@plus,  Vertices_sm, center);

%     % Compute normals
%     VertNormals = tess_normals(Vertices_sm, Faces, VertConn);
%     Vertices_sm = Vertices_sm + 0.00001 * nIterations * VertNormals;
end


text:
%%*******************************************************************
%% NTrhsfun: compute the right-hand side vector of the
%%           Schur complement equation for the NT direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [rhs,EinvRc,hRd] = NTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ)

spdensity = par.spdensity;
m = length(rp);
if (nargin > 8)
    corrector = 1;
else
    corrector = 0;
    hRd = zeros(m,1);
end
hEinvRc = zeros(m,1);
EinvRc  = cell(size(blk,1),1);
rhsfree = [];
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});  numblk = length(pblk{2});
    if strcmp(pblk{1},'l')
        if iscell(sigmu)
            EinvRc{p} = sigmu{p}./Z{p} -X{p};
        else
            EinvRc{p} = sigmu./Z{p} -X{p};
        end
        Rq = sparse(n,1);
        if (corrector) && (norm(par.parbarrier{p})==0)
            Rq = dX{p}.*dZ{p}./Z{p};
        else
            tmp  = par.dd{p}.*Rd{p};
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = EinvRc{p} - Rq;
        tmp2 = mexMatvec(At{p},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'q')
        if iscell(sigmu)
            EinvRc{p} = qops(pblk,-sigmu{p}./(par.gamz{p}.*par.gamz{p}),Z{p},4) -X{p};
        else
            EinvRc{p} = qops(pblk,-sigmu./(par.gamz{p}.*par.gamz{p}),Z{p},4) -X{p};
        end
        Rq = sparse(n,1);
        if (corrector) && (norm(par.parbarrier{p})==0)
            w = sqrt(par.gamz{p}./par.gamx{p});
            hdx = qops(pblk,w,par.ff{p},5,dX{p});
            hdz = qops(pblk,w,par.ff{p},6,dZ{p});
            hdxdz = Arrow(pblk,hdx,hdz);
            vv = qops(pblk,w,par.ff{p},5,X{p});
            Vihdxdz = Arrow(pblk,vv,hdxdz,1);
            Rq = qops(pblk,w,par.ff{p},6,Vihdxdz);
        else
            tmp  = par.dd{p}.*Rd{p} + qops(pblk,qops(pblk,Rd{p},par.ee{p},1),par.ee{p},3);
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = EinvRc{p} - Rq;
        tmp2 = mexMatvec(At{p},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'s')
        n2 = pblk{2}.*(pblk{2}+1)/2;
        if iscell(sigmu)
            %%ss = [0,cumsum(pblk{2})];
            %%sigmuvec = zeros(n,1);
            %%for k = 1:length(pblk{2});
            %%   sigmuvec(ss(k)+1:ss(k+1)) = sigmu{p}(k)*ones(pblk{2}(k),1);
            %%end
            sigmuvec = mexexpand(pblk{2},sigmu{p});
            tmp = spdiags(sigmuvec./par.sv{p} -par.sv{p},0,n,n);
        else
            tmp = spdiags(sigmu./par.sv{p} -par.sv{p},0,n,n);
        end
        EinvRc{p} = Prod3(pblk,par.G{p}',tmp,par.G{p},1);
        Rq = sparse(n,n);
        if (corrector) && (norm(par.parbarrier{p})==0)
            hdZ = Prod3(pblk,par.G{p},dZ{p},par.G{p}',1);
            hdX = spdiags(qops(pblk,par.parbarrier{p}',1./par.sv{p},3)-par.sv{p},0,n,n)-hdZ;
            tmp = Prod2(pblk,hdX,hdZ,0);
            tmp = 0.5*(tmp+tmp');
            if (numblk == 1)
                d = par.sv{p};
                e = ones(pblk{2},1);
                Rq = 2*tmp./(d*e'+e*d');
                if (nnz(Rq) <= spdensity*n2); Rq = sparse(Rq); end
            else
                Rq = sparse(n,n);
                ss = [0, cumsum(pblk{2})];
                for i = 1:numblk
                    pos = ss(i)+1 : ss(i+1);
                    d = par.sv{p}(pos); e = ones(length(pos),1);
                    Rq(pos,pos) = 2*tmp(pos,pos)./(d*e' + e*d'); %#ok
                end
            end
            Rq = Prod3(pblk,par.G{p}',Rq,par.G{p},1);
        else
            tmp = Prod3(pblk,par.W{p},Rd{p},par.W{p},1,par.nzlistAy{p});
            tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),{tmp});
            hRd = hRd + tmp2;
        end
        EinvRc{p} = EinvRc{p} - Rq;
        tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'u')
        rhsfree = [rhsfree; Rd{p}]; %#ok
    end
end
%%
rhs = rp + hRd - hEinvRc;
rhs = full([rhs; rhsfree]);
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function sample_runscript_matlab(doParallelize,saveIncrement,fileName)
% SAMPLE_RUNSCRIPT_MATLAB  Template function for looping over a hctsa analysis when
%                           using local Matlab files (initializing with TS_init).
%
%---INPUTS:
% doParallelize, (binary) whether to parallelize computations over available
%                   CPUs using Matlab's Parallel Computing Toolbox (default: 1).
% saveIncrement, saves the results back to the hctsa file after computing for
%               this many time series (default: 5)
% fileName, the source of the HCTSA dataset (a .mat file)
%
%---EXAMPLE USAGE:
% sample_runscript_matlab(1,'HCTSA.mat');

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Set parameters for computation:
%-------------------------------------------------------------------------------
if nargin < 1
    % Set to 1 to
    doParallelize = true;
end

if nargin < 2
    saveIncrement = 5;
end

if nargin < 3
    fileName = 'HCTSA.mat';
end

% ------------------------------------------------------------------------------
%% Preparations:
% ------------------------------------------------------------------------------
% Get IDs of all time series in the hctsa datafile:
[~,TimeSeries,Operations] = TS_LoadData(fileName);
tsIDs = TimeSeries.ID;

% Set the IDs to compute at each iteration:
ID_inc = 1:saveIncrement:length(tsIDs)+1;

% Provide a quick message:
fprintf(1,['About to calculate across %u time series and %u operations ' ...
                'saving back to file every %u time series\n'], ...
                height(TimeSeries),height(Operations),saveIncrement);

%-------------------------------------------------------------------------------
% Running calculations:
%-------------------------------------------------------------------------------
% Loop across time series:
for i = 1:length(ID_inc)-1
	fprintf(1,'\n\n\nWe''re looking at time series with IDs from %u--%u\n\n\n', ...
                            	tsIDs(ID_inc(i)),tsIDs(ID_inc(i+1)-1));

    % Compute any missing values for this range of time series, then save back:
    TS_Compute(doParallelize,tsIDs(ID_inc(i):ID_inc(i+1)-1),[],'missing',fileName);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function R = sw_han_converse(a, n_vals, eps)
% Compute the converse sum rate (at the symmetrical rate point) from Han's SW converse
% for the joint source distribution: [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)].

% Input:
% a: parameter to define the joint source distribution, should be in range [0.25, 1)
% n_vals: a vector containing the blocklengths
% eps: target error probability
% Output:
% R: the rates for the blocklengths specified in n_vals.

    p = [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)];                                 % Stores joint probabilities
    H = sum(sum(p.*log(1./p)));                                             % Joint entropy
    b = [a/(a+1/3*(1-a)) 1/3*(1-a)/(2*1/3*(1-a)) 1/3*(1-a)/(a+1/3*(1-a))];  % Stores conditional probabilities

    % Compute log binomial coefficients
    log_b = binomial_coeff(max(n_vals));

    R = zeros(length(n_vals),1);

    for i = 1:length(n_vals)
        n = n_vals(i);
        m = 0:n;
        log_Pr = m*log(a)+(n-m)*log(1/3*(1-a));

        % Search the range of gamma specified in r_vals
        % The range of r_vals could be set smaller for large n to save computation time.
        if n <= 500
            r_vals = 0.0001:0.00005:2;
        else
            r_vals = 0.0001:0.00005:0.1;
        end

        R_r = zeros(1,length(r_vals));
        for ri = 1:length(r_vals)
            r = r_vals(ri);
            x = 0;
            y = H + 0.1;

            %%%%%%%%%% Use this part to check if the range is valid for bisection algorithm %%%%%%%%%%
            c = x;
            m0 = max(0, floor(n*(log(1/3*(1-a))+c+r)/log(1/3*(1-a)/a)));
            if m0 >= n
                err = 1 - 3*exp(-n*r);
            else
                log_err_m1 = my_logsumexp(log_Pr(1:m0+1)+log_b{n+1}(1:m0+1)+(n-(0:m0))*log(3));
                log_err_m2 = [];
                for m = m0+1:n
                    B_m = ceil((m*log(b(1))+(n-m)*log(b(2))+n*(c/2+r))/(log(b(2)/b(3))));
                    if B_m <= 0
                        C_m = log_Pr(m+1) + log_b{n+1}(m+1) + (n-m)*log(3);
                        log_err_m2 = [log_err_m2 C_m];
                    elseif B_m <= n-m
                        C_m1 = my_logsumexp(log_b{n-m+1}(B_m+1:n-m+1)+(n-m-(B_m:n-m))*log(2));
                        log_err_k = [];
                        for k = 0:B_m-1
                            if B_m <= n-m-k
                                D_m = log_b{n-m+1}(k+1) + my_logsumexp(log_b{n-m-k+1}(B_m+1:n-m-k+1));
                                log_err_k = [log_err_k D_m];
                            end
                        end
                        C_m2 = my_logsumexp(log_err_k);
                        C_m = log_Pr(m+1) + log_b{n+1}(m+1) + my_logsumexp([C_m1 C_m2]);
                        log_err_m2 = [log_err_m2 C_m];
                    end
                end
                log_err = my_logsumexp([log_err_m1 log_err_m2]);
                err = exp(log_err) - 3*exp(-n*r);
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % Start of bisection iteration
            iter = 0;
            while abs(x - y) >= 0.00001 && iter <= 30
                c = (x+y)/2;
                %%%%%%%%%%%%%%%%%%% body of iteration %%%%%%%%%%%%%%%%%%%%
                m0 = max(0, floor(n*(log(1/3*(1-a))+c+r)/log(1/3*(1-a)/a)));
                if m0 >= n
                    err = 1 - 3*exp(-n*r);
                else
                    log_err_m1 = my_logsumexp(log_Pr(1:m0+1)+log_b{n+1}(1:m0+1)+(n-(0:m0))*log(3));
                    log_err_m2 = [];
                    for m = m0+1:n
                        B_m = ceil((m*log(b(1))+(n-m)*log(b(2))+n*(c/2+r))/(log(b(2)/b(3))));
                        if B_m <= 0
                            C_m = log_Pr(m+1) + log_b{n+1}(m+1) + (n-m)*log(3);
                            log_err_m2 = [log_err_m2 C_m];
                        elseif B_m <= n-m
                            C_m1 = my_logsumexp(log_b{n-m+1}(B_m+1:n-m+1)+(n-m-(B_m:n-m))*log(2));
                            log_err_k = [];
                            for k = 0:B_m-1
                                if B_m <= n-m-k
                                    D_m = log_b{n-m+1}(k+1) + my_logsumexp(log_b{n-m-k+1}(B_m+1:n-m-k+1));
                                    log_err_k = [log_err_k D_m];
                                end
                            end
                            C_m2 = my_logsumexp(log_err_k);
                            C_m = log_Pr(m+1) + log_b{n+1}(m+1) + my_logsumexp([C_m1 C_m2]);
                            log_err_m2 = [log_err_m2 C_m];
                        end
                    end
                    log_err = my_logsumexp([log_err_m1 log_err_m2]);
                    err = exp(log_err) - 3*exp(-n*r);
                end
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                if err - eps < 0
                    y = c;
                else
                    x = c;
                end
                iter = iter + 1;
            end
            R_r(ri) = c;
        end
        % Optimize with respect to gamma
        R(i) = max(R_r);
    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function [y,fs]=readflac(filename,mode)
%READWAV  Read a .FLAC format sound file [Y,FS]=(FILENAME,MODE)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .WAV extension) or alternatively
%                 can be the FIDX output from a previous call to READWAV
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1
%             'r'    Raw unscaled data (integer values)
%             'q'    Scaled to make 0dBm0 be unity mean square
%             'p' *	 Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%
% FLAC (Free lossless audio codec) is a compressed audio file format described here:
% http://flac.sourceforge.net/

%      Copyright (C) Mike Brookes 2008
%      Version: $Id: readflac.m,v 1.1 2009/09/13 15:52:23 dmb Exp $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2
    mode='p';
else
    mode = [mode(:).' 'p'];
end
dirt=voicebox('dir_temp');
[fnp,fnn,fne,fnv]=fileparts(filename);
filetemp=fullfile(dirt,[fnn '.wav' fnv]);
doscom=['"' voicebox('flac') '"' ' -d -f -o "' filetemp '" "' filename '"'];
%                     fprintf(1,'Executing: %s\n',doscom);
[doss,dosr]=dos(doscom); % run the program
if doss % test for errors
    error(sprintf('Error running DOS command: %s',doscom));
end
if exist(filetemp)~=2
    error(sprintf('No output file from: %s',doscom));
end
[y,fs]=readwav(filetemp,mode);
doscom=['del /f "' filetemp '"'];
if dos(doscom) % run the program
    error(sprintf('Error running DOS command: %s',doscom));
end

author: decouples


text:
%% Machine Learning Online Class
%  Exercise 8 | Anomaly Detection and Collaborative Filtering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     estimateGaussian.m
%     selectThreshold.m
%     cofiCostFunc.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ================== Part 1: Load Example Dataset  ===================
%  We start this exercise by using a small dataset that is easy to
%  visualize.
%
%  Our example case consists of 2 network server statistics across
%  several machines: the latency and throughput of each machine.
%  This exercise will help us find possibly faulty (or very fast) machines.
%

fprintf('Visualizing example dataset for outlier detection.\n\n');

%  The following command loads the dataset. You should now have the
%  variables X, Xval, yval in your environment
load('ex8data1.mat');

%  Visualize the example dataset
plot(X(:, 1), X(:, 2), 'bx');
axis([0 30 0 30]);
xlabel('Latency (ms)');
ylabel('Throughput (mb/s)');

fprintf('Program paused. Press enter to continue.\n');
pause

%% ================== Part 2: Estimate the dataset statistics ===================
%  For this exercise, we assume a Gaussian distribution for the dataset.
%
%  We first estimate the parameters of our assumed Gaussian distribution,
%  then compute the probabilities for each of the points and then visualize
%  both the overall distribution and where each of the points falls in
%  terms of that distribution.
%
fprintf('Visualizing Gaussian fit.\n\n');

%  Estimate my and sigma2
[mu sigma2] = estimateGaussian(X);

%  Returns the density of the multivariate normal at each data point (row)
%  of X
p = multivariateGaussian(X, mu, sigma2);

%  Visualize the fit
visualizeFit(X,  mu, sigma2);
xlabel('Latency (ms)');
ylabel('Throughput (mb/s)');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================== Part 3: Find Outliers ===================
%  Now you will find a good epsilon threshold using a cross-validation set
%  probabilities given the estimated Gaussian distribution
%

pval = multivariateGaussian(Xval, mu, sigma2);

[epsilon F1] = selectThreshold(yval, pval);
fprintf('Best epsilon found using cross-validation: %e\n', epsilon);
fprintf('Best F1 on Cross Validation Set:  %f\n', F1);
fprintf('   (you should see a value epsilon of about 8.99e-05)\n');
fprintf('   (you should see a Best F1 value of  0.875000)\n\n');

%  Find the outliers in the training set and plot the
outliers = find(p < epsilon);

%  Draw a red circle around those outliers
hold on
plot(X(outliers, 1), X(outliers, 2), 'ro', 'LineWidth', 2, 'MarkerSize', 10);
hold off

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================== Part 4: Multidimensional Outliers ===================
%  We will now use the code from the previous part and apply it to a
%  harder problem in which more features describe each datapoint and only
%  some features indicate whether a point is an outlier.
%

%  Loads the second dataset. You should now have the
%  variables X, Xval, yval in your environment
load('ex8data2.mat');

%  Apply the same steps to the larger dataset
[mu sigma2] = estimateGaussian(X);

%  Training set
p = multivariateGaussian(X, mu, sigma2);

%  Cross-validation set
pval = multivariateGaussian(Xval, mu, sigma2);

%  Find the best threshold
[epsilon F1] = selectThreshold(yval, pval);

fprintf('Best epsilon found using cross-validation: %e\n', epsilon);
fprintf('Best F1 on Cross Validation Set:  %f\n', F1);
fprintf('   (you should see a value epsilon of about 1.38e-18)\n');
fprintf('   (you should see a Best F1 value of 0.615385)\n');
fprintf('# Outliers found: %d\n\n', sum(p < epsilon));

author: Ayatans

repo: Machine-Learning-homework


text:
classdef UF1 < PROBLEM
% <multi> <real> <large/none>
% Unconstrained benchmark MOP

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,
% Multiobjective optimization test instances for the CEC 2009 special
% session and competition, School of CS & EE, University of Essex, Working
% Report CES-487, 2009.
%------------------------------- 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 = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = [0,zeros(1,obj.D-1)-1];
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 3 : 2 : obj.D;
            J2 = 2 : 2 : obj.D;
            Y  = X - sin(6*pi*repmat(X(:,1),1,obj.D)+repmat(1:obj.D,size(X,1),1)*pi/obj.D);
            PopObj(:,1) = X(:,1)         + 2*mean(Y(:,J1).^2,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + 2*mean(Y(:,J2).^2,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
        	R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function confidences = test_boosted_dt_mc(classifier, features)
% confidences = test_boosted_dt_mc(classifier, features)
%
% Returns a log likelihod ratio for each class in the classifier
%
% Input:
%  classifier: boosted decision tree classifier
%  features:   classifier features (ndata, nvariables)
% Output:
%   confidences(ndata, nclasses):
%      P(class=k|features) \propto 1./(1+exp(-confidences(k)))

npred = classifier.wcs(1).dt.npred;
if size(features, 2)~=npred
    error('Incorrect number of attributes')
end

wcs = classifier.wcs;
nclasses = size(wcs, 2);

ntrees = size(wcs, 1);

confidences = zeros(size(features, 1), nclasses);
for c = 1:nclasses
    for t = 1:ntrees
        if ~isempty(wcs(t,c).dt)
            if 1
                dt = wcs(t,c).dt;
%                 dt = tree_getNewVersion(dt);
                [var, cut, children, catsplit] = tree_getParameters(dt);
                nodes = treevalc(int32(var), cut, int32(children(:, 1)), ...
                        int32(children(:, 2)), catsplit(:, 1), features');
                %disp(num2str(nodes));
            else
                [class_indices, nodes, classes] = treeval(wcs(t, c).dt, features);
            end
            confidences(:, c) = confidences(:, c) + wcs(t, c).confidences(nodes);
        end
    end
    confidences(:, c) = confidences(:, c) + classifier.h0(c);
end

author: zouchuhang

repo: LayoutNet


text:
% Diagram of van Albada limiter

% Theory in Section 10.8 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program makes Fig. 10.27 in the book

figure(1), clf, hold on

x = 0.0:0.01:4.0; y = (x.*x + x)./(1 + x.*x);	% Graph of van Albada limiter
plot(x,y)

x = 0.0:0.01:1.1; y = 0.5*(1+sqrt(2))*x; plot(x,y)
x = 0.0:0.01:4.0; y = 0.5*(1+sqrt(2))*ones(size(x)); plot(x,y)
y = ones(size(x)); plot(x,y)
x = 0.0:0.01:1.3; plot(x,x)
y = 0.0:0.01:1; x = ones(size(y)); plot(x,y,'--')

repo: Aero-matlab


text:
function [oHDev] = calculateHDEV(tau,sPeriod,readings)
%Overlapping Hadamard (ADEV) function that uses phase error or time error values
%input argument is readings already grouped as tau values. tau is the
%desired gate or tau time needed for calculating overlap.
%sPeriod is used to determine the overlap

N = numel(readings); %get the reading count
n = tau / sPeriod; %calculate averaging factor
n = floor(n); %make sure this is an integer
const = 1/(6*(N - (3*n))*tau^2); %calculate the const mult 1/(2*(N - (2*n))*tau^2)
%sum from i=1 to N-(3*n) (Xi+3m - 3Xi+2m + 3Xi+m - Xi)^2
sum = 0; %variable to store summation calculation

%loop for performing summation
for i = 1:(N-(3*n))
   sum = sum + (readings(i+(3*n)) - (3*readings(i+(2*n))) + (3*readings(i+n)) - readings(i))^2;
end

oHDev = sqrt(const*sum); %square root for Allan Dev


text:
function [u, flag, normres, Q] = gmres(varargin)
%GMRES   Iterative solution of chebfun operator equations.
%   U = GMRES(A, F) attempts to solve the operator equation A(U) = F, where F
%   and U are CHEBFUNs and A is a function handle defining a linear operator
%   on CHEBFUN.
%
%   U = GMRES(A, F, RESTART) chooses a restart parameter. Use Inf or [] for no
%   restarts. The default is Inf.
%
%   U = GMRES(A, F, RESTART, TOL) tries to find a solution for which the
%   residual norm is less than TOL relative to the norm of F. The default
%   value is 1e-10.
%
%   U = GMRES(A, F, RESTART, TOL, MAXIT) stops after MAXIT total iterations.
%   This defaults to 36.
%
%   U = GMRES(A, F, RESTART, TOL, MAXIT, M1INV, M2INV) applies the left and
%   right preconditioners M1INV and M2INV, defined as functions. They should
%   approximate the inverse of A. Note that this usage of preconditioners
%   differs from that in the built-in GMRES.
%
%   [U, FLAG] = GMRES(A,F,...) also returns a convergence FLAG:
%     0 GMRES converged to the desired tolerance TOL within MAXIT iterations.
%     1 GMRES iterated MAXIT times but did not converge.
%
%   [U, FLAG, NORMRES] = GMRES(A, F, ...) also returns a vector of the relative
%   residual norms for all iterations. Note the output ordering is not the same
%   as for built-in GMRES. This calling sequence will also print out updates on
%   the progress of the iteration.
%
% Example:
%   % To solve a simple Volterra integral equation:
%   f = chebfun('exp(-4*x.^2)', [-1 1]);
%   A = @(u) cumsum(u) + 20*u;
%   u = gmres(A, f, Inf, 1e-14);
%
% See also GMRES, CHEBOP/MLDIVIDE.

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

% Parse inputs and supply defaults:
defaults = {[], [], Inf, 1e-10, 36, [], []};
idx = nargin+1:length(defaults);
args = [varargin, defaults(idx)];
[L, f, m, tol, maxiter, M1inv, M2inv] = deal( args{:} );
showtrace = (nargout > 2);

if ( m == 0 )
    % No restarts:
    m = Inf;
end
% Avoid warning about Inf in for loop:
m = min(m, maxiter + 1);

% Initialise:
u = chebfun(0, domain(f));
normb = norm(f);
r = f;
if ( ~isempty(M1inv) )
    r = M1inv(r);
end
normres(1) = norm(r)/normb;

j = 1;                                             % total iterations
while ( (normres(j) > tol) && (j < maxiter) )      % outer iterations

    H = [];
    QTb = norm(r);
    Q = r/QTb ;                                    % Krylov basis

    for n = 1:m                                    % inner iterations

        % Next Krylov vector, with preconditioners.
        q = Q(:,n);
        if ( ~isempty(M2inv) )
            q = M2inv(q);
        end
        v = L(q);

        if ( ~isempty(M1inv) )
            v = M1inv(v);
        end

        % Modified Gram-Schmidt iteration.
        for k = 1:n
            H(k,n) = Q(:,k)' * v ;                   %#ok<AGROW>
            v = v - H(k,n)*Q(:,k);
        end
        H(n+1,n) = norm(v);                          %#ok<AGROW>

        % Use QR factorization to find the residual norm.
        % TODO: This could be made more efficient (worthwhile?).
        QTb(n+1,1) = 0;                              % by orthogonality
        j = j + 1;
        [P, R] = qr(H);
        normres(j) = abs( P(1,n+1)*QTb(1) ) / normb; %#ok<AGROW>

        % Done?
        if ( normres(j) < tol )
            showtrace = false;
            flag = 0;
            break
        end

        % Trace iteration--for debugging, not an official option.
        if ( showtrace && (rem(j - 1, 5) == 0) )
            fprintf('Iteration %2i: relative residual norm = %.3e\n', ...
                j - 1, normres(j));
        end

        % Give up?
        if ( j == maxiter + 1 )
            warning('CHEBFUN:CHEBFUN:gmres:maxiter', ...
                'Maximum number of iterations reached.')
            showtrace = false;
            break
        end

        % New basis vector:
        Q(:,n+1) = v / H(n+1,n);

        % Reorthogonalize (for research only--not an official option).
        if ( rem(n, Inf) == 0 )
            [Q, ignored] = qr(Q, 0);
        end

    end   % end inner iterations

    y = R(1:n,1:n)\(P(:,1:n)'*QTb);                % least squares soln
    u0 = Q(:,1:n)*y;                               % new part of solution
    if ( ~isempty(M2inv) )
        u0 = M2inv(u0);
    end
    u = u + u0;                                    % solution
    v = L(u0);
    if ( ~isempty(M1inv) )
        v = M1inv(v);
    end
    r = r - v;                                     % new residual
    if ( showtrace )
        fprintf('  (restart)\n')
    end

end   % end outer iterations

if ( j >= maxiter )
    flag = 1;
end
if ( nargout < 2 )
    fprintf('\n  Final relative residual: %.3e\n\n', normres(end))
end

end   % end main function

author: chebfun