text:
function browse_mdl(bpath)
% browse_mdl(bpath), Browse models
%
% A quick and dirty simulink model file browser
% Copies pictures of the simulink models files specified by bpath into a new word document with the model name
% If bpath is not specified it lists all the simulink files in the current directory.

if nargin < 1
    files=dir('*.mdl');
else
    [pathstr,name] = fileparts(bpath);
    if exist(bpath,'dir')
        name = [name '\*'];
    end
    ext = '.mdl';
    files=dir(fullfile(pathstr,[name ext]));
end

if ~isempty(files)
    whandle = actxserver('word.application');
    whandle.Visible = 1;
    Add(whandle.Documents);
    sel=whandle.Selection;
    set(sel.ParagraphFormat,'Alignment','wdAlignParagraphLeft');
    set(sel.PageSetup,'LeftMargin',28.35,'RightMargin',28.35,'BottomMargin',28.35*1.43);
    for index=1:size(files,1)
        title=files(index).name;
        eval(title(1:end-4));
        print('-dmeta',['-s', title(1:end-4)]);
        bdclose;
        close('all');
        TypeText(sel,title);
        set(sel.ParagraphFormat,'KeepWithNext',-1)
        TypeParagraph(sel);
        Paste(sel);
        set(sel.ParagraphFormat,'KeepWithNext',0)
        TypeParagraph(sel);
    end
    release(whandle);
end

return


text:
function phase = ph(params,freq)
% This function computes the phase of an analog allpass network with 2 legs and N sections.
% The phase is evalutated at the set of frequencies given in the vector: freq.
% function phase = ph(params,freq)
%
% Inputs:
% params (Nx2)	Matrix of parameters (1/RC for 2 legs and N sections)
% freq	(1xM)	frequency vector (in Hz)
%
% Output:
% phase	(2xM)	output phase vector in deg.

[N,L] = size(params);
[dummy,M] = size(freq);

temp = zeros(N,2*M);		% set array size
%temp(:) = 180 - (360/pi)*atan(2*pi*(1./(10.^params(:)))*freq);	% compute phase matrix
temp(:) = (360/pi)*atan(2*pi*(1./(10.^params(:)))*freq);	% compute phase matrix (ignor the fact that each stage inverts signal)
phase = zeros(2,M);			% set array size
phase(:) = sum(temp,1);		% sum up the phase over N cascaded sections
return

author: Sable

repo: mcbench-benchmarks


text:
function []=plot_sb_baselines(ix)
%PLOT_SB_BASELINES plot the small baselines in small_baselines.list
%Optional an input argument ix can be specified containing a vector of the
%small baseline interferograms to keep which will be plotted in the baseline plot.
%
%   Andy Hooper, June 2007
%
%   ======================================================================
%   09/2010 AH: Add option to plot in MERGED directory
%   09/2010 AH: For SMALL_BASELINES/MERGED don't plot dropped ifgs
%   12/2012 DB: Added meaning of ix to the syntax of the code
%   04/2013 DB: Command variable
%   03/2014 DB: Suppress command line output
%   01/2017 DB: Relax the name for SB such it could have a prefix or suffix
%   01/2017 DB: add fontsize
%   ======================================================================

% 05/2019 MG: Modified to plot colorize network based on interferogram standard
% deviation

if nargin <1
   ix=[];
end

fontsize = 16;

currdir=pwd;
dirs=strread(currdir,'%s','delimiter','/');

if ~isempty(findstr(dirs{end},'SMALL_BASELINES'))

    if strcmp(dirs{end},'SMALL_BASELINES')
        [a,b] = system(['\ls -d [1,2]* | sed ''' 's/_/ /''' ' > small_baselines.list']);
    else
        cd ../SMALL_BASELINES
        [a,b] = system(['\ls -d [1,2]* | sed ''' 's/_/ /''' ' > ' currdir filesep 'small_baselines.list']);
        cd(currdir)
    end
    load ../psver
    psname=['../ps',num2str(psver)];
    small_baseline_flag='y';

elseif strcmp(dirs{end},'MERGED')
    cd ../SMALL_BASELINES
    [a,b] = system(['\ls -d [1,2]* | sed ''' 's/_/ /''' ' > ../MERGED/small_baselines.list']);
    cd ../MERGED
    load ../psver
    psname=['../ps',num2str(psver)];
    small_baseline_flag='y';
else
    load psver
    psname=['ps',num2str(psver)];
    small_baseline_flag='n';
end

sb=load('small_baselines.list');
n_ifg=size(sb,1);
if small_baseline_flag=='y' & isempty(ix) & exist('./parms.mat','file')
    drop_ifg_index=getparm('drop_ifg_index');
    if ~isempty(drop_ifg_index)
       ix=setdiff([1:n_ifg],drop_ifg_index);
    end
end

if ~isempty(ix)
    sb=sb(ix,:);
else
    ix=1:size(sb,1);
end

ps=load(psname);

% loading the std information
    ifgstdname=['ifgstd',num2str(2)];
    if exist([ifgstdname,'.mat'],'file')
       stdin=load(ifgstdname);
       ifg_std=stdin.ifg_std;
       clear stdin
    else
    ifg_std=zeros(ps.n_ifg,1);
    end

    if ~isempty(ix)
    ifg_std = ifg_std(ix,:);
    else
    ix=1:size(ifg_std,1);
    end

n_ifg=size(sb,1);
[yyyymmdd,I,J]=unique(sb);
ifg_ix=reshape(J,n_ifg,2);
x=ifg_ix(:,1);
y=ifg_ix(:,2);

day=str2num(datestr(ps.day,'yyyymmdd'));
[B,I]=intersect(day,yyyymmdd);

x=I(x);
y=I(y);

%create colormap
%colors
C1 = [1 0 0]; %red
C2 =[1 1 1];  %white
C3 = [0 0 1];  %blue

invert_flag = 0;

%invert colormap
if invert_flag ==1
  CM1=[linspace(C1(1),C2(1),100)',linspace(C1(2),C2(2),100)', linspace(C1(3),C2(3),100)'];
  CM2=[linspace(C2(1),C3(1),100)',linspace(C2(2),C3(2),100)', linspace(C2(3),C3(3),100)'];
  CM = [CM1;CM2];
elseif invert_flag ==0
  CM1=[linspace(C2(1),C1(1),100)',linspace(C2(2),C1(2),100)', linspace(C2(3),C1(3),100)'];
  CM2=[linspace(C3(1),C2(1),100)',linspace(C3(2),C2(2),100)', linspace(C3(3),C2(3),100)'];
  CM = [CM2;CM1];
else
end

%cm = colormap('parula'); % take your pick (doc colormap)
%cm = CM;
%cm = interp1(linspace(min(ifg_std),max(ifg_std),length(cm)),cm,ifg_std); % map color to y values
%cm = uint8(cm'*255); % need a 4xN uint8 array
%cm(4,:) = 255; % last column is transparency

h_fig=figure('name','Baseline plot','position',[ 440   170   821   628]);
for i=1:length(x)
    hold on
    x1 = [ps.day(x(i)),ps.day(y(i))]';
    y1 = [ps.bperp(x(i)),ps.bperp(y(i))]';
    xx = [x1 x1];
    yy = [y1 y1];
    zz = zeros(size(xx));
    cc = [ifg_std(i) ifg_std(i)];
    l=surf(xx,yy,zz,[cc;cc],'EdgeColor','interp','FaceColor','none');
    %l=plot([ps.day(x(i)),ps.day(y(i))],[ps.bperp(x(i)),ps.bperp(y(i))],'-');
    text((ps.day(x(i))+ps.day(y(i)))/2,(ps.bperp(x(i))+ps.bperp(y(i)))/2,num2str(ix(i)),'fontsize',fontsize-7);
    set(l,'linewidth',2)
    %set(l.Edge,'ColorBinding','interpolated','ColorData',cm)
    %drawnow i
end

if length(ix) ~= 1
colormap(CM);
cb = colorbar;
ylabel(cb,'Ifg std [deg]');
else
end
shading flat
view(2)
p=plot(ps.day,ps.bperp,'ko');
set(p,'markersize',6,'linewidth',2)
hold off
box on
set(gca,'fontsize',fontsize)
datetick('x',12)
ylabel('B_{perp}','fontsize',fontsize)
try
    print(h_fig,'-dpng','baseline_plot.png')
catch
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% Copyright (C) 1993-2013, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for MATLAB (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% We will show how to acquire images from files, movies, the web and attached
% video cameras.

% We start by opening a movie file, a standard MPG format movie file
camera = Movie('traffic_sequence.mpg');
about camera
% and we see that the result is a Movie object

% The displayed value of this object
camera
% shows the frame dimension, frame rate, number of frames etc.

% To grab a frame from the movie
im = camera.grab();
% and then display it
idisp(im)

% We can display subsequent frames from the movie
idisp( camera.grab() );
idisp( camera.grab() );
idisp( camera.grab() );
idisp( camera.grab() );

% Finally we destroy the camera object which deallocates resources
clear camera

% We can grab an image from an internet webcam if we know its URL
camera = AxisWebCamera('http://wc2.dartmouth.edu')
about camera
% and we see that the result is an AxisWebCamera  object

% We can display a frame from that image source
idisp( camera.grab() );
clear camera

% We can grab an image from Google Maps
ev = EarthView();
% which does require you to have a Google key (available for free) and given
% in the environment GOOGLE_KEY
idisp( ev.grab(-27.475722,153.0285, 17) )
% Alternatively we can a geographic lookup service to find a place by name
% rather than latitude/longitude
idisp( ev.grab('brisbane', 14) )

% Finally, if you have a camera in, or attached to, your computer you can grab
% a live image frame.
% First create a VideoCamera object
camera = VideoCamera();
idisp( camera.grab() )
% and close the camera
clear camera

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [ y, j ] = j_carry_hebrew ( y, j )

%*****************************************************************************80
%
%% J_CARRY_HEBREW carries year-days to years in a Hebrew date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, a YJ date.
%
  while ( 1 )

    days = year_length_hebrew ( y );

    if ( j < days )
      break
    end

    j = j - days;
    y = y + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Please cite the following article when using this source code:
%

%
% Author:
% - Mehrtash Harandi (mehrtash.harandi at gmail dot com)
%
% This file is provided without any warranty of
% fitness for any purpose. You can redistribute
% this file and/or modify it under the terms of
% the GNU General Public License (GPL) as published
% by the Free Software Foundation, either version 3
% of the License or (at your option) any later version.

clear;
clc;

Beta = 1e1;
SPARSE_SOLVER = 'SPAMS'; %options are CVX and SPAMS

load toy_data; %The toy data is indeed a split of Brodatz-16c dataset.

TestSet;
TrainSet;

%Generate Kernel
galleryKernel = exp(-Beta*Stein_Divergence(TrainSet.X,TrainSet.X));
probeKernel = exp(-Beta*Stein_Divergence(TestSet.X,TrainSet.X))';

%Normalize dictionary
KD = galleryKernel./repmat(sqrt(sum(galleryKernel.^2)),[size(galleryKernel,1) 1]);

%
KX = probeKernel./repmat(sqrt(sum(probeKernel.^2)),[size(probeKernel,1) 1]);
L1 = size(probeKernel,2);
if (strcmpi(SPARSE_SOLVER,'CVX'))
    SR_Lambda = 5e-1;
    A = KD;
    Number_Visual_Codes = size(A,2);
    scX = zeros(Number_Visual_Codes,L1);
    for tmpC1=1:L1
        if (mod(tmpC1,10) == 1)
            fprintf('Working on sample %d/%d\n',tmpC1,L1);
        end
        cvx_begin quiet
        variable tmp_Sparse_X(Number_Visual_Codes);
        minimize( -2*tmp_Sparse_X'*KX(:,tmpC1)+tmp_Sparse_X'*KD*tmp_Sparse_X+SR_Lambda*norm(tmp_Sparse_X,1));
        % if you are interested in positive sparse codes, uncomment the following two lines
        %             subject to
        %                 tmp_Sparse_X >= 0;
        cvx_end
        scX(:,tmpC1) = tmp_Sparse_X;
    end %endfor
elseif (strcmpi(SPARSE_SOLVER,'SPAMS'))
    SR_Lambda = 1e-2;
    [KD_U,KD_D,~] = svd(KD);
    A = diag(sqrt(diag(KD_D)))*KD_U';
    D_Inv = KD_U*diag(1./sqrt(diag(KD_D)));
    KX = D_Inv'*KX;
    scX = full(mexLasso(KX,A,struct('mode',2,'lambda',SR_Lambda,'lambda2',0)));
else
    error('Sparse solver is not defined. This version accepts CVX or SPAMS');
end

%% Classification
qLabel = Classify_SRC(A,TrainSet.y,scX,KX);

CRR = sum((qLabel - TestSet.y == 0))/L1;

fprintf('Clasiffication Ratio: %f\n', CRR);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% TEST_cartPoleDynamics.m
%
% A simple script to run a forward simulation of the cart-pole dynamics, to
% make sure that everything looks reasonable.
clc; clear;

p.m1 = 1;
p.m2 = 1;
p.g = 9.81;
p.l = 1;

z0 = [...
    0;              %Initial horizontal position
    (pi/180)*170;   %Initial angle of the pendulum (from stable config)
    0;              %Initial velocity of cart
    0];             %initial rotation rate of pendulum

%Simple PD controller
u = @(z)( -[ 10*z(1) + 2*z(3) ; 10*z(2) + 2*z(4) ]  );

tSpan = [0;2];   %Time span for simulation

dynFun = @(t,z)( cartPoleDynamics(z,u(z),p) );   %Function handle to give ode45
sol = ode45(dynFun,tSpan,z0);  %Run simulation

%Extract uniformly spaced samples via high-order interpolation
t = linspace(tSpan(1),tSpan(2),250);
z = deval(sol,t);

% Make some simple plots:
figure(1); clf;

subplot(2,2,1);
plot(t,z(1,:));
xlabel('t')
ylabel('x')

subplot(2,2,2);
plot(t,z(2,:));
xlabel('t')
ylabel('q')

subplot(2,2,3);
plot(t,z(3,:));
xlabel('t')
ylabel('dx')

subplot(2,2,4);
plot(t,z(4,:));
xlabel('t')
ylabel('dq')

author: MatthewPeterKelly

repo: dscTutorials


text:
function islands=bwislands(img)
%
% islands=bwislands(img)
%
% return the indices of non-zero elements in a 2D or 3D image
% grouped by connected regions in a cell array
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%	 img: a 2D or 3D array
% output:
%	 islands: a cell array, each cell records the indices
%		  of the non-zero elements in img for a connected
%		  region (or an island)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

img=logical(1-img);
idx=find(1-img(:));
islands={};

count=1;
while(length(idx))
        [I,J]=ind2sub(size(img),idx(1));
	imgnew=imfill(img,[I,J]);
	islands{count}=find(imgnew~=img);
	count=count+1;
	img=imgnew;
	idx=find(1-img(:));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%This Matlab script can be used to reproduce Figure 4.11 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Number of BS antennas
M = 100;

%Define the pilot reuse factor
f = 2;

%Select the number of setups with random UE locations
nbrOfSetups = 10;

%Select the number of channel realizations per setup
nbrOfRealizations = 1000;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Noise figure at the BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the exact Gaussian local scattering model
accuracy = 1;

%Range of angular standard deviation in the local scattering model (in degrees)
ASDdegRange = [1e-4 5:5:50];

%Prepare to save simulation results
sumSE_MR_theorem = zeros(length(ASDdegRange),nbrOfSetups);
sumSE_MR_UatF_corollary = zeros(length(ASDdegRange),nbrOfSetups);
sumSE_MR_UatF_alt1 = zeros(length(ASDdegRange),nbrOfSetups);
sumSE_MR_UatF_alt2 = zeros(length(ASDdegRange),nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup and many ASDs
    [R,channelGaindB] = functionExampleSetup(L,K,M,accuracy,ASDdegRange);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Go through all ASDs
    for spr = 1:length(ASDdegRange)

        %Output simulation progress
        disp([num2str(spr) ' ASDs out of ' num2str(length(ASDdegRange))]);

        %Generate channel realizations with estimates and estimation error
        %correlation matrices
        [Hhat,C,tau_p,~,H] = functionChannelEstimates(R(:,:,:,:,:,spr),channelGainOverNoise,nbrOfRealizations,M,K,L,p,f);

        %Compute SE using Monte Carlo simulations
        [signal,interference,prelogFactor,SE_MR] = functionComputeSE_UL_MR_MC(Hhat,H,C,tau_c,tau_p,nbrOfRealizations,M,K,L,p);

        %Save results
        sumSE_MR_theorem(spr,n) = mean(sum(SE_MR,1));

        sumSE_MR_UatF_corollary(spr,n) = prelogFactor*mean(sum(log2(1+signal(:,:,1)./interference(:,:,1)),1));

        sumSE_MR_UatF_alt1(spr,n) = prelogFactor*mean(sum(log2(1+signal(:,:,2)./interference(:,:,2)),1));

        sumSE_MR_UatF_alt2(spr,n) = prelogFactor*mean(sum(log2(1+signal(:,:,3)./interference(:,:,3)),1));

        %Delete large matrices
        clear Hhat H C;

    end

    %Delete large matrices
    clear R;

end

%Compute the average sum SE over setups
sumSE_MR_theorem_opt = mean(sumSE_MR_theorem,2);
sumSE_MR_UatF_corollary_opt = mean(sumSE_MR_UatF_corollary,2);
sumSE_MR_UatF_alt1_opt = mean(sumSE_MR_UatF_alt1,2);
sumSE_MR_UatF_alt2_opt = mean(sumSE_MR_UatF_alt2,2);

%% Plot the simulation results
figure(1);
hold on; box on;

plot(ASDdegRange,sumSE_MR_theorem_opt,'rs--','LineWidth',1);
plot(ASDdegRange,sumSE_MR_UatF_alt2_opt,'b*-','LineWidth',1);
plot(ASDdegRange,sumSE_MR_UatF_alt1_opt,'kd:','LineWidth',1);
plot(ASDdegRange,sumSE_MR_UatF_corollary_opt,'kp--','LineWidth',1);

xlabel('ASD $$(\sigma_{\varphi})$$ [degree]','Interpreter','latex');
ylabel('Average sum SE [bit/s/Hz/cell]');
legend('Original bound (Th. 4.1)','UatF: h/||h||^2','UatF: h/||h||','UatF: h (Corr. 4.5)','Location','NorthEast');
ylim([0 50]);

author: emilbjornson

repo: massivemimobook


text:
function face = compute_delaunay(vertex)

% compute_delaunay - compute faces of a Delaunay triangulation
%
%   face = compute_delaunay(vertex);
%
%   Copyright (c) 2008 Gabriel Peyre

if size(vertex,1)>size(vertex,2)
    vertex = vertex';
end

face = delaunay(vertex(1,:), vertex(2,:))';

author: gpeyre

repo: matlab-toolboxes


text:
function [soln,eqn,info] = StokesP2P1(node,elem,bdFlag,pde,option)
%% STOKESP2P1 Stokes equation: P2-P1 Taylor-Hood elements.
%
%   [soln,eqn,info] = STOKESP2P1(node,elem,pde,bdFlag) use quadratic and piceswise
%   linear elements to approximate velocity u and pressure p, repectively.
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%   with
%       Dirichlet boundary condition        u = g_D  on \Gamma_D,
%       Neumann boundary condition du/dn - np = g_N  on \Gamma_N.
%
%   It is a choice of option.fem in Stokes. Please read Stokes for more
%   information on the input and output.
%
% See also Stokes, Poisson, StokesP2P1
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end

t = cputime;
%% Construct Data Structure
[elem2dof,edge,bdDof] = dofP2(elem);
N = size(node,1);  NT = size(elem,1);  Nu = N+size(edge,1);   Np = N;

%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix for Laplace operator
% generate sparse pattern
ii = zeros(21*NT,1); jj = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        index = index + NT;
    end
end
% quadrature points
if ~isfield(pde,'nu'), pde.nu = []; end
if ~isfield(option,'quadorder')
    % constant viscosity
    option.quadorder = 2;        % default order
    if ~isempty(pde.nu) && isnumeric(pde.nu) % numerical viscosity
        option.quadorder = 3;    % exact for linear diffusion coefficient
    end
end
[lambda, w] = quadpts(option.quadorder);
nQuad = size(lambda,1);
% compute non-zeros
sA = zeros(21*NT,nQuad);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    index = 0;
    for i = 1:6
        for j = i:6
            Aij = 0;
            if isempty(pde.nu) || isnumeric(pde.nu)
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
            else
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
            end
            if ~isempty(pde.nu) && (pde.nu~=1)
                Aij = pde.nu*Aij;
            end
            Aij = Aij.*area;
            sA(index+1:index+NT,p) = Aij;
            index = index + NT;
        end
    end
end
sA = sum(sA,2);
% assemble the matrix
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Nu,Nu);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Nu,Nu);
A = A + AU + AU';
A = blkdiag(A,A);
clear Aij ii jj sA

%% Assemble the matrix for divergence operator
Dx = sparse(Np,Nu);
Dy = sparse(Np,Nu);
[lambda, w] = quadpts(2); % (div(P2), P1) is P2
nQuad = size(lambda,1);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    for i = 1:6
        for j = 1:3
            Dxij = 0;
            Dyij = 0;
            Dxij = Dxij + w(p)*Dphip(:,1,i).*lambda(p,j);
            Dyij = Dyij + w(p)*Dphip(:,2,i).*lambda(p,j);
            Dx = Dx + sparse(elem(:,j),double(elem2dof(:,i)),Dxij.*area,Np,Nu);
            Dy = Dy + sparse(elem(:,j),double(elem2dof(:,i)),Dyij.*area,Np,Nu);
        end
    end
end
B = -[Dx Dy];
clear Dxij Dyij Dx Dy Dphip

%% Assemble right hand side
f1 = zeros(Nu,1);
f2 = zeros(Nu,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 4;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,weight] = quadpts(option.fquadorder);
    % basis values at quadrature points
    phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
    phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
    phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
    phi(:,4) = 4*lambda(:,2).*lambda(:,3);
    phi(:,5) = 4*lambda(:,3).*lambda(:,1);
    phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    nQuad = size(lambda,1);
    ft1 = zeros(NT,6);
    ft2 = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        % function values at quadrature points
        fp = pde.f(pxy);
        % evaluate fp outside.
        for j = 1:6
            ft1(:,j) = ft1(:,j) + fp(:,1).*phi(p,j)*weight(p);
            ft2(:,j) = ft2(:,j) + fp(:,2).*phi(p,j)*weight(p);
        end
    end
    ft1 = ft1.*repmat(area,1,6);
    ft2 = ft2.*repmat(area,1,6);
    f1 = accumarray(elem2dof(:),ft1(:),[Nu 1]);
    f2 = accumarray(elem2dof(:),ft2(:),[Nu 1]);
end
clear phi ft1 ft2

%% Boundary Conditions
[AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesP2P1;

%% Record assembeling time
assembleTime = cputime - t;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(ufreeDof), return; end
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if length(f)+length(g) <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else          % Multigrid-type  solver for large size systems
        option.solver = 'asmg';
    end
end
solver = option.solver;

% solve
switch solver
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'direct'
        t = cputime;
        bigA = [AD, B'; ...
                B, sparse(Np,Np)];
        bigF = [f; g];
        bigu = [u; p];
        bigFreeDof = [ufreeDof; 2*Nu+pDof];
        bigu(bigFreeDof) = bigA(bigFreeDof,bigFreeDof)\bigF(bigFreeDof);
        u = bigu(1:2*Nu);
        p = bigu(2*Nu+1:end);
        residual = norm(bigF - bigA*bigu);
        info = struct('solverTime',cputime - t,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
%         option.tol = Np^(-2);
        option.solver  = 'WCYCLE';
        [u(ufreeDof),p,info] = mgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                        u(ufreeDof),p,elem,ufreeDof,option);
    case 'asmg'
        [u(ufreeDof),p,info] = asmgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                          u,p,node,elem,bdFlag,ufreeDof,option);
end

%% Post-process
if length(pDof) ~= Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(mean(p(elem),2).*area)/sum(area);
    p = p - c;
end

%% Output
soln = struct('u',u,'p',p);
eqn = struct('A',AD,'B',BD,'Lap',A,'f',f,'g',g,...
             'edge',edge,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdStokesP2P1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesP2P1
    %% Boundary condition of Stokes equation: P2-P0 elements

    %% Initial set up
%     f = [f1; f2];    % set in Neumann boundary condition
    g = zeros(Np,1);
    u = zeros(2*Nu,1);
    p = zeros(Np,1);
    ufreeDof = (1:Nu)';
    pDof = (1:Np)';

    if ~exist('bdFlag','var'), bdFlag = []; end
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Find Dirichlet dof and modify the matrix
    % Find Dirichlet boundary dof: fixedDof and pDof
    isFixedDof = false(Nu,1);
    if ~isempty(bdFlag)       % case: bdFlag is not empty
        elem2edge = elem2dof(:,4:6)-N;
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;   % nodes of all D-edges
        isFixedDof(N + find(isDirichlet')) = true;% dof on D-edges
        fixedDof = find(isFixedDof);
        ufreeDof = find(~isFixedDof);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        fixedDof = bdDof;
        isFixedDof(fixedDof) = true;
        ufreeDof = find(~isFixedDof);
    end
    if isempty(fixedDof) % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        fixedDof = 1;
        ufreeDof = (2:Nu)';    % eliminate the kernel by enforcing u(1) = 0;
    end

    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedDof,fixedDof)=I, AD(fixedDof,ufreeDof)=0, AD(ufreeDof,fixedDof)=0.
    % BD(:,fixedDof) = 0 and thus BD'(fixedDof,:) = 0.
    bdidx = zeros(2*Nu,1);
    bdidx([fixedDof; Nu+fixedDof]) = 1;
    Tbd = spdiags(bdidx,0,2*Nu,2*Nu);
    T = spdiags(1-bdidx,0,2*Nu,2*Nu);
    AD = T*A*T + Tbd;
    BD = B*T;

    %% Part 2: Find boundary edges and modify the right hand side f and g
    % Find boundary edges: Neumann and Robin
    Neumann = []; Robin = []; %#ok<*NASGU>
    if ~isempty(bdFlag)
        isNeumann(elem2edge((bdFlag(:)==2)|(bdFlag(:) == 3))) = true;
        isRobin(elem2edge(bdFlag(:)==3)) = true;
        Neumannidx = find(isNeumann);
        Neumann   = edge(isNeumann,:);
        Robin     = edge(isRobin,:);
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        Neumann = edge(bdDof>N,:);
        if ~isempty(pde.g_R)
            Robin = Neumann;
        end
    end

    % Neumann boundary condition
    if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        [lambda,w] = quadpts1(3);
        nQuad = size(lambda,1);
        % quadratic bases (1---3---2)
        bdphi(:,1) = (2*lambda(:,1)-1).*lambda(:,1);
        bdphi(:,2) = (2*lambda(:,2)-1).*lambda(:,2);
        bdphi(:,3) = 4*lambda(:,1).*lambda(:,2);
        % length of edge
        ve = node(Neumann(:,1),:) - node(Neumann(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        % update RHS
        gex = zeros(size(Neumann,1),2);   % x-component
        gey = zeros(size(Neumann,1),2);   % y-component
        for pp = 1:nQuad
            pxy = lambda(pp,1)*node(Neumann(:,1),:)+lambda(pp,2)*node(Neumann(:,2),:);
            gp = pde.g_N(pxy);
            gex(:,1) = gex(:,1) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,1);
            gex(:,2) = gex(:,2) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,2);
            gey(:,1) = gey(:,1) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,1);
            gey(:,2) = gey(:,2) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,2);
            f1(N+Neumannidx) = f1(N+Neumannidx) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,3); % interior bubble
            f2(N+Neumannidx) = f2(N+Neumannidx) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,3); % interior bubble
        end
        f1(1:N) = f1(1:N) + accumarray(Neumann(:), gex(:),[N,1]);
        f2(1:N) = f2(1:N) + accumarray(Neumann(:), gey(:),[N,1]);
    end
    f = [f1; f2];
    % The case non-empty Neumann but g_N=[] corresponds to the zero flux
    % boundary condition on Neumann edges and no modification is needed.

    % Dirichlet boundary conditions
    if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
        u1 = zeros(Nu,1);
        u2 = zeros(Nu,1);
        idx = (fixedDof > N);              % index of edge dof
        uD = pde.g_D(node(fixedDof(~idx),:));  % bd value at vertex dofs
        u1(fixedDof(~idx)) = uD(:,1);
        u2(fixedDof(~idx)) = uD(:,2);
        bdEdgeIdx = fixedDof(idx)-N;
        bdEdgeMid = (node(edge(bdEdgeIdx,1),:)+node(edge(bdEdgeIdx,2),:))/2;
        uD = pde.g_D(bdEdgeMid);         % bd values at middle points of edges
        u1(fixedDof(idx)) = uD(:,1);
        u2(fixedDof(idx)) = uD(:,2);
        u = [u1; u2]; % Dirichlet bd condition is built into u
        f = f - A*u;  % bring affect of nonhomgenous Dirichlet bd condition to
        g = g - B*u;  % the right hand side
        g = g - mean(g);
        f(fixedDof) = u1(fixedDof);
        f(fixedDof+Nu) = u2(fixedDof);
    end
    % The case non-empty Dirichlet but g_D=[] corresponds to the zero Dirichlet
    % boundary condition and no modification is needed.

    % modfiy pressure dof for pure Dirichlet
    if isempty(Neumann)
        pDof = (1:Np-1)';
    end

    ufreeDof = [ufreeDof; Nu+ufreeDof];
    end % end of function getbdStokesP2P1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

repo: ifem


text:
function [fMshift, fStretch, fMTransFit, fARate] = calc_Shiftstretch(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, nCalculateMC, fMcOrg, fBValueOrg)
% function [fMshift, fStretch, fMTransFit] = plot_Shiftstretch(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, nCalculateMC, fMcOrg, fBValueOrg)
% ----------------------------------------------------------------------------------------------------------------------------------------------
% Function to calculate shift and stretch: Mnew = c*Mold+ dM;
%
% Incoming variables:
% mCatalog     : current earthquake catalog
% bTimePeriod   : Use catalog from beginning to end (0), use time periods (1)
% fSplitTime   : Splittime of catalog
% fTimePeriod  : Time period in decimal years
% nCalculateMC : Method to determine Mc (1-5 see help calc_Mc)
% fMcOrg       : Mc of entire catalog
%
% Outgoing variable:
% fMshift    : magnitude shift
% fStretch   : stretch factor
% fMTransFit : Goodness of fit by modeling second period manipulating the first period
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 09.10.02

% Track of changes:
% 21.08.02: Solved stability  problem (Error: vFMD(1,1) index exceeds matrix dimension) by adding
%           ~isempty(vFMD) & ~isempty(vFMDSecond) into if statement
% 03.10.02: Still little problems with mRes

% Track variables
% fProbability: log likelihood score

global fProbability

[mFirstCatalog, mSecondCatalog, fFirstPeriodExact, fSecondPeriodExact, fFirstPeriod,...
        result.fSecondPeriod] = ex_SplitCatalog(mCatalog, fSplitTime, bTimePeriod, fTimePeriod, bTimePeriod, fTimePeriod);
% Time periods for normalization
if bTimePeriod == 0
    fPeriod1 = (max(mFirstCatalog(:,3))-min(mFirstCatalog(:,3)));
    fPeriod2 = (max(mSecondCatalog(:,3))-min(mSecondCatalog(:,3)));
else
    fPeriod1 = fTimePeriod;
    fPeriod2 = fTimePeriod;
end

% Create the frequency magnitude distribution vectors for the two time
% periods and entire catalog
[vFMD, vNonCFMD] = calc_FMD(mFirstCatalog);
[vFMDSecond, vNonCFMDSecond] = calc_FMD(mSecondCatalog);

% Calculate magnitude of completeness
fMc = calc_Mc(mFirstCatalog, nCalculateMC);
fMcSecond = calc_Mc(mSecondCatalog, nCalculateMC);

if (~isempty(fMc) & ~isempty(fMcSecond) & ~isnan(fMc) & ~isnan(fMcSecond) & ~isempty(vFMD) & ~isempty(vFMDSecond))
    % First period
    [nIndexLo, fMagHi, vSel, vMagnitudes] = fMagToFitBValue(mFirstCatalog, vFMD, fMc);
    % Calculate the b-value etc. for M > Mc
    [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemagMag(mFirstCatalog(vSel,:));
    vPoly = [-1*fBValue fAValue];
    fBFunc = 10.^(polyval(vPoly, vMagnitudes));
    %% Second period
    [nIndexLoSecond, fMagHiSecond, vSelSecond, vMagnitudesSecond] = fMagToFitBValue(mSecondCatalog, vFMDSecond, fMcSecond);
    % Calculate the b-value etc. for M > Mc
    [fMeanMagSecond, fBValueSecond, fStdDevSecond, fAValueSecond] = calc_bmemagMag(mSecondCatalog(vSelSecond,:));
    vPolySecond = [-1*fBValueSecond fAValueSecond];
    fBFuncSecond = 10.^(polyval(vPolySecond, vMagnitudesSecond));
    % Determine magnitude shift
    fMintercept = 1/fBValueSecond*(fAValueSecond-log10(vFMD(2,nIndexLo)));
    fMshift = fMintercept - vFMD(1,nIndexLo);

    sPer1 = ['Period 1: a = ' num2str(fAValue) ' const. b = ' num2str(fBValue) ' Mc = ' num2str(fMc)];
    sPer2 = ['Period 2: a = ' num2str(fAValueSecond) ' const. b = ' num2str(fBValueSecond)...
            ' Mc = ' num2str(fMcSecond)];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Max. likelihood estimations for a and b using different methods for the
    %% SECOND time period

    % Model 1: Variable a and b-value
    mControlM1 = [min(mSecondCatalog(:,3)) fMcSecond fMcSecond 0.1];
    [fBValue2M1, fAValue2M1] = calc_MaxLikelihoodABCombined(mSecondCatalog, mControlM1, 0);
    fP_modelM1 = fProbability;
    nDegFreeM1 = 2;
    sMod1 = ['Model 1: a = ' num2str(fAValue2M1) ' b = ' num2str(fBValue2M1)];

    % Model 2: Variable a, const. b-value from max. likelihood estimation of
    % first period
    vSel = (mSecondCatalog(:,6) >= fMcSecond);
    mSecondComp = mSecondCatalog(vSel,:);
    if isempty(mSecondComp)
        disp('oo');
    end
    fBValue2M2 = fBValue;
    mControlM2 = [min(mSecondCatalog(:,3)) fMcSecond fMcSecond 0.1];
    [fAValue2M2] = calc_MaxLikelihoodA_org(mSecondComp, mControlM2, fBValue2M2, 0);
    % [fAValue2M2] = calc_MaxLikelihoodA(mSecondComp, fBValue2M2)
    fP_modelM2 = fProbability;
    nDegFreeM2 = 1;
    sMod2 = ['Model 2: a = ' num2str(fAValue2M2) ' const. b = ' num2str(fBValue)];

    % Model 3:

    % Model 4: Variable a, const. b-value from max. likelihood estimation of of
    % entire catalog
    vSel = (mSecondCatalog(:,6) >= fMcSecond);
    mSecondComp = mSecondCatalog(vSel,:);
    mControlM4 = [min(mSecondCatalog(:,3)) fMcSecond fMcSecond 0.1];
    [fAValue2M4] = calc_MaxLikelihoodA_org(mSecondComp, mControlM4, fBValueOrg, 0);
    % [fAValue2M4] = calc_MaxLikelihoodA(mSecondComp, fBValueOrg);
    fBValue2M4 = fBValueOrg;
    fP_modelM4 = fProbability;
    nDegFreeM4 = 1;
    sMod4 = ['Model 4: a = ' num2str(fAValue2M4) ' const. b = ' num2str(fBValueOrg)];

    % Model 5: Variable a, const. b-value from max. likelihood estimation of b
    % from second period
    vSel = (mSecondCatalog(:,6) >= fMcSecond);
    mSecondComp = mSecondCatalog(vSel,:);
    mControlM5 = [min(mSecondCatalog(:,3)) fMcSecond fMcSecond 0.1];
    [fAValue2M5] = calc_MaxLikelihoodA_org(mSecondComp, mControlM5, fBValueSecond, 0);
    % [fAValue2M4] = calc_MaxLikelihoodA(mSecondComp, fBValueOrg);
    fBValue2M5 = fBValueSecond;
    fP_modelM5 = fProbability;
    nDegFreeM5 = 1;
    sMod5 = ['Model 5: a = ' num2str(fAValue2M5) ' const. b = ' num2str(fBValue2M5)];

    % Bayesian Information Criteria (BIC) for model decision
    % Take the one with the highest BIC value
    %n_samples = length(fMcSecond:0.1:ceil(max(mSecondCatalog(:,6))));
    n_samples = length(mSecondCatalog(:,6));
    fBIC_1 = 2*fP_modelM1 + 2*log(n_samples)*nDegFreeM1;
    fBIC_2 = 2*fP_modelM2 + 2*log(n_samples)*nDegFreeM2;
    %     fBIC_3 = 2*fP_modelM3 + 2*log(n_samples)*nDegFreeM3;
    fBIC_4 = 2*fP_modelM4 + 2*log(n_samples)*nDegFreeM4;
    fBIC_5 = 2*fP_modelM5 + 2*log(n_samples)*nDegFreeM5;
    sBIC = ['Model comparison: Model 1 BIC: ' num2str(fBIC_1) ';    Model 2 BIC: ' num2str(fBIC_2) '; '...
            ' Model 4 BIC: ' num2str(fBIC_4) ' Model 5 BIC: ' num2str(fBIC_5)];

    %% Stretch factors from b-value ratios
    fStretchM1 = fBValue/fBValue2M1;
    fStretchM2 = fBValue/fBValue;
    %     fStretchM3 = fBValue/fBValue2M3;
    fStretchM4 = fBValue/fBValueOrg;
    fStretchM5 = fBValue/fBValue2M5;

    %% Result - Matrix
    mRes = [];
    mRes = [mRes; fAValue2M1 fBValue2M1 fBIC_1 fStretchM1; fAValue2M2 fBValue2M2 fBIC_2 fStretchM2];
    mRes = [mRes; fAValue2M4 fBValue2M4 fBIC_4 fStretchM4; fAValue2M5 fBValue2M5 fBIC_5 fStretchM5];

    vSel2 = min(find(mRes(:,3) == min(mRes(:,3))));
    vBestModel = mRes(vSel2,:)

    %% Magnitude Shift and stretch for the best model
    fMshift = vBestModel(:,4)*( vBestModel(:,1)- fAValue)/fBValue;
    fStretch = vBestModel(:,4);
    %%% Rate factor froma-value ratio
    fARate = vBestModel(1,1)/fAValue;

else
    disp('fMc, fMcSecond or vFMD / vFMDSecond not derivable');
    fMshift=NaN;
    fStretch = NaN;
end

%%% Goodness of fit determination using formula from Wiemer & Wyss (BSSA, 2000)
mCatModel = mFirstCatalog;
mCatModel(:,6) = fStretch.*mCatModel(:,6)+fMshift;
[vEv_val2 vMags2 vEventsum2 vEvsum_rev2,  vMags_rev2] =calc_cumulsum(mSecondCatalog);
[vEv_valMod vMagsMod vEventsumMod vEvsum_revMod,  vMags_revMod] =calc_cumulsum(mCatModel);
fMTransFit  = sum(abs(vEventsum2-vEventsumMod))/sum(vEventsum2);
fMTransFit = 100-fMTransFit*100;

author: CelsoReyes


text:
clear;
load input.mat; %read the input
load input_adi.mat;
fs=10;

%compute the Allan variance curve
[t_axis, allan, Nrw]=cp_allan_var_v000(input,fs,1,60);
[t_axis_adi, allan_adi, Nrw_adi]=cp_allan_var_v000(input_adi,fs,2,60);

%approximation and rw parameter
apx_lvl=16;
[inp_appx, boun_appx]=lvl_filter_v000(input, apx_lvl,0,0);
rw_coef=cp_wn_pow_v000(input,apx_lvl+1,0);
rw_pow=diag(rw_coef'*rw_coef).^0.5;

apx_lvl_adi=16;
[inp_appx_adi, boun_appx_adi]=lvl_filter_v000(input_adi, apx_lvl_adi,0,0);
rw_coef_adi=cp_wn_pow_v000(input_adi,apx_lvl_adi+1,0);
rw_pow_adi=diag(rw_coef_adi'*rw_coef_adi).^0.5;

%White noise level
whi_coef=cp_wn_pow_v000(input,3,1);
whi_pow=diag(whi_coef'*whi_coef).^0.5;
whi_coef_adi=cp_wn_pow_v000(input_adi,3,1);
whi_pow_adi=diag(whi_coef_adi'*whi_coef_adi).^0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Individual IMU Results %%%%%%%%%%
%%%%% AMS
flt_lvl=9;
[inp_filt, boun_filt]=lvl_filter_v000(input,flt_lvl,0,0);
inputx=inp_filt-inp_appx;
cor_sz=2^16;
inp_corr_x=xcorr(inputx(1,:),cor_sz,'biased');
inp_corr_y=xcorr(inputx(2,:),cor_sz,'biased');
inp_corr_z=xcorr(inputx(3,:),cor_sz,'biased');

plot(0:cor_sz,inp_corr_x(cor_sz+1:2*cor_sz+1));grid;
hold on;
plot(0:cor_sz,inp_corr_y(cor_sz+1:2*cor_sz+1),'r');
plot(0:cor_sz,inp_corr_z(cor_sz+1:2*cor_sz+1),'k');

inp_corr=[inp_corr_x(:,cor_sz+1:2*cor_sz+1);inp_corr_y(:,cor_sz+1:2*cor_sz+1);inp_corr_z(:,cor_sz+1:2*cor_sz+1)];
[mod_h exp_h]=expo_fit_v000([(2^flt_lvl),8000;(2^flt_lvl),15000;(2^flt_lvl),20000], inp_corr,1,1);

flt_lvl=4;
[inp_filt, boun_filt]=lvl_filter_v000(input,flt_lvl,0,0);
inputx=inp_filt-inp_appx;
cor_sz=2^11;
inp_corr_x=xcorr(inputx(1,:),cor_sz,'biased');
inp_corr_y=xcorr(inputx(2,:),cor_sz,'biased');
inp_corr_z=xcorr(inputx(3,:),cor_sz,'biased');

prev_cor_x=exp_h(1,2)*exp_h(1,1).^[0:cor_sz]+exp_h(1,3);
prev_cor_y=exp_h(2,2)*exp_h(2,1).^[0:cor_sz]+exp_h(2,3);
prev_cor_z=exp_h(3,2)*exp_h(3,1).^[0:cor_sz]+exp_h(3,3);
figure;
plot(0:cor_sz,inp_corr_x(cor_sz+1:2*cor_sz+1)-prev_cor_x);grid;
hold on;
plot(0:cor_sz,inp_corr_y(cor_sz+1:2*cor_sz+1)-prev_cor_y,'r');
plot(0:cor_sz,inp_corr_z(cor_sz+1:2*cor_sz+1)-prev_cor_z,'k');

inp_corr=[inp_corr_x(:,cor_sz+1:2*cor_sz+1)-prev_cor_x;inp_corr_y(:,cor_sz+1:2*cor_sz+1)-prev_cor_y;inp_corr_z(:,cor_sz+1:2*cor_sz+1)-prev_cor_z];
[mod_l exp_l]=expo_fit_v000([16,500;16, 1200;16, 1200], inp_corr,1,1);

%compare the allan results
allan_th=arma2allan_v000([1 -mod_h(1,1);1 -mod_l(1,1);1 0;1 -1],[mod_h(1,2);mod_l(1,2);whi_pow(1);rw_pow(1)],fs,t_axis,allan(1,:));
allan_th=arma2allan_v000([1 -mod_h(2,1);1 -mod_l(2,1);1 0;1 -1],[mod_h(2,2);mod_l(2,2);whi_pow(2);rw_pow(2)],fs,t_axis,allan(2,:));
allan_th=arma2allan_v000([1 -mod_h(3,1);1 -mod_l(3,1);1 0;1 -1],[mod_h(3,2);mod_l(3,2);whi_pow(3);rw_pow(3)],fs,t_axis,allan(3,:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% ADI
flt_lvl=9;
[inp_filt, boun_filt]=lvl_filter_v000(input_adi,flt_lvl,0,0);
inputx=inp_filt-inp_appx_adi;
cor_sz=2^16;
inp_corr_x=xcorr(inputx(1,:),cor_sz,'biased');
inp_corr_y=xcorr(inputx(2,:),cor_sz,'biased');
inp_corr_z=xcorr(inputx(3,:),cor_sz,'biased');

figure;
plot(0:cor_sz,inp_corr_x(cor_sz+1:2*cor_sz+1));grid;
hold on;
plot(0:cor_sz,inp_corr_y(cor_sz+1:2*cor_sz+1),'r');
plot(0:cor_sz,inp_corr_z(cor_sz+1:2*cor_sz+1),'k');

inp_corr=[inp_corr_x(:,cor_sz+1:2*cor_sz+1);inp_corr_y(:,cor_sz+1:2*cor_sz+1);inp_corr_z(:,cor_sz+1:2*cor_sz+1)];
[mod_h_adi exp_h_adi]=expo_fit_v000([(2^flt_lvl),18000;(2^flt_lvl),18000;(2^flt_lvl),18000], inp_corr,1,1);

%compare the allan results
allan_th=arma2allan_v000([1 -mod_h_adi(1,1);1 0;1 -1],[mod_h_adi(1,2);whi_pow_adi(1);rw_pow_adi(1)],fs,t_axis,allan_adi(1,:));
allan_th=arma2allan_v000([1 -mod_h_adi(2,1);1 0;1 -1],[mod_h_adi(2,2);whi_pow_adi(2);rw_pow_adi(2)],fs,t_axis,allan_adi(2,:));
allan_th=arma2allan_v000([1 -mod_h_adi(3,1);1 0;1 -1],[mod_h_adi(3,2);whi_pow_adi(3);rw_pow_adi(3)],fs,t_axis,allan_adi(3,:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Combined (Multi) IMU Results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% combine outputs
R=diag([diag(whi_coef'*whi_coef);diag(whi_coef_adi'*whi_coef_adi)]);
H=[sin(45*pi/180) cos(45*pi/180);
   sin(90*pi/180) cos(90*pi/180);
   sin(70*pi/180) cos(70*pi/180)
   sin(10*pi/180) cos(10*pi/180)
   sin(0*pi/180) cos(0*pi/180)
   sin(-10*pi/180) cos(-10*pi/180)];
%standard method
T2=inv(H'*inv(R)*H)*H'*inv(R);
output=T2*[input;input_adi];

%compute the Allan variance curve
[t_axis, allan_out, Nrw]=cp_allan_var_v000(output,fs,3,60);

%approximation and rw parameter
apx_lvl=16;
[out_appx, boun_appx]=lvl_filter_v000(output, apx_lvl,0,0);
rw_coef_out=cp_wn_pow_v000(output,apx_lvl+1,0);
rw_pow_out=diag(rw_coef_out'*rw_coef_out).^0.5;

%White noise level
whi_coef_out=cp_wn_pow_v000(output,3,1);
whi_pow_out=diag(whi_coef_out'*whi_coef_out).^0.5;

%%%%%%% Start modeling
flt_lvl=9;
llim=2^flt_lvl; %autocorr is valid for k>=llim;
[out_filt, boun_filt]=lvl_filter_v000(output,flt_lvl,0,0);
outputx=out_filt-out_appx;
cor_sz=2^15;
out_corr_x=xcorr(outputx(1,:),cor_sz,'biased');
out_corr_y=xcorr(outputx(2,:),cor_sz,'biased');
out_corr_xy=xcorr(outputx(1,:),outputx(2,:),cor_sz,'biased');

figure;
plot(0:cor_sz,out_corr_x(cor_sz+1:2*cor_sz+1));grid;
hold on;
plot(0:cor_sz,out_corr_y(cor_sz+1:2*cor_sz+1),'r');
plot(0:cor_sz,out_corr_xy(cor_sz+1:2*cor_sz+1),'k');
plot(-cor_sz:0,out_corr_xy(1:cor_sz+1),'k');

out_corr=[out_corr_x(:,cor_sz+1:2*cor_sz+1);out_corr_y(:,cor_sz+1:2*cor_sz+1);out_corr_xy(:,cor_sz+1:2*cor_sz+1);fliplr(out_corr_xy(:,1:cor_sz+1));];
[mod_h_out exp_h_out]=expo_fit_v000([ones(3,1)*[llim,20000];[llim 10000]], out_corr,1,1);

%%%lower level
flt_lvl=3;
[out_filt, boun_filt]=lvl_filter_v000(output,flt_lvl,0,0);
outputx=out_filt-out_appx;
cor_sz=2^10;
out_corr_x=xcorr(outputx(1,:),cor_sz,'biased');
out_corr_y=xcorr(outputx(2,:),cor_sz,'biased');
out_corr_xy=xcorr(outputx(1,:),outputx(2,:),cor_sz,'biased');

prev_cor_x=exp_h_out(1,2)*exp_h_out(1,1).^[0:cor_sz]+exp_h_out(1,3);
prev_cor_y=exp_h_out(2,2)*exp_h_out(2,1).^[0:cor_sz]+exp_h_out(2,3);
prev_cor_xy_pos=exp_h_out(3,2)*exp_h_out(3,1).^[0:cor_sz]+exp_h_out(3,3);
prev_cor_xy_neg=exp_h_out(4,2)*exp_h_out(4,1).^abs([-cor_sz:0])+exp_h_out(4,3);

plot(0:cor_sz,out_corr_x(cor_sz+1:2*cor_sz+1)-prev_cor_x);grid;
hold on;
plot(0:cor_sz,out_corr_y(cor_sz+1:2*cor_sz+1)-prev_cor_y,'r');
plot(0:cor_sz,out_corr_xy(cor_sz+1:2*cor_sz+1)-prev_cor_xy_pos,'k');
plot(-cor_sz:0,out_corr_xy(1:cor_sz+1)-prev_cor_xy_neg,'k');

out_corr=[out_corr_x(:,cor_sz+1:2*cor_sz+1)-prev_cor_x;out_corr_y(:,cor_sz+1:2*cor_sz+1)-prev_cor_y;out_corr_xy(:,cor_sz+1:2*cor_sz+1)-prev_cor_xy_pos;fliplr(out_corr_xy(:,1:cor_sz+1)-prev_cor_xy_neg)];
[mod_l_out exp_l_out]=expo_fit_v000(ones(4,1)*[8,1000], out_corr,1,1);

allan_th=arma2allan_v000([1 -mod_h_out(1,1);1 -mod_l_out(1,1);1 0;1 -1],[mod_h_out(1,2);mod_l_out(1,2);whi_pow_out(1);rw_pow_out(1)],fs,t_axis,allan_out(1,:));
allan_th=arma2allan_v000([1 -mod_h_out(2,1);1 -mod_l_out(2,1);1 0;1 -1],[mod_h_out(2,2);mod_l_out(2,2);whi_pow_out(2);rw_pow_out(2)],fs,t_axis,allan_out(2,:));

%%%Readjust parameters for multi model
%% 1st order MW appx
[exp_h_apx]=exp_approx_v000(ones(4,1)*[2^9,15000], exp_h_out,1);
[exp_l_apx]=exp_approx_v000(ones(4,1)*[1,1000], exp_l_out,1);

mod_h_apx=exp2markov_v000(exp_h_apx(1:3,:),1);
mod_l_apx=exp2markov_v000(exp_l_apx(1:3,:),1);

allan_th=arma2allan_v000([1 -mod_h_apx(1,1);1 -mod_l_apx(1,1);1 0;1 -1],[mod_h_apx(1,2);mod_l_apx(1,2);whi_pow_out(1);rw_pow_out(1)],fs,t_axis,allan_out(1,:));
allan_th=arma2allan_v000([1 -mod_h_apx(2,1);1 -mod_l_apx(2,1);1 0;1 -1],[mod_h_apx(2,2);mod_l_apx(2,2);whi_pow_out(2);rw_pow_out(2)],fs,t_axis,allan_out(2,:));

%% 2nd method (mod_l will be approximated as second order process)
exp_h2_apx=exp_approx_v000(ones(4,1)*[2^9,15000], exp_h_out,1);
exp_l2a_apx=exp_approx_v000(ones(2,1)*[10,1000], exp_l_out([1 3],:),1);
exp_l2b_apx=[exp_l2a_apx(1,1) (exp_l2a_apx(2,2)^2)/exp_l2a_apx(1,2)];
[mod_l2bres_apx exp_l2bres_apx]=exp_approx2_v000([0,1000], exp_l_out(2,:), exp_l2b_apx,1);

[mod_h2_apx mx_a]=exp2markov_v000(exp_h2_apx(1:3,:),1);
[mod_l2a_apx mx_a]=exp2markov_v000([exp_l2a_apx(1,:);exp_l2b_apx;exp_l2a_apx(2,:)],1);

allan_th=arma2allan_v000([1 -mod_h2_apx(1,1);1 -mod_l2a_apx(1,1);1 0;1 -1],[mod_h2_apx(1,2);mod_l2a_apx(1,2);whi_pow_out(1);rw_pow_out(1)],fs,t_axis,allan_out(1,:));
allan_th=arma2allan_v000([1 -mod_h2_apx(2,1) 0;1 -mod_l2a_apx(2,1) 0; 1 -mod_l2bres_apx(1,1) -mod_l2bres_apx(1,2);1 0 0;1 -1 0],[mod_h2_apx(2,2);mod_l2a_apx(2,2);mod_l2bres_apx(1,3);whi_pow_out(2);rw_pow_out(2)],fs,t_axis,allan_out(2,:));

repo: ekfmonoslam


text:
% Li's ellipsoid specific fitting algorithm:

function v=FindFit4(S)

% Create constraint matrix C:

  C(6,6)=0;

  C(1,1)=-1; C(2,2)=-1; C(3,3)=-1;

  C(4,4)=-4; C(5,5)=-4; C(6,6)=-4;

  C(1,2)=1; C(2,1)=1;
  C(1,3)=1; C(3,1)=1;
  C(2,3)=1; C(3,2)=1;

% Solve generalized eigensystem

  S11=S(1:6, 1:6);  S12=S(1:6, 7:10);
  S22=S(7:10,7:10);

  A=S11-S12*pinv(S22)*S12';

  %[gevec, geval]=eig(A, C);

  CA=inv(C)*A;
  [gevec, geval]=eig(CA);

% Find the largest eigenvalue(the only positive eigenvalue)

  In=1;
  maxVal=geval(1,1);
  for i=2:6
    if (geval(i,i)>maxVal)
        maxVal=geval(i,i);
        In=i;
    end;
  end;

% Find the fitting

  v1=gevec(:, In);
  v2=-pinv(S22)*S12'*v1;
  v=[v1; v2];


text:
function D = driving_function_mono_sdm_ls(x0,nk,f,conf)
%DRIVING_FUNCTION_MONO_SDM_LS driving signal for a line source in SDM
%
%   Usage: D = driving_function_mono_sdm_ls(x0,nk,f,conf)
%
%   Input parameters:
%       x0          - position of the secondary sources / m [nx3]
%       nk          - position of virtual line source / m [nx3]
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   See also: driving_function_mono_sdm

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

%% ===== Checking of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);
isargmatrix(x0,nk);
isargpositivescalar(f);
isargstruct(conf);

%% ===== Configuration ==================================================
xref = conf.xref;
c = conf.c;
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation ====================================================
% Calculate the driving function in time-frequency domain

% Frequency
omega = 2*pi*f;

if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    % Ensure 2D
    x0 = x0(:,1:2);
    nk = nk(:,1:2);

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    % Reference point
    xref = repmat(xref,[size(x0,1) 1]);

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2.5D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('3D',dimension)

    % === 3-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 3D line source.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [B]=cellTrace(A)

% function [B]=cellTrace(A)
% ------------------------------------------------------------------------
% Computes the trace for each matrix contained in the cell array A, i.e.
% b=trace(a) is executed for each cell entry.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 28/11/2013
%------------------------------------------------------------------------

[B]=cellfun(@trace,A,'UniformOutput',0);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% Main SDRSAC Algorithm
% Input:
%       M, B: Input point clouds, each is a 3xN matrix
%       config: Contains all parameters required to run the algoirthm.
%               See readConfig.m for more info
% Output: a structure out containing output variables:
%         out.bestR: Best rotation matrix
%         out.bestT: Best translation vector

function [out] = SDRSAC(M, B, config)
    add_dependencies;
    ps = 0.99;

    iter=0;
    T_max = 1e10;

    % Using KDTree for quick computing of consensus size
    B_Tree = KDTreeSearcher(B');

    % Prepare sampling
    n = config.pointPerSample;  % N_sample
    maxInls = 0;
    bestR = [];
    bestT = [];

    % Start the sampling iterations
    stop = false;
    while (iter < config.maxIter && ~stop)

        idxM = randsample(size(M,2), n);
        m = M(:, idxM);
        B_to_sample = B;
        scount = 0;
        while (size(B_to_sample,2) > n*2 && scount < 2)
            scount = scount + 1;
            fprintf('Current B size: %d\n', size(B_to_sample,2));

            idxB = randsample(size(B_to_sample,2), n);
            b = B_to_sample(:, idxB);
            B_to_sample(:,idxB) = [];

            % Solve SDP
            [Rs,ts, ~, corrB] = sdpReg(m, b, config);

            if ~isempty(corrB)
                TM = Rs*M + repmat(ts, 1, size(M,2));

                % Conducting ICP
                [Ricp, Ticp] = icp(B, TM, 'Matching', 'kDtree', 'WorstRejection', 0.1, 'iter', 100);
                TMICP = Ricp*Rs*M + repmat(Ricp*ts + Ticp, 1, size(M, 2));
                inls_icp = countCorrespondences(TMICP, B_Tree, config.epsilon);

                if inls_icp > maxInls
                    maxInls = inls_icp;
                    bestR = Ricp*Rs;
                    bestT = Ricp*ts + Ticp;

                    fprintf('Best-so-far consensus size: %d\n', maxInls);%
                    fprintf('--------------');

                    % For debugging purpose:
                     %close all; plotPointClouds(B, TMICP, 'b.','r.');

                    % Compute stopping criterion
                    pI = maxInls./size(M,2);
                    T_max = log(1-ps)./log(1-pI^config.k);
                end
            end

            iter = iter+1;

            % Stopping criterion is satisfied (minimum of 5 iterations)
            if iter >= T_max && iter >=5
                 stop = true;
            end
        end

    end

    out.inls = maxInls;
    out.iter = iter;
    out.R = bestR;
    out.T = bestT;
    plotPointClouds(B, bestR*M + bestT, 'b.','r.')
end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function cliques = cliques_from_engine(engine)
% CLIQUES_FROM_ENGINE Return the cliques stored inside the inf. engine (jtree)
% cliques = cliques_from_engine(engine)

cliques = engine.cliques;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function f=iwpfbt(c,par,varargin)
%IWPFBT   Inverse Wavelet Packet Filterbank Tree
%   Usage:  f=iwpfbt(c,info);
%           f=iwpfbt(c,wt,Ls);
%
%   Input parameters:
%         c       : Coefficients stored in a cell-array.
%         info,wt : Transform parameters struct/Wavelet Filterbank tree.
%         Ls      : Length of the reconstructed signal.
%
%   Output parameters:
%         f     : Reconstructed data.
%
%   `f = iwpfbt(c,info)` reconstructs signal *f* from the coefficients *c*
%   using parameters from `info` struct. both returned by |wfbt| function.
%
%   `f = iwpfbt(c,wt,Ls)` reconstructs signal *f* from the coefficients *c*
%   using filter bank tree defined by *wt*. Plese see |wfbt| function for
%   possible formats of *wt*. The *Ls* parameter is mandatory due to the
%   ambiguity of reconstruction lengths introduced by the subsampling
%   operation and by boundary treatment methods.
%
%   Please see help for |wfbt| description of possible formats of *wt* and
%   of the additional flags.
%
%   Scaling of intermediate outputs:
%   --------------------------------
%
%   The following flags control scaling of the intermediate coefficients.
%   The intermediate coefficients are outputs of nodes which ale also
%   inputs to nodes further in the tree.
%
%      'intsqrt'
%               Each intermediate output is scaled by `1/sqrt(2)`.
%               If the filterbank in each node is orthonormal, the overall
%               undecimated transform is a tight frame.
%               This is the default.
%
%      'intnoscale'
%               No scaling of intermediate results is used.
%
%      'intscale'
%               Each intermediate output is scaled by `1/2`.
%
%   If 'intnoscale' is used, 'intscale' must have been used in |wpfbt|
%   (and vice versa) in order to obtain a perfect reconstruction.
%
%   Examples:
%   ---------
%
%   A simple example showing perfect reconstruction using the "full
%   decomposition" wavelet tree:::
%
%     f = gspi;
%     J = 7;
%     wt = {'db10',J,'full'};
%     c = wpfbt(f,wt);
%     fhat = iwpfbt(c,wt,length(f));
%     % The following should give (almost) zero
%     norm(f-fhat)
%
%   See also: wpfbt, wfbtinit
%
complainif_notenoughargs(nargin,2,'IWPFBT');

if(~iscell(c))
    error('%s: Unrecognized coefficient format.',upper(mfilename));
end

if(isstruct(par)&&isfield(par,'fname'))
   complainif_toomanyargs(nargin,2,'IWPFBT');

   if ~strcmpi(par.fname,'wpfbt')
      error('%s: Wrong func name in info struct. The info parameter was created by %s.',upper(mfilename),par.fname);
   end

   wt = wfbtinit({'dual',par.wt},par.fOrder);
   Ls = par.Ls;
   ext = par.ext;
   interscaling = par.interscaling;
   % Use the "oposite" scaling
   if strcmp(interscaling,'intscale')
      interscaling = 'intnoscale';
   elseif strcmp(interscaling,'intnoscale')
      interscaling = 'intscale';
   end

   % Determine next legal input data length.
   L = wfbtlength(Ls,wt,ext);
else
   if nargin<3
      error('%s: Too few input parameters.',upper(mfilename));
   end
   %% PARSE INPUT
   definput.keyvals.Ls=[];
   definput.import = {'fwt','wfbtcommon'};
   definput.flags.interscaling = {'intsqrt', 'intscale', 'intnoscale'};
   [flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);
   complainif_notposint(Ls,'Ls');

   ext = flags.ext;
   interscaling = flags.interscaling;
   % Initialize the wavelet tree structure
   wt = wfbtinit(par,flags.forder);

   [Lc,L]=wpfbtclength(Ls,wt,ext);

   % Do a sanity check
   if ~isequal(Lc,cellfun(@(cEl) size(cEl,1),c))
      error(['%s: The coefficients subband lengths do not comply with the'...
             ' signal length *Ls*.'],upper(mfilename));
   end

end

wtPath = fliplr(nodeBForder(0,wt));
[pOutIdxs,chOutIdxs] = treeWpBFrange(wt);
f = comp_iwpfbt(c,wt.nodes(wtPath),pOutIdxs,chOutIdxs,L,ext,interscaling);
f = postpad(f,Ls);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% Copyright (C) 2018 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym arg (@var{x})
%% @defmethodx @@sym angle (@var{x})
%% Symbolic polar angle.
%%
%% Example:
%% @example
%% @group
%% x = sym(2+3*i);
%% y = arg(x)
%%   @result{} y = (sym) atan(3/2)
%% @end group
%% @end example
%% @seealso{arg, @@sym/abs}
%% @end defmethod

function y = arg (x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('arg', x);
end

%!test
%! syms x
%! assert (isequal (angle (x), arg (x)));

repo: octsympy


text: function sarima(y,c,p,q,d,seasonal,x,startingVals,options,holdBack,sigma2)

author: bashtage

repo: mfe-toolbox


text:
function model=replaceMets(model,metabolite,replacement,verbose)
% replaceMets
%   Replaces metabolite names and annotation with replacement metabolite
%   that is already in the model. If this results in duplicate metabolites,
%   the replacement metabolite will be kept, while the S matrix is updated
%   to use the replacement metabolite instead.
%
%   model               a model structure
%   metabolite          string with name of metabolite to be replace
%   replacement         string with name of replacement metabolite
%   verbose             logical whether to print the ids of reactions that
%                       involve the replaced metabolite (opt, default
%                       false)
%
%   This function is useful when the model contains both 'oxygen' and 'o2'
%   as metabolites.
%
%   Usage: model=replaceMets(model,metabolite,replacement,verbose)

metabolite=char(metabolite);
replacement=char(replacement);

if nargin<4
    verbose=false;
end

% Find occurence of replacement metabolites. Annotation will be taken from
% first metabolite found. Metabolite ID from replacement will be used where
% possible.
repIdx = find(strcmp(replacement,model.metNames));
if isempty(repIdx)
    error('The replacement metabolite name cannot be found in the model.');
end

% Change name and information from metabolite to replacement metabolite
metIdx = find(strcmp(metabolite,model.metNames));
if isempty(metIdx)
    error('The to-be-replaced metabolite name cannot be found in the model.');
end
if verbose==true
    fprintf('\n\nThe following reactions contain the replaced metabolite as reactant:\n')
    fprintf(strjoin(model.rxns(find(model.S(metIdx,:))),'\n'))
    fprintf('\n')
end
model.metNames(metIdx) = model.metNames(repIdx(1));
if isfield(model,'metFormulas')
    model.metFormulas(metIdx) = model.metFormulas(repIdx(1));
end
if isfield(model,'metMiriams')
    model.metMiriams(metIdx) = model.metMiriams(repIdx(1));
end
if isfield(model,'metCharges')
    model.metCharges(metIdx) = model.metCharges(repIdx(1));
end
if isfield(model,'inchis')
    model.inchis(metIdx) = model.inchis(repIdx(1));
end
if isfield(model,'metSmiles')
    model.metSmiles(metIdx) = model.metSmiles(repIdx(1));
end
% Run through replacement metabolites and their compartments. If any of the
% to-be-replaced metabolites is already present (checked by
% metaboliteName[compartment], then the replacement metabolite is kept and
% the to-be-replace metabolite ID deleted.

% Build list of metaboliteName[compartment]
metCompsN =cellstr(num2str(model.metComps));
map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
metCompsN = map.values(metCompsN);
metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');

idxDelete=[];
for i = 1:length(repIdx)
    metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
    if length(metCompsNidx)>1
        for j = 2:length(metCompsNidx)
            model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
            idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
        end
    end
end

if ~isempty(idxDelete)
    model.S(idxDelete,:) =[];
    model.mets(idxDelete) = [];
    model.metNames(idxDelete) = [];
    model.metComps(idxDelete) = [];
    model.b(idxDelete) = [];
    if isfield(model,'metFormulas')
        model.metFormulas(idxDelete) = [];
    end
    if isfield(model,'unconstrained')
        model.unconstrained(idxDelete) = [];
    end
    if isfield(model,'metMiriams')
        model.metMiriams(idxDelete) = [];
    end
    if isfield(model,'metCharges')
        model.metCharges(idxDelete) = [];
    end
    if isfield(model,'inchis')
        model.inchis(idxDelete) = [];
    end
    if isfield(model,'metSmiles')
        model.metSmiles(idxDelete) = [];
    end
    if isfield(model,'metFrom')
        model.metFrom(idxDelete) = [];
    end
end

% This could now have created duplicate reactions. Contract model.
model=contractModel(model);
end

author: SysBioChalmers

repo: RAVEN


text:
function tissueModel = GIMME(model, expressionRxns, threshold, obj_frac)
% Use the GIMME algorithm (`Becker and Palsson, 2008`) to extract a context
% specific model using data. GIMME minimizes usage of low-expression
% reactions while keeping the objective (e.g., biomass) above a certain
% value. Note that this algorithm does not favor the inclusion of reactions
% not related to the objective.
%
% USAGE:
%
%    tissueModel = GIMME(model, expressionRxns, threshold)
%
% INPUTS:
%
%    model:               input model (COBRA model structure)
%    expressionRxns:      expression data, corresponding to model.rxns (see
%                         mapGeneToRxn.m)
%    threshold:           expression threshold, reactions below this are
%                         minimized
%
% OPTIONAL INPUT:
%    obj_frac:            minimum fraction of the objective(s) of model
%                         (default value - 0.9)
%
% OUTPUTS:
%    tissueModel:         extracted model
%
% `Becker and Palsson (2008). Context-specific metabolic networks are
% consistent with experiments. PLoS Comput. Biol. 4, e1000082.`
%
% .. Author: - Originally written by Becker and Palsson, adapted by S. Opdam and A. Richelle - May 2017

    if nargin < 4 || isempty(obj_frac)
        obj_frac =0.9;
    end

    objectiveCol = [find(model.c) obj_frac];

    nRxns = size(model.S,2);

    %first make model irreversible
    [modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);

    nbExpressionRxns = size(expressionRxns,1);
    if (nbExpressionRxns < nRxns)
        display('Warning: Fewer expression data inputs than reactions');
        expressionRxns(nbExpressionRxns+1:nRxns,:) = zeros(nRxns-nbExpressionRxns, size(expressionRxns,2));
    end

    nIrrevRxns = size(irrev2rev,1);
    expressionRxnsIrrev = zeros(nIrrevRxns,1);
    for i=1:nIrrevRxns
        expressionRxnsIrrev(i,1) = expressionRxns(irrev2rev(i,1),1);
    end

    nObjectives = size(objectiveCol,1);
    for i=1:nObjectives
        objectiveColIrrev(i,:) = [rev2irrev{objectiveCol(i,1),1}(1,1) objectiveCol(i,2)];
    end

    %Solve initially to get max for each objective
    for i=1:size(objectiveCol)
        %define parameters for initial solution
        modelIrrev.c=zeros(nIrrevRxns,1);
        modelIrrev.c(objectiveColIrrev(i,1),1)=1;

        %find max objective
        FBAsolution = optimizeCbModel(modelIrrev);
        if (FBAsolution.stat ~= 1)
            not_solved=1;
            error('Failed to solve initial FBA problem');
        end
        maxObjective(i)=FBAsolution.f;
    end

    model2gimme = modelIrrev;
    model2gimme.c = zeros(nIrrevRxns,1);

    for i=1:nIrrevRxns
        if (expressionRxnsIrrev(i,1) > -1)   %if not absent reaction
            if (expressionRxnsIrrev(i,1) < threshold)
                model2gimme.c(i,1) = threshold-expressionRxnsIrrev(i,1); %FIX: use expression level as weight
            end
        end
    end

    for i=1:size(objectiveColIrrev,1)
        model2gimme.lb(objectiveColIrrev(i,1),1) = objectiveColIrrev(i,2) * maxObjective(i);
    end

    gimmeSolution = optimizeCbModel(model2gimme,'min');

    if (gimmeSolution.stat ~= 1)
    %No solution for the problem
        display('Failed to solve GIMME problem');
        gimmeSolution.x = zeros(nIrrevRxns,1);
    end

    reactionActivityIrrev = zeros(nIrrevRxns,1);
    for i=1:nIrrevRxns
        if ((expressionRxnsIrrev(i,1) > threshold) || (expressionRxnsIrrev(i,1) == -1))
            reactionActivityIrrev(i,1)=1;
        elseif (gimmeSolution.x(i,1) > 0)
            reactionActivityIrrev(i,1)=2;
        end
    end

    %Translate reactionActivity to reversible model
    reactionActivity = zeros(nRxns,1);
    for i=1:nRxns
        for j=1:size(rev2irrev{i,1},2)
            if (reactionActivityIrrev(rev2irrev{i,1}(1,j)) > reactionActivity(i,1))
                reactionActivity(i,1) = reactionActivityIrrev(rev2irrev{i,1}(1,j));
            end
        end
    end

    remove = model.rxns(reactionActivity == 0);
    tissueModel = removeRxns(model,remove);
end

author: opencobra

repo: cobratoolbox


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Laplace Transform properties

% Integration in the complex frequency

x=sin(t);

L=laplace(x/t,s);
ezplot(L,[0 100])

figure
syms u
X=laplace(x,u);
R=int(X,u,s,inf);
ezplot(R,[0 100])

author: Sable

repo: mcbench-benchmarks


text:
function [lb,ub,A,b] = remove_bounds_from_Ab(A,b,lb,ub);

if size(A,1)>0
    cand_rows_lp = find(sum(A~=0,2)==1);
    if ~isempty(cand_rows_lp)
        [ii,jj,kk] = find(A(cand_rows_lp,:));
        s_pos = find(kk>0);
        s_neg = find(kk<=0);
        if ~isempty(s_pos)
            for s = 1:length(s_pos)
                ub(jj(s_pos(s)),1) = full(min(ub(jj(s_pos(s))),b(cand_rows_lp(ii(s_pos(s))))./kk(s_pos(s))));
            end
        end
        if ~isempty(s_neg)
            for s = 1:length(s_neg)
                lb(jj(s_neg(s)),1) = full(max(lb(jj(s_neg(s))),b(cand_rows_lp(ii(s_neg(s))))./kk(s_neg(s))));
            end
        end
    end
    A(cand_rows_lp,:) = [];
    b(cand_rows_lp,:) = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP