%==============================================================================
%
% function [yc,His] = lBFGS(fctn,yc,varargin)
%
% limited BFGS optimizer
%
% Input:
% ------
%   fctn           function handle
%   yc          starting guess
%   varargin    optional parameter, see below
%
% Output:
% -------
%   yc          numerical optimizer (current iterate)
%   His         iteration history
%==============================================================================

function [yc,His] = lBFGS(fctn,yc,varargin)

if nargin ==0, % help and minimal example
  help(mfilename);
  E9_MRIhead_MLIRlBFGS_NGF_mbElas;
  yc = 'endOfMinimalExample';
  His = [];
  return;
end;

% parameter initialization -----------------------------------------------
maxIter      = 10;              % maximum number of iterations
tolJ         = 1e-3;            % for stopping, objective function
tolY         = 1e-2;            %   - " -     , current value
tolG         = 1e-2;            %   - " -     , norm of gradient
solver       = regularizer;
maxLBFGS     = 5;              % maximum number of BFGS vectors
LSmaxIter    = 10;              % maximum number of line search iterations
LSreduction  = 1e-4;            % minimal reduction in line search
vecNorm      = @norm;           % norm to be used for dJ and dy
yStop        = [];              % used for stopping in multi-level framework
Jstop        = [];              %
Plots        = @(iter,para) []; % for plots;

for k=1:2:length(varargin),     % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;
objFctn = fctn([]);

if isempty(yStop), yStop  = yc; end; % yStop used for stopping only
% -- end parameter set-up   ----------------------------------------------

% some output
FAIRmessage([mfilename '(JM 2009/01/31)']);
fprintf('[ maxIter=%d / maxLBFGS=%d / tolJ=%s / tolY=%s / tolG=%s / length(yc)=%d ]\n',...
  maxIter,maxLBFGS,num2str(tolJ),num2str(tolY),num2str(tolG),length(yc));

% -- initialize  ---------------------------------------------------------
STOP = zeros(5,1);

zBFGS = [];     % memory for BFGS gradient directions
sBFGS = [];     % memory for BFGS directions
nBFGS  = 0;     % counter for number of limited BFGS directions

if isempty(Jstop),
  % evaluate objective function for stopping values and plots
  [Jstop,para] = fctn(yStop); Jstop = abs(Jstop) + (Jstop == 0);
  Plots('stop',para);
end;

% evaluate objective function for starting values and plots
[Jc,para,dJ,H0] = fctn(yc);
Plots('start',para);
iter = 0; yOld = 0*yc; Jold = Jc; y0 = yc;

hisStr    = {'iter','J','Jold-J','|\nabla J|','|dy|','LS'};
his        = zeros(maxIter+2,6);
his(1,1:3) = [-1,Jstop,Jstop-Jc];
his(2,:)   = [0,Jc,Jstop-Jc,vecNorm(dJ),vecNorm(yc-yStop),0];

% some output
fprintf('%4s %-12s %-12s %-12s %-12s %4s\n%s\n',...
  hisStr{:},char(ones(1,64)*'-'));
dispHis = @(var) ...
  fprintf('%4d %-12.4e %-12.3e %-12.3e %-12.3e %4d\n',var);
dispHis(his(1,:));
% -- end initialization   ------------------------------------------------

% ==============================================================================
% MAIN LOOP
% ==============================================================================
while 1,
  % check stopping rules
  STOP(1) = abs(Jold-Jc)  <= tolJ*(1+abs(Jstop));
  STOP(2) = (iter>0) && (norm(yc-yOld) <= tolY*(1+norm(y0)));
  STOP(3) = norm(dJ)      <= tolG*(1+abs(Jstop));
  STOP(4) = norm(dJ)      <= 1e3*eps;
  STOP(5) = (iter >= maxIter);
  if all(STOP(1:3)) || any(STOP(4:5)), break;  end;

  iter      = iter + 1;

  % update at most maxLBFGS BFGS directions
  if iter > 1,
    zz = (dJ - dJold)';
    ss =  yc - yOld;
    if zz'*ss > 0,
      start = 2-(nBFGS<maxLBFGS);
      zBFGS = [zBFGS(:,start:end),zz];
      sBFGS = [sBFGS(:,start:end),ss];
      nBFGS = min(maxLBFGS,nBFGS+1);
    else
      warning('   y''*s < 0, skip update');
    end;
  end;

  % compute search direction using recursive and limited BFGS
  dy = bfgsrec(solver,nBFGS,sBFGS,zBFGS,H0,-dJ');
  if iter == 1,
    if strcmp(objFctn,'NPIR'),
      m     = para.m;
      omega = para.omega;
      dim   = length(para.omega)/2;
      maxdy = norm(sqrt(sum(reshape(center(dy,m),[],dim).^2,2)),'inf');
      fac   = 0.5*max((omega(2:2:end)-omega(1:2:end)));
      if maxdy>fac,
        dy = (fac/maxdy)*dy;
        %H0 = fac*H0;%speye(length(dy),length(dy));
      end;
    else
      H0 = norm(dJ) * speye(length(dy),length(dy));
    end;
  end;
  % perform Armijo line-search
  [t,yt,LSiter] = Armijo(fctn,yc,dy,Jc,dJ,...
        'LSmaxIter',LSmaxIter,'LSreduction',LSreduction);
  if (t == 0),
    break;
  end; % break if line-search fails

  % update variables
  yOld = yc; Jold = Jc; dJold = dJ; yc = yt;
  [Jc,para,dJ] = fctn(yc);  % evaluate objective function

  % some output
  his(iter+2,:) = [iter,Jc,Jold-Jc,vecNorm(dJ),vecNorm(yc-yOld),LSiter];
  dispHis(his(iter+1,:));
  para.normdY = vecNorm(yc - yOld);
  Plots(iter,para);

end;%while; % end of iteration loop
% ==============================================================================

% clean up
His.str = hisStr;
His.his = his(1:iter+2,:);
fprintf('STOPPING:\n');
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(1),...
  '(Jold-Jc)',(Jold-Jc),'tolJ*(1+|Jstop|)',tolJ*(1+abs(Jstop)));
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(2),...
  '|yc-yOld|',norm(yc-yOld),'tolY*(1+norm(yc)) ',tolY*(1+norm(yc)));
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(3),...
  '|dJ|',norm(dJ),'tolG*(1+abs(Jstop)',tolG*(1+abs(Jstop)));
fprintf('%d[ %-10s=%16.8e <= %-25s=%16.8e]\n',STOP(4),...
  'norm(dJ)',norm(dJ),'eps',1e3*eps);
fprintf('%d[ %-10s=  %-14d >= %-25s=  %-14d]\n',STOP(5),...
  'iter',iter,'maxIter',maxIter);

FAIRmessage([mfilename,' : done !']);

%------------------------------------------------------------------------------

function[d] = bfgsrec(solver,n,S,Z,H,d)

maxIterCG = 10; tolCG = 1e-2;
if isempty(solver)
  if isstruct(H),
    if isfield(H,'solver'),
      solver = H.solver;
    elseif isfield(H,'d2S') && isfield(H.d2S,'solver'),
      solver = H.d2S.solver;
    else
      error('solver has not been defined')
    end;

  end
end

if n == 0,
  switch solver,
    case {'MG-elastic'}
        d = MGsolver(d,H);
    case {'pcg'}
      L   = tril(H); % Symmetric Gauss Seidel Preconditioning,
      D   = diag(H); % L is lower, D is diagonal, U = L'
      SGS = @(x) L\(D.*(L'\x));
      d   = pcg(H,d,tolCG,maxIterCG,SGS);
    case 'PCG-curvature',
      M     = @(x) H.d2D.M;
      Afctn = @(x) M(x) + H.d2S.d2S(x,H.omega,H.m);
      % preconditioner
      Ddiag = H.d2D.M;
      D = Ddiag  +  H.d2S.diag(H.omega,H.m);
      PC = @(x) D.\x; % Jacobi preconditioner
      [d,flag,relres,iter] = pcg(Afctn,d,tolCG,maxIterCG,PC);

    otherwise,
      if isnumeric(H),
        % if H is a matrix, solve the linear system using MATLAB's backslash
        d = H\d;
      else
        error('nyi - solver %s', solver)
      end;
  end;
else
  alpha = (S(:,n)'*d)/(Z(:,n)'*S(:,n));
  d     = d - alpha*Z(:,n);
  d     = bfgsrec(solver,n-1,S,Z,H,d);
  d     = d + (alpha - (Z(:,n)'*d)/(Z(:,n)'*S(:,n)))*S(:,n);
end;
%==============================================================================


text:
function plot_boundary_orange(V,F)
%PLOT_BOUNDARY_ORANGE  plots the boundary vertices of a triangle mesh in
%orange, and its interior in blue.
%
% plot_boundary_orange(V,F);
%
% Input:
%  V,F  mesh whose boundary is to be plotted

%Compute a function u that is 0 on all interior vertices and 1 on all
%boundary vertices.
b = unique(outline(F)); %find all boundary vertices
u = zeros(size(V,1),1); %construct a constant zero function
u(b) = 1; %set all boundary vertices to 1

%Plot the function u with the right colors.
t = tsurf(F,V, 'CData',u);
shading interp;
axis equal;
axis off;
cm = flipud(cbrewer('RdYlBu', 500));
colormap(cm(100:450,:));
light('Position',[-1.5 1 1],'Style','local');
lights = camlight;
set(t, 'FaceLighting','gouraud', 'FaceColor','interp');
set(t, 'DiffuseStrength',0.5, 'SpecularStrength',0.2, 'AmbientStrength',0.3);
camproj('perspective');
add_shadow([t],lights);

end

author: odedstein

repo: sgi-introduction-course


text:
function varargout = min2(varargin)
%MIN2   Global minimum of a DISKFUN.
%   M = MIN2(F) returns the global minimum of F over its domain.
%
%   [M, LOC] = MIN2(F) returns the global minimum in M and its location in
%   LOC.
%
%  This command may be faster if the OPTIMIZATION TOOLBOX is installed.
%
% See also DISKFUN/MAX2, DISKFUN/MINANDMAX2.

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

f = varargin{1};

if ( isempty(f) )
 varargout= [];
 return
end

% Give cdr to separableApprox.
f = cart2pol(f, 'cdr');
[out{1:nargout}] = min2@separableApprox(f);
if numel(out) == 1
    varargout = out;
else
    % Return loc in Cartesian coords.
    [val, loc] = out{:};
    t = loc(1);
    loc(1) = loc(2).*cos(t);
    loc(2) = loc(2).*sin(t);
    varargout = {val, loc};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef MeshSymmetrizerTest < testShowingError

    properties (Access = protected)
        tol = 1e-14;
    end

    properties (Access = private)
        mesh
        symmetricMesh
    end

    properties (Access = private)
        alpha
        nx
        ny
        x1Max
        x1Min
        x2Max
        x2Min
    end

    methods (Access = public)

        function obj = MeshSymmetrizerTest()
            obj.init();
            obj.createMesh();
            obj.createSymmetricMesh();
            obj.plotMeshes();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.nx = 20;
            obj.ny = 30;
            obj.x1Max = 3;
            obj.x1Min = 0;
            obj.x2Max = 2;
            obj.x2Min = 0;
            obj.alpha = 0.3;
        end

        function createMesh(obj)
            [x1,x2] = obj.createCoord();
            x3 = zeros(size(x2));
            [F,V] = mesh2tri(x1,x2,x3,'x');
            s.coord  = V(:,1:2);
            s.connec = F;
            m = Mesh(s);
            obj.mesh = m;
        end

        function [x1,x2] = createCoord(obj)
            [x1,x2] = obj.createHorizontalVerticalCoords();
            [x1,x2] = obj.rotateCoord(x1,x2);
        end

        function [x1v,x2v] = createHorizontalVerticalCoords(obj)
            x1 = linspace(obj.x1Min,obj.x1Max,obj.nx);
            x2 = linspace(obj.x2Min,obj.x2Max,obj.ny);
            [x1v,x2v] = meshgrid(x1,x2);
        end

        function [x1r,x2r] = rotateCoord(obj,x1,x2)
            ca = cos(obj.alpha);
            sa = sin(obj.alpha);
            x1r = ca*x1 - sa*x2;
            x2r = sa*x1 + ca*x2;
        end

        function createSymmetricMesh(obj)
            s.mesh = obj.mesh;
            s.symmetricLine.vector = [cos(obj.alpha);sin(obj.alpha)];
            s.symmetricLine.point = [0;0];
            mS = Symmetrizer(s);
            m = mS.computeSymmetricMesh();
            obj.symmetricMesh = m;
        end

        function plotMeshes(obj)
            figure()
            obj.mesh.plot();
            figure()
            obj.symmetricMesh.plot();
        end

    end

    methods (Access = protected)

        function computeError(obj)
            obj.error = 0;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ ] = display_graph(x_category, y_category, algorithm_list, w_list, info_list)
% SHow graphs of optimizations
%
% Inputs:
%       x_category          "numofgrad" or "iter" or "epoch" or "grad_calc_count"
%       y_category          "cost" or "optimality_gap" or "gnorm"
%       algorithms_list     algorithms to be evaluated
%       w_list              solution produced by each algorithm
%       info_list           statistics produced by each algorithm
%
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Oct. 23, 2016
% Modified by H.Kasai on Nov. 02, 2016

    % for plotting
    linetype = {'r','b','c','g','m','y','r--','b--','c--','g--','m--','y--','r:','b:','c:','g:','m:','y:','r.','b.','c.','g.','m.','y.'};
    fontsize = 16;
    markersize = 5;
    linewidth = 2;

    % initialize
    legend_str = cell(1);
    alg_num = 0;

    % plot
    figure;
    for alg_idx=1:length(algorithm_list)
        if ~isempty(info_list{alg_idx})
            alg_num = alg_num + 1;

            if strcmp(x_category, 'numofgrad')
                x_plot_data = info_list{alg_idx}.grad_calc_count;
            elseif strcmp(x_category, 'iter')
                x_plot_data = info_list{alg_idx}.iter;
            elseif strcmp(x_category, 'epoch')
                x_plot_data = info_list{alg_idx}.epoch;
            elseif strcmp(x_category, 'grad_calc_count')
                x_plot_data = info_list{alg_idx}.grad_calc_count;
            elseif strcmp(x_category, 'time')
                x_plot_data = info_list{alg_idx}.time;
            else
            end

            if strcmp(y_category, 'cost')
                y_plot_data = info_list{alg_idx}.cost;
            elseif strcmp(y_category, 'optimality_gap')
                y_plot_data = info_list{alg_idx}.optgap;
            elseif strcmp(y_category, 'gnorm')
                y_plot_data = info_list{alg_idx}.gnorm;
            elseif strcmp(y_category, 'K')
                y_plot_data = info_list{alg_idx}.K;
            elseif strcmp(y_category, 'orth')
                y_plot_data = info_list{alg_idx}.orth;
            elseif strcmp(y_category, 'symmetry')
                y_plot_data = [info_list{alg_idx}.symmetry];
            elseif strcmp(y_category, 'clustering_acc')
                y_plot_data = [info_list{alg_idx}.clustering_acc.acc];
            elseif strcmp(y_category, 'clustering_nmi')
                y_plot_data = [info_list{alg_idx}.clustering_acc.nmi];
            elseif strcmp(y_category, 'clustering_purity')
                y_plot_data = [info_list{alg_idx}.clustering_acc.purity];
            end

            semilogy(x_plot_data, y_plot_data, linetype{alg_num}, 'MarkerSize', markersize, 'Linewidth', linewidth); hold on;
            %plot(x_plot_data, y_plot_data, linetype{alg_num}, 'MarkerSize', markersize, 'Linewidth', linewidth); hold on;

            legend_str{alg_num} = algorithm_list{alg_idx};
        else
            %
        end
    end
    hold off;

    % X label
    if strcmp(x_category, 'numofgrad')
        xlabel('Number of gradient evaluations', 'FontSize', fontsize);
    elseif strcmp(x_category, 'iter')
        xlabel('Iteration', 'FontSize', fontsize);
    elseif strcmp(x_category, 'epoch')
        xlabel('Epoch', 'FontSize', fontsize);
    elseif strcmp(x_category, 'grad_calc_count')
        xlabel('# of grad', 'FontSize', fontsize);
    elseif strcmp(x_category, 'time')
        xlabel('Time', 'FontSize', fontsize);
    end

    % Y label
    if strcmp(y_category, 'cost')
        ylabel('Cost', 'FontSize', fontsize);
    elseif strcmp(y_category, 'optimality_gap')
        ylabel('Optimality gap', 'FontSize', fontsize);
    elseif strcmp(y_category, 'gnorm')
        ylabel('Norm of gradient', 'FontSize', fontsize);
    elseif strcmp(y_category, 'K')
        ylabel('Batch size', 'FontSize', fontsize);
    elseif strcmp(y_category, 'orth')
        ylabel('Orthogonality', 'FontSize', fontsize);
    elseif strcmp(y_category, 'symmetry')
        ylabel('norm(W-Wt)', 'FontSize', fontsize);
    elseif strcmp(y_category, 'clustering_acc')
        ylabel('ACC', 'FontSize', fontsize);
    elseif strcmp(y_category, 'clustering_nmi')
        ylabel('NMI', 'FontSize', fontsize);
    elseif strcmp(y_category, 'clustering_purity')
        ylabel('Purity', 'FontSize', fontsize);
    end
    legend(legend_str);
    set(gca, 'FontSize', fontsize);
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function []=uw_sb_unwrap_space_time(day,ifgday_ix,unwrap_method,time_win,la_flag,bperp,n_trial_wraps,prefilt_win,scf_flag,temp,n_temp_wraps,max_bperp_for_temp_est)
%UW_SB_UNWRAP_SPACE_TIME smooth and unwrap phase diffs between neighboring data points in time
%
%   Andy Hooper, June 2007
%
%   ======================================================================
%   07/2008 AH: Allow isolated images, not included in any interferogram
%   01/2010 AH: Allow for non-positive definite inverse var/covar matrix
%   01/2012 AH: New option la_flag that estimates LA error for each arc
%   02/2012 AH: New method 3D_NEW
%   02/2012 AH: New method 3D_FULL
%   03/2012 AH: version updated for subsets of full network
%   10/2012 AH: Bug corrected in 3D_FULL
%   11/2012 AH: Add temperature option and scf_flag
%   09/2013 DB: Fix compatibility issue with matlab 2009
%   10/2013 DB: Fix introduced bug on them variable
%   09/2014 DB: Fix the order in which isolated dates are removed. Start
%               from big to small.
%   02/2014 AH: Add predefined ph_uw option
%   ======================================================================
%
%
disp('Unwrapping in time-space...')

uw=load('uw_grid');
ui=load('uw_interp');

n_ifg=uw.n_ifg;
n_ps=uw.n_ps;
nzix=uw.nzix;
ij=uw.ij;

if isempty(uw.ph_uw_predef)
    predef_flag='n';
else
    predef_flag='y';
end

n_image=size(day,1);
master_ix=find(day==0);
[nrow,ncol]=size(ui.Z);

day_pos_ix=find(day>0);
[tempdummy,I]=min(day(day_pos_ix));
dph_space=((uw.ph(ui.edgs(:,3),:).*conj(uw.ph(ui.edgs(:,2),:))));
if predef_flag=='y'
    dph_space_uw=uw.ph_uw_predef(ui.edgs(:,3),:)-uw.ph_uw_predef(ui.edgs(:,2),:);
    predef_ix=~isnan(dph_space_uw);
    dph_space_uw=dph_space_uw(predef_ix);
else
    predef_ix=[];
end
clear uw tempdummy

dph_space=dph_space./abs(dph_space);

ifreq_ij=[];
jfreq_ij=[];

G=zeros(n_ifg,n_image);
for i=1:n_ifg
    G(i,ifgday_ix(i,1))=-1;
    G(i,ifgday_ix(i,2))=1;
end

nzc_ix=sum(abs(G))~=0; % non-zero column index
day=day(nzc_ix);
n_image=size(day,1);
G=G(:,nzc_ix);
zc_ix=find(nzc_ix==0);
zc_ix = sort(zc_ix,'descend');              % fix to make sure the order goes from largest to smallest
for i=1:length(zc_ix)
    ifgday_ix(ifgday_ix>zc_ix(i))=ifgday_ix(ifgday_ix>zc_ix(i))-1;
end
n=size(G,2);

if ~isempty(temp)
    temp_flag='y';
else
    temp_flag='n';
end

if strcmpi(temp_flag,'y')
    fprintf('   Estimating temperature correlation (elapsed time=%ds)\n',round(toc))
    ix=abs(bperp)<max_bperp_for_temp_est;
    temp_sub=temp(ix);
    temp_range=max(temp)-min(temp);
    temp_range_sub=max(temp_sub)-min(temp_sub);
    dph_sub=dph_space(:,ix); % only ifgs using ith image
    n_temp_wraps=n_temp_wraps*(temp_range_sub/temp_range);

    trial_mult=[-ceil(8*n_temp_wraps):ceil(8*n_temp_wraps)];
    n_trials=length(trial_mult);
    trial_phase=temp_sub/temp_range_sub*pi/4;
    trial_phase_mat=exp(-j*trial_phase*trial_mult);
    Kt=zeros(ui.n_edge,1,'single');
    coh=zeros(ui.n_edge,1,'single');
    for i=1:ui.n_edge
        cpxphase=dph_sub(i,:).';
        cpxphase_mat=repmat(cpxphase,1,n_trials);
        phaser=trial_phase_mat.*cpxphase_mat;
        phaser_sum=sum(phaser);
        coh_trial=abs(phaser_sum)/sum(abs(cpxphase));
        [coh_max,coh_max_ix]=max(coh_trial);
        falling_ix=find(diff(coh_trial(1:coh_max_ix))<0); % segemnts prior to peak where falling
        if ~isempty(falling_ix)
            peak_start_ix=falling_ix(end)+1;
        else
            peak_start_ix=1;
        end
        rising_ix=find(diff(coh_trial(coh_max_ix:end))>0); % segemnts after peak where rising
        if ~isempty(rising_ix)
            peak_end_ix=rising_ix(1)+coh_max_ix-1;
        else
            peak_end_ix=n_trials;
        end
        coh_trial(peak_start_ix:peak_end_ix)=0;
        if coh_max-max(coh_trial)>0.1 % diff between peak and next peak at least 0.1
            K0=pi/4/temp_range_sub*trial_mult(coh_max_ix);
            resphase=cpxphase.*exp(-1i*(K0*temp_sub)); % subtract approximate fit
            offset_phase=sum(resphase);
            resphase=angle(resphase*conj(offset_phase)); % subtract offset, take angle (unweighted)
            weighting=abs(cpxphase);
            mopt=double(weighting.*temp_sub)\double(weighting.*resphase);
            Kt(i)=K0+mopt;
            phase_residual=cpxphase.*exp(-1i*(Kt(i)*temp_sub));
            mean_phase_residual=sum(phase_residual);
            coh(i)=abs(mean_phase_residual)/sum(abs(phase_residual));
        end
    end

    clear cpxphase_mat trial_phase_mat phaser
    Kt(coh<0.31)=0; % not to be trusted;
    dph_space=dph_space.*exp(-1i*Kt*temp');
    if predef_flag=='y'
        dph_temp=Kt*temp';
        dph_space_uw=dph_space_uw-dph_temp(predef_ix);
        clear dph_temp
    end
    dph_sub=dph_sub.*exp(-1i*Kt*temp_sub');

end

if strcmpi(la_flag,'y')

    fprintf('   Estimating look angle error (elapsed time=%ds)\n',round(toc))

    bperp_range=max(bperp)-min(bperp);
    ix=find(abs(diff(ifgday_ix,[],2))==1);

    if length(ix)>=length(day)-1 % if almost full chain is present
        fprintf('   using sequential daisy chain of interferograms\n')
        dph_sub=dph_space(:,ix); % only ifgs using ith image
        bperp_sub=bperp(ix);
        bperp_range_sub=max(bperp_sub)-min(bperp_sub);
        n_trial_wraps=n_trial_wraps*(bperp_range_sub/bperp_range);
    else
        ifgs_per_image=sum(abs(G));
        [max_ifgs_per_image,max_ix]=max(ifgs_per_image);
        if max_ifgs_per_image>=length(day)-2; % can create chain from (almost) complete single master series
            fprintf('   Using sequential daisy chain of interferograms\n')
            ix=G(:,max_ix)~=0;
            gsub=G(ix,max_ix);
            sign_ix=-sign(single(gsub'));
            dph_sub=dph_space(:,ix); % only ifgs using chosen master image
            bperp_sub=[bperp(ix)];
            bperp_sub(sign_ix==-1)=-bperp_sub(sign_ix==-1);
            bperp_sub=[bperp_sub;0]; % add master
            sign_ix=repmat(sign_ix,ui.n_edge,1);
            dph_sub(sign_ix==-1)=conj(dph_sub(sign_ix==-1)); % flip sign if necessary to make ith image master
            dph_sub=[dph_sub,mean(abs(dph_sub),2)]; % add zero phase master
            slave_ix=sum(ifgday_ix(ix,:),2)-max_ix;
            day_sub=day([slave_ix;max_ix]); % extract days for subset
            [day_sub,sort_ix]=sort(day_sub); % sort ascending day
            dph_sub=dph_sub(:,sort_ix); % sort ascending day
            bperp_sub=bperp_sub(sort_ix);
            bperp_sub=diff(bperp_sub);
            bperp_range_sub=max(bperp_sub)-min(bperp_sub);
            n_trial_wraps=n_trial_wraps*(bperp_range_sub/bperp_range);
            n_sub=length(day_sub);
            dph_sub=dph_sub(:,[2:end]).*conj(dph_sub(:,1:end-1)); % sequential dph, to reduce influence of defo
            dph_sub=dph_sub./abs(dph_sub); % normalise
        else % just use all available
            dph_sub=dph_space;
            bperp_sub=bperp;
            bperp_range_sub=bperp_range;
        end
    end

    trial_mult=[-ceil(8*n_trial_wraps):ceil(8*n_trial_wraps)];
    n_trials=length(trial_mult);
    trial_phase=bperp_sub/bperp_range_sub*pi/4;
    trial_phase_mat=exp(-j*trial_phase*trial_mult);
    K=zeros(ui.n_edge,1,'single');
    coh=zeros(ui.n_edge,1,'single');
    for i=1:ui.n_edge
        cpxphase=dph_sub(i,:).';
        cpxphase_mat=repmat(cpxphase,1,n_trials);
        phaser=trial_phase_mat.*cpxphase_mat;
        phaser_sum=sum(phaser);
        coh_trial=abs(phaser_sum)/sum(abs(cpxphase));
        [coh_max,coh_max_ix]=max(coh_trial);
        falling_ix=find(diff(coh_trial(1:coh_max_ix))<0); % segemnts prior to peak where falling
        if ~isempty(falling_ix)
            peak_start_ix=falling_ix(end)+1;
        else
            peak_start_ix=1;
        end
        rising_ix=find(diff(coh_trial(coh_max_ix:end))>0); % segemnts after peak where rising
        if ~isempty(rising_ix)
            peak_end_ix=rising_ix(1)+coh_max_ix-1;
        else
            peak_end_ix=n_trials;
        end
        coh_trial(peak_start_ix:peak_end_ix)=0;
        if coh_max-max(coh_trial)>0.1 % diff between peak and next peak at least 0.1
            K0=pi/4/bperp_range_sub*trial_mult(coh_max_ix);
            resphase=cpxphase.*exp(-1i*(K0*bperp_sub)); % subtract approximate fit
            offset_phase=sum(resphase);
            resphase=angle(resphase*conj(offset_phase)); % subtract offset, take angle (unweighted)
            weighting=abs(cpxphase);
            mopt=double(weighting.*bperp_sub)\double(weighting.*resphase);
            K(i)=K0+mopt;
            phase_residual=cpxphase.*exp(-1i*(K(i)*bperp_sub));
            mean_phase_residual=sum(phase_residual);
            coh(i)=abs(mean_phase_residual)/sum(abs(phase_residual));
        end
    end

    clear cpxphase_mat trial_phase_mat phaser dph_sub
    K(coh<0.31)=0; % not to be trusted;
    if strcmpi(temp_flag,'y')
        dph_space(K==0,:)=dph_space(K==0,:).*exp(1i*Kt(K==0)*temp'); % add back temp correction if not able to estimate DEM reliably
        Kt(K==0)=0;
        K(Kt==0)=0;
    end
    dph_space=dph_space.*exp(-1i*K*bperp');
    if predef_flag=='y'
        dph_scla=K*bperp';
        dph_space_uw=dph_space_uw-dph_scla(predef_ix);
        clear dph_scla
    end

end

spread=sparse(zeros(ui.n_edge,n_ifg));

if strcmpi(unwrap_method,'2D')
    dph_space_uw=angle(dph_space);
    if strcmpi(la_flag,'y')
        dph_space_uw=dph_space_uw+K*bperp';   % equal to dph_space + integer cycles
    end
    if strcmpi(temp_flag,'y')
        dph_space_uw=dph_space_uw+Kt*temp';   % equal to dph_space + integer cycles
    end
    dph_noise=[];
    save('uw_space_time','dph_space_uw','spread','dph_noise');
elseif strcmpi(unwrap_method,'3D_NO_DEF')
    dph_noise=angle(dph_space);
    dph_space_uw=angle(dph_space);
    if strcmpi(la_flag,'y')
        dph_space_uw=dph_space_uw+K*bperp';   % equal to dph_space + integer cycles
    end
    if strcmpi(temp_flag,'y')
        dph_space_uw=dph_space_uw+Kt*temp';   % equal to dph_space + integer cycles
    end
    save('uw_space_time','dph_space_uw','dph_noise','spread');
else
    fprintf('   Smoothing in time (elapsed time=%ds)\n',round(toc))

    if strcmpi(unwrap_method,'3D_FULL')

       dph_smooth_ifg=nan(size(dph_space),'single');

       for i=1:n_image
         ix=G(:,i)~=0;
         if sum(ix)>=n_image-2 % allow only 1 missing ifg for each image as master
           gsub=G(ix,i);
           dph_sub=dph_space(:,ix); % only ifgs using ith image
           sign_ix=repmat(-sign(single(gsub')),ui.n_edge,1);
           dph_sub(sign_ix==-1)=conj(dph_sub(sign_ix==-1)); % flip sign if necessary to make ith image master
           slave_ix=sum(ifgday_ix(ix,:),2)-i;
           day_sub=day(slave_ix); % extract days for subset
           [day_sub,sort_ix]=sort(day_sub); % sort ascending day
           dph_sub=dph_sub(:,sort_ix); % sort ascending day
           dph_sub_angle=angle(dph_sub);
           n_sub=length(day_sub);
           dph_smooth=zeros(ui.n_edge,n_sub,'single');
           for i1=1:n_sub
               time_diff=(day_sub(i1)-day_sub)';
               weight_factor=exp(-(time_diff.^2)/2/time_win^2);
               weight_factor=weight_factor/sum(weight_factor);

               dph_mean=sum(dph_sub.*repmat(weight_factor,ui.n_edge,1),2);
               %dph_mean_adj=angle(dph_sub.*repmat(conj(dph_mean),1,n_sub)); % subtract weighted mean
               dph_mean_adj=mod(dph_sub_angle-repmat(angle(dph_mean),1,n_sub)+pi,2*pi)-pi;
               GG=[ones(n_sub,1),time_diff'];
               if size(GG,1)>1
                   m=lscov(GG,double(dph_mean_adj)',weight_factor);
               else
                   m=zeros(size(GG,1),ui.n_edge);
               end
               %dph_mean_adj=mod(dph_mean_adj-(GG*m)'+pi,2*pi)-pi; % subtract first estimate
               %m2=lscov(GG,double(dph_mean_adj)',weight_factor);
               %dph_smooth(:,i1)=dph_mean.*exp(1i*(m(1,:)'+m2(1,:)')); % add back weighted mean
               dph_smooth(:,i1)=dph_mean.*exp(1i*(m(1,:)')); % add back weighted mean
            end
            dph_smooth_sub=cumsum([angle(dph_smooth(:,1)),angle(dph_smooth(:,2:end).*conj(dph_smooth(:,1:end-1)))],2);
            close_master_ix=find(slave_ix-i>0);
            if isempty(close_master_ix)
                close_master_ix=n_sub;
            else
                close_master_ix=close_master_ix(1);
                if close_master_ix>1
                    close_master_ix=[close_master_ix-1;close_master_ix];
                end
            end
            dph_close_master=mean(dph_smooth_sub(:,close_master_ix),2);
            dph_smooth_sub=dph_smooth_sub-repmat(dph_close_master-angle(exp(j*dph_close_master)),1,n_sub);
            dph_smooth_sub=dph_smooth_sub.*sign_ix;
            already_sub_ix=find(~isnan(dph_smooth_ifg(1,ix))); % already unwrapped index
            ix=find(ix);
            already_ix=ix(already_sub_ix);
            std_noise1=std(angle(dph_space(:,already_ix).*exp(-1i*dph_smooth_ifg(:,already_ix))));
            std_noise2=std(angle(dph_space(:,already_ix).*exp(-1i*dph_smooth_sub(:,already_sub_ix))));
            keep_ix=true(n_sub,1);
            keep_ix(already_sub_ix(std_noise1<std_noise2))=false; % keep least noisy
            dph_smooth_ifg(:,ix(keep_ix))=dph_smooth_sub(:,keep_ix);
          end
        end

        dph_noise=angle(dph_space.*exp(-1i*dph_smooth_ifg));
        dph_noise(std(dph_noise,0,2)>1.2,:)=nan;

    else
        x=(day-day(1))*(n-1)/(day(end)-day(1)); % use dates for smoothing
        if predef_flag=='y'
            n_dph=size(dph_space,1);
            dph_space_angle=double(angle(dph_space));
            dph_space_angle(predef_ix)=dph_space_uw;
            dph_space_series=zeros(n,n_dph);
            for i=1:n_dph
                 W=predef_ix(i,:)+0.01; % give more weight to the predfined unwrapped
                 dph_space_series(2:end,i)=lscov(double(G(:,2:end)),dph_space_angle(i,:)',W);
            end
        else
            dph_space_series=[zeros(1,ui.n_edge);double(G(:,2:end))\double(angle(dph_space))'];
        end

        dph_smooth_series=zeros(size(G,2),ui.n_edge,'single');

        for i1=1:n
            time_diff_sq=(day(i1)-day).^2;
            weight_factor=exp(-time_diff_sq/2/time_win^2);
            weight_factor=weight_factor/sum(weight_factor);
            dph_smooth_series(i1,:)=sum(dph_space_series.*repmat(weight_factor,1,ui.n_edge));
        end

        dph_smooth_ifg=(G*dph_smooth_series)';
        dph_noise=angle(dph_space.*exp(-1i*dph_smooth_ifg));

        if strcmpi(unwrap_method,'3D_SMALL_DEF')|...
           strcmpi(unwrap_method,'3D_QUICK')
          not_small_ix=find(std(dph_noise,0,2)>1.3)';
          fprintf('   %d edges with high std dev in time (elapsed time=%ds)\n',length(not_small_ix),round(toc))
          dph_noise(not_small_ix,:)=nan;
        else % 3D
          uw=load('uw_grid');
          ph_noise=angle(uw.ph.*conj(uw.ph_lowpass));
          clear uw
          dph_noise_sf=((ph_noise(ui.edgs(:,3),:)-(ph_noise(ui.edgs(:,2),:))));
          m_minmax=repmat([-pi,pi],5,1).*repmat([0.5;0.25;1;0.25;1],1,2);
          anneal_opts=[1;15;0;0;0;0;0];
          covm=cov((dph_noise_sf)); % estimate of covariance
          [W,P]=chol(inv(covm)); % weighting matrix
          if P~=0
              W=diag(1./sqrt(diag(covm)));
          end
          not_small_ix=find(std(dph_noise,0,2)>1)';
          fprintf('   Performing complex smoothing on %d edges (elapsed time=%ds)\n',length(not_small_ix),round(toc))

          n_proc=0;
          for i=not_small_ix
            dph=angle(dph_space(i,:))';
            dph_smooth_series(:,i)=uw_sb_smooth_unwrap(m_minmax,anneal_opts,G,W,dph,x);

            n_proc=n_proc+1;
            if round(n_proc/1000)==n_proc/1000
                save('uw_unwrap_time','G','dph_space','dph_smooth_series');
                fprintf('%d edges of %d reprocessed (elapsed time=%ds)\n',n_proc,length(not_small_ix),round(toc))
            end
          end
          dph_smooth_ifg=(G*dph_smooth_series)';
          dph_noise=angle(dph_space.*exp(-1i*dph_smooth_ifg));
        end
    end
    clear dph_space
    dph_space_uw=dph_smooth_ifg+dph_noise;
    clear dph_smooth_ifg

    if strcmpi(la_flag,'y')
        dph_space_uw=dph_space_uw+K*bperp';   % equal to dph_space + integer cycles
    end
    if strcmpi(temp_flag,'y')
        dph_space_uw=dph_space_uw+Kt*temp';   % equal to dph_space + integer cycles
    end

    if strcmpi(scf_flag,'y')

        fprintf('   Calculating local phase gradients (elapsed time=%ds)\n',round(toc))
        ifreq_ij=nan(n_ps,n_ifg,'single');
        jfreq_ij=nan(n_ps,n_ifg,'single');
        ifgw=zeros(nrow,ncol);
        uw=load('uw_grid');
        for i=1:n_ifg
            ifgw(nzix)=uw.ph(:,i);
            [ifreq,jfreq,grad_ij,Hmag]=gradient_filt(ifgw,prefilt_win);
            ix=~isnan(ifreq)&Hmag./((abs(ifreq))+1)>3;
            if sum(ix(:))>2
                ifreq_ij(:,i)=griddata(grad_ij(ix,2),grad_ij(ix,1),ifreq(ix),ij(:,2),ij(:,1),'linear');
            end
            ix=~isnan(jfreq)&Hmag./((abs(jfreq))+1)>3;
            if sum(ix(:))>2
                jfreq_ij(:,i)=griddata(grad_ij(ix,2),grad_ij(ix,1),jfreq(ix),ij(:,2),ij(:,1),'linear');
            end
        end
        clear uw

        spread2=zeros(size(spread),'single');
        dph_smooth_uw2=nan(ui.n_edge,n_ifg,'single');

        fprintf('   Smoothing using local phase gradients (elapsed time=%ds)\n',round(toc))
        for i=1:ui.n_edge
            nodes_ix=ui.edgs(i,[2:3]);
            ifreq_edge=mean(ifreq_ij(nodes_ix,:));
            jfreq_edge=mean(jfreq_ij(nodes_ix,:));
            diff_i=diff(ij(nodes_ix,1));
            diff_j=diff(ij(nodes_ix,2));
            dph_smooth_uw2(i,:)=diff_i*ifreq_edge+diff_j*jfreq_edge;
%           spread2(i,:)=diff_i*diff(ifreq_ij(nodes_ix,:))+diff_j*diff(jfreq_ij(nodes_ix,:));
            spread2(i,:)=diff(ifreq_ij(nodes_ix,:))+diff(jfreq_ij(nodes_ix,:));
        end
        fprintf('   Choosing between time and phase gradient smoothing (elapsed time=%ds)\n',round(toc))
        std_noise=std(dph_noise,0,2);
        dph_noise2=angle(exp(j*(dph_space_uw-dph_smooth_uw2)));
        std_noise2=std(dph_noise2,0,2);
        dph_noise2(std_noise2>1.3,:)=nan;
        shaky_ix=isnan(std_noise) | std_noise>std_noise2; % spatial smoothing works better index

        fprintf('   %d arcs smoothed in time, %d in space (elapsed time=%ds)\n',ui.n_edge-sum(shaky_ix),sum(shaky_ix),round(toc))

        dph_noise(shaky_ix,:)=dph_noise2(shaky_ix,:);
        dph_space_uw(shaky_ix,:)=dph_smooth_uw2(shaky_ix,:)+dph_noise2(shaky_ix,:);
        spread(shaky_ix,:)=spread2(shaky_ix,:);
    else
        shaky_ix=[];
    end

    save('uw_space_time','dph_space_uw','dph_noise','G','spread','ifreq_ij','jfreq_ij','shaky_ix','predef_ix');
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function cinh_test ( )

%*****************************************************************************80
%
%% CINH_TEST tests R4_CINH and R8_CINH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CINH_TEST:\n' );
  fprintf ( 1, '  Test CINH_VALUES, R4_CINH, R8_CINH.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X         CINH(X)\n' );
  fprintf ( 1, '                    R4_CINH(X)        Diff\n' );
  fprintf ( 1, '                    R8_CINH(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = cinh_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_cinh ( single ( x ) );
    fx3 = r8_cinh ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% How to load the data, deconvolve it
clear all;
DataFolder = '\\feevault\data0-shared\shared\EmilyShijieShared\CaExtraction\IsolatePreTut\6991';
cnmfeFilePath = fullfile('\\feevault\data0-shared\shared\EmilyShijieShared\CaExtraction\IsolatePreTut\6991', 'CNMFE_BatchVer.mat');
load(fullfile(DataFolder, 'singinginfo.mat'))
load(cnmfeFilePath)
%% compile
L =1; % units of seconds
SONGfs = 1/diff(singinginfo(1).SpecTime(1:2));
VIDEOfs = singinginfo(1).VIDEOfs;
Lsong = L*SONGfs;
Lneural = L*VIDEOfs;

% for each file
% drop long gaps from spectrogram and neural
% concatenate song and neural, and keep track of new seg times
% do extraction

Fnums = 1:length(neuron_batch);
compSONG = [];
compNEURO = [];
FirstOfFileNeuro = [];
FirstOfFileSong = [];
for fi = 1:length(Fnums)
    filei = Fnums(fi);
    RawSong = singinginfo(filei).SongSpec;
    RawNeural =  neuron_batch(filei).DeconvSpiketrain(:,:);
    RawNeural(isnan(RawNeural)) = 0; % deconv sometimes makes nans
    [N,T] = size(RawNeural);

    if length(singinginfo(filei).segs)>0
        % take out long gaps in neural
        bouts = SegsToBouts(singinginfo(filei).segs,.05*singinginfo(filei).SOUNDfs);
        onsetTimes = floor(bouts(:,1)*singinginfo(filei).VIDEOfs/singinginfo(filei).SOUNDfs);
        offsetTimes = ceil(bouts(:,2)*singinginfo(filei).VIDEOfs/singinginfo(filei).SOUNDfs);
        % BUFFER = nan(size(NEURAL,1),Lneural);%.*mean(SPEC,2);
        NEURAL1 = [];
        for segi = 1:size(bouts,1)
            onsetTime =max(onsetTimes(segi),1);
            offsetTime =min(offsetTimes(segi),size(RawNeural,2));
            % expand spectrogram at syll onset
            NEURAL1 = [NEURAL1 RawNeural(:,onsetTime:offsetTime)];
        end
        NEURAL = NEURAL1;
%         NEURAL(NEURAL>=1) = 1; % binarize, and scale for plotting
%         NEURAL(:,end-L:end) = 0; % can't fit last L timepoints, so pad with zero
        FirstOfFileNeuro(fi) = size(compNEURO,2)+1;
        compNEURO = [compNEURO NEURAL];

        % take out long gaps in song
        onsetTimes = ceil((onsetTimes)/singinginfo(filei).VIDEOfs/diff(singinginfo(filei).SpecTime(1:2)));
        offsetTimes = floor((offsetTimes+1)/singinginfo(filei).VIDEOfs/diff(singinginfo(filei).SpecTime(1:2)));
        SONG1 = [];
        for segi = 1:size(bouts,1)
            onsetTime =max(onsetTimes(segi),1);
            offsetTime =min(offsetTimes(segi), size(RawSong,2));
            % expand spectrogram at syll onset
            SONG1 = [SONG1  RawSong(:,onsetTime:offsetTime)];
        end
        SONG = SONG1;
        FirstOfFileSong(fi) = size(compSONG,2)+1;
        compSONG = [compSONG SONG];
        figure(5); imagesc(compSONG); drawnow; shg
    end
end
NEURAL = compNEURO;
%
NEURAL = NEURAL./(prctile(NEURAL,100,2)+prctile(NEURAL(:),95));
NEURAL(max(NEURAL,[],2)<.5,:) = [];
SONG = compSONG;
SONG = SONG-prctile(SONG(:),70); SONG(SONG<0)=0;
SONG = SONG/max(SONG(:));
%% break data into training set and test set
splitN = floor(size(NEURAL,2)*.75);
splitS = floor(size(SONG,2)*.75);
trainNEURAL = NEURAL(:,1:splitN);
trainSONG = SONG(:,1:splitS);
testNEURAL = NEURAL(:,(splitN+1):end);
testSONG = SONG(:,(splitS+1):end);

%% choosing lambda
K = 10;
X = trainNEURAL;
lambdas = sort([logspace(-1,-5,100)], 'ascend');
loadings = [];
regularization = [];
cost = [];
for li = 1:length(lambdas)
    X = trainNEURAL;
    [N,T] = size(X);
    [W, H, ~,loadings(li,:),power]= seqNMF(X,'K',K,'L',Lneural,...
        'lambdaL1W', .1, 'lambda', lambdas(li), 'maxiter', 100, 'showPlot', 0);
    [cost(li),regularization(li),~] = helper.get_seqNMF_cost(X,W,H);
    li
end
%% plot costs as a function of lambda
R = (regularization-min(regularization)); R = R./mean(R);
C = (cost -min(cost)); C = C./mean(C);

clf; hold on
scatter(lambdas, R, 'b', 'markerfacecolor', 'flat');
scatter(lambdas, C, 'r', 'markerfacecolor', 'flat');
xlabel('Lambda'); ylabel('Cost (au)')
set(legend('Regularization cost', 'Reconstruction cost'), 'Box', 'on')
set(gca, 'xscale', 'log', 'ytick', [], 'color', 'none')

savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'IsolateSequences_Feb13';
thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [4 3];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_raw_costvslambda.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end
readytosave = 0
%% chose lambda=.004; run multiple times, see loadings
loadings = [];
pvals = [];
is_significant = [];
X = trainNEURAL;
for iteri = 1:100
    [W, H, ~,loadings(iteri,:),power]= seqNMF(X,'K',K,'L',Lneural,...
            'lambdaL1W', .1, 'lambda', .004, 'maxiter', 100, 'showPlot', 0);
    p = .01;
    [pvals(iteri,:),is_significant(iteri,:)] = test_significance(testNEURAL,W,p, 5000)
    W = W(:,is_significant(iteri,:)==1,:);
    H = H(is_significant(iteri,:)==1,:);
    [max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
    indSort = hybrid(:,3);
    tstart = 300;
    clf; WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
    iteri
end
%% plot significant loadings
clf; hold on
h = histogram(sum(is_significant,2), 'edgecolor', 'w', 'facecolor', .7*[1 1 1]);
h.BinCounts = h.BinCounts/sum(h.BinCounts)*100;
xlim([0 10]);
xlabel('# significant factors')
ylabel('% seqNMF runs')
% %%

% XX = repmat(1:size(loadings,1),size(loadings,1),1)+.25*rand(size(loadings,1))-.125;
% XX = XX(:);
% LL = loadings(:);
%
% scatter(XX(is_significant(:)==1),LL(is_significant(:)==1), 'k*')
% scatter(XX(is_significant(:)==0),LL(is_significant(:)==0), 'ko')
% xlabel('Sorted factor #'); ylabel('Loading')

savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'IsolateSequences_Feb13';
thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [3 2];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_NumSigFacs.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end
readytosave = 0

%% plot one example factorization
rng(122); %
X = trainNEURAL;
[W, H, ~,loadings(iteri,:),power]= seqNMF(X,'K',K,'L',Lneural,...
            'lambdaL1W', .1, 'lambda', .004, 'maxiter', 100, 'showPlot', 1);
p = .01;
[pvals,is_significant] = test_significance(testNEURAL,W,p, 5000)
W = W(:,is_significant,:);
H = H(is_significant,:);
[max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
indSort = hybrid(:,3);
tstart = 300;
clf; WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
%
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'IsolateSequences_Feb13';
thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_RawSorted.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end

clf; WHPlot(W(indSort,:,:),H(:,tstart:end), [], 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'IsolateSequences_Feb13';
thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_Reconstruction.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end

indperm = indSort(randperm(N));
clf; WHPlot(W(indperm,:,:),H(:,tstart:end), X(indperm,tstart:end), 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
saveTitle = 'IsolateSequences_Feb13';
thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [8 6];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_RawUnsorted.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end
readytosave = 0
% %% run seqNMF on neural data training set
% rng(111)
% X = trainNEURAL;
% [N,T] = size(X);
% K = 3;
% clf; shg
% [W,H,errneural] = seqNMF(X,'K',K,'L',Lneural,...
%     'lambdaL1W', .1, 'lambda', .01, 'maxiter', 150);
% [max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
% indSort = hybrid(:,3);
%
% tstart = 300;
% clf; WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
% % clf; WHPlot(W(indSort,:,:),H, X(indSort,:), 0,trainSONG)
% savedir = 'C:\Users\emackev\Dropbox (MIT)\SeqNMF\Figures';
% saveTitle = 'IsolateSequences_Feb13';
% thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
% papersize = [6 4];
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
% if readytosave
%     set(gcf, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_raw_sorted.pdf']));
%     copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
% end
% readytosave = 0

% grayscaleWHPlot(WTWO,HTWO, double(NEURAL>0), SONG, VIDEOfs/SONGfs)

% %% testing significance
% p = .05;
% [pvals,is_significant] = test_significance(testNEURAL,W,p, 5000)

%% Factor-triggered song examples and rastors
HTriggeredSpec(H,trainSONG,VIDEOfs,SONGfs,Lsong);

thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [3 2];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_HtriggeredSpec.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end

HTriggeredRaster(H,trainNEURAL(indSort,:),Lneural)
thisFile = 'C:\Users\emackev\Documents\MATLAB\seqnmf\TestbedFigs_elm.m';
papersize = [3 3];
set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
if readytosave
    set(gcf, 'color', 'none')
    saveas(gcf, fullfile(savedir, [saveTitle '_HtriggeredRaster.pdf']));
    copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
end

readytosave = 0

%
%
%
%
%
% HTriggeredRaster(HTWO,NEURAL(indSort,:),Lneural);
% papersize = [3 7]; %[185.5 116.17*3/4]/30;
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
% if readytosave
%     set(gcf, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_rasters.pdf']));
%     copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
% end
% clf
% HTriggeredSpec(HTWO,SONG,VIDEOfs,SONGfs,Lsong)
% papersize = [3 4]; %[185.5 116.17*3/4]/30;
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
% if readytosave
%     set(gcf, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_specs.pdf']));
%     copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
%     set(gcf, 'color', 'w')
% end
% readytosave = 0;
% %%
% helper.AutoSelectLambda(NEURAL, 6, Lneural)
% papersize = [4 3];
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
% if readytosave
%     set(gcf, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_choosinglambda.pdf']));
%     copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
%     set(gcf, 'color', 'w')
% end
% readytosave = 0;
% %% compute loadings for multiple runs
% rng(3333)
%
% [N,T] = size(NEURAL);
% K = 15;
% config = struct();
% config.ploton = 0;
% config.maxiter = 100;
% config.lambdaW = .002;
% config.competeW = 0;
% config.lambdaL1W = .25;
% config.lambdaL1H = 0;
% config.shift = 1;
% clf; shg
% for iteri = 1:10
%     config.lambdaW = .002;
%     config.lambdaL1W = .25;
% %     [Wneural,Hneural,errneural] = algorithms.seq_NMF(NEURAL,K,Lneural,config);
% %     LOADING(iteri,:) = helper.computeLoadingPercentPower(NEURAL,Wneural,Hneural);
%     config.lambdaW = 0;
%     config.lambdaL1W = 0;
%     [Wneural,Hneural,errneural] = algorithms.seq_NMF(NEURAL,K,Lneural,config);
%     LOADING_cnmf(iteri,:) = helper.computeLoadingPercentPower(NEURAL,Wneural,Hneural);
%     iteri
% end
% %%
% clf
% X = repmat(1:K,iteri,1) + .5*rand(iteri,K)-.25;
% scatter(X(:), LOADING_cnmf(:), 'k', 'markerfacecolor', 'flat')
% hold on
% scatter(X(:), LOADING(:), 'r', 'markerfacecolor', 'flat')
% legend('CNMF', 'seqNMF')
% xlabel('Factor #')
% ylabel('Loading (% power explained)');
% %
% papersize = [3 2.5];
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
% if readytosave
%     set(gcf, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_Loadings.pdf']));
%     copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
%     set(gcf, 'color', 'w')
% end
% readytosave = 0;
% %% compute rmse & pow explained for multiple K's
% rng(3333)
%
% [N,T] = size(NEURAL);
% K = 10;
% config = struct();
% config.ploton = 0;
% config.maxiter = 100;
% config.lambdaW = .002;
% config.competeW = 0;
% config.lambdaL1W = .25;
% config.lambdaL1H = 0;
% config.shift = 1;
% clf; shg
% RMSE = [];
% RMSE_cnmf = [];
% POWEXP = [];
% POWEXP_cnmf = [];
% W = {};
% H = {};
% W_cnmf = {};
% H_cnmf = {};
% Score = [];
% Score_cnmf = [];
%
%
% for Ki = 1:K
%     for iteri = 1:10
%         config.lambdaW = .002;
%         config.lambdaL1W = .25;
%         [W{iteri,Ki},H{iteri,Ki},RMSE(iteri,Ki),~,POWEXP(iteri,Ki)] = algorithms.seq_NMF(NEURAL,Ki,Lneural,config);
%         config.lambdaW = 0;
%         config.lambdaL1W = 0;
%         [W_cnmf{iteri,Ki},H_cnmf{iteri,Ki},RMSE_cnmf(iteri,Ki),~,POWEXP_cnmf(iteri,Ki)] = algorithms.seq_NMF(NEURAL,Ki,Lneural,config);
%         iteri
%         RMSE
%         RMSE_cnmf
%     end
%     Score(:,Ki) = helper.consistency(W(:,Ki), H(:,Ki));
%     Score_cnmf(:,Ki) = helper.consistency(W_cnmf(:,Ki), H_cnmf(:,Ki));
% end
% %% plot it
% clf
% subplot(1,3,1); cla
% X = repmat(1:K,iteri,1) + .5*rand(iteri,K)-.25;
% scatter(X(:), RMSE_cnmf(:), 5,'k', 'markerfacecolor', 'flat')
% hold on;
% scatter(X(:), RMSE(:), 5,'r', 'markerfacecolor', 'flat')
% % legend('CNMF', 'seqNMF')
% axis tight;
% ylims = ylim;
% ylim([0 ylims(2)])
% xlabel('K')
% ylabel('RMSE');
%
% subplot(1,3,2); cla
% X = repmat(1:K,iteri,1) + .5*rand(iteri,K)-.25;
% scatter(X(:), POWEXP_cnmf(:), 5,'k', 'markerfacecolor', 'flat')
% hold on
% scatter(X(:), POWEXP(:), 5,'r', 'markerfacecolor', 'flat')
% % legend('CNMF', 'seqNMF')
% axis tight;
% ylims = ylim;
% ylim([0 ylims(2)])
% xlabel('K')
% ylabel('Power explained (%)');
% %
% subplot(1,3,3); cla
% patches = plot.violin(Score_cnmf+rand(size(Score))*.00001, 'facecolor','k');
% for ki = 1:K
%     patches(ki).Vertices(:,1) = patches(ki).Vertices(:,1) -.25;
% end
% hold on
% plot.violin(Score+rand(size(Score))*.00001, 'facecolor','r');
% ylims = ylim;
% ylim([0 ylims(2)])
% % X = repmat(1:K,size(Score,1),1) + .5*rand(size(Score,1),K)-.25;
% % scatter(X(:), Score_cnmf(:), 2,'k', 'markerfacecolor', 'flat')
% % hold on
% % scatter(X(:), Score(:), 2,'r', 'markerfacecolor', 'flat')
% % legend('CNMF', 'seqNMF')
% legend off; box off
% xlabel('K')
% ylabel('Consistency');
%
% papersize = [8 2];
% set(gcf, 'papersize', papersize, 'paperposition', [0 0 papersize]);%, 'color', 'none')
% if readytosave
%     set(gcf, 'color', 'none')
%     saveas(gcf, fullfile(savedir, [saveTitle '_chooseK.pdf']));
%     copyfile(thisFile, fullfile(savedir, [saveTitle '.m']));
%     set(gcf, 'color', 'w')
% end
% readytosave = 0;

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function mono_between_next_grlex_test ( )

%*****************************************************************************80
%
%% MONO_BETWEEN_NEXT_GRLEX_TEST tests MONO_BETWEEN_NEXT_GRLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_BETWEEN_NEXT_GRLEX_TEST\n' );
  fprintf ( 1, '  MONO_BETWEEN_NEXT_GRLEX can list the monomials\n' );
  fprintf ( 1, '  in M variables, of total degree N between N1 and N2,\n' );
  fprintf ( 1, '  one at a time, in graded lexicographic order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We start the process with (0,0,...,0,N1).\n' );
  fprintf ( 1, '  The process ends with (N2,0,...,0,0)\n' );

  n1 = 2;
  n2 = 3;
  m = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Let M =  %d\n', m );
  fprintf ( 1, '      N1 = %d\n', n1 );
  fprintf ( 1, '      N2 = %d\n', n2 );
  fprintf ( 1, '\n' );

  x = [ 0, 0, n1 ];
  i = 1;

  while ( 1 )

    fprintf ( 1, '  %2d:', i );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );

    if ( x(1) == n2 )
      break
    end

    x = mono_between_next_grlex ( m, n1, n2, x );
    i = i + 1;

  end

  return
end


text:
function vid_bic = interp_bicubic(vid_l, upscale)

N = length(vid_l);

for i = 1: N
    vid_bic{i} = imresize(vid_l{i}, upscale, 'bicubic');
    %vid_bic{i} = imresize(vid_l{i}, upscale, 'nearest');
end

author: thomas-koehler

repo: SupER


text:
function I = haze_linear(R, t, L)
%HAZE_LINEAR  Generate hazy image from clean image using the linear haze model
%corresponding to Lambert-Beer law
%   Inputs:
%       -|R|: H-by-W-by-|image_channels| clean image representing true radiance
%       of scene.
%       -|t|: H-by-W transmission map.
%       -|L|: 1-by-1-by-|image_channels| homogeneous atmospheric light.
%
%   Outputs:
%       -|I|: synthetic hazy image, with same size as input clean image |R|.

image_channels = size(L, 3);

% Auxiliary matrix with replicates of transmission map for all channels,
% allowing to express hazy image conveniently.
t_replicated = repmat(t, 1, 1, image_channels);

% Apply linear haze model.
I = t_replicated .* R + (1 - t_replicated) .* repmat(L, size(t));

end

author: sakaridis


text:
function power_method_test ( )

%*****************************************************************************80
%
%% POWER_METHOD_TEST tests the POWER_METHOD library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_METHOD_TEST\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Test the POWER_METHOD library.\n' );

  power_method_test01 ( );
  power_method_test02 ( );
  power_method_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_METHOD_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function power_method_test01 ( )

%*****************************************************************************80
%
%% POWER_METHOD_TEST01 uses POWER_METHOD on the Fibonacci2 matrix.
%
%  Discussion:
%
%    This matrix, despite having a single dominant eigenvalue, will generally
%    converge only very slowly under the power method.  This has to do with
%    the fact that the matrix has only 3 eigenvectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2008
%
%  Author:
%
%    John Burkardt
%
  n = 50;

  a = fibonacci2 ( n );

  seed = 123456789;

  [ x, seed ] = r8vec_uniform_01 ( n, seed );

  it_max = 300;
  tol = 0.000001;

  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_METHOD_TEST01\n' );
  fprintf ( 1, '  Use the power method on the Fibonacci2 matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N       = %d\n', n );
  fprintf ( 1, '  Maximum iterations   = %d\n', it_max );
  fprintf ( 1, '  Error tolerance      = %g\n', tol );

  ctime1 = cputime;

  [ x, lambda, it_num ] = power_method ( n, a, x, it_max, tol );

  ctime2 = cputime;
  ctime = ctime2 - ctime1;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', it_num );
  fprintf ( 1, '  CPU time             = %f\n', ctime );
  fprintf ( 1, '  Estimated eigenvalue = %f\n', lambda );
  fprintf ( 1, '  Correct value        = %f\n', phi );
  fprintf ( 1, '  Error                = %e\n', abs ( lambda - phi ) );
%
%  X2 is the exact eigenvector.
%
  x2(1:n,1) = phi.^(0:n-1);
  x2 = x2 / norm ( x2 );
%
%  The sine of the angle between X and X2 is a measure of error.
%
  cos_x1x2 = x' * x2;
  sin_x1x2 = sqrt ( ( 1.0 - cos_x1x2 ) * ( 1.0 + cos_x1x2 ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sine of angle between true and estimated vectors = %e\n', ...
    sin_x1x2 );

  return
end
function power_method_test02 ( )

%*****************************************************************************80
%
%% POWER_METHOD_TEST02 uses POWER_METHOD2 on the Fibonacci2 matrix.
%
%  Discussion:
%
%    This matrix, despite having a single dominant eigenvalue, will generally
%    converge only very slowly under the power method.  This has to do with
%    the fact that the matrix has only 3 eigenvectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2008
%
%  Author:
%
%    John Burkardt
%
  n = 50;

  a = fibonacci2 ( n );

  seed = 123456789;

  [ x, seed ] = r8vec_uniform_01 ( n, seed );

  it_max = 300;
  tol = 0.000001;

  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_METHOD_TEST02\n' );
  fprintf ( 1, '  Use the power method2 on the Fibonacci2 matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N       = %d\n', n );
  fprintf ( 1, '  Maximum iterations   = %d\n', it_max );
  fprintf ( 1, '  Error tolerance      = %g\n', tol );

  ctime1 = cputime;

  [ lambda, v, it_num ] = power_method2 ( n, a, x, it_max, tol );

  ctime2 = cputime;
  ctime = ctime2 - ctime1;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', it_num );
  fprintf ( 1, '  CPU time             = %f\n', ctime );
  fprintf ( 1, '  Estimated eigenvalue = %f  %f\n', real ( lambda ), imag ( lambda ) );
  fprintf ( 1, '  Correct value        = %f\n', phi );
  fprintf ( 1, '  Error                = %e\n', abs ( lambda - phi ) );

  return
end
function power_method_test03 ( )

%*****************************************************************************80
%
%% POWER_METHOD_TEST03 uses POWER_METHOD2 on the TRIS matrix.
%
%  Discussion:
%
%    This matrix, despite having a single dominant eigenvalue, will generally
%    converge only very slowly under the power method.  This has to do with
%    the fact that the matrix has only 3 eigenvectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2008
%
%  Author:
%
%    John Burkardt
%
  n = 50;

  alpha = -1.0;
  beta = 10.0;
  gamma = 8.0;

  a = tris ( n, n, alpha, beta, gamma );

  seed = 123456789;

  [ x, seed ] = r8vec_uniform_01 ( n, seed );

  it_max = 4000;
  tol = 0.000001;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_METHOD_TEST03\n' );
  fprintf ( 1, '  Use the power method2 on the TRIS matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N       = %d\n', n );
  fprintf ( 1, '  Maximum iterations   = %d\n', it_max );
  fprintf ( 1, '  Error tolerance      = %g\n', tol );

  ctime1 = cputime;

  [ lambda, v, it_num ] = power_method2 ( n, a, x, it_max, tol );

  ctime2 = cputime;
  ctime = ctime2 - ctime1;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', it_num );
  fprintf ( 1, '  CPU time             = %f\n', ctime );
  fprintf ( 1, '  Estimated eigenvalue = %f  %f\n', ...
    real ( lambda ), imag ( lambda ) );

  lambda_vec = tris_eigenvalues ( n, alpha, beta, gamma );

  lambda_max = lambda_vec(1);

  for i = 2 : n
    if ( abs ( lambda_max ) < abs ( lambda_vec(i) ) )
      lambda_max = lambda_vec(i);
    end
  end

  fprintf ( 1, '  Correct value        = %f  %f\n', ...
    real ( lambda_max ), imag ( lambda_max ) );
  fprintf ( 1, '  Error                = %e\n', abs ( lambda - lambda_max ) );

  return
end
function a = fibonacci2 ( n )

%*****************************************************************************80
%
%% FIBONACCI2 returns the Fibonacci2 matrix.
%
%  Example:
%
%    N = 5
%
%    0 1 0 0 0
%    1 1 0 0 0
%    0 1 1 0 0
%    0 0 1 1 0
%    0 0 0 1 1
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is banded, with bandwidth 3.
%
%    A is integral: int ( A ) = A.
%
%    A is a zero/one matrix.
%
%    If N = 1 then
%      det ( A ) = 0
%    else
%      det ( A ) = (-1)**(N-1)
%
%    If 1 < N, then A is unimodular.
%
%    For 2 <= N, A has the eigenvalues:
%
%      PHI   (once),
%      1     (N-2) times,
%      1-PHI (once).
%
%    When applied to a Fibonacci1 matrix B, the Fibonacci2 matrix
%    A produces the "next" Fibonacci1 matrix C = A*B.
%
%    Let PHI be the golden ratio (1+sqrt(5))/2.
%
%    For 2 <= N, the eigenvalues and eigenvectors are:
%
%    LAMBDA(1)     = PHI,     vector = (1,PHI,PHI^2,...PHI^(N-1));
%    LAMBDA(2:N-1) = 1        vector = (0,0,0,...,0,1);
%    LAMBDA(N)     = 1 - PHI. vector = ((-PHI)^(N-1),(-PHI)^(N-2),...,1)
%
%    Note that there is only one eigenvector corresponding to 1.
%    Hence, for 3 < N, the matrix is defective.  This fact means,
%    for instance, that the convergence of the eigenvector in the power
%    method will be very slow.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  for i = 1 : n
    for j = 1 : n

      if ( i == 1 )

        if ( j == 2 )
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end

      else

        if ( j == i-1 | j == i )
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end

      end

    end
  end

  return
end
function [ r, seed ] = r8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real R(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM_01 - Fatal error!' );
  end

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + i4_huge;
    end

    r(i) = seed * 4.656612875E-10;

  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
function a = tris ( m, n, x, y, z )

%*****************************************************************************80
%
%% TRIS returns the tridiagonal scalar matrix.
%
%  Formula:
%
%    if ( J = I-1 )
%      A(I,J) = X
%    elseif ( J = I )
%      A(I,J) = Y
%    elseif ( J = I + 1 )
%      A(I,J) = Z
%    else
%      A(I,J) = 0
%
%  Example:
%
%    M = 5, N = 5, X = 1, Y = 2, Z = 3
%
%    2 3 0 0 0
%    1 2 3 0 0
%    0 1 2 3 0
%    0 0 1 2 3
%    0 0 0 1 2
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is banded, with bandwidth 3.
%
%    A is Toeplitz: constant along diagonals.
%
%    If Y is not zero, then for A to be singular, it must be the case that
%
%      0.5 * Y / sqrt ( X * Z ) < 1
%
%    and
%
%      cos (K*PI/(N+1)) = - 0.5 * Y / sqrt ( X * Z ) for some 1 <= K <= N.
%
%    If Y is zero, then A is singular when N is odd, or if X or Z is zero.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A has eigenvalues
%
%      LAMBDA(I) = Y + 2 * sqrt(X*Z) * COS(I*PI/(N+1))
%
%    The eigenvalues will be complex if X * Z < 0.
%
%    If X = Z, the matrix is symmetric.
%
%    As long as X and Z are nonzero, the matrix is irreducible.
%
%    If X = Z = -1, and Y = 2, the matrix is a symmetric, positive
%    definite M matrix, the negative of the second difference matrix.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Todd,
%    Basic Numerical Mathematics,
%    Volume 2: Numerical Algebra,
%    Academic Press, 1978, page 155.
%
%  Parameters:
%
%    Input, integer M, N, the order of A.
%
%    Input, real X, Y, Z, the scalars that define A.
%
%    Output, real A(M,N), the matrix.
%
  for i = 1 : m
    for j = 1 : n

      if ( j == i - 1 )
        a(i,j) = x;
      elseif ( j == i )
        a(i,j) = y;
      elseif ( j == i + 1 )
        a(i,j) = z;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end
function lambda = tris_eigenvalues ( n, x, y, z )

%*****************************************************************************80
%
%% TRIS_EIGENVALUES returns the eigenvalues of the tridiagonal scalar matrix.
%
%  Discussion:
%
%    The eigenvalues will be complex if X * Z < 0.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X, Y, Z, the scalars that define A.
%
%    Output, complex LAMBDA(N), the eigenvalues.
%
  for i = 1 : n
    angle = i * pi / ( n + 1 );
    lambda(i) = y + 2.0 * sqrt ( x * z ) * cos ( angle );
  end

  return
end


text:
function cadj=ref_spreadadj_3(coef)
%REF_SPREADADJ_3  Symbol of adjoint spreading function.
%   Usage: cadj=ref_spreadadj_3(c,number);
%
%   Development version by FJ for comparison of different implementations
%   cadj=SPREADADJ(c) will compute the symbol cadj of the spreading
%   operator that is the adjoint of the spreading operator with symbol c.
%
%   This is an improved implementation of the direct formula with low memory
%   needs
%
%   This implementation uses the direct formula given in case 2 with
%   the following Optimizations :
%
%-     Avoiding mod : In the loop of case 2, we see that
%      mod(L-ii,L)~=L-ii only for ii==0 (idem for jj), so we can
%      remove the mod by processing separetly the cases ii==0 or
%      jj==0.
%
%-     Precomputation of exp : In the loop of case 2, we see that we
%      compute many time complex exponential terms with the same
%      values. Using precomputation and modulo, we can reduce the
%      computation time
%-
%

L=size(coef,1);

cadj=zeros(L);

% Proceesing for ii==0 or jj==0
cadj(1,1)=conj(coef(1,1));
cadj(2:end,1)=conj(coef(end:-1:2,1));
cadj(1,2:end,1)=conj(coef(1,end:-1:2));

% Proceesing for ii~=0 and jj~=0

% Precomputation for exponential term
% Optimization note : here we are computing the Lth root of unity
% which have many known special properties and symetries that
% could be exploited to highly reduce the computation for the
% following line (this is not a critical part but some
% improvements Optimizations are easy to identify)
temp=exp((-i*2*pi/L)*(0:L-1));

for ii=1:L-1
  for jj=1:L-1
    cadj(ii+1,jj+1)=conj(coef(L-ii+1,L-jj+1))...
        *temp(mod(ii*jj,L)+1);
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%DOUBLEM Datafile mapping for conversion to double
%
%	B = DOUBLEM(A)
%	B = A*DOUBLEM
%
% For datasets B = A, but the data inside A is converted to double.
% For datafiles B = A, but the command makes sure that the data in B is
% first converted to double when B is converted to a dataset.
% In all other cases A itself is converted to double.
% This is useful to convert INT8 and other formats inside a dataset.

function a = doublem(a)

	if nargin < 1 || isempty(a)
		a = prmapping(mfilename,'fixed');
		a = setname(a,'double');
	elseif isdataset(a)
		a = setdat(a,double(+a));
	elseif isdatafile(a)
		a = a*filtm([],'double');
	else
		a = double(a);
	end

return

author: marianux

repo: ecg-kit

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


text:
function setup_cnnimageretrieval()
% SETUP_CNNIMAGERETRIEVAL Setup the toolbox.
%
%   setup_cnnimageretrieval()  Adds the toolbox to MATLAB path.
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

[root] = fileparts(mfilename('fullpath')) ;

% add paths from this package
addpath(fullfile(root, 'cnnblocks'));
addpath(fullfile(root, 'cnninit'));
addpath(fullfile(root, 'cnntrain'));
addpath(fullfile(root, 'cnnvecs'));
addpath(fullfile(root, 'examples'));
addpath(fullfile(root, 'whiten'));
addpath(fullfile(root, 'utils'));

author: filipradenovic

repo: cnnimageretrieval


text:
function value = rhs_test03 ( x_num, x, t )

%*****************************************************************************80
%
%% RHS_TEST03 evaluates the right hand side for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM,1), the node coordinates.
%
%    Input, real T, the current time.
%
%    Output, real VALUE(X_NUM,1), the source term.
%
  value(1:x_num,1) = 0.0;

  return
end


text:
function mm_header_print ( input_filename, id, type, rep, field, symm )

%*****************************************************************************80
%
%% MM_HEADER_PRINT prints header information from a Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ( len = * ) INPUT_FILENAME, the input file name.
%
%    Input, character ( len = 14 ) ID, the Matrix Market identifier.
%    This value must be '%%MatrixMarket'.
%
%    Input, character ( len = 6 ) TYPE, the Matrix Market type.
%    This value must be 'matrix'.
%
%    Input, character ( len = 10 ) REP, the Matrix Market 'representation'
%    indicator.  Possible values include:
%    'coordinate'   (for sparse data)
%    'array'        (for dense data)
%    'elemental'    (to be added)
%
%    Input, character ( len = 7 ) FIELD, the Matrix Market 'field'.
%    Possible values include:
%    'real'
%    'double'
%    'complex'
%    'integer'
%    'pattern'
%
%    Input, character ( len = 19 ) SYMM, the Matrix Market symmetry.
%    Possible values include:
%    'symmetric'
%    'hermitian'
%    'skew-symmetric'
%    'general'
%
  fprintf ( 1, '\n');
  fprintf ( 1, 'MM_HEADER_PRINT:\n' );
  fprintf ( 1, '  Header information from Matrix Market file "%s".\n', ...
    input_filename );
  fprintf ( 1, '\n');
  fprintf ( 1, '  Matrix Market ID = "%s".\n', id );
  fprintf ( 1, '    "%%%%MatrixMarket" is only allowed value.\n' );
  fprintf ( 1, '\n');
  fprintf ( 1, '  Matrix Market TYPE = "%s".\n', type );
  fprintf ( 1, '    "matrix" is only allowed value.\n' );
  fprintf ( 1, '\n');
  fprintf ( 1, '  Representation type REP = "%s".\n', rep );
  fprintf ( 1, '    "coordinate" for sparse matrices,\n' );
  fprintf ( 1, '    "array"      for dense matrices,\n' );
  fprintf ( 1, '    "elemental"  for unassembled finite element matrices.\n' );
  fprintf ( 1, '\n');
  fprintf ( 1, '  Numeric FIELD = "%s".\n', field );
  fprintf ( 1, '    "integer" for integer values,\n' );
  fprintf ( 1, '    "real"    for real values,\n' );
  fprintf ( 1, '    "double"  for double precision real values,\n' );
  fprintf ( 1, '    "complex" for complex values,\n' );
  fprintf ( 1, '    "pattern" for nonzero pattern only.\n' );
  fprintf ( 1, '\n');
  fprintf ( 1, '  Symmetry type SYMM = "%s".\n', symm );
  fprintf ( 1, '    "general"         no symmetry,\n' );
  fprintf ( 1, '    "symmetric"       A(I,J) = A(J,I),\n' );
  fprintf ( 1, '                      input only lower triangle.\n' );
  fprintf ( 1, '    "skew-symmetric"  A(I,J) = - A(J,I),\n' );
  fprintf ( 1, '                      input only strict lower triangle.\n' );
  fprintf ( 1, '    "Hermitian"       A(I,J) = A*(J,I),\n' );
  fprintf ( 1, '                      input only lower triangle.\n' );

  return
end


text:
function [nodes2, edges2] = grMergeNodesMedian(nodes, edges, mnodes)
%GRMERGENODESMEDIAN Replace several nodes by their median coordinate
%
%   [NODES2, EDGES2] = grMergeNodesMedian(NODES, EDGES, NODES2MERGE)
%   NODES ans EDGES are the graph structure, and NODES2MERGE is the list of
%   indices of nodes to be merged.
%   The median coordinate of merged nodes is computed, and all nodes are
%   merged to this new node.
%
%

%   -----
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 13/08/2003.
%

%   HISTORY
%   10/02/2004 : documentation

% coordinates of reference node
x = median(nodes(mnodes, 1));
y = median(nodes(mnodes, 2));

% index of reference node
refNode = findPoint([x y], nodes);
mnodes = sort(mnodes(mnodes ~= refNode));

for n = 1:length(mnodes)
    node = mnodes(n);

    % process each neighbor of the current node
    neighbors = grAdjacentNodes(edges, node);
    for e = 1:length(neighbors)
        edge = neighbors(e);

        if edges(edge, 1) == refNode || edges(edge, 2) == refNode
            continue;
        end

        % find if the node is referenced as 1 or 2 in the edge,
        % and replace it with the reference node.
        if edges(edge, 1) == node
            edges(edge, 1) = refNode;
        else
            edges(edge, 2) = refNode;
        end

    end
end

% remove nodes from the list, except the reference node.
for n = 1:length(mnodes)
    [nodes, edges] = grRemoveNode(nodes, edges, mnodes(n)-n+1);
end

nodes2 = nodes;
edges2 = edges;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% 15.7.02
%
% Unscented Kalman Filter (UKF) applied to FitzHugh-Nagumo neuron dynamics.
% Voltage observed, currents and inputs estimated.
%
% FitzHughNagumo() is the main program and calls the other programs.
%
% A detailed description is provided in
% H.U. Voss, J. Timmer & J. Kurths, Nonlinear dynamical system identification from uncertain and indirect measurements, Int. J. Bifurcation and Chaos 14, 1905-1933 (2004).
% I will be happy to email this paper on request. It contains a tutorial about the estimation of hidden states and unscented Kalman filtering.
%
% For commercial use and questions, please contact me.
%
% ++++++++++++++++++++++++++++++++++++++++++++++
% Henning U. Voss, Ph.D.
% Associate Professor of Physics in Radiology
% Citigroup Biomedical Imaging Center
% Weill Medical College of Cornell University
% 516 E 72nd Street
% New York, NY 10021
% Tel. 001-212 746-5216, Fax. 001-212 746-6681
% Email: hev2006@med.cornell.edu
% ++++++++++++++++++++++++++++++++++++++++++++++

function FitzHughNagumo()

global dT

orient tall

% Dimensions: dq for param. vector, dx augmented state, dy observation
dq=1; dx=dq+2; dy=1;

fct='FitzHughNagumo_fct'; % model function F(x)
obsfct='FitzHughNagumo_obsfct';   % observation function G(x)

ll=800; % number of data samples
dT=0.2; % sampling time step (global variable)

% Simulating data: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x0=zeros(2,ll); x0(:,1)=[0; 0]; % true trajectory
dt=.1*dT; nn=fix(dT/dt); % the integration time step is smaller than dT

% External input, estimated as parameter p later on:
z=[1:ll]/250*2*pi; z=-.4-1.01*abs(sin(z/2));

% 4th order Runge-Kutta integrator:

for n=1:ll-1;
  xx=x0(:,n);
  for i=1:nn
    k1=dt*FitzHughNagumo_int(xx,z(n));
    k2=dt*FitzHughNagumo_int(xx+k1/2,z(n));
    k3=dt*FitzHughNagumo_int(xx+k2/2,z(n));
    k4=dt*FitzHughNagumo_int(xx+k3,z(n));
    xx=xx+k1/6+k2/3+k3/3+k4/6;
  end;
  x0(:,n+1)=xx;
end;

x=[z; x0]; % augmented state vector (notation a bit different to paper)

R=.2^2*var(FitzHughNagumo_obsfct(x))*eye(dy,dy); % observation noise covariance matrix
randn('state',0); y=feval(obsfct,x)+sqrtm(R)*randn(dy,ll); % noisy data

% Initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xhat=zeros(dx,ll);
xhat(:,2)=x(:,2); % first guess of x_1 set to observation

Q=.015; % process noise covariance matrix

Pxx=zeros(dx,dx,ll);
Pxx(:,:,1)=blkdiag(Q,R,R);

errors=zeros(dx,ll); % not so important
Ks=zeros(dx,dy,ll);  % Kalman gains

% Main loop for recursive estimation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=2:ll

[xhat(:,k),Pxx(:,:,k),Ks(:,:,k)]=ut(xhat(:,k-1),Pxx(:,:,k-1),y(:,k),fct,obsfct,dq,dx,dy,R);

Pxx(1,1,k)=Q;

errors(:,k)=sqrt(diag(Pxx(:,:,k)));

end; % k

% Results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

chisq=mean((x(1,:)-xhat(1,:)).^2+(x(2,:)-xhat(2,:)).^2+(x(3,:)-xhat(3,:)).^2)
est=xhat(1:dq,ll)'
error=errors(1:dq,ll)'

figure(1)

subplot(2,1,1)
plot(y,'bd','MarkerEdgeColor','blue', 'MarkerFaceColor','blue','MarkerSize',3);
hold on;
plot(x(dq+1,:),'black','LineWidth',2);
%plot(xhat(dq+1,:),'r','LineWidth',2);
xlabel(texlabel('t'));
ylabel(texlabel('x_1, y'));
hold off;
axis tight
title('(a)')
drawnow

subplot(2,1,2)
plot(x(dq+2,:),'black','LineWidth',2);
hold on
plot(xhat(dq+2,:),'r','LineWidth',2);
plot(x(1,:),'black','LineWidth',2);
for i=1:dq; plot(xhat(i,:),'m','LineWidth',2); end;
for i=1:dq; plot(xhat(i,:)+errors(i,:),'m'); end;
for i=1:dq; plot(xhat(i,:)-errors(i,:),'m'); end;
xlabel(texlabel('t'));
ylabel(texlabel('z, estimated z, x_2, estimated x_2'));
hold off
axis tight
title('(b)')

function r=FitzHughNagumo_int(x,z)
% Function for modeling data, not for UKF execution
a=.7; b=.8; c=3.;
r=[c*(x(2)+x(1)-x(1)^3/3+z); -(x(1)-a+b*x(2))/c];

function r=FitzHughNagumo_obsfct(x)
% Observation function
r=x(2,:);

function r=FitzHughNagumo_fct(dq,x)
% Runge-Kutta integrator for FitzHugh-Nagumo system with parameters

global dT
dt=.02; % local integration step
nn=fix(dT/dt);

p=x(1:dq,:);
xnl=x(dq+1:size(x(:,1)),:);
for n=1:nn
k1=dt*fc(xnl,p);
k2=dt*fc(xnl+k1/2,p);
k3=dt*fc(xnl+k2/2,p);
k4=dt*fc(xnl+k3,p);
xnl=xnl+k1/6+k2/3+k3/3+k4/6;
end
r=[x(1:dq,:); xnl];

function r=fc(x,p);
a=.7; b=.8; c=3.;
r=[c*(x(2,:)+x(1,:)-x(1,:).^3/3+p); -(x(1,:)-a+b*x(2,:))/c];

function [xhat,Pxx,K]=ut(xhat,Pxx,y,fct,obsfct,dq,dx,dy,R);
% Unscented transformation. Not specific to FitzHugh-Nagumo model

  N=2*dx;

  xsigma=chol( dx*Pxx )'; % Pxx=root*root', but Pxx=chol'*chol
  Xa=xhat*ones(1,N)+[xsigma, -xsigma];
  X=feval(fct,dq,Xa);

  xtilde=sum(X')'/N;

  Pxx=zeros(dx,dx);
  for i=1:N;
    Pxx=Pxx+(X(:,i)-xtilde)*(X(:,i)-xtilde)'/N;
  end;

  Y=feval(obsfct,X);

  ytilde=sum(Y')'/N;
  Pyy=R;
  for i=1:N;
    Pyy=Pyy+(Y(:,i)-ytilde)*(Y(:,i)-ytilde)'/N;
  end;
  Pxy=zeros(dx,dy);
  for i=1:N;
    Pxy=Pxy+(X(:,i)-xtilde)*(Y(:,i)-ytilde)'/N;
  end;

  K=Pxy*inv(Pyy);
  xhat=xtilde+K*(y-ytilde);
  Pxx=Pxx-K*Pxy';

author: Sable

repo: mcbench-benchmarks


text:
function [R,eff] = randmio_und_connected(R, ITER)
%RANDMIO_UND_CONNECTED     Random graph with preserved degree distribution
%
%   R = randmio_und_connected(W,ITER);
%   [R eff] = randmio_und_connected(W, ITER);
%
%   This function randomizes an undirected network, while preserving the
%   degree distribution. The function does not preserve the strength
%   distribution in weighted networks. The function also ensures that the
%   randomized network maintains connectedness, the ability for every node
%   to reach every other node in the network. The input network for this
%   function must be connected.
%
%   Input:      W,      undirected (binary/weighted) connection matrix
%               ITER,   rewiring parameter
%                       (each edge is rewired approximately ITER times)
%
%   Output:     R,      randomized network
%               eff,    number of actual rewirings carried out
%
%   References: Maslov and Sneppen (2002) Science 296:910
%
%
%   2007-2012
%   Mika Rubinov, UNSW
%   Jonathan Power, WUSTL
%   Olaf Sporns, IU

%   Modification History:
%   Jun 2007: Original (Mika Rubinov)
%   Apr 2008: Edge c-d is flipped with 50% probability, allowing to explore
%             all potential rewirings (Jonathan Power)
%   Mar 2012: Limit number of rewiring attempts, count number of successful
%             rewirings (Olaf Sporns)

n=size(R,1);
[i,j]=find(tril(R));
K=length(i);
ITER=K*ITER;

% maximal number of rewiring attempts per 'iter'
maxAttempts= round(n*K/(n*(n-1)));
% actual number of successful rewirings
eff = 0;

for iter=1:ITER
    att=0;
    while (att<=maxAttempts)                                     %while not rewired
        rewire=1;
        while 1
            e1=ceil(K*rand);
            e2=ceil(K*rand);
            while (e2==e1),
                e2=ceil(K*rand);
            end
            a=i(e1); b=j(e1);
            c=i(e2); d=j(e2);

            if all(a~=[c d]) && all(b~=[c d]);
                break           %all four vertices must be different
            end
        end

        if rand>0.5
            i(e2)=d; j(e2)=c; 	%flip edge c-d with 50% probability
            c=i(e2); d=j(e2); 	%to explore all potential rewirings
        end

        %rewiring condition
        if ~(R(a,d) || R(c,b))
            %connectedness condition
            if ~(R(a,c) || R(b,d))
                P=R([a d],:);
                P(1,b)=0; P(2,c)=0;
                PN=P;
                PN(:,d)=1; PN(:,a)=1;

                while 1
                    P(1,:)=any(R(P(1,:)~=0,:),1);
                    P(2,:)=any(R(P(2,:)~=0,:),1);
                    P=P.*(~PN);
                    if ~all(any(P,2))
                        rewire=0;
                        break
                    elseif any(any(P(:,[b c])))
                        break
                    end
                    PN=PN+P;
                end
            end %connectedness testing

            if rewire               %reassign edges
                R(a,d)=R(a,b); R(a,b)=0;
                R(d,a)=R(b,a); R(b,a)=0;
                R(c,b)=R(c,d); R(c,d)=0;
                R(b,c)=R(d,c); R(d,c)=0;

                j(e1) = d;          %reassign edge indices
                j(e2) = b;
                eff = eff+1;
                break;
            end %edge reassignment
        end %rewiring condition
        att=att+1;
    end %while not rewired
end %iterations

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Plane Wave Absorption Example
%
% This example illustrates the characteristics of the Kelvin-Voigt
% absorption model used in the k-Wave simulation functions pstdElastic2D,
% pstdElastic3D. It builds on the Explosive Source In A Layered Medium
% Example.
%
% author: Bradley Treeby
% date: 17th January 2014
% last update: 14th February 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SET GRID PARAMETERS
% =========================================================================

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 32;            % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed_compression = 1800;  % [m/s]
medium.sound_speed_shear       = 1200;  % [m/s]
medium.density                 = 1000;  % [kg/m^3]

% set the absorption properties
medium.alpha_coeff_compression = 1;     % [dB/(MHz^2 cm)]
medium.alpha_coeff_shear       = 1;     % [dB/(MHz^2 cm)]

% define binary sensor mask with two sensor positions
sensor.mask = zeros(Nx, Ny);
pos1 = 45;                              % [grid points]
pos2 = 65;                              % [grid points]
sensor.mask(pos1, Ny/2) = 1;
sensor.mask(pos2, Ny/2) = 1;

% calculate the distance between the sensor positions
d_cm = (pos2 - pos1)*dx*100;            % [cm]

% set sensor to record to particle velocity
sensor.record = {'u'};

% define source mask
source_mask = ones(Nx, Ny);
source_pos = 35;                        % [grid points]

% set the CFL
cfl = 0.05;

% define the properties of the PML to allow plane wave propagation
pml_alpha = 0;
pml_size  = [30, 2];

% set the input arguments
input_args = {'PlotScale', 'auto', 'PMLSize', pml_size,...
    'PMLAlpha', pml_alpha, 'PlotPML', false};

% =========================================================================
% COMPRESSIONAL PLANE WAVE SIMULATION
% =========================================================================

% define source
source.u_mask = source_mask;
source.ux = zeros(Nx, Ny);
source.ux(source_pos, :) = 1;
source.ux = smooth(kgrid, source.ux, true);
source.ux = 1e-6.*reshape(source.ux, [], 1);

% set end time
t_end = 3.5e-6;

% create the time array
kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression, medium.sound_speed_shear), cfl, t_end);

% run the simulation
sensor_data_comp = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});

% calculate the amplitude spectrum at the two sensor positions
[~, as1] = spect(sensor_data_comp.ux(1, :), 1/kgrid.dt);
[f_comp, as2] = spect(sensor_data_comp.ux(2, :), 1/kgrid.dt);

% calculate the attenuation from the amplitude spectrums
attenuation_comp = -20*log10(as2./as1)./d_cm;

% calculate the corresponding theoretical attenuation in dB/cm
attenuation_th_comp = medium.alpha_coeff_compression .* (f_comp./1e6).^2;

% calculate the maximum supported frequency
f_max_comp = medium.sound_speed_compression / (2*dx);

% find the maximum frequency in the frequency vector
[~, f_max_comp_index] = findClosest(f_comp, f_max_comp);

% =========================================================================
% SHEAR PLANE WAVE SIMULATION
% =========================================================================

% define source
clear source
source.u_mask = source_mask;
source.uy = zeros(Nx, Ny);
source.uy(source_pos, :) = 1;
source.uy = smooth(kgrid, source.uy, true);
source.uy = 1e-6.*reshape(source.uy, [], 1);

% set end time
t_end = 4e-6;

% create the time array
kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression, medium.sound_speed_shear), cfl, t_end);

% run the simulation
sensor_data_shear = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});

% calculate the amplitude at the two sensor positions
[~, as1] = spect(sensor_data_shear.uy(1, :), 1/kgrid.dt);
[f_shear, as2] = spect(sensor_data_shear.uy(2, :), 1/kgrid.dt);

% calculate the attenuation from the amplitude spectrums
attenuation_shear = -20*log10(as2./as1)./d_cm;

% calculate the corresponding theoretical attenuation in dB/cm
attenuation_th_shear = medium.alpha_coeff_shear .* (f_shear./1e6).^2;

% calculate the maximum supported frequency
f_max_shear = medium.sound_speed_shear / (2*dx);

% find the maximum frequency in the frequency vector
[~, f_max_shear_index] = findClosest(f_shear, f_max_shear);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot compressional wave traces
figure;
subplot(4, 1, 1)
t_axis = (0:length(sensor_data_comp.ux)-1)*kgrid.dt*1e6;
plot(t_axis, sensor_data_comp.ux.', 'k-');
axis tight;
xlabel('Time [\mus]');
ylabel('Particle Velocity');
title('Compressional Wave');

% plot compressional wave absorption
subplot(4, 1, 2)
plot(f_comp./1e6, attenuation_comp, 'ko', f_comp./1e6, attenuation_th_comp, 'k-');
set(gca, 'XLim', [0 f_max_comp/1e6], 'YLim', [0 attenuation_th_comp(f_max_comp_index)*1.1]);
box on;
xlabel('Frequency [MHz]');
ylabel('\alpha [dB/cm]');

% plot shear wave traces
subplot(4, 1, 3)
t_axis = (0:length(sensor_data_shear.uy)-1)*kgrid.dt*1e6;
plot(t_axis, sensor_data_shear.uy.', 'k-');
axis tight;
xlabel('Time [\mus]');
ylabel('Particle Velocity');
title('Shear Wave');

% plot shear wave absorption
subplot(4, 1, 4)
plot(f_shear./1e6, attenuation_shear, 'ko', f_shear./1e6, attenuation_th_shear, 'k-');
set(gca, 'XLim', [0 f_max_shear/1e6], 'YLim', [0 attenuation_th_shear(f_max_shear_index)*1.1]);
box on;
xlabel('Frequency [MHz]');
ylabel('\alpha [dB/cm]');

scaleFig(1, 1.5);

repo: Aero-matlab


text:
function [pca,c] = spm_vpca_update (T,S,pca,c,m)
% Update VPCA parameters
% FORMAT [pca,c] = spm_vpca_update (T,S,pca,c,m)
%
% T     [d x N] matrix containing N d-dimensional data vectors
%       The nth data sample, t_n, is nth column of T
% S
% pca   data structure (see eg. spm_vpca.m)
% c     information about single component
% m     cluster number (used for mixtures of VPCA model)
%
% pca,c updated info
%__________________________________________________________________________
% Copyright (C) 2012-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_vpca_update.m 5962 2014-04-17 12:47:43Z spm $

q=c(m).q;
d=pca.d;
N=pca.N;
if isfield(pca,'update_mean')
    update_mean=pca.update_mean;
else
    update_mean=1;
end

% Update latent variables
c(m).sxn2=zeros(q,q);
c(m).Sigma_x=inv(eye(q)+pca.mean_tau*c(m).avg_WtW);
for n=1:N,
    c(m).M_x(:,n)=pca.mean_tau*c(m).Sigma_x*c(m).M_w'*(T(:,n)-c(m).mean_mu);
    c(m).xn2(:,:,n)=c(m).Sigma_x+c(m).M_x(:,n)*c(m).M_x(:,n)';
    c(m).sxn2=c(m).sxn2+S(m,n)*c(m).xn2(:,:,n);
end

% Update factors
c(m).Sigma_w = inv(diag(c(m).mean_alpha)+pca.mean_tau*c(m).sxn2);
for k=1:d,
    obs_err=T(k,:)-c(m).mean_mu(k)*ones(1,N);
    obs_err=S(m,:).*obs_err;
    x_err=c(m).M_x*obs_err';
    mwk=pca.mean_tau*c(m).Sigma_w*x_err;
    c(m).M_w(k,:) = mwk';
end
c(m).avg_WtW=c(m).M_w'*c(m).M_w+d*c(m).Sigma_w;

% Update mean
if update_mean
    c(m).Sigma_mu=1/(pca.beta+sum(S(m,:))*pca.mean_tau);
    t_err=zeros(d,1);
    for n=1:N,
        t_err=t_err+S(m,n)*(T(:,n)-c(m).M_w*c(m).M_x(:,n));
    end
    c(m).mean_mu=pca.mean_tau*c(m).Sigma_mu*t_err;
end

% Update alphas - shrinkage priors on factor columns
pca.qa_alpha=pca.a_alpha+d/2;
for i=1:q,
    c(m).qb_alpha(i)=pca.b_alpha+0.5*c(m).avg_WtW(i,i);
    c(m).mean_alpha(i)=pca.qa_alpha/c(m).qb_alpha(i);
end

c(m).sum_quad_exp=0;
T3=0.5*spm_logdet(c(m).Sigma_x);
quad_term=c(m).mean_mu'*c(m).mean_mu+trace(c(m).Sigma_mu);
trSx=trace(c(m).Sigma_x);
mutW=2*c(m).mean_mu'*c(m).M_w;
for n=1:N,
    term1=T(:,n)'*T(:,n)+quad_term;
    term2=trace(c(m).avg_WtW*c(m).xn2(:,:,n));
    term3=mutW*c(m).M_x(:,n);
    term4=-2*T(:,n)'*c(m).M_w*c(m).M_x(:,n)-2*T(:,n)'*c(m).mean_mu;
    c(m).sum_quad_exp=c(m).sum_quad_exp+S(m,n)*(term1+term2+term3+term4);
    c(m).T1b(n)=-0.5*(c(m).M_x(:,n)'*c(m).M_x(:,n)+trSx);
    T2=-0.5*pca.mean_tau*(term1+term2+term3+term4);
    pca.log_QS(m,n)=c(m).T1b(n)+T2+T3;
end

if pca.M==1
    % Update tau - precision of observation noise
    pca.qa_tau=pca.a_tau+0.5*N*d;
    pca.qb_tau=pca.b_tau+0.5*c(m).sum_quad_exp;
    pca.mean_tau=pca.qa_tau/pca.qb_tau;
end

author: spm


text:
function [nodes2, edges2] = grRemoveNodes(nodes, edges, rmNodes)
%GRREMOVENODES Remove several nodes in a graph.
%
%   usage:
%   [NODES2 EDGES2] = grRemoveNodes(NODES, EDGES, NODES2REMOVE)
%   remove the nodes with indices NODE2REMOVE from array NODES, and also
%   remove edges containing the nodes NODE2REMOVE.
%
%   Example
%     nodes = [...
%         10 10; 20 10; 30 10; ...
%         10 20; 20 20; 30 20];
%     edges = [...
%         1 2; 1 4; 1 5; ...
%         2 3; 2 5; 2 6; ...
%         3 6; 4 5; 5 6];
%     toRemove = [3 4];
%     [nodes2 edges2] = grRemoveNodes(nodes, edges, toRemove);
%     drawGraph(nodes2, edges2);
%     axis equal; axis([0 40 0 30]);
%
%   See also
%     grRemoveEdges
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-08-13
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

%% edges processing

% remove all edges connected to one of the nodes to remove
edges2 = edges(~any(ismember(edges, rmNodes), 2), :);

% change edges information, due to the node index shift
for i = 1:length(rmNodes)
    inds = edges2 > (rmNodes(i) - i + 1);
    edges2(inds) = edges2(inds) - 1;
end

%% nodes processing

% number of nodes
N   = size(nodes, 1);
NR  = length(rmNodes);
N2  = N-NR;

% allocate memory
nodes2 = zeros(N2, 2);

% process the first node
nodes2(1:rmNodes(1)-1,:) = nodes(1:rmNodes(1)-1,:);

for i = 2:NR
    inds = rmNodes(i-1)+1:rmNodes(i)-1;
    if isempty(inds)
        continue;
    end
    nodes2(inds - i + 1, :) = nodes(inds, :);
end

% process the last node
nodes2(rmNodes(NR)-NR+1:N2, :) = nodes(rmNodes(NR)+1:N, :);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function[] = Param_Space(Cat)

rdc= ReasenbergDeclusterClass(Cat, ....
    'AutoShowPlots',false,'DelayProcessing',true,'InteractiveMode',false);

lTaumin     = 1;
lTaumax     = 10;
lP          = .95;
lXk         = .4;
lXmeff      = 3;
lRfact      = 10;
lErr        = 2;
lDerr       = 5;

rdc.taumin  = lTaumin;
rdc.taumax  = lTaumax;
rdc.P       = lP;
rdc.xk      = lXk;
rdc.xmeff   = lXmeff;
rdc.rfact   = lRfact;
rdc.err     = lErr;
rdc.derr    = lDerr;

%% loop for P

resFileOut = 'decResP';
parmFileOut = 'decParmP';
P_Sim = .8:.01:1;
numSim = length(P_Sim);

for simNum = 1:numSim
    rdc.P = P_sim(simNum);
    [declusCat] = rdc.ReasenbergDeclus();
    decResult(simNum) = {declusCat};
    save(resFileOut,'decResult');

    monteParms(simNum) = {[rdc.taumin;rdc.taumax;rdc.P;rdc.xk;rdc.xmeff;rdc.rfact;rdc.err;rdc.derr]};
    save(parmFileOut,'monteParms');
    disp(num2str(simNum));
end

rdc.P = lTauP;

%% loop for Taumin

resFileOut = 'decResTmin';
parmFileOut = 'decParmTmin';
T_Sim = .2:.1:5;
numSim = length(T_Sim);

for simNum = 1:numSim
    rdc.taumin = T_Sim(simNum);
    [declusCat] = rdc.ReasenbergDeclus();
    decResult(simNum) = {declusCat};
    save(resFileOut,'decResult');

    monteParms(simNum) = {[rdc.taumin;rdc.taumax;rdc.P;rdc.xk;rdc.xmeff;rdc.rfact;rdc.err;rdc.derr]};
    save(parmFileOut,'monteParms');
    disp(num2str(simNum));
end

rdc.taumin = lTaumin;

%% loop for Taumax

resFileOut = 'decResTmax';
parmFileOut = 'decParmTmax';
T_Sim = 1:1:20;
numSim = length(T_Sim);

for simNum = 1:numSim
    rdc.taumax = T_Sim(simNum);
    [declusCat] = rdc.ReasenbergDeclus();
    decResult(simNum) = {declusCat};
    save(resFileOut,'decResult');

    monteParms(simNum) = {[rdc.taumin;rdc.taumax;rdc.P;rdc.xk;rdc.xmeff;rdc.rfact;rdc.err;rdc.derr]};
    save(parmFileOut,'monteParms');
    disp(num2str(simNum));
end
rdc.taumax  = lTaumax;

%% loop for rFact

resFileOut = 'decRes_rFact';
parmFileOut = 'decParm_rFact';
R_Sim = 0:1:40;
numSim = length(R_Sim);

for simNum = 1:numSim
    rdc.rfact = R_Sim(simNum);
    [declusCat] = rdc.ReasenbergDeclus();
    decResult(simNum) = {declusCat};
    save(resFileOut,'decResult');

    monteParms(simNum) = {[rdc.taumin;rdc.taumax;rdc.P;rdc.xk;rdc.xmeff;rdc.rfact;rdc.err;rdc.derr]};
    save(parmFileOut,'monteParms');
    disp(num2str(simNum));
end

%% loop for xK

resFileOut = 'decRes_xk';
parmFileOut = 'decParm_xk';
Xk_Sim = 0:.1:1;
numSim = length(Xk_Sim);

for simNum = 1:numSim
    rdc.xk = Xk_Sim(simNum);
    [declusCat] = rdc.ReasenbergDeclus();
    decResult(simNum) = {declusCat};
    save(resFileOut,'decResult');

    monteParms(simNum) = {[rdc.taumin;rdc.taumax;rdc.P;rdc.xk;rdc.xmeff;rdc.rfact;rdc.err;rdc.derr]};
    save(parmFileOut,'monteParms');
    disp(num2str(simNum));
end

author: CelsoReyes


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
%   - data                 Hand, Omega=(0,20)x(0,25),
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             SSD
%   - regularizer          mbElastic
%   - optimizer            Gauss-Newton
% ===============================================================================

close all, help(mfilename)

% load data, initialize image viewer, interpolator, transformation, distance
setup2DhandData
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-2);
trafo('reset','trafo','affine2D');
distance('reset','distance','SSD');
regularizer('reset','regularizer','mbElastic','alpha',1e3,'mu',1,'lambda',0);

PIRobj = @PIRobjFctn;
PIRopt = optPara('GN');

NPIRobj = @NPIRobjFctn;
NPIRopt = optPara('NPIR-GN')

[yc,wc,his] = MLIR(ML,'PIRobj',PIRobj,'PIRopt',PIRopt,'NPIRobj',NPIRobj,...
  'NPIRopt',NPIRopt,'pause',1,'parametric',1);

% ===============================================================================