text:
function rf=v_lpcla2rf(la)
%V_LPCLA2RF Convert log areas to reflection coefficients RF=(LA)

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

[nf,p2]=size(la);
rf=-tanh((la(:,1:p2-1)-la(:,2:p2))/2);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function val = get_params_kernel(K, name)
% GET_PARAMS_KERNEL Accessor function for a field (tabular_kernel)
% val = get_params_kernel(K, name)
%
% e.g., get_params_kernel(K, 'table')

switch name
 case 'table', val = K.table;
 otherwise,
  error(['invalid field name ' name]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [Pro,Res] = interpolationAMGsa(A,node2agg,omega,smoothingstep)
%% INTERPOLATIONAMGSA construct prolongation using smoothed aggregation
%
% [Pro,Res] = INTERPOLATIONAMGSA(A,node2agg) construct prolongation and
% restriction matrices using smoothed aggregation.
%
% In the input, A is a SPD matrix and node2agg assigns an aggregate for
% each node in A. In the output Pro and Res are prolongation and
% restriction matrices satisfying Res = Pro'.
%
% The prolongation operator is first set to be piecewise constant and then
% using matrix A to smooth out this simple prolongation using weighted
% Jacobi iteration for several steps. The default choice of the weight is
% 0.35 and the smoothingstep is 2. Larger smoothing steps will result a
% denser prolongation operator.
%
% Example
%   load lakemesh
%   A = assemblematrix(node,elem);
%   [node2agg,As] = coarsenAMGa(A);
%   [Pro,Res] = interpolationAMGsa(As,node2agg);
%
% See also: coarsenAMGc, amg
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('omega','var'), omega = 0.35; end
if ~exist('smoothingstep','var'), smoothingstep = 2; end

%% A simple prolongation
N = size(A,1);
Nc = max(node2agg);
idx = find(node2agg~=0);
Pro = sparse(idx,node2agg(idx),1,N,Nc);

%% Smooth the piecewise constant prolongation
for k = 1:smoothingstep
    Pro = Pro - omega*(A*Pro);
end

%% Normalize the prolongation such that the row sum is one
rowsum = sum(Pro,2);
D = spdiags(1./rowsum,0,N,N);
Pro = D*Pro;
Res = Pro';

repo: ifem


text:
classdef ConstantDetectionProbabilityX < DetectionModelX
% ConstantDetectionProbabilityX class
%
% Summary of ConstantDetectionProbabilityX
% This is a class implementation of a detection model, described by a Poisson
% distributed false alarm rate and a Uniform spatial distribution.
%
% ConstantDetectionProbabilityX Properties:
%   + DetectionProbability - The detection probability.
%
% ConstantDetectionProbabilityX Methods:
%   + pdf       - Function to evaluate the instensity p(b_t) at a given point in
%                 the state space.
%
% (+) denotes puplic properties/methods
%
%  February 2018 Lyudmil Vladimirov, University of Liverpool

    properties
        DetectionProbability
    end

    methods (Access = protected)
        function initialise_(this, config)
            if(isfield(config,'DetectionProbability'))
                this.DetectionProbability = config.DetectionProbability;
            end
        end
    end

    methods
        function this = ConstantDetectionProbabilityX(varargin)
        % ConstantDetectionProbabilityX Constructor method
        %
        % Parameters
        % ----------
        % DetectionProbability: scalar
        %   The detection probability over the search space.
        %
        % See also pdf.

            % Call SuperClass method
            this@DetectionModelX();

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

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

        end

        function int = pdf(this, varargin)
        % pdf Evaluates the detection intensity at given points xk of the state space
        %
        % Parameters
        % ----------
        % xk: matrix
        %   A (NumStateDims x Ns) matrix, whose columns correspond to Ns individual state vectors.
        %
        % Returns
        % -------
        % prob: matrix
        %   A (1 x Ns) matrix, where :math:`prob(i,j)=p(D_k^i|x_k^j)

            Ns = 1;
            if(nargin>1)
                xk = varargin{1};
                Ns = size(xk,2);
            end
            int = this.DetectionProbability*ones(1,Ns);
            x = [-3422.85261310741,-3191.58441229017,-2993.55608122595,-3228.08558459284, -3422.85261310741];
            y = [1.472966334316646e+03,1.123217741935625e+03,1.243952901078573e+03,1.579667558371468e+03, 1.472966334316646e+03];
            a = find(inpolygon(xk(1,:),xk(3,:),x,y));
            if numel(a)>0
                int(a) = 0.1;
            end
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function str = VBA_summaryMFX(out)
% writes a summary string from output of VBA_MFX model inversion
% function str = VBA_summaryMFX(out)
% IN:
%   - out: the 'out' structure of VBA_MFX inversion routine
% OUT:
%   - str: a cell array of strings, which summarize the VBA inversion

str{1} = sprintf(['Date: ',datestr(out.date)]);
s0 = ['VB converged in ',num2str(out.it),' iterations'];
try
    if floor(out.dt./60) == 0
        timeString = [num2str(floor(out.dt)),' sec'];
    else
        timeString = [num2str(floor(out.dt./60)),' min'];
    end
    str{2} = sprintf([s0,' (took ~',timeString,')']);
catch
    str{2} = sprintf(s0);
end
str{3} = sprintf(['Dimensions of the MFX analysis:','\n ',...
    '    - subjects: ns=',num2str(out.options.dim.ns),'\n ',...
    '    - data: p=',num2str(out.options.dim.p),'\n ',...
    '    - hidden states: n=',num2str(out.options.dim.n),'\n ',...
    '    - evolution parameters: n_theta=',num2str(out.options.dim.n_theta),'\n ',...
    '    - observation parameters: n_phi=',num2str(out.options.dim.n_phi)]);
str{4} = sprintf('\n     Within-subject generative model:\n');
if isa(out.options.g_fname,'function_handle')
    gfn = func2str(out.options.g_fname);
else
    gfn = out.options.g_fname;
end
if out.options.dim.n >= 1
    if isa(out.options.f_fname,'function_handle')
        ffn = func2str(out.options.f_fname);
    else
        ffn = out.options.f_fname;
    end
    str{4} = sprintf([str{4},...
        '       - observation function: ',gfn,'\n',...
        '       - evolution function: ',ffn]);
else
    str{4} = sprintf([str{4},'     - observation function: ',gfn]);
end
mF = mean(out.within_fit.F(:));
sF = std(out.within_fit.F(:));
mF0 = mean(out.within_fit.LLH0(:));
sF0 = std(out.within_fit.LLH0(:));
str{5} = sprintf([...
    '\n     - Bayesian log model evidences: <log p(y|m)> = ',num2str(mF,'%4.3e'),' +/- ',num2str(sF,'%4.3e'),'\n',...
    '     - Bayesian log evidences under the null: <log p(y|H0)> = ',num2str(mF0,'%4.3e'),' +/- ',num2str(sF0,'%4.3e')]);

if any([out.options.sources.type]==0)
    R2str = 'coefficient of determination: <R2>';
    mR = mean(out.within_fit.R2);
    sR = std(out.within_fit.R2);
    str{6} = sprintf([...
        '\n     - ',R2str,' = ',num2str(mR,'%4.3f '),' +/- ',num2str(sR,'%4.3e ')]);
end

if any([out.options.sources.type]>0)
    R2str = 'balanced classification accuracy <Acc>';
    mR = mean(out.within_fit.R2);
    sR = std(out.within_fit.R2);
    str{7} = sprintf([...
        '\n     - ',R2str,' = ',num2str(mR,'%4.3f '),' +/- ',num2str(sR,'%4.3f ')]);
end

author: MBB-team

repo: VBA-toolbox


text:
% This example demonstrates how to use the auxilliary routine
% (postprocess.m) to calculate the four remaining field
% components (Hz, Ex, Ey, Ez) once the transverse magnetic
% field components (Hx, Hy) are known.  The waveguide
% considered here is a uniaxial channel waveguide with the
% c-axis oriented at pi/4 relative to the x-y axes.

n1 = 1.55;
n2x = 2.156;
n2y = 2.232;
n2z = 2.232;
theta = pi/4;

e2xx = n2x^2*cos(theta)^2 + n2y^2*sin(theta)^2;
e2yy = n2y^2*cos(theta)^2 + n2x^2*sin(theta)^2;
e2xy = cos(theta)*sin(theta)*(n2x^2-n2y^2);
e2yx = e2xy;

Rx = 0.30;
Ry = 0.20;
side = 0.20;

dx = 0.0025;        % grid size (x)
dy = dx;            % grid size (y)

lambda = 1.00;      % wavelength
nmodes = 1;         % number of modes to compute

fprintf (1,'generating index mesh...\n');

[x,y,xc,yc,nx,ny,epsxx,edges] = ...
    waveguidemeshfull([n1,sqrt(e2xx),n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epsxy,edges] = ...
    waveguidemeshfull([0,sqrt(e2xy),0],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epsyx,edges] = ...
    waveguidemeshfull([0,sqrt(e2yx),0],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epsyy,edges] = ...
    waveguidemeshfull([n1,sqrt(e2yy),n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);
[x,y,xc,yc,nx,ny,epszz,edges] = ...
    waveguidemeshfull([n1,n2z,n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy);

% Now we stretch out the mesh at the boundaries:
[x,y,xc,yc,dx,dy] = stretchmesh(x,y,[80,80,80,80],[4,4,4,4]);

[Hx,Hy,neff] = wgmodes (lambda, n2y, nmodes, dx, dy, epsxx, epsxy, epsyx, epsyy, epszz, '0000');

fprintf(1,'neff = %7.5f\n',neff);

[Hz,Ex,Ey,Ez] = postprocess (lambda, neff, Hx, Hy, dx, dy, epsxx, epsxy, epsyx, epsyy, epszz, '0000');

figure(1);

subplot(231);
contourmode(x,y,Hx);
title('Hx');
for v = edges, line(v{:}); end

subplot(232);
contourmode(x,y,Hy);
title('Hy');
for v = edges, line(v{:}); end

subplot(233);
contourmode(x,y,Hz);
title('Hz');
for v = edges, line(v{:}); end

subplot(234);
contourmode(x,y,Ex);
title('Ex');
for v = edges, line(v{:}); end

subplot(235);
contourmode(x,y,Ey);
title('Ey');
for v = edges, line(v{:}); end

subplot(236);
contourmode(x,y,Ez);
title('Ez');
for v = edges, line(v{:}); end

figure(2);

ii = 1;
colormap(jet(256));
hn = abs(interp2(y,x,Hy,side+Ry,0));
en = abs(interp2(yc,xc,Ex,side+Ry,0));

subplot(231);
imagemode(x,y,Hx/hn);
title(sprintf('Hx (mode %d)',ii));
for v = edges, line(v{:}); end

subplot(232);
imagemode(x,y,Hy/hn);
title(sprintf('Hy (mode %d)',ii));
for v = edges, line(v{:}); end

subplot(233);
imagemode(x,y,Hz/hn);
title(sprintf('Hz (mode %d)',ii));
for v = edges, line(v{:}); end

subplot(234);
imagemode(x,y,Ex/en);
title(sprintf('Ex (mode %d)',ii));
for v = edges, line(v{:}); end

subplot(235);
imagemode(x,y,Ey/en);
title(sprintf('Ey (mode %d)',ii));
for v = edges, line(v{:}); end

subplot(236);
imagemode(x,y,Ez/en);
title(sprintf('Ez (mode %d)',ii));
for v = edges, line(v{:}); end

author: Sable

repo: mcbench-benchmarks


text:
function report(dir1, suffix1, showcls, do_auc_ap)
% Print scores for all classes.

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

if nargin < 3
  showcls = true;
end

if nargin < 4
  do_auc_ap = false;
end

files = dir([dir1 '/*' suffix1 '.mat']);
classes = cell(length(files), 1);
for i = 1:length(files)
  j = strfind(files(i).name, '_');
  classes{i} = files(i).name(1:j(1)-1);
end

count1 = 0;
for i=1:length(classes)
  cls = classes{i};
  try
    load([dir1 cls suffix1]);
    ap1 = ap;
    if do_auc_ap
      ap2 = xVOCap(recall, prec);
      score2(i) = ap2;
    end
    if showcls
      if do_auc_ap
        fprintf('%12s %.4f\n', cls, ap2);
      else
        fprintf('%12s %.4f\n', cls, ap1);
      end
    else
      if do_auc_ap
        fprintf('%.4f\n', ap2);
      else
        fprintf('%.4f\n', ap1);
      end
    end
    score1(i) = ap1;
  catch
    score1(i) = nan;
    score2(i) = nan;
    if showcls
      fprintf('%12s -\n', cls);
    else
      fprintf('-\n');
    end
  end
end

a1 = nanmean(score1);
if do_auc_ap
  a2 = nanmean(score2);
end
if showcls
  fprintf('%s\n', repmat('-', [1 12]));
  if do_auc_ap
    fprintf('%12s %.4f\n', 'mAP', a2);
  else
    fprintf('%12s %.4f\n', 'mAP', a1);
  end
else
  if do_auc_ap
    fprintf('%.4f\n', a2);
  else
    fprintf('%.4f\n', a1);
  end
end

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
% function Cpsi = cwt_adm(type, opt)
%
% Calculate cwt admissibility constant int(|f(w)|^2/w, w=0..inf) as
% per Eq. (4.67) of [1].
%
% 1. Mallat, S., Wavelet Tour of Signal Processing 3rd ed.
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function Cpsi = cwt_adm(type, opt)
    if nargin<2, opt=struct(); end
    switch type
      case 'sombrero',
        if ~isfield(opt,'s'), s = 1; else s = opt.s; end
        Cpsi = (4/3)*s*sqrt(pi);
      case 'shannon',
        Cpsi = log(2);
      otherwise
        psihfn = wfiltfn(type, opt);
        Cpsi = quadgk(@(x) (conj(psihfn(x)).*psihfn(x))./x, 0, Inf);
    end

    % Normalize
    Cpsi = Cpsi / (4*pi);
end

author: HeLiangHIT


text:
%BAGC Bag classifier for classifying sets of object instances
%
%		[WBAG,WOBJ] = BAGC(A,OBJCLASSF,BAGINDEX,BAGCOMBC,BAGCLASSF,BAGLAB)
%		 WBAG       = A*BAGC([],OBJCLASSF,BAGINDEX,BAGCOMBC,BAGCLASSF,BAGLAB)
%		 WBAG       = A*BAGC(OBJCLASSF,BAGINDEX,BAGCOMBC,BAGCLASSF,BAGLAB)
%		 D          = B*WBAG
%
% INPUT
%   A          Training dataset with object labels and bag indices
%   B          Test Dataset with index list of bags, stored as label list
%   OBJCLASSF  Trained or untrained object classifier, default QDC
%   BAGINDEX   String or a label_list_index defining in which label list
%              of A the bag indices are stored, default 2.
%   BAGCOMBC   Combiner for objects in a bag, default VOTEC
%   BAGCLASSF  Untrained classifier for bags, default FISHERC
%   BAGLAB     String or a label_list_index defining in which label list
%              of A the bag labels are stored. Objects with the same
%              bag index should have the same bag label.
%              Default is the current labeling of A
%
% OUTPUT
%   WBAG       Trained bag classifier
%   WOBJ       Trained object classifier
%   D          Classification matrix of bags in B
%
% DESCRIPTION
% This routine offers a classifier for bags (e.g. images) of objects
% (e.g. pixels) stored in a single dataset. The objects in the training
% set A should have at least two labels: bag labels (the class of their bag)
% and bag indices, defining which objects belong to the same bag. These two
% label sets should be stored by the ADDLABELS command in the dataset A.
% Refer to the multi-labeling system (see MULTI_LABELING) offered by
% PRTools. The current object labels of A can be the bag labels, but may
% also be different, e.g. true object labels.
%
% BAGINDEX should be a label_list_name or a label_list_index defining the
% label list used for storing the bag indices that refer to the bag an
% object belongs to. The same label_list_name or label_list_index should be
% used for defining the bags of the test objects in B.
%
% All objects in A are used to train the object classifier OBJCLASSF if it
% is untrained. The current object labels are used for that. Classification
% results of the objects in the same bag are combined by BAGCOMBC, which
% can be any of the fixed combiners MEANC, PRODC, PERC, MAXC, etcetera.
% This results for every bag in a single confidence vector for the classes.
%
% If an untrained bag classifier BAGCLASSF is supplied, the bag confidence
% vectors are used to train a bag classifier.
%
% New bags, organised in a dataset like B, with the proper bag indices per
% object stored in a label list with the same name or label_list_index as
% used in A, can be classified by the bag classifier WBAG.
%
% If no bag classifier BAGCLASSF was defined during training, just the
% results of the object classifier WOBJ are returned combined by BAGCOMBC
% over the objects in the same bag in B. In this case the final result is
% identical to B*(A*WOBJ)*BAGCC([],BAGCOMBC), provided that A has class
% labels and B is labeled by its bag indices.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, MULTI_LABELING, BAGCC, LOSO,
% DATASET/ADDLABELS, DATASET/CHANGELABLIST

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function [out1,out2] = bagc(varargin)

	argin = shiftargin(varargin,'prmapping');
  argin = setdefaults(argin,[],qdc,2,votec,[],[]);
  if mapping_task(argin,'definition')
    out1 = define_mapping(argin,'untrained');
    out1 = setname(out1,'Set classifier');
  else
    [a,objclassf,bagindex,bagcombc,bagclassf,baglab] = deal(argin{:});
    if isuntrained(objclassf) | nargin > 2 | ~strcmp(getmapping_file(objclassf),mfilename)
      % train the mapping (classifier)

      % we need datasets with at least 1 object per class and 2 classes
      isvaldset(a,1,2);

      % if the object classifier is untrained, train it, else use it
      if isuntrained(objclassf)
        wobj = a*objclassf;
      else
        wobj = objclassf;
      end

      if ismapping(bagclassf) & isuntrained(bagclassf)

        % if the bag labels are not given,
        % use the objcts labels for the bags too
        if isempty(setlab), setlab = curlablist(a); end

        % classifiy the dataset and change labeling to bag index
        x = changelablist(a*wobj,setindex);

        % avoid empty bags
        x = setlablist(x);

        % combine object results to bag results
        d = bagcc(x,bagcombc);

        % change to bag labels
        d = changelablist(d,baglab);

        % train bag classifier
        bagclassf = d*bagclassf;

        % get outputlabels
        labels_out = getlabels(bagclassf);

      else
        labels_out = getlabels(wobj);
      end

      % store all what is needed for execution in the mapping
      out1 = prmapping(mfilename,'trained',{wobj,bagcombc,bagclassf,bagindex,baglab}, ...
        labels_out, size(a,2),size(labels_out,1));

      % prevent batch execution
      out1 = setbatch(out1,0);

      % return the object classifier as well
      out2 = wobj;

    else % here we are for execution

      % the mapping is stored in objclassf
      w = getdata(objclassf);

      % save current lablist
      curlist = curlablist(a);

      % use the bag index for the test set if supplied
      if ~isempty(w{4}), testset = changelablist(a,w{4}); end

      % classify test set by the object classifier
      d = testset*w{1};

      % avoid empty bags
      d = setlablist(d);

      % combine objects in the same bag
      d = bagcc(d,w{2});

      % reset lablist for classification matrix
      d = changelablist(d,curlist);

      % apply the set classifier, if defined
      if ~isempty(w{3}), d = d*w{3}; end

      % that is it, define class labels as feature labels
      out1 = setfeatlab(d,getlabels(objclassf));

    end

  end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% test the main formulation

addpath('../MALSAR/c_files/largescale_ops/');
addpath('../MALSAR/functions/pacifier/');

clear; clc;
rng(1985)
n = 12;
k = 40;
d = 1000;
t_min = 150;
t_max = 200;
density = 0.001;

X0 = cell(n, 1);
for i = 1 : n
    X0{i} = sprand(d, t_min + randi(t_max-t_min, 1), density);
end

reg_l1 = 1e-10;
reg_l2 = 1e-10;

reg_smooth = 10;

[ U1, V1, Ss1, fv1 ] = pacifier_iba( X0, k+5, reg_l1, reg_l2, reg_smooth);
[ U2, V2, Ss2, fv2 ] = pacifier_sba( X0, k+5, reg_l1, reg_l2, reg_smooth);

figure
plot(fv1)
title('Pacifier IBA Objective Value')

figure
plot(fv2)
title('Pacifier SBA Objective Value')

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [x, hist] = chordsol(x, f, tola, tolr)
%CHORDSOL  Chord method code for scalar equations.
%
% This code terminates on small relative-absolute errors
%
% function [X, HIST] = CHORDSOL(X, F, TOLA, TOLR)
%
% Inputs: X = initial iterate
%         F = function
%         TOLA = absolute error tolerance
%         TOLR = relative error tolerance
%
% Output: x=approximate result
%      HIST = array of iteration history, useful for tables and plots
%                The two columns are iteration number and residual.
%
%                If you leave this out by calling chordsol as
%                Z = CHORDSOL(X, F, TOLA, TOLR)
%                the storage for hist is not allocated and
%                the iteration history is not stored.
%
% Uses: fordiff.m, Tim's numerical derivative
%
%
maxit=20;    % maximum iterations
alpha=1.d-4;
small=1.d-8;
if nargout == 2
    histmp=zeros(maxit,2);
end
%
% Initialize
%

itc=0;
fc=feval(f,x);
fp=fordiff(f,x,fc);
tol=tolr*abs(fc)+tola;
%
%    Store iteration history?
%
if nargout == 2
    histmp(itc+1,:)=[itc, abs(fc)];
end
%
%    Main chord loop.
%
while(abs(fc) > tol)
%
%    Iteration becoming unbounded?
%
    if(abs(x) > 1.d9)
        disp(' iterate too large');
        if nargout == 2
            histmp(itc+1,:)=[itc,abs(fc)];
            hist=histmp(1:itc,:);
        end
    return
    end
%
%    Check for small relative derivatives.
%
    if(abs(fp) <= small*abs(fc))
        disp(' small derivative error');
    return;
    end
    s=-fc/fp;
    xt=x+s;
    ft=feval(f,xt);
    x=xt;
    itc=itc+1;
%
%    Too many iterations?
%
        if(itc > maxit)
            disp(' maxit reached');
            if nargout == 2
                histmp(itc+1,:)=[itc,abs(fc)];
                hist=histmp(1:itc,:);
            end
            return
        end
%
        fc=feval(f,x);
    if nargout == 2
        histmp(itc+1,:)=[itc,abs(fc)];
        end
    end
        if nargout == 2
        hist=histmp(1:itc+1,:);
        end
%
% A simple forward difference in MATLAB.
%
function fp = fordiff(f,x,fc)
        h=1.d-7;
%
%       Notice how f is evaluated using the feval command.
%
        fright=feval(f,x+h);
        fp=(fright-fc)/h;

repo: Aero-matlab


text:
classdef(Abstract) AbstractGeometricVectorConstraint < AbstractConstraint
    %AbstractGeometricVectorConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        vector AbstractGeometricVector
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        type(1,:) char

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            if(not(isempty(obj.frame)))
                frame = obj.frame;
            else
                frame = stateLogEntry.centralBody.getBodyCenteredInertialFrame();
            end

            value = lvd_GeometricVectorTasks(stateLogEntry, obj.type, obj.vector, frame);

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                valueStateComp = lvd_GeometricVectorTasks(stateLogEntryStateComp, obj.type, obj.vector, frame);
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(~, ~)
            tf = false;
        end

        function tf = usesExtremum(~, ~)
            tf = false;
        end

        function tf = usesGroundObj(~, ~)
            tf = false;
        end

        function tf = usesGeometricVector(obj, vector)
            tf = obj.vector == vector;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = '';
            lbLim = 0;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
            if(lvdData.geometry.vectors.getNumVectors() >= 1)
%                 addConstraintTf = lvd_EditGeometricVectorConstraintGUI(obj, lvdData);

                output = AppDesignerGUIOutput({false});
                lvd_EditGeometricVectorConstraintGUI_App(obj, lvdData, output);
                addConstraintTf = output.output{1};
            else
                errordlg('There are currently no geometric vectors in this scenario.  Add at least one new vector first.');

                addConstraintTf = false;
            end
        end

        function vector = selectConstraintObj(obj, lvdData)
            [listBoxStr, vectors] = lvdData.geometry.vectors.getListboxStr();

            vector = [];
            if(isempty(vectors))
                warndlg('Cannot create vector value object: no vectors have been created.  Create a vector first.','Vector Value Constraint','modal');
            else
                [Selection,ok] = listdlg('PromptString',{'Select a vector:'},...
                                'SelectionMode','single',...
                                'Name','Vectors',...
                                'ListString',listBoxStr);

                if(ok == 0)
                    vector = [];
                else
                    vector = vectors(Selection);
                end
            end
        end

        function useObjFcn = setupForUseAsObjectiveFcn(obj,lvdData)
            vectorSel = obj.selectConstraintObj(lvdData);

            if(not(isempty(vectorSel)))
                obj.vector = vectorSel;
                useObjFcn = true;
            else
                useObjFcn = false;
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function d = dell2nose_pair(trx,fly1,fly2,istry)

nsamples = 20;

% initialize
d = nan(1,trx(fly1).nframes);

% get start and end frames of overlap
t0 = max(trx(fly1).firstframe,trx(fly2).firstframe);
t1 = min(trx(fly1).endframe,trx(fly2).endframe);

% no overlap
if t1 < t0,
  return;
end

% position of nose2
xnose = trx(fly2).x_mm + 2*trx(fly2).a_mm.*cos(trx(fly2).theta_mm);
ynose = trx(fly2).y_mm + 2*trx(fly2).a_mm.*sin(trx(fly2).theta_mm);

% ellipse 1
x_mm1 = trx(fly1).x_mm;
y_mm1 = trx(fly1).y_mm;
a_mm1 = trx(fly1).a_mm;
b_mm1 = trx(fly1).b_mm;
theta_mm1 = trx(fly1).theta_mm;

off1 = trx(fly1).off;
off2 = trx(fly2).off;

if nargin < 4,
  tstry = t0:t1;
else
  tstry = istry(:)' - off1;
end

for t = tstry,
  i = t + off1;
  j = t + off2;
  d(i) = ellipsedist_hack(x_mm1(i),y_mm1(i),...
    2*a_mm1(i),2*b_mm1(i),theta_mm1(i),...
    xnose(j),ynose(j),nsamples);
end

author: kristinbranson

repo: JAABA


text:
function varargout = dim(varargin)
% Format
% For getting the value
% dat = dim(obj)
%
% For setting the value
% obj = dim(obj,dat)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

if nargin==2,
    varargout{1} = asgn(varargin{:});
elseif nargin==1,
    varargout{1} = ref(varargin{:});
else
    error('Wrong number of arguments.');
end;
return;

function dat = ref(obj)
dat = obj.dim;
return;

function obj = asgn(obj,dat)
if isnumeric(dat) && all(dat>=0) && all(rem(dat,1)==0),
    dat = [double(dat(:)') 1 1];
    lim = max([2 find(dat~=1)]);
    dat = dat(1:lim);
    obj.dim = dat;
else
    error('"dim" must be a vector of positive integers.');
end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x,flag,relres,iter,resvec] = Pgmres_deflation(A,b,restart,tol,maxit,M1,M2,x,varargin)
%GMRES   Generalized Minimum Residual Method.
%   X = GMRES(A,B) attempts to solve the system of linear equations A*X = B
%   for X.  The N-by-N coefficient matrix A must be square and the right
%   hand side column vector B must have length N. This uses the unrestarted
%   method with MIN(N,10) total iterations.
%
%   X = GMRES(AFUN,B) accepts a function handle AFUN instead of the matrix
%   A. AFUN(X) accepts a vector input X and returns the matrix-vector
%   product A*X. In all of the following syntaxes, you can replace A by
%   AFUN.
%
%   X = GMRES(A,B,RESTART) restarts the method every RESTART iterations.
%   If RESTART is N or [] then GMRES uses the unrestarted method as above.
%
%   X = GMRES(A,B,RESTART,TOL) specifies the tolerance of the method.  If
%   TOL is [] then GMRES uses the default, 1e-6.
%
%   X = GMRES(A,B,RESTART,TOL,MAXIT) specifies the maximum number of outer
%   iterations. Note: the total number of iterations is RESTART*MAXIT. If
%   MAXIT is [] then GMRES uses the default, MIN(N/RESTART,10). If RESTART
%   is N or [] then the total number of iterations is MAXIT.
%
%   X = GMRES(A,B,RESTART,TOL,MAXIT,M) and
%   X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2) use preconditioner M or M=M1*M2
%   and effectively solve the system inv(M)*A*X = inv(M)*B for X. If M is
%   [] then a preconditioner is not applied.  M may be a function handle
%   returning M\X.
%
%   X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) specifies the first initial
%   guess. If X0 is [] then GMRES uses the default, an all zero vector.
%
%   [X,FLAG] = GMRES(A,B,...) also returns a convergence FLAG:
%    0 GMRES converged to the desired tolerance TOL within MAXIT iterations.
%    1 GMRES iterated MAXIT times but did not converge.
%    2 preconditioner M was ill-conditioned.
%    3 GMRES stagnated (two consecutive iterates were the same).
%
%   [X,FLAG,RELRES] = GMRES(A,B,...) also returns the relative residual
%   NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL. Note with
%   preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).
%
%   [X,FLAG,RELRES,ITER] = GMRES(A,B,...) also returns both the outer and
%   inner iteration numbers at which X was computed: 0 <= ITER(1) <= MAXIT
%   and 0 <= ITER(2) <= RESTART.
%
%   [X,FLAG,RELRES,ITER,RESVEC] = GMRES(A,B,...) also returns a vector of
%   the residual norms at each inner iteration, including NORM(B-A*X0).
%   Note with preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).
%
%   Example:
%      n = 21; A = gallery('wilk',n);  b = sum(A,2);
%      tol = 1e-12;  maxit = 15; M = diag([10:-1:1 1 1:10]);
%      x = gmres(A,b,10,tol,maxit,M);
%   Or, use this matrix-vector product function
%      %-----------------------------------------------------------------%
%      function y = afun(x,n)
%      y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0];
%      %-----------------------------------------------------------------%
%   and this preconditioner backsolve function
%      %------------------------------------------%
%      function y = mfun(r,n)
%      y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
%      %------------------------------------------%
%   as inputs to GMRES:
%      x1 = gmres(@(x)afun(x,n),b,10,tol,maxit,@(x)mfun(x,n));
%
%   Class support for inputs A,B,M1,M2,X0 and the output of AFUN:
%      float: double
%
%   See also BICG, BICGSTAB, BICGSTABL, CGS, LSQR, MINRES, PCG, QMR, SYMMLQ,
%   TFQMR, ILU, FUNCTION_HANDLE.

%   References
%   H.F. Walker, "Implementation of the GMRES Method Using Householder
%   Transformations", SIAM J. Sci. Comp. Vol 9. No 1. January 1988.

%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.21.4.15 $ $Date: 2011/05/17 02:33:07 $
global ZH AH
if (nargin < 2)
    error(message('MATLAB:gmres:NumInputs'));
end

% Determine whether A is a matrix or a function.
[atype,afun,afcnstr] = iterchk(A);
if strcmp(atype,'matrix')
    % Check matrix and right hand side vector inputs have appropriate sizes
    [m,n] = size(A);
    if (m ~= n)
        error(message('MATLAB:gmres:SquareMatrix'));
    end
    if ~isequal(size(b),[m,1])
        error(message('MATLAB:gmres:VectorSize', m));
    end
else
    m = size(b,1);
    n = m;
    if ~iscolumn(b)
        error(message('MATLAB:gmres:Vector'));
    end
end

% Assign default values to unspecified parameters
if (nargin < 3) || isempty(restart) || (restart == n)
    restarted = false;
else
    restarted = true;
end
if (nargin < 4) || isempty(tol)
    tol = 1e-6;
end
warned = 0;
if tol < eps
    warning(message('MATLAB:gmres:tooSmallTolerance'));
    warned = 1;
    tol = eps;
elseif tol >= 1
    warning(message('MATLAB:gmres:tooBigTolerance'));
    warned = 1;
    tol = 1-eps;
end
if (nargin < 5) || isempty(maxit)
    if restarted
        maxit = min(ceil(n/restart),10);
    else
        maxit = min(n,10);
    end
end

if restarted
    outer = maxit;
    if restart > n
        warning(message('MATLAB:gmres:tooManyInnerItsRestart',restart, n));
        restart = n;
    end
    inner = restart;
else
    outer = 1;
    if maxit > n
        warning(message('MATLAB:gmres:tooManyInnerItsMaxit',maxit, n));
        maxit = n;
    end
    inner = maxit;
end

% Check for all zero right hand side vector => all zero solution
n2b = norm(b);                   % Norm of rhs vector, b
if (n2b == 0)                    % if    rhs vector is all zeros
    x = zeros(n,1);              % then  solution is all zeros
    flag = 0;                    % a valid solution has been obtained
    relres = 0;                  % the relative residual is actually 0/0
    iter = [0 0];                % no iterations need be performed
    resvec = 0;                  % resvec(1) = norm(b-A*x) = norm(0)
    if (nargout < 2)
        itermsg('gmres',tol,maxit,0,flag,iter,NaN);
    end
    return
end

if ((nargin >= 6) && ~isempty(M1))
    existM1 = 1;
    [m1type,m1fun,m1fcnstr] = iterchk(M1);  %Check the size of Precondioner.
    if strcmp(m1type,'matrix')
        if ~isequal(size(M1),[m,m])
            error(message('MATLAB:gmres:PreConditioner1Size', m));
        end
    end
else
    existM1 = 0;
    m1type = 'matrix';
end

if ((nargin >= 7) && ~isempty(M2))
    existM2 = 1;
    [m2type,m2fun,m2fcnstr] = iterchk(M2);
    if strcmp(m2type,'matrix')
        if ~isequal(size(M2),[m,m])
            error(message('MATLAB:gmres:PreConditioner2Size', m));
        end
    end
else
    existM2 = 0;
    m2type = 'matrix';
end

if ((nargin >= 8) && ~isempty(x))
    if ~isequal(size(x),[n,1])
        error(message('MATLAB:gmres:XoSize', n));
    end
else
    x = zeros(n,1);
end

if ((nargin > 8) && strcmp(atype,'matrix') && ...
        strcmp(m1type,'matrix') && strcmp(m2type,'matrix'))
    error(message('MATLAB:gmres:TooManyInputs'));
end

% Set up for the method
flag = 1;
xmin = x;                        % Iterate which has minimal residual so far
imin = 0;                        % "Outer" iteration at which xmin was computed
jmin = 0;                        % "Inner" iteration at which xmin was computed
tolb = tol * n2b;                % Relative tolerance
evalxm = 0;
stag = 0;
moresteps = 0;
maxmsteps = min([floor(n/50),5,n-maxit]);
maxstagsteps = 3;
minupdated = 0;

x0iszero = (norm(x) == 0);
r = b - iterapp('mtimes',afun,atype,afcnstr,x,varargin{:});
normr = norm(r);                 % Norm of initial residual
if (normr <= tolb)               % Initial guess is a good enough solution
    flag = 0;
    relres = normr / n2b;
    iter = [0 0];
    resvec = normr;
    if (nargout < 2)
        itermsg('gmres',tol,maxit,[0 0],flag,iter,relres);
    end
    return
end
minv_b = b;

if existM1

%    r = P*r;
    r = r-A*(ZH*(AH\(ZH'*r)));
    r = iterapp('mldivide',m1fun,m1type,m1fcnstr,r,varargin{:});
    if ~x0iszero

        b0 = b-A*(ZH*(AH\(ZH'*b)));
        minv_b = iterapp('mldivide',m1fun,m1type,m1fcnstr,b0,varargin{:});
    else
        minv_b = r;
    end
    if ~all(isfinite(r)) || ~all(isfinite(minv_b))
        flag = 2;
        x = xmin;
        relres = normr / n2b;
        iter = [0 0];
        resvec = normr;
        return
    end
end

if existM2
    r = iterapp('mldivide',m2fun,m2type,m2fcnstr,r,varargin{:});
    if ~x0iszero
        minv_b = iterapp('mldivide',m2fun,m2type,m2fcnstr,minv_b,varargin{:});
    else
        minv_b = r;
    end
    if ~all(isfinite(r)) || ~all(isfinite(minv_b))
        flag = 2;
        x = xmin;
        relres = normr / n2b;
        iter = [0 0];
        resvec = normr;
        return
    end
end

normr = norm(r);                 % norm of the preconditioned residual
n2minv_b = norm(minv_b);         % norm of the preconditioned rhs
clear minv_b;
tolb = tol * n2minv_b;
if (normr <= tolb)               % Initial guess is a good enough solution
    flag = 0;
    relres = normr / n2minv_b;
    iter = [0 0];
    resvec = n2minv_b;
    if (nargout < 2)
        itermsg('gmres',tol,maxit,[0 0],flag,iter,relres);
    end
    return
end

resvec = zeros(inner*outer+1,1);  % Preallocate vector for norm of residuals
resvec(1) = normr;                % resvec(1) = norm(b-A*x0)
normrmin = normr;                 % Norm of residual from xmin

%  Preallocate J to hold the Given's rotation constants.
J = sparse(2,inner);

U = sparse(n,inner);
R = sparse(inner,inner);
w = sparse(inner+1,1);

for outiter = 1 : outer
    %  Construct u for Householder reflector.
    %  u = r + sign(r(1))*||r||*e1
    u = r;
    normr = norm(r);
    beta = scalarsign(r(1))*normr;
    u(1) = u(1) + beta;
    u = u / norm(u);

    U(:,1) = u;

    %  Apply Householder projection to r.
    %  w = r - 2*u*u'*r;
    w(1) = -beta;

    for initer = 1 : inner
        %  Form P1*P2*P3...Pj*ej.
        %  v = Pj*ej = ej - 2*u*u'*ej
        v = -2*(u(initer)')*u;
        v(initer) = v(initer) + 1;
        %  v = P1*P2*...Pjm1*(Pj*ej)
        for k = (initer-1):-1:1
            v = v - U(:,k)*(2*(U(:,k)'*v));
        end
        %  Explicitly normalize v to reduce the effects of round-off.
        v = v/norm(v);

        %  Apply A to v.
        v = iterapp('mtimes',afun,atype,afcnstr,v,varargin{:});
        %  Apply Preconditioner.
        if existM1

            v = v-A*(ZH*(AH\(ZH'*v)));
            v = iterapp('mldivide',m1fun,m1type,m1fcnstr,v,varargin{:});
            %v = P*v;
            if ~all(isfinite(v))
                flag = 2;
                break
            end
        end

        if existM2
            v = iterapp('mldivide',m2fun,m2type,m2fcnstr,v,varargin{:});
            if ~all(isfinite(v))
                flag = 2;
                break
            end
        end
        %  Form Pj*Pj-1*...P1*Av.
        for k = 1:initer
            v = v - U(:,k)*(2*(U(:,k)'*v));
        end

        %  Determine Pj+1.
        if (initer ~= length(v))
            %  Construct u for Householder reflector Pj+1.
            u = [zeros(initer,1); v(initer+1:end)];
            alpha = norm(u);
            if (alpha ~= 0)
                alpha = scalarsign(v(initer+1))*alpha;
                %  u = v(initer+1:end) +
                %        sign(v(initer+1))*||v(initer+1:end)||*e_{initer+1)
                u(initer+1) = u(initer+1) + alpha;
                u = u / norm(u);
                U(:,initer+1) = u;

                %  Apply Pj+1 to v.
                %  v = v - 2*u*(u'*v);
                v(initer+2:end) = 0;
                v(initer+1) = -alpha;
            end
        end

        %  Apply Given's rotations to the newly formed v.
        for colJ = 1:initer-1
            tmpv = v(colJ);
            v(colJ)   = conj(J(1,colJ))*v(colJ) + conj(J(2,colJ))*v(colJ+1);
            v(colJ+1) = -J(2,colJ)*tmpv + J(1,colJ)*v(colJ+1);
        end

        %  Compute Given's rotation Jm.
        if ~(initer==length(v))
            rho = norm(v(initer:initer+1));
            J(:,initer) = v(initer:initer+1)./rho;
            w(initer+1) = -J(2,initer).*w(initer);
            w(initer) = conj(J(1,initer)).*w(initer);
            v(initer) = rho;
            v(initer+1) = 0;
        end

        R(:,initer) = v(1:inner);

        normr = abs(w(initer+1));
        resvec((outiter-1)*inner+initer+1) = normr;
        normr_act = normr;

        if (normr <= tolb || stag >= maxstagsteps || moresteps)
            if evalxm == 0
                ytmp = R(1:initer,1:initer) \ w(1:initer);
                additive = U(:,initer)*(-2*ytmp(initer)*conj(U(initer,initer)));
                additive(initer) = additive(initer) + ytmp(initer);
                for k = initer-1 : -1 : 1
                    additive(k) = additive(k) + ytmp(k);
                    additive = additive - U(:,k)*(2*(U(:,k)'*additive));
                end
                if norm(additive) < eps*norm(x)
                    stag = stag + 1;
                else
                    stag = 0;
                end
                xm = x + additive;
                evalxm = 1;
            elseif evalxm == 1
                addvc = [-(R(1:initer-1,1:initer-1)\R(1:initer-1,initer))*...
                    (w(initer)/R(initer,initer)); w(initer)/R(initer,initer)];
                if norm(addvc) < eps*norm(xm)
                    stag = stag + 1;
                else
                    stag = 0;
                end
                additive = U(:,initer)*(-2*addvc(initer)*conj(U(initer,initer)));
                additive(initer) = additive(initer) + addvc(initer);
                for k = initer-1 : -1 : 1
                    additive(k) = additive(k) + addvc(k);
                    additive = additive - U(:,k)*(2*(U(:,k)'*additive));
                end
                xm = xm + additive;
            end
            r = b - iterapp('mtimes',afun,atype,afcnstr,xm,varargin{:});
            if norm(r) <= tol*n2b
                x = xm;
                flag = 0;
                iter = [outiter, initer];
                break
            end
            minv_r = r;
            if existM1

                %minv_r = P*minv_r;
                minv_r = r-A*(ZH*(AH\(ZH'*r)));
                minv_r = iterapp('mldivide',m1fun,m1type,m1fcnstr,minv_r,varargin{:});
                if ~all(isfinite(minv_r))
                    flag = 2;
                    break
                end
            end
            if existM2
                minv_r = iterapp('mldivide',m2fun,m2type,m2fcnstr,minv_r,varargin{:});
                if ~all(isfinite(minv_r))
                    flag = 2;
                    break
                end
            end

            normr_act = norm(minv_r);
            resvec((outiter-1)*inner+initer+1) = normr_act;

            if normr_act <= normrmin
                normrmin = normr_act;
                imin = outiter;
                jmin = initer;
                xmin = xm;
                minupdated = 1;
            end

            if normr_act <= tolb
                x = xm;
                flag = 0;
                iter = [outiter, initer];
                break
            else
                if stag >= maxstagsteps && moresteps == 0
                    stag = 0;
                end
                moresteps = moresteps + 1;
                if moresteps >= maxmsteps
                    if ~warned
                        warning(message('MATLAB:gmres:tooSmallTolerance'));
                    end
                    flag = 3;
                    iter = [outiter, initer];
                    break;
                end
            end
        end

        if normr_act <= normrmin
            normrmin = normr_act;
            imin = outiter;
            jmin = initer;
            minupdated = 1;
        end

        if stag >= maxstagsteps
            flag = 3;
            break;
        end
    end         % ends inner loop

    evalxm = 0;

    if flag ~= 0
        if minupdated
            idx = jmin;
        else
            idx = initer;
        end
        y = R(1:idx,1:idx) \ w(1:idx);
        additive = U(:,idx)*(-2*y(idx)*conj(U(idx,idx)));
        additive(idx) = additive(idx) + y(idx);
        for k = idx-1 : -1 : 1
            additive(k) = additive(k) + y(k);
            additive = additive - U(:,k)*(2*(U(:,k)'*additive));
        end
        x = x + additive;
        xmin = x;
        r = b - iterapp('mtimes',afun,atype,afcnstr,x,varargin{:});
        minv_r = r;
        if existM1

            minv_r = r-A*(ZH*(AH\(ZH'*r)));
            minv_r = iterapp('mldivide',m1fun,m1type,m1fcnstr,minv_r,varargin{:});
            if ~all(isfinite(minv_r))
                flag = 2;
                break
            end
        end
        if existM2
            minv_r = iterapp('mldivide',m2fun,m2type,m2fcnstr,minv_r,varargin{:});
            if ~all(isfinite(minv_r))
                flag = 2;
                break
            end
        end
        normr_act = norm(minv_r);
        r = minv_r;
    end

    if normr_act <= normrmin
        xmin = x;
        normrmin = normr_act;
        imin = outiter;
        jmin = initer;
    end

    if flag == 3
        break;
    end
    if normr_act <= tolb
        flag = 0;
        iter = [outiter, initer];
        break;
    end
    minupdated = 0;
end         % ends outer loop

% returned solution is that with minimum residual
if flag == 0
    relres = normr_act / n2minv_b;
else
    x = xmin;
    iter = [imin jmin];
    relres = normr_act / n2minv_b;
end

% truncate the zeros from resvec
if flag <= 1 || flag == 3
    resvec = resvec(1:(outiter-1)*inner+initer+1);
    indices = resvec==0;
    resvec = resvec(~indices);
else
    if initer == 0
        resvec = resvec(1:(outiter-1)*inner+1);
    else
        resvec = resvec(1:(outiter-1)*inner+initer);
    end
end

% only display a message if the output flag is not used
% if nargout < 2
%     if restarted
%         itermsg(sprintf('gmres(%d)',restart),tol,maxit,[outiter initer],flag,iter,relres);
%     else
%         itermsg(sprintf('gmres'),tol,maxit,initer,flag,iter(2),relres);
%     end
% end

function sgn = scalarsign(d)
sgn = sign(d);
if (sgn == 0)
    sgn = 1;
end

function y = iterapp(op,afun,atype,afcnstr,x,varargin)
%ITERAPP   Apply matrix operator to vector and error gracefully.
%   ITERAPP(OP,AFUN,ATYPE,AFCNSTR,X) applies matrix operator AFUN to vector
%   X. If ATYPE is 'matrix, then AFUN is a matrix and the OP is applied
%   directly. OP is either 'mtimes' or 'mldivide'.
%   ATYPE and AFCNSTR are used in case of error.
%   ITERAPP(OP,AFUN,ATYPE,AFCNSTR,X,P1,P2,...) allows extra arguments to
%   AFUN(X,P1,P2,...) although this usage is now discouraged in favor of
%   using anonymous functions.
%   AFUN(X,P1,P2,...,PN,TFLAG) should accept a TFLAG as its final input if
%   the calling function is BICG, LSQR or QMR. TFLAG is either 'transp' or
%   'notransp' depending on whether A' OP X or A OP X is required.
%   ITERAPP is designed for use by iterative methods like PCG which
%   require matrix operators AFUN representing matrices A to operate on
%   vectors X and return A*X and may also have operators MFUN representing
%   preconditioning matrices M operate on vectors X and return M\X.
%
%   See also BICG, BICGSTAB, CGS, GMRES, LSQR, MINRES, PCG, QMR, SYMMLQ.

%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.7.4.8 $ $Date: 2011/05/17 02:33:16 $

if strcmp(atype,'matrix')
    switch lower(op)
        case 'mtimes'
            if (nargin >= 6) && isequal(varargin{end},'transp')
                y = afun' * x;
            else
                y = afun * x;
            end
        case 'mldivide'
            if (nargin >= 6) && isequal(varargin{end},'transp')
                y = afun' \ x;
            else
                y = afun \ x;
            end
        otherwise
            error(message('MATLAB:iterapp:InvalidOp'))
    end
else
    try
        if (nargin >= 6) && isequal(varargin{end},'notransp')
            % A request for A*x coming from BICG, LSQR and QMR
            try
                % New syntax: we now request afun(x,P1,P2,...,PN,'notransp')
                y = afun(x,varargin{:});
            catch
                % Old syntax: we used to accept afun(x,P1,P2,...,PN)
                y = afun(x,varargin{1:end-1});
            end
        else
            % A request for A*x
            % coming from BICGSTAB, CGS, GMRES, MINRES, PCG or SYMMLQ
            % with the call always afun(P1,P2,...,PN)
            % or a request for A'*x coming from
            % BICG, LSQR and QMR in the afun(x,P1,P2,...,PN,'transp') case
            y = afun(x,varargin{:});
        end
    catch ME
        error(message('MATLAB:iterapp:InvalidInput', atype,afcnstr, ME.message));
    end

    if ~iscolumn(y)
        error(message('MATLAB:iterapp:MustReturnColumn', atype, afcnstr));
    end
end

function [atype,afun,afcnstr] = iterchk(A)
%ITERCHK  Checks arguments to iterative methods.
%   [ATYPE,AFUN,AFCNSTR] = ITERCHK(A) returns the following:
%   ATYPE is either 'matrix', 'function', 'expression' or 'inline object'.
%   AFUN is the function name or inline object.
%   AFUN is '' if ATYPE is 'matrix'.
%   AFCNSTR is the function name if ATYPE is 'function'.
%   AFCNSTR is the formula of the function if ATYPE is 'expression' or
%   'inline object'.  AFCNSTR is '' if ATYPE is 'matrix'.
%
%   See also BICG, BICGSTAB, CGS, GMRES, LSQR, MINRES, PCG, QMR, SYMMLQ.

%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 1.8.4.3 $ $Date: 2010/08/23 23:12:57 $

[afun,afunmsg] = fcnchk(A);
if isempty(afunmsg)
   if isa(afun,'inline')
      if isa(A,'inline')
         atype = 'inline object';
      else
         atype = 'expression';
      end
      afcnstr = formula(afun);
   else % both function_handles @fun and function names 'fun'
      atype = 'function';
      if isa(A,'function_handle')
          afcnstr = func2str(A);
      else
          afcnstr = A;
      end
   end
elseif isa(A,'float')
   afun = A;
   atype = 'matrix';
   afcnstr = '';
else
   error(message('MATLAB:iterchk:InvalidInput'));
end


text:
function varargout = pm_restore_ramp(varargin)
%
% Restores linear phase-ramps in the x-, y- and z-direction
% that has previously been removed from pm by pm_estimate_ramp.
% FORMAT: pm = pm_estimate_ramp(pm,mask,ramps)
%
% Input:
% pm     : 2 or 3D phasemap that has been unwrapped and
%          that has had its ramps removed by pm_remove_ramp
% mask   : Mask that indicates which voxels are worth
%          bothering with and which are not.
% ramps  : 3x1 vector signifying the slope of the ramps in
%          the x-, y- and z-directions. This SHOULD be the
%          values returned by a previous call to pm_estimate_ramp.
%
% Output:
% pm     : Same as pm in, but with linear ramps restored.
%
% This routine was written on the suggestion of Mark J, and will
% potentially improve performance of subsequent phase-unwrapping.
% I haven't actually found it particularly helpful, and it may
% simply have been a sneaky fMRIB attempt to delay the SPM
% phasemap toolbox.
%__________________________________________________________________
% Jesper Andersson 30/9-03
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: pm_restore_ramp.m 1317 2008-04-08 16:16:38Z chloe $

error('mex-function pm_restore_ramp.c not compiled');

author: spm


text:
% Gianni Schena  July 2005, schena@units.it
% Lattice Boltzmann LBE, geometry: D2Q9, model: BGK
% Application to permeability in porous media

Restart=false % to restart from an earlier convergence
logical(Restart);

if Restart==false;
close all, clear all % start from scratch and clean ...
Restart=false;
% type of channel geometry ;
% one of the flollowing flags == true
Pois_test=true, % no obstacles in the 2D channel
% porous systems
obs_regolare=false %
obs_irregolare=false %
tic
%   IN
% |vvvv|    + y
% |vvvv|     ^
% |vvvv|     | -> + x
%  OUT

% Pores in 2D : Wet and Dry locations (Wet ==1 , Dry ==0 )
wXh_Dry=[3,1];wXh_Wet=[3,4];

if obs_regolare, % with internal obstacles

A=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]);A=[A,zeros(wXh_Dry)];
B=ones(size(A));
C=[A;B]  ; D=repmat(C,4,1);
D=[B;D]
end

if obs_irregolare, % with int obstacles
A1=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]);
A1=[A1,zeros(wXh_Dry)]  ;
B=ones(size(A1));
C1=repmat([ones(wXh_Wet),zeros(wXh_Dry)],[1,3]); C1=[C1,ones(wXh_Dry)];
E=[A1;B;C1;B];
D=repmat(E,2,1);
D=[B;D]
end

if ~Pois_test
figure,imshow(D,[])
Channel2D=D;
Len_Channel_2D=size(Channel2D,1); % Length
Width=size(Channel2D,2); % should not be hod
Channel_2D_half_Width=Width/2,
end

% test without obstacles (i.e. 2D channel & no obstacles)

if Pois_test
%over-writes the definition of the pore space
clear Channel2D
Len_Channel_2D=36, % lunghezza canale 2d
Channel_2D_half_Width=8; Width=Channel_2D_half_Width*2;
Channel2D=ones(Len_Channel_2D,Width); % define wet area
%Channel2D(6:12,6:8)=0; % put fluid obstacle
imshow(Channel2D,[]);
end

[Nr Mc]=size(Channel2D); % Number rows and Munber columns

% porosity
porosity=nnz(Channel2D==1)/(Nr*Mc)

% FLUID PROPERTIES
% physical properties
cs2=1/3; %
cP_visco=0.5; % [cP] 1 CP Dinamic water viscosity 20 C
density=1.; % fluid density
Lky_visco=cP_visco/density; % lattice kinematic viscosity
omega=(Lky_visco/cs2+0.5).^-1; %  omega: relaxation frequency
%Lky_visco=cs2*(1/omega - 0.5) , % lattice kinematic viscosity
%dPdL= Pressure / dL;% External pressure gradient [atm/cm]

uy_fin_max=-0.2;
%dPdL = abs( 2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2) );
dPdL=-0.0125;
uy_fin_max=dPdL*(Channel_2D_half_Width.^2)/(2*Lky_visco); % Poiseuille Gradient;
% max poiseuille final  velocity on the flow profile
uy0=-0.001; ux0=0.0001; %  linear vel .. inizialization

%
% uy_fin_max=-0.2; % max poiseuille final  velocity on the flow profile
% omega=0.5, cs2=1/3; % omega: relaxation frequency
% Lky_visco=cs2*(1/omega - 0.5) , % lattice kinematic viscosity
% dPdL = abs( 2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2) ); % Poiseuille Gradient;
%

uyf_av=uy_fin_max*(2/3);; % average fluid velocity on the profile

x_profile=([-Channel_2D_half_Width:+Channel_2D_half_Width-1]+0.5);
uy_analy_profile=uy_fin_max.*(1-  ( x_profile /Channel_2D_half_Width).^2 ); % analytical velocity profile

av_vel_t=1.e+10; % inizialization (t=0)
%PixelSize= 5; % [Microns]
%dL=(Nr*PixelSize*1.0E-4); % sample hight [cm]

%
% EXPERIMENTAL SET-UP
% inlet and outlet buffers
inb=2, oub=2; % inlet and outlet buffers thickness
% add fluid at the inlet (top) and outlet (down)
inlet=ones(inb,Mc); outlet=ones(oub,Mc);
Channel2D=[ [inlet]; Channel2D ;[outlet] ] ; % add flux in and down (E to W)
[Nr Mc]=size(Channel2D); % update size
% boundaries related to the experimental set up
wb=2; % wall thickness
Channel2D=[zeros(Nr,wb), Channel2D , zeros(Nr,wb)]; % add walls (no fluid leak)
[Nr Mc]=size(Channel2D); % update size
uy_analy_profile=[zeros(1,wb), uy_analy_profile, zeros(1,wb) ] ; % take into account walls
x_pro_fig=[[x_profile(1)-[wb:-1:1]], [x_profile, [1:wb]+x_profile(end)] ];

% Figure plots analytical parabolic profile
figure(20), plot(x_pro_fig,uy_analy_profile,'-'), grid on,
title('Analytical parab. profile for Poiseuille planar flow in a channel')

% VISUALIZE PORE SPACE & FLUID OSTACLES & MEDIAL AXIS
figure, imshow(Channel2D); title('Vassel geometry');
Channel2D=logical(Channel2D);
% obstacles for Bounce Back ( in front of the grain)
Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
border=logical(ones(Nr,Mc));
border([1:inb,Nr-oub:Nr],[wb+2:Mc-wb-1])=0;
Obstacles=Obstacles.*(border);
figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );
%
Medial_axis=bwmorph(Channel2D,'thin',Inf); %
figure, imshow(Medial_axis); title('Medial axis');
figure(10) % used to visualize evolution of rho
figure(11) % used to visualize ux
figure(12) % used to visualize uy (i.e. top -> down)

% INDICES
% Wet locations etc.
[iabw1 jabw1]=find(Channel2D==1); % indices i,j, of active lattice locations i.e. pore
lena=length(iabw1); % number of active location i.e. of pore space lattice cells
ija= (jabw1-1)*Nr+iabw1; % equivalent single index (i,j)->> ija for active locations
% absolute (single index) position of the obstacles in for bounce back in Channel2D
% Obstacles
[iobs jobs]=find(Obstacles);lenobs=length(iobs); ijobs= (jobs-1)*Nr+iobs; % as above
% Medial axis of the pore space
[ima jma]=find(Medial_axis); lenma=length(ima);  ijma= (jma-1)*Nr+ima; % as above
% Internal wet locations : wet & ~obstables
% (i.e. internal wet lattice location non in contact with dray locations)
[iawint jawint]=find(( Channel2D==1 & ~Obstacles)); % indices i,j, of active lattice locations
lenwint=length(iawint); % number of internal (i.e. not border) wet locations
ijaint= (jawint-1)*Nr+iawint; % equivalent singl
NxM=Nr*Mc;

% DIRECTIONS: E N W S NE NW SW SE ZERO (ZERO:Rest Particle)
%    y^
%  6 2 5           ^         NW  N  NE
%  3 9 1 ... +x-> +y         W   RP  E
%  7 4 8                     SW  S  SE
%   -y
% x & y components of velocities , +x is to est , +y is to nord
East=1; North=2; West=3; South=4; NE=5; NW=6; SW=7; SE=8; RP=9;
N_c=9 ; % number of directions
% versors D2Q9
C_x=[1 0 -1  0 1 -1 -1  1 0];
C_y=[0 1  0 -1 1  1 -1 -1 0]; C=[C_x;C_y]

% BOUNCE BACK SCHEME
% after collision the fluid elements densities f are sent back to the
% lattice node they come from with opposite direction
% indices opposite to 1:8 for fast inversion after bounce
ic_op = [3 4 1 2 7 8 5 6]; %   i.e. 4 is opposite to 2 etc.

% PERIODIC BOUNDARY CONDITIONS - reinjection rules
yi2=[Nr , 1:Nr , 1]; % this definition allows implemening Period Bound Cond
%yi2=[1, Nr , 2:Nr-1 , 1,Nr]; % re-inj the second last to as first
% directional weights (density weights)
w0=16/36. ; w1=4/36. ; w2=1/36.;
W=[ w1 w1 w1 w1 w2 w2 w2 w2 w0];
%c constants (sound speed related)
cs2=1/3; cs2x2=2*cs2; cs4x2=2*cs2.^2;
f1=1/cs2; f2=1/cs2x2; f3=1/cs4x2;
f1=3., f2=4.5; f3=1.5; % coef. of the f equil.

% declarative statemets
f=zeros(Nr,Mc,N_c); % array of fluid density distribution
feq=zeros(Nr,Mc,N_c); % f at equilibrium
rho=ones(Nr,Mc); % macro-scopic density
temp1=zeros(Nr,Mc);
ux=zeros(Nr,Mc);   uy=zeros(Nr,Mc); uyout=zeros(Nr,Mc);  % dimensionless velocities
uxsq=zeros(Nr,Mc); uysq=zeros(Nr,Mc);   usq=zeros(Nr,Mc);  % higher degree velocities

% initialization arrays : start values in the wet area
for ia=1:lena % stat values in the active cells only ; 0 outside
    i=iabw1(ia);  j=jabw1(ia);
    f(i,j,:)=1/9; % uniform density distribution for a start
end
uy(ija)=uy0; ux(ija)=ux0; % initialize fluid velocities
rho(ija)=density;

% EXTERNAL (Body) FORCES e.g. inlet pressure or inlet-outlet gradient
% directions: E N W S NE NW SW SE ZERO
force = -dPdL*(1/6)*1*[0 -1 0 1 -1 -1 1  1  0]'; %;
%...                   E  N E S NE NW SW SE RP ...
% the pressure pushes the fluid down i.e. N to S

% While .. MAIN TIME EVOLUTION LOOP
StopFlag=false; % i.e. logical(0)
Max_Iter=3000; % max allowed number of iteration
Check_Iter=1; Output_Every=20; % frequency of check & output
Cur_Iter=0; % current iteration counter inizialization
toler=1.0e-8; % tollerance to declare convegence
Cond_path=[]; % recording values of the convergence criterium
density_path=[]; % recording aver. density values for convergence
end % ends if restart

if(Restart==true)
 StopFlag=false;  Max_Iter=Max_Iter+3000; toler=1.0e-12;
end

while(~StopFlag)
    Cur_Iter=Cur_Iter+1 % iteration counter update

    % density and moments
    rho=sum(f,3); % density

    if Cur_Iter >1 % use inizialization ux uy to start
        % Moments ... Note:C_x(9)=C_y(9)=0
        ux=zeros(Nr,Mc); uy=zeros(Nr,Mc);
        for ic=1:N_c-1;
            ux = ux + C_x(ic).*f(:,:,ic) ; uy = uy + C_y(ic).*f(:,:,ic)  ;
        end
       % uy=f(:,:,2) +f(:,:,5)+f(:,:,6)-f(:,:,4)-f(:,:,7)-f(:,:,8); % in short !
       % ux=f(:,:,1) +f(:,:,5)+f(:,:,8)-f(:,:,3)-f(:,:,6)-f(:,:,7); % in short !
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ux(ija)=ux(ija)./rho(ija); uy(ija)=uy(ija)./rho(ija);
    uxsq(ija)=ux(ija).^2; uysq(ija)=uy(ija).^2;
    usq(ija)=uxsq(ija)+uysq(ija); %

    % weighted densities : rest particle, principal axis, diagonals
    rt0 = w0.*rho; rt1 = w1.*rho; rt2 = w2.*rho;

    % Equilibrium distribution
    % main  directions ( + cross)
    feq(ija)= rt1(ija) .*(1 +f1*ux(ija) +f2*uxsq(ija) -f3*usq(ija));
    feq(ija+NxM*(2-1))= rt1(ija) .*(1 +f1*uy(ija) +f2*uysq(ija) -f3*usq(ija));
    feq(ija+NxM*(3-1))= rt1(ija) .*(1 -f1*ux(ija) +f2*uxsq(ija) -f3*usq(ija));
    %feq(ija+NxM*(3)=f(ija)-2*rt1(ija)*f1.*ux(ija); % much faster... !!
    feq(ija+NxM*(4-1))= rt1(ija) .*(1 -f1*uy(ija) +f2*uysq(ija) -f3*usq(ija));

    % diagonals (X diagonals) (ic-1)
    feq(ija+NxM*(5-1))= rt2(ija) .*(1 +f1*(+ux(ija)+uy(ija)) +f2*(+ux(ija)+uy(ija)).^2 -f3.*usq(ija));
    feq(ija+NxM*(6-1))= rt2(ija) .*(1 +f1*(-ux(ija)+uy(ija)) +f2*(-ux(ija)+uy(ija)).^2 -f3.*usq(ija));
    feq(ija+NxM*(7-1))= rt2(ija) .*(1 +f1*(-ux(ija)-uy(ija)) +f2*(-ux(ija)-uy(ija)).^2 -f3.*usq(ija));
    feq(ija+NxM*(8-1))= rt2(ija) .*(1 +f1*(+ux(ija)-uy(ija)) +f2*(+ux(ija)-uy(ija)).^2 -f3.*usq(ija));
    % rest particle (.) ic=9
    feq(ija+NxM*(9-1))= rt0(ija) .*(1 - f3*usq(ija));

    %Collision (between fluid elements)omega=relaxation frequency
    f=(1.-omega).*f + omega.*feq;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %add external body force due to the pressure gradient prop. to dPdL
    for ic=1:N_c;%-1
        for ia=1:lena
            i=iabw1(ia);  j=jabw1(ia);
            % if Obstacles(i,j)==0 % the i,j is not aderent to the boundaries
            % if ( f(i,j,ic) + force(ic) ) >0; %! avoid negative distributions
            %i=1 ;% force only on the first row !
            f(i,j,ic)= f(i,j,ic) + force(ic);
            % end
            % end
        end
    end

    % % STREAM
    % Forward Propagation step & % Bounce Back (collision fluid with obstacles)
    %f(:,:,9) = f(:,:,9); % Rest element do not move

    feq = f; % temp storage of f in feq
        for ic=1:1:N_c-1, % select velocity layer

        ic2=ic_op(ic); % selects the layer of the velocity opposite to ic for BB
        temp1=feq(:,:,ic); %

        % from wet location that are NOT on the border to other wet locations
        for ia=1:1:lenwint % number of internal (i.e. not border) wet locations
            i=iawint(ia);  j=jawint(ia);  % so that we care for the wet space only !
            i2 = i+C_y(ic); j2 = j+C_x(ic); % Expected final locations to move
            i2=yi2(i2+1); % i2 corrected for PBC when necessary (flow out re-fed to inlet)
            % i.e the new position (i2,j2) is sure another wet location
            % therefore normal propagation from (i,j) to (i2,j2) on layer ic
            f(i2,j2,ic)=temp1(i,j); % see circshift(..) fnct for circularly shifts
        end ; % i and j single loop

        % from wet locations that ARE on the border of obstacles
        for ia=1:1:lenobs % wet border locations
            i=iobs(ia);  j=jobs(ia);  % so that we care for the wet space only !
            i2 = i+C_y(ic); j2 = j+C_x(ic); % Expected final locations to move
            i2=yi2(i2+1); % i2 corrected for PBC

            if( Channel2D(i2,j2) ==0 ) % i.e the new position (i2,j2) is dry
                f(i,j,ic2) =temp1(i,j); % invert direction: bounce-back in the opposite direction ic2
            else % otherwise, normal propagation from (i,j) to (i2,j2) on layer ic
                f(i2,j2,ic)=temp1(i,j); % see circshift(..) fnct for circularly shifts
            end ; % b.b. and propagations

        end ; % i and j single loop
        % special treatment for Corners
        %   f(1,wb+1,ic)=temp1(Nr,Mc-wb);      f(1,Mc-wb,ic)=temp1(Nr,wb+1);
        %   f(Nr,wb+1,ic)=temp1(1,Mc-wb);      f(Nr,Mc-wb,ic)=temp1(1,wb+1);

    end ; %  for ic direction

    % ends of Forward Propagation step &  Bounce Back Sections

    % re-calculate  uy as uyout for convergence
    rho=sum(f,3); % density
    % check velocity
    uyout= zeros(Nr,Mc);
    for ic=1:N_c-1;
        uyout= uyout + C_y(ic).*f(:,:,ic) ; % flow dim.less velocity out
    end
   % uyout(ija)=uyout(ija)./rho(ija); % from momentum to velocity

    % Convergence check on velocity values
    if (mod(Cur_Iter,Check_Iter)==0) ; % check for convergence every 'Check_Iter' iterations

        % variables monitored
        % mean density and
        vect=rho(ija); vect=vect(:);
        cur_density=mean(vect);
        % mean 'interstitial' velocity
        % uy(ija)=uy(ija)/rho(ija); ?
        vect=uy(ija); av_vel_int= mean(vect)  ; % seepage velocity (in the wet area)
        % on the whole cross-sectional area of flow (wet + dry)
        av_vel_int=av_vel_int*porosity, % av. vel. on the wet + dry area
        %av_vel_int=mean2(uy),
        av_vel_tp1 = av_vel_int;
        Condition=abs( abs(av_vel_t/av_vel_tp1 )-1), % should --> 0

        Cond_path=[Cond_path, Condition]; % records the convergence path (value)
        density_path=[density_path, cur_density];
        %
        av_vel_t=av_vel_tp1; % time t & t+1

        if (Condition < toler) | (Cur_Iter > Max_Iter)
            StopFlag=true;
            display( 'Stop iteration: Convergence met or iteration exeeding the max allowed' )
            display( ['Current iteration: ',num2str(Cur_Iter),...
                ' Max Number of iter: ',num2str(Max_Iter)] )
            break % Terminate execution of WHILE .. exit the time evolution loop.

        end    % if(Condition < toler

    end

    if (mod(Cur_Iter,Output_Every)==0) ;  % Output from loop every ...
        %if (Cur_Iter>60) ;  % Output from loop every ...

        rho=sum(f,3); % density
        figure(10); imshow(rho,[0.1 0.9]); title(' rho'); % visualize density evolution
        figure(11); imshow(ux,[ ]); title(' ux' ); % visualize fluid velocity horizontal
        figure(12); imshow(-uy,[ ]); title(' uy' ); % visualize fluid velocity down
        figure(14), imshow(-uyout,[]), title('uyout'); % vis vel flow out
        up=2; % linear section to visualize up from the lower row
        figure(15), hold off, feather(ux(Nr-up,:),uy(Nr-up,:)),
        figure(15), hold on , plot(uy_analy_profile,'r-')
        title('Analytical vs LB calculated, fluid velocity parabolic profile')
        pause(3); % time given to visualize properly

    end % every

   % pause(1);

end %  End main time Evolution Loop

% Output & Draw after the end of the time evolution

figure, plot(Cond_path(2:end)); title('convergence path')
%figure, plot(density_path(2:end)); title('density convergence path')
figure, plot( [uy(Nr-up,:)-uy_analy_profile] ); title('difference : LB - Analytical solution')

toc

% Permeability K

K_Darcy_Porous_Sys= (av_vel_int*porosity)/dPdL*Lky_visco ,

K_Analy_2D_Channel=(Width^2)/12


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Forward FFT w.r.t. the second variable %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fs=fty(s);
 fs=fftshift(fft(fftshift(s.'))).';


text:
if ~exist('calib_name')|~exist('format_image'),
   data_calib_no_read;
   return;
end;

check_directory;

if ~exist('n_ima'),
   data_calib_no_read;
   return;
end;

check_active_images;

images_read = active_images;

if exist('image_numbers'),
   first_num = image_numbers(1);
end;

% Just to fix a minor bug:
if ~exist('first_num'),
   first_num = image_numbers(1);
end;

image_numbers = first_num:n_ima-1+first_num;

no_image_file = 0;

%% Step used to clean the memory if a previous atttempt has been made to read the entire set of images into memory:
for kk = 1:n_ima,
    if (exist(['I_' num2str(kk)])==1),
        clear(['I_' num2str(kk)]);
    end;
end;

fprintf(1,'\nChecking directory content for the calibration images (no global image loading in memory efficient mode)\n');

one_image_read = 0;

i = 1;

while (i <= n_ima), % & (~no_image_file),

    if active_images(i),

        %fprintf(1,'Loading image %d...\n',i);

        if ~type_numbering,
            number_ext =  num2str(image_numbers(i));
        else
            number_ext = sprintf(['%.' num2str(N_slots) 'd'],image_numbers(i));
        end;

        ima_name = [calib_name  number_ext '.' format_image];

        if i == ind_active(1),
            fprintf(1,'Found images: ');
        end;

        if exist(ima_name),

            fprintf(1,'%d...',i);

            if ~one_image_read

                if format_image(1) == 'p',
                    if format_image(2) == 'p',
                        I = double(loadppm(ima_name));
                    else
                        I = double(loadpgm(ima_name));
                    end;
                else
                    if format_image(1) == 'r',
                        I = readras(ima_name);
                    else
                        I = double(imread(ima_name));
                    end;
                end;

                if size(I,3)>1,
                    I = 0.299 * I(:,:,1) + 0.5870 * I(:,:,2) + 0.114 * I(:,:,3);
                end;

                if size(I,1)~=480,
                    small_calib_image = 1;
                else
                    small_calib_image = 0;
                end;

                [Hcal,Wcal] = size(I); 	% size of the calibration image

                [ny,nx] = size(I);

                one_image_read = 1;

            end;

        else

            images_read(i) = 0;

        end;

    end;

    i = i+1;

end;

ind_read = find(images_read);

if ~(exist('map')==1), map = gray(256); end;

active_images = images_read;

fprintf(1,'\ndone\n');

repo: ekfmonoslam


text:
function varargout = metabolites_in_database(varargin)
% m-file for metabolites_in_database.fig
% metabolites_in_database, by itself, creates a new object or raises the existing singleton*.
%
% USAGE:
%
%    varargout = metabolites_in_database(varargin)
%
% INPUTS:
%    varargin:    various input arguments
%
% OUTPUTS:
%    varargout:   various output arguments
%
% EXAMPLE:
%
%    H = metabolites_in_database() returns the handle to a new metabolites_in_database or the handle to
%    the existing singleton*.
%
%    metabolites_in_database('CALLBACK',hObject,eventData,handles,...) calls the local
%    function named CALLBACK in metabolites_in_database.M with the given input arguments.
%
%    metabolites_in_database('Property','Value',...) creates a new metabolites_in_database or raises the
%    existing singleton*.  Starting from the left, property value pairs are
%    applied to the GUI before `metabolites_in_database_OpeningFcn` gets called.  An
%    unrecognized property name or invalid value makes property application
%    stop.  All inputs are passed to `metabolites_in_database_OpeningFcn` via `varargin`.
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% NOTE:
%    See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%    instance to run (singleton)".
%    See also: GUIDE, GUIDATA, GUIHANDLES
% .. Edit the above text to modify the response to help metabolites_in_database
% .. Last Modified by GUIDE v2.5 31-Aug-2010 18:09:02
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is
gui_Singleton = 1; % Begin initialization code - DO NOT EDIT
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @metabolites_in_database_OpeningFcn, ...
                   'gui_OutputFcn',  @metabolites_in_database_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

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

% --- Executes just before metabolites_in_database is made visible.
function metabolites_in_database_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to metabolites_in_database (see VARARGIN)
data = varargin{1};
set(handles.uitable1,'data',data);
% Choose default command line output for metabolites_in_database
handles.output = 0;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes metabolites_in_database wait for user response (see UIRESUME)
 uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = metabolites_in_database_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;
delete(gcf)

% --- Executes on button press in act_continue.
function act_continue_Callback(hObject, eventdata, handles)
% hObject    handle to act_continue (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.output = 1;
guidata(hObject,handles)
uiresume;

% --- Executes on button press in act_cancel.
function act_cancel_Callback(hObject, eventdata, handles)
% hObject    handle to act_cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
uiresume;

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
uiresume

author: opencobra

repo: cobratoolbox


text:
function lam = line_loop_adj_eigenvalues ( n )

%*****************************************************************************80
%
%% LINE_LOOP_ADJ_EIGENVALUES: the eigenvalues of the LINE_LOOP_ADJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAM(N,1), the eigenvalues.
%
  lam = zeros ( n, 1 );

  for i = 1 : n
    angle = i * pi / ( n + 1 );
    lam(i) = 1.0 + 2.0 * cos ( angle );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [G,model] = penlab_callback_dcon2(x,k,model)

global latest_x_g
global latest_G
global latest_g
x = x(:);

G = spalloc(length(model.c),length(model.c),0);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function handles = select_plane_from_point_cloud(handles,preview_filenames,ind_preview)

fig = figure;
dcm_obj = datacursormode(fig);

h = uicontrol('Position',[20 20 100 40],'String','Continue',...
              'Callback',@continue_callback);
h1 = uicontrol('Position',[220 20 100 40],'String','back',...
              'Enable','off','Callback',@back_callback);

for i=1:numel(preview_filenames)
    disp(preview_filenames{i});
    done = 0;
    pc = import_point_cloud_from_file(preview_filenames{i});
    pcshow(pc);
    title(preview_filenames{i}, 'interpreter', 'none');

    uiwait(gcf);

end
close(gcf);

function continue_callback(objectHandle , eventData )
    if done == 0
        cur_info = getCursorInfo(dcm_obj);
        if isfield(cur_info,'Position')
            arrayPos=find(handles.active_image_numbers==ind_preview(i));

            [theta, alpha, minError, inliers] = ObtainPlanePoint(pc.Location,cur_info.Position);
%             [theta,alpha,minError] = tls_robust(inliers);
            handles.user_selected_planes{arrayPos}.theta = theta;
            handles.user_selected_planes{arrayPos}.alpha = alpha;
            handles.user_selected_planes{arrayPos}.e = minError;
            handles.user_selected_planes{arrayPos}.inliers = inliers;

            pcshow(pointCloud(inliers));
            h1.Enable = 'on';
            done = 1;
        end
    else
        h1.Enable = 'off';
        uiresume(gcbf);
    end

end
function back_callback(objectHandle , eventData )

    pcshow(pc);
    done = 0;
    objectHandle.Enable = 'off';
end

end

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
% KM_DEMO_KDE Kernel density estimation demo.
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2013
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

close all
clear

%% PARAMETERS

N1 = 100; % number of data in data set 1
m1 = -1; % mean value
s1 = 0.1; % variance

N2 = 500; % number of data in data set 2
m2 = 2; % mean value
s2 = 0.5; % variance

sigma = 0.1; % bandwidth
nsteps = 100; % number of abscis points in kde

%% PROGRAM
tic

x1 = m1 + sqrt(s1)*randn(N1,1);
x2 = m2 + sqrt(s2)*randn(N2,1);
x = [x1;x2];

[xx,pp] = km_kde(x,sigma,nsteps,1.2);

toc
%% OUTPUT

figure;

[n,h] = hist(x,100);
bar(h,n/sum(n),'b','EdgeColor','b')

hold on
plot(xx,pp,'r','LineWidth',3);

legend('histogram','KDE')

author: epfl-lasa


text:
function N = tangentspacefactory(M, x)
% Returns a manifold structure representing the tangent space to M at x.
%
% N = tangentspacefactory(M, x)
%
% N defines a (linear) manifold that is the tangent space to M at x. Points
% are represented as tangent vectors to M at x. Tangent vectors are also
% represented as tangent vectors to M at x.
%
% This is chiefly useful to solve optimization problems involving tangent
% vectors to M at x, which notably comes up when solving linear systems
% involving, for example, the Hessian of the cost on M at x (think of the
% Newton equations.) The Riemannian (actually, Euclidean) structure on N is
% that of the tangent space to M, that is, the inner product is inherited.
%
% See also: preconhessiansolve

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 9, 2015.
% Contributors:
% Change log:
%
%   Jan. 25, 2017 (NB):
%       Following a comment by Jesus Briales on the Manopt forum, the
%       functions N.egrad2rgrad, N.ehess2rhess and N.tangent now include a
%       projection (they were formerly identities.)
%
%   Feb. 2, 2017 (NB):
%       Following a comment by Jesus Briales on the Manopt forum, the
%       function N.proj now calls M.proj(x, .) instead of M.proj(y, .).
%       Furthermore, N.ehess2rhess was corrected in the same way.
%
%   Dec. 14, 2019 (NB):
%       Fixed N.tangent so that it should now work with factories that
%       have a non-identity tangent2ambient, e.g., fixedrankembeddedfactory
%       and rotationsfactory.

    % N is the manifold we build. y will be a point on N, thus also a
    % tangent vector to M at x. This is a typical Euclidean space, hence it
    % will be easy to describe in terms of the tools available for M.
    N = struct();

    % u, u1 and u2 will be tangent vectors to N at y. The tangent space to
    % N at y is the tangent space to M at x, thus u, u1 and u2 are also
    % tangent vectors to M at x.

    if isfield(M, 'name')
        N.name  = @() ['Tangent space to ' M.name()];
    end
    N.dim   = @() M.dim();
    N.inner = @(y, u1, u2) M.inner(x, u1, u2);
    N.norm  = @(y, u) M.norm(x, u);
    N.proj  = @(y, u) M.proj(x, u);
    N.typicaldist = @() sqrt(N.dim());
    if isfield(M, 'tangent2ambient')
        N.tangent = @(y, u) M.proj(x, M.tangent2ambient(x, u));
    else
        N.tangent = N.proj;
    end

    N.egrad2rgrad = N.proj;
    N.ehess2rhess = @(y, eg, eh, d) M.proj(x, eh);
    N.exp = @exponential;
    N.retr = @exponential;
    N.log = @(y1, y2) M.lincomb(x, 1, y2, -1, y1);
    N.pairmean = @(y1, y2) M.lincomb(x, 0.5, y1, 0.5, y2);
    N.rand = @() M.randvec(x);
    N.randvec = @(y) M.randvec(x);
    N.zerovec = M.zerovec;
    N.lincomb = M.lincomb;
    N.transp = @(y1, y2, u) u;
    N.hash = @(y) ['z' hashmd5(M.vec(x, y))];

    % In a Euclidean space, the exponential is merely the sum: y + tu.
    function yy = exponential(y, u, t)
        if nargin == 2
            t = 1;
        end
        yy = M.lincomb(x, 1, y, t, u);
    end

end

author: NicolasBoumal

repo: manopt


text:
function res = matching( des, varargin )
%MATCHING Run the image matching task
%  RES = MATCHING(DESCRIPTOR) Compute the results for descriptor
%  DESCRIPTOR. Descriptor can either be an object of desc.memdesc or a
%  string, in which case a descriptor from
%
%     `<HBROOT>/data/descriptors/DESCRIPTOR`
%
%  are loaded.
%  Stores the results in csv files in:
%
%      `<HBROOT>/matlab/results/default/DESCRIPTOR/matching.csv`
%
%  If the scores file already exists, skips the computation. Use
%  `'override', true` to change this behaviour.
%
%  Additionally accepts the following 'OptionName', OptionValue arguments:
%
%  'scoresroot' :: '<HBROOT>/matlab/results/default/'
%     Change for a different target path for scores.
%
%  'sequences' :: all of the present for the descriptor
%     Limit the sequences for which the scores are computed (cell array of
%     strings).
%
%  'geom_noise' :: {'easy', 'hard', 'tough'}
%     Limit the geometry noise computed.
%
%  'impairs' :: [1 2; 1 3; 1 4; 1 5; 1 6]
%     Change the image pairs for which the scores are computed.
%
%  'override' :: false
%     Do not overwrite existing score files.

% Copyright (C) 2017 Karel Lenc
% 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).
if ischar(des), [des, varargin] = desc.memdesc(des, varargin{:}); end;
opts.scoresroot = fullfile(hb_path, 'matlab', 'scores', 'default');
[opts, varargin] = vl_argparse(opts, varargin);
opts.scorespath = fullfile(opts.scoresroot, des.name, 'matching.csv');
opts.sequences = utls.listdirs(des.path);
opts.geom_noise = {'easy', 'hard', 'tough'};
opts.impairs = [1 2; 1 3; 1 4; 1 5; 1 6];
opts.override = false;
opts.verbose = false;
[opts, ~] = vl_argparse(opts, varargin);

if opts.verbose, display(opts); end;
if ~iscell(opts.sequences), opts.sequences  = {opts.sequences}; end;
if ~iscell(opts.geom_noise), opts.geom_noise  = {opts.geom_noise}; end;
vl_xmkdir(fileparts(opts.scorespath));
if ~opts.override && exist(opts.scorespath, 'file')
  res = readtable(opts.scorespath);
  fprintf('Results loaded from %s.\n', opts.scorespath);
  return;
end;

numtasks = numel(opts.sequences) * numel(opts.geom_noise) * size(opts.impairs, 1);
status = utls.textprogressbar(numtasks, 'updatestep', 1, ...
  'startmsg', sprintf('Evaluating "matching" for %s', des.name));
stepi = 1;
res = cell(numtasks, 1);
for gni = 1:numel(opts.geom_noise)
  geom_noise = opts.geom_noise{gni};
  for si = 1:numel(opts.sequences)
    sequence = opts.sequences{si};
    for pi = 1:size(opts.impairs, 1)
      pair = opts.impairs(pi, :);
      descA = des.getimdescs(des, sequence, geom_noise, pair(1));
      descB = des.getimdescs(des, sequence, geom_noise, pair(2));
      res_s = bench.matching.eval(descA, descB);
      res{stepi} = struct(...
        'descriptor', des.name, 'sequence', sequence, ...
        'geom_noise', geom_noise, 'ima', pair(1)-1, 'imb', pair(2)-1, ...
        'ap', res_s.ap, 'auc', res_s.auc, 'sr', res_s.sr);
      status(stepi); stepi = stepi+1;
    end
  end
end
res = struct2table(cell2mat(res), 'AsArray', true);
writetable(res, opts.scorespath);
fprintf('\n');

author: hpatches

repo: hpatches-benchmark


text:
function shape_scaled = scaleshape(shape, scale)
%SCALESHAPE Summary of this function goes here
%   Function: scale input shape using a scale ratio
%   Detailed explanation goes here
%   Input:
%        bbox: the bbox of current sample
%        shape: input shape
%        scale: scale ratio
%   Output:
%        shape_scaled: scaled shape

shape_scaled = bsxfun(@plus, scale*(bsxfun(@minus, shape, mean(shape))), mean(shape));

end

author: jwyang

repo: face-alignment


text:
function face = meshFace(faces, index)
%MESHFACE Return the vertex indices of a face in a mesh
%
%   FACE = meshFace(FACES, INDEX)
%   Return the vertex indices of the i-th face in the face array. This is
%   mainly an utilitary function that manages faces stored either as int
%   array (when all faces have same number of sides) or cell array (when
%   faces may have different number of edges).
%
%   Example
%   meshFace
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-10-06,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% process mesh given as structure
if isstruct(faces)
    if isfield(faces, 'faces')
        faces = faces.faces;
    else
        error('Mesh structure should contains a field ''faces''');
    end
end

% switch between numeric or cell array
if isnumeric(faces)
    face = faces(index, :);
else
    face = faces{index};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function fft_plot(samples, Fs)
    N = length(samples);
    df = Fs/N;
    f = -Fs/2:df:Fs/2-df;
    Y = fftshift(abs(fft(samples)));
    subplot(121);
    plot(f, Y);
    subplot(122);
    nfft = 128;
    spectrogram(samples, hanning(nfft), round(nfft*0.9), nfft, Fs);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
% Configuration function for some constants used in the library
%
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef Config < handle
    properties (Constant)
        EFFECTIVE_SUPPORT_THRESHOLD = double(5e-4)
%         LOG_FLOOR_VALUE = double(1e-5)
        LOG_FLOOR_VALUE = double(1e-10)
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function predictAllRF_Aerts(pathRF)

startpath = pwd;

cd(pathRF), load('training'), load('testing')
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
fSetNames = {'PETsign','CTsign','CTorigsign'}; nFset = numel(fSetNames); % The two sets tested for the radiomic signature

for o = 1:nOutcomes
    if strcmp(nameOutcomes{o},'DeathSign')
        outcome = testing.outcomes.Death;
        time = testing.timeToEvents.Death;
    else
        outcome = testing.outcomes.(nameOutcomes{o});
        time = testing.timeToEvents.(nameOutcomes{o});
    end
    censoring = 1 - outcome;
    for f = 1:nFset
        results = struct;
        RF = load(['RF_',[fSetNames{f},'Clinic'],'_',nameOutcomes{o}]); RF = struct2cell(RF); RF = RF{1};
        indClinic = training.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f});
        text = testing.textures.(nameOutcomes{o}).(fSetNames{f}); nText = size(text,2);
        tableTest = [text,testing.clinical.table(:,indClinic)];
        [prob] = predictRF(tableTest,RF);
        results.probResponse = prob;
        [AUC,sensitivity,specificity,accuracy] = calcPerformMetrics(prob,outcome,0.5);
        CI = calcCI(prob,time,censoring);
        results.AUC = AUC; results.Sensitivity = sensitivity; results.Specificity = specificity; results.Accuracy = accuracy;
        results.CI = CI;
        save(['testResultsRF_',[fSetNames{f},'Clinic'],'_',nameOutcomes{o}],'results'), clear results
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics