text:
function [opts,rem_opts] = vgg_argparse(opts,varargin)

%VGG_ARGPARSE  Parse variable arguments into a structure
%  opts = vgg_argparse(inopts,varargin)
%    inopts: structure (cells array) listing valid members and default values
%    varargin: variable arguments of form '<name>',<value>,...
%    opts: opts modified by varargin
%
%  Example:
%    function f = foo(varargin)
%    opts = vgg_argparse(struct('maxiters',10,'verbose',0), varargin)
%    ...
%
%  An unknown option (ie, present in varargin but absent in inopts)
%  causes an error. Calling the function as
%  [opts,rem_opts] = vgg_argparse(inopts,varargin) returns the unknown
%  option(s) in rem_opts for later use rather than causes an error.
%
%  May also use OPTS = VGG_ARGPARSE(OPTS, ASTRUCT) where ASTRUCT is a struct
%  of options.

% Author: Mark Everingham <me@robots.ox.ac.uk>
% modified by werner, Jan 03
% Date: 16 Jan 02

if iscell(opts)
  opts=struct(opts{:});
end

if length(varargin) & iscell(varargin{1})
    if isempty(varargin{1})
        inopts = struct([]);
    else
        inopts=struct(varargin{1}{:});
    end
else
    if isempty(varargin)
        inopts = struct([]);
    elseif isstruct(varargin{1})
        inopts = varargin{1};
    else
        inopts=struct(varargin{:});
    end
end

rem_opts = [];
fn = fieldnames(inopts);
for i=1:length(fn)
    if isfield(opts,fn{i})
        %opts.(fn{i})=inopts.(fn{i});
        opts = setfield(opts,fn{i},getfield(inopts,fn{i}));
    else
        if nargout < 2
            error(sprintf('bad argument: ''%s''', fn{i}));
        else
            %rem_opts.(fn{i}) = inopts.(fn{i});
            rem_opts = setfield(rem_opts,fn{i},getfield(inopts,fn{i}));
        end
    end
end

repo: SelfExSR


text:
function v = mean2( f )
%MEAN2   Mean of a SEPARABLEAPPROX.
%   V = MEAN2(F) returns the mean of a SEPARABLEAPPROX:
%
%   V = 1/A*integral2( f )
%
% 	where the A is the area of the domain of F.
%
% See also MEAN, STD2.

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

% Empty check:
if ( isempty( f ) )
    return
end

% Apply the formula:
v = sum2( f ) / domainarea( f );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = a123 ( )

%*****************************************************************************80
%
%% A123 returns the A123 matrix.
%
%  Example:
%
%    1 2 3
%    4 5 6
%    7 8 9
%
%  Properties:
%
%    A is integral.
%
%    A is not symmetric.
%
%    A is singular.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(3,3), the matrix.
%
  a = zeros ( 3, 3 );

  k = 0;
  for i = 1 : 3
    for j = 1 : 3
      k = k + 1;
      a(i,j) = k;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function T45=T45(h)
T45=2.78e-3*h+139.1;

author: Sable

repo: mcbench-benchmarks


text:
function elec_write_3dspace(p)

% elec_write_3dspace - Write electrode coordinate file (ascii).
%
% Useage: elec_write_3dspace(p)
%
% Write out the electrode labels, type, and Cartesian (x,y,z)
% coordinates in the format of Neuroscan 3Dspace ascii files.
%
% Each row of the file comprises an electrode label, an
% electrode type code (see below), and the x,y,z
% coordinates (cm). Each field is separated by spaces.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no express or implied warranties
% History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\nELEC_WRITE_3DSPACE...\n'); tic;

if ~exist('p', 'var'),
    msg = sprintf('...no input p struct.\n\n');
    error(msg);
end

[path,name,ext] = fileparts(strcat(p.elec.path,filesep,p.elec.file));
ext = '.dat';
file = fullfile(path,[name ext]);

fprintf('...writing electrode data to:\n\t%s\n', file);

write_elec(file,p.elec.data);

t=toc; fprintf('...done (%5.2f sec).\n\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function write_elec(file,elec);

fid = fopen(file,'w','ieee-le');

if(fid == -1),
    fprintf('...could not open file:\n\t%s',file);
    return;
else

    % Write out the electrode coordinates

    fprintf('...converting from meters to cm.\n');
    % Convert to Neuroscan 3Dspace coordinate
    % metric (cm) from meters.

    nasion = elec.nasion .* 100;
    rpa    = elec.rpa .* 100;
    lpa    = elec.lpa .* 100;

    x = elec.x .* 100;
    y = elec.y .* 100;
    z = elec.z .* 100;

    ref = elec.ref .* 100;

    centroid = elec.centroid .* 100;

    fprintf('...writing output coordinates.\n');

    fprintf(fid,'%+10s\t%3d\t%+12.6f\t%+12.6f\t%+12.6f\n','Nasion',  110, nasion);
    fprintf(fid,'%+10s\t%3d\t%+12.6f\t%+12.6f\t%+12.6f\n','Left',    108, lpa);
    fprintf(fid,'%+10s\t%3d\t%+12.6f\t%+12.6f\t%+12.6f\n','Right',   114, rpa);

    Nelec = size(x,1);
    type = 69;
    for e = 1:Nelec,
        fprintf(fid,'%+10s\t%3d\t%+12.6f\t%+12.6f\t%+12.6f\n',char(elec.label(e)), type, x(e), y(e), z(e));
    end

    fprintf(fid,'%+10s\t%3d\t%+12.6f\t%+12.6f\t%+12.6f\n','Centroid', 99, centroid(1),centroid(2),centroid(3));
    fprintf(fid,'%+10s\t%3d\t%+12.6f\t%+12.6f\t%+12.6f\n','Ref',     120, ref(1),ref(2),ref(3));

    fclose(fid);

end

return

author: PatternRecognition

repo: OpenBMI


text:
function c8vec_norm_l1_test ( )

%*****************************************************************************80
%
%% C8VEC_NORM_L1_TEST tests C8VEC_NORM_L1;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2015
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8VEC_NORM_L1_TEST\n' );
  fprintf ( 1, '  C8VEC_NORM_L1 computes the L1 norm of a C8VEC.\n' );

  a = c8vec_indicator ( n );

  c8vec_print ( n, a, '  The "indicator" vector:' );

  value = c8vec_norm_l1 ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L1 norm = %g\n', value );

  return
end


text:
function c=ref_dfft(f,order)
%REF_DFFT  Reference Discrete Fractional Fourier Transform
%   Usage:  c=ref_dfft(f,order);
%

L=size(f,1);

% Create matrix representation of the DFT
F=idft(eye(L));

c=(F^order)*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function time = compute_wait(plaza)
[a,b] = size(plaza);
time = 0;
for i = 1:a
for j = 1:b
time = time + (plaza(i,j) > 0);
end
end

author: Eurus-Holmes

path: github-repos/MATLAB/Eurus-Holmes-Mathematical_Modeling/Mathematical_Modeling-cb9dd53af84ffbd455ec62ab89886351e0ec98d9/Cellular automaton/元胞自动机/1月13日课件（元胞自动机）/程序/The Booth Tolls for Thee/compute_wait.m


text:
function [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains5(tpl)
% create_restrictions_and_markov_chains5 -- creates restrictions and
% markov chains for the SVAR model in which both coefficients and variance
% for the monetary policy equation are changing with two independent Markov
% processes
%
% ::
%
%
%   [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains5(tpl)
%
% Args:
%
%    - **tpl** [struct]: template created for SVAR objects
%
% Returns:
%    :
%
%    - **lin_restr** [cell]: two column-cell (see below). The first column
%    contains COEF objects or linear combinations of COEF objects, which are
%    themselves COEF objects.
%
%    - **nonlin_restr** [cell]: one column-cell containing inequality or
%    nonlinear restrictions
%
%    - **tpl** [struct]: modified template
%
% Note:
%
%    - The syntax to construct an advanced COEF object is
%    a=coef(eqtn,vname,lag,chain_name,state)
%      - **eqtn** [integer|char]: integer or variable name
%      - **vname** [integer|char]: integer or variable name
%      - **lag** [integer]: integer or variable name
%      - **chain_name** [char]: name of the markov chain
%      - **state** [integer]: state number
%
%    - RISE sorts the endogenous variables alphabetically and use this order
%    to tag each equation in SVAR and RFVAR models.
%
%    - The lag coefficients are labelled a0, a1, a2,...,ak, for a model with k
%    lags. Obviously, a0 denotes the contemporaneous coefficients.
%
%    - The constant terms labelled c_1_1, c_2_2,...,c_n_n, for a model with n
%    endogenous variables.
%
%    - The standard deviations labelled sig_1_1, sig_2_2,...,sig_n_n, for a
%    model with n endogenous variables.
%
% Example:
%    coef('pi','ygap',0,'policy',1)
%
%    coef(2,3,0,'policy',1)
%
%    coef(2,'ygap',0,'policy',1)
%
%    coef('pi',3,0,'policy',1)
%
%    See also:

% The parameter restrictions are identical to those in the model
% with regime switching in the policy coefficients only. Hence the mp_coef
% markov chain will also be common to those two models.
[lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains3(tpl);

% We add the volatility Markov chain from the model in which only the
% volatility of the monetary policy equation changes. N.B: In the process,
% we want to make sure we do not over-write the restrictions above!
%--------------------------------------------------------------------------
[~,~,tpl]=create_restrictions_and_markov_chains4(tpl);

end

author: jmaih


text:
function TS_ConvertToTable(hctsaFile)
% TS_ConvertToTable     Converts a file using the old structure array format -> table
%
%---INPUTS:
% hctsaFile, an hctsa file to convert

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

if nargin < 1
    hctsaFile = 'HCTSA.mat';
end

[~,TimeSeries,Operations,whatDataFile] = TS_LoadData(hctsaFile,false);
MasterOperations = TS_GetFromData(whatDataFile,'MasterOperations');

% Convert the three structure arrays to tables:
if isstruct(TimeSeries)
    TimeSeries = struct2table(TimeSeries);
end
if isstruct(Operations)
    Operations = struct2table(Operations);
end
if isstruct(MasterOperations)
    MasterOperations = struct2table(MasterOperations);
end

% Save back to file
fprintf(1,'Converted TimeSeries, Operations, and MasterOperations to table, saving back to %s...',whatDataFile);
save(whatDataFile,'TimeSeries','Operations','MasterOperations','-append')
fprintf(1,'Saved.\n');

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function z = uplus( x )

%   Disciplined convex/geometric programming information for UPLUS (+):
%      Unary plus may be used in DCPs and DGPs without restrictions.

z = x;

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% batch run full-clustering on all fish

% f.pushbutton_autoclus_Callback

global VAR;

data_masterdir = GetCurrentDataDir();

% range_fish = [5,6,7];
% M_ClusGroup = [2,2,2,2];
% M_Cluster = [1,1,1,1];
range_fish = 2:2;
M_ClusGroup = 1;
M_Cluster = 1;
M_stim = [1];
% M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2];

%%
M_param = 5:25:30%0.3:0.1:0.8;

ParamScores = cell(length(M_param),length(range_fish));
ParamScores_raw = cell(length(M_param),length(range_fish));

for k_param = 1:length(M_param),
    numK2 = M_param(k_param);
    %thres_split = M_param(k_param);
    setappdata(hfig,'thres_split',thres_split);

    for k_fish = 1:length(range_fish),
        i_fish = range_fish(k_fish);
        disp(i_fish);
        LoadFullFish(hfig,i_fish);
        absIX = getappdata(hfig,'absIX');

        %% partitions for CV
        timelists = getappdata(hfig,'timelists');
        timelists_names = getappdata(hfig,'timelists_names');
        periods = getappdata(hfig,'periods');
        if length(periods)>1,
            timelistsCV = cell(length(M_stim),2);
            for k_stim = 1:length(M_stim), % :3
                i_stim = M_stim(k_stim);
                TL = timelists{i_stim};
                period = periods(i_stim);
                nrep = size(TL,2)/periods(i_stim); % integer
                n = floor(nrep/2);
                timelistsCV{k_stim,1} = TL(1):TL(n*period);
                timelistsCV{k_stim,2} = TL(1+n*period):TL(2*n*period);
            end
        end

        %% CV loop: auto-clustering with the partitions
        Score = zeros(length(M_stim),2);

        for k_stim = 1:length(M_stim), % :3
            i_stim = M_stim(k_stim);
            NumClus = zeros(1,2);
            CIX = cell(1,2);
            GIX = cell(1,2);
            for k = 1:2,
                %% Cluster to start auto-clustering
                i_ClusGroup = M_ClusGroup(k_fish);
                i_Cluster = M_Cluster(k_fish);
                ClusGroup = VAR(i_fish).ClusGroup{i_ClusGroup};
                numK = ClusGroup(i_Cluster).numK;
                gIX = ClusGroup(i_Cluster).gIX;
                cIX_abs = ClusGroup(i_Cluster).cIX_abs; % convert absolute index to index used for this dataset
                [~,cIX] = ismember(cIX_abs,absIX);

                % ~UpdateTimeIndex
                tIX = timelistsCV{k_stim,k};
                M_0 = GetTimeIndexedData_Default_Direct(hfig,cIX,tIX,'isAllCells');

                isWkmeans = 1;
                [cIX,gIX] = AutoClustering(cIX,gIX,absIX,i_fish,M_0,isWkmeans,numK2);

                NumClus(k) = length(unique(gIX));
                CIX{k} = cIX;
                GIX{k} = gIX;
            end
            % plot cell-matching figure
            Score(k_stim,1) = HungarianCV(NumClus(1),NumClus(2),CIX{1},CIX{2},GIX{1},GIX{2},timelists_names{i_stim});
            Score(k_stim,2) = HungarianCV(NumClus(2),NumClus(1),CIX{2},CIX{1},GIX{2},GIX{1},timelists_names{i_stim});
        end

        ParamScores{k_param,k_fish} = mean(mean(Score));
        ParamScores_raw{k_param,k_fish} = Score;
    end

end

author: xiuyechen

repo: FishExplorer


text:
function f=loglikGaPExpWfix(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

[Hkt_r, cx, cy, Wxk] = loglikGaPreadparamWfix(x,varargin);
Hkt=exp(Hkt_r); %nonnegativity constrains

[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
P=Wxkbg*Hktbg; %current approximation

%Poisson contribution
t1=Vxt.*log(P) - P;
%Gamma contribution
% t2=(alpha-1)*log(Hkt)-1/beta*Hkt-alpha*log(beta)-log(gamma(alpha));
t2=-1/beta*Hkt;

f=sum(t1(:));%+sum(t2(:)));
f=-f; %conjugate gradient is mimimizing!
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% LFToolboxVersion - returns a string describing the current toolbox version

% Copyright (c) 2013-2020 Donald G. Dansereau

function VersionStr = LFToolboxVersion
VersionStr = 'v0.5.1';

repo: LFToolbox


text:
function[a]=morseafun(varargin)
%MORSEAFUN  Returns the generalized Morse wavelet amplitude or a-function.
%
%   MORSEAFUN is a low-level function called by many a number of the Morse
%   wavelet functions.
%
%   A=MORSEAFUN(GAMMA,BETA) returns the generalized Morse wavelet
%   amplitude, called "A_{BETA,GAMMA}" by Lilly and Olhede (2009).
%
%   By default, A is chosen such that the maximum of the frequency-
%   domain wavelet is equal to 2, the ``bandpass normalization.''
%
%   A=MORSEAFUN(GAMMA,BETA,'energy') instead returns the coefficient
%   giving the wavelet unit energy.
%
%   A=MORSEAFUN(K,GAMMA,BETA,'energy') returns the unit energy coefficient
%   appropriate for the Kth-order wavelet.  The default choice is K=1.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2006--2016 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1},'--t')
      morseafun_test;return
end

str='band';
if ischar(varargin{end})
    str=varargin{end};
    varargin=varargin(1:end-1);
end

k=1;
if length(varargin)==3
    k=varargin{1};
    varargin=varargin(2:end);
end

ga=varargin{1};
be=varargin{2};
if strcmpi(str(1:3),'ban')
    om=morsefreq(ga,be);
%    a=frac(2,(om.^be).*exp(-om.^ga));
    a=frac(2,exp(be.*log(om)-om.^ga));
    a(be==0)=2;
elseif strcmpi(str(1:3),'tes')
    a=sqrt(frac(2 *pi*ga.*2.^frac(2*be+1,ga),gamma(frac(2*be+1,ga))));
elseif strcmpi(str(1:3),'ene')
    r=frac(2*be+1,ga);
    a=double((2*pi*ga.*(2.^r).*exp(gammaln(k)-gammaln(k+r-1))).^(1/2));
end

function[]=morseafun_test

ga1=(2:1:9);
be1=(1:1:10);
[ga,be]=meshgrid(ga1,be1);
om=reshape(morsefreq(ga,be),10,8);

dom=0.01;
omgrid=permute((0:dom:20)',[3 2 1]);
omgrid=vrep(omgrid,length(ga1),2);
omgrid=vrep(omgrid,length(be1),1);

omgrid=omgrid.*vrep(om,size(omgrid,3),3);
a=morseafun(ga,be,'energy');

begrid=vrep(be,size(omgrid,3),3);
gagrid=vrep(ga,size(omgrid,3),3);
agrid=vrep(a,size(omgrid,3),3);

psi=agrid.*omgrid.^begrid.*exp(-omgrid.^gagrid);
psiint=vsum(psi.^2,3).*dom.*om./(2*pi);

reporttest('MORSEAFUN unit energy', allall(abs(psiint-1)<1e-2))

a2=morseafun(ga,be,'test');
reporttest('MORSEAFUN unit energy, alternate formulation', aresame(a,a2,1e-6));

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [energy, jumpPenalty, dataError] = energyL2Potts( u, f, gamma, A, isotropic )
%energyL2Potts Computes the energy of the L2 Potts functional

if exist('A', 'var') && not(isempty(A))
    Au = A * u;
    dataError = sum((Au(:) - f(:)).^2);
else
    dataError = sum((u(:) - f(:)).^2);
end

if isvector(u)
    % vectors
    jumpPenalty = sum(diff(u(:)) ~= 0);
else
    % matrices
    nJumpComp = 0; % jumps in compass directions
    nJumpDiag = 0; % jumps in diagonal directions

    % count jumps
    for i = 1:size(u, 1)
        for j = 1:size(u, 2)-1
            if ~all((u(i,j,:) == u(i,j+1,:)))
                nJumpComp = nJumpComp + 1;
            end
        end
    end
    for i = 1:size(u, 1)-1
        for j = 1:size(u, 2)
            if ~all((u(i,j,:) == u(i+1,j,:)))
                nJumpComp = nJumpComp + 1;
            end
        end
    end
    for i = 1:size(u, 1)-1
        for j = 1:size(u, 2)-1
            if ~all((u(i,j,:) == u(i+1,j+1,:)))
                nJumpDiag = nJumpDiag + 1;
            end
        end
    end
    for i = 1:size(u, 1)-1
        for j = 2:size(u, 2)
            if ~all((u(i,j,:) == u(i+1,j-1,:)))
                nJumpDiag = nJumpDiag + 1;
            end
        end
    end

    % set weights (isotropic by default)
    if ~exist('isotropic', 'var') || isotropic
        omega1 = sqrt(2) - 1;
        omega2 = 1 - sqrt(2)/2;
    else
        omega1 = 1;
        omega2 = 0;
    end

    % compute energy
    jumpPenalty = (omega1 * nJumpComp + omega2* nJumpDiag);

end

energy = gamma * jumpPenalty + dataError;

end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function [X, R, S] = snapshot_gen_sym(design, doas, wavelength, stype, t, ncov, sp)
%SNAPSHOT_GEN_SYM Generates snapshots from symbols.
%Syntax:
%   X = SNAPSHOT_GEN_SYM(design, doas, wavelength, stype[, t, ncov, sp]);
%   [X, R] = SNAPSHOT_GEN_SYM(design, doas, wavelength, stype[, t, ncov, sp]);
%   [X, R, S] = SNAPSHOT_GEN_SYM(design, doas, wavelength, stype[, t, ncov, sp]);
%Inputs:
%   design - Array design.
%   doas - DOA vector. For 2D DOAs, each column represents a DOA pair.
%   wavelength - Wavelength.
%   stype - Symbol type.
%   t - Number of snapshots.
%   ncov - Covariance matrix of the additive complex circular-symmetric
%          Gaussian noise. Can be a scalar, vector (for uncorrelated noise
%          with different powers), or a matrix.
%   sp - Sources powers.
%Outputs:
%   X - Snapshots, where each columns is a single snapshot.
%   R - Sample covariance matrix (averaged by the number of snapshots).
%   S - A source_count x snapshot_count matrix consists of source signal
%       vectors.
if nargin <= 6
    sp = 1;
end
if nargin <= 5
    ncov = 1;
end
if nargin <= 4
    t = 1;
end
A = steering_matrix(design, wavelength, doas);
[m, k] = size(A);
S_internal = gen_symbols(k, t, stype, sp);
X = A * S_internal + gen_ccsg(m, t, ncov);
if nargout >= 2
    R = (X*X')/t;
    if nargout == 3
        S = S_internal;
    end
end
end

function S = gen_symbols(m, n, stype, sp)
switch lower(stype)
    case 'bpsk'
        symbols = [-1 1];
    case 'qpsk'
        symbols = exp(1j*[1 3 5 7]/4*pi);
    otherwise
        error('Unsupported symbol type.')
end
S = symbols(randi([1 length(symbols)], m, n));
if isscalar(sp)
    S = sqrt(sp) * S;
else
    S = bsxfun(@times, sqrt(sp(:)), S);
end
end

function X = gen_ccsg(m, n, cov)
X0 = randn(m, n) + 1j * randn(m, n);
if isscalar(cov)
    X = sqrt(cov/2) * X0;
elseif isvector(cov)
    X = bsxfun(@times, X0, cov(:));
else
    C = sqrtm(cov/2);
    X = C*X0;
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
%%*****************************************************************
%% convertcmpsdp: convert SDP with complex data into one
%%                with real data by converting
%%
%%  C - sum_{k=1}^m yk*Ak psd
%%  to
%%  [CR,-CI] - sum ykR*[AkR,-AkI] psd
%%  [CI, CR]           [AkI, AkR]
%%
%%  ykI = 0 for k = 1:m
%%
%% [bblk,AAt,CC,bb] = convertcmpsdp(blk,A,C,b);
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [bblk,AAt,CC,bb,iscmp] = convertcmpsdp(blk,A,C,b);

  m = length(b);
  [pp,mm] = size(A);
  if (pp ~= size(blk,1))
     error('blk and A not compatible');
  end
  numblk = size(blk,1);
  iscmp = zeros(numblk,m+1);
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     len = size(A(p),2);
     for k = 1:len
        if ~isempty(A{p,k})
           iscmp(p,k) = 1-isreal(A{p,k});
        end
     end
     iscmp(p,m+1) = 1-isreal(C{p});
  end
  iscmp = norm(iscmp,'fro');
%%
  if (iscmp == 0)
     %% data is real
     bblk = blk; AAt = A; CC = C; bb = b;
     return;
  end
%%
  bb = [real(b)];
  bblk = cell(size(blk,1),2);
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     if (size(pblk{2},1) > size(pblk{2},2))
        pblk{2} = pblk{2}';
     end
     if strcmp(pblk{1},'s')
        ss = [0,cumsum(pblk{2})];
        ss2 = [0,cumsum(2*pblk{2})];
        n = sum(pblk{2});
        n2 = sum(pblk{2}.*(pblk{2}+1))/2;
        AR = cell(1,m); Ctmp = sparse(2*n,2*n);
        if (size(A{p},1)==n2 & size(A{p},2)==m);
           Atype = 1;
        elseif (size(A(p),1)==1 & size(A(p),2)==1);
           Atype = 2;
        else
           error('convertcmp: At is not properly coded');
        end
        for k = 0:m
           if (k == 0)
              Ak = C{p};
           else
              if (Atype == 1)
                 Ak = smat(pblk,A{p}(:,k),1);
              elseif (Atype == 2)
                 Ak = A{p,k};
              end
           end
           Atmp = sparse(2*n,2*n);
           if (length(pblk{2}) == 1)
              tmp = [real(Ak),-imag(Ak);  imag(Ak), real(Ak)];
              if (k==0)
                 Ctmp = tmp;
              else
                 Atmp = tmp;
              end
           else
              for j = 1:length(pblk{2})
                 idx = [ss(j)+1: ss(j+1)];
                 Akj = Ak(idx,idx);
                 tmp = [real(Akj),-imag(Akj);  imag(Akj), real(Akj)];
                 idx2 = [ss2(j)+1: ss2(j+1)];
                 if (k==0)
                    Ctmp(idx2,idx2) = tmp;
                 else
   	            Atmp(idx2,idx2) = tmp;
                 end
              end
           end
           if (k==0);
              CC{p,1} = Ctmp;
           else
              AR{k} = Atmp;
           end
        end
        bblk{p,1} = 's'; bblk{p,2} = 2*pblk{2};
        AAt(p,1) = svec(bblk(p,:),AR);
     elseif strcmp(pblk{1},'q');
        error('SOCP block with complex data is currently not allowed');
     elseif strcmp(pblk{1},'l');
        if isreal(A{p}) & isreal(C{p})
           bblk(p,:) = blk(p,:);
           AAt{p,1} = A{p}; CC{p,1} = C{p};
        else
           error('data for linear block must be real');
        end
     elseif strcmp(pblk{1},'u');
        if isreal(A{p}) & isreal(C{p})
           bblk(p,:) = blk(p,:);
           AAt{p,1} = A{p}; CC{p,1} = C{p};
        else
           error('data for unrestricted block must be real');
        end
     end
  end
%%*********************************************************


text:
function status = exportTable(table, pth);
% Write out a table to a text file.
%
% status = exportTable(table, [pth=dialog]);
%
% table should be a cell array of strings (1D or 2D).
% THis function was written for tables produces by the MATLAB
% anova functions.
%
% ras, 11/2007.
if notDefined('table'), error('Need an ANOVA table .'); end
if notDefined('pth'),
	pth = mrvSelectFile('w', 'txt', 'Specify ANOVA Table file');
end

status = -1;

fid = fopen(pth, 'w');
if fid < 1
	error( sprintf('Couldn''t open file %s for writing.', pth) )
end

for i = 1:size(table, 1)
	for j = 1:size(table, 2)
		val = table{i,j};
		if isnumeric(val)
			fprintf(fid, '%s\t', num2str(val));
		else
			fprintf(fid, '%s\t', val);
		end
	end

	fprintf(fid, '\n');
end

st = fclose(fid);

if st==0  % ok flag from fclose
	status = 1;
	fprintf('Saved table as %s.\n', pth);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ convertedstr ] = convertSBMLID( str, toSBML )
%CONVERTSBMLID converts the given str to a valid SBML ID
% USAGE:
%
%       convertedstr = convertSBMLID(str,toSBML)
%
% INPUT:
%
%       str:        The String to convert
%
% OPTIONAL INPUT:
%       toSBML:     Whether to convert to SBML format (or undo a
%                   conversion) (default true, ie convert to SBML)
%
% OUTPUT:
%
%       convertedstr:   The converted String.
%
% .. Authors:
%       - Thomas Pfau May 2017

if ~exist('toSBML','var')
    toSBML = true;
end
if toSBML
    convertedstr = regexprep(str,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
else
    convertedstr = regexprep(str,'__([0-9]+)__','${char(str2num($1))}');
    %Biomodels specific conversions for Input of Biomodels
    convertedstr = strrep(convertedstr,'-DASH-','-');
    convertedstr = strrep(convertedstr,'_DASH_','-');
    convertedstr = strrep(convertedstr,'_FSLASH_','/');
    convertedstr = strrep(convertedstr,'_BSLASH_','\');
    convertedstr = strrep(convertedstr,'_LPAREN_','(');
    convertedstr = strrep(convertedstr,'_LSQBKT_','[');
    convertedstr = strrep(convertedstr,'_RSQBKT_',']');
    convertedstr = strrep(convertedstr,'_RPAREN_',')');
    convertedstr = strrep(convertedstr,'_COMMA_',',');
    convertedstr = strrep(convertedstr,'_PERIOD_','.');
    convertedstr = strrep(convertedstr,'_APOS_','''');
end

end

author: opencobra

repo: cobratoolbox


text:
function net = cnn_initNet_pairwiseModel( networkFile, structuredLoss, numPairwiseClusters, extraLayers, meanImage )
%cnn_initNet_pairwiseModel initializes the pairwise model
%
% net = cnn_initNet_pairwiseModel( networkFile, structuredLoss, numPairwiseClusters, extraLayers )
%
% Input:
%   networkFile - file with the pretrained network (the dropout layers should already be inside)
%   structuredLoss - type of structure loss to use: 'svmStructCompact' or 'logisticScoresCompact'
%   numPairwiseClusters - number of clusters in the pairwise potentials
%   extraLayers - structure containing information about extra layers (default: [] - no extra layers). Fields:
%       unary - vector containing numbers of nodes in the extra layers for the unary potentials
%       pairwise - vector containing numbers of nodes in the extra layers for the pairwise potentials
%       useDropout - flag showing whether to use dropout in the extra layers
%       cutoffConvolutions - parameter to show if there is need to cut off more layers with weights from the pretrained net
%           One such layer is cut off anyway (connected to the classes from the pretraining task)
%   meanImage - if provided, sets the new mean image for the normalization
%
% Output:
%   net - structure containing the pairwise model

if ~exist('extraLayers', 'var') || isempty(extraLayers)
    extraLayers = struct;
    extraLayers.unary = [];
    extraLayers.pairwise = [];
    extraLayers.useDropout = true;
    extraLayers.cutoffConvolutions = 0;
end

% the initialization parameters
scal = 1 ;
init_bias = 0.1;
dropout_rate = 0.5;

%% read the network file
pretrainedNet = load( networkFile, '-mat' );
if ~isfield(pretrainedNet, 'layers')
    if ~isfield(pretrainedNet, 'net')
        error('cnn_initNet_pairwiseModel:incorrectInitialization', 'File to initialize the network is of incorrect format')
    else
        pretrainedNet = pretrainedNet.net;
    end
end
pretrainedNet.layers(end) = []; % cut off the loss layer
if ~exist('meanImage', 'var')
    meanImage = pretrainedNet.normalization.averageImage;
end

%% determine the number of output features
% find the last convolutional layer
lastConvolutionalLayer = length(pretrainedNet.layers);
cutoffConvolutions = 0;
while (~isequal( pretrainedNet.layers{ lastConvolutionalLayer }.type, 'conv' ) || cutoffConvolutions < extraLayers.cutoffConvolutions) ...
        && lastConvolutionalLayer > 1
    if isequal( pretrainedNet.layers{ lastConvolutionalLayer }.type, 'conv' )
        cutoffConvolutions = cutoffConvolutions + 1;
    end
    lastConvolutionalLayer = lastConvolutionalLayer - 1;
end
if isequal( pretrainedNet.layers{ lastConvolutionalLayer }.type, 'conv' )
    numFeatures = size( pretrainedNet.layers{ lastConvolutionalLayer }.weights{1}, 3 );
else
    error('cnn_initNet_pairwiseModel:wrongNetwork', 'Cannot determine the number of features from the pretrained network');
end

%% create the feature extractor network
net = struct;
net.featureExtractor = struct;
net.featureExtractor.layers = pretrainedNet.layers;
net.featureExtractor.layers(lastConvolutionalLayer : 1 : end) = []; % cut off the last layer

%% add the loss layer
net.lossLayer = struct;
net.lossLayer.type = structuredLoss;

numClasses = 2;
switch net.lossLayer.type
    case 'svmStructCompact'
        numUnaryOutputs = 1;
        numPairwiseOutputs = numPairwiseClusters;
    case 'logisticScoresCompact'
        numUnaryOutputs = 1;
        numPairwiseOutputs = numPairwiseClusters;
    otherwise
        error('cnn_initNet_pairwiseModel:unknownLoss',['The structured loss is not recognised: ', net.lossLayer.type]);
end

%% add the unary network
net.unaryNetwork = struct;
net.unaryNetwork.layers = cell(0, 0);

numExtraLayers = length(extraLayers.unary);
numNodesExtraLayersUnary = [ numFeatures; extraLayers.unary(:)];
for iLayer = 1 : numExtraLayers
    net.unaryNetwork.layers{end+1} = struct('type', 'conv', ...
        'weights', {{ 0.01/scal * randn(1, 1, numNodesExtraLayersUnary(iLayer), numNodesExtraLayersUnary(iLayer+1),'single'),... % filters
        init_bias*ones(1,numNodesExtraLayersUnary(iLayer+1),'single') }}, ... % biases
        'stride', 1, ...
        'pad', 0, ...
        'learningRate', [1, 2], ...
        'weightDecay', [1, 0], ...
        'name', ['unaryNetwork_conv', num2str(iLayer)]) ;
    net.unaryNetwork.layers{end+1} = struct('type', 'relu') ;
    if extraLayers.useDropout
        net.unaryNetwork.layers{end+1} = struct('type', 'dropout', 'rate', dropout_rate) ;
    end
end
net.unaryNetwork.layers{end+1} = struct('type', 'conv', ...
    'weights', {{ 0.01/scal * randn(1, 1, numNodesExtraLayersUnary(end), numUnaryOutputs,'single'),... % weights
    init_bias * ones(1, numUnaryOutputs, 'single') }}, ... % biases
    'stride', 1, ...
    'pad', 0, ...
    'learningRate', [1, 2], ...
    'weightDecay', [1, 0], ...
    'name', 'unaryNetwork_convLast') ;

%% add the pairwise network
net.pairwiseNetwork = struct;
net.pairwiseNetwork.layers = cell(0, 0);

numExtraLayers = length(extraLayers.pairwise);
numNodesExtraLayersPairwise = [ 2 * numFeatures; extraLayers.pairwise(:)];
for iLayer = 1 : numExtraLayers
    net.pairwiseNetwork.layers{end+1} = struct('type', 'conv', ...
        'weights', {{ 0.01/scal * randn(1, 1, numNodesExtraLayersPairwise(iLayer), numNodesExtraLayersPairwise(iLayer+1),'single'), ... % filters
        init_bias*ones(1,numNodesExtraLayersPairwise(iLayer+1),'single') }}, ... % biases
        'stride', 1, ...
        'pad', 0, ...
        'learningRate', [1, 2], ...
        'weightDecay', [1, 0], ...
        'name', ['pairwiseNetwork_conv', num2str(iLayer)]) ;
    net.pairwiseNetwork.layers{end+1} = struct('type', 'relu') ;
    if extraLayers.useDropout
        net.pairwiseNetwork.layers{end+1} = struct('type', 'dropout', 'rate', dropout_rate) ;
    end
end
net.pairwiseNetwork.layers{end+1} = struct('type', 'conv', ...
    'weights', {{ 0.01/scal * randn(1, 1, numNodesExtraLayersPairwise(end), numPairwiseOutputs,'single'), ... % filters
    init_bias * ones(1, numPairwiseOutputs, 'single') }}, ... % biases
    'stride', 1, ...
    'pad', 0, ...
    'learningRate', [1, 2], ...
    'weightDecay', [1, 0], ...
    'name', 'pairwiseNetwork_convLast') ;

%% collect all trainable layers
net.layers = cell(0,0);

for iLayer = 1 : length(net.featureExtractor.layers)
    if isequal( net.featureExtractor.layers{iLayer}.type, 'conv' )
        net.layers{end+1} = net.featureExtractor.layers{iLayer};
        net.featureExtractor.layers{iLayer} = struct( 'type', 'convPtr', 'index', length(net.layers) );
    end
end
for iLayer = 1 : length(net.unaryNetwork.layers)
    if isequal( net.unaryNetwork.layers{iLayer}.type, 'conv' )
        net.layers{end+1} = net.unaryNetwork.layers{iLayer};
        net.unaryNetwork.layers{iLayer} = struct( 'type', 'convPtr', 'index', length(net.layers) );
    end
end
for iLayer = 1 : length(net.pairwiseNetwork.layers)
    if isequal( net.pairwiseNetwork.layers{iLayer}.type, 'conv' )
        net.layers{end+1} = net.pairwiseNetwork.layers{iLayer};
        net.pairwiseNetwork.layers{iLayer} = struct( 'type', 'convPtr', 'index', length(net.layers) );
    end
end

%% mean image normalization
% set the new mean image
net.normalization.averageImage = meanImage;

%% class order
net.classes = struct;
net.classes.name = {'head', 'background'};
net.classes.description = {'head', 'background'};

end

author: aosokin


text:
classdef billiard

    % this are the dimensions of the table
    properties(Constant = true)
        XMAX = 100;
        XMIN = 0;
        YMAX = 50;
        YMIN = 0;
    end

    methods
        function whatever = billiard()
            % step 1. Draw the billiard table
            whatever.drawtable();

            % step 2. Draw the balls
          % ===========================
          % draw a single ball
          x = 10;   % Initial x-position
          y = 25;   % Initial y-position
          u = 1;    % Initial x-velocity (moves in x # pixels per iteration)
          v = 1;    % Initial y-velcotiy (moves in y # pixels per iteration)
          d = 5;    % Diameter of the ball
          m = 1;    % The mass of the ball

          ball = rectangle('position',[x,y,d,d],...
                    'FaceColor','b',...
                    'curvature',1.0 ...
                    );
          % ===========================
            % Update position
          % ===========================
           % the initial position
           x_new = x;
           y_new = y;

           % Update position

           dt = 1; % the time step is constant

           for iteration = 1:950;

               disp(x_new)

               if (abs(x_new-whatever.XMAX+d)<eps) || (abs(x_new-whatever.XMIN)<eps)
                   u = -u;
               end

               if (abs(y_new-whatever.YMAX+d)<eps) || (abs(y_new-whatever.YMIN)<eps)
                   v = -v;
               end

           x_new = x_new + u*dt;
           y_new = y_new + v*dt;

           ball.Position = [x_new,y_new,d,d];

           drawnow

           end

           % ===========================
        end

        function drawtable(whatever)
            % properties of the table
            w = whatever.XMAX-whatever.XMIN;
            h = whatever.YMAX-whatever.YMIN;
            x0 = whatever.XMIN;
            y0 = whatever.YMIN;
            % draw the table
            figure( 'name'   ,'TheCrazyBilliard2000',...
                    'menubar','none',...
                    'units'  ,'pixels'...
                    );
            rectangle('position',[x0,y0,w,h],...
                      'FaceColor','w',...
                      'curvature',0.1...
                      );
            axis('equal','off');
        end

    end
end

repo: Aero-matlab


text:
report_this_filefun(mfilename('fullpath'));

str = [];
[newmatfile, newpath] = uiputfile([ hodo '*.dat'], 'Save As'); %Syntax change Matlab Version 7, no window positioning on macs

s = [a.Longitude   a.Latitude  a.Date  a.Date.Month  a.Date.Day  a.Magnitude  a.Depth a.Date.Hour a.Date.Minute a(:,10) a(:,11) a(:,12) a(:,13) a(:,14)   ];
fid = fopen([newpath newmatfile],'w') ;
fprintf(fid,'%7.3f  %7.3f %6.2f %6.0f  %6.0f %6.1f %6.2f  %6.0f  %6.0f  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f\n',s');
fclose(fid);
clear s
return

author: CelsoReyes


text:
% Crude approximations for the W function (used by bisect)
%
% Syntax:
% crude(xx,nb)
%
% xx is the argument, W(xx)
%
% nb is the branch of the W function needed:
% nb = 0 - upper branch, Wp
% nb <> 0 - lower branch, Wm
%
function c=crude(xx,nb)
%
% D. A. Barry, 23 June 2003 (d.a.barry@ed.ac.uk)
%
% Various constants
%
em=-exp(-1);
em9=-exp(-9);
c13=1/3;
em2=2/em;
s2=sqrt(2);
s21=2*s2-3;
s22=4-3*s2;
s23=s2-2;
if nb == 0
%
%  Calculations for crude Wp
%
   if xx <= 20
      reta=s2*sqrt(1-xx/em);
      an2=4.612634277343749*sqrt(sqrt(reta+1.09556884765625));
      c=reta/(1+reta/(3+(s21*an2+s22)*reta/(s23*(an2+reta))))-1;
   else
      zl=log(xx);
      c=log(xx/log(xx/zl^exp(-1.124491989777808/(.4225028202459761+zl))));
   end
else
%
%  Calculations for crude Wm
%
   if xx <= em9
      zl=log(-xx);
      t=-1-zl;
      ts=sqrt(t);
      c=zl-(2*ts)/(s2+(c13-t/(270+ts*127.0471381349219))*ts);
   else
      zl=log(-xx);
      eta=2-em2*xx;
      c=log(xx/log(-xx/((1-.5043921323068457*(zl+1))*(sqrt(eta)+eta/3)+1)));
   end
end


text:
%% Test Problems for CLP
clc
clear

%% CLP Options Function
clc
clpset

a = clpset

%% LP1 [-31.4]
clc
f = -[6 5]';
A = sparse([1,4; 6,4; 2, -5]);
rl = -Inf(3,1);
ru = [16;28;6];
lb = [0;0];
ub = [10;10];

opts = [];
opts.display = 1;
opts.maxiter = 15;
opts.algorithm = 'automatic';
opts.numThreads = 4;
[x,ff,e,i,lambda] = clp([],f,A,rl,ru,lb,ub,opts);

%% LP1 with redundant row [-31.4]
clc
f = -[6 5]';
A = sparse([1,4; 2 8; 6,4; 2, -5]);
rl = -Inf(4,1);
ru = [16;32;28;6];
lb = [0;0];
ub = [10;10];

opts = [];
opts.display = 1;
opts.maxiter = 15;
opts.algorithm = 3;
% opts.writeprob = 'prob.dat-s';
% opts.writesol = 'sol.dat-s';
[x,ff,e,i,lambda] = clp([],f,A,rl,ru,lb,ub,opts)

%% LP2 [2]
clc
f = [8,1]';
A = sparse([-1,-2;1,-4;3,-1;1,5;-1,1;-1,0;0,-1]);
rl = -Inf(7,1);
ru = [-4,2,21,39,3,0,0]';

[x,ff,e,i,lambda] = clp([],f,A,rl,ru)

%% LP3 [-3.75]
clc
f = -[-1, 2]';
A = sparse([2, 1;-4, 4]);
rl = -Inf(2,1);
ru = [5, 5]';

[x,ff,e,i,lambda] = clp([],f,A,rl,ru)

%% LP4 [-97.5]
clc
f = -[1 2 3]';
A = sparse([-1,1,1; 1,-3,1]);
b = [20,30]';
Aeq = sparse([1 1 1]);
A = [A;Aeq;-Aeq];
beq = 40;
ru = [b;beq;-beq];
rl = -Inf(size(ru));
lb = [0 0 0]';
ub =[40 inf inf]';

[x,ff,e,i,lambda] = clp([],f,A,rl,ru,lb,ub,opts)

%% LP5 [-30, no linear constraints]
clc
f = -[1, 2]';
lb = [0;0];
ub = [10;10];

[x,ff,e,i,lambda] = clp([],f,[],[],[],lb,ub,opts)

%% LP6 infeasible
clc
%Objective & Constraints
f = -[6 5]';
A = sparse([1,4; 6,4; 2, -5]);
rl = -Inf(3,1);
ru = [16;28;-30];
lb=[0;0];
ub=[10;10];

opts = [];
opts.display = 1;
opts.maxiter = 150;
opts.algorithm =5;

[x,ff,e,i,lambda] = clp([],f,A,rl,ru,lb,ub,opts)

%% LP7 unbounded
clc
f = -[1, 2]';
lb = [0;0];
ub = [10;Inf];

opts = [];
opts.display = 1;
opts.maxiter = 150;
opts.algorithm = 0;

[x,ff,e,i,lambda] = clp([],f,[],[],[],lb,ub,opts)

%% LARGE LP
clc
clear clp
prob = coinRead('maros-r7.mps');

opts = [];
opts.display = 1;
opts.maxiter = 10000;
opts.doPresolve = 1;
opts.algorithm = 5;
opts.numThreads = 4;

[~,ff,e,i,lambda] = clp(prob.H,prob.f,prob.A,prob.rl,prob.ru,prob.lb,prob.ub,opts)

%% AMPL Problem LP
clc
clear clp
prob = coinRead('prod.mps');

opts = [];
opts.display = 1;
opts.maxiter = 10000;
opts.algorithm = 0;
opts.doPresolve = 1;
opts.numThreads = 1;

[~,ff,e,i,lambda] = clp(prob.H,prob.f,prob.A,prob.rl,prob.ru,prob.lb,prob.ub,opts)

%% LP No Presolve Dual
clc
f = -[-1, 2]';
A = sparse([2, 1;-4, 4]);
rl = -Inf(2,1);
ru = [5, 5]';

opts = [];
opts.display = 1;
opts.maxiter = 15;
opts.algorithm = 0;
opts.doPresolve = 0;
opts.numThreads = 1;
[x,ff,e,i,lambda] = clp([],f,A,rl,ru,[],[],opts);

%% LP Empty after Presolve Abc
clc
f = -[-1, 2]';
A = sparse([2, 1;-4, 4]);
rl = -Inf(2,1);
ru = [5, 5]';

opts = [];
opts.display = 1;
opts.algorithm = 0;
opts.numThreads = 2;
[x,ff,e,i,lambda] = clp([],f,A,rl,ru,[],[],opts);

%% QP1 -2.83333333301227;
clc
H = speye(3);
f = -[2 3 1]';
A = sparse([1 1 1;3 -2 -3; 1 -3 2]);
b = [1;1;1];

opts = [];
opts.display = 2;
opts.algorithm = 5;
opts.numThreads = 1;

% prob=optiprob('qp',H,f,'ineq',A,b);
% coinWrite(prob,'testqp.mps')

[x,fval,e,i,l] = clp(H,f,A,-Inf(size(b)),b,[],[],opts)

%% QP2 -8.22222220552525
clc
H = (sparse([1 -1; -1 2]));
f = -[2 6]';
A = sparse([1 1; -1 2; 2 1]);
b = [2; 2; 3];
lb = [0;0];

opts = [];
opts.display = 1;
opts.algorithm = 3;
opts.objbias = 10;

[x,fval,e,i,l] = clp(tril(H),f,A,-Inf(size(b)),b,[],[],opts)

%% QP3 -6.41379310344827
H = tril(sparse([1 -1; -1 2]));
f = -[2 6]';
A = sparse([1 1; -1 2; 2 1]);
b = [2; 2; 3];
Aeq = sparse([1 1.5]);
beq = 2;
A = [A;Aeq]; rl = [zeros(size(b));beq]; ru = [b;beq];
lb = [0;0];
ub = [10;10];

opts = [];
opts.display = 2;
[x,fval,e,i,l] = clp(H,f,A,rl,ru,lb,ub,opts)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


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

%%
clear; close all; clc;

%% Syntax
% |p=csapsPar(varargin);|

%% Description
% This function aims to compute the cubic smoothing spline parameter p for
% the MATLAB csaps function. The input features the curve vertices V and
% the parameter pw. The smoothing parameter p is derived using:
%
% p=1./(1+(((h.^3)/6)*f));
%
% Where h is the point spacing (derived from V) and f is defined as:
%
% f=(1/pw)-1;
%
% See also MATLAB's csaps documentation:
%
% The interesting range for p is close to 1./(1+((h.^3)/6)). The following
% form is used introducing the factor f: p=1./(1+(((h.^3)/6)*f)). By using
% f=10 we obtain p=1./(1+((h.^3)/60)) which should result in a close
% following of the data. If instead f=0.1 is used, leading to
% p=1./(1+((h.^3)/0.6)), a smoother result is obtained.

%% Examples
%

%%
% Create saw-tooth example curve
V=[0 0 0; 1 2 0; 2 0 0; 3 -2 0; 4 0 0];

n1=2;
n2=2*n1;
[V1]=subCurve(V,n1);
[V2]=subCurve(V,n2);

%%
% Computing p parameters for the cubic smoothing spline based resampling

pw=0.5;

p1=csapsPar(V1,pw)
p2=csapsPar(V2,pw)

% Resample curves using smoothing parameters
n=100;
V1f=evenlySampleCurve(V1,n,p1,0);
V2f=evenlySampleCurve(V2,n,p2,0);

%%
% Visualizing curves

cFigure;
subplot(1,2,1); hold on;
title(['Input curve with ',num2str(size(V1,1)),' points'])
plotV(V1,'k.-','LineWidth',1,'MarkerSize',25);
plotV(V1f,'b-','LineWidth',3);
axis tight; axis equal; grid on; box on;

subplot(1,2,2); hold on;
title(['Input curve with ',num2str(size(V2,1)),' points'])
plotV(V2,'k.-','LineWidth',1,'MarkerSize',25);
plotV(V2f,'r-','LineWidth',3);
axis tight; axis equal; grid on; box on;

drawnow;

%%
%
% <<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:
function SO3VF = plus(SO3VF1, SO3VF2)
% overloads |SO3VF1 + SO3VF2|
%
% Syntax
%   SO3VF = SO3VF1 + SO3VF2
%   SO3VF = a + SO3VF1
%   SO3VF = SO3VF1 + a
%   SO3VF = SO3VF1 + SO3F;
%   SO3VF = SO3F + SO3VF2;
%
% Input
%  SO3VF1, SO3VF2 - @SO3VectorField
%  a - double, @vector3d
%  SO3F - @SO3Fun
%
% Output
%  SO3VF - @SO3VectorField
%

if isa(SO3VF2,'vector3d')
  SO3VF2 = SO3VF2.xyz.';
end

if isnumeric(SO3VF1) || isa(SO3VF1,'SO3Fun')
  SO3VF = SO3VF2;
  SO3VF.SO3F = SO3VF1 + SO3VF2.SO3F;
  return
end
if isnumeric(SO3VF2) || isa(SO3VF2,'SO3Fun')
  SO3VF = SO3VF2 + SO3VF1;
  return
end

ensureCompatibleSymmetries(SO3VF1,SO3VF2)
SO3VF1 = SO3VectorFieldHarmonic(SO3VF1);
SO3VF2 = SO3VectorFieldHarmonic(SO3VF2);
SO3VF = SO3VF1;
SO3VF.SO3F = SO3VF1.SO3F+SO3VF2.SO3F;

end

author: mtex-toolbox

repo: mtex

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


text:
function [Image,x0,y0] = move_circle2(Image,I_ellipse,xlo,xhi,ylo,yhi,ix,iy);

image_size = length(Image);

while max(max(Image(ix,iy))) > 1

    Image(ix,iy) = Image(ix,iy) - I_ellipse;

    x0 = ceil(rand*image_size);
    y0 = ceil(rand*image_size);

    nlo = x0 - xlo; nhi = x0 + xhi;
    ix = mod(nlo:nhi,length(Image));ix(ix==0)=length(Image);
    nlo = y0 - ylo; nhi = y0 + yhi;
    iy = mod(nlo:nhi,length(Image));iy(iy==0)=length(Image);

    Image(ix,iy) = Image(ix,iy) + I_ellipse;

end

author: Sable

repo: mcbench-benchmarks


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

% modified by Yu Xiang

function tracker = LK_tracking(frame_id, dres_image, dres_det, tracker)

% current frame + motion
J = dres_image.Igray{frame_id};
ctrack = apply_motion_prediction(frame_id, tracker);
w = tracker.dres.w(end);
h = tracker.dres.h(end);
BB3 = [ctrack(1)-w/2; ctrack(2)-h/2; ctrack(1)+w/2; ctrack(2)+h/2];
[J_crop, BB3_crop, bb_crop, s] = LK_crop_image_box(J, BB3, tracker);

num_det = numel(dres_det.x);
for i = 1:tracker.num
    BB1 = [tracker.x1(i); tracker.y1(i); tracker.x2(i); tracker.y2(i)];
    I_crop = tracker.Is{i};
    BB1_crop = tracker.BBs{i};

    % LK tracking
    [BB2, xFJ, flag, medFB, medNCC, medFB_left, medFB_right, medFB_up, medFB_down] = LK(I_crop, J_crop, ...
        BB1_crop, BB3_crop, tracker.margin_box, tracker.level_track);

    BB2 = bb_shift_absolute(BB2, [bb_crop(1) bb_crop(2)]);
    BB2 = [BB2(1)/s(1); BB2(2)/s(2); BB2(3)/s(1); BB2(4)/s(2)];

    ratio = (BB2(4)-BB2(2)) / (BB1(4)-BB1(2));
    ratio = min(ratio, 1/ratio);

    if isnan(medFB) || isnan(medFB_left) || isnan(medFB_right) || isnan(medFB_up) || isnan(medFB_down) ...
            || isnan(medNCC) || ~bb_isdef(BB2) || ratio < tracker.max_ratio
        medFB = inf;
        medFB_left = inf;
        medFB_right = inf;
        medFB_up = inf;
        medFB_down = inf;
        medNCC = 0;
        o = 0;
        score = 0;
        ind = 1;
        angle = -1;
        flag = 2;
        BB2 = [NaN; NaN; NaN; NaN];
    else
        % compute overlap
        dres.x = BB2(1);
        dres.y = BB2(2);
        dres.w = BB2(3) - BB2(1);
        dres.h = BB2(4) - BB2(2);
        if isempty(dres_det.fr) == 0
            overlap = calc_overlap(dres, 1, dres_det, 1:num_det);
            [o, ind] = max(overlap);
            score = dres_det.r(ind);
        else
            o = 0;
            score = -1;
            ind = 0;
        end

        % compute angle
        centerI = [(BB1(1)+BB1(3))/2 (BB1(2)+BB1(4))/2];
        centerJ = [(BB2(1)+BB2(3))/2 (BB2(2)+BB2(4))/2];
        v = compute_velocity(tracker);
        v_new = [centerJ(1)-centerI(1), centerJ(2)-centerI(2)] / double(frame_id - tracker.frame_ids(i));
        if norm(v) > tracker.min_vnorm && norm(v_new) > tracker.min_vnorm
            angle = dot(v, v_new) / (norm(v) * norm(v_new));
        else
            angle = 1;
        end
    end

    tracker.bbs{i} = BB2;
    tracker.points{i} = xFJ;
    tracker.flags(i) = flag;
    tracker.medFBs(i) = medFB;
    tracker.medFBs_left(i) = medFB_left;
    tracker.medFBs_right(i) = medFB_right;
    tracker.medFBs_up(i) = medFB_up;
    tracker.medFBs_down(i) = medFB_down;
    tracker.medNCCs(i) = medNCC;
    tracker.overlaps(i) = o;
    tracker.scores(i) = score;
    tracker.indexes(i) = ind;
    tracker.angles(i) = angle;
    tracker.ratios(i) = ratio;
end

% combine tracking and detection results
% [~, ind] = min(tracker.medFBs);
ind = tracker.anchor;
if tracker.overlaps(ind) > tracker.overlap_box
    index = tracker.indexes(ind);
    bb_det = [dres_det.x(index); dres_det.y(index); ...
        dres_det.x(index)+dres_det.w(index); dres_det.y(index)+dres_det.h(index)];
    tracker.bb = mean([repmat(tracker.bbs{ind}, 1, tracker.weight_tracking) bb_det], 2);
else
    tracker.bb = tracker.bbs{ind};
end

% compute pattern similarity
if bb_isdef(tracker.bb)
    pattern = generate_pattern(dres_image.Igray{frame_id}, tracker.bb, tracker.patchsize);
    nccs = distance(pattern, tracker.patterns, 1); % measure NCC to positive examples
    tracker.nccs = nccs';
else
    tracker.nccs = zeros(tracker.num, 1);
end

if tracker.is_show
    fprintf('\ntarget %d: frame ids ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%d ', tracker.frame_ids(i))
    end
    fprintf('\n');
    fprintf('target %d: medFB ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.medFBs(i))
    end
    fprintf('\n');

    fprintf('target %d: medFB left ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.medFBs_left(i))
    end
    fprintf('\n');

    fprintf('target %d: medFB right ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.medFBs_right(i))
    end
    fprintf('\n');

    fprintf('target %d: medFB up ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.medFBs_up(i))
    end
    fprintf('\n');

    fprintf('target %d: medFB down ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.medFBs_down(i))
    end
    fprintf('\n');

    fprintf('target %d: medNCC ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.medNCCs(i))
    end
    fprintf('\n');

    fprintf('target %d: overlap ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.overlaps(i))
    end
    fprintf('\n');
    fprintf('target %d: detection score ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.scores(i))
    end
    fprintf('\n');
    fprintf('target %d: flag ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%d ', tracker.flags(i))
    end
    fprintf('\n');
    fprintf('target %d: angle ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.angles(i))
    end
    fprintf('\n');
    fprintf('target %d: ncc ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.nccs(i))
    end
    fprintf('\n\n');
    fprintf('target %d: bb overlaps ', tracker.target_id);
    for i = 1:tracker.num
        fprintf('%.2f ', tracker.bb_overlaps(i))
    end
    fprintf('\n\n');

    if tracker.flags(ind) == 2
        fprintf('target %d: bounding box out of image\n', tracker.target_id);
    elseif tracker.flags(ind) == 3
        fprintf('target %d: too unstable predictions\n', tracker.target_id);
    end
end

author: yuxng


text:
function ldQ = VBA_logDet(Q,indIn,t)
% computes the log-determinant of matrix Q
% function ldQ = VBA_logDet(Q,indIn,t)
% IN:
%   - Q: nxn matrix
%   - indIn: vector of indices that defines a submatrix of Q, whose
%   determinant is computed {1:n}
%   - t: threshold on the eigenvalues of the submatrix {eps}
% OUT:
%   - ldQ: log-determinant of the submatrix of Q

if nargin < 3
    t = eps;
end
if nargin < 2 || isempty(indIn)
    dq = diag(Q);
    indIn = find(~isinf(dq)&dq~=0);
end
subQ = full(Q(indIn,indIn));
if isequal(subQ,eye(length(indIn))) % identity matrix
    ldQ = 0;
elseif isequal(subQ,diag(diag(subQ))) % diagonal matrix
    dQ = diag(subQ);
    ldQ = sum(log(dQ(abs(dQ)>t)));
else % full matrix
    if ~ VBA_isWeird (subQ)
        ev = eig(subQ);
    else
        ev = 0;
    end
    ev = ev(abs(ev)>t);
    ldQ = sum(log(ev));
end
ldQ = real(ldQ);

author: MBB-team

repo: VBA-toolbox


text:
function SetAppData(varargin)
for k=1:nargin; setappdata(0, inputname(k), varargin{k}); end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% feat is a cell array of matrix, each with TxD size. T is the number of
% frames and D is the number of dimensions.

function writeKaldiArchiveText(file_name, uttID, feat)

FID = fopen(file_name, 'W');

for i=1:length(uttID)
    PrintProgress(i, length(uttID), 500, 'Writing Kaldi archive');
    fprintf(FID, '%s [ ', uttID{i});
    write_matrix(FID, feat{i}, 1);
    fprintf(FID, ' ]\n');
end

fclose(FID);

author: singaxiong

repo: SignalGraph


text:
classdef  VideoMaker_TopOpt_density3DXmirror < VideoMaker_TopOpt
  properties

  end

  methods (Access = public)

        function Make_video_design_variable(obj,output_video_name)
            post = Postprocess_TopOpt_density();
            field2print = post.density_name;
            componentfield = post.density_name_component;
            obj.Make_video_density_field(field2print,componentfield,output_video_name)
        end

        function Make_video_design_variable_reg(obj,output_video_name)
            post = Postprocess_TopOpt_density();
            field2print = post.density_name_reg;
            componentfield = post.density_name_component_reg;
            obj.Make_video_density_field(field2print,componentfield,output_video_name)
        end

  end

  methods (Access = private)
        function Make_video_density_field(obj,field2print,componentfield,output_video_name)
            file_tcl_name = 'tcl_gid.tcl';
            file_list = obj.create_file_list(obj.iterations_to_print,obj.file_name,obj.files_folder);
            [output_video_name] = obj.replace_special_character(output_video_name);
            [output_photo]=obj.replace_special_character(strcat(obj.files_folder,'\',obj.file_name,'.png'));
            [file_list] = obj.replace_special_character(file_list);

            file_tcl_name_with_path = fullfile(obj.files_folder,file_tcl_name);
            file_path_in = fullfile(pwd,'FEM','PostProcess','Make_Video_density3DXmirror.tcl');
            filepath = obj.replace_special_character(file_path_in);

            fid = fopen(file_tcl_name_with_path,'w+');

            fprintf(fid,'GiD_Process PostProcess \n');
            fprintf(fid,'%s\n',['set arg1 "',file_list,'"']);
            fprintf(fid,'%s\n',['set arg2 "',output_video_name,'"']);
            fprintf(fid,'%s\n',['set arg3 "',field2print,'"']);
            fprintf(fid,'%s\n',['set arg4 "',componentfield,'"']);
            fprintf(fid,'%s\n',['set arg5 "',output_photo,'"']);

            fprintf(fid,'%s\n',['source "',filepath,'"']);
            fprintf(fid,'%s\n',['Make_Video_density $arg1 $arg2 $arg3 $arg4 $arg5']);
            fprintf(fid,'%s\n',['GiD_Process Mescape Quit']);
            fclose(fid);
            obj.execute_tcl_files(obj.gidPath,file_tcl_name_with_path)
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [fx,tx,pv,fv]=fxpefac(s,fs,tinc,m,pp)
%FXPEFAC PEFAC pitch tracker [FX,TT,PV,FV]=(S,FS,TINC,M,PP)
%
% Input:   s(ns)      Speech signal
%          fs         Sample frequency (Hz)
%          tinc       Time increment between frames (s) [0.01]
%                     or [start increment end]
%          m          mode
%                     'g' plot graph showing waveform and pitch
%                     'G' plot spectrogram with superimposed pitch using
%                         options pp.sopt [default: 'ilcwpf']
%                     'x' use external files for algorithm parameter
%                         initialization: fxpefac_g and fxpefac_w
%          pp         structure containing algorithm parameters
%
% Outputs: fx(nframe)     Estimated pitch (Hz)
%          tx(nframe)     Time at the centre of each frame (seconds).
%          pv(nframe)     Probability of the frame of being voiced
%          fv             structure containing feature vectors
%                           fv.vuvfea(nframe,2) = voiced/unvoiced GMM features

% References
%  [1]  S.Gonzalez and M. Brookes,
%       A pitch estimation filter robust to high levels of noise (PEFAC), Proc EUSIPCO,Aug 2011.

% Bugs/Suggestions
% (1) do long files in chunks
% (2) option of n-best DP

%	   Copyright (C) Sira Gonzalez and Mike Brookes 2011
%      Version: $Id: fxpefac.m 3601 2013-10-11 15:27:30Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent w_u m_u v_u w_v m_v v_v dpwtdef
% initialize persistent variables
if ~numel(w_u)

    % voiced/unvoiced decision based on 2-element feature vector
    % (a) mean power of the frame's log-freq spectrum (normalized so its short-term average is LTASS)
    % (b) sum of the power in the first three peaks
    %===== VUV
    if nargin>3 && any(m=='x')
        fxpefac_g;     % read in GMM parameters
        fxpefac_w;     % read in Weights parameters
    else
        w_u=[0.1461799 0.3269458 0.2632178 0.02331986 0.06360947 0.1767271 ]';

        m_u=[13.38533 0.4199435 ;
             12.23505 0.1496836 ;
             12.76646 0.2581733 ;
             13.69822 0.6893078 ;
             9.804372 0.02786567 ;
             11.03848 0.07711229 ];

        v_u=reshape([0.4575519 0.002619074 0.002619074 0.01262138 ;
             0.7547719 0.008568089 0.008568089 0.001933864 ;
             0.5770533 0.003561592 0.003561592 0.00527957 ;
             0.3576287 0.01388739 0.01388739 0.04742106 ;
             0.9049906 0.01033191 0.01033191 0.0001887114 ;
             0.637969 0.009936445 0.009936445 0.0007082946 ]',[2 2 6]);

        w_v=[0.1391365 0.221577 0.2214025 0.1375109 0.1995124 0.08086066 ]';

        m_v=[15.36667 0.8961554 ;
             13.52718 0.4809653 ;
             13.95531 0.8901121 ;
             14.56318 0.6767258 ;
             14.59449 1.190709 ;
             13.11096 0.2861982 ];

        v_v=reshape([0.196497 -0.002605404 -0.002605404 0.05495016 ;
             0.6054919 0.007776652 0.007776652 0.01899244 ;
             0.5944617 0.0485788 0.0485788 0.03511229 ;
             0.3871268 0.0292966 0.0292966 0.02046839 ;
             0.3377683 0.02839657 0.02839657 0.04756354 ;
             1.00439 0.03595795 0.03595795 0.006737475 ]',[2 2 6]);
    end
    %===== PDP
    %     dfm = -0.4238; % df mean
    %     dfv = 3.8968; % df variance (although treated as std dev here)
    %     delta = 0.15;
    %     dflpso=[dfm 0.5/(log(10)*dfv^2) -log(2*delta/(dfv*sqrt(2*pi)))/log(10)]; % scale factor & offset for df pdf
    %     dpwtdef=[1.0000, 0.8250, 1.3064, 1.9863]; % default DP weights
    dpwtdef=[1.0000, 0.8250, 0.01868, 0.006773, 98.9, -0.4238]; % default DP weights
    %===== END

end
% Algorithm parameter defaults

p.fstep=5;              % frequency resolution of initial spectrogram (Hz)
p.fmax=4000;            % maximum frequency of initial spectrogram (Hz)
p.fres = 20;            % bandwidth of initial spectrogram (Hz)
p.fbanklo = 10;         % low frequency limit of log filterbank (Hz)
p.mpsmooth = 21;       % width of smoothing filter for mean power
% p.maxtranf = 1000;      % maximum value of tranf cost term
p.shortut = 7;          % max utterance length to average power of entire utterance
p.pefact = 1.8;         % shape factor in PEFAC filter
p.numopt = 3;           % number of possible frequencies per frame
p.flim = [60 400];      % range of feasible fundamental frequencies (Hz)
p.w = dpwtdef;          % DP weights
% p.rampk = 1.1;          % constant for relative-amplitude cost term
% p.rampcz = 100;         % relative amplitude cost for missing peak
p.tmf = 2;              % median frequency smoothing interval (s)
p.tinc = 0.01;          % default frame increment (s)
p.sopt = 'ilcwpf';      % spectrogram options

% update parameters from pp argument

if nargin>=5 && isstruct(pp)
    fnq=fieldnames(pp);
    for i=1:length(fnq)
        if isfield(p,fnq{i})
            p.(fnq{i})=pp.(fnq{i});
        end
    end
end

% Sort out input arguments
if nargin>=3  && numel(tinc)>0
    p.tinc = tinc;   % 0.01 s between consecutive time frames
end
if nargin<4
    m='';
end

% Spectrogram of the mixture
fmin = 0; fstep = p.fstep; fmax = p.fmax;
fres = p.fres;  % Frequency resolution (Hz)
[tx,f,MIX]=spgrambw(s,fs,fres,[fmin fstep fmax],[],p.tinc);
nframes=length(tx);
txinc=tx(2)-tx(1);  % actual frame increment
%  ==== we could combine spgrambw and filtbankm into a single call to spgrambw or use fft directly ====
% Log-frequency scale
[trans,cf]=filtbankm(length(f),2*length(f)-1,2*f(end),p.fbanklo,f(end),'usl');
O = MIX*trans'; % Original spectrum in Log-frequency scale

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Amplitude Compression

% Calculate alpha based on LTASS ratios
ltass = stdspectrum(6,'p',cf);
auxf = [cf(1),(cf(1:end-1)+cf(2:end))./2,cf(end)];
ltass = ltass.*diff(auxf);                  % weight by bin width

% estimated ltass
O = O.*repmat(diff(auxf),nframes,1);     % weight spectrum by bin width
O1 = O;

if tx(end)<p.shortut                        % if it is a short utterance
    eltass = mean(O,1);                     % mean power per each frequency band
    eltass = smooth(eltass,p.mpsmooth);     % smooth in log frequency
    eltass= eltass(:).';                    % force a row vector

    % Linear AC
    alpha = (ltass)./(eltass);
    alpha = alpha(:).';
    alpha = repmat(alpha,nframes,1);
    O = O.*alpha;                           % force O to have an average LTASS spectrum

    % ==== should perhaps exclude the silent portions ***
else                                        % long utterance

    tsmo = 3; % time smoothing over 3 sec
    stt = round(tsmo/txinc);
    eltass = timesm(O,stt);
    eltass = smooth(eltass,p.mpsmooth);     % filter in time and log frequency

    % Linear AC
    alpha = repmat(ltass,nframes,1)./(eltass);
    O = O.*alpha;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the filter to detect the harmonics
ini = find(cf>3*cf(1));
sca = cf/cf(ini(1)); % bin frequencies start at approximately 0.33 with sca(ini(1))=1 exactly

% Middle
sca = sca(sca<10.5 & sca>0.5);  % restrict to 0.5 - 10.5 times fundamental

sca1 = sca;
filh = 1./(p.pefact-cos(2*pi*sca1));
filh = filh - sum(filh(1:end).*diff([sca1(1),(sca1(1:end-1)+sca1(2:end))./2,sca1(end)]))/sum(diff([sca1(1),(sca1(1:end-1)+sca1(2:end))./2,sca1(end)]));

posit = find(sca>=1);  % ==== this should just equal ini(1) ====
negat = find(sca<1);
numz = length(posit)-1-length(negat);
filh = filh./max(filh);
filh = [zeros(1,numz) filh]; % length is always odd with central value = 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filter the log-frequency scaled spectrogram
B = imfilter(O,filh);  % does a convolution with zero lag at centre of filh

% Feasible frequency range
numopt = p.numopt; % Number of possible fundamental frequencies per frame
flim = p.flim;
pfreq = find(cf>flim(1) & cf<flim(2)); % flim = permitted fx range = [60 400]
ff = zeros(nframes,numopt);
amp = zeros(nframes,numopt);
for i=1:nframes
    [pos,peak]=v_findpeaks(B(i,pfreq),[],5/(cf(pfreq(2))-cf(pfreq(1)))); % min separation = 5Hz @ fx=flim(1) (could pre-calculate) ====
    if numel(pos)
        [peak,ind]=sort(peak,'descend');
        pos = pos(ind);                     % indices of peaks in the B array
        posff = cf(pfreq(pos));             % frequencies of peaks
        fin = min(numopt,length(posff));
        ff(i,1:fin)=posff(1:fin);           % save both frequency and amplitudes
        amp(i,1:fin)=peak(1:fin);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Probabilitly of the frame of being voiced

% voiced/unvoiced decision based on 2-element feature vector
% (a) mean power of the frame's log-freq spectrum (normalized so its short-term average is LTASS)
% (b) sum of the power in the first three peaks

pow = mean(O,2);

vuvfea = [log(pow) 1e-3*sum(amp,2)./(pow+1.75*1e5)];

% %%%%%%%%%%%%%%%%%%%%%

pru=gaussmixp(vuvfea,m_u,v_u,w_u);  % Probability of being unvoiced
prv=gaussmixp(vuvfea,m_v,v_v,w_v);  % Probability of being voiced

pv=(1+exp(pru-prv)).^(-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dynamic programming

% w(1): relative amp, voiced local cost
% w(2): median pitch deviation cost
% w(3): df cost weight
% w(4): max df cost
% w(5): relative amp cost for missing peaks (very high)
% w(6): df mean

w = p.w;

% Relative amplitude
camp = -amp./repmat(max(amp,[],2),1,numopt);  % relative amplitude used as cost
camp(amp==0)=w(5); % If no frequency found

% Time interval for the median frequency
tmf = p.tmf; % in sec
inmf = round(tmf/txinc);

%--------------------------------------------------------------------------
% FORWARDS
% Initialize values
cost = zeros(nframes,numopt);
prev = zeros(nframes,numopt);
medfx = zeros(nframes,1);
dffact=2/txinc;

% First time frame
% cost(1,:) = w(1)*ramp(1,:);
cost(1,:) = w(1)*camp(1,:);  % only one cost term for first frame
fpos = ff(1:min(inmf,end),1);
mf=median(fpos(pv(1:min(inmf,end))>0.6));   % calculate median frequency of first 2 seconds
if isnan(mf)
    mf=median(fpos(pv(1:min(inmf,end))>0.5));
    if isnan(mf)
        mf=median(fpos(pv(1:min(inmf,end))>0.4));
        if isnan(mf)
            mf=median(fpos(pv(1:min(inmf,end))>0.3)); % ==== clumsy way of ensuring that we take the best frames ====
            if isnan(mf)
                mf=0;
            end
        end
    end
end
medfx(1)=mf;

for i=2:nframes              % main dynamic programming loop
    if i>inmf
        fpos = ff(i-inmf:i,1);  % fpos is the highest peak in each frame
        mf=median(fpos(pv(1:inmf)>0.6));  % find median frequency over past 2 seconds
        if isnan(mf)
            mf=median(fpos(pv(1:inmf)>0.5));
            if isnan(mf)
                mf=median(fpos(pv(1:inmf)>0.4));
                if isnan(mf)
                    mf=median(fpos(pv(1:inmf)>0.3));% ==== clumsy way of ensuring that we take the best frames ====
                    if isnan(mf)
                        mf=0;
                    end
                end
            end
        end
    end
    medfx(i)=mf;
    % Frequency difference between candidates and cost
    df = dffact*(repmat(ff(i,:).',1,numopt) - repmat(ff(i-1,:),numopt,1))./(repmat(ff(i,:).',1,numopt) + repmat(ff(i-1,:),numopt,1));
    costdf=w(3)*min((df-w(6)).^2,w(4));

    % Cost related to the median pitch
    if mf==0                                   % this test was inverted in the original version
        costf = zeros(1,numopt);
    else
        costf = abs(ff(i,:) - mf)./mf;
    end
    [cost(i,:),prev(i,:)]=min(costdf + repmat(cost(i-1,:),numopt,1),[],2); % ==== should we allow the possibility of skipping frames ? ====
    cost(i,:)=cost(i,:)+w(2)*costf + w(1)*camp(i,:);  % add on costs that are independent of previous path

end

% Traceback

fx=zeros(nframes,1);
ax=zeros(nframes,1);
best = zeros(nframes,1);

nose=find(cost(end,:)==min(cost(end,:))); % ==== bad method (dangerous) ===
best(end)=nose(1);
fx(end)=ff(end,best(end));
ax(end)=amp(end,best(end));
for i=nframes:-1:2
    best(i-1)=prev(i,best(i));
    fx(i-1)=ff(i-1,best(i-1));
    ax(i-1)=amp(i-1,best(i-1));
end

if nargout>=4
    fv.vuvfea=vuvfea;  % voiced-unvoiced features
    fv.best=best;  % selected path
    fv.ff=ff;  % pitch candidates
    fv.amp=amp;  % pitch candidate amplitudes
    fv.medfx=medfx;  % median pitch
    fv.w=w;  % DP weights
    fv.dffact=dffact;  % df scale factor
    fv.hist = [log(mean(O,2)) sum(amp,2)./((mean(O,2)))];
end

if ~nargout || any(m=='g') || any(m=='G')
    nax=0;  % number of axes sets to link
    msk=pv>0.5; % find voiced frames as a mask
    fxg=fx;
    fxg(~msk)=NaN; % allow only good frames
    fxb=fx;
    fxb(msk)=NaN; % allow only bad frames
    if any(m=='G') || ~nargout && ~any(m=='g')
        clf;
        spgrambw(s,fs,p.sopt); % draw spectrogram with log axes
        hold on
        plot(tx,log10(fxg),'-b',tx,log10(fxb),'-r'); % fx track
        yy=get(gca,'ylim');
        plot(tx,yy(1)+yy*[-1;1]*(0.02+0.05*pv),'-k'); % P(V) track
        hold off
        nax=nax+1;
        axh(nax)=gca;
        if any(m=='g')
            figure;   % need a new figure if plotting two graphs
        end
    end
    if any(m=='g')
        ns=length(s);
        [tsr,ix]=sort([(1:ns)/fs 0.5*(tx(1:end-1)+tx(2:end))']); % intermingle speech and frame boundaries
        jx(ix)=1:length(ix); % create inverse index
        sp2fr=jx(1:ns)-(0:ns-1);  % speech sample to frame number
        spmsk=msk(sp2fr);   % speech sample voiced mask
        sg=s;
        sg(~spmsk)=NaN;   % good speech samples only
        sb=s;
        sb(spmsk)=NaN;    % bad speech samples only
        clf;
        subplot(5,1,1);
        plot(tx,pv,'-b',(1:ns)/fs,0.5*mod(cumsum(fx(sp2fr)/fs),1)-0.6,'-b');
        nax=nax+1;
        axh(nax)=gca;
        ylabel('\phi(t), P(V)');
        set(gca,'ylim',[-0.65 1.05]);
        subplot(5,1,2:3);
        plot((1:ns)/fs,sg,'-b',(1:ns)/fs,sb,'-r');
        nax=nax+1;
        axh(nax)=gca;
        subplot(5,1,4:5);
        plot(tx,fxg,'-b',tx,fxb,'-r');
        ylabel('Pitch (Hz)');
        %         semilogy(tx,fxg,'-b',tx,fxb,'-r');
        %         ylabel(['Pitch (' yticksi 'Hz)']);
        set(gca,'ylim',[min(fxg)-30 max(fxg)+30]);
        nax=nax+1;
        axh(nax)=gca;
    end
    if nax>1
        linkaxes(axh,'x');
    end
end

function y=smooth(x,n)
nx=size(x,2);
nf=size(x,1);
c=cumsum(x,2);
y=[c(:,1:2:n)./repmat(1:2:n,nf,1) (c(:,n+1:end)-c(:,1:end-n))/n (repmat(c(:,end),1,floor(n/2))-c(:,end-n+2:2:end-1))./repmat(n-2:-2:1,nf,1)];

function y=timesm(x,n)
if ~mod(n,2)
    n = n+1;
end
nx=size(x,2);
nf=size(x,1);
c=cumsum(x,1);
mid = round(n/2);
y=[c(mid:n,:)./repmat((mid:n).',1,nx); ...
    (c(n+1:end,:)-c(1:end-n,:))/n; ...
    (repmat(c(end,:),mid-1,1) - c(end-n+1:end-mid,:))./repmat((n-1:-1:mid).',1,nx)];

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function option = mgoptions(option,N)
%% MGOPTIONS default options for multigrid solver
%
% x0 = 0;	tol = 1e-8;   solver = 'CG';    preconditioner = 'Vcycle';
% N0 = 500; coarsegridsolver = 'direct';    mu = 1;
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if isfield(option,'x0') && strcmpi(option.x0,'rand')
    option.x0 = rand(N,1);
end
if ~isfield(option,'x0')
    option.x0 = zeros(N,1);
end
if ~isfield(option,'tol')
    option.tol = 1e-8;
end
if ~isfield(option,'refType')
    option.refType = 'red';
end
if isfield(option,'solvermaxIt')
    option.solvermaxit = option.solvermaxIt;
elseif isfield(option,'solvermaxit')
    option.solvermaxit = option.solvermaxit;
else
    option.solvermaxit = min(N,200);
end
if isfield(option,'solver')
    option.solver = upper(option.solver);
else
    option.solver = 'CG';
end
if isfield(option,'outsolver')
    option.outsolver = upper(option.outsolver);
else
    option.outsolver = 'CG';
end
if isfield(option,'smoother')
    option.smoother = upper(option.smoother);
else
    option.smoother = 'GS';  % defacult one is Gauss-Seidel
    % option.solver = 'JAC' % Jacobi preconditioner
end
if ~isfield(option,'smoothingstep')  % smoothing steps
    option.smoothingstep = 1;
end
if ~isfield(option,'smoothingratio')  % ratio of variable smoothing
    option.smoothingratio = 1;
end
if ~isfield(option,'smoothingparameter')  % smoothing parameter
    option.smoothingparameter = 1;
end
if isfield(option,'preconditioner')
    option.preconditioner = upper(option.preconditioner);
elseif ~strcmp(option.solver,'Vcycle') && ~strcmp(option.solver,'Wcycle') ...
        && ~strcmp(option.solver,'Fcycle')
    option.preconditioner = 'V';
else
    option.preconditioner = 'NO';
end
if ~isfield(option,'N0')
    option.N0 = 50;
end
if ~isfield(option,'coarsegridsolver')
    option.coarsegridsolver = 'direct';    % solver in the coarsest grid
end
if ~isfield(option,'printlevel')
    option.printlevel = 1;
end
if ~isfield(option,'setupflag')
    option.setupflag = true;
end

repo: ifem


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function currentPose = inferenceReg(images,model,Pr,level,regsInfo)

m = length(images);
n_pts = size(model{level}.tpt,2) / 2;
tpt = model{max(1,level-1)}.tpt;
mt = size(tpt,1);

% 61. Sampling
index_inference = zeros(m,regsInfo.samplingTot(level));
currentPose_inference = zeros(m,2*n_pts,regsInfo.samplingTot(level));
for i = 1:m
    if level > 1
        [~,tmp] = sort(Pr(i,:),'descend');
        index_inference(i,:) = tmp(1:regsInfo.samplingTot(level));
    else
        index_inference(i,:) = randSampleNoReplace(mt, regsInfo.samplingTot(level), Pr(i,:));
    end;
    for j = 1:regsInfo.samplingTot(level)
        currentPose_inference(i,:,j) = tpt(index_inference(i,j),:);
    end;
end;

% 62. Inference Iteration
featInfo.scale = regsInfo.SIFTscale;
for iter = 1:regsInfo.iterTot(level)
    featOri_batch = extractSIFTs_toosimple_samples(images,currentPose_inference,iter,featInfo);
    for j = 1:regsInfo.samplingTot(level)
        currentPose_inference(:,:,j) = currentPose_inference(:,:,j) + ...
            (featOri_batch(:,:,j) - repmat(model{level}.reg{iter}.mu,m,1)) * model{level}.reg{iter}.A ...
             + repmat(model{level}.reg{iter}.b,m,1);
    end;
end;

currentPose = poseVoting(currentPose_inference,regsInfo.dominantIterTot(level));

end

author: zhusz


text:
 function x = eml_psca(x, G, yi, ci, ri, nsubiter)
%function x = eml_psca(x, G, yi, ci, ri, nsubiter)
% Runs one iteration of the ML-PSCA algorithm for emission Poisson problem
% (paraboloidal surrogates coordinate ascent)
% model: Y_i ~ Poisson(c_i [G x]_i + r_i)
% in
%	see em_fbp.m for model, G, yi, ci, ri
% out
%	x [np,1]	updated image
%
% This is for TESTING ONLY.
% IT IS NOT USEFUL because it is unregularized and slow
%
% Copyright Mar 2000, Jeff Fessler, The University of Michigan

if nargin < 3, ir_usage, end

[nb na] = size(yi);

if ~isvar('ci') || isempty(ci)
	ci = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end
if ~isvar('nsubiter') || isempty(nsubiter)
	nsubiter = 1;
end

%
% compute surrogate curvatures
%
li0 = reshape(G * x(:), size(yi));		% l=G*x "line integrals"
yb = ci .* li0 + ri;				% predicted measurement means

doth = ci .* (yi ./ yb - 1);			% first derivative
ni = eml_curvature(yi, ci, ri, li0, yb, 'oc');	% curvatures
ni = ni(:);
%disp(range(ni)')

% backproject curvatures for denominator
dj = (G').^2 * ni(:);

%
% now we are maximizing the quadratic surrogate function:
% Q(x) == P(G*x ; G*x0)
%	where P(l;l0) = doth' (l-l0) - 1/2 (l-l0)' D(n) (l-l0)
% d/dli P(l;l0) = dothi - ni (li-l0i)
% Q(x) == doth' G (x-x0) - 1/2 (x-x0)' G' D(n) G (x-x0)
% ??? d/dx_j Q(x) = \sumi \gij \dothi - \sumi \gij \ni [G(x-x0)]_i
% d/dx_j Q(x) = \sumi \gij [d/dli P(l;l0)]
% -d^2/dx_j^2 Q(x) = \sumi \gij^2 \ni
%
dqi = doth(:);	% initial surrogate derivatives

%xs = zeros([length(x) nsubiter]);	% save and return the subiterations too

%x0 = x;		% save initial image

like0 = eql_obj(x, G, yi(:), ci(:), ri(:));

%
% loop over subiterations of paraboloid
%
for is=1:nsubiter

	%
	% loop over pixels
	%
	for jj=1:numel(x)
		g = G(:,jj);
		x0j = x(jj);
		x(jj) = x0j + (dqi' * g) / dj(jj);
		x(jj) = max(x(jj),0);
		dqi = dqi - ni .* (g * (x(jj) - x0j));

	end
%	xs(:,is) = x;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% works fine, timestepping can be improved (ODE solver?)
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
%% define the geometry
Nx = 100; % number of cells in x direction
Ny = 30; % number of cells in y direction
W = 300; % [m] length of the domain in x direction
H = 30; % [m] length of the domain in y direction
m = createMesh1D(Nx, W); % creates a 1D mesh
%% define the physical parametrs
krw0 = 1.0;
kro0 = 0.76;
nw = 2.4;
no = 2.0;
sor=0.12;
swc=0.09;
sws=@(sw)((sw>swc).*(sw<1-sor).*(sw-swc)/(1-sor-swc)+(sw>=1-sor).*ones(size(sw)));
kro=@(sw)((sw>=swc).*kro0.*(1-sws(sw)).^no+(sw<swc).*(1+(kro0-1)/swc*sw));
krw=@(sw)((sw<=1-sor).*krw0.*sws(sw).^nw+(sw>1-sor).*(-(1-krw0)/sor.*(1.0-sw)+1.0));
dkrwdsw=@(sw)((sw<=1-sor).*nw.*krw0.*(1/(1-sor-swc)).*sws(sw).^(nw-1)+(sw>1-sor)*((1-krw0)/sor));
dkrodsw=@(sw)((sw>=swc).*(-kro0*no*(1-sws(sw)).^(no-1))/(-swc-sor+1)+(sw<swc).*((kro0-1)/swc));
p0 = 100e5; % [bar] pressure
pin = 150e5; % [bar] injection pressure at the left boundary
u_in= 1.0/(24*3600); % [m/s] equal to 1 m/day
sw0 = swc+0.1; % initial water saturation
sw_in = 1;
mu_oil = 2e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 2e-12; % [m^2] average reservoir permeability
phi0 = 0.2; % average porosity
clx=1.2;
cly=0.2;
V_dp=0.7; % Dykstra-Parsons coef.
perm_val= k0; %field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
%% Define the boundaries
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% left boundary pressure gradient
BCp.left.a(:)=(krw(sw_in)*lw.xvalue(1,:)+kro(sw_in)*lo.xvalue(1,:)); BCp.left.b(:)=0; BCp.left.c(:)=-u_in;
% change the right boandary to constant pressure (Dirichlet)
% BCp.left.a(:)=0; BCp.left.b(:)=1; BCp.left.c(:)=pin;
BCp.right.a(:)=0; BCp.right.b(:)=1; BCp.right.c(:)=p0;
% change the left boundary to constant saturation (Dirichlet)
BCs.left.a(:)=0; BCs.left.b(:)=1; BCs.left.c(:)=1;
%% define the time step and solver properties
% dt = 1000; % [s] time step
dt=(W/Nx)/u_in/10; % [s]
t_end = 1000*dt; % [s] final time
eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-5; % saturation accuracy
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
p = p_old;
uw = -gradientTerm(p_old); % an estimation of the water velocity
%% start the main loop
% generate intial pressure profile (necessary to initialize the fully
% implicit solver)
dp_alwd= 100.0; % Pa
dsw_alwd= 0.05;
t = 0;
% fprintf(1, 'progress (%%):  ');
while (t<t_end)
    error_p = 1e5;
    error_sw = 1e5;
    % Implicit loop
    loop_count=0;
    while ((error_p>eps_p) || (error_sw>eps_sw))
        loop_count=loop_count+1;
        if loop_count>10
            break
        end
        % calculate parameters
        pgrad = gradientTerm(p);
        sw_face = upwindMean(sw, -pgrad); % average value of water saturation
        labdao = lo.*funceval(kro, sw_face);
        labdaw = lw.*funceval(krw, sw_face);
        dlabdaodsw = lo.*funceval(dkrodsw, sw_face);
        dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
        labda = labdao+labdaw;
        dlabdadsw = dlabdaodsw+dlabdawdsw;
        % compute [Jacobian] matrices
        Mdiffp1 = diffusionTerm(-labda);
        Mdiffp2 = diffusionTerm(-labdaw);
        Mconvsw1 = convectionUpwindTerm(-dlabdadsw.*pgrad);
        Mconvsw2 = convectionUpwindTerm(-dlabdawdsw.*pgrad);
        [Mtranssw2, RHStrans2] = transientTerm(sw_old, dt, phi);
        % Compute RHS values
        RHS1 = divergenceTerm(-dlabdadsw.*sw_face.*pgrad);
        RHS2 = divergenceTerm(-dlabdawdsw.*sw_face.*pgrad);
        % include boundary conditions
        [Mbcp, RHSbcp] = boundaryCondition(BCp);
        [Mbcsw, RHSbcsw] = boundaryCondition(BCs);
        % Couple the equations; BC goes into the block on the main diagonal
        M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw];
        RHS = [RHS1+RHSbcp; RHS2+RHStrans2+RHSbcsw];
        % solve the linear system of equations
        x = M\RHS;
        % x = agmg(M, RHS, [], 1e-10, 500, [], [p.value(:); sw.value(:)]);
        % separate the variables from the solution
        p_new = reshapeCell(m,full(x(1:(Nx+2))));
        sw_new = reshapeCell(m,full(x((Nx+2)+1:end)));
        % calculate error values
        error_p = max(abs((p_new(:)-p.value(:))./p_new(:)));
        error_sw = max(abs(sw_new(:)-sw.value(:)));
        % assign new values of p and sw
        p.value = p_new;
        sw.value = sw_new;
    end
    if loop_count>10
      p=p_old;
      sw=sw_old;
      dt=dt/5;
      continue
    end
    dsw=max(abs(sw_new(:)-sw_old.value(:))./sw_new(:));
    t=t+dt;
    fprintf(1,'\b\b\b\b\b\b\b\b\b\b\b\b\b\bProgress: %d %%',floor(t/t_end*100));
    dt=min([dt*(dsw_alwd/dsw), 2*dt, t_end-t]);
    p_old = p;
    sw_old = sw;
    figure(1);visualizeCells(sw); drawnow;
end
fprintf('\n')

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
% absolute difference in orientation between a fly and the closest fly
% according to type
function [data,units] = compute_absthetadiff(trx,n,type)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);

for i1 = 1:nflies,
  fly1 = flies(i1);

  % fly closest to fly1 according to type
  closestfly = trx(fly1).(['closestfly_',type]);

  % orientation of fly1
  theta_mm1 = trx(fly1).theta_mm;

  % loop over all flies
  for i2 = 1:nflies,

    fly2 = flies(i2);
    if i1 == i2, continue; end

    % frames where this fly is closest
    idx = find(closestfly == fly2);
    if isempty(idx), continue; end

    % orientation of fly2
    off = trx(fly1).firstframe - trx(fly2).firstframe;
    theta_mm2 = trx(fly2).theta_mm(off+idx);

    % absolute difference in orientation
    data{i1}(idx) = abs(modrange(theta_mm2 - theta_mm1(idx),-pi,pi));

  end
end

units = parseunits('rad');

author: kristinbranson

repo: JAABA


text:
function OUT = test_pattern_on_kragel_2018_n270_data(varargin)
%
% test a multivariate signature pattern stored in an fmri_data object on Kragel et al. 2018 datasets, N = 270
%
% -   requires kragel_2018_nat_neurosci_270_subjects_test_images.mat
%     if not found, will attempt to download from Neurovault using
%     retrieve_neurovault_collection(). See load_image_set('kragel18_alldata')
%
% -   Dataset is from Kragel et al. 2018, Nature Neuroscience.
%     Uses load_image_set to load it:
%     test_images = load_image_set('kragel18_alldata', 'noverbose');
%
% -   If not found, load_image_set uses the utility retrieve_neurovault_collection
%     to pull the images from Neurovault.org if you don't have them.
%     [files_on_disk, url_on_neurovault, mycollection, myimages] = retrieve_neurovault_collection(504);
%     image_obj = fmri_data(files_on_disk);
%     It also adds meta-data to label the images.
%
% :Usage:
% ::
%
%     [list outputs here] = test_pattern_on_kragel_2018_n270_data(obj)
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2019 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **obj:**
%        An fmri_data object storing a single image of multivariate pattern weights
%
% :Optional Inputs:
%   **'noverbose':**
%        Turn off verbose output
%
%   **'noplots':**
%       Turn off plots
%
%   **'threshold_type'
%       Followed by keyphrase, e.g., 'Optimal balanced error rate'
%       Default is 'Optimal overall accuracy'
%
%
% :Outputs:
%
%   **OUT:**
%        A structure with a summary table of sensivity, specificity, PPV,
%        and effect size for each classification (pain vs. other, cog
%        control vs. other, negative emotion vs. other.
%        Also includes roc_plot output for all three classifications, with
%        more stats and function handles.
%
%
% :Examples:
% ::
%
%    % give examples of code here
%    param1 = abc();
%    param2 = xyz();
%    [out1,out2] = func_call(param1, param2)
%
% :References:
%   CITATION(s) HERE
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
% ..

% ----------------------------------------------------------------------
% Defaults
% ----------------------------------------------------------------------
% Must define ALL inputs here. All these can be changed by passing in arguments.

obj = [];
num_cols = 18;      % this is for the test dataset
doverbose = true;
doplot = true;
fontsize = 14;
threshold_type = 'Optimal balanced error rate'; % 'Optimal overall accuracy'; %

colors = [repmat({[1 .2 0]}, 1, 6) repmat({[.4 .5 .2]}, 1, 6) repmat({[.1 0 .9]}, 1, 6)];

sim_string = 'cosine_similarity';  % or 'dotproduct' or 'correlation'
weight_string = 'none';  % or 'weighted'

canlab_parse_inputs_subfcn();

% ----------------------------------------------------------------------
% Display helper functions: Called by later scripts
% ----------------------------------------------------------------------

dashes = '----------------------------------------------';
printstr = @(dashes) disp(dashes);
printhdr = @(str) fprintf('%s\n%s\n%s\n', dashes, str, dashes);

disp('Process: test_pattern_on_kragel_2018_n270_data');
printstr(dashes)

disp(' ');

% ----------------------------------------------------------------------
% Load and prep test data
% ----------------------------------------------------------------------
test_images = load_image_set('kragel18_alldata', 'noverbose');

test_images.Y = test_images.metadata_table.Studynumber;

% vector_data = apply_mask(test_images, obj, 'pattern_expression', sim_string);
vector_data = apply_mask(test_images, obj, 'pattern_expression', sim_string,weight_string); %weight by overall cosine similarity

for s = 1:length(test_images.additional_info)
    pattern_response{s} = vector_data(test_images.Y==s);
end

% clear exp
% labels = test_images.additional_info;
labels = format_strings_for_legend(test_images.additional_info);

% ----------------------------------------------------------------------
% ROC and stats
% ----------------------------------------------------------------------

% [ispain, iscog, isemo] = deal(false(270, 1));
%
% ispain(1:6*15) = true;
% iscog(6*15 + 1 : (12*15)) = true;
% isemo(12*15+1:18*15) = true;

ispain = strcmp(test_images.metadata_table.Domain, 'Pain');
iscog = strcmp(test_images.metadata_table.Domain, 'Cog_control');
isemo = strcmp(test_images.metadata_table.Domain, 'Neg_Emotion');

wh_isnan=isnan(vector_data);
vector_data=vector_data(~wh_isnan);

ispain=ispain(~wh_isnan);
iscog=iscog(~wh_isnan);
isemo=isemo(~wh_isnan);

% Boostrap confidence intervals for effect size
fprintf('Bootstrapping CIs\n')

bci_pain = bootci(5000, @cohens_d_2sample, vector_data, ispain);

bci_cog = bootci(5000, @cohens_d_2sample, vector_data, iscog);

bci_emo = bootci(5000, @cohens_d_2sample, vector_data, isemo);

% ----------------------------------------------------------------------
% ROC calculation and classification output
% ----------------------------------------------------------------------

printhdr('Pattern classification results')
disp('Applied pattern to test data from 18 studies, n = 270 participants')
disp('From Kragel et al. 2018, Nature Neuroscience')
disp(' ')

disp('Classification of pain vs. other')

roc_pain = roc_plot(vector_data,ispain, 'noplot', 'threshold_type', threshold_type);

fprintf('Threshold for cosine_sim = %3.4f, Cohen''s d(pain vs no) = %3.2f, Bootstrapped CI is [%3.2f %3.2f]\n', ...
    roc_pain.class_threshold, cohens_d_2sample(vector_data, ispain), bci_pain(1), bci_pain(2));

disp(' ')

disp('Classification of cognitive control vs. other')

roc_cog = roc_plot(vector_data,iscog, 'noplot', 'threshold_type', threshold_type);

fprintf('Threshold for cosine_sim = %3.4f, Cohen''s d(cog vs no) = %3.2f, Bootstrapped CI is [%3.2f %3.2f]\n', ...
    roc_cog.class_threshold, cohens_d_2sample(vector_data, iscog), bci_cog(1), bci_cog(2));

disp(' ')

disp('Classification of emotion vs. other')

roc_emo = roc_plot(vector_data,isemo, 'noplot', 'threshold_type', threshold_type);

fprintf('Threshold for cosine_sim = %3.4f, Cohen''s d(emo vs no) = %3.2f, Bootstrapped CI is [%3.2f %3.2f]\n', ...
    roc_emo.class_threshold, cohens_d_2sample(vector_data, isemo), bci_emo(1), bci_emo(2));
disp(' ')

% ----------------------------------------------------------------------
% Output table
% ----------------------------------------------------------------------

tabledata(1, :) = [roc_pain.class_threshold cohens_d_2sample(vector_data, ispain) roc_pain.sensitivity roc_pain.specificity roc_pain.PPV ];
tabledata(2, :) = [roc_cog.class_threshold cohens_d_2sample(vector_data, iscog) roc_cog.sensitivity roc_cog.specificity roc_cog.PPV ];
tabledata(3, :) = [roc_emo.class_threshold cohens_d_2sample(vector_data, isemo) roc_emo.sensitivity roc_emo.specificity roc_emo.PPV ];

T = array2table(tabledata, 'VariableNames', {'Cos_Sim_Thresh' 'Cohens_d' 'Sensitivity' 'Specificity' 'PPV'}, 'RowNames', {'Pain' 'Cog Control' 'Neg Emotion'});

disp(T);

OUT.summary_table = T;
OUT.roc_pain = roc_pain;
OUT.roc_cog = roc_cog;
OUT.roc_emo = roc_emo;

% ----------------------------------------------------------------------
% Plot
% ----------------------------------------------------------------------
if doplot

    create_figure('plot', 1, 3);

    clear hh

    subplot(1,3,3)

    roc = roc_plot(vector_data,ispain, 'threshold_type', threshold_type, 'nooutput');
    set(roc.line_handle, 'Color', 'r');
    hh(1) = roc.line_handle(2);

    roc = roc_plot(vector_data,iscog, 'threshold_type', threshold_type, 'nooutput');
    set(roc.line_handle, 'Color', 'g');
    hh(2) = roc.line_handle(2);

    roc = roc_plot(vector_data,isemo, 'threshold_type', threshold_type, 'nooutput');
    set(roc.line_handle, 'Color', 'b');
    hh(3) = roc.line_handle(2);

    set(gca,'fontsize', fontsize)
    legend(hh, {'Pain vs Other' 'Cog vs. Other' 'Emo vs. Other'});

    printhdr('Stats for individual studies')

    subplot(1,3,1:2)

    set(gca, 'FontSize', fontsize);
    barplot_columns(pattern_response,'names',labels,'colors',colors,'nofig');
    ylabel 'Pattern Response (Cosine Similarity)'
    hold on;

    if roc_pain.accuracy_p < .05
        plot([0 19],[roc_pain.class_threshold roc_pain.class_threshold],'r--')
    end

    if roc_cog.accuracy_p < .05
        plot([0 19],[roc_cog.class_threshold roc_cog.class_threshold],'g--')
    end

    if roc_emo.accuracy_p < .05
        plot([0 19],[roc_emo.class_threshold roc_emo.class_threshold],'b--')
    end

else
    % no plot, but print stats

    printhdr('Stats for individual studies')
    barplot_columns(pattern_response,'names',labels, 'colors',colors, 'skipallplots');

end % plot

% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
% Nested functions
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------

    function canlab_parse_inputs_subfcn()
        % INPUT PARSER: This is a nested function.
        % It modifies variables in the calling workspace directly.
        % It validates inputs.
        % It uses inputParser, but also allows for special keywords that change
        % variables, which are not handled by inputParser.

        % ----------------------------------------------------------------------
        % Parse inputs
        % ----------------------------------------------------------------------
        p = inputParser;

        % Validation functions - customized for each type of input
        % ----------------------------------------------------------------------

        valfcn_scalar = @(x) validateattributes(x, {'numeric'}, {'nonempty', 'scalar'});
        valfcn_number = @(x) validateattributes(x, {'numeric'}, {'nonempty'}); % scalar or vector
        valfcn_cell = @(x) validateattributes(x, {'cell'}, {'nonempty'}, {'numel'}, num_cols);

        valfcn_image_vector = @(x) isa(x, 'image_vector');

        % Required inputs
        % ----------------------------------------------------------------------
        p.addRequired('obj', valfcn_image_vector);

        % Optional inputs - key-value pairs
        % ----------------------------------------------------------------------
        % Pattern: keyword, value, validation function handle

        p.addParameter('fontsize', fontsize);
        p.addParameter('colors', colors);
        p.addParameter('sim_string', sim_string);
        p.addParameter('weight_string',weight_string);
        p.addParameter('threshold_type', threshold_type);

        % Optional inputs - Logical flags and keywords
        % ----------------------------------------------------------------------
        % Identify and remove these before parsing, because they do not follow name-value pair pattern
        wh = strcmp(varargin, 'noverbose');
        if any(wh), doverbose = false; varargin(wh) = []; end

        wh = strcmp(varargin, 'noplots');
        if any(wh), doplot = false; varargin(wh) = []; end

        wh = strcmp(varargin, 'noplot');
        if any(wh), doplot = false; varargin(wh) = []; end

        % Parse inputs and distribute out to variable names in workspace
        % ----------------------------------------------------------------------
        p.parse(varargin{:});

        IN = p.Results;
        fn = fieldnames(IN);

        for i = 1:length(fn)
            str = sprintf('%s = IN.(''%s'');', fn{i}, fn{i});
            eval(str)
        end

    end  % Parse subfunction

end % main function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore