text:
function H = get_shelve_lagrange(f,H,FlagSub,fSub,FlagAliasing,fAliasing,Bandwidth_in_Oct)
%GET_SHELVE_LAGRANGE Lagrange interpolation towards shelving filter
%
%   Usage: H = get_shelve_lagrange(f,H,FlagSub,fSub,FlagAliasing, ...
%                                  fAliasing,Bandwidth_in_Oct)
%
%   Input parameters:
%       f                -  frequency vector in Hz, typical 0 to half
%                           sampling frequency, equidistant sampling is
%                           assumed (i.e. DFT frequencies) but not required
%       H                -  complex spectrum (mag/phase) at specified frequencies
%                           with meaningful slope, typically used for +3dB/oct.
%       FlagSub          -  use low-shelf part, 0 or 1
%       fSub             -  cut-frequency in Hz for low-shelf,
%                           check e.g. fLow=100Hz in [Fig2, Sch13]
%       FlagAliasing     -  use high-sehlf part, 0 or 1
%       fAliasing        -  cut-frequency in Hz for high-shelf,
%                           check e.g. fAliasing=2kHz in [Fig2, Sch13]
%       Bandwidth_in_Oct -  interpolation bandwidth, i.e. shelf knee range
%                           in octaves, allowed: 0.5, 1, 2, 3, 4
%
%   Output parameters:
%       H                -  return complex spectrum at specified frequencies
%                           after low-/ and high-shelf interpolation              -
%
%   GET_SHELVE_LAGRANGE(f,H,FlagSub,fSub,FlagAliasing,fAliasing,Bandwidth_in_Oct)
%   does an Lagrange interplation to get a shelving filter. This function is used
%   in wfs_iir_prefilter(conf).
%
%   Note 1: there is a analytical expression for arbitrary Bandwidth_in_Oct
%   which is not implemented yet
%   Note 2: the Offset_dB specified for Bandwidth_in_Oct work only for the
%   3dB/oct. slope, other driving functions may have other slopes and
%   different interpolation offsets may be required
%
%   See also: wfs_iir_prefilter
%
%   References:
%       Schultz, Erbes, Spors, Weinzierl (2013) - "Derivation of IIR prefilters
%       for soundfield synthesis using linear secondary source distributions",
%       International Conference on Acoustics (AIA-DAGA), p.2372-2375,
%       http://pub.dega-akustik.de/AIA_DAGA_2013/data/articles/000604.pdf

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
% Revision: 07/02/2013 frank.schultz@uni-rostock.de initial development      *
%*****************************************************************************
    Hphase = unwrap(angle(H));  %save original phase
    H = 20*log10(abs(H));       %interpolation in dB

    if Bandwidth_in_Oct==4
        Offset_dB = 1.5053;
    elseif Bandwidth_in_Oct==3
        Offset_dB = 1.1286;
    elseif Bandwidth_in_Oct==2
        Offset_dB = 0.7528;
    elseif Bandwidth_in_Oct==1
        Offset_dB = 0.3766;
    elseif Bandwidth_in_Oct==0.5
        Offset_dB = 0.1857;
    else
        disp('error, Bandwidth_in_Oct must be 1,2,3 or 4')
    end

    if FlagSub
        %Sub-Bass limitation:
        tmp=(find(f<=fSub)); fSub_i = tmp(end);
        H_Sub = H*0+H(fSub_i);
        H_Sub(1) = H_Sub(2);
    end

    if FlagAliasing
        % Aliasing frequency limitation:
        tmp=(find(f<=fAliasing)); falias_i = tmp(end);
        H_Alias = H*0+H(falias_i);
        H_Alias(1) = H_Alias(2);
    end

    if FlagSub
        % Lagrange interpolation for SUB:
        fl = f(fSub_i)*2^(-Bandwidth_in_Oct/2);
        fh = f(fSub_i)*2^(+Bandwidth_in_Oct/2);
        tmp=(find(f<=fl)); fl_i = tmp(end);
        tmp=(find(f<=fh)); fh_i = tmp(end);

        P1x = log10(f(fl_i));
        P1y = H_Sub(fl_i);
        P2x = log10(f(fSub_i));
        P2y = H_Sub(fSub_i)+Offset_dB;
        P3x = log10(f(fh_i));
        P3y = H(fh_i);

        % Lagrange interpolation tmp variables:
        ipol1=P1y/((P1x-P2x)*(P1x-P3x));
        ipol2=P2y/((P2x-P1x)*(P2x-P3x));
        ipol3=P3y/((P3x-P1x)*(P3x-P2x));

        for fi=1:fl_i-1
        H(fi) = H_Sub(fi);
        end

        for fi=fl_i:fh_i
            H(fi)=          (ipol1*(log10(f(fi))-P2x)*(log10(f(fi))-P3x))+...
                            (ipol2*(log10(f(fi))-P1x)*(log10(f(fi))-P3x))+...
                            (ipol3*(log10(f(fi))-P1x)*(log10(f(fi))-P2x));
        end
    end

    if FlagAliasing
        % Lagrange interpolation for ALIASING:
        fl = f(falias_i)*2^(-Bandwidth_in_Oct/2);
        fh = f(falias_i)*2^(+Bandwidth_in_Oct/2);
        tmp=(find(f<=fl)); fl_i = tmp(end);
        tmp=(find(f<=fh)); fh_i = tmp(end);

        P1x = log10(f(fl_i));
        P1y = H(fl_i);
        P2x = log10(f(falias_i));
        P2y = H_Alias(falias_i)-Offset_dB;
        P3x = log10(f(fh_i));
        P3y = H_Alias(fh_i);

        % Lagrange interpolation tmp variables:
        ipol1=P1y/((P1x-P2x)*(P1x-P3x));
        ipol2=P2y/((P2x-P1x)*(P2x-P3x));
        ipol3=P3y/((P3x-P1x)*(P3x-P2x));

        for fi=fl_i:fh_i
            H(fi)=      (ipol1*(log10(f(fi))-P2x)*(log10(f(fi))-P3x))+...
                        (ipol2*(log10(f(fi))-P1x)*(log10(f(fi))-P3x))+...
                        (ipol3*(log10(f(fi))-P1x)*(log10(f(fi))-P2x));
        end

        for fi=fh_i+1:length(H)
            H(fi) = H_Alias(fi);
        end
    end

    %apply original phase and delog
    H = 10.^(H/20).*exp(1i*Hphase);
    H(1) = abs(H(2));
    %we ignore the special treatment for DC and fs/2 as would be usual for
    %DFT data, since we need a 'classical' frequency response
end

author: sfstoolbox

repo: sfs-matlab


text:
function io_save(fname, varargin)
% Like MATLAB's builtin save(), except that the fname can be platform-independent.
% io_save(Filename, Arguments...)
%
% Also, additional -options are supported, which offer additional fault-tolerance and automation,
% especially for long-running batch scripts.
%
% In:
%   Filename    :   Platform-Independent file name. If no ending is given, .mat will be assumed.
%                   If the ending is .sto, then BCILAB's Storage file format will be used (10-50x as
%                   fast as .mat).
%
%   Arguments...:   optional arguments to save; in addition to what is supported by save(), the following options are allowed:
%                   -makedirs   : make directories (recursively), if possible
%                   -retryinput : ask the user whether he/she wants to retry creating the directory or the file, if unsuccessful (using input())
%                   -attributes attribs : specify file/directory attributes; these are the second or second and third argument to be passed to
%                                         the fileattrib() function, separated by a comma, to set attributes for the current user or specific user groups,
%                                         e.g., '+w' or '+w','a'.
%                   -prunehandles : prune unreferenced variables from anonymous function handles; these are normally not used but saved anyway since cases
%                                   can be constructed in which they are actually needed by the function (e.g. using evalin('caller',...))
%                                   since these can be extremely large, they can easily cause save() to fail
%                   -serialized : compress the data using hlp_serialize
%                                 (when loading with io_load, the data will be automatically uncompressed,
%                                 not so when using load directly)
%                   -dataformat format : specify the data format to use for saving; the default is
%                                        'v7.3', but for interop with some software, 'v7' is needed
%
% Example:
%   % save the variables a,b, and c to a .mat file
%   io_save('store:/myoutput.mat', 'a','b','c');
%
%   % as before, but try to create directories if necessary
%   io_save('store:/myoutput.mat', 'a','b','c','-makedirs');
%
%   % like before, but this time prompt the user to re-try if an error occurs
%   io_save('store:/myoutput.mat', 'a','b','c','-makedirs','-retryinput');
%
%   % like before, but this time make the created file writable for all users
%   io_save('store:/myoutput.mat', 'a','b','c','-attributes','''+w'',''a''');
%
% See also:
%   save, io_load
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-06-21

% translate the file name
fname = env_translatepath(fname);
if length(fname)<4 || ~any(fname(end-3:end-1)=='.')
    fname = [fname '.mat']; end

% parse arguments...
saveargs = {};        % default args to save()
variables = {};       % names of variables to save (subset of saveargs)
makedirs = false;     % whether to create directories
retryinput = false;   % whether to ask the user for a retry
prunehandles = false; % whether to prune unreferenced workspace variables from function handles in arguments
serialized = false;   % whether to write variables as serialized
nooverwrite = false;  % whether to disable overwriting existing files
dataformat = 'v7.3';  % data format to use
fileattriblist = [];  % the list of attributes for files & directories
i = 1;
while i <= length(varargin)
    arg = varargin{i};
    if ischar(arg) && ~isempty(arg) && arg(1) == '-'
        % we are dealing with an option argument
        switch(strtrim(arg))
            case '-makedirs'
                makedirs = true;
            case '-retryinput'
                retryinput = true;
            case '-attributes'
                fileattriblist = varargin{i+1};
                i = i+1; % skip next arg, too
            case '-prunehandles'
                prunehandles = true;
            case {'-serialized','-serialize'}
                serialized = true;
            case '-nooverwrite'
                nooverwrite = true;
            case '-dataformat'
                nextarg = varargin{i+1};
                i = i+1;
                if ~strcmp(nextarg, 'default')
                    if nextarg(1) ~= 'v'
                        disp('WARNING: io_save dataformat is expected to begin with a "v" character.'); end
                    dataformat = nextarg;
                end
            otherwise
                % regular argument; append it
                saveargs{end+1} = arg;
        end
    elseif ~isempty(arg)
        % append the argument
        saveargs{end+1} = arg;
        variables{end+1} = arg;
    end
    i = i+1; % go to next argument
end

% prepend save format option
saveargs = [{['-' dataformat]} saveargs];

% move variables into a local workspace for optional processing
if isempty(variables)
    % no variables listed: obtain all
    varinfos = evalin('caller','whos');
    variables = {varinfos.name};
    saveargs = [saveargs variables];
end
for k=1:length(variables)
    wkspace.(variables{k}) = evalin('caller',variables{k}); end

% prune their handles
if prunehandles
    for k=1:length(variables)
        wkspace.(variables{k}) = utl_prune_handles(wkspace.(variables{k})); end
end

% serialize the variables
if serialized
    for k=1:length(variables)
        wkspace.(variables{k}) = hlp_serialize(wkspace.(variables{k})); end
    wkspace.is_serialized__ = ['variables were serialized using ' hlp_funcversion(@hlp_serialize) '; use hlp_deserialize (see File Exchange) to return them to their original form.'];
    saveargs{end+1} = 'is_serialized__';
end

% reformat the save-args into a string
tmp = [saveargs' repmat({' '},length(saveargs),1)]';
saveargs = [tmp{:}];

% reformat the fileattriblist into a cell array...
fileattriblist = eval(['{' fileattriblist '}']);

while 1
    if makedirs
        % first ensure that the directory exists
        while 1
            try
                io_mkdirs(fname,fileattriblist);
                break;
            catch
                warning('BCILAB:io_save:cannot_create_directory',['Error creating target directory ' fname '. Please review your permissions.']);
                if retryinput
                    if strcmp(input('Would you like to retry? [yes/no]: ','s'),'no')
                        disp('Did not save the data.');
                        break;
                    end
                else
                    break;
                end
            end
        end
    end

    % we have copied the variables to be saved into the 'wkspace' variable; now save
    try
        if nooverwrite && exist(fname,'file')
            return; end
        if length(fname)>4 && strcmpi(fname(end-3:end),'.sto')
            % save as .sto file (fast)
            bytes = hlp_serialize(wkspace);
            f = fopen(fname,'w');
            fwrite(f,bytes);
            fclose(f);
        else
            % save as .mat file (slow)
            save_workspace(fname,saveargs,wkspace);
        end
        if ~isempty(fileattriblist)
            % try to set file fileattriblist...
            if ~isempty(fileattriblist)
                warning off MATLAB:FILEATTRIB:SyntaxWarning
                fileattrib(fname,fileattriblist{:});
            end
        end
        break;
    catch e
        warning('BCILAB:io_save:cannot_save_file',['Error saving target file ' fname '. Please review your permissions: ' e.message]);
        if retryinput
            if strcmp(input('Would you like to retry? [yes/no]','s'),'no')
                disp('Did not save the data.');
                break;
            end
        else
            break;
        end
    end
end

% save, given a workspace of variables
function save_workspace(fname__,saveargs__,workspace__)
for fn__=fieldnames(workspace__)'
    eval([fn__{1} ' = workspace__.(fn__{1});']); end
clear workspace__ fn__;
eval(['save ''' fname__ ''' ' saveargs__ ';']);

repo: REST


text:
function r = invwishrand(S,nu);
%INVWISHRND Random matrices from inverse Wishard distribution.
%
%   R = INVWISHRAND(S,N) returns a matrix of random numbers chosen
%   from the central inverse Wishard distribution with parameters S and NU.
%
%   S is a symmetric positive definite scale matrix
%   NU is degrees of freedom
%
%   Note: E[R]=S*nu/(nu-k-1)
%
%	See also WISHRAND
%
% Copyright (c) 1999 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 2
  error('Requires two input arguments.');
end;

[d d2] = size(S);
if d ~= d2
  error('Matrix S must be square');
end

[t,p]=chol(S);
if p < 0
  error('Matrix S must be positive definite.');
end

r=inv(wishrand(inv(S),nu));

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [fx,tt]=v_fxrapt(s,fs,mode,q)
%V_FXRAPT RAPT pitch tracker [FX,VUV]=(S,FS,M,Q)
%
% Input:   s(ns)      Speech signal
%          fs         Sample frequency (Hz)
%          mode       'g' will plot a graph [default if no output arguments]
%                     'u' will include unvoiced fames (with fx=NaN)
%          q          stucture with parameter values (e.g. q.f0min=40); see below for a list
%
% Outputs: fx(nframe)     Larynx frequency for each fram,e (or NaN for silent/unvoiced)
%          tt(nframe,3)  Start and end samples of each frame. tt(*,3)=1 at the start of each talk spurt
%
% Plots a graph if no outputs are specified showing lag candidates and selected path
%

% Bugs/Suggestions:
%   (1) Include backward DP pass and output the true cost for each candidate.
%   (2) Add an extra state to distinguish between voiceless and silent
%   (3) N-best DP to allow longer term penalties (e.g. for frequent pitch doubling/halving)

% The algorithm is taken from [1] with the following differences:
%
%      (a)  the factor AFACT which in the Talkin algorithm corresponds roughly
%           to the absolute level of harmonic noise in the correlation window. This value
%           is here calculated as the maximum of three figures:
%                   (i) an absolute floor set by p.absnoise
%                  (ii) a multiple of the peak signal set by p.signoise
%                 (iii) a multiple of the noise floor set by p.relnoise
%      (b) The LPC used in calculating the Itakura distance uses a Hamming window rather than
%          a Hanning window.
%
% A C implementation of this algorithm by Derek Lin and David Talkin is included as  "get_f0.c"
% in the esps.zip package available from http://www.speech.kth.se/esps/esps.zip under the BSD
% license.
%
% Refs:
%      [1]   D. Talkin, "A Robust Algorithm for Pitch Tracking (RAPT)"
%            in "Speech Coding & Synthesis", W B Kleijn, K K Paliwal eds,
%            Elsevier ISBN 0444821694, 1995

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

s=s(:); % force s to be a column
if nargin<4
    q=[];
    if nargin<3
        mode=' ';
    end
end
doback=0;   % don't do backwards DP for now

% set default parameters

p0.f0min=50;           % Min F0 (Hz)
p0.f0max=500;          % Max F0 (Hz)
p0.tframe=0.01;        % frame size (s)
p0.tlpw=0.005;         % low pass filter window size (s)
p0.tcorw=0.0075;       % correlation window size (s)
p0.candtr=0.3;         % minimum peak in NCCF
p0.lagwt=0.3;          % linear lag taper factor
p0.freqwt=0.02;        % cost factor for F0 change
p0.vtranc=0.005;       % fixed voice-state transition cost
p0.vtrac=0.5;          % delta amplitude modulated transition cost
p0.vtrsc=0.5;          % delta spectrum modulated transition cost
p0.vobias=0.0;         % bias to encourage voiced hypotheses
p0.doublec=0.35;       % cost of exact doubling or halving
p0.absnoise=0;         % absolute rms noise level
p0.relnoise=2;         % rms noise level relative to noise floor
p0.signoise=0.001;     % ratio of peak signal rms to noise floor (0.001 = 60dB)
p0.ncands=20;          % max hypotheses at each frame
p0.trms=0.03;                      % window length for rms measurement
p0.dtrms=0.02;                     % window spacing for rms measurement
p0.preemph=-7000;                  % s-plane position of preemphasis zero
p0.nfullag=7;                      % number of full lags to try (must be odd)
p=v_paramsetch(p0,q);

% redefine commonly used parameters

candtr=p.candtr;          % minimum peak in NCCF                      [0.3]
vtranc=p.vtranc;          % fixed voice-state transition cost         [0.005]
vtrac=p.vtrac;            % delta amplitude modulated transition cost [0.5]
vtrsc=p.vtrsc;            % delta spectrum modulated transition cost  [0.5]
vobias=p.vobias;          % bias to encourage voiced hypotheses       [0.0]
doublec=p.doublec;        % cost of exact doubling or halving         [0.35]
ncands=p.ncands;          % max hypotheses at each frame              [20]
nfullag=p.nfullag;        % number of full lags to try (must be odd)  [7]

% derived parameters (mostly dependent on sample rate fs)

krms=round(p.trms*fs);            % window length for rms measurement
kdrms=round(p.dtrms*fs);          % window spacing for rms measurement
% rmswin=hanning(krms).^2;
rmswin=(0.5-0.5*cos(2*pi*(1:krms)'/(krms+1))).^2; % squared Hanning window
kdsmp=round(0.25*fs/p.f0max);
hlpw=round(p.tlpw*fs/2);          % force window to be an odd length
blp=sinc((-hlpw:hlpw)/kdsmp).*hamming(2*hlpw+1).';
fsd=fs/kdsmp;
kframed=round(fsd*p.tframe);      % downsampled frame length
kframe=kframed*kdsmp;           % frame increment at full rate
rmsix=(1:krms)+floor((kdrms-kframe)/2); % rms index according to Talkin; better=(1:krms)+floor((kdrms-krms+1)/2)
minlag=ceil(fsd/p.f0max);
maxlag=round(fsd/p.f0min);        % use round() only because that is what Talkin does
kcorwd=round(fsd*p.tcorw);        % downsampled correlation window
kcorw=kcorwd*kdsmp;             % full rate correlation window
spoff=max(hlpw-floor(kdsmp/2),1+kdrms-rmsix(1)-kframe);  % offset for first speech frame at full rate
sfoff=spoff-hlpw+floor(kdsmp/2); % offset for downsampling filter
sfi=1:kcorwd;                   % initial decimated correlation window index array
sfhi=1:kcorw;                   % initial correlation window index array
sfj=1:kcorwd+maxlag;
sfmi=repmat((minlag:maxlag)',1,kcorwd)+repmat(sfi,maxlag-minlag+1,1);
lagoff=(minlag-1)*kdsmp;        % lag offset when converting to high sample rate
beta=p.lagwt*p.f0min/fs;            % bias towards low lags
log2=log(2);
lpcord=2+round(fs/1000);        % lpc order for itakura distance
hnfullag=floor(nfullag/2);
jumprat=exp((doublec+log2)/2);  % lag ratio at which octave jump cost is lowest
ssq=s.^2;
csssq=cumsum(ssq);
sqrt(min(csssq(kcorw+1:end)-csssq(1:end-kcorw))/kcorw);
afact=max([p.absnoise^2,max(ssq)*p.signoise^2,min(csssq(kcorw+1:end)-csssq(1:end-kcorw))*(p.relnoise/kcorw)^2])^2*kcorw^2;

% downsample signal to approx 2 kHz to speed up autocorrelation calculation
% kdsmp is the downsample factor

sf=filter(blp/sum(blp),1,s(sfoff+1:end));
sp=filter([1 exp(p.preemph/fs)],1,s); % preemphasised speech for LPC calculation
sf(1:length(blp)-1)=[];         % remove startup transient
sf=sf(1:kdsmp:end);             % downsample to =~2kHz
nsf=length(sf);                 % length of downsampled speech
ns=length(s);                   % length of full rate speech

% Calculate the frame limit to ensure we don't run off the end of the speech or decimated speech:
%   (a) For decimated autocorrelation when calculating sff():  (nframe-1)*kframed+kcorwd+maxlag <= nsf
%   (b) For full rate autocorrelation when calculating sfh():  max(fho)+kcorw+maxlag*kdsamp+hnfllag <= ns
%   (c) For rms ratio window when calculating rr            :  max(fho)+rmsix(end) <= ns
% where max(fho) = (nframe-1)*kframe + spoff

nframe=floor(1+min((nsf-kcorwd-maxlag)/kframed,(ns-spoff-max(kcorw-maxlag*kdsmp-hnfullag,rmsix(end)))/kframe));

% now search for autocorrelation peaks in the downsampled signal

cost=zeros(nframe,ncands);      % cumulative cost
prev=zeros(nframe,ncands);      % traceback pointer
mcands=zeros(nframe,1);         % number of actual candidates excluding voiceless
lagval=repmat(NaN,nframe,ncands-1);    % lag of each voiced candidate
tv=zeros(nframe,3);             % diagnostics: 1=voiceless cost, 2=min voiced cost, 3:cumulative voiceless-min voiced
if doback
    costms=cell(nframe,1);
end

% Main processing loop for each 10 ms frame

for iframe=1:nframe       % loop for each frame (~10 ms)

    % Find peaks in the normalized autocorrelation of subsampled (2Khz) speech
    % only keep peaks that are > 30% of highest peak

    fho=(iframe-1)*kframe+spoff;
    sff=sf((iframe-1)*kframed+sfj);
    sffdc=mean(sff(sfi));       % mean of initial correlation window length
    sff=sff-sffdc;              % subtract off the mean
    nccfd=normxcor(sff(1:kcorwd),sff(minlag+1:end));
    [ipkd,vpkd]=v_findpeaks(nccfd,'q');

    % Debugging: execute the line below to plot the autocorrelation peaks.
    % v_findpeaks(nccfd,'q'); xlabel(sprintf('Lag = (x+%d)*%g ms',minlag-1,1000*kdsmp/fs)); ylabel('Normalized Cross Correlation'); title (sprintf('Frame %d/%d',iframe,nframe));

    vipkd=[vpkd ipkd];
    vipkd(vpkd<max(vpkd)*candtr,:)=[];          % eliminate peaks that are small
    if size(vipkd,1)
        if size(vipkd,1)>ncands-1
            vipkd=sortrows(vipkd);
            vipkd(1:size(vipkd,1)-ncands+1,:)=[];   % eliminate lowest to leave only ncands-1
        end
        lagcan=round(vipkd(:,2)*kdsmp+lagoff);        % convert the lag candidate values to the full sample rate
        nlcan=length(lagcan);
    else
        nlcan=0;
    end

    % If there are any candidate lag values (nlcan>0) then refine their accuracy at the full sample rate

    if nlcan
        laglist=reshape(repmat(lagcan(:)',nfullag,1)+repmat((-hnfullag:hnfullag)',1,nlcan),nfullag*nlcan,1);
        sfh=s(fho+(1:kcorw+max(lagcan)+hnfullag));
        sfhdc=mean(sfh(sfhi));
        sfh=sfh-sfhdc;
        e0=sum(sfh(sfhi).^2);                     % energy of initial correlation window (only needed to store in tv(:,6)
        lagl2=repmat(lagcan(:)',nfullag+kcorw-1,1)+repmat((1-hnfullag:hnfullag+kcorw)',1,nlcan);
        nccf=normxcor(sfh(1:kcorw),sfh(lagl2),afact);

        [maxcc,maxcci]=max(nccf,[],1);
        vipk=[maxcc(:) lagcan(:)+maxcci(:)-hnfullag-1];
        vipk=vipk(:,[1 2 2]);
        maxccj=maxcci(:)'+nfullag*(0:nlcan-1);    % vector index into nccf array
        msk=mod(maxcci,nfullag-1)~=1 & 2*nccf(maxccj)-nccf(mod(maxccj-2,nfullag*nlcan)+1)-nccf(mod(maxccj,nfullag*nlcan)+1)>0;  % don't do quadratic interpolation for the end ones
        if any(msk)
            maxccj=maxccj(msk);
            vipk(msk,3)=vipk(msk,3)+(nccf(maxccj+1)-nccf(maxccj-1))'./(2*(2*nccf(maxccj)-nccf(maxccj-1)-nccf(maxccj+1)))';
        end
        vipk(maxcc<max(maxcc)*candtr,:)=[];          % eliminate peaks that are small
        if size(vipk,1)>ncands-1
            vipk=sortrows(vipk);
            vipk(1:size(vipk,1)-ncands+1,:)=[];   % eliminate lowest to leave only ncands-1
        end

        % vipk(:,1) has NCCF value, vipk(:,2) has integer peak position, vipk(:,3) has refined peak position

        mc=size(vipk,1);
    else
        mc=0;
    end

    % We now have mc lag candidates at the full sample rate

    mc1=mc+1;               % total number of candidates including "unvoiced" possibility
    mcands(iframe)=mc;      % save number of lag candidates (needed for pitch consistency cost calculation)
    if mc
        lagval(iframe,1:mc)=vipk(:,3)';
        cost(iframe,1)=vobias+max(vipk(:,1));   % voiceless cost
        cost(iframe,2:mc1)=1-vipk(:,1)'.*(1-beta*vipk(:,3)');   % local voiced costs
        tv(iframe,2)=min(cost(iframe,2:mc1));
    else
        cost(iframe,1)=vobias;          % if no lag candidates (mc=0), then the voiceless case is the only possibility
    end
    tv(iframe,1)=cost(iframe,1);
    if iframe>1                         % if it is not the first frame, then calculate pitch consistency and v/uv transition costs
        mcp=mcands(iframe-1);
        costm=zeros(mcp+1,mc1);         % cost matrix: rows and cols correspond to candidates in previous and current frames (incl voiceless)

        % if both frames have at least one lag candidate, then calculate a pitch consistency cost

        if mc*mcp
            lrat=abs(log(repmat(lagval(iframe,1:mc),mcp,1)./repmat(lagval(iframe-1,1:mcp)',1,mc)));
            costm(2:end,2:end)=p.freqwt*min(lrat,doublec+abs(lrat-log2));  % allow pitch doubling/halving
        end

        % if either frame has a lag candidate, then calculate the cost of voiced/voiceless transition and vice versa

        if mc+mcp
            rr=sqrt((rmswin'*s(fho+rmsix).^2)/(rmswin'*s(fho+rmsix-kdrms).^2)); % amplitude "gradient"
            ss=0.2/(v_distitar(v_lpcauto(sp(fho+rmsix),lpcord),v_lpcauto(sp(fho+rmsix-kdrms),lpcord),'e')-0.8);   % Spectral stationarity: note: Talkin uses Hanning instead of Hamming windows for LPC
            costm(1,2:end)= vtranc+vtrsc*ss+vtrac/rr;   % voiceless -> voiced cost
            costm(2:end,1)= vtranc+vtrsc*ss+vtrac*rr;
            tv(iframe,4:5)=[costm(1,mc1) costm(mcp+1,1)];
        end
        costm=costm+repmat(cost(iframe-1,1:mcp+1)',1,mc1);  % add in cumulative costs
        [costi,previ]=min(costm,[],1);
        cost(iframe,1:mc1)=cost(iframe,1:mc1)+costi;
        prev(iframe,1:mc1)=previ;
    else                            % first ever frame
        costm=zeros(1,mc1); % create a cost matrix in case doing a backward recursion
    end
    if mc
        tv(iframe,3)=cost(iframe,1)-min(cost(iframe,2:mc1));
        tv(iframe,6)=5*log10(e0*e0/afact);
    end
    if doback
        costms{iframe}=costm; % need to add repmatted cost into this
    end
end

% now do traceback

best=zeros(nframe,1);
[cbest,best(nframe)]=min(cost(nframe,1:mcands(nframe)+1));
for i=nframe:-1:2
    best(i-1)=prev(i,best(i));
end
vix=find(best>1);
fx=repmat(NaN,nframe,1);                        % unvoiced frames will be NaN
fx(vix)=fs*lagval(vix+nframe*(best(vix)-2)).^(-1); % leave as NaN if unvoiced
tt=zeros(nframe,3);
tt(:,1)=(1:nframe)'*kframe+spoff;       % find frame times
tt(:,2)=tt(:,1)+kframe-1;
jratm=(jumprat+1/jumprat)/2;
tt(2:end,3)=abs(fx(2:end)./fx(1:end-1)-jratm)>jumprat-jratm;    % new spurt if frequency ratio is outside (1/jumprat,jumprat)
tt(1,3)=1;           % first frame always starts a spurt
tt(1+find(isnan(fx(1:end-1))),3)=1; % NaN always forces a new spurt

% plot results if there are no output arguments of if the 'g' mode option is specified

if ~nargout || any(mode=='g')
    tf=spoff+(0:nframe-1)'*kframe;      % one sample before start of each frame
    blag=repmat(NaN,nframe,1);                        % unvoiced frames will be NaN
    blag(vix)=lagval(vix+nframe*(best(vix)-2)); % leave as NaN if unvoiced
    ts=(1:ns)/fs;                       % time scale for speech samples
    tsa=[1:tf(1) tf(end)+kframe+1:ns];  % indexes for unprocessed speech [-1 term is an error methinks]
    sup=repmat(NaN,ns,1);               % unprocessed speech - plot in black
    sup(tsa)=s(tsa);
    sv=reshape(s(tf(1)+1:tf(end)+kframe),kframe,nframe);               % processed speech
    su=sv;
    su(:,best>1)=NaN;                   % delete all voiced samples
    sv(:,best==1)=NaN;                  % delete all unvoiced samples
    tsuv=(tf(1)+1:tf(end)+kframe)/fs;
    su=su(:);
    sv=sv(:);
    ax=zeros(2,1);
    ax(1)=subplot(211);
    plot(ts,sup,'-k',tsuv,su,'r-',tsuv,sv,'b-');
    title('Speech');
    ax(2)=subplot(212);
    plot((tf+(kframe+1)/2)/fs,lagval*1000/fs,'xr',(tf+(kframe+1)/2)/fs,blag*1000/fs,'-b')
    xlabel('Time (s)');
    ylabel('Period (ms)');
    title('Lag Candidates');
    linkaxes(ax,'x');
end
if ~any(mode=='u')
    tt(isnan(fx),:)=[];    % remove NaN spurts
    fx(isnan(fx),:)=[];
end

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

function v=normxcor(x,y,d)
% Calculate the normalized cross correlation of column vectors x and y
% we can calculate this in two ways but fft is much faster even for nx small
% We must have nx<=ny and the output length is ny-nx+1
% note that this routine does not do mean subtraction even though this is normally a good idea
% if y is a matrix, we correlate with each column
% d is a constant added onto the normalization factor
% v(j)=x'*yj/sqrt(d + x'*x * yj'*yj) where yj=y(j:j+nx-1) for j=1:ny-nx+1

if nargin<3
    d=0;
end
nx=length(x);
[ny,my]=size(y);
nv=1+ny-nx;
if nx>ny
    error('second argument is shorter than the first');
end

nf=pow2(nextpow2(ny));
w=v_irfft(repmat(conj(v_rfft(x,nf,1)),1,my).*v_rfft(y,nf,1));
s=zeros(ny+1,my);
s(2:end,:)=cumsum(y.^2,1);
v=w(1:nv,:)./sqrt(d+(x'*x).*(s(nx+1:end,:)-s(1:end-nx,:)));

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function int_exactness_gen_hermite ( quad_filename, degree_max, alpha, option )

%*****************************************************************************80
%
%% MAIN is the main program for INT_EXACTNESS_GEN_HERMITE.
%
%  Discussion:
%
%    This program investigates a generalized Gauss-Hermite quadrature rule
%    by using it to integrate monomials over (-oo,+oo), and comparing the
%    approximate result to the known exact value.
%
%    The user specifies:
%    * the "root" name of the R, W and X files that specify the rule;
%    * DEGREE_MAX, the maximum monomial degree to be checked;
%    * ALPHA, the power of X in the weighting function;
%    * OPTION, whether the rule is for |x|^alpha*exp(-x*x)*f(x) or f(x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Investigate the polynomial exactness of a generalized Gauss-Hermite\n' );
  fprintf ( 1, '  quadrature rule by integrating exponentially weighted\n' );
  fprintf ( 1, '  monomials up to a given degree over the (-oo,+oo) interval.\n' );
%
%  Get the quadrature file root name:
%
  if ( 1 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE:\n' );

    quad_filename = input ( '  Enter the "root" name of the quadrature files.' );

  end
%
%  Create the names of:
%    the quadrature X file;
%    the quadrature W file;
%    the quadrature R file;
%
  quad_x_filename = strcat ( quad_filename, '_x.txt' );
  quad_w_filename = strcat ( quad_filename, '_w.txt' );
  quad_r_filename = strcat ( quad_filename, '_r.txt' );
%
%  The second command line argument is the maximum degree.
%
  if ( 2 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE:\n' );

    degree_max = input ( '  Please enter the maximum degree to check.' );

  end
%
%  The third command line argument is ALPHA.
%
  if ( 3 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE:\n' );
    fprintf ( 1, '  ALPHA is the power of |X| in the weighting function.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ALPHA is a real number greater than -1.0.\n' );
    fprintf ( 1, '\n' );

    alpha = input ( '  Please enter ALPHA.' );

  end
%
%  The fourth command line argument is OPTION.
%  0 for the standard rule for integrating |x|^alpha*exp(-x*x)*f(x),
%  1 for a rule for integrating f(x).
%
  if ( 4 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE:\n' );
    fprintf ( 1, '  OPTION chooses the standard or modified rule:\n' );
    fprintf ( 1, '  0: standard rule for integrating |x|^alpha*exp(-x*x)*f(x);\n' );
    fprintf ( 1, '  1: modified rule for integrating                     f(x).\n' );
    fprintf ( 1, '\n' );
    option = input ( '  Please enter OPTION' );

  end
%
%  Summarize the input.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE: User input:\n' );
  fprintf ( 1, '  Quadrature rule X file = "%s".\n', quad_x_filename );
  fprintf ( 1, '  Quadrature rule W file = "%s".\n', quad_w_filename );
  fprintf ( 1, '  Quadrature rule R file = "%s".\n', quad_r_filename );
  fprintf ( 1, '  Maximum degree to check = %d\n', degree_max );
  fprintf ( 1, '  Weighting function exponent ALPHA = %f\n', alpha );
  if ( option == 0 )
    fprintf ( 1, '  OPTION = 0, integrate |x|^alpha*exp(-x*x)*f(x).\n' );
  else
    fprintf ( 1, '  OPTION = 1, integrate                     f(x).\n' );
  end
%
%  Read the X file.
%
  [ dim_num, order ] = r8mat_header_read ( quad_x_filename );

  if ( dim_num ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE - Fatal error!\n' );
    fprintf ( 1, '  The spatial dimension should be 1.\n');
    fprintf ( 1, '  The spatial dimension in the X file is %d\n', dim_num );
    error ( 'INT_EXACTNESS_GEN_HERMITE - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', dim_num );
  fprintf ( 1, '  Number of points  = %d\n', order );

  x = r8mat_data_read ( quad_x_filename, dim_num, order );
%
%  Read the W file.
%
  [ dim_num2, point_num ] = r8mat_header_read ( quad_w_filename );

  if ( dim_num2 ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE - Fatal error!\n' );
    fprintf ( 1, '  The quadrature weight file should have exactly\n');
    fprintf ( 1, '  one value on each line.\n' );
    error ( 'INT_EXACTNESS_GEN_HERMITE - Fatal error!' );
  end

  if ( point_num ~= order )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE - Fatal error!\n' );
    fprintf ( 1, '  The quadrature weight file should have exactly\n' );
    fprintf ( 1, '  the same number of lines as the abscissa file.\n' );
    error ( 'INT_EXACTNESS_GEN_HERMITE - Fatal error!' );
  end

  w = r8mat_data_read ( quad_w_filename, 1, order );
%
%  Read the R file.
%
  [ dim_num2, point_num ] = r8mat_header_read ( quad_r_filename );

  if ( dim_num2 ~= dim_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE - Fatal error!\n' );
    fprintf ( 1, '  The quadrature region file should have the same\n' );
    fprintf ( 1, '  number of values on each line as the abscissa file\n' );
    fprintf ( 1, '  does.\n' );
    error ( 'INT_EXACTNESS_GEN_HERMITE - Fatal error!' );
  end

  if ( point_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE - Fatal error!\n' );
    fprintf ( 1, '  The quadrature region file should have two lines.\n' );
    error ( 'INT_EXACTNESS_GEN_HERMITE - Fatal error!' );
  end

  r = r8mat_data_read ( quad_r_filename, dim_num, 2 );
%
%  Print the input quadrature rule.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The quadrature rule to be tested is\n' );
  fprintf ( 1, '  a generalized Gauss-Hermite rule\n' );
  fprintf ( 1, '  ORDER = %d\n', order );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '\n' );
  if ( option == 0 )
    fprintf ( 1, '  OPTION = 0, standard rule:\n' );
    fprintf ( 1, '    Integral ( -oo < x < +oo ) |x|^alpha exp(-x*x) f(x) dx\n' );
    fprintf ( 1, '    is to be approximated by\n' );
    fprintf ( 1, '    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).\n' );
  else
    fprintf ( 1, '  OPTION = 1, modified rule:\n' );
    fprintf ( 1, '    Integral ( -oo < x < +oo ) f(x) dx\n' );
    fprintf ( 1, '    is to be approximated by\n' );
    fprintf ( 1, '    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weights W:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : order
    fprintf ( 1, '  w(%d) = %24.16f\n', i, w(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Abscissas X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : order
    fprintf ( 1, '  x(%d) = %24.16f\n', i, x(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Region R:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 2
    fprintf ( 1, '  r(%d) = %e\n', i, r(i) );
  end
%
%  Explore the monomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A generalized Gauss-Hermite rule would be able to exactly\n' );
  fprintf ( 1, '  integrate monomials up to and including \n' );
  fprintf ( 1, '  degree = %d\n', 2 * order - 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Error    Degree\n' );
  fprintf ( 1, '\n' );

  for degree = 0 : degree_max

    quad_error = monomial_quadrature_gen_hermite ( degree, alpha, order, option, w, x );

    fprintf ( 1, '  %24.16f   %2d\n', quad_error, degree );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'INT_EXACTNESS_GEN_HERMITE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function value = gen_hermite_integral ( expon, alpha )

%*****************************************************************************80
%
%% GEN_HERMITE_INTEGRAL evaluates a monomial generalized Hermite integral.
%
%  Discussion:
%
%    H(n,alpha) = Integral ( -oo < x < +oo ) x^n |x|^alpha exp(-x^2) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, int EXPON, the exponent of the monomial.
%
%    Input, real ALPHA, the exponent of |X| in the integral.
%    -1.0 < ALPHA.
%
%    Output, real VALUE, the value of the integral.
%
  if ( mod ( expon, 2 ) == 1 )

    value = 0.0;

  else

    a = alpha + expon;

    if ( a <= -1.0 )

      value = - r8_huge ( );

    else

      value = r8_gamma ( ( a + 1.0 ) / 2.0 );

    end

  end

  return
end
function quad_error = monomial_quadrature_gen_hermite ( expon, alpha, order, ...
  option, w, x )

%*****************************************************************************80
%
%% MONOMIAL_QUADRATURE_GEN_HERMITE applies a quadrature rule to a monomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer EXPON, the exponent.
%
%    Input, real ALPHA, the exponent of X in the weight factor.
%
%    Input, intege ORDER, the number of points in the rule.
%
%    Input, integer OPTION, indicates standard or modified rule.
%    0, standard generalized Gauss-Hermite rule for
%       integrand |x|^alpha*exp(-x*x)*f(x).
%    1, modified generalized Gauss-Laguerre rule for
%       integrand                     f(x).
%
%    Input, real W(ORDER), the quadrature weights.
%
%    Input, real X(ORDER), the quadrature points.
%
%    Output, real QUAD_ERROR, the quadrature error.
%

%
%  Get the exact value of the integral of the monomial.
%
   exact = gen_hermite_integral ( expon, alpha );
%
%  Evaluate the unweighted monomial at the quadrature points.
%
  if ( option == 0 )
    value(1:order) = x(1:order).^expon;
  else
    value(1:order) = ( abs ( x(1:order) ) ).^alpha ...
    .* exp ( - x(1:order).^2 ) .* x(1:order).^expon;
  end
%
%  Compute the weighted sum.
%
  quad = w(1:order) * value(1:order)';
%
%  Error:
%
  if ( exact == 0.0 )
    quad_error = abs ( quad );
  else
    quad_error = abs ( ( quad - exact ) / exact );
  end

  return
end
function value = r8_gamma ( x )

%*****************************************************************************80
%
%% R8_GAMMA evaluates Gamma(X) for a real argument.
%
%  Discussion:
%
%    This routine calculates the gamma function for a real argument X.
%
%    Computation is based on an algorithm outlined in reference 1.
%    The program uses rational functions that approximate the gamma
%    function to at least 20 significant decimal digits.  Coefficients
%    for the approximation over the interval (1,2) are unpublished.
%    Those for the approximation for 12 <= X are from reference 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by William Cody, Laura Stoltz.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    William Cody,
%    An Overview of Software Development for Special Functions,
%    in Numerical Analysis Dundee, 1975,
%    edited by GA Watson,
%    Lecture Notes in Mathematics 506,
%    Springer, 1976.
%
%    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
%    Charles Mesztenyi, John Rice, Henry Thatcher,
%    Christoph Witzgall,
%    Computer Approximations,
%    Wiley, 1968,
%    LC: QA297.C64.
%
%  Parameters:
%
%    Input, real X, the argument of the function.
%
%    Output, real VALUE, the value of the function.
%

%
%  Coefficients for minimax approximation over (12, INF).
%
  c = [ ...
   -1.910444077728E-03, ...
    8.4171387781295E-04, ...
   -5.952379913043012E-04, ...
    7.93650793500350248E-04, ...
   -2.777777777777681622553E-03, ...
    8.333333333333333331554247E-02, ...
    5.7083835261E-03 ];
%
%  Mathematical constants
%
  one = 1.0;
  half = 0.5;
  twelve = 12.0;
  two = 2.0;
  zero = 0.0;
  sqrtpi = 0.9189385332046727417803297;
%
%  Machine dependent parameters
%
  xbig = 171.624E+00;
  xminin = 2.23E-308;
  eps = 2.22E-16;
  xinf = 1.79E+308;
%
%  Numerator and denominator coefficients for rational minimax
%  approximation over (1,2).
%
  p = [ ...
   -1.71618513886549492533811E+00, ...
    2.47656508055759199108314E+01, ...
   -3.79804256470945635097577E+02, ...
    6.29331155312818442661052E+02, ...
    8.66966202790413211295064E+02, ...
   -3.14512729688483675254357E+04, ...
   -3.61444134186911729807069E+04, ...
    6.64561438202405440627855E+04 ];

  q = [ ...
   -3.08402300119738975254353E+01, ...
    3.15350626979604161529144E+02, ...
   -1.01515636749021914166146E+03, ...
   -3.10777167157231109440444E+03, ...
    2.25381184209801510330112E+04, ...
    4.75584627752788110767815E+03, ...
   -1.34659959864969306392456E+05, ...
   -1.15132259675553483497211E+05 ];

  parity = 0;
  fact = one;
  n = 0;
  y = x;
%
%  Argument is negative.
%
  if ( y <= zero )

    y = - x;
    y1 = floor ( y );
    res = y - y1;

    if ( res ~= zero )

      if ( y1 ~= floor ( y1 * half ) * two )
        parity = 1;
      end

      fact = - pi / sin ( pi * res );
      y = y + one;

    else

      res = xinf;
      value = res;
      return

    end

  end
%
%  Argument is positive.
%
  if ( y < eps )
%
%  Argument < EPS.
%
    if ( xminin <= y )
      res = one / y;
    else
      res = xinf;
      value = res;
      return
    end

  elseif ( y < twelve )

    y1 = y;
%
%  0.0 < argument < 1.0.
%
    if ( y < one )

      z = y;
      y = y + one;
%
%  1.0 < argument < 12.0.
%  Reduce argument if necessary.
%
    else

      n = floor ( y ) - 1;
      y = y - n;
      z = y - one;

    end
%
%  Evaluate approximation for 1.0 < argument < 2.0.
%
    xnum = zero;
    xden = one;
    for i = 1 : 8
      xnum = ( xnum + p(i) ) * z;
      xden = xden * z + q(i);
    end

    res = xnum / xden + one;
%
%  Adjust result for case  0.0 < argument < 1.0.
%
    if ( y1 < y )

      res = res / y1;
%
%  Adjust result for case 2.0 < argument < 12.0.
%
    elseif ( y < y1 )

      for i = 1 : n
        res = res * y;
        y = y + one;
      end

    end

  else
%
%  Evaluate for 12.0 <= argument.
%
    if ( y <= xbig )

      ysq = y * y;
      sum = c(7);
      for i = 1 : 6
        sum = sum / ysq + c(i);
      end
      sum = sum / y - y + sqrtpi;
      sum = sum + ( y - half ) * log ( y );
      res = exp ( sum );

    else

      res = xinf;
      value = res;
      return

    end

  end
%
%  Final adjustments and return.
%
  if ( parity )
    res = - res;
  end

  if ( fact ~= one )
    res = fact / res;
  end

  value = res;

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [featvec] = MR8fast(im)
%computes MR8 filterbank using recursive Gaussian filters
%input: intensity image
%output: MR8 feature vector

persistent MR8filterNorm;
obtainFilterNorm = 0;

% problematic to do L1 normalization for recursive filters,
% hence solve it this way...
if (isempty(MR8filterNorm))
    MR8filterNorm = 1;
    obtainFilterNorm = 1;
    a=zeros(256,256);
    a(128,128)=1;
end;

in = double(im) - mean(mean(im));
in = in ./ sqrt(mean(mean(in .^ 2)));

ims = cell(1, 8);
i=1;
n=1;

sfac = 0.25;% 1.0;
mulfac = 2.0;

s1 = 3*sfac; s2 = 1*sfac;
for j=0:2,
    for k=0:5,
        phi = (k/6.0)*180.0;

        if (obtainFilterNorm)
            % this should be done only once....
            im1 = s2 .* anigauss(a, s1, s2, phi, 0, 1);
            im2 = (s2*s2) .* anigauss(a, s1, s2, phi, 0, 2);
            n1 = 1.0/sum(sum(abs(im1)));
            n2 = 1.0/sum(sum(abs(im2)));
            MR8filterNorm = [MR8filterNorm, n1, n2];
        else
            n1 = MR8filterNorm(n);
            n = n+1;
            n2 = MR8filterNorm(n);
            n = n+1;
        end;

        im1 = n1 .* anigauss(in, s1, s2, phi, 0, 1);
        im2 = n2 .* anigauss(in, s1, s2, phi, 0, 2);

        % take max of abs response for first order derivative
        % Varma&Zisserman also take abs max of second order...
        im1 = abs(im1);
        %im2 = abs(im2);
        if (k==0)
            maxim1 = im1;
            maxim2 = im2;
        else
            maxim1 = max(maxim1, im1);
            maxim2 = max(maxim2, im2);
        end
    end

    ims{i} = maxim1; i=i+1;
    ims{i} = maxim2; i=i+1;

    % next octave
    s1 = s1*mulfac; s2 = s2*mulfac;
end

sigma = 10.0*sfac;

if (obtainFilterNorm)
    im1 = anigauss(a, sigma, sigma, 0.0, 2, 0);
    im2 = anigauss(a, sigma, sigma, 0.0, 0, 2);
    im1 = (s2*s2) .* (im1+im2);
    im2 = anigauss(a, sigma, sigma);
    n1 = 1.0/sum(sum(abs(im1)));
    n2 = 1.0/sum(sum(abs(im2))); % this one normally should be positive
    MR8filterNorm = [MR8filterNorm, n1, n2];
else
    n1 = MR8filterNorm(n);
    n = n+1;
    n2 = MR8filterNorm(n);
end;

im1 = anigauss(in, sigma, sigma, 0.0, 2, 0);
im2 = anigauss(in, sigma, sigma, 0.0, 0, 2);
ims{i} = n1 .* (im1+im2);
i=i+1;
ims{i} = n2 .* anigauss(in, sigma, sigma);

% just throw away 25 pixel border...(half support of sigma=10 filter)
if 0
[R,C] = size(ims{1});
for j=1:i,
    ims{j} = ims{j}(26:R-25,26:C-25);
end
end

featvec = [ims{8}(:) ims{7}(:) ims{1}(:) ims{3}(:) ims{5}(:) ims{2}(:) ims{4}(:) ims{6}(:)]';

author: Cloud-CV

repo: object-proposals


text:
% ------------------------------------------------------------------------
%      Analyze non-linear ROIs that were clicked as 'object points'
% ------------------------------------------------------------------------

%% ENTER PARAMETERS AND FILE LOCATION

% file location of object points
save_folder = 'C:\Drive\Histology\brainX\processed';

% directory of reference atlas files
annotation_volume_location = 'C:\Drive\Histology\for tutorial\annotation_volume_10um_by_index.npy';
structure_tree_location = 'C:\Drive\Histology\for tutorial\structure_tree_safe_2017.csv';

% name of the saved object points
object_save_name_suffix = '';

% either set to 'all' or a list of indices from the clicked objects in this file, e.g. [2,3]
objects_to_analyze = 'all';

% plane used to view when points were clicked ('coronal' -- most common, 'sagittal', 'transverse')
plane = 'coronal';

% brain figure black or white
black_brain = true;

%% LOAD THE REFERENCE ANNOTATIONS AND PROBE POINTS

% load the reference brain annotations
if ~exist('av','var') || ~exist('st','var')
    disp('loading reference atlas...')
    av = readNPY(annotation_volume_location);
    st = loadStructureTree(structure_tree_location);
end

% load object points
objectPoints = load(fullfile(save_folder, ['probe_points' object_save_name_suffix]));

% determine which objects to analyze
if strcmp(objects_to_analyze,'all')
    objects = 1:size(objectPoints.pointList.pointList,1);
else
    objects = objects_to_analyze;
end

%% BRING UP THE RELEVANT DATA FOR EACH PROBE POINTS, FOR FURTHER ANALYSIS

% initialize cell array containing info on each clicked point
if length(objects) > 1
    roi_annotation = cell(length(objects),1);
    roi_location = cell(length(objects),1);
end

% generate needed values
bregma = allenCCFbregma(); % bregma position in reference data space
atlas_resolution = 0.010; % mm

% plot brain grid
ProbeColors = [1 1 1; 1 .75 0;  .3 1 1; .4 .6 .2; 1 .35 .65; .7 .7 1; .65 .4 .25; .7 .95 .3; .7 0 0; .6 0 .7; 1 .6 0];
% order of colors: {'white','gold','turquoise','fern','bubble gum','overcast sky','rawhide', 'green apple','purple','orange','red'};
fwireframe = plotBrainGrid([], [], [], black_brain); hold on;
fwireframe.InvertHardcopy = 'off';

for object_num = objects

    selected_object = objects(object_num);

    % get the object points for the currently analyzed object
    if strcmp(plane,'coronal')
        curr_objectPoints = objectPoints.pointList.pointList{selected_object,1}(:, [3 2 1]);
    elseif strcmp(plane,'sagittal')
        curr_objectPoints = objectPoints.pointList.pointList{selected_object,1}(:, [1 2 3]);
    elseif strcmp(plane,'transverse')
        curr_objectPoints = objectPoints.pointList.pointList{selected_object,1}(:, [1 3 2]);
    end

    % plot points on the wire frame brain
    figure(fwireframe); hold on
    hp = plot3(curr_objectPoints(:,1), curr_objectPoints(:,3), curr_objectPoints(:,2), '.','linewidth',2, 'color',[ProbeColors(object_num,:) .2],'markers',10);

    % use the point's position in the atlas to get the AP, DV, and ML coordinates
    ap = -(curr_objectPoints(:,1)-bregma(1))*atlas_resolution;
    dv = (curr_objectPoints(:,2)-bregma(2))*atlas_resolution;
    ml = (curr_objectPoints(:,3)-bregma(3))*atlas_resolution;

    roi_location_curr = [ap dv ml];

    % initialize array of region annotations
    roi_annotation_curr = cell(size(curr_objectPoints,1),3);

    % loop through every point to get ROI locations and region annotations
    for point = 1:size(curr_objectPoints,1)

        % find the annotation, name, and acronym of the current ROI pixel
        ann = av(curr_objectPoints(point,1),curr_objectPoints(point,2),curr_objectPoints(point,3));
        name = st.safe_name{ann};
        acr = st.acronym{ann};

        roi_annotation_curr{point,1} = ann;
        roi_annotation_curr{point,2} = name;
        roi_annotation_curr{point,3} = acr;

    end

    % save results in cell array
    if length(objects) > 1
        roi_annotation{object_num} = roi_annotation_curr;
        roi_location{object_num} = roi_location_curr;
    else
        roi_annotation = roi_annotation_curr;
        roi_location = roi_location_curr;
    end

    % display results in a table
    disp(['Clicked points for object ' num2str(selected_object)])
    roi_table = table(roi_annotation_curr(:,2),roi_annotation_curr(:,3), ...
                        roi_location_curr(:,1),roi_location_curr(:,2),roi_location_curr(:,3), roi_annotation_curr(:,1), ...
         'VariableNames', {'name', 'acronym', 'AP_location', 'DV_location', 'ML_location', 'avIndex'});
     disp(roi_table)

end

% now, use roi_location and roi_annotation for your further analyses

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function write_html_dir(obj, output_dir)
% Write HTML coverage report for m-file collection
%
% write_html_dir(obj, output_dir)
%
% Inputs:
%   obj                 MOcovMFileCollection instance
%   output_dir          HTML output directory
%
% Notes:
%   - this function writes a file index.html, as well as node*.html files
%     for each individual MOcovMFile.

    if ~mocov_util_isfolder(output_dir)
        mkdir(output_dir);
    end

    monitor=obj.monitor;

    index_rel_fn='index.html';
    index_fn=fullfile(output_dir, index_rel_fn);

    mfiles=obj.mfiles;
    n=numel(mfiles);
    mfile_node_fns=cell(n,1);

    % write single HTML file for each m-file
    notify(monitor,sprintf('Writing html files in %s', output_dir));

    for k=1:n
        node_rel_fn=sprintf('node%d.html',k);
        mfile_node_fns{k}=node_rel_fn;

        mfile=mfiles{k};
        node_fn=fullfile(output_dir, node_rel_fn);
        write_html(mfile, node_fn, index_rel_fn);

        msg=sprintf('Written to %s', node_fn);
        notify(monitor,'.',msg);
    end

    % write index HTML file
    write_index_html(index_fn, mfiles, mfile_node_fns);
    msg=sprintf('Index written to %s\n', index_fn);
    notify(monitor,msg);

function write_index_html(output_fn, mfiles, mfile_node_fns)
    % build an index html file

    fid=fopen(output_fn,'w');
    cleaner=onCleanup(@()fclose(fid));

    fprintf(fid,['<!DOCTYPE html>\n'...
                        '<title>Coverage</title>'...
                        '
\n']);
    fprintf(fid,['
Lines	Executable	'...
                    'Executed	Coverage	'...
                    'File
\n']);

    stats=get_stats(mfiles);

    [unused,i]=sort(stats(:,4));

    n_files=numel(mfiles);

    % pattern with placeholders for:
    % - number of lines
    % - number of executable lines
    % - number of executed lines
    % - string describing the row
    row_pat=['
%d	'...
                        '%d	'...
                        '%d	'...
                        '%.1f%%	'...
                        '%s
\n'];

    % add summary line at top
    overall_stats=sum(stats,1);
    overall_stats(4)=100*overall_stats(3)/overall_stats(2);
    summary_str=sprintf('%d files',n_files);

    fprintf(fid,row_pat,overall_stats,summary_str);

    % add row for each file
    for k=1:numel(mfiles)
        mfile=mfiles{i(k)};
        stat=stats(i(k),:);
        node_fn=mfile_node_fns{i(k)};

        mfile_name=get_filename(mfile);
        mfile_anchor=sprintf('%s',node_fn,mfile_name);

        fprintf(fid,row_pat,stat,mfile_anchor);
    end

    % add summary line at bottom
    fprintf(fid,row_pat,overall_stats,summary_str);

    fprintf(fid,'
');

function stats=get_stats(mfiles)
    n=numel(mfiles);
    stats=zeros(n,4);
    for k=1:n
        mfile=mfiles{k};

        executable=get_lines_executable(mfile);
        executed=get_lines_executed(mfile);

        executed(~executable)=false;

        if sum(executable)==0
            coverage=100;
        else
            coverage=100*sum(executed)/sum(executable);
        end

        stat=[numel(executable),...
                     sum(executable),...
                     sum(executed),...
                     coverage];
        stats(k,:)=stat;
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function vectorfield(F,xval,yval)
% vectorfield([F1,F2],a:dx:b,c:dy:d)
% Plot 2D vectorfield [F1,F2] for
% using x-values from a to b with spacing of dx
%       y-values from c to d with spacing of dy

[xg,yg] = meshgrid(xval,yval);             % values x,y on a grid
F1f = inline(vectorize(F(1)),'x','y');
F2f = inline(vectorize(F(2)),'x','y');
F1g = F1f(xg,yg);   % values of F1 on this grid
F2g = F2f(xg,yg);   % values of F1 on this grid
quiver(xg,yg,F1g,F2g,'k')
axis equal; axis tight

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function test_suite = test_publish_run_scripts()
% smoke tests for publish_run_scripts (using dry-run);
% currently no files are actually generated by this test function
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_publish_run_scripts_dry_run_all_files()
    if cosmo_skip_test_if_no_external('!publish')
        return;
    end

    [srcpat,names_to_publish]=helper_find_files(true);
    helper_run_publish_run_scripts(names_to_publish,srcpat,',''-dry''');

function test_publish_run_scripts_dry_run_single_file()
    if cosmo_skip_test_if_no_external('!publish')
        return;
    end

    [srcpat,names_to_publish]=helper_find_files(false);
    assert(numel(names_to_publish)==1);

    helper_run_publish_run_scripts(names_to_publish,srcpat,',''-dry''');

function test_publish_run_scripts_run_single_file_pass()
    cosmo_notify_test_skipped(['publish functionality tests '...
                                'currently disabled because coverage '...
                                'functionality does not work']);
    return;
    result=helper_publish_run_scripts_single_file('disp(''ok'')');
    assertTrue(result);

function test_publish_run_scripts_run_single_file_error()
    cosmo_notify_test_skipped(['publish functionality tests '...
                                'currently disabled because coverage '...
                                'functionality does not work']);
    return;

    result=helper_publish_run_scripts_single_file('error(''wrong'')');
    assertFalse(result);

function result=helper_publish_run_scripts_single_file(cmd)
    if cosmo_skip_test_if_no_external('!publish')
        return;
    end

    warning_state=cosmo_warning();
    cleaner_warning=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    tmp_dir=tempname();
    cleaner_tmp=onCleanup(@()remove_dir_helper(tmp_dir));
    mkdir(tmp_dir);

    fn=cosmo_make_temp_filename(fullfile(tmp_dir,'run_'),'.m');
    fid=fopen(fn,'w');
    cleaner_fid=onCleanup(@()fclose(fid));
    fprintf(fid,'%s\n',cmd);
    clear cleaner_fid;

    [unused,nm,e]=fileparts(fn);
    args=sprintf(',''-o'',''%s''',tmp_dir);

    result=helper_run_publish_run_scripts({[nm e]},fn,args);

    html_fn=fullfile(tmp_dir,[nm '.html']);
    assertTrue(logical(exist(html_fn,'file')));

function [srcpat,names_to_publish]=helper_find_files(entire_directory)
    if cosmo_skip_test_if_no_external('!evalc')
        return;
    end

    cosmomvpa_rootdir=fileparts(fileparts(mfilename('fullpath')));
    example_dir=fullfile(cosmomvpa_rootdir,'examples');

    to_publish_runs=dir(fullfile(example_dir,'run_*.m'));
    to_publish_demos=dir(fullfile(example_dir,'demo_*.m'));

    to_publish_all=cat(1,to_publish_runs,to_publish_demos);

    names_to_publish={to_publish_all.name};

    if entire_directory
        srcpat=example_dir;
    else
        rp=randperm(numel(names_to_publish));
        names_to_publish=names_to_publish(rp(1));

        srcpat=fullfile(example_dir,names_to_publish{1});
    end

function result=helper_run_publish_run_scripts(names_to_publish,...
                                        srcpat,args)
    assertTrue(numel(names_to_publish)>0,'example directory not found');
    orig_warning_state=cosmo_warning();
    warning_state_resetter=onCleanup(@()cosmo_warning(orig_warning_state));
    cosmo_warning('off');

    expr=sprintf('cosmo_publish_run_scripts(''%s''%s);',...
                            srcpat,args);
    [s,result]=run_evalc_helper(expr);

    lines=cosmo_strsplit(s,'\n');
    n_lines=numel(lines);
    files_processed=cell(n_lines,1);

    for k=1:n_lines
        sp=cosmo_strsplit(lines{k},':');
        if numel(sp)==2
            fn_without_whitespace=regexprep(sp{2},'\s*','');
            mfile_fn=regexprep(fn_without_whitespace,'\.html$','.m');

            files_processed{k}=mfile_fn;
        end
    end

    files_processed_msk=~cellfun(@isempty,files_processed);
    files_processed=files_processed(files_processed_msk);

    assert(all(cosmo_match(names_to_publish,files_processed)));

function [s,result]=run_evalc_helper(expr)
    orig_path=path();
    cleaner=onCleanup(@()path(orig_path));

    s=evalc(['result=' expr]);

function remove_dir_helper(tmp_dir)
    if cosmo_wtf('is_octave')
        rmdir_state=confirm_recursive_rmdir();
        state_resetter=onCleanup(@()confirm_recursive_rmdir(rmdir_state));
        confirm_recursive_rmdir(false,'local');
    end
    rmdir(tmp_dir,'s');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [trainX, testX, minimums, ranges] = featNormalize(trainX, testX)

	minimums = min(trainX, [], 1);
	ranges = max(trainX, [], 1) - minimums + eps;
%
	trainX = (trainX - repmat(minimums, size(trainX, 1), 1)) ./ repmat(ranges, size(trainX, 1), 1);

	testX = (testX - repmat(minimums, size(testX, 1), 1)) ./ repmat(ranges, size(testX, 1), 1);

%     trainX = bsxfun(@rdivide,trainX,sqrt(sum(trainX.^2,2)));
% 	testX = bsxfun(@rdivide,testX,sqrt(sum(testX.^2,2)));
end

author: HuiZeng


text:
function xyz_io_test04 ( )

%*****************************************************************************80
%
%% XYZ_IO_TEST04 tests XYZL_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2009
%
%  Author:
%
%    John Burkardt
%
  xyz_filename = 'cube.xyz';
  xyzl_filename = 'cube.xyzl';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'XYZ_IO_TEST04\n' );
  fprintf ( 1, '  XYZ_HEADER_READ  reads the header of an XYZ  file.\n' );
  fprintf ( 1, '  XYZ_DATA_READ    reads the data   of an XYZ  file.\n' );
  fprintf ( 1, '  XYZL_HEADER_READ reads the header of an XYZL file.\n' );
  fprintf ( 1, '  XYZL_DATA_READ   reads the data   of an XYZL file.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Examine XYZ file "%s".\n', xyz_filename );

  point_num = xyz_header_read ( xyz_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points = %d\n', point_num );

  xyz = xyz_data_read ( xyz_filename, point_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Point data:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : point_num
    fprintf ( 1, '  %4d  %10f  %10f  %10f\n', i, xyz(1:3,i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Examine XYZL file "%s".\n', xyzl_filename );

  [ line_num, line_data_num ] = xyzl_header_read ( xyzl_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of lines      = %d\n', line_num );
  fprintf ( 1, '  Number of line items = %d\n', line_data_num );

  [ line_pointer, line_data ] = xyzl_data_read ( xyzl_filename, line_num, ...
    line_data_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Line pointers:\n' );
  fprintf ( 1, '\n' );

  for line = 1 : line_num
    fprintf ( 1, '  %4d  %8d  %8d\n', line, line_pointer(line), ...
      line_pointer(line+1) - 1 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Line data:\n' );
  fprintf ( 1, '\n' );

  for line = 1 : line_num
    fprintf ( 1, '  %4d    ', line );
    for j = line_pointer(line) : line_pointer(line+1) - 1
      fprintf ( 1, '  %8d', line_data(j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function n = square_arbq_size ( degree )

%*****************************************************************************80
%
%% SQUARE_ARBQ_SIZE: size of quadrature rule for the symmetric square.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2014
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired degree of exactness.
%    1 <= DEGREE <= 30.
%
%    Output, integer N, the number of points in the
%    corresponding rule.
%
  n_save = [ ...
      1,   3,   4,   6,   7, ...
     10,  12,  16,  17,  22, ...
     24,  31,  33,  41,  44, ...
     52,  55,  64,  68,  78, ...
     82,  93,  98, 109, 115, ...
    127, 132, 147, 152, 167 ];

  if ( degree < 1 | 30 < degree )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SQUARE_ARBQ_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of DEGREE.\n' );
    error ( 'SQUARE_ARBQ_SIZE - Fatal error!' );
  end

  n = n_save ( degree );

  return
end


text:
function [ y ] = lp_thresh_real( x,threshold,p )
% generalized softthresholding (GST) based on the paper by Zuo et al.
%
% (c) Marc Fischer, Thomas Kuestner
% ---------------------------------------------------------------------

threshold_nonconvex = (2*threshold*(1-p))^(1/(2-p)) + threshold*p*(2*threshold*(1-p))^((p-1)/(2-p));
y = zeros(size(x));
i0 = find(abs(x)>threshold_nonconvex);

if p < 1
    J = 3;
else
    J = 1;
end;
% if length(i0) > 1
    x0 = x(i0);
    x_temp = abs(x0);
    for j = 1:J
        x_temp = abs(x0) - threshold*p*(x_temp.^(p-1));
    end;
    y(i0) = sign(x0).*x_temp;
% end

author: thomaskuestner


text:
function [FigData, iDS, iFig] = gui_figure_data( hFig, varargin )
% GUI_FIGURE_DATA: Retrieve the time series displayed in figure hFig.
%
% USAGE:  FigData = gui_figure_data( hFig )
%         FigData = gui_figure_data( ..., 'SelectedChannels' )
%         FigData = gui_figure_data( ..., 'SelectedTime' )
%
% INPUT:
%    - hFig : handle to a brainstorm figure
%    - 'SelectedChannels' : if some channels are highlighted in the figure, return only them
%    - 'SelectedTime'     : if a time segment is highlighted in the figure, return only this segment
% OUTPUT:
%    - FigData : structure

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

% === INITIALIZATIONS ===
global GlobalData;
% Parse inputs
isSelectedChannels = any(strcmpi(varargin, 'SelectedChannels'));
isSelectedTime     = any(strcmpi(varargin, 'SelectedTime'));
% Initialize save TimeSeries matrix
FigData = struct('F',          [], ...
                 'Time',       [], ...
                 'FigTitle',   get(hFig, 'Name'), ...
                 'FigType',    [], ...
                 'AxesTitle',  [], ...
                 'AxesLegend', [], ...
                 'Events',     repmat(db_template('Event'),0));
FigData.F = {};
FigData.AxesTitle = {};
FigData.AxesLegend = {};

% === GET FIGURE INFORMATION ===
% Get figure description
[hFig, iFig, iDS] = bst_figures('GetFigure', hFig);
Handles = GlobalData.DataSet(iDS).Figure(iFig).Handles;
% Check downsampling factor
if (Handles.DownsampleFactor > 1)
    error(['The signals were downsampled before being displayed and cannot be exported from the figure.' 10 ...
        'For continuous files: import the recordings first.' 10 ...
        'For imported time series: read directly the corresponding .mat files.' 10 ...
        'For scouts: use the proces "Extract > Scout time series".' ]);
end
% Get y-factor and offsets
Factor = Handles(1).DisplayFactor;
if isfield(Handles(1), 'ChannelOffsets')
    ChannelOffsets = Handles(1).ChannelOffsets;
else
    ChannelOffsets = [];
end
FigData.FigType = GlobalData.DataSet(iDS).Figure(iFig).Id.Type;

% === FIND CHANNELS INDICES ===
% Only valid for data time series (real channels)
if strcmpi(GlobalData.DataSet(iDS).Figure(iFig).Id.Type, 'DataTimeSeries')
    iChannel = GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels;
    Channel = GlobalData.DataSet(iDS).Channel(iChannel);
else
    Channel = [];
end

% Find selected indices
if isSelectedTime
    timeSelection = getappdata(hFig, 'GraphSelection');
    if ~isempty(timeSelection) && ~isinf(timeSelection(2))
        selTime = [min(timeSelection), max(timeSelection)];
        % Get time vector
        if ~isempty(GlobalData.DataSet(iDS).Measures.Time)
            TimeVector = bst_memory('GetTimeVector', iDS, [], 'UserTimeWindow');
        % If no time window defined: use the X values from the lines in axes
        elseif ~isempty(Handles(1).hLines)
            iValidLine = find(Handles(1).hLines > 0, 1);
            hLine = Handles(1).hLines(iValidLine);
            TimeVector = get(hLine, 'XData');
        else
            error('You must define a channel file to use this feature.');
        end
        % Get time indices to save
        iTimeBounds = bst_closest(selTime, TimeVector);
        iTime = iTimeBounds(1):iTimeBounds(2);
    else
        iTime = []; % Select all time indices
    end
else
    iTime = [];
end

% Process axes : one cell by axes displayed
for iAxes = 1:length(Handles)
    % === DATA ===
    % Get lines handles
    iDispChannels = find(ishandle(Handles(iAxes).hLines));
    hLines = Handles(iAxes).hLines(iDispChannels);
    % If no data in these axes => next ones
    if isempty(hLines)
        continue
    end

    % Keep only selected channels
    iSel = [];
    if isSelectedChannels
        % Get display width of all lines
        linesWidth = get(hLines, 'LineWidth');
        if iscell(linesWidth)
            linesWidth = [linesWidth{:}];
        end
        % Selected ones have a width > 1
        iSel = find(linesWidth > 1);
    end
    % If no channel selected: use all the channels
    if isempty(iSel)
        iSel = 1:length(hLines);
    end

    % Extract F (Y-values)
    linesY = get(hLines(iSel), 'YData');
    if iscell(linesY)
        FigData.F{iAxes} = cat(1,linesY{:});
    else
        FigData.F{iAxes} = linesY;
    end

    % Remove offsets
    if ~isempty(ChannelOffsets)
        chOffset = ChannelOffsets(iDispChannels(iSel));
        % Channel with negative offsets: ignore them (set them to zero)
        FigData.F{iAxes}(chOffset < 0, :) = 0;
        % Positive offsets: remove them
        iPosOffset = find(chOffset > 0);
        if ~isempty(iPosOffset)
            % Apply offsets
            FigData.F{iAxes}(iPosOffset,:) = FigData.F{iAxes}(iPosOffset,:) - repmat(chOffset, 1, size(FigData.F{iAxes},2));
        end
    end
    % Remove display factor
    FigData.F{iAxes} = FigData.F{iAxes} ./ Factor;
    % Undo the magnetometer / gradiometer scaling (CTF + Neuromag)
    if ~isempty(Channel)
        FigData.F{iAxes} = bst_scale_gradmag( FigData.F{iAxes}, Channel(iDispChannels(iSel)), 'reverse' );
    end

    % Keep only selected time indices
    if isSelectedTime && ~isempty(iTime)
        FigData.F{iAxes} = FigData.F{iAxes}(:,iTime);
    end

    % If some sensors are not displayed: recreate a full matrix
    if (length(iDispChannels) ~= length(Handles(iAxes).hLines)) && (length(iSel) == length(hLines))
        Ffull = zeros(length(Handles(iAxes).hLines), size(FigData.F{iAxes},2));
        Ffull(iDispChannels(iSel),:) = FigData.F{iAxes};
        FigData.F{iAxes} = Ffull;
    end

    % === COMMENTS ===
    % Get title for these axes
    FigData.AxesTitle{iAxes} = get(get(Handles(iAxes).hAxes, 'Title'), 'String');
    % Get legend of these axes (newer matlab create legend if does not exist: avoided with this test)
    if (bst_get('MatlabVersion') >= 903) && isempty(Handles(iAxes).hAxes.Legend)
        hLegend = [];
    else
        hLegend = legend(Handles(iAxes).hAxes);
    end
    if ~isempty(hLegend)
        if (bst_get('MatlabVersion') < 804)
            legendText = get(findobj(hLegend, 'type','text'), 'string');
            FigData.AxesLegend{iAxes} = bst_flip(legendText,1);
        else
            FigData.AxesLegend{iAxes} = get(hLegend, 'String');
        end
    % Legend are not accessible in the figure: need to reconstruct them
    elseif ~isempty(Channel)
        FigData.AxesLegend{iAxes} = {Channel(iDispChannels).Name};
    else
        FigData.AxesLegend{iAxes} = {};
    end

    % === TIME ===
    % Extract Time (X-values)
    if (iAxes == 1)
        FigData.Time = get(hLines(iSel(1)), 'XData');
        % If subselection of the time indices
        if isSelectedTime && ~isempty(iTime)
            FigData.Time = FigData.Time(iTime);
        end
    end
end
% If no data was extracted : exit
if isempty(FigData.F)
    FigData = [];
    return
end

% ===== GET EVENTS =====
hEventsBar = findobj(hFig, '-depth', 1, 'Tag', 'AxesEventsBar');
if (iAxes == 1) && ~isempty(hEventsBar) && ~isempty(GlobalData.DataSet(iDS).Measures.sFile) && ~isempty(GlobalData.DataSet(iDS).Measures.sFile.events)
    % Get events in current time window
    events = panel_record('GetEventsInTimeWindow', hFig);
    % Change the time vector
    FigData.Events = panel_record('ChangeTimeVector', events, 1./GlobalData.DataSet(iDS).Measures.SamplingRate, FigData.Time);
end

author: brainstorm-tools


text:
function [CL, u_U, v_U, w_U, l_t, l_s] = LiftCoeff(gamma, panel, geo, Fu_bar, Fv_bar, Fw_bar)
%Determine the lift coefficient given the vortex strengths
%and the influence coefficients from the vortex lattice method

u_U.c = 1/(4*pi)*Fu_bar*[gamma.c]';  %Eqn 25 in NASA paper
v_U.c = 1/(4*pi)*Fv_bar*[gamma.c]';    %Eqn 23 in NASA paper
w_U.c = 1/(4*pi)*Fw_bar*[gamma.c]';

u_U.alpha = 1/(4*pi)*Fu_bar*[gamma.alpha]';
v_U.alpha = 1/(4*pi)*Fv_bar*[gamma.alpha]';
w_U.alpha = 1/(4*pi)*Fw_bar*[gamma.alpha]';

for i = 1:geo.ns
    for j = 1:geo.nc
        if i == geo.ns  %wingtip
            DeltaGamma(i,j).c = gamma(i,j).c;
            DeltaGamma(i,j).alpha = gamma(i,j).alpha;
        elseif j == 1 %Leading edge
            DeltaGamma(i,j).c = gamma(i,j).c-gamma(i+1,j).c;
            DeltaGamma(i,j).alpha = gamma(i,j).alpha-gamma(i+1,j).alpha;
        else
            DeltaGamma(i,j).c = DeltaGamma(i,j-1).c + gamma(i,j).c - gamma(i+1,j).c;
            DeltaGamma(i,j).alpha = DeltaGamma(i,j-1).alpha + gamma(i,j).alpha - gamma(i+1,j).alpha;
        end
    end
end

l_t.c = 2/geo.S*[DeltaGamma.c]'.*[panel.cc]'.*[v_U.c];  %Eqn 24 in NASA paper
l_t.alpha = 2/geo.S*[DeltaGamma.alpha]'.*[panel.cc]'.*[v_U.alpha];  %Eqn 24 in NASA paper

l_s.c = 2/geo.S*[gamma.c]'*2.*[panel.s]'.*((ones(size(u_U))-[u_U.c])+[v_U.c].*tan([panel.sweep]'))*cos(geo.dih);  %Eqn 28 in NASA paper
l_s.alpha = 2/geo.S*[gamma.alpha]'*2.*[panel.s]'.*((ones(size(u_U))-[u_U.alpha])+[v_U.alpha].*tan([panel.sweep]'))*cos(geo.dih);  %Eqn 28 in NASA paper

% reshape([l_t.c],geo.ns,geo.nc)
% reshape([l_t.alpha],geo.ns,geo.nc)
% reshape([l_s.c],geo.ns,geo.nc)
% reshape(((ones(size(u_U))-[u_U.alpha])+[v_U.alpha].*tan([panel.sweep]')),geo.ns,geo.nc)
% reshape([panel.s],geo.ns,geo.nc)
% reshape([l_s.alpha],geo.ns,geo.nc)
% reshape([u_U.alpha],geo.ns,geo.nc)
% reshape([v_U.alpha],geo.ns,geo.nc)
% reshape([gamma.alpha],geo.ns,geo.nc)
% sum([gamma.c])
% geo.S

CL.c = 2*(sum([l_t.c])+sum([l_s.c]));
CL.alpha = 2*(sum([l_t.alpha])+sum([l_s.alpha]));

author: Sable

repo: mcbench-benchmarks


text:
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% For whatever reason, the code is very slow and does not converge
% DOES NOT WORK WELL
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc; clear;
%% define the geometry
Nx = 50; % number of cells in x direction
Ny = 30; % number of cells in y direction
W = 50; % [m] length of the domain in x direction
H = 30; % [m] length of the domain in y direction
x1=linspace(0,W, Nx);
x2=x1+0.001;
x=zeros(2*Nx,1);
j=0;
for i=1:Nx
    j=j+1;
    x(j)=x1(i);
    j=j+1;
    x(j)=x2(i);
end
y1=linspace(0,H, Ny);
y2=y1+0.01;
y=zeros(2*Ny,1);
j=0;
for i=1:Ny
    j=j+1;
    y(j)=y1(i);
    j=j+1;
    y(j)=y2(i);
end
m = createMesh2D(x, y); % creates a 2D mesh
[X,Y]=ndgrid(m.cellsize.x, m.cellsize.y);
%% define the physical parametrs
p0 = 1e5; % [bar] pressure
pin = 50e5; % [bar] injection pressure at the left boundary
q_in=5e-5; % [m^3/s] water injection
sw0 = 0; % initial water saturation
swin = 1;
mu_oil = 1e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 2e-12; % [m^2] average reservoir permeability
k_frac=100e-12; % [m^2] fracture permeability
phi0 = 0.2; % average porosity
clx=0.05;
cly=0.05;
V_dp=0.1; % Dykstra-Parsons coef.
perm_val= field2d(2*Nx,2*Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, k0);
for i=4:2:2*Nx
    k.value(i,3:end-2)=k_frac;
end
for j=4:2:2*Ny
    k.value(3:end-2,j)=k_frac;
end
k.value(1:2,:)=k_frac;
k.value(end-1:end,:)=k_frac;
k.value(:,1:2)=k_frac;
k.value(:,end-1:end)=k_frac;
%k.value(1:2,:)=k_frac;
phi=createCellVariable(m, phi0);
krw0 = 0.4;
kro0 = 0.8;
nw = 1;
no = 1;
krw = @(sw)(krw0*sw.^nw);
dkrwdsw = @(sw)(krw0*nw*sw.^(nw-1));
kro = @(sw)(kro0*(1-sw).^no);
dkrodsw = @(sw)(-kro0*no*(1-sw).^(no-1));
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
% change the right boandary to constant pressure (Dirichlet)
BCp.right.a(end-5:end)=0; BCp.right.b(end-5:end)=1; BCp.right.c(end-5:end)=p0;
% change the left boundary to constant flow
BCp.left.a(1:5)=-k_frac/mu_water; BCp.left.b(1:5)=0; BCp.left.c(1:5)=q_in;
% change the left boundary to constant saturation (Dirichlet)
BCs.left.a(1:5)=0; BCs.left.b(1:5)=1; BCs.left.c(1:5)=1;
%% define the time step and solver properties
dt = 10; % [s] time step
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)

t = 0;
while (t<t_end)
    error_p = 1e5;
    error_sw = 1e5;
    % estimation loop: sequential
%     for i = 1:2
%         pgrad = gradientTerm(m, p.value);
%         sw_face = upwindMean(m, -pgrad, sw.value); % average value of water saturation
%         labdao = lo.*funceval(kro, sw_face);
%         labdaw = lw.*funceval(krw, sw_face);
%         labda = labdao+labdaw;
%         Mdiffp1 = diffusionTerm(m, -labda);
%         [Mbcp, RHSbcp] = boundaryCondition(m, BCp);
%         p.value = reshapeCell(m, full((Mdiffp1+Mbcp)\RHSbcp));
%         for j = 1:2
%             pgrad = gradientTerm(m, p.value);
%             sw_face = upwindMean(m, -pgrad, sw.value);
%             labdaw = lw.*funceval(krw, sw_face);
%             dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
%             [Mbcsw, RHSbcsw] = boundaryCondition(m, BCs);
%             [Mtranssw2, RHStrans2] = transientTerm(m, phi, dt, sw);
%             Mconvsw2 = convectionUpwindTerm(m, -dlabdawdsw.*pgrad);
%             RHSs = RHSbcsw+RHStrans2+divergenceTerm(m, (labdaw-dlabdawdsw.*sw_face).*pgrad);
%             Ms = Mbcsw+Mconvsw2+Mtranssw2;
%             sw.value = reshapeCell(m, full(Ms\RHSs));
%         end
%     end
    % Implicit loop
    hold all
    while ((error_p>eps_p) || (error_sw>eps_sw))
        % 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:(2*Nx+1)*(2*Ny+1))));
        sw_new = reshapeCell(m,full(x((2*Nx+1)*(2*Ny+1)+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
    t=t+dt;
    p_old = p;
    sw_old = sw;
    figure(1);visualizeCells(sw);shading interp; drawnow;
    sw_tot=sw.value.*X.*Y/(W*H);
    oil_prod=(sum(sum(sw_tot(2:end-1,2:end-1))));
    figure(2);semilogx(t, p_new(2,2), 'o'); drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function info = get_performance_q4s(Xopt,yopt,Sopt,SDP,POP,path)

if iscell(path)
    for i = 1:length(path)
        addpath(genpath(path{i}))
    end
else
    addpath(genpath(path))
end

blk         = SDP.blk;
At          = SDP.At;
b           = SDP.b;
C           = svec(blk,SDP.C);
Xoptmat     = Xopt;
Xopt        = svec(blk,Xopt);
Sopt        = svec(blk,Sopt);
Amap  = @(X) AXfun(blk,At,X);
ATmap = @(y) Atyfun(blk,At,y);

%% standard SDP KKT residuals
Rp          = Fnorm(b - Amap(Xopt))/(1+Fnorm(b));
Rd          = Fnorm(ops(ops(C,'-',ATmap(yopt)),'-',Sopt))/(1+Fnorm(C));
pobj        = blktrace(blk,C,Xopt);
dobj        = b'*yopt;
gap         = abs(pobj - dobj)/(1+abs(pobj) + abs(dobj));

%% round a rank-one solution
Xmom        = Xoptmat{1};
[V,~]       = sorteig(Xmom);
v           = V(:,1)/V(1,1);
x_est       = v(2:1+POP.d);
x_est       = x_est / norm(x_est);
f_est       = POP.f.coefficient * ( prod(x_est.^(POP.f.degmat),1) )';

S_mineig    = mineig( smat(blk, ops(C,'-',ATmap(yopt)) ) );
S_mineig_1  = min(S_mineig,0);
M           = SDP.M;
f_lb        = SDP.b'*yopt + M'*S_mineig_1;
eta         = abs(f_est - f_lb)/(1+abs(f_est)+abs(f_lb));

info.Rp     = Rp;
info.Rd     = Rd;
info.Rg     = gap;
info.Rs     = eta;
info.pobj   = pobj;
info.dobj   = dobj;
info.f_est      = f_est;
info.S_mineig   = S_mineig;
info.f_lb       = f_lb;
info.x_est      = x_est;

fprintf('\n========== Performance of Quartic Opt Over Sphere ==========\n')
fprintf('Rp: %3.2e, Rd: %3.2e, Rg: %3.2e, Rs: %3.2e.\n',Rp,Rd,gap,eta);
fprintf('f_est: %3.4e, f_lb: %3.4e, pobj: %3.4e, dobj: %3.4e.\n',f_est,f_lb,pobj,dobj);
fprintf('==============================================================\n')

if iscell(path)
    for i = 1:length(path)
        rmpath(genpath(path{i}))
    end
else
    rmpath(genpath(path))
end
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function points = intersectLineCircle(line, circle)
%INTERSECTLINECIRCLE Intersection point(s) of a line and a circle
%
%   INTERS = intersectLineCircle(LINE, CIRCLE);
%   Returns a 2-by-2-by-N array, containing on each row the coordinates of
%   an intersection point for each line-circle pair, i.e. INTERS(:,:,k)
%   contains the intersections between LINE(k,:) and CIRCLE(k,:).
%
%   If a line-circle pair does not intersect, the corresponding results are
%   set to NaN.
%
%   Example
%     % base point
%     center = [10 0];
%     % create vertical line
%     l1 = [center 0 1];
%     % circle
%     c1 = [center 5];
%     pts = intersectLineCircle(l1, c1)
%     pts =
%     10   -5
%     10    5
%     % draw the result
%     figure; clf; hold on;
%     axis([0 20 -10 10]);
%     drawLine(l1);
%     drawCircle(c1);
%     drawPoint(pts, 'rx');
%     axis equal;
%
%   See also
%   lines2d, circles2d, intersectLines, intersectCircles
%
%   References
%   http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
%   http://mathworld.wolfram.com/Circle-LineIntersection.html
%

% ------
% Author: David Legland, david.legland@inra.fr
% Author: JuanPi Carbajal <ajuanpi+dev@gmail.com>
% Created: 2011-01-14,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% HISTORY
% 2011-06-06 fix bug in delta test
% 2017-05-05 included some suggestions from code by JuanPi Carbajal <ajuanpi+dev@gmail.com>
% 2017-08-08 update doc

% check size of inputs
nLines = size(line, 1);
nCircles = size(circle, 1);
if nLines ~= nCircles
  error ('matGeom:geom3d:invalidArguments', ...
      'Requires same number of lines and circles');
end

% center parameters
center = circle(:, 1:2);
radius = circle(:, 3);

% line parameters
dp = line(:, 1:2) - center;
vl = line(:, 3:4);

% coefficient of second order equation
a = sum(line(:, 3:4).^2, 2);
b = 2*sum(dp .* vl, 2);
c = sum(dp.^2, 2) - radius.^2;

% discriminant
delta = b .^ 2 - 4 * a .* c;

points = nan(2, 2, nCircles);

valid = delta >= 0;

if any(valid)
    % compute roots
    u = bsxfun(@plus, -b(valid), bsxfun(@times, [-1 1], sqrt(delta(valid))));
    u = bsxfun(@rdivide, u, a(valid)) / 2;

    if nCircles == 1
        points = [...
            line(1:2) + u(:,1) .* line(3:4); ...
            line(1:2) + u(:,2) .* line(3:4)];
    else
        tmp = [...
            line(valid, 1:2) + u(:,1) .* line(valid, 3:4) ...
            line(valid, 1:2) + u(:,2) .* line(valid, 3:4)].';
	    points(:, :, valid) = permute(reshape(tmp, [2, 2, nCircles]), [2 1 3]);
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips