text:
function gX = mlpardKernDiagGradX(kern, X)

% MLPARDKERNDIAGGRADX Gradient of MLPARD kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the automatic relevance determination multi-layer perceptron kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : mlpardKernParamInit, kernDiagGradX, mlpardkernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

gX = zeros(size(X));
for i = 1:size(X, 1);
  gX(i, :) = mlpardKernDiagGradXpoint(kern, X(i, :));
end

function gX = mlpardKernDiagGradXpoint(kern, x)

% MLPARDKERNDIAGGRADXPOINT Diagonal gradient with respect to one point of x.

innerProd = x*sparse(diag(kern.inputScales))*x';
numer = innerProd*kern.weightVariance + kern.biasVariance;
denom = numer + 1;
arg = numer./denom;
gX = zeros(size(x));
twooverpi = 2/pi;
for j = 1:size(x, 2)
  gX(:, j)=1./denom...
           - numer./denom.^2;
  gX(:, j) = twooverpi*2*kern.inputScales(j)*x(:, j)*kern.weightVariance*kern.variance*gX(:, j)./sqrt(1-arg.*arg);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [sz, maxSize] = h5size(filepath, dataset, dim)
%H5SIZE Returns the size of the specified dataset.
% Usage:
%   [size, maxSize] = h5size(filepath, dataset)
%   [size, maxSize] = h5size(filepath, dataset, dim)

if nargin < 3; dim = []; end

if nargout > 1
    info = h5info(filepath, dataset);

    sz = info.Dataspace.Size;
    maxSize = info.Dataspace.MaxSize;

    if ~isempty(dim)
        sz = sz(dim);
        maxSize = maxSize(dim);
    end
else
    sz = size(h5file(filepath, dataset));
    if ~isempty(dim)
        sz = sz(dim);
    end
end
end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function [rows, cols, vals] = rcvize_matrix(M, i, j)
%function [rows, cols, vals] = rcvize_matrix(M, i, j)
%
% Given a DxD matrix M that we wish to insert into the (i,j)th DxD
% block of a sparse matrix, this function computes and returns the
% corresponding {row, col, val} vectors describing this block

% Copyright (C) 2016 by David M. Rosen

D = size(M, 1);

vals = reshape(M, [1, D^2]);  %Vectorize M by concatenating its columns

rows = repmat( [D*(i-1) + 1 : D*(i-1) + D], [1,D]);
cols = kron( [D*(j-1) + 1 : D*(j-1) + D], ones(1, D));

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function y = cvx_id( x )
y = cvx_id( struct2cell( x ) );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [y] = spm_fx_lz(x,u,P)
% flow for Lorenz attractor
% FORMAT [y] = spm_fx_lz(x,u,P)
% x - state
% u - input
% P - parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_lz.m 1961 2008-07-26 09:38:46Z karl $

% flow for Lorenz attractor
%--------------------------------------------------------------------------
try
    P(3) = P(3)*(1 + u);
end
J    = [-P(1) P(1) 0; (P(3) - x(3)) -1 -x(1); x(2) x(1) P(2)];
y    = J*x;

author: spm


text:
% ========================================================================
% Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion, TIP,2020
% algorithm Version 1.0
% Copyright(c) 2020, Hui Li, Kede Ma, Yongwei Yong and Lei Zhang
% All Rights Reserved.
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereby
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
% This is a dynamic scene implementation of "Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion"
% Please refer to the following paper:
% H. Li et al., "Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion, 2020" In press
% IEEE Transactions on Image Processing
% Please kindly report any suggestions or corrections to xiaohui102788@126.com
%----------------------------------------------------------------------

clear ;
close all;
addpath(genpath(pwd));

    imgSeqColor= loadImg('ArchSequence'); % [0,1]
    imgSeqColor = downSample(imgSeqColor, 1024);

tic
%% motion detection
r=10;
    imgSeqColor = detection_color(imgSeqColor,r);

%% fusion as static
r1=4;
[ D1,i_mean1,aa1,N1] = scale_fine(imgSeqColor,r1);

%% the intermediate  scale
[w,h,~,~]=size(imgSeqColor);
nlev = floor(log(min(w,h)) / log(2))-5;

D2 = cell(nlev,1);
aa2= cell(nlev,1);
N2= cell(nlev,1);

r2=4;
for ii=1:nlev
    [ D2{ii},i_mean2,aa2{ii},N2{ii}] = scale_interm(i_mean1,r2);
    i_mean1=i_mean2;
end

%% the coarsest  scale
r3=4;
[fI3,i_mean3,aa3,N3] = scale_coarse(i_mean2,r3);

%% reconstruct
%% Intermediate layers
for ii=nlev:-1:1
    temp=aa2{ii};
    fI=zeros(size(temp));
    fI(1:2:size(temp,1),1:2:size(temp,2))=fI3;
    B2=boxfilter(fI, r2)./ N2{ii}+D2{ii};

    fI3=B2;
end
%% finest layers
fI=zeros(size(aa1));
fI(1:2:size(aa1,1),1:2:size(aa1,2))=B2;
B1=boxfilter(fI, r1)./ N1;
C_out=repmat(B1,[1 1 3])+D1;
toc

figure,imshow(C_out)

author: thfylsty


text:
function logistic_cdf_values_test ( )

%*****************************************************************************80
%
%% LOGISTIC_CDF_VALUES_TEST demonstrates the use of LOGISTIC_CDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LOGISTIC_CDF_VALUES_TEST:\n' );
  fprintf ( 1, '  LOGISTIC_CDF_VALUES stores values of\n' );
  fprintf ( 1, '  the Logistic Cumulative Density Function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Mu        Beta          X                  CDF(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, mu, beta, x, fx ] = logistic_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %12f  %24.16f\n', mu, beta, x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% SCRIPT TEST FOR THE KINEMATIC PROBLEM FOR THE SAWYER ROBOT
% LOAD THE SAWYER ROBOT. SIMILAR RESULTS CAN BE ACHIEVED WITH OTHER
% REDUNDANT ROBOTS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
close all
fprintf('\nSimple test: try to reach T with the Sawyer robot')
robot = load_robot('practicals', 'SAWYER');
adjust_view(robot)

T1 = [1   0  0  0.5;
     0   1  0  0.5;
     0   0  1  0.5;
     0   0  0   1];

T2 =[0.0000    0.0000    1.0000    1.0148;
   -0.0000   -1.0000    0.0000    0.1603;
    1.0000   -0.0000   -0.0000    0.3170;
         0         0         0    1.0000];

T3 =[0.3164    0.0369    0.9479    0.9796;
   -0.2628   -0.9567    0.1250    0.2655;
    0.9115   -0.2886   -0.2931    0.1686;
         0         0         0    1.0000];

T = T1;

% try this initial seed for the inverse kinematicdifferent seeds
% do allow to obtain different solutions
q = [0.2 -0.2 -0.4 -0.4 -0.2 -0.2 -0.2]';
% q=[0.0 0.0 0.0 0.0 0.0 0.0 0.0]';

drawrobot3d(robot, q)
qinv = inversekinematic(robot, T, q)
T_reach = directkinematic(robot, qinv)

'diff T-Treach'
T-T_reach

repo: ARTE


text:
function B = xrepmat(A,varargin)
%XREPMAT Replicate and tile an array.
%
%   XREPMAT(A,M,N) and XREPMAT(A,[M N]) replicates and tiles the matrix
%   A to produce a M-by-N block matrix.
%
%   XREPMAT(A,(M,N,P,...) and XREPMAT(A,[M N P ...]) tiles the array A
%   to produce a M-by-N-by-P-by-... block array.  A can be N-D.
%
%   XREPMAT(A,M,N) when A is a scalar is commonly used to produce an
%   M-by-N matrix filled with A's value.  This can be much faster than
%   A*ONES(M,N) when M and/or N are large.
%
%   Example:
%       xrepmat(magic(2),2,3)
%       xrepmat(NaN,2,3)
%
%   See also REPMAT, MESHGRID, NDGRID.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:50:27 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

if nargin < 2
   error( 'Not enough input arguments.' );
elseif nargin == 2
   if length(varargin{1}) == 1          % XREPMAT(A,M)
      Rsiz = [ varargin{1} varargin{1} ];
   else                                 % XREPMAT(A,[M N P ...])
      Rsiz = varargin{1};
   end
else                                    % XREPMAT(A,[M N P ...])
   Rsiz = [ varargin{:} ];
end

if length(A) == 1
   nelems = prod(Rsiz);
   if nelems > 0
      % Since B doesn't exist, the first statement creates a B with the
      % right size and type.  Then use scalar expansion to fill the
      % array.  Finally reshape to the specified size.
      B(nelems) = A;
      B(:) = A;
      B = reshape(B,Rsiz);
   else
      B = A(ones(Rsiz));
   end
else
   Asiz = size(A);
   Adim = length(Asiz);
   Rdim = length(Rsiz);
   if ( Adim == 2 ) & ( Rdim == 2 )
      mind = (1:Asiz(1))';
      nind = (1:Asiz(2))';
      mind = mind(:,ones(1,Rsiz(1)));
      nind = nind(:,ones(1,Rsiz(2)));
      B = A(mind,nind);
   else
      Bdim = max(Adim,Rdim);
      Asiz = [ Asiz ones(1,Bdim-Adim) ];
      Rsiz = [ Rsiz ones(1,Bdim-Rdim) ];

      % This might be `clever', but this is usually slower than the
      % for-loop below.
%      v = reshape( reshape( 1:2*Bdim, Bdim, 2 )', 1, 2*Bdim );
%      A = A(:);
%      B = reshape( A(:,ones(1,prod(Rsiz))), [ Asiz Rsiz ] );
%      B = reshape( permute( B, v ), Asiz.*Rsiz );

      for i = 1:Bdim
         ind = (1:Asiz(i))';
         subs{i} = ind(:,ones(1,Rsiz(i)));
      end
      B = A(subs{:});

   end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function wfa = sst2wfa(sst,w)
%
%SST2WFA: Convert events from SST form to WFA form.
%      WFA - WaveForm Array form, events separated in nx1 waveform array
%          - (or) into 1xl cell array, each cell containing a nx1 waveform
%          - array of events from a separate station/channel.
%      NAN - NaN form, (1xl waveform) for l station/channels with events in
%            one time series w/ non-events = NaN.
%      SST - Start/Stop Times form, nx2 array of matlab start/stop times
%            (or) 1xl cell array, each cell containing nx2 start/stop
%            times from a separate station/channel.
%
%USAGE: wfa = sst2wfa(sst)
%
%INPUTS: w    - Waveform object that contains events
%        sst  - Events (WFA form)
%
%OUTPUTS: wfa - Events (WFA form)
%
% See also NAN2SST, NAN2WFA, SSD2SST, SST2NAN, SST2SSD, SST2VAL,
%          WFA2NAN, WFA2SST
%
% Author: Dane Ketner, Alaska Volcano Observatory
% $Date$
% $Revision$

if iscell(sst)
   if numel(sst) == numel(w)
      for l=1:numel(sst)
         for n=1:size(sst{l},1)
            wfa{l}(n) = extract(w,'time',sst{l}(n,1),sst{l}(n,2));
         end
      end
   else
      error(['sst2nan:ArgumentDimensions - Number of elements in',...
             ' waveform input ''w'' and cell input ''sst'' must match']);
   end
else
   for n=1:size(sst,1)
      wfa(n) = extract(w,'time',sst(n,1),sst(n,2));
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function [ fea, out ] = ex_linearelasticity6( varargin )
%EX_LINEARELASTICITY6 NAFEMS LE6 skew plate benchmark
%
%   [ FEA, OUT ] = EX_LINEARELASTICITY6( VARARGIN ) Skew plate under
%   normal pressure (NAFEMS LE6 Benchmark).
%
%   Reference:
%
%   [1] National Agency for Finite Element Methods and Standards. The
%   Standard NAFEMS Benchmarks. Rev. 3. United Kingdom: NAFEMS,
%   October 1990.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       hmax        string {0.05}          Grid size
%       sfun        string {sflag2}        Shape function
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'hmax',    0.05;
            'sfun',     'sflag2';
            'iplot',    1;
            'tol',      0.1;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

E    = 210e9;
nu   = 0.3;
rho  = 7800;
t    = 0.01;
load = -0.7e3;

% Grid definition.
fea.sdim = {'x', 'y', 'z'};
a = cos(pi*30/180);
p = [0 1 1+a a;0 0 0.5 0.5]';
geom_2d.objects{1} = gobj_polygon(p);
grid_2d = gridgen( geom_2d, 'hmax', opt.hmax, 'fid', fid );
fea.grid = gridextrude( grid_2d, max(3,ceil(t/opt.hmax)), t );

% Equations and problem definition.
fea = addphys( fea, @linearelasticity );
fea.phys.el.eqn.coef{1,end} = { nu   };
fea.phys.el.eqn.coef{2,end} = { E    };
fea.phys.el.eqn.coef{3,end} = { rho  };
fea.phys.el.sfun            = { opt.sfun, opt.sfun, opt.sfun };

% Set/constrain w = 0 on sides (boundaries 1-4).
n_bdr = 6;
bctype = mat2cell( zeros(3,n_bdr), [1 1 1], ones(1,n_bdr) );
[bctype{3,[1:4]}] = deal(1);
fea.phys.el.bdr.coef{1,5} = bctype;

% Set/constrain u, v = 0 on edge 2, x = y = 0.
edg(1).index = 2;
edg(1).type = 'constraint';
edg(1).dvar = 1;
edg(1).expr = 0;

edg(2).index = 2;
edg(2).type = 'constraint';
edg(2).dvar = 2;
edg(2).expr = 0;

% Set/constrain v = 0 on edge 1, x = 0, y = 1.
edg(3).index = 1;
edg(3).type = 'constraint';
edg(3).dvar = 2;
edg(3).expr = 0;
fea.edg = edg;

% Apply vertical load to top boundary.
bccoef = mat2cell( zeros(3,n_bdr), [1 1 1], ones(1,n_bdr) );
bccoef{3,6} = load;
fea.phys.el.bdr.coef{1,end} = bccoef;

% Solve problem.
fea = parsephys( fea );
fea = parseprob( fea );

fea.sol.u = solvestat( fea, 'fid', fid );

% Postprocessing.
if( opt.iplot>0 )
  postplot( fea, 'surfexpr', 'sqrt(u^2+v^2+w^2)', ...
            'deformexpr', {'u', 'v', 'w'} )
end

% Error checking.
out = [];
p_E = [0.933012701892219; 0.25; 0];
ps1_E = evalexpr( fea.phys.el.eqn.vars{12,2}, p_E, fea );
ps2_E = evalexpr( fea.phys.el.eqn.vars{13,2}, p_E, fea );
ps3_E = evalexpr( fea.phys.el.eqn.vars{14,2}, p_E, fea );
ps_E_max = max([ps1_E,ps2_E,ps3_E]);
ps_E_ref = 0.802e6;
out.w_E = evalexpr( 'w', p_E, fea );
out.ps_E = [ps1_E, ps2_E, ps3_E];
out.err  = abs(ps_E_max-ps_E_ref)/ps_E_ref;
out.pass = out.err < opt.tol;

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function test_example_sourcemodel_aligned2mni

% MEM 4gb
% WALLTIME 00:10:00

%
%% Create MNI-aligned grids in individual head-space
%
% When combining the source-level estimates of activity in multiple subjects, the activity can first be interpolated on the individuals MRI (using **[ft_sourceinterpolate](https://github.com/fieldtrip/fieldtrip/blob/release/ft_sourceinterpolate.m)**) and then spatially normalized to a template brain (using **[ft_volumenormalise](https://github.com/fieldtrip/fieldtrip/blob/release/ft_volumenormalise.m)**).
%
% It is also possible to define the source-reconstruction grid for each individual subject in such a way that all these grids are already aligned in MNI-space. The combination or statistic of source-level data across subjects can then directly be computed within the source-structure _without the need_ to interpolate and normalize each volume.
%
% The content of this page has moved to another location, and is dealt with **in context** in the tutorial about [the construction of different types of source models](/tutorial/sourcemodel#subject-specific_grids_that_are_equivalent_across_subjects_in_normalized_space).

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function BF_PlotCorrMat(D_corr,rangeHow,makeFigure)
% BF_PlotCorrMat    Visualization of a pairwise similarity matrix
%
% Attempts to determine a set of smaller clusters of objects showing similar
% behavior.
%
%---INPUTS:
% D_corr, a pairwise distance matrix
% rangeHow, the colorbar extent: (i) 'auto' (automatic)
%                                (ii) '-1to1' (range from -1 to 1)
%                                (iii) '0to1' (range from 0 to 1)
%                                (iv) 'balanced' (range from -X to X, for max possible X)
%                                (v) 'positive' (range from 0 to X, for max possible X)
% makeFigure, whether to generate a new figure

% ------------------------------------------------------------------------------
% 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.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs, set defaults:
%-------------------------------------------------------------------------------

if nargin < 2 || isempty(rangeHow)
    rangeHow = 'auto';
end
if nargin < 3
    makeFigure = false;
end

%-------------------------------------------------------------------------------
% Make a matrix
if any(size(D_corr)==1)
    D_corr = squareform(D_corr);
end

% Open a new figure window:
if makeFigure
    f = figure('color','w');
end

% Plot the image:
imagesc(D_corr)

axis square

%-------------------------------------------------------------------------------
% Set color limits and colormap

% If auto setting, can set to balanced if both positive and negative
if strcmp(rangeHow,'auto')
    if any(D_corr(:) > 0) && any(D_corr(:) < 0)
        rangeHow = 'balanced';
    end
end

switch rangeHow
case 'positive'
    maxDev = max(D_corr(:));
    caxis([0,maxDev])
    colormap(BF_GetColorMap('greenblue',9,0))
case 'balanced'
    maxDev = max(abs(D_corr(:)));
    caxis([-maxDev,maxDev])
    colormap([flipud(BF_GetColorMap('blues',9,0));[1,1,1],;BF_GetColorMap('reds',9,0)])
case 'n1to1'
    caxis([-1,1])
    colormap([flipud(BF_GetColorMap('blues',9,0));BF_GetColorMap('reds',9,0)])
case '0to1' % assume 0-1 range (a normalized distance metric)
    caxis([0,1])
    colormap([1,1,1;BF_GetColorMap('reds',9,0)])
case 'auto'
    colormap(flipud(gray(64)))
otherwise
    error('Unknown colormap range option: ''%s''',rangeHow);
end

% ------------------------------------------------------------------------------
% Superimpose black rectangles over NaN values

if any(isnan(D_corr(:)))
    [theNaNs_i,theNaNs_j] = find(isnan(D_corr));
    for i = 1:length(theNaNs_i)
        rectangle('Position',[theNaNs_j(i)-0.5,theNaNs_i(i)-0.5,1,1],'FaceColor','k', ...
                        'EdgeColor','k')
    end
end

% Remove ticks:
set(gca,'TickLength',[0,0])

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% This is a script to run a mesh unfold from the console. It It really an
% example of how such a script should work - it is not particularly
% beautiful. Or big. Or clever.
% ARW 021701

% Basically, we just collect all the date we need to pass to
% unfoldMeshFromGUI. The busyHandle and statusHandle are zero

% Test assignToNearest

saveExtra=1;
showFigures=1;
truePerimDist=1;
busyHandle=0;
statusHandle=0;
scaleFactor=[240/256 240/256 1.2];
perimDist=input('Unfold size: ');
startCoords=input('Start position: ');
fprintf('Browsing for input mesh file:');
[whiteFileName,whitePathName]=uigetfile('*.class','Select class file:');
whiteFileName=[whitePathName,whiteFileName];
[grayFileName,grayPathName]=uigetfile('*.gray','Select gray file:');
grayFileName=[grayPathName,grayFileName];
[flatFileName,flatPathName]=uiputfile('*.mat','Select output file:');
flatFileName=[flatPathName,flatFileName];

dummy=unfoldMeshFromGUI(whiteFileName,grayFileName,flatFileName,startCoords,scaleFactor,perimDist,statusHandle,busyHandle,showFigures,saveExtra,truePerimDist);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function rgb = computeSliceRGB(slice, displayRange, lut)
%COMPUTESLICERGB Convert slice data to renderable slice
%
%   RGB = computeSliceRGB(SLICE)
%
%   RGB = computeSliceRGB(SLICE, DISPLAYRANGE, LUT)
%
%   Example
%     computeSliceRGB
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-11-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% RGB slices are returned directly
if ndims(slice) > 2
    rgb = slice;
    return;
end

% forces computation of extreme values
if isempty(displayRange) || islogical(slice)
    displayRange = [0 max(slice(:))];
end

% converts to uint8, rescaling data between 0 and max value
displayRange = double(displayRange);
extent = displayRange(2) - displayRange(1);
slice = uint8((double(slice) - displayRange(1)) * 255 / extent);

% eventually apply a LUT
if ~isempty(lut)
    if ischar(lut)
        lut = feval(lut, 256);
    end

    lutMax = max(lut(:));
    dim = size(slice);
    rgb = zeros([dim 3], 'uint8');

    % compute each channel
    for c = 1:size(lut, 2)
        res = zeros(size(slice));
        for i = 0:size(lut,1)-1
            res(slice==i) = lut(i+1, c);
        end
        rgb(:,:,c) = uint8(res * 255 / lutMax);
    end

else
    % if no LUT, simply use gray equivalent
    rgb = repmat(slice, [1 1 3]);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [f, lb, ub, A_x, be_u_abs, be_l_abs, s_current_m, error_state] = ...
    prepareOptimizationProblem(VehicleDynamicState, dot_d_numerical, PathPos, v_traj, ...
    ax_diff_traj, ax_traj, ay_traj, ax_lim_mps2, ay_lim_mps2, ...
    d_lim_ub_m, d_lim_lb_m, d_Target_m, dot_d_Target_mps, ...
    vx_lin, kappa_lin, UncertaintyTube, v_terminal_mps, ax_traj_old, ay_traj_old, solution_old, ...
    sys, P_VDC_VirtualController, drag_coefficient, roh_air, vehiclemass_kg, ...
    P_VDC_PositiveAxLimScale, P_VDC_EnableNumLatErrorDer_b, P_VDC_MaxTightening, ...
    P_VDC_MinVelSlipCalc_mps)

% function prepareOptimizationProblem
% Authors:       Martin Euler
%                Salih Guemues
%                Alexander Wischnewski
%
% Description:
%   function used to generate the constraint matrices for the TMPC-Controller
% Inputs/parameters:
%   VehicleDynamicState:    Dynamic state of the racecar
%   PathPos:                Vehicle position on the target path
%   v_traj:                 Target trajectory for the PH (prediction horizon)
%   ax_diff_traj:           Derivative of the target trajectory velocity in the PH
%   ax_traj:                Longitudinal acceleration for the target trajectory in the PH
%   ay_traj:                Lateral acceleration for the target trajectory in the PH
%   ax_lim_mps2             Longitudinal acceleration limits of the target trajectory in the PH
%   ay_lim_mps2             Lateral acceleration limits of the target trajectory in the PH
%   d_lim_ub_m              Size of the driving tube to the right side
%   d_lim_lb_m              Size of the driving tube to the left side
%   d_Target_m              Target trajectory projected into the curvilinear coordinate frame
%   dot_d_Target_mps        Derivative of the target trajectory projected into the c.l. coordinates
%   vx_lin                  Linearization velocity profile
%   kappa_lin               Linearization curvature profile
%   UncertaintyTube         Uncertainy ellipsoids for predicted uncertainty tube
%   v_terminal_set          Terminal set speed
%   ax_traj_old             Previous iteration target trajectory long. acceleration
%   ay_traj_old             Previous iteration target trajectory lat. acceleration
%   solution_old            Previous iteration solution
%   P_VDC_VirtualController Virtual feedback controller used for uncertainty tube calculation
%   drag_coefficient        Vehicle drag coefficient
%   roh_air                 Air density
%   vehiclemass_kg          Vehicle mass
%   P_VDC_PositiveAxLimScale Scaling factor for positive accelerations to consider RWD
%
% Outputs:
%
%   f:                      linear weight matrix for QP
%   lb:                     lower bound vector for inequality constraints
%   ub:                     upper bound vector for inequality constraints
%   A_x:                    nonzero-value-vector of inequality constraint coefficient matrix A_ineq
%   be_u_abs:               upper terminal set bounds in absolute coordinates
%   be_l_abs:               lower terminal set bounds in absolute coordinates
%   s_current:              Vehicle position in raceline coordinates
%   errorState:             current errorState of the system

%% ------- define optimization problem variables ----------------------- %%
% get UNSCALED (in terms of preconditioning) constraint matrices/vectors
A_ineq = zeros(sys.osqp_m, sys.osqp_n);
A_ineq(sys.A_i_lin) = sys.A_x_par;
lb = sys.l_par;
ub = sys.u_par;

%% ------------------------- calculate current error state ------------------------------- %%
% set s_0 to current distance on racetrack
s_current_m = PathPos.s_m;
% calculate assumed velocity. limit to lower value to prevent linearization issues.
v_mps = max(VehicleDynamicState.v_mps, 0.5*P_VDC_MinVelSlipCalc_mps);
% velocity difference: delta_vx = v_path - v_car;
delta_vx = v_traj(1) - v_mps;
% lateral error delta_y as output from local_path_matching with added initial offset
d = PathPos.d_m;
% Delta between vehicle and path heading by interpolation
d_psi = PathPos.psi_rad + VehicleDynamicState.beta_rad;
% lateral error dynamics (either numerical or analytical)
if(P_VDC_EnableNumLatErrorDer_b)
    d_dot = dot_d_numerical;
else
    d_dot = v_mps * sin(d_psi);
end
% current state vector for controller (state vector of error diff. equat.)
error_state = [delta_vx; d; d_dot];
% calculate system response to initial state
x0_resp = sys.Ax0_MPC*error_state;

%% ----------- calc linear weight matrix for QP-Solver ------------%
f = sys.osqp_qpar + 2*(x0_resp'*sys.f_x0)' - 2*(d_Target_m'*sys.f_d_m)' ...
    - 2*(dot_d_Target_mps'*sys.f_dot_d_mps)' ...
    - 2*([ax_traj_old(1); ax_traj]'*sys.f_Dax)' ...
    + 2*([ay_traj_old(1); ay_traj]'*sys.f_Day)' ...
    + 2*(solution_old(1)*sys.f_D_deltaax)' + 2*(solution_old(2)*sys.f_D_deltaay)';
% its minus the ax term since the definitions of delta ax and delta ay are changed

%% - calculation of input and state constraints within prediction horizon - %%
for j = 1:1:sys.N_hor
    % get uncertainty matrix for the inputs
    M_states_current = UncertaintyTube(:, (j-1)*3+1:(j-1)*3+3);
    M_inputs_current = P_VDC_VirtualController*M_states_current*P_VDC_VirtualController';
    % calculate tightenings for all state variables and inputs
    tight_delta_v_mps = sqrt(M_states_current(1, 1));
    tight_d_m = sqrt(M_states_current(2, 2));
    tight_dot_d_mps = sqrt(M_states_current(3, 3));
    % the factors here are swapped with respect to the inputs following this scheme
    % such that only upper bounds are given as the ax_lim has to be modified for the upper bounds
    % ay_lim*ax_util + P_VDC_PositiveAxLimScale*ax_lim*ay_util <= P_VDC_PositiveAxLimScale*ax_lim*ay_lim
    tight_tire1 = sqrt([ay_lim_mps2(j), P_VDC_PositiveAxLimScale*ax_lim_mps2(j)]*M_inputs_current*[ay_lim_mps2(j); P_VDC_PositiveAxLimScale*ax_lim_mps2(j)]);
    % -ay_lim*ax_util + -ax_lim*ay_util <= ax_lim*ay_lim
    tight_tire2 = sqrt([-ay_lim_mps2(j), -ax_lim_mps2(j)]*M_inputs_current*[-ay_lim_mps2(j); -ax_lim_mps2(j)]);
    % -ay_lim*ax_util + ax_lim*ay_util <= ax_lim*ay_lim
    tight_tire3 = sqrt([-ay_lim_mps2(j), ax_lim_mps2(j)]*M_inputs_current*[-ay_lim_mps2(j); ax_lim_mps2(j)]);
    % ay_lim*ax_util + -P_VDC_PositiveAxLimScale*ax_lim*ay_util <= P_VDC_PositiveAxLimScale*ax_lim*ay_lim
    tight_tire4 = sqrt([ay_lim_mps2(j), -P_VDC_PositiveAxLimScale*ax_lim_mps2(j)]*M_inputs_current*[ay_lim_mps2(j); -P_VDC_PositiveAxLimScale*ax_lim_mps2(j)]);
    % calculate normalization factor for acceleration constraints
    norm_acc_con_upper = 1/(P_VDC_PositiveAxLimScale*ax_lim_mps2(j)*ay_lim_mps2(j));
    norm_acc_con_lower = 1/(ax_lim_mps2(j)*ay_lim_mps2(j));

    % limit tigthening to given factor
    tight_d_m = min(d_lim_ub_m(j)*P_VDC_MaxTightening, tight_d_m);
    % get tightened admissible lateral error for current step
    ub(sys.n_constr*j) = d_lim_ub_m(j) - tight_d_m - x0_resp(2+(j-1)*sys.n_sys);
    lb(sys.n_constr*j) = d_lim_lb_m(j) + tight_d_m - x0_resp(2+(j-1)*sys.n_sys);

    % prepare linearization around d = 0, dot_d = 0 and vx = vx_pred for longitudinal acceleration
    ax_grad_d = ax_diff_traj(j)*kappa_lin(j)*vx_lin(j);
    ax_grad_v = ax_diff_traj(j) + roh_air*drag_coefficient*vx_lin(j)/vehiclemass_kg;
    % prepare full gradient for easier specification of constraints
    ax_grad_states = ay_lim_mps2(j)*[ax_grad_v, ax_grad_d, 0];
    ax_grad_inputs = ay_lim_mps2(j)*[-1, 0]; % this is minus one since delta_ax is different sign
    ax_op = ax_diff_traj(j)*(2*vx_lin(j) - v_traj(j)) + 0.5*roh_air*drag_coefficient*vx_lin(j)^2/vehiclemass_kg;

    % prepare linearization around d = 0, dot_d = 0 and vx = vx_pred for lateral acceleration
    ay_grad_d = kappa_lin(j)^2*vx_lin(j)^2;
    ay_grad_dot_d = -ax_diff_traj(j);
    ay_grad_v = -2*kappa_lin(j)*vx_lin(j);
    % prepare full gradient for easier specification of constraints
    ay_grad_states_upper = P_VDC_PositiveAxLimScale*ax_lim_mps2(j)*[ay_grad_v, ay_grad_d, ay_grad_dot_d];
    ay_grad_inputs_upper = P_VDC_PositiveAxLimScale*ax_lim_mps2(j)*[0, 1];
    ay_grad_states_lower = ax_lim_mps2(j)*[ay_grad_v, ay_grad_d, ay_grad_dot_d];
    ay_grad_inputs_lower = ax_lim_mps2(j)*[0, 1];
    ay_op = kappa_lin(j)*vx_lin(j)^2 + 2*kappa_lin(j)*vx_lin(j)*(v_traj(j) - vx_lin(j));

    % update gradients for states and inputs in constraint matrix
    % ay_lim*ax_util + P_VDC_PositiveAxLimScale*ax_lim*ay_util <= P_VDC_PositiveAxLimScale*ax_lim*ay_lim
    A_ineq(sys.n_constr*(j-1)+1, 1:sys.N_hor*sys.m_sys) = ...
        (ax_grad_states+ay_grad_states_upper)*sys.ABK_MPC(1+(j-1)*sys.n_sys:j*sys.n_sys, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+1, sys.m_sys*(j-1)+1:sys.m_sys*j) = ...
        (A_ineq(sys.n_constr*(j-1)+1, sys.m_sys*(j-1)+1:sys.m_sys*j) + ax_grad_inputs + ay_grad_inputs_upper);
    % -ay_lim*ax_util + -ax_lim*ay_util <= ax_lim*ay_lim
    A_ineq(sys.n_constr*(j-1)+2, 1:sys.N_hor*sys.m_sys) = ...
        -(ax_grad_states+ay_grad_states_lower)*sys.ABK_MPC(1+(j-1)*sys.n_sys:j*sys.n_sys, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+2, sys.m_sys*(j-1)+1:sys.m_sys*j) = ...
        (A_ineq(sys.n_constr*(j-1)+2, sys.m_sys*(j-1)+1:sys.m_sys*j) - ax_grad_inputs - ay_grad_inputs_lower);
    % -ay_lim*ax_util + ax_lim*ay_util <= ax_lim*ay_lim
    A_ineq(sys.n_constr*(j-1)+3, 1:sys.N_hor*sys.m_sys) = ...
        (-ax_grad_states+ay_grad_states_lower)*sys.ABK_MPC(1+(j-1)*sys.n_sys:j*sys.n_sys, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+3, sys.m_sys*(j-1)+1:sys.m_sys*j) = ...
        (A_ineq(sys.n_constr*(j-1)+3, sys.m_sys*(j-1)+1:sys.m_sys*j) - ax_grad_inputs + ay_grad_inputs_lower);
    % ay_lim*ax_util + -P_VDC_PositiveAxLimScale*ax_lim*ay_util <= P_VDC_PositiveAxLimScale*ax_lim*ay_lim
    A_ineq(sys.n_constr*(j-1)+4, 1:sys.N_hor*sys.m_sys) = ...
        -(-ax_grad_states+ay_grad_states_upper)*sys.ABK_MPC(1+(j-1)*sys.n_sys:j*sys.n_sys, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+4, sys.m_sys*(j-1)+1:sys.m_sys*j) = ...
        (A_ineq(sys.n_constr*(j-1)+4, sys.m_sys*(j-1)+1:sys.m_sys*j) + ax_grad_inputs - ay_grad_inputs_upper);

    % scale matrices appropriately
    A_ineq(sys.n_constr*(j-1)+1, 1:sys.N_hor*sys.m_sys) = ...
        norm_acc_con_upper*A_ineq(sys.n_constr*(j-1)+1, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+2, 1:sys.N_hor*sys.m_sys) = ...
        norm_acc_con_lower*A_ineq(sys.n_constr*(j-1)+2, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+3, 1:sys.N_hor*sys.m_sys) = ...
        norm_acc_con_lower*A_ineq(sys.n_constr*(j-1)+3, 1:sys.N_hor*sys.m_sys);
    A_ineq(sys.n_constr*(j-1)+4, 1:sys.N_hor*sys.m_sys) = ...
        norm_acc_con_upper*A_ineq(sys.n_constr*(j-1)+4, 1:sys.N_hor*sys.m_sys);

    % limit tire tightening to resonable values
    tight_tire1 = min(P_VDC_PositiveAxLimScale*ax_lim_mps2(j)*ay_lim_mps2(j)*P_VDC_MaxTightening, ...
                        tight_tire1);
    tight_tire2 = min(ax_lim_mps2(j)*ay_lim_mps2(j)*P_VDC_MaxTightening, ...
                        tight_tire2);
    tight_tire3 = min(ax_lim_mps2(j)*ay_lim_mps2(j)*P_VDC_MaxTightening, ...
                        tight_tire3);
    tight_tire4 = min(P_VDC_PositiveAxLimScale*ax_lim_mps2(j)*ay_lim_mps2(j)*P_VDC_MaxTightening, ...
                        tight_tire4);

    % ay_lim*ax_util + P_VDC_PositiveAxLimScale*ax_lim*ay_util <= P_VDC_PositiveAxLimScale*ax_lim*ay_lim
    ub(sys.n_constr*(j-1)+1) = norm_acc_con_upper*(-ax_op*ay_lim_mps2(j) ...
        - ay_op*P_VDC_PositiveAxLimScale*ax_lim_mps2(j) ...
        - ax_grad_states*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        - ay_grad_states_upper*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        - tight_tire1 + P_VDC_PositiveAxLimScale*ax_lim_mps2(j) * ay_lim_mps2(j));
    % -ay_lim*ax_util + -ax_lim*ay_util <= ax_lim*ay_lim
    ub(sys.n_constr*(j-1)+2) = -norm_acc_con_lower*(-ax_op*ay_lim_mps2(j) ...
        - ay_op*ax_lim_mps2(j) ...
        - ax_grad_states*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        - ay_grad_states_lower*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        + tight_tire2 - ax_lim_mps2(j) * ay_lim_mps2(j));
    % -ay_lim*ax_util + ax_lim*ay_util <= ax_lim*ay_lim
    ub(sys.n_constr*(j-1)+3) = norm_acc_con_lower*(ax_op*ay_lim_mps2(j) ...
        - ay_op*ax_lim_mps2(j) ...
        + ax_grad_states*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        - ay_grad_states_lower*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        - tight_tire3 + ax_lim_mps2(j) * ay_lim_mps2(j));
    % ay_lim*ax_util + -P_VDC_PositiveAxLimScale*ax_lim*ay_util <= P_VDC_PositiveAxLimScale*ax_lim*ay_lim
    ub(sys.n_constr*(j-1)+4) = -norm_acc_con_upper*(ax_op*ay_lim_mps2(j) ...
        - ay_op*P_VDC_PositiveAxLimScale*ax_lim_mps2(j) ...
        + ax_grad_states*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        - ay_grad_states_upper*x0_resp(1+(j-1)*sys.n_sys:j*sys.n_sys) ...
        + tight_tire4 - P_VDC_PositiveAxLimScale*ax_lim_mps2(j) * ay_lim_mps2(j));
end

%% -------------------- terminal state constraints ----------------------- %%
be_u = zeros(3, 1);
be_l = zeros(3, 1);
be_u(1) = inf;
be_l(1) = v_traj(end) - v_terminal_mps;

% lateral error is constrained to be
tight_terminal_d_m = sqrt(UncertaintyTube(2, end-1));
tight_terminal_d_m = min(d_lim_ub_m(end)*P_VDC_MaxTightening, tight_terminal_d_m);
be_u(2) = d_lim_ub_m(end) - tight_terminal_d_m;
be_l(2) = d_lim_lb_m(end) + tight_terminal_d_m;
% calculate lateral erorr speed as average of the last five values as they might be noisy
be_u(3) = mean(dot_d_Target_mps(end-5:end));
be_l(3) = mean(dot_d_Target_mps(end-5:end));
% write terminal constraints to logs
be_u_abs = be_u;
be_l_abs = be_l;
% add offset velocity to make logging meaningful
be_u_abs(1) = v_traj(end) - be_u_abs(1);
be_l_abs(1) = v_traj(end) - be_l_abs(1);
% calculate ub for terminal constraint on states
ub(sys.n_constr*sys.N_hor+1:sys.n_constr*sys.N_hor+sys.n_sys) = be_u - ...
    x0_resp(sys.N_hor*sys.n_sys+1:(sys.N_hor+1)*sys.n_sys);
% calculate lb for terminal constraint on states
lb(sys.n_constr*sys.N_hor+1:sys.n_constr*sys.N_hor+sys.n_sys) = be_l - ...
    x0_resp(sys.N_hor*sys.n_sys+1:sys.n_sys*(sys.N_hor+1));
% no terminal constraints on longitudinal acceleration
ub(sys.n_constr*sys.N_hor+sys.n_sys+1) = inf;
lb(sys.n_constr*sys.N_hor+sys.n_sys+1) = -inf;
% lateral acceleration is constrained to be close to target
ub(sys.n_constr*sys.N_hor+sys.n_sys+2) = 0;
lb(sys.n_constr*sys.N_hor+sys.n_sys+2) = 0;

%% --------------------- setup slack constraints ------------------------------------- %%
% limit slacks to be strictly positive
lb(sys.N_hor*sys.n_constr+sys.n_sys+sys.m_sys+1:end) = 0;
% upper limit of tire slacks set to certain percentage of maximum available accelerations
ub(sys.N_hor*sys.n_constr+sys.n_sys+sys.m_sys+1:sys.n_slacks:end) = sys.slack_lim_rel;
ub(sys.N_hor*sys.n_constr+sys.n_sys+sys.m_sys+2:sys.n_slacks:end) = sys.slack_lim_rel;
% upper limit of lateral error slacks set to equal amount of allowed deviation
% which allows to go for the double amount if required
ub(sys.N_hor*sys.n_constr+sys.n_sys+sys.m_sys+3:sys.n_slacks:end) = 0.5*(d_lim_ub_m(1)-d_lim_lb_m(1));
ub(sys.N_hor*sys.n_constr+sys.n_sys+sys.m_sys+4:sys.n_slacks:end) = 0.5*(d_lim_ub_m(1)-d_lim_lb_m(1));

%% ----------- transform A_ineq and P to CSC compatible format using sparsity pattern ----------- %%
A_x = A_ineq(sys.A_i_lin);

end

author: TUMFTM


text:
function [new, fittype]=dfcreatecopy(original);

%   $Revision: 1.1.6.3 $  $Date: 2004/01/24 09:35:23 $
%   Copyright 2003-2004 The MathWorks, Inc.

fittype = original.fittype;

new = copyfit(original);
new = java(new);

author: zouchuhang

repo: LayoutNet


text:
function [X,err,iter] = lrtcR_snn(M,omega,alpha,opts)

% Solve the Noisy Low-Rank Tensor Completion (LRTC) based on Sum of Nuclear Norm (SNN) problem by M-ADMM
%
% min_{X,E} \sum_i \alpha_i*||X_{i(i)}||_* + loss(E),
% s.t. P_Omega(X) + E = M.
% loss(E) = ||E||_1 or 0.5*||E||_F^2
%
% ---------------------------------------------
% Input:
%       M       -    d1*d2*...dk tensor
%       omega   -    index of the observed entries
%       alpha   -    k*1 vector, parameters
%       opts    -    Structure value in Matlab. The fields are
%           opts.loss       -   'l1' (default): loss(E) = ||E||_1
%                               'l2': loss(E) = 0.5*||E||_F^2
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       X       -    d1*d2*...*dk tensor
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 24/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;
loss = 'l1';

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'loss');        loss = opts.loss;            end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

dim = size(M);
k = length(dim);

omegac = setdiff(1:prod(dim),omega);

X = zeros(dim);
Y = cell(k,1);
Z = Y;
E = X;
Y2 = E;
for i = 1 : k
    Y{i} = X;
    Z{i} = X;
end

iter = 0;
for iter = 1 : max_iter
    Xk = X;
    Ek = E;
    Zk = Z;
    % first super block {Z_i,E}
    sumtemp = zeros(dim);
    for i = 1 : k
        Z{i} = Fold(prox_nuclear(Unfold(X+Y{i}/mu,dim,i), alpha(i)/mu),dim,i);
        sumtemp = sumtemp + Z{i} - Y{i}/mu;
    end
    if strcmp(loss,'l1')
        E = prox_l1(-X+M-Y2/mu,1/mu);
    elseif strcmp(loss,'l2')
        E = (-X+M-Y2/mu)*(mu/(1+mu));
    else
        error('not supported loss function');
    end
    % second super block {X}
    X(omega) = (sumtemp(omega)-Y2(omega)/mu-E(omega)+M(omega))/(k+1);
    X(omegac) = sumtemp(omegac)/k;

    chg = max([max(abs(Xk(:)-X(:))), max(abs(Ek(:)-E(:))) ]);
    err = 0;
    for i = 1 : k
        dY = X-Z{i};
        err = err+norm(dY(:))^2;
        Y{i} = Y{i}+mu*dY;
        chg = max([chg,max(abs(dY(:))), max(abs((Zk{i}(:)-Z{i}(:))))]);
    end
    dY = E-M;
    dY(omega) = dY(omega)+X(omega);
    chg = max(chg,max(abs(dY(:))));
    Y2 = Y2 + mu*dY;
    err = sqrt(err+norm(dY(:))^2);

    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', err=' num2str(chg)]);
        end
    end
    if chg < tol
        break;
    end
    mu = min(rho*mu,max_mu);
end

author: canyilu

repo: LibADMM-toolbox


text:
%% sampleCurveEvenly
% Below is a basic demonstration of the features of the |sampleCurveEvenly| function.

%%
clear; close all; clc;

% PLOT SETTINGS
fig_color='w'; fig_colordef='white';
markerSize=15;
lineWidth=2;

%% EXAMPLE USING DEFAULT SETTINGS FOR RESAMPLING A CURVE EVENLY

%Simulating the case of an unevenly sampled loop curve
ns=50;
x=linspace(0,1,ns); %evenly spaced x data
y=sin(x*2*pi); %resulting y data
V=[x(:) y(:)]; %nonlinearity causing uneven point spacing

cPar=[]; %If the control structure is empty these are the settings used:

% cPar.nd=size(V,1); %resample with the same about of points
% cPar.typeOpt='num'; %curve steps are based on the desired number of points
% cPar.interpMethod='pchip'; %use pchip interpolation
% cPar.closeLoopOpt=0; %the curve is not a closed loop

[Vg]=sampleCurveEvenly(V,cPar);

cFigure;
subplot(1,2,1); hold on;
title('Unevenly sampled');
plotV(V,'r.-','MarkerSize',markerSize);
view(2); grid on; axis equal; axis tight;
subplot(1,2,2); hold on;
title('Evenly sampled using default settings');
plotV(Vg,'g.-','MarkerSize',markerSize);
plotV(Vg(1,:),'r.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
plotV(Vg(end,:),'b.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
view(2); grid on; axis equal; axis tight;
drawnow;

%% EXAMPLE USING CUSTOM SETTINGS FOR RESAMPLING A CURVE EVENLY

%Simulating the case of an unevenly sampled loop curve

cPar.nd=size(V,1)*2; %upsample twice
cPar.typeOpt='num'; %curve steps are based on the desired number of points
cPar.interpMethod='linear'; %use linear interpolation
cPar.closeLoopOpt=0; %the curve is not a closed loop

[Vg]=sampleCurveEvenly(V,cPar);

cFigure;
subplot(1,2,1); hold on;
title('Unevenly sampled');
plotV(V,'r.-','MarkerSize',markerSize);
view(2); grid on; axis equal; axis tight;
subplot(1,2,2); hold on;
title('Evenly sampled allong curve');
plotV(Vg,'g.-','MarkerSize',markerSize);
plotV(Vg(1,:),'r.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
plotV(Vg(end,:),'b.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
view(2); grid on; axis equal; axis tight;
drawnow;

%% EXAMPLE USING CUSTOM SETTINGS FOR RESAMPLING A CURVE EVENLY

%Simulating the case of an unevenly sampled loop curve

cPar.nd=0.1; %resmaple based on curve length step size
cPar.typeOpt='dist'; %curve steps are based on the desired number of points
cPar.interpMethod='linear'; %use linear interpolation
cPar.closeLoopOpt=0; %the curve is not a closed loop

[Vg]=sampleCurveEvenly(V,cPar);

cFigure;
subplot(1,2,1); hold on;
title('Unevenly sampled');
plotV(V,'r.-','MarkerSize',markerSize);
view(2); grid on; axis equal; axis tight;
subplot(1,2,2); hold on;
title('Evenly sampled allong curve');
plotV(Vg,'g.-','MarkerSize',markerSize);
plotV(Vg(1,:),'r.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
plotV(Vg(end,:),'b.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
view(2); grid on; axis equal; axis tight;
drawnow;

V(:,2)=V(:,2)*2; %Increasing amplitude increases point spacing and  hence increases number of point used for curve step size resampling

cPar.nd=0.1; %resample based on curve length step size
cPar.typeOpt='dist'; %curve steps are based on the desired number of points
cPar.interpMethod='linear'; %use linear interpolation
cPar.closeLoopOpt=0; %the curve is not a closed loop

[Vg]=sampleCurveEvenly(V,cPar);

cFigure;
subplot(1,2,1); hold on;
title('Unevenly sampled');
plotV(V,'r.-','MarkerSize',markerSize);
view(2); grid on; axis equal; axis tight;
subplot(1,2,2); hold on;
title('Evenly sampled allong curve');
plotV(Vg,'g.-','MarkerSize',markerSize);
plotV(Vg(1,:),'r.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
plotV(Vg(end,:),'b.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
view(2); grid on; axis equal; axis tight;
drawnow;

%% EXAMPLE USING CUSTOM SETTINGS FOR RESAMPLING A CLOSED POLYGON EVENLY

%Simulating the case of an unevenly sampled loop curve
ns=150;
t=sort(linspace(0,2*pi,ns)+pi/10*rand(1,ns));
t=unique(t); %removing double points
t=t(t<2*pi);%Removing 2*pi points since they are the same as the 0 point
r=3+2.*sin(5*t);
[x,y] = pol2cart(t,r);
z=y;
V=[x(:) y(:) z(:)];

cPar.nd=200; %Desired number of pointsupsample twice
cPar.typeOpt='num'; %curve steps are based on the desired number of points
cPar.interpMethod='pchip'; %use pchip interpolation
cPar.closeLoopOpt=1; %the curve is close so the end is considered equal to the start

[Vg]=sampleCurveEvenly(V,cPar);

cFigure;
subplot(1,2,1); hold on;
title('Unevenly sampled');
plot3(V(:,1),V(:,2),V(:,3),'r.-','MarkerSize',markerSize);
drawnow; view(3); grid on; axis equal; axis tight;
subplot(1,2,2); hold on;
title('Evenly sampled allong curve');
plot3(Vg(:,1),Vg(:,2),Vg(:,3),'g.-','MarkerSize',markerSize);
plot3(Vg(1,1),Vg(1,2),Vg(1,3),'r.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
plot3(Vg(end,1),Vg(end,2),Vg(end,3),'b.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
drawnow; view(3); grid on; axis equal; axis tight;

%% EXAMPLE USING CURVE SMOOTHENING

%Adding noise
V=V+0.2.*randn(size(V));

cPar.nd=200;
cPar.typeOpt='num';
cPar.interpMethod=0.5;
cPar.closeLoopOpt=1;

[Vg]=sampleCurveEvenly(V,cPar);

cFigure;
subplot(1,2,1); hold on;
title('Unevenly sampled');
plot3(V(:,1),V(:,2),V(:,3),'r.-','MarkerSize',markerSize);
drawnow; view(3); grid on; axis equal; axis tight;
subplot(1,2,2); hold on;
title('Evenly sampled allong curve and smoothened');
plot3(Vg(:,1),Vg(:,2),Vg(:,3),'g.-','MarkerSize',markerSize,'lineWidth',lineWidth);
plot3(Vg(1,1),Vg(1,2),Vg(1,3),'r.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
plot3(Vg(end,1),Vg(end,2),Vg(end,3),'b.','MarkerSize',2*markerSize,'lineWidth',lineWidth);
drawnow; view(3); grid on; axis equal; axis tight;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% Description errorb(x,y,varargin)
% errorb(Y,E) plots Y and draws an error bar at each element of Y. The
% error bar is a distance of E(i) above and below the curve so that each
% bar is symmetric and 2*E(i) long.
% If Y and E are a matrices, errob groups the bars produced by the elements
% in each row and plots the error bars in their appropriate place above the
% bars.
%
% errorb(X,Y,E) plots Y versus X with
% symmetric error bars 2*E(i) long. X, Y, E must
% be the same size. When they are vectors, each error bar is a distance of E(i) above
% and below the point defined by (X(i),Y(i)).
%
% errorb(X,Y,'Parameter','Value',...) see below
%
%% Optional Parameters
%    horizontal: will plot the error bars horizontally rather than vertically
%    top: plot only the top half of the error bars (or right half for horizontal)
%    barwidth: the width of the little hats on the bars (default scales with the data!)
%              barwidth is a scale factor not an absolute value.
%    linewidth: the width of the lines the bars are made of (default is 2)
%    points: will plot the points as well, in the same colors.
%    color: specify a particular color for all the bars to be (default is black, this can be anything like 'blue' or [.5 .5 .5])
%    multicolor: will plot all the bars a different color (thanks to my linespecer function)
%    fill: this will plot error bounds in shaded color.
%                colormap: in the case that multicolor is specified, one
%                           may also specify a particular colormap to
%                           choose the colors from.
%% Examples
% y=rand(1,5)+1; e=rand(1,5)/4;
% hold off; bar(y,'facecolor',[.8 .8 .8]); hold on;
% errorb(y,e);
%
% defining x and y
% x=linspace(0,2*pi,8); y=sin(x); e=rand(1,8)/4;
% hold off; plot(x,y,'k','linewidth',2); hold on;
% errorb(x,y,e)
%
% group plot:
% values=abs(randn(2,3))+2; errors=abs(randn(2,3)/1.5+.5)/2;
% errorb(values,errors);
% errorb(values,errors,'top');
%
% % motivation for the function
% It is possible to plot nice error bars on top of a bar plot with Matlab's
% built in errorbar function by setting tons of different parameters to be
% various things.
% This function plots what I would consider to be nice error bars as the
% default, with no modifications necessary.
% It also plots, only the error bars, and in black. There are some very
% useful abilities that this function has over the matlab one, see below:
%
%% Acknowledgments
% Thank you to the AP-Lab at Boston University for funding me while I
% developed these functions. Thank you to the AP-Lab, Avi and Eli for help
% with designing and testing them.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Jonathan Lansey, last update Oct 2011,                                  %
%                   questions to Lansey at gmail.com                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function errorb(x,y,varargin)
%% first things first
%save the initial hold state of the figure.
hold_state = ishold;
if ~hold_state
    cla;
end

h=[];

%% If you are plotting errobars on Matlabs grouped bar plot
if size(x,1)>1 && size(x,2)>1 % if you need to do a group plot
% Plot bars
    num2Plot=size(x,2);
    e=y; y=x;
	handles.bars = bar(x, 'edgecolor','k', 'linewidth', 2);
	hold on
	for i = 1:num2Plot
		x =get(get(handles.bars(i),'children'), 'xdata');
		x = mean(x([1 3],:));
%       Now the recursive call!
		errorb(x,y(:,i), e(:,i),'barwidth',1/(num2Plot),varargin{:}); %errorb(x,y(:,i), e(:,i),'barwidth',1/(4*num2Plot),varargin{:});
    end
    if ~hold_state
        hold off;
    end
    return; % no need to see the rest of the function
else
    x=x(:)';
    y=y(:)';
    num2Plot=length(x);
end
%% Check if X and Y were passed or just X

if ~isempty(varargin)
    if ischar(varargin{1})
        justOneInputFlag=1;
        e=y; y=x; x=1:num2Plot;
    else
        justOneInputFlag=0;
        e=varargin{1}(:)';
    end
else % not text arguments, not even separate 'x' argument
    e=y; y=x; x=1:length(e);
    justOneInputFlag=0;
end

hold on; % axis is already cleared if hold was off
%% Check that your vectors are the proper length
if num2Plot~=length(e) || num2Plot~=length(y)
    error('your data must be vectors of all the same length')
end

%% Check that the errors are all positive
signCheck=min(0,min(e(:)));
if signCheck<0
    error('your error values must be greater than zero')
end

%% In case you don't specify color:
color2Plot = [ 0 0 0];
% set all the colors for each plot to be the same one ... hope you like black
for kk=1:num2Plot
    lineStyleOrder{kk}=color2Plot;
end

%% Initialize some things before accepting user parameters
horizontalFlag=0;
topFlag=0;
pointsFlag=0;
barFactor=1;
linewidth=2;
colormapper='jet';
multicolorFlag=0;
fillFlag=0;

%% User entered parameters
%  if there is just one input, then start at 1,
%  but if X and Y were passed then we need to start at 2
k = 1 + 1 - justOneInputFlag; %
%
while k <= length(varargin) && ischar(varargin{k})
    switch (lower(varargin{k}))
      case 'horizontal'
        horizontalFlag=1;
        if justOneInputFlag % need to switch x and y now
            x=y; y=1:num2Plot; % e is the same
        end
      case 'color' %
        color2Plot = varargin{k + 1};
%       set all the colors for each plot to be the same one
        for kk=1:num2Plot
            lineStyleOrder{kk}=color2Plot;
        end
        k = k + 1;
      case 'linewidth'
        linewidth = varargin{k + 1};
        k = k + 1;
      case {'barwidth','width','thickness'}
        barFactor = varargin{k + 1};
%         barWidthFlag=1;
        k = k + 1;
      case 'points'
        pointsFlag=1;
      case 'multicolor'
          multicolorFlag=1;
      case 'colormap' % used only if multicolor
        colormapper = varargin{k+1};
        k = k + 1;
      case 'top'
          topFlag=1;
        case 'fill'
            fillFlag=1;
      otherwise
        warning('Dude, you put in the wrong argument');
    end
    k = k + 1;
end

if ~fillFlag % plot errobars normally.

if multicolorFlag
    lineStyleOrder=linspecer(num2Plot,colormapper);
end

%% Set the bar's width if not set earlier
if num2Plot==1
%   defaultBarFactor=how much of the screen the default bar will take up if
%   there is only one number to work with.
    defaultBarFactor=20;
    p=axis;
    if horizontalFlag
        barWidth=barFactor*(p(4)-p(3))/defaultBarFactor;
    else
        barWidth=barFactor*(p(2)-p(1))/defaultBarFactor;
    end
else % is more than one datum
    if horizontalFlag
        barWidth=barFactor*(y(2)-y(1))/4;
    else
        barWidth=barFactor*(x(2)-x(1))/4;
    end
end

%% Plot the bars
for k=1:num2Plot
    if horizontalFlag
        ex=e(k);
        esy=barWidth/2;
%       the main line
        if ~topFlag || x(k)>=0  %also plot the bottom half.
            h(end+1) = plot([x(k)+ex x(k)],[y(k) y(k)],'color',lineStyleOrder{k},'linewidth',linewidth);
    %       the hat
            h(end+1) = plot([x(k)+ex x(k)+ex],[y(k)+esy y(k)-esy],'color',lineStyleOrder{k},'linewidth',linewidth);
        end
        if ~topFlag || x(k)<0  %also plot the bottom half.
            h(end+1) = plot([x(k) x(k)-ex],[y(k) y(k)],'color',lineStyleOrder{k},'linewidth',linewidth);
            h(end+1) = plot([x(k)-ex x(k)-ex],[y(k)+esy y(k)-esy],'color',lineStyleOrder{k},'linewidth',linewidth);
            %rest?
        end
    else %plot then vertically
        ey=e(k);
        esx=barWidth/2;
%         the main line
        if ~topFlag || y(k)>=0 %also plot the bottom half.
            h(end+1) = plot([x(k) x(k)],[y(k)+ey y(k)],'color',lineStyleOrder{k},'linewidth',linewidth);
    %       the hat
            h(end+1) = plot([x(k)+esx x(k)-esx],[y(k)+ey y(k)+ey],'color',lineStyleOrder{k},'linewidth',linewidth);
        end
        if ~topFlag || y(k)<0 %also plot the bottom half.
            h(end+1) = plot([x(k) x(k)],[y(k) y(k)-ey],'color',lineStyleOrder{k},'linewidth',linewidth);
            h(end+1) = plot([x(k)+esx x(k)-esx],[y(k)-ey y(k)-ey],'color',lineStyleOrder{k},'linewidth',linewidth);
        end
    end
end
%
%% plot the points, very simple

if pointsFlag
    for k=1:num2Plot
        h(end+1) = plot(x(k),y(k),'o','markersize',8,'color',lineStyleOrder{k},'MarkerFaceColor',lineStyleOrder{k});
    end
end

else % plot error bars by filling them in.

%     Plot the mean:
    plot(x,y,'linewidth',2,'color',color2Plot);
%     Make the polygon:

    xPoly = [x  fliplr(x) x(1)];
    yPoly = [y+e fliplr(y-e) y(1)+e(1)];

    h(end+1) = plot(xPoly,yPoly,'color',color2Plot)
    hReg = fill(xPoly,yPoly,color2Plot); % draw region
    set(get(get(hReg,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); % Exclude line from legend

end % fillFlag check over

for hLoop = 1:length(h)
    set(get(get(h(hLoop),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end

drawnow;
% return the hold state of the figure
if ~hold_state
    hold off;
end

end

%%

function lineStyleOrder=linspecer(N,varargin)

% default colormap
colormap hot; A=colormap;

%% interperet varagin
if ~isempty(varargin)>0
    colormap (varargin{1}); A=colormap;
end

%%

if N<=0
    lineStyleOrder={};
    return;
end

values=round(linspace(1,50,N));
for n=1:N
    lineStyleOrder(n) = {A(values(n),:)};
end

end

author: Sable

repo: mcbench-benchmarks


text:
%% kwJet
% Below is a demonstration of the features of the |kwJet| function

%%
clear; close all; clc;

%% Syntax
% |cMap=kwJet(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function OUT = MovCorrCent(X, window, i)
% =======================================================================
% Computes correlation between X_{i} and X_j for all j different
% of i, from windowdowSize to the number of observations.  The
% first windowdowSize rows are NaN.
% =======================================================================
% OUT = MovCorr(X, windowdowSize, i)
% -----------------------------------------------------------------------
% INPUT
%   - X: panel time series T observations x N variables
%   - windowdow: size of the moving windowdow
%   - i: the variable X_i against which correlations should be returned
% -----------------------------------------------------------------------
% OUTPUT
%   - OUT: matrix with moving correlation T observations x N-1
%       variables. The first windowdowSize-1 rows are NaN.
% =======================================================================
% EXAMPLE
% X = rand(100,5);
% OUT = MovCorrCent(X,10,1)
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015

[nobs,nvars] = size(X);
OUT = nan(nobs,nvars-1);

for tt = window+1:nobs-window
    C = corrcoef(X(tt-window:tt+window, :));
    idx = setdiff(1:nvars, i);
    OUT(tt, :) = C(i, idx);
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [x, y, z] = geodetic_to_ecf(lat, lon, alt)
%GEODETIC_TO_ECF Convert geodetic coordinates to ECF
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////
%
%
% Convert geodetic latitude, longitude, and altitude to ECF (Earth
% Centered Fixed) coordinates.
%
% USAGE:
%   pos_ecf = geodetic_to_ecf(pos_lla)
%   [pos_ecf_x, pos_ecf_y, pos_ecf_z] = geodetic_to_ecf(lat, lon, alt)
%
% INPUTS:
%   pos_lla - required : geodetic latitude, longitude, and altitude   [deg, deg, m]
%
% OUTPUTS:
%   pos_ecf - required : ecf x, y, z coordinates                      [m, m, m]
%
% NOTES:
%   Zhu, J. Conversion of Earth-centered, Earth-fixed coordinates to
%   geodetic coordinates. IEEE Transactions on Aerospace and Electronic
%   Systems, 30, 3 (July 1994), 957-962.
%
% VERSION:
%   1.0
%     - Sean Hatch 20070911
%     - initial version
%   1.1
%     - Wade Schwartzkopf 20130708
%     - vectorized and componentwise data handling
%
% TODO:

% define constants
e2 = 6.6943799901377997e-3; % eccentricity squared of Earth (WGS 84 value)
a = 6378137.0;              % semimajor radius of the Earth (WGS 84 value)

% Handle different forms in input arguments
if nargin==3 % Componentwise inputs, separate arguments for lat,lon,alt
    % Nothing to do.  Processing uses this form.
elseif size(lat,1)==3 % Array of 3-element vectors
    alt = lat(3,:);
    lon = lat(2,:);
    lat = lat(1,:);
elseif numel(lat)==3 % Horizontal 3-element vector
    alt = lat(3);
    lon = lat(2);
    lat = lat(1);
else
    error('WGS_84_NORM:INVALID_INPUTS', 'Invalid inputs.');
end

% calculate distance to surface of ellipsoid
R = a ./ sqrt(1.0 - e2 .* sind(lat) .* sind(lat));

% calculate coordinates
x = (R + alt) .* cosd(lat) .* cosd(lon);
y = (R + alt) .* cosd(lat) .* sind(lon);
z = (R + alt - e2 .* R) .* sind(lat);

if nargout < 2 % Matrix/vector form, rather than componentwise form, was requested
    x = [x; y; z];
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [thresh,cntR,sumR,cntP,sumP,miHist] = affinityPRMI(pairs,segs,nthresh)
% function [thresh,cntR,sumR,cntP,sumP,miHist] = affinityPRMI(pairs,segs,nthresh)
%
% Calcualte precision/recall and mutual information data.
%
% INPUT
%	pairs		3xN array of (i,j,wij) triples.
%	segs		Array of segmentations.
%	[nthresh]	Number of points in PR curve.
%
% OUTPUT
%	thresh		Vector of threshold values.
%	cntR,sumR	Ratio gives recall.
%	cntP,sumP	Ratio gives precision.
%	miHist		Joint histogram of (sameseg,Wij) for MI.
%
% For pairs, the i,j indices are 1-based 1D indices into the
% segmentations.
%
% See also calcMI.
%
% David Martin <dmartin@eecs.berkeley.edu>
% January 2003

if nargin<3, nthresh = 100; end
nthresh = max(1,nthresh);

if size(pairs,1)~=3,
  error('pairs must be 3xN');
end
n = size(pairs,2);

if pairs(3,:)<0 | pairs(3,:)>1,
  error('illegal wij value in pairs(3,:); wij must be in [0,1]');
end

nsegs = length(segs);
if nsegs==0,
  error('segs is empty');
end

[height,width] = size(segs{1});
thresh = linspace(1/(nthresh+1),1-1/(nthresh+1),nthresh)';

% For the mutual information, we need the joint distribution of the
% same segment indicator (given by the segmentations) and Wij.  We
% will bin the Wij values, and so get a 2D histogram estimate of the
% joint.  This histogram is also sufficient information to compute
% precision and recall.
miHist = zeros(nthresh,2);
for index = 1:n,
  % groundtruth is same-segment iff all humans give same-segment
  i = pairs(1,index);
  j = pairs(2,index);
  same = 1;
  for s = 1:nsegs,
    same = same & (segs{s}(i)==segs{s}(j));
  end
  % bin the wij value
  wij = pairs(3,index);
  bin = 1+round(wij*(nthresh-1));
  % increment histogram
  miHist(bin,1+same) = miHist(bin,1+same) + 1;
end

% compute precision and recall
cumHist = flipud(cumsum(flipud(miHist)));
cntR = cumHist(:,2);
sumR = zeros(size(cntR))+cumHist(1,2);
cntP = cumHist(:,2);
sumP = sum(cumHist,2);

author: Cloud-CV

repo: object-proposals


text:
classdef Reshape < dagnn.ElementWise
  properties
    shape = {} ;
  end

  properties (Transient)
    inputSizes = {}
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnreshape(inputs{1}, obj.shape) ;
      obj.inputSizes = cellfun(@size, inputs, 'UniformOutput', false) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nnreshape(inputs{1}, obj.shape, derOutputs{1}) ;
      derParams = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes{1} = [0,0,0,0] ;
    end

    function rfs = getReceptiveFields(obj)
    end

    function load(obj, varargin)
      s = dagnn.Layer.argsToStruct(varargin{:}) ;
      load@dagnn.Layer(obj, s) ;
    end

    function obj = Reshape(varargin)
      obj.load(varargin{:}) ;
    end
  end
end

repo: MFT


text:
function [data smpr meta] = read(db,index)
%READ - returns waveform from ADT object.
%[data smpr meta] = read(ADTobj,index).
% See query function for the description of the
% output arguments.
%Exemples:
%  wave  = read(ADTobj);  -  returns wave data in the form of cell array.
%  wave = read(ADTobj,5); - returns wave data only from sentence 5;
%
%See also query, filterdb, read, play.

if ~isempty(db.enteries)
    Odata=struct();
    count = 0;
    if nargin<2 %if index not enterd reads from all sentances
        from = 1;
        to = length(db.enteries);
    else         %else read only form sentence number index
        from = index;
        to = index;
    end

    fprintf('Reading Data... \n');

    for  index=from:to  %reads data according to the database kind

        if db.kind==1 %reads sentences
            count = count + 1;
            [Odata(count).data smpr] = sentence_read(db,index);
            Odata(count).info =getmeta(db,index);
        end

        if db.kind==2 %reads words
            [data smpr] = sentence_read(db,index);
            for ii=1:length(db.enteries(index).word.b)
                if db.enteries(index).word.flag(ii)==1
                    count = count + 1;
                    b = db.enteries(index).word.b(ii);
                    e = db.enteries(index).word.e(ii);
                    Odata(count).data = data(b:e);
                    Odata(count).info = getmeta(db,index);
                    Odata(count).info.word = strrep(db.enteries(index).word.name(ii,1:end),' ','');
                    Odata(count).info.wordNumber = ii;
                end
            end
        end

        if db.kind==3 %reads phomemes
            [data smpr] = sentence_read(db,index);
            for ii=1:length(db.enteries(index).phoneme.b)
                b = db.enteries(index).phoneme.b(ii);
                e = db.enteries(index).phoneme.e(ii);
                in_word = db.enteries(index).word.flag(db.enteries(index).phoneme.from(ii));
                if(db.enteries(index).phoneme.flag(ii)==1)&&(in_word==1)
                    count = count + 1;
                    Odata(count).data = data(b:e);
                    Odata(count).info = getmeta(db,index);
                    Odata(count).info.word = strrep(db.enteries(index).word.name(db.enteries(index).phoneme.from(ii),1:end),' ','');
                    Odata(count).info.phoneme = strrep(db.enteries(index).phoneme.name(ii,1:end),' ','');
                end
            end
        end
    end
    %converts resalt in the from of struct to cellarray
    data = cell(length(Odata),1);
    meta = cell(length(Odata),1);
    for ii=1:length(Odata)
        data{ii} = Odata(ii).data;
        meta{ii} = Odata(ii).info;
    end
    fprintf('Enteries: %d\n', length(Odata));
else
    fprintf('An empty Database!\n');
    data = [];
    meta =[];
    smpr = [];
end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function Y=reshape(varargin)
%RESHAPE (overloaded)

Y = varargin{1};
% simple hack to figure out new sizes
temp = reshape(ones(size(Y.basis,1),1),varargin{2:end});
Y.conicinfo = [0 0];
Y.dim = size(temp);
Y = clean(Y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% READELP - read electrode locations from an .elp (electrode positions)
%             file as generated, for example, by a Polhemus tracking device
% Usage:
%   >> [eloc, elocnames, X, Y, Z] = readelp(filename);
%
% Inputs:
%   filename - name of the .elp file containing cartesian (XYZ) electrode
%              locations
% Outputs:
%   eloc      - structure containing the names and locations of the channels
%   elocnames - cell array containing the names of the electrodes
%   X,Y,Z     - vectors containing the xyz locations of the electrodes
%
% Author: Arnaud Delorme, CNL / Salk Institute, 28 Feb 2002
%
% Note: ignores comments and the sensor type field
% Note: convert output XYZ locations to polar coordinates using CART2POL
%
% See also: READLOCS, SNAPREAD, FLOATREAD, CART2POL

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 28 Feb 2002
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [eloc, names, x, y, z] = readelp( filename );

if nargin < 1
	help readelp;
	return;
end

% open file
% ---------
fid = fopen(filename, 'r');
if fid == -1
  disp('Cannot open file'); return;
end

index = 1;
countfid  = 1;
fidlabels = { 'Nz' 'LPA' 'RPA' };

needToReadNumbers = 0;
tmpstr = fgetl(fid);

while 1
    if needToReadNumbers==1 % Has sparsed $N.
        if (~isempty(str2num(tmpstr)))
            tmp = sscanf(tmpstr, '%f');

            eloc(index).X  = tmp(1); x(index) = tmp(1);
            eloc(index).Y  = tmp(2); y(index) = tmp(2);
            eloc(index).Z  = tmp(3); z(index) = tmp(3);
            eloc(index).type = 'EEG';

            index = index + 1;
            needToReadNumbers = 0;
        end
    elseif tmpstr(1) == '%'
        if tmpstr(2) == 'F' % fiducial
            tmp = sscanf(tmpstr(3:end), '%f');

            eloc(index).labels = fidlabels{countfid};
            eloc(index).X  = tmp(1); x(index) = tmp(1);
            eloc(index).Y  = tmp(2); y(index) = tmp(2);
            eloc(index).Z  = tmp(3); z(index) = tmp(3);
            eloc(index).type = 'FID';

            index     = index    + 1;
            countfid  = countfid + 1;

        elseif tmpstr(2) == 'N' % regular channel
            nm = strtok(tmpstr(3:end));
            if ~(strcmp(nm, 'Name') || strcmp(nm, 'Status'))
                eloc(index).labels = nm;
                needToReadNumbers = 1;
            end
        end
    end

    % Get the next line
    tmpstr = fgetl(fid);
    while isempty(tmpstr)
        tmpstr = fgetl(fid);
        if ~ischar(tmpstr) && tmpstr == -1
            break;
        end;
    end

    if ~ischar(tmpstr) && tmpstr == -1
        break;
    end;
end;

names = { eloc.labels };

fclose(fid);    % close the file

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% Reduce polygon to a given number of vertices
%
% poly = reduce_poly(poly, num)
%
% Inputs: poly        Polygon (2 rows, n columns)
%         num         Required number of vertices
%
% Outputs: poly       Final polygon
%
% Description: This code reduces the number of vertices in a closed polygon
% to the number specified by 'num'. It does this by calculating the
% importance of each vertex based on angle and segment length and then
% removing the least important. The process is repeated until the desired
% number of vertices is reached.
%
% Example:
% t = 0:0.1:2*pi;
% poly1 = [sin(t); cos(t)];
% poly2 = reduce_poly(poly1, 21);
% poly_draw = [poly2 poly2(:,1)];
% plot(poly_draw(1,:), poly_draw(2,:), '.-')
% axis equal
%
% Coded by: Peter Bone (peterbone@hotmail.com)
%------------------------------------------------------------------------
function poly = reduce_poly(poly, num)

numv = length(poly);

% Calculate initial importance of each vertex
imp = zeros(1,numv);
for v = 1 : numv
    imp(v) = vertex_importance(v, poly, numv);
end

% Iterate until desired number of vertices is reached
while numv > num

    [~, i] = min(imp(1:numv));

    % Remove vertex with least importance
    if i < numv
        poly(:,i:numv-1) = poly(:,i+1:numv);
        imp(i:numv-1) = imp(i+1:numv);
        vp = i;
    else
        vp = 1;
    end
    numv = numv - 1;

    % Recalculate importance for vertices neighbouring the removed one
    vm = 1 + mod(i - 2, numv);
    imp(vp) = vertex_importance(vp, poly, numv);
    imp(vm) = vertex_importance(vm, poly, numv);

end

% Clip polygon to the final length
poly = poly(:,1:num);

function a = vertex_importance(v, poly, numv)

% Find adjacent vertices
vp = 1 + mod(v, numv);
vm = 1 + mod(v - 2, numv);

% Obtain adjacent line segments and their lengths
dir1 = poly(:,v) - poly(:,vm);
dir2 = poly(:,vp) - poly(:,v);
len1 = norm(dir1);
len2 = norm(dir2);

% Calculate angle between vectors and multiply by segment lengths
% This is the importance of the vertex.
% Vertices with large angle and large segments attached are less
% likely to be removed
len1len2 = len1 * len2;
a = abs(acos((dir1' * dir2) / len1len2)) * len1len2;
%a = abs(1 - ((dir1' * dir2) / len1len2)) * len1len2;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [net,stats] = cnn_train_dag(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.averageImage = [];
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.gamma = 0.9;
opts.constraint = 10;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.saveMomentum = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train',opts) ;
    [net, state] = processEpoch(net, state, params, 'val',opts) ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train',opts) ;
      [net, state] = processEpoch(net, state, params, 'val',opts) ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode,opts)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.momentum)
  state.r = num2cell(zeros(1, numel(net.params))) ;
  state.momentum = num2cell(zeros(1, numel(net.params))) ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  state.momentum = cellfun(@gpuArray, state.momentum, 'uniformoutput', false) ;
  %state.r = cellfun(@gpuArray, state.r, 'uniformoutput', false) ;
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end
    opts.epoch = params.epoch;
    inputs = params.getBatch(params.imdb, batch,opts) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch,opts) ;
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (s ~= 1) ;
      net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    state = accumulateGradients(net, state, params, batchSize, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveMomentum
  state.momentum = [] ;
else
  state.momentum = cellfun(@gather, state.momentum, 'uniformoutput', false) ;
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
    %parDer(parDer>params.constraint) = params.constraint;
    %parDer(parDer<-params.constraint) = -params.constraint;
  end

  % skip
  thisLR = net.params(p).learningRate ;
  if(thisLR ==0 )
      continue;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.params(p).value) ;

        % Update momentum.
        state.momentum{p} = vl_taccum(...
          params.momentum, state.momentum{p}, ...
          -1, parDer) ;

        % Nesterov update (aka one step ahead).
        if params.nesterovUpdate
          delta = vl_taccum(...
            params.momentum, state.momentum{p}, ...
            -1, parDer) ;
        else
          delta = state.momentum{p} ;
        end

        % Update parameters.
        net.params(p).value = vl_taccum(...
          1,  net.params(p).value, thisLR, delta) ;
      end
    case 'rmsprop'
      opts = params;
      thisDecay = opts.weightDecay * net.params(p).weightDecay ;
      thisLR = opts.learningRate * net.params(p).learningRate ;
      state.r{p} = (1 - opts.gamma) * (net.params(p).der).^2 + opts.gamma * state.r{p};
      der_constraint = (1 / batchSize) * net.params(p).der./ (sqrt(state.r{p})+1e-8);
      der_constraint(der_constraint>opts.constraint) = opts.constraint;
      der_constraint(der_constraint<-opts.constraint) = -opts.constraint;
      state.momentum{p} = opts.momentum * state.momentum{p} ...
          - thisDecay * net.params(p).value ...
          - der_constraint ;
      net.params(p).value = net.params(p).value  + thisLR * state.momentum{p} ;
    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
state.r = [];
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end

author: layumi

repo: Image-Text-Embedding


text:
% created: Zoya Bylinskii, March 6
% based on: Kummerer et al.
% (http://www.pnas.org/content/112/52/16054.abstract)

% This finds the information-gain of the saliencyMap over a baselineMap

function score = InfoGain(saliencyMap, fixationMap, baselineMap)
% saliencyMap is the saliency map
% fixationMap is the human fixation map (binary matrix)
% baselineMap is another saliency map (e.g. all fixations from other images)

%%
map1 = imresize(saliencyMap,size(fixationMap));
mapb = imresize(baselineMap,size(fixationMap));
%%
% normalize and vectorize saliency maps
map1 = (map1(:) - min(map1(:)))/(max(map1(:))-min(map1(:)));
mapb = (mapb(:) - min(mapb(:)))/(max(mapb(:))-min(mapb(:)));
%%

% turn into distributions
map1 = map1./sum(map1);
mapb = mapb./sum(mapb);

locs = logical(fixationMap(:));
%%
score = mean(log2(eps+map1(locs))-log2(eps+mapb(locs)));

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function [out_im,status ] = open_bitfield_bmp( filename )
% OPEN_BITFIELD_BMP open a bitfield compressed bitmap image.
%
%   IM = OPEN_BITFIELD_BMP( FILENAME ) opens a 16-bit bitmapped compressed file
%   named FILENAME. The output is formatted as three planes of Red, Green
%   and Blue data.
%
%   See also imread, imfinfo

% Get information about the image
info = imfinfo(filename);

% if we have a bitfield compressed R-G-B bitmap image
if ( strcmp(info.Format , 'bmp') && ...
     strcmp(info.CompressionType , 'bitfields') &&  ...
     (info.BitDepth == 16 || info.BitDepth == 32) )

    % indicate the input bmp file is not of 8-bit
    status = 0;

    % Open the file for reading
    fid = fopen(filename, 'r');

    % Extract relevvant image info
    data_offset = info.ImageDataOffset;
    width = info.Width;
    height = info.Height;

    % Create space for output image
    out_im = zeros(height, width, 3);

    % Seek to where the image data begins (i.e. skip the file header
    fseek(fid, data_offset, 'bof');

    % Read in the image data and format it into a matrix
    if (info.BitDepth == 16)
        % compressed_image = (fread(fid, [width + 1, height] , 'uint16'))';
        compressed_image = (fread(fid, [width, height] , 'uint16'))';
    else
        compressed_image = (fread(fid, [width, height] , 'uint32'))';
    end

    % Eliminate last column of junk data (scanline row terminators) in 16
    % bit mode
    if (info.BitDepth == 16)
        compressed_image = compressed_image(:, 1:(width));
    end;

    % Invert Row Order since it is flipped
    new_image = flipud(compressed_image);

    % Extract color bitmasks to decompress
    red_mask = info.RedMask;
    blue_mask = info.BlueMask;
    green_mask = info.GreenMask;

    % Extract color components and form output image
    out_im(:,:,1) = (bitand(new_image, red_mask) / red_mask * 255);
    out_im(:,:,2) = (bitand(new_image, green_mask) / green_mask * 255);
    out_im(:,:,3) = (bitand(new_image, blue_mask) / blue_mask * 255);

    % typecast
    out_im = uint8(out_im);

    % display image
    %imshow(out_im);

else

    % indicate the input bmp file is 8-bit
    status = 1;

    % general imread for all other cases
    out_im = imread(filename);
    %imshow(out_im);

end

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


text:
function v = ToVector(im)
sz = size(im);
if size(im,3)==3
    v = reshape(im, [prod(sz(1:2)) 3]);
else
    v = reshape(im, [prod(sz(1:2)) 1]);
end


text:
function [pID,pN] = FDR(p,q)
% :Usage:
% ::
%
%     pt = FDR(p,q)
%
% :Inputs:
%
%   **p:**
%        vector of p-values
%
%   **q:**
%        False Discovery Rate level
%
% :Outputs:
%
%   **pID:**
%        p-value threshold based on independence or positive dependence
%
%   **pN:**
%        Nonparametric p-value threshold
%
% ..
%    % @(#)FDR.m	1.3 Tom Nichols 02/01/18
%
%    The checking code below was added by Tor Wager
% ..

    p(isnan(p)) = [];

    if any(p == 0)
        disp('******************************************')
        disp('Warning! Some p-values are zero.')
        disp('FDR.m will interpret these as ineligible voxels.')
        disp('If these are valid p-values, they should have some not-exactly-zero value.')
        p(p == 0) = [];
        disp('******************************************')

    end

    p = sort(p(:));
    V = length(p);
    I = (1:V)';

    cVID = 1;
    cVN = sum(1./(1:V));

    pID = p(max(find(p<=I/V*q/cVID)));
    pN = p(max(find(p<=I/V*q/cVN)));

    return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore