text:
function [globalOpt, options, optionsDyn, initXOptions] = ...
    hsvargplvmRegressionInitX(globalOpt, options, optionsDyn, inpX, Ytr, stackedOpt)

if nargin < 6 || isempty(stackedOpt)
    stackedOpt = [];
end

% ---- Potential special initialisations for X -----
if ~iscell(globalOpt.initX) && strcmp(globalOpt.initX, 'inputs')
    options = rmfield(options, 'initX');
    for i=1:options.H
        options.initX{i} = inpX;
    end
    optionsDyn.initX = inpX;
    globalOpt.initX = options.initX;
end

% Initialise half of the latent spaces with inputs, half with PCA on outputs
if ~iscell(globalOpt.initX) && strcmp(globalOpt.initX, 'inputsOutputs')
    options = rmfield(options, 'initX');
    if iscell(globalOpt.Q)
        oldQ = globalOpt.Q{end};
    else
        oldQ = globalOpt.Q;
    end
    for i=options.H:-1:floor(options.H/2)+1
        options.initX{i} = inpX;
        Q{i} = size(inpX,2);
    end
    optionsDyn.initX = inpX;

    YtrScaled = scaleData(Ytr{1}, options.scale2var1);
    Xpca  = ppcaEmbed(YtrScaled, oldQ);
    for i=1:floor(options.H/2)
        options.initX{i} = Xpca;
        Q{i} = oldQ;
    end
    options.Q = Q;
    globalOpt.Q = Q;
    globalOpt.initX = options.initX;
end

% Just rewrite all options into a struct of cells
optionsAll = hsvargplvmCreateOptions(Ytr, options, globalOpt);
% Don't mind the following for loop... it just gives the extra possibility
% of initialising the latent space with Bayesian GPLVM or GPLVM (see
% hsvargplvm_init on how to activate this).
initXOptions = cell(1, options.H);
for h=1:options.H
    if strcmp(optionsAll.initX, 'vargplvm') | strcmp(optionsAll.initX, 'fgplvm')
        initXOptions{h}{1} = optionsAll;
        % DOn't allow the D >> N trick for layers > 1
        if h~=1
            if isfield(initXOptions{h}{1}, 'enableDgtN')
                initXOptions{h}{1}.enableDgtN = false;
            end
        end
        initXOptions{h}{1}.latentDim = optionsAll.Q{h};
        initXOptions{h}{1}.numActive = optionsAll.K{h}{1};
        initXOptions{h}{1}.kern = optionsAll.kern{h}{1};
        initXOptions{h}{1}.initX = 'ppca';
        initXOptions{h}{1}.initSNR = 90;
        initXOptions{h}{1}.numActive = 50;
        initXOptions{h}{2} = 160;
        initXOptions{h}{3} = 30;
        if isfield(stackedOpt, 'stackedInitVardistIters'),  initXOptions{h}{2} = stackedOpt.stackedInitVardistIters;   end
        if isfield(stackedOpt, 'stackedInitIters'), initXOptions{h}{3} = stackedOpt.stackedInitIters;   end
        if isfield(stackedOpt, 'stackedInitSNR'), initXOptions{h}{1}.initSNR = stackedOpt.stackedInitSNR; end
        if isfield(stackedOpt, 'stackedInitK'), initXOptions{h}{1}.numActive = stackedOpt.stackedInitK; end
    else
        initXOptions{h} = {};
    end
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
%% LIPSOL TESTING

%% LP1 [-31.4]
clc
%Objective & Constraints
f = -[6 5]';
A = ([1,4; 6,4; 2, -5]);
b = [16;28;6];
lb = [0;0];
ub = [10;10];
%Solve
[x,fval,ef,info] = opti_lipsol(f,A,b,[],[],lb,ub)

%% LP2
clc
%Objective & Constraints
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
Aeq = [1 1 1];
beq = 40;
lb = [0;0;0];
ub = [40;inf;inf];
%Solve
[x,fval,ef,info] = opti_lipsol(f,A,b,Aeq,beq,lb,ub)

%% LP3 (needs lower bounds...)
% clc
% %Objective & Constraints
% f = [8,1]';
% A = [-1,-2;1,-4;3,-1;1,5;-1,1;-1,0;0,-1];
% b = [-4,2,21,39,3,0,0]';
% %Solve
% [x,fval,ef,info] = opti_lipsol(f,A,b)

%% Raw LIPSOL Test
clc
%Problem
A = sparse([ 1     3     4     1     0
             5    -1     1     0    -1
             2     1    -1     0     0]);
b = [8;0;4];
c = [2;1;-2.5;0;0];
lb = [2;0;0;0;0];
ub = [1e32;1e32;9;1e32;1e32];
[xsol,fp,fd,info,msg,times] = lipsol(A,b,c,lb,ub)

%% MPS Files
clc
prob = coinRead('maros-r7.mps');
opts = optiset('solver','lipsol','display','iter','maxiter',1e4);
Opt = opti(prob,opts);
%Solve
[~,f,e,i] = solve(Opt)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function B = LDRtoLDR(A, expA, expB)

Radiance = LDRtoHDR(A, expA);
B = HDRtoLDR(Radiance, expB);

author: qingsenyangit

repo: AHDRNet


text:
function [Iwhite,lilI,lilmask] = maskCropper(I,mask)
%function [Iwhite,lilI,lilmask] = maskCropper(I,mask)
%Given an Image and a Binary Segmentation Mask create a small icon
%showing only the image content inside the segment
%Iwhite is the resulting icon with background colored white
%lilI is the rectangular subimage defined by extent of mask
%lilmask is the small mask
%Tomasz Malisiewicz (tomasz@cmu.edu)

sss = [size(I,1) size(I,2)];

[uu,vv] = find(mask);
minu = min(uu);
maxu = max(uu);
minv = min(vv);
maxv = max(vv);

lilmask = mask(minu:maxu,minv:maxv);
I = I(minu:maxu,minv:maxv,:);

Iwhite = (I.*repmat(lilmask,[1 1 3])) + repmat(~lilmask,[1 1 3]);
lilI = I;

author: quantombone

repo: exemplarsvm


text:
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function sample = SampleMultinomial(probabilities)

dice = rand(1,1);
accumulate = 0;
for i=1:length(probabilities)
    accumulate = accumulate + probabilities(i);
   if accumulate/sum(probabilities) > dice
       break
   end
end
sample = i;

author: anhncs

repo: Probabilistic-Graphical-Models


text:
%%*****************************************************************
%% linsysolvefun: Solve H*x = b
%%
%% x = linsysolvefun(L,b)
%% where L contains the triangular factors of H.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function x = linsysolvefun(L,b)

  x = zeros(size(b));
  for k=1:size(b,2)
     if strcmp(L.matfct_options,'chol')
        x(L.perm,k) = mextriang(L.R, mextriang(L.R,b(L.perm,k),2) ,1);
        %% x(L.perm,k) = L.R \ (b(L.perm,k)' / L.R)';
     elseif strcmp(L.matfct_options,'spchol')
        x(L.perm,k) = mextriangsp(L.Rt,mextriangsp(L.R,b(L.perm,k),2),1);
     elseif strcmp(L.matfct_options,'ldl')
	x(L.p,k) = ((L.D\ (L.L \ b(L.p,k)))' / L.L)';
     elseif strcmp(L.matfct_options,'spldl')
        btmp = b(:,k).*L.s;
        xtmp(L.p,1) = L.Lt\ (L.D\ (L.L \ btmp(L.p)));
	x(:,k) = xtmp.*L.s;
     elseif strcmp(L.matfct_options,'lu')
        x(:,k) = L.U \ (L.L \ b(L.p,k));
     elseif strcmp(L.matfct_options,'splu')
	btmp = b(:,k)./L.s;
        x(L.q,k) = L.U \ (L.L \ (btmp(L.p)));
     end
  end
%%*************************************************************************

author: yu-jiang

repo: radpbook

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


text:
function C = ProdKL(A,B,K,L)
%PRODKL       Matrix product in approximately K-fold precision stored in L results
%
%   C = ProdKL(A,B,K,L)
%
% Input A or B or both may be cell arrays, 1 <= L <= K.
% Output C is cell array iff L>1, default is L=1.
% Simple application of SumKL. All input must be real.
% Relative error of the result is of order  eps^L + eps^K cond(product) , see
%   S.M. Rump: Inversion of extremely ill-conditioned matrices in floating-point,
%      Japan J. Indust. Appl. Math. (JJIAM), 26:249-277, 2009.
%
%Reference implementation! Slow due to interpretation!
%

% written  02/17/08     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, default L=1
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin==3
    L = 1;
  end

  if iscell(A)              % input A cell array
    n = size(A{1},1);
    lenA = length(A);
    AA = zeros(n,lenA*n);
    for i=1:lenA
      AA(:,(i-1)*n+1:i*n) = A{i};
    end
    if iscell(B)            % both A and B cell arrays
      lenB = length(B);
      m = lenA*n;
      AAA = zeros(n,lenB*m);
      BBB = zeros(lenB*m,n);
      for i=1:lenB
        AAA(:,(i-1)*m+1:i*m) = AA;
        BBB((i-1)*m+1:i*m,:) = repmat(B{i},lenA,1);
      end
      A = AAA;
      B = BBB;
    else
      A = AA;
      B = repmat(B,lenA,1);
    end
  else
    n = size(A,1);
  end
  if iscell(B)              % input B cell array
    lenB = length(B);
    A = repmat(A,1,lenB);
    BB = zeros(lenB*n,n);
    for i=1:lenB
      BB((i-1)*n+1:i*n,:) = B{i};
    end
    B = BB;
  end
  if L==1
    C = zeros(n);
  else
    for i=1:L
      C{i} = zeros(n);
    end
  end
  for i=1:n
    for j=1:n
      [x,y] = TwoProduct(A(i,:),B(:,j)');
      res = SumKL([x y],K,L);
      if L==1
        C(i,j) = res;
      else
        for k=1:L
          C{k}(i,j) = res(k);
        end
      end
    end
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function hc = binHistIndep(mu)
% hc = binHistIndep(mu)
% 	Computes expected histogram under independence assumption
%   P(X)=PROD(P(x_i))
%
% Code from the paper: 'Generating spike-trains with specified
% correlations', Macke et al., submitted to Neural Computation
%
% www.kyb.mpg.de/bethgegroup/code/efficientsampling

% generate all possible binary patterns
n = size(mu,1);
c = 0:2^n-1;
pattern = zeros(n,size(c,2));

for i=n:-1:1
    idx = c>=2^(i-1);
    pattern(i,idx)=1;
    c(idx) = c(idx) - 2^(i-1);
end

pattern = flipud(pattern);

% transform to probabilities
mu = mu/2+.5;

% find relevant probabilities for independent model
pMat = (repmat(mu,1,size(pattern,2)).*pattern) + (repmat(1-mu,1,size(pattern,2)).* (~pattern));

% calculate histogram
hc = prod(pMat);

author: Sable

repo: mcbench-benchmarks


text:
function DEMO_GROUP_PEB
% Demonstration routine for empirical Bayes and Bayesian model reduction
%--------------------------------------------------------------------------
% This routine illustrates the use of Bayesian model reduction when
% inverting hierarchical (linear) models - it is essentially a software
% validation demo and proof of concept. It uses a parametric empirical
% Bayesian model (i.e., nested linear models) to eschew local minima issues
% and to assure the Laplace assumption is correct. In brief, the data are
% generated for multiple subjects, under a linear model with subject
% specific parameters at the first level and group specific parameters at
% the second. These model a group effect common to all subjects in a subset
% of parameters and differences in a further subset. In this demo, we
% consider the full hierarchical inversion of a multisubject  study by
% updating the priors at the first level, using the empirical priors from
% the second level. Crucially, this is done during the optimisation  at the
% first level (i.e., after  every iteration - or small number of iterations
% - at the first level.
%
% This provides a generic scheme for the hierarchical inversion of
% nonlinear and possibly dynamic models in which the first level
% optimisation is informed by the sufficient statistics of the second level
% (namely the empirical priors). This should be contrasted with the summary
% statistic approach, in which the second level optimisation, based upon
% the sufficient statistics of the first level (posteriors and priors) are
% computed after convergence of the first level. The results of inversion
% are compared in terms of the second level posteriors (and the second
% level free energy over iterations). Specifically, we compare a gold
% standard (PEB) inversion, with the summary statistic approach to
% empirical Bayes and the hierarchical inversion demonstrated in this
% routine.
%
% The parameterisation of the models uses the format of DCM. This means
% parameters are specified as a structure with key parameters being in the
% fields A, B and C.
%
% See also: spm_dcm_bmr, spm_dcm_peb and spm_dcm_peb_bma
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston, Peter Zeidman
% $Id: DEMO_GROUP_PEB.m 6737 2016-03-03 12:05:51Z karl $

% set up
%==========================================================================
rng('default')

% model space - defined in terms of combinations of some parameters
%--------------------------------------------------------------------------
k     = spm_perm_mtx(3);
for i = 1:8;
    B{i} = k(i,:);
end

% model space
%--------------------------------------------------------------------------
mw  = 3;                              % true model (within)
mx  = 4;                              % true model (between)
Ns  = 16;                             % number of subjects
C   = 32;                             % within:between [co]variance ratio

% create subject-specifc GLM
%==========================================================================

% within subject effects:  condition specific effects 'B' (2 s.d.)
%--------------------------------------------------------------------------
pC          = 1/8;
sd          = sqrt(pC/C);
DCM.Ep.A    = randn(4,1)*sd;
DCM.Ep.B{1} = B{mw}*2*sd;
Np          = spm_length(DCM.Ep);
DCM.M.pE    = spm_zeros(DCM.Ep);
DCM.M.pC    = eye(Np,Np)*pC;

% between subject effects: constant and group difference
%--------------------------------------------------------------------------
X           = [ones(Ns,1) kron([-1;1],ones(Ns/2,1))];
DCM.Ex      = spm_zeros(DCM.Ep);
DCM.Ex.B{1} = B{mx}*2*sd;

% create subject-specifc DCM
%--------------------------------------------------------------------------
Ex    = spm_vec(DCM.Ex);
Ep    = spm_vec(DCM.Ep);
pC    = DCM.M.pC;
Cp    = sd*diag(~~spm_vec(Ep));
Ny    = 16;
for i = 1:Ns

    % generate data
    %----------------------------------------------------------------------
    Pp    = X(i,1)*Ep + X(i,2)*Ex + Cp*randn(Np,1);

    % generate data
    %----------------------------------------------------------------------
    Z{i,i} = randn(Ny,Np);
    y{i,1} = Z{i,i}*Pp + randn(Ny,1)/8;

    % design matrix and data
    %----------------------------------------------------------------------
    GCM{i,1}.xU    = Z{i,i};
    GCM{i,1}.xY.y  = y{i,1};
    GCM{i,1}.xY.X0 = [];

    % likelihood model and priors
    %----------------------------------------------------------------------
    GCM{i,1}.M.IS  = @(P,M,U) U*spm_vec(P);
    GCM{i,1}.M.pE  = DCM.M.pE;
    GCM{i,1}.M.pC  = pC;
    GCM{i,1}.Tp    = Pp;

end

% PEB (GLM) for inversion to provide a reference for BMR
%==========================================================================
Nx    = size(X,2);
Q     = spm_Ce(ones(1,Np));
for i = 1:Np
    Q{i} = kron(eye(Ns,Ns),Q{i})/128;
end
P{1}.X = spm_cat(Z);
P{1}.C = spm_Ce(ones(1,Ns)*Ny);
P{2}.X = kron(X,eye(Np,Np));
P{2}.C = Q;
P{3}.X = kron(zeros(Nx,1),zeros(Np,1));
P{3}.C = kron(eye(Nx,Nx),pC);

% Full hierarchical parametric empirical Bayes inversion
%--------------------------------------------------------------------------
[qP,~,F] = spm_PEB(spm_cat(y),P,1);

% record estimates as a reference
%--------------------------------------------------------------------------
PB.F  = F;
PB.Ep = qP{3}.E;
PB.Cp = qP{3}.C;

% repeat using non-linear empirical Bayes
%==========================================================================

% second level model
%--------------------------------------------------------------------------
M.X   = X;
M.pE  = DCM.M.pE;
M.pC  = DCM.M.pC;

% hierarchical inversion
%--------------------------------------------------------------------------
[gcm,peb,M] = spm_dcm_peb_fit(GCM,M);

% repeated using Bayesian model reduction summary statistic approach
%==========================================================================
GCM  = spm_dcm_fit(GCM);
PEB  = spm_dcm_peb(GCM,M);

% second level parameter estimates
%==========================================================================
spm_figure('GetWin','Figure 1'); clf

% estimated and true second level parameters
%--------------------------------------------------------------------------
subplot(2,2,1), spm_plot_ci(PB.Ep,PB.Cp), hold on, bar([Ep;Ex],1/2), hold off
xlabel('parameters'), ylabel('expectation'),
title('Parametric Bayes','FontSize',16), axis square

subplot(2,2,2), spm_plot_ci(PEB.Ep(:),PEB.Cp), hold on, bar([Ep;Ex],1/2), hold off
xlabel('parameters'), ylabel('expectation'),
title('Sufficient statistics','FontSize',16), axis square

subplot(2,2,3), spm_plot_ci(peb.Ep(:),peb.Cp), hold on, bar([Ep;Ex],1/2), hold off
xlabel('parameters'), ylabel('expectation'),
title('Hierarchical inversion','FontSize',16), axis square

subplot(2,2,4), bar(gcm{1}.FEB)
xlabel('iterations'), ylabel('free energy'),
title('second level free energy','FontSize',16), axis square

return

% Bayesian model reduction with and without hierarchical inversion
%==========================================================================

% define the model space in terms of a matrix
%--------------------------------------------------------------------------
K     = ones(length(B),spm_length(DCM.Ep));
k     = spm_fieldindices(DCM.M.pE,'B');
for i = 1:length(B)
    K(i,k) = spm_vec(B{i})';
end

% defined model in terms of prior covariance
%--------------------------------------------------------------------------
Nm    = size(K,2);
for i = 1:Ns
    for j = 1:Nm
        gcm{i,j}      = gcm{i,1};
        GCM{i,j}      = GCM{i,1};
        gcm{i,j}.M.pC = diag(K(j,:))*M.pC*diag(K(j,:));
        GCM{i,j}.M.pC = diag(K(j,:))*M.pC*diag(K(j,:));
    end
end

rcm   = spm_dcm_bmr(gcm);
RCM   = spm_dcm_bmr(GCM);

% Free energies
%--------------------------------------------------------------------------
for i = 1:Ns
    for j = 1:Nm
        G(i,j,1) = rcm{i,j}.F - rcm{i,1}.F;
        G(i,j,2) = RCM{i,j}.F - RCM{i,1}.F;
    end
end

%  free energy model comparison
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2');clf

p  = spm_softmax(sum(G(:,:,1))'); [m i] = max(p);
subplot(2,2,1), bar(p)
text(i - 1/4,m/2,sprintf('%-2.0f%%',m*100),'Color','w','FontSize',8)
xlabel('model'), ylabel('probability'), title('Iterative PEB','FontSize',16)
axis([0 (length(p) + 1) 0 1]), axis square

p  = spm_softmax(sum(G(:,:,1))'); [m i] = max(p);
subplot(2,2,2), bar(p)
text(i - 1/4,m/2,sprintf('%-2.0f%%',m*100),'Color','w','FontSize',8)
xlabel('model'), ylabel('probability'), title('Bayesian model reduction','FontSize',16)
axis([0 (length(p) + 1) 0 1]), axis square

author: spm


text:
function pgs(hf)
ss2=get(hf,'UserData');
status=ss2(1);
status2=ss2(2);

ssz=get(0,'ScreenSize');

set(status,'FontUnits','points','FontSize',8*ssz(4)/768);
set(status2,'FontUnits','points','FontSize',8*ssz(4)/768);
drawnow;
pagesetupdlg(hf);

set(status,'FontUnits','points','FontSize',10*ssz(4)/768);
set(status2,'FontUnits','points','FontSize',10*ssz(4)/768);


text:
classdef PTKVentilationMask < PTKPlugin
    % PTKVentilationMask. Plugin for segmenting ventilation values from
    %     hyperpolarised gas MRI
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Ventilation <BR>Mask'
        ToolTip = ''
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            original_image = dataset.GetResult('PTKOriginalImage');
            noise_roi = original_image.RawImage(1, 20:220, 10:60); % Kaushik et al
            roi = dataset.GetResult('PTKLungROI');
            std_noise = std(double(noise_roi(:)));
            mean_signal = mean(double(noise_roi(:)));
            threshold = mean_signal + 2*std_noise;

            threshold_image = logical(roi.RawImage >= threshold);
            results = roi.BlankCopy;
            results.ChangeRawImage(threshold_image);

            results.BinaryMorph(@imopen, 6);

            results.ChangeRawImage(3*uint8(results.RawImage));

            results.ImageType = PTKImageType.Colormap;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function BC=betweenness_bin(G)
%BETWEENNESS_BIN    Node betweenness centrality
%
%   BC = betweenness_bin(A);
%
%   Node betweenness centrality is the fraction of all shortest paths in
%   the network that contain a given node. Nodes with high values of
%   betweenness centrality participate in a large number of shortest paths.
%
%   Input:      A,      binary (directed/undirected) connection matrix.
%
%   Output:     BC,     node betweenness centrality vector.
%
%   Note: Betweenness centrality may be normalised to the range [0,1] as
%   BC/[(N-1)(N-2)], where N is the number of nodes in the network.
%
%   Reference: Kintali (2008) arXiv:0809.1906v2 [cs.DS]
%              (generalization to directed and disconnected graphs)
%
%
%   Mika Rubinov, UNSW/U Cambridge, 2007-2012

n=length(G);                %number of nodes
I=eye(n)~=0;                %logical identity matrix
d=1;                     	%path length
NPd=G;                      %number of paths of length |d|
NSPd=NPd;                  	%number of shortest paths of length |d|
NSP=NSPd; NSP(I)=1;        	%number of shortest paths of any length
L=NSPd; L(I)=1;           	%length of shortest paths

%calculate NSP and L
while find(NSPd,1)
    d=d+1;
    NPd=NPd*G;
    NSPd=NPd.*(L==0);
    NSP=NSP+NSPd;
    L=L+d.*(NSPd~=0);
end
L(~L)=inf; L(I)=0;          %L for disconnected vertices is inf
NSP(~NSP)=1;                %NSP for disconnected vertices is 1

Gt=G.';
DP=zeros(n);            	%vertex on vertex dependency
diam=d-1;                  	%graph diameter

%calculate DP
for d=diam:-1:2
    DPd1=(((L==d).*(1+DP)./NSP)*Gt).*((L==(d-1)).*NSP);
    DP=DP + DPd1;           %DPd1: dependencies on vertices |d-1| from source
end

BC=sum(DP,1);               %compute betweenness

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [irf_record,D_record,gamma_record,struct_irf_record,irf_estimates,D_estimates,gamma_estimates,strshocks_record,strshocks_estimates]=...
    panel6irf(y,Xtilde,theta_gibbs,sigma_gibbs,B_gibbs,Xi,It,Bu,IRFperiods,IRFband,IRFt,rho,thetabar,N,n,m,p,T,d,favar)

% create first the cells storing the results
irf_record={};
irf_estimates={};

% deal with shocks in turn (there are N*n shocks with this model)
for ii=1:N*n

   % start iterating
   for jj=1:It-Bu

   % draw theta from its posterior distribution
   theta=theta_gibbs(:,jj,T);

   % draw B from its posterior distribution
   B=reshape(B_gibbs(:,jj),d,d);
   % obtain its Choleski factor: as B is diagonal, it is simply the square roots of its diagonal entries
   cholB=sparse(diag((diag(B).^0.5)));

   % create a matrix of zeros of dimension p*(N*n)
   Ysim=zeros(p,N*n);
   %  step 2: set the value of the last row, column i, equal to 1
   Ysim(p,ii)=1;

      % step 4: for each iteration kk, repeat the algorithm for periods T+1 to T+h
      for kk=1:IRFperiods-1

      % recover theta for period T+kk
      theta=(1-rho)*thetabar+rho*theta+cholB*mvnrnd(zeros(d,1),eye(d))';

      % use the function lagx to obtain a matrix temp, containing the endogenous regressors
      temp=bear.lagx(Ysim,p-1);

      % define the vector X for the period
      Xsim_t=[temp(end,:) zeros(1,m)];

      % obtain Xbar_t et Xtilde_t
      Xbar_t=kron(speye(N*n),Xsim_t);
      Xtilde_t=Xbar_t*Xi;

      % obtain the predicted value for T+kk
      yp=Xtilde_t*theta;

      % concatenate yp at the top of Y
      Ysim=[Ysim;yp'];

      % repeat until values are obtained for T+h
      end

      % record the results from current iteration in cell irf_record
      % loop over variables
      for kk=1:N*n
      % consider column kk of matrix Ysim and trim the (p-1) initial periods: what remains is the series of IRFs for period T to period T+h-1, for variable kk
      temp=Ysim(p:end,kk);
      % record these values in the corresponding matrix of irf_record
      irf_record{kk,ii}(jj,:)=temp';
      end

   % then go for next iteration
   end

% conduct the same process with shocks in other variables
end

% then apply the structural decomposition (if applicable)
% if IRFs have been set to an unrestricted VAR (IRFt=1):
if IRFt==1
% generate empty structural IRFs (just to be consistent with subsequent parts of the code
struct_irf_record=[];
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(irf_record,N*n,IRFperiods,IRFband,IRFt,[],[],favar);
% if IRFs have been set to an SVAR with Choleski identification (IRFt=2):
elseif IRFt==2
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record,D_record,gamma_record]=bear.irfchol(sigma_gibbs(:,:,T),irf_record,It,Bu,IRFperiods,N*n,favar);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,N*n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);
% if IRFs have been set to an SVAR with triangular factorisation (IRFt=3):
elseif IRFt==3
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record,D_record,gamma_record]=bear.irftrig(sigma_gibbs(:,:,T),irf_record,It,Bu,IRFperiods,N*n,favar);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,N*n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);
end

% record the structural decomposition matrices (the ones obtained from the structural decomposition function only apply to the final period, but incoming functions will require records for all the periods)
% also, if a structural identification was implemented, compute structural shocks
strshocks_record={};
strshocks_estimates={};
% in the case of no structural decomposition
if IRFt==1
% run a pseudo Gibbs sampler to obtain records for D and gamma (for the trivial SVAR)
D_record=repmat(reshape(eye(N*n),(N*n)^2,1),[1 It-Bu T]);
gamma_record=sigma_gibbs;
% obtain estimates
D_estimates=repmat(reshape(eye(N*n),(N*n)^2,1),[1 1 T]);
   for ii=1:T
   gamma_estimates(:,1,ii)=quantile(gamma_record(:,:,ii),0.5,2);
   end
elseif IRFt==2
% run the Gibbs sampler
[strshocks_record,D_record,gamma_record]=bear.strshockspan6(theta_gibbs,sigma_gibbs,y,Xtilde,N,n,T,It,Bu,IRFt);
% obtain point estimates and credibility intervals
[strshocks_estimates]=bear.strsestimates(strshocks_record,N*n,T,IRFband);
% reshape
strshocks_estimates=reshape(strshocks_estimates,[n 1 N]);
% obtain estimates
   for ii=1:T
   D_estimates(:,1,ii)=quantile(D_record(:,:,ii),0.5,2);
   end
gamma_estimates=repmat(reshape(eye(N*n),(N*n)^2,1),[1 1 T]);
elseif IRFt==3
% run the Gibbs sampler
[strshocks_record,D_record,gamma_record]=bear.strshockspan6(theta_gibbs,sigma_gibbs,y,Xtilde,N,n,T,It,Bu,IRFt);
% obtain point estimates and credibility intervals
[strshocks_estimates]=bear.strsestimates(strshocks_record,N*n,T,IRFband);
% reshape
strshocks_estimates=reshape(strshocks_estimates,[n 1 N]);
% obtain estimates
   for ii=1:T
   D_estimates(:,1,ii)=quantile(D_record(:,:,ii),0.5,2);
   gamma_estimates(:,1,ii)=quantile(gamma_record(:,:,ii),0.5,2);
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
% eeg_helpsigproc() - Help file for EEGLAB

function noname();

command = { ...
'pophelp(''acsobiro.m'');' ...
'pophelp(''adjustlocs.m'');' ...
'pophelp(''angtimewarp.m'');' ...
'pophelp(''axcopy.m'');' ...
'pophelp(''binica.m'');' ...
'pophelp(''biosig2eeglabevent.m'');' ...
'pophelp(''blockave.m'');' ...
'pophelp(''cart2topo.m'');' ...
'pophelp(''cbar.m'');' ...
'pophelp(''celltomat.m'');' ...
'pophelp(''chancenter.m'');' ...
'pophelp(''changeunits.m'');' ...
'pophelp(''compvar.m'');' ...
'pophelp(''condstat.m'');' ...
'pophelp(''convertlocs.m'');' ...
'pophelp(''copyaxis.m'');' ...
'pophelp(''coregister.m'');' ...
'pophelp(''crossf.m'');' ...
'pophelp(''dftfilt.m'');' ...
'pophelp(''eegfilt.m'');' ...
'pophelp(''eegfiltfft.m'');' ...
'pophelp(''eegplot.m'');' ...
'pophelp(''eegplot2event.m'');' ...
'pophelp(''eegplot2trial.m'');' ...
'pophelp(''eegrej.m'');' ...
'pophelp(''eegthresh.m'');' ...
'pophelp(''entropy_rej.m'');' ...
'pophelp(''env.m'');' ...
'pophelp(''envtopo.m'');' ...
'pophelp(''epoch.m'');' ...
'pophelp(''erpimage.m'');' ...
'pophelp(''eventalign.m'');' ...
'pophelp(''eventlock.m'');' ...
'pophelp(''eyelike.m'');' ...
'pophelp(''fastif.m'');' ...
'pophelp(''floatread.m'');' ...
'pophelp(''floatwrite.m'');' ...
'pophelp(''forcelocs.m'');' ...
'pophelp(''headplot.m'');' ...
'pophelp(''icaact.m'');' ...
'pophelp(''icadefs.m'');' ...
'pophelp(''icaproj.m'');' ...
'pophelp(''icavar.m'');' ...
'pophelp(''imagesctc.m'');' ...
'pophelp(''jader.m'');' ...
'pophelp(''jointprob.m'');' ...
'pophelp(''kurt.m'');' ...
'pophelp(''loadavg.m'');' ...
'pophelp(''loadcnt.m'');' ...
'pophelp(''loaddat.m'');' ...
'pophelp(''loadeeg.m'');' ...
'pophelp(''loadtxt.m'');' ...
'pophelp(''matsel.m'');' ...
'pophelp(''mattocell.m'');' ...
'pophelp(''movav.m'');' ...
'pophelp(''moveaxes.m'');' ...
'pophelp(''nan_mean.m'');' ...
'pophelp(''openbdf.m'');' ...
'pophelp(''parsetxt.m'');' ...
'pophelp(''phasecoher.m'');' ...
'pophelp(''plotchans3d.m'');' ...
'pophelp(''plotdata.m'');' ...
'pophelp(''ploterp.m'');' ...
'pophelp(''plotmesh.m'');' ...
'pophelp(''plotsphere.m'');' ...
'pophelp(''plottopo.m'');' ...
'pophelp(''posact.m'');' ...
'pophelp(''projtopo.m'');' ...
'pophelp(''qqdiagram.m'');' ...
'pophelp(''quantile.m'');' ...
'pophelp(''read_erpss.m'');' ...
'pophelp(''readbdf.m'');' ...
'pophelp(''readedf.m'');' ...
'pophelp(''readeetraklocs.m'');' ...
'pophelp(''readegi.m'');' ...
'pophelp(''readegihdr.m'');' ...
'pophelp(''readelp.m'');' ...
'pophelp(''readlocs.m'');' ...
'pophelp(''readneurodat.m'');' ...
'pophelp(''readneurolocs.m'');' ...
'pophelp(''realproba.m'');' ...
'pophelp(''rejkurt.m'');' ...
'pophelp(''rejstatepoch.m'');' ...
'pophelp(''rejtrend.m'');' ...
'pophelp(''reref.m'');' ...
'pophelp(''rmbase.m'');' ...
'pophelp(''runica.m'');' ...
'pophelp(''runica_ml.m'');' ...
'pophelp(''runica_ml2.m'');' ...
'pophelp(''runica_mlb.m'');' ...
'pophelp(''sbplot.m'');' ...
'pophelp(''sdfopen.m'');' ...
'pophelp(''sdfread.m'');' ...
'pophelp(''shuffle.m'');' ...
'pophelp(''signalstat.m'');' ...
'pophelp(''slider.m'');' ...
'pophelp(''snapread.m'');' ...
'pophelp(''sobi.m'');' ...
'pophelp(''spec.m'');' ...
'pophelp(''spectopo.m'');' ...
'pophelp(''sph2topo.m'');' ...
'pophelp(''spher.m'');' ...
'pophelp(''spherror.m'');' ...
'pophelp(''strmultiline.m'');' ...
'pophelp(''textsc.m'');' ...
'pophelp(''timef.m'');' ...
'pophelp(''timefdetails.m'');' ...
'pophelp(''timewarp.m'');' ...
'pophelp(''timtopo.m'');' ...
'pophelp(''topo2sph.m'');' ...
'pophelp(''topoplot.m'');' ...
'pophelp(''transformcoords.m'');' ...
'pophelp(''trial2eegplot.m'');' ...
'pophelp(''uigetfile2.m'');' ...
'pophelp(''uiputfile2.m'');' ...
'pophelp(''writecnt.m'');' ...
'pophelp(''writelocs.m'');' ...
};
vartext = { ...
'acsobiro.m' ...
'adjustlocs.m' ...
'angtimewarp.m' ...
'axcopy.m' ...
'binica.m' ...
'biosig2eeglabevent.m' ...
'blockave.m' ...
'cart2topo.m' ...
'cbar.m' ...
'celltomat.m' ...
'chancenter.m' ...
'changeunits.m' ...
'compvar.m' ...
'condstat.m' ...
'convertlocs.m' ...
'copyaxis.m' ...
'coregister.m' ...
'crossf.m' ...
'dftfilt.m' ...
'eegfilt.m' ...
'eegfiltfft.m' ...
'eegplot.m' ...
'eegplot2event.m' ...
'eegplot2trial.m' ...
'eegrej.m' ...
'eegthresh.m' ...
'entropy_rej.m' ...
'env.m' ...
'envtopo.m' ...
'epoch.m' ...
'erpimage.m' ...
'eventalign.m' ...
'eventlock.m' ...
'eyelike.m' ...
'fastif.m' ...
'floatread.m' ...
'floatwrite.m' ...
'forcelocs.m' ...
'headplot.m' ...
'icaact.m' ...
'icadefs.m' ...
'icaproj.m' ...
'icavar.m' ...
'imagesctc.m' ...
'jader.m' ...
'jointprob.m' ...
'kurt.m' ...
'loadavg.m' ...
'loadcnt.m' ...
'loaddat.m' ...
'loadeeg.m' ...
'loadtxt.m' ...
'matsel.m' ...
'mattocell.m' ...
'movav.m' ...
'moveaxes.m' ...
'nan_mean.m' ...
'openbdf.m' ...
'parsetxt.m' ...
'phasecoher.m' ...
'plotchans3d.m' ...
'plotdata.m' ...
'ploterp.m' ...
'plotmesh.m' ...
'plotsphere.m' ...
'plottopo.m' ...
'posact.m' ...
'projtopo.m' ...
'qqdiagram.m' ...
'quantile.m' ...
'read_erpss.m' ...
'readbdf.m' ...
'readedf.m' ...
'readeetraklocs.m' ...
'readegi.m' ...
'readegihdr.m' ...
'readelp.m' ...
'readlocs.m' ...
'readneurodat.m' ...
'readneurolocs.m' ...
'realproba.m' ...
'rejkurt.m' ...
'rejstatepoch.m' ...
'rejtrend.m' ...
'reref.m' ...
'rmbase.m' ...
'runica.m' ...
'runica_ml.m' ...
'runica_ml2.m' ...
'runica_mlb.m' ...
'sbplot.m' ...
'sdfopen.m' ...
'sdfread.m' ...
'shuffle.m' ...
'signalstat.m' ...
'slider.m' ...
'snapread.m' ...
'sobi.m' ...
'spec.m' ...
'spectopo.m' ...
'sph2topo.m' ...
'spher.m' ...
'spherror.m' ...
'strmultiline.m' ...
'textsc.m' ...
'timef.m' ...
'timefdetails.m' ...
'timewarp.m' ...
'timtopo.m' ...
'topo2sph.m' ...
'topoplot.m' ...
'transformcoords.m' ...
'trial2eegplot.m' ...
'uigetfile2.m' ...
'uiputfile2.m' ...
'writecnt.m' ...
'writelocs.m' ...
};
textgui( vartext, command,'fontsize', 15, 'fontname', 'times', 'linesperpage', 18, 'title',strvcat( 'Signal processing functions', '(Click on blue text for help)'));
icadefs; set(gcf, 'COLOR', BACKCOLOR);h = findobj('parent', gcf, 'style', 'slider');set(h, 'backgroundcolor', GUIBACKCOLOR);return;

author: PatternRecognition

repo: OpenBMI


text:
%% Piggyback image
%  Changed: Dec 31st, 2011
%
function [I,lim] = piggyback(I,scale)

    if nargin<2; scale = 2; end; % default, piggybacked image twice as big

    if size(size(I),2)>3; [I,lim]=piggyback_multichannel(I,scale); return; end; % check if I is multichannel

    Ip  = zeros(ceil(size(I)*scale));
    lim = bsxfun(@plus, floor(size(I)*(scale-1)/2), [[1 1 1]; size(I)]); % image limits
    Ip(lim(1):lim(2),lim(3):lim(4),lim(5):lim(6)) = I;                   % piggybacked image
    I = Ip;

end

%% Piggyback image
%  Changed: Dec 31st, 2011
%
function [I,lim] = piggyback_multichannel(I,scale)

    if nargin<2; scale = 2; end; % default, piggybacked image twice as big
    nchannels = size(I,4);

    for i=1:nchannels
        Ip(:,:,:,i) = piggyback(I(:,:,:,i),scale);
    end

    imagesize = [size(I,1) size(I,2) size(I,3)];
    lim = bsxfun(@plus, floor(imagesize*(scale-1)/2), [[1 1 1]; imagesize]); % image limits
    I = Ip;

end


text:
function [f, g] = gpsimMapGradFuncWrapper(param, model)

% GPSIMMAPGRADFUNCWRAPPER wraps the log-likelihood function and the gradient function
% together for minimize optimisation.

% FORMAT
% ARG param : a vector of parameters from the model.
% ARG model : the model.
% RETURN f : the log likelihood of the data set.
% RETURN g : gradient function of the model with respect to the parameters.

% SEEALSO : gpsimMapCreate, gpsimMapLogLikelihood, gpsimMapLogLikeGradients
%
% COPYRIGHT : Pei Gao, Magnus Rattray and Neil D. Lawrence, 2008

% SHEFFIELDML

g = zeros(1,length(param));
if isfield(model, 'comp')
    Nrep = length(model.comp);
    for rep=1:Nrep  %Work out likelihood gradient for each replicate
        options = defaultOptions;
        model.comp{rep} = gpsimMapExpandParam(model.comp{rep}, param);
        model.comp{rep} = gpsimMapUpdateF(model.comp{rep}, options);
        ll(rep) = gpsimMapLogLikelihood(model.comp{rep});
        if nargout > 1
          dg{rep} = gpsimMapLogLikeGradients(model.comp{rep});
          g = g - dg{rep};
        end
        fprintf('log-likelihood %2.4f;\t', ll(rep));
    end
else
    Nrep = 1;
    options = defaultOptions;
    model = gpsimMapExpandParam(model, param);
    model = gpsimMapUpdateF(model, options);
    ll = gpsimMapLogLikelihood(model);
    dg = gpsimMapLogLikeGradients(model);
    g = -dg;
    fprintf('log-likelihood %2.4f;\t', ll);
end

fprintf('\n');
f = - sum(ll);

if nargout>1
  g = g';
end

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Example script showcasing the functionality of slocstat
%
% @author B. Schauerte
% @date   2012

%%
% You can call it on a filename
call = 'slocstat.m';
[sloc,stat] = slocstat(call);
fprintf('Total number of lines %d for ''%s''\n',sloc,call);

%%
% You can call it on a function name
call = 'slocstat';
[sloc,stat] = slocstat(call);
fprintf('Total number of lines %d for ''%s''\n',sloc,call);

% ... and of course of any (non-builtin) function in Matlab's search path,
% e.g., "dct"
call = 'dct';
[sloc,stat] = slocstat(call);
fprintf('Total number of lines %d for ''%s''\n',sloc,call);

%%
% You can call it on a folder name and then it also makes sense to use
% printslocstat
call = '.';
[sloc,stat] = slocstat(call);
fprintf('Total number of lines %d for ''%s''\n',sloc,call);
printslocstat(stat);

author: opencobra

repo: cobratoolbox


text:
function varargout = cluster_tool_draw_sphere(varargin)
% CLUSTER_SRC_TOOL_DRAW_SPHERE M-file for cluster_src_tool_draw_sphere.fig
%      CLUSTER_SRC_TOOL_DRAW_SPHERE, by itself, creates a new CLUSTER_SRC_TOOL_DRAW_SPHERE or raises the existing
%      singleton*.
%
%      H = CLUSTER_SRC_TOOL_DRAW_SPHERE returns the handle to a new CLUSTER_SRC_TOOL_DRAW_SPHERE or the handle to
%      the existing singleton*.
%
%      CLUSTER_SRC_TOOL_DRAW_SPHERE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in CLUSTER_SRC_TOOL_DRAW_SPHERE.M with the given input arguments.
%
%      CLUSTER_SRC_TOOL_DRAW_SPHERE('Property','Value',...) creates a new CLUSTER_SRC_TOOL_DRAW_SPHERE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before cluster_src_tool_draw_sphere_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to cluster_src_tool_draw_sphere_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help cluster_src_tool_draw_sphere

% Last Modified by GUIDE v2.5 01-Aug-2006 16:11:27

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @cluster_tool_draw_sphere_OpeningFcn, ...
                   'gui_OutputFcn',  @cluster_tool_draw_sphere_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

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

% --- Executes just before cluster_src_tool_draw_sphere is made visible.
function cluster_tool_draw_sphere_OpeningFcn(hObject, eventdata, handles, varargin)
movegui(hObject,'center');
set(handles.title,'string',cluster_tool_getver);
handles.src='Image';
handles.intype='Manual';
handles.pwd=pwd;handles.fdir=pwd;
for k=1:size(varargin(:),1)
    if strcmp(varargin{k},'cl'),handles.cl=varargin{k+1};
    elseif strcmp(varargin{k},'fdir'),handles.fdir=varargin{k+1};
    elseif strcmp(varargin{k},'img'),handles.img=varargin{k+1};
        set(handles.sourcespace,'string',['Image: ' handles.img]);
    end
end
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to cluster_src_tool_draw_sphere (see VARARGIN)

% Choose default command line output for cluster_src_tool_draw_sphere
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

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

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

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

function radius_Callback(hObject, eventdata, handles)
try handles.r=str2double(get(hObject,'string'));catch set(hObject,'string','');handles.r=[];end
guidata(hObject,handles);
% Hints: get(hObject,'String') returns contents of radius as text
%        str2double(get(hObject,'String')) returns contents of radius as a double

% --- Executes during object creation, after setting all properties.
function radius_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in selector.
function selector_Callback(hObject, eventdata, handles)
if strcmp(handles.src,'Image')
    if isfield(handles,'img'),msgbox('Selected image will be used instead of previously selected image','','warn');end
    handles.pwd=pwd;cd(handles.fdir);
    [FileName,PathName]=uigetfile('*.img','Select image file.');
    cd(handles.pwd);
    if ~PathName,return,end
    handles.fdir=PathName;
    set(handles.sourcespace,'String',['Source: ' PathName FileName]);
    handles.img=[PathName FileName];
else
    handles.pwd=pwd;cd(handles.fdir);
    [FileName,PathName]=uigetfile('*.mat','Select mat file.');
    cd(handles.pwd);
    if ~PathName,return,end
    handles.fdir=PathName;
    load([PathName FileName]);
    try
        handles.cl_space=cl;
        set(handles.sourcespace,'String',['Source: ' PathName FileName]);
        handles.cl_space(1).fname=[PathName FileName];
    catch
        beep,disp('Error: .mat file likely does not contain a cl structure'),return
    end
end
guidata(hObject,handles);

% --- Executes on button press in done.
function done_Callback(hObject, eventdata, handles)
try
    if get(handles.image_src,'Value')
        cl_space=iimg_read_img(handles.img,1);
    else
        cl_space=handles.cl_space(1);
    end
    if ~isfield(cl_space,'M')
        if ~isfield(cl_space,'mat'),beep,disp('Error: no .M or .mat field in source space.'),return,end
        M=cl_space.mat;
    else
        M=cl_space.M;
    end
    if get(handles.Man_in,'Value')
        XYZmm=sphere_3d([handles.x_val handles.y_val handles.z_val]',handles.r,diag(M(1:3,1:3))');
    else
        XYZmm=sphere_3d(handles.sphere_centers,handles.r,diag(M(1:3,1:3))');
    end
    XYZ=mmToVoxel(XYZmm,M,'valid');
    cl_ind=spm_clusters(XYZ);
    for k=1:max(cl_ind)
        cl(k).XYZ=XYZ(:,find(cl_ind==k));
        cl(k).XYZmm=voxelToMm(cl(k).XYZ,M);
        cl(k).numVox=length(find(cl_ind==k));
        cl(k).M=M;
        cl(k).voxSize=diag(cl(k).M(1:3,1:3));
        cl(k).name='';
        cl(k).Z(1:cl(k).numVox)=1;
        cl(k).mm_center=center_of_mass(cl(k).XYZmm,cl(k).Z);
        cl(k).dim=cl_space.dim;
        cl(k).threshold=[1 NaN];
        cl(k).title=['Cluster of ' num2str(cl(k).numVox) ' voxels, created by cluster_tool_draw_sphere GUI. Image space from: ' cl_space.fname];
    end
    if get(handles.domask,'Value')
        if get(handles.image_src,'Value')
            cl=cl_img_intersect(cl,handles.img);
        else
            cl=cl_cl_intersect(cl,handles.cl_space);
        end
    end
    if isfield(handles,'cl')
        if isfield(handles,'cl'),append=cluster_tool_callgui('cluster_tool_append');else handles.append=0;end
        if append
            handles.cl=[handles.cl(:);cl(:)];
        else
            handles.cl=cl;
        end
    else
        handles.cl=cl;
    end
catch beep,disp('Error: Not all required data is specified. Check that you have selected a file, specified a radius, and specified valid sphere centers.'),return
end
guidata(handles.figure1,handles);

% --- Executes on button press in ret.
function ret_Callback(hObject, eventdata, handles)
if isfield(handles,'cl'),cluster_tool('cl',handles.cl,'fdir',handles.fdir);
else cluster_tool('fdir',handles.fdir);end
close(handles.figure1);

function x_Callback(hObject, eventdata, handles)
try handles.x_val=str2double(get(hObject,'string'));catch set(hObject,'string','');handles.x_val=[];end
guidata(hObject,handles);
% Hints: get(hObject,'String') returns contents of x as text
%        str2double(get(hObject,'String')) returns contents of x as a double

% --- Executes during object creation, after setting all properties.
function x_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function y_Callback(hObject, eventdata, handles)
try handles.y_val=str2double(get(hObject,'string'));catch set(hObject,'string','');handles.y_val=[];end
guidata(hObject,handles);

% Hints: get(hObject,'String') returns contents of y as text
%        str2double(get(hObject,'String')) returns contents of y as a double

% --- Executes during object creation, after setting all properties.
function y_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function z_Callback(hObject, eventdata, handles)
try handles.z_val=str2double(get(hObject,'string'));catch set(hObject,'string','');handles.z_val=[];end
guidata(hObject,handles);

% Hints: get(hObject,'String') returns contents of z as text
%        str2double(get(hObject,'String')) returns contents of z as a double

% --- Executes during object creation, after setting all properties.
function z_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function var_name_Callback(hObject, eventdata, handles)
try
    handles.sphere_centers=evalin('base',get(hObject,'String'));
catch beep,disp('Error: Could not get data from workspace variable'),return,end
if size(size(handles.sphere_centers))~=2||~isnumeric(handles.sphere_centers)||(min(size(handles.sphere_centers))~=3&max(size(handles.sphere_centers))~=3)
    beep,disp('Error: Specified variable does not appear to be a valid n x 3 or 3 x n array of center coordinates'),return
end
if size(handles.sphere_centers,1)~=3
    handles.sphere_centers=handles.sphere_centers';
elseif size(handles.sphere_centers,2)==3
    beep,disp('Warning: Assuming that X Y Z values are in ROWS, not COLUMNS (that is, that each COLUMN is a point).'),disp('If this is not the case, transpose the matrix and try again')
end
guidata(hObject,handles);

% Hints: get(hObject,'String') returns contents of var_name as text
%        str2double(get(hObject,'String')) returns contents of var_name as a double

% --- Executes during object creation, after setting all properties.
function var_name_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function source_selector_SelectionChangeFcn(hObject, eventdata, handles)
handles.src=get(hObject,'String');
if strcmp(handles.src,'Image')
    set(handles.selector,'String','Select Image');
    if isfield(handles,'img')
        set(handles.sourcespace,'String',['Source: ' handles.img]);
    else
        set(handles.sourcespace,'string','Source: (none selected)');
    end
else
    set(handles.selector,'String','Select File');
    if isfield(handles,'cl_space')
        set(handles.sourcespace,'String',['Source: ' handles.cl_space.fname]);
    else
        set(handles.sourcespace,'String','Source: (none selected)');
    end
end
guidata(handles.figure1,handles);

% --------------------------------------------------------------------
function center_input_selector_SelectionChangeFcn(hObject, eventdata, handles)
handles.intype=get(hObject,'String');
if strcmp(handles.intype,'Manual')
    set(handles.x,'Visible','on');
    set(handles.y,'Visible','on');
    set(handles.z,'Visible','on');
    set(handles.x_text,'Visible','on');
    set(handles.y_text,'Visible','on');
    set(handles.z_text,'Visible','on');
    set(handles.var_name,'Visible','off');
    set(handles.center_text,'String','Sphere Center:');
    if isfield(handles,'wrkspc_var'),handles=rmfield(handles,'wrkspc_var');end
else
    set(handles.x,'Visible','off');
    set(handles.y,'Visible','off');
    set(handles.z,'Visible','off');
    set(handles.x_text,'Visible','off');
    set(handles.y_text,'Visible','off');
    set(handles.z_text,'Visible','off');
    set(handles.var_name,'Visible','on');
    set(handles.center_text,'String','Variable Name:');
    if isfield(handles,'x_val'),handles=rmfield(handles,'x_val');end
    if isfield(handles,'y_val'),handles=rmfield(handles,'y_val');end
    if isfield(handles,'z_val'),handles=rmfield(handles,'z_val');end
end
guidata(handles.figure1,handles);


text:
function [varargout] = VBA_spm_unvec(vX,varargin)
% Unvectorise a vectorised array - a compiled routine
% FORMAT [varargout] = spm_unvec(vX,varargin)
% varargin  - numeric, cell or structure array
% vX        - spm_vec(X)
%
% i.e. X           = spm_unvec(spm_vec(X),X)
%      [X1,X2,...] = spm_unvec(spm_vec(X1,X2,...),X1,X2,...)
%
% See spm_vec
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_unvec.m 6151 2014-09-03 18:09:33Z guillaume $

%error('spm_unvec.c not compiled - see Makefile')

% deal to multiple outputs if necessary
%--------------------------------------------------------------------------
if nargout > 1
    varargout = VBA_spm_unvec(vX,varargin);
    return
end
if numel(varargin) == 1
    X  = varargin{1};
else
    X  = varargin;
end

% reshape numerical arrays
%--------------------------------------------------------------------------
if isnumeric(X) || islogical(X)
    if ndims(X) > 2
        X(:)  = full(vX);
    else
        X(:)  = vX;
    end
    varargout = {X};
    return
end

% fill in structure arrays
%--------------------------------------------------------------------------
if isstruct(X)
    f = fieldnames(X);
    for i = 1:numel(f)
        c          = {X.(f{i})};
        if isnumeric(c)
            n      = numel(c);
        else
            n      = VBA_spm_length(c);
        end
        c          = VBA_spm_unvec(vX(1:n),c);
        [X.(f{i})] = deal(c{:});
        vX         = vX(n + 1:end);
    end
    varargout      = {X};
    return
end

% fill in cell arrays
%--------------------------------------------------------------------------
if iscell(X)
    for i = 1:numel(X)
        if isnumeric(X{i})
            n      = numel(X{i});
        else
            n      = VBA_spm_length(X{i});
        end
        X{i}  = VBA_spm_unvec(vX(1:n),X{i});
        vX    = vX(n + 1:end);
    end
    varargout = {X};
    return
end

% else
%--------------------------------------------------------------------------
X         = [];
varargout = {X};

author: MBB-team

repo: VBA-toolbox


text:
%% fnc_GetIndex: give the index of the element in the vector that
%%               corresponds to the set of inputs
%
% Usage:
%   i = fnc_GetIndex(C)
%
% Inputs:
%    C                array of input indexes
%
% Output:
%     i               index in the vector that contains all the variances
%
% ------------------------------------------------------------------------
% See also
%
% Author : Flavio Cannavo'
% e-mail: flavio(dot)cannavo(at)gmail(dot)com
% Release: 1.0
% Date   : 29-01-2011
%
% History:
% 1.0  29-01-2011  First release.
%%

function i = fnc_GetIndex(C)

i = sum(2.^(C-1));


text:
function pts = trianglePoints(dt, idx)
    trIndices = dt.ConnectivityList(idx, :);
    pts = zeros(length(trIndices), 2);
    for i = 1:length(trIndices)
        pts(i, :) = dt.Points(trIndices(i), :);
    end
end

author: kavli-ntnu


text:
function T = convert_to_table(CPD, domain, evidence)
% CONVERT_TO_TABLE Convert a discrete CPD to a table
% T = convert_to_table(CPD, domain, evidence)
%
% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.
% The resulting table can easily be converted to a potential.

domain = domain(:);
CPT = CPD_to_CPT(CPD);
odom = domain(~isemptycell(evidence(domain)));
vals = cat(1, evidence{odom});
map = find_equiv_posns(odom, domain);
index = mk_multi_index(length(domain), map, vals);
T = CPT(index{:});
T = T(:);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% summarize conserved clusters as graph

% P1 = vertcat(Pairs_AllClusAllFish{:,1});
% P2 = vertcat(Pairs_AllClusAllFish{:,2});
% Nodes = unique([P1;P2],'rows');
% numNodes = length(Nodes);
%
% [~,IX1] = ismember(P1,Nodes,'rows');
% [~,IX2] = ismember(P2,Nodes,'rows');
%
% G = graph(IX1,IX2); figure;plot(G)
% bins = conncomp(G);
% %%
% thres_fishcount = round(length(range_fish)/2);
% range_bins = unique(bins);
% fishCountForG = zeros(length(range_bins),1);
% for i = range_bins,
%     IX = find(bins==range_bins(i));
%     fishCountForG(i) = length(unique(Nodes(IX,1)));
% end
%
%
% %% Plot anat
% anat_yx_norm = getappdata(hfig,'anat_yx_norm');
%
% figure;
% hold on;
% image(anat_yx_norm)
% axis equal
% axis ij
% axis off
%
% clrmap = hsv(round(max(bins)*1.1));
%
% for i_node = 1:length(Nodes),
%     i_fish = Nodes(i_node,1);
%     clusIX = Nodes(i_node,2);
%
%     GIX = VAR(i_fish).ClusGroup{3}.gIX;
%     IX = find(GIX == clusIX);
%     CIX_abs = VAR(i_fish).ClusGroup{3}.cIX_abs;
%     cIX_abs = CIX_abs(IX);
%     LoadFishDataWithoutTS(hfig,i_fish);
%     CellXYZ_norm = getappdata(hfig,'CellXYZ_norm');
%     xyz_norm = CellXYZ_norm(cIX_abs,:);
%     %    cmap = repmat(clrmap(bins(i_node),:)',1,length(cIX_abs));
%     scatter(xyz_norm(:,2),xyz_norm(:,1),2,clrmap(bins(i_node),:));%,20,clrmap(IX,:))
% end
%
% %%
% anat_yx_norm = getappdata(hfig,'anat_yx_norm');
%
% figure;
% hold on;
% image(anat_yx_norm)
% axis equal
% axis ij
% axis off
%
% clrmap = hsv(round(max(bins)*1.1));
% range_fish_plot = unique(Nodes(:,1))';
% for i_fish = range_fish_plot,
%     [cIX,gIX,M_xyz_norm] = GetDefaultClustersFromLoad(hfig,i_fish);
%     IX_node = find(Nodes(:,1)==i_fish);
%     U = Nodes(IX_node,2);
%     clr = bins(IX_node);
%     for i = 1:length(U),
%         IX = find(gIX==U(i));
%         scatter(M_xyz_norm(IX,2),M_xyz_norm(IX,1),2,clrmap(clr(i),:));%,20,clrmap(IX,:))
%     end
% end

%% JUST SCREEN NORMALLY

% k_consrv_ratio = 1/3;% or just set int number of fish......
k_consrv = 6; % at least in 5/6 more fish within this set

range_fish = 1:18;
for i_fish = range_fish,

    m = TF_fishrange{i_fish};
    m1 = sum(logical(m),2);
%     m1 = sum(m,2);
    U = find(m1>=k_consrv);
    [cIX_in,gIX_in,~,~,~,absIX] = GetDefaultClustersFromLoad(hfig,i_fish);

    cIX = [];gIX = [];
    for i = 1:length(U),
        IX = find(gIX_in==U(i));
        cIX = [cIX;cIX_in(IX)];
        gIX = [gIX;gIX_in(IX)];
    end

    % save cluster
    name = 'D12screen_6/18fish';
    clusgroupID = 10;
    clusIDoverride = 2; % 1
    SaveCluster_Direct(cIX,gIX,absIX,i_fish,name,clusgroupID,clusIDoverride);
end
%%
% SaveVARwithBackup();
% %%
%
% % fish6_clus47:right abducence
% i_fish = 6;
% ix = 47;
%
% m = TF_fishrange{i_fish};
% fIX = find(m(ix,:));
%

author: xiuyechen

repo: FishExplorer


text:
function jacobi_eigenvalue_test02 ( )

%*****************************************************************************80
%
%% JACOBI_EIGENVALUE_TEST02 uses a 4x4 test matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'JACOBI_EIGENVALUE_TEST02\n' );
  fprintf ( 1, '  For a symmetric matrix A,\n' );
  fprintf ( 1, '  JACOBI_EIGENVALUE computes the eigenvalues D\n' );
  fprintf ( 1, '  and eigenvectors V so that A * V = D * V.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  As a "sanity check", input a diagonal matrix.\n' );

  n = 4;

  a = [ ...
      4.0, 0.0, 0.0, 0.0;
      0.0, 1.0, 0.0, 0.0;
      0.0, 0.0, 3.0, 0.0;
      0.0, 0.0, 0.0, 2.0 ];

  r8mat_print ( n, n, a, '  Input matrix A:' );

  it_max = 100;

  [ v, d, it_num, rot_num ] = jacobi_eigenvalue ( n, a, it_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', it_num );
  fprintf ( 1, '  Number of rotations  = %d\n', rot_num );

  r8vec_print ( n, d, '  Eigenvalues D:' );

  r8mat_print ( n, n, v, '  Eigenvector matrix V:' );
%
%  Compute eigentest.
%
  error_frobenius = r8mat_is_eigen_right ( n, n, a, v, d );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm error in eigensystem A*V-D*V = %g\n', ...
    error_frobenius );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%TEST_DESCI_CASSI_BIRD Test decompress snapshot compressive imaging (DeSCI)
%for simulated coded aperture snapshot spectral imaging (CASSI) `bird`
%dataset.
% Reference
%   [1] Y. Liu, X. Yuan, J. Suo, D.J. Brady, and Q. Dai, Rank Minimization
%       for Snapshot Compressive Imaging, IEEE Trans. Pattern Anal. Mach.
%       Intell. (TPAMI), DOI:10.1109/TPAMI.2018.2873587, 2018.
%   [2] X. Yuan, Generalized alternating projection based total variation
%       minimization for compressive sensing, in Proc. IEEE Int. Conf.
%       Image Process. (ICIP), pp. 2539-2543, 2016.
% Dataset
%   `toy` from the CAVE multispectral image database
%     (http://www1.cs.columbia.edu/CAVE/databases/multispectral/).
%   `bird` from the multiframe CASSI system [3] captured in [4].
%     [3] D. Kittle, K. Choi, A. Wagadarikar, and D. J. Brady, Multiframe
%         image estimation for coded aperture snapshot spectral imagers,
%         Appl. Opt., vol. 49, no. 36, pp. 6824-6833, 2010.
%     [4] A. Rajwade, D. Kittle, T.-H. Tsai, D. Brady, and L. Carin, Coded
%         Hyperspectral Imaging and Blind Compressive Sensing, SIAM J. on
%         Imag. Sci., vol. 6, no. 2, pp. 782-812, 2013.
% Contact
%   Xin Yuan, Bell Labs, xyuan@bell-labs.com, initial version Jul 2, 2015.
%   Yang Liu, Tsinghua University, y-liu16@mails.tsinghua.edu.cn, last
%     update Dec 17, 2018.
%   See also GAPDENOISE_CACTI, GAPDENOISE.
clear; clc;
% [0] environment configuration
addpath(genpath('../algorithms')); % algorithms
addpath(genpath('../packages')); % packages
addpath(genpath('../utils')); % utilities

datasetdir = '../dataset'; % dataset
resultdir  = '../results'; % results

% [1] load dataset
para.type   = 'cassi'; % type of dataset, cassi or cacti
para.name   = 'bird'; % name of dataset
para.number = 24; % number of frames in the dataset

datapath = sprintf('%s/%s%d_%s.mat',datasetdir,para.name,...
    para.number,para.type);

if exist(datapath,'file')
    load(datapath); % mask, meas, orig (and para)
else
    error('File %s does not exist, please check dataset directory!',...
        datapath);
end

para.nframe = 1; % number of coded frames in this test
para.MAXB   = 255;

[nrow,ncol,nmask] = size(mask);
nframe = para.nframe; % number of coded frames in this test
MAXB = para.MAXB;

% [2] apply GAP-Denoise for reconstruction
% yall = meas/MAXB;

para.Mfunc  = @(z) A_xy(z,mask);
para.Mtfunc = @(z) At_xy_nonorm(z,mask);

para.Phisum = sum(mask.^2,3);
para.Phisum(para.Phisum==0) = 1;
% common parameters
para.lambda   =     1; % correction coefficiency
para.acc      =     1; % enable GAP-acceleration
para.flag_iqa = false; % disable image quality assessments in iterations

%% [2.1] GAP-TV, ICIP'16
para.denoiser = 'tv'; % TV denoising
  para.maxiter  = 250; % maximum iteration
  para.tvweight = 5; % weight for TV denoising
  para.tviter   = 5; % number of iteration for TV denoising

[vgaptv,psnr_gaptv,ssim_gaptv,tgaptv] = ...
    gapdenoise_cacti(mask,meas,orig,[],para);

fprintf('GAP-%s mean PSNR %2.2f dB, mean SSIM %.4f, total time % 4.1f s.\n',...
    upper(para.denoiser),mean(psnr_gaptv),mean(ssim_gaptv),tgaptv);

%% [2.2] DeSCI (with GAP-TV for initialization), TPAMI'18
para.denoiser = 'wnnm'; % WNNM denoising
  para.wnnm_int_fwise = true; % enable GAP-WNNM integrated (with frame-wise denoising)
    para.blockmatch_period = 20; % period of block matching
  para.sigma   = [12]/MAXB; % noise deviation (to be estimated and adapted)
  para.vrange  = 1; % value range
  para.maxiter = [40];
  para.patchsize = 32; % patch size
  para.iternum = 1; % iteration number in WNNM
  para.enparfor = true; % enable parfor
  if para.enparfor % if parfor is enabled, start parpool in advance
      delete(gcp('nocreate')); % delete current parpool
      mycluster = parcluster('local');
      div = 1;
      while nmask/div > mycluster.NumWorkers
          div = div+1;
      end
      parnum = max(min(ceil(nmask/div),mycluster.NumWorkers),1);
      parpool(mycluster,parnum);
  end

[vdesci,psnr_desci,ssim_desci,tdesci,psnrall] = ...
    gapdenoise_cacti(mask,meas,orig,vgaptv,para); % vgaptv as initialization

  tdesci = tdesci + tgaptv;
  delete(gcp('nocreate')); % delete parpool

fprintf('DeSCI mean PSNR %2.2f dB, mean SSIM %.4f, total time % 4.1f s.\n',...
    mean(psnr_desci),mean(ssim_desci),tdesci);

%% [3] save results as mat file
matdir = [resultdir '/savedmat'];
if ~exist(matdir,'dir')
    mkdir(matdir);
end
save([matdir '/desci_' para.type '_' para.name '.mat']);

repo: DeSCI


text:
function geometry_test174 ( )

%*****************************************************************************80
%
%% TEST174 tests RADEC_TO_XYZ, XYZ_TO_RADEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 6;

  ptest = [ ...
     1.0,  0.0,  0.0; ...
     0.0,  1.0,  0.0; ...
     0.0,  0.0,  1.0; ...
     1.0,  1.0,  1.0; ...
     5.0, -2.0, -1.0; ...
    -2.0, -2.0, -2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST174\n' );
  fprintf ( 1, '  RADEC_TO_XYZ converts XYZ to RADEC coordinates.\n' );
  fprintf ( 1, '  XYZ_TO_RADEC converts RADEC to XYZ coordinates.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          P1          RA     DEC           P2\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    p1(1:3,1) = ptest(1:3,i);

    [ ra, dec ] = xyz_to_radec ( p1 );
    p2 = radec_to_xyz ( ra, dec );

    fprintf ( 1, '  %10f  %10f  %10f  %10f  %10f  %10f  %10f  %10f\n', ...
      p1(1:3,1), ra, dec, p2(1:3,1) );

  end

  return
end


text:
function [source] = fixsource(source)

% FIXSOURCE converts old style source structures into new style source structures
%
% Use as
%   output = fixsource(input)
% where input is a structure representing source data
%
% Typically, old style source structures contain source.avg.XXX or source.trial.XXX.
% Furthermore. old style source structrures do not contain a dimord field.
%
% The new style source structure contains:
%   source.pos       Nx3 list with source positions
%   source.dim       optional, if the list of positions describes a 3D volume
%   source.XXX       the old style subfields in avg/trial
%   source.XXXdimord string how to interpret the respective XXX field:
%
% For example
%   source.pow             = zeros(Npos,Ntime)
%   source.powdimord       = 'pos_time'
%
%   source.mom             = cell(1,Npos)
%   source.mom{1}          = zeros(Nori,Nrpttap)
%   source.momdimord       = '{pos}_ori_rpttap'
%
%   source.leadfield       = cell(1,Npos)
%   source.leadfield{1}    = zeros(Nchan,Nori)
%   source.leadfielddimord = '{pos}_chan_ori'
%
% See also FT_CHECKDATA, FIXVOLUME

npos = size(source.pos, 1);

if ~isfield(source, 'inside')
  source.inside = 1:npos;
  source.outside = [];
end

% old data structures may use latency/frequency instead of time/freq
if isfield(data, 'frequency'),
  data.freq = data.frequency;
  data      = rmfield(data, 'frequency');
end
if isfield(data, 'latency'),
  data.time = data.latency;
  data      = rmfield(data, 'latency');
end

% copy the fields from the average into the main structure
if isfield(source, 'avg')
  avg    = source.avg;
  source = rmfield(source,  'avg');

  fn = fieldnames(avg);
  for i=1:length(fn)
    source.(fn{i}) = avg.(fn{i});
  end
end

% convering the trials has precedence over the average
if isfield(source, 'trial')
  trial  = source.trial;
  source = rmfield(source,  'trial');
  nrpt   = length(trial);
  for i=1:nrpt
    % FIXME
    keyboard
  end
end

hastime = isfield(source, 'time');
hasfreq = isfield(source, 'freq');

if hastime
  ntime = length(source.time);
else
  ntime = 1;
end

if hasfreq
  nfreq = length(source.freq);
else
  nfreq = 1;
end

% FIXME the code bellow should be repatterned
% first get the dimord string
% then determine the data dimensions
% then determine the dimensions according to the dimord
% if consistent, continue
% if inconsistent but can be reshaped, reshape and give warning
% if inconsistent and cannot can be reshaped, give error

% check the dimord and dimensions for each of the fields
fn = fieldnames(source);
for i=1:length(fn)
  element = source.(fn{i});

  if ischar(element) || isstruct(element)
    % the dimord only applies to cell-arrays or numeric arrays
    continue
  end

  if any(strcmp(fn{i}, {'inside', 'outside'}))
    % these don't need a dimord
    continue
  end

  % determine the dimensions of this data element
  if iscell(element)
    dim = numel(element);
    dim = [dim size(element{source.inside(1)})];
  else
    dim = size(element);
  end

  if isfield(source, [fn{i} 'dimord'])
    % check whether the specific dimord applies
    dimord = source.([fn{i} 'dimord']);

    dimtok = tokenize(dimord, '_');
    descr  = [];
    for j=1:length(dimtok)
      if strcmp(dimtok{j}, 'pos')
        descr = [descr size(source.pos,1)];
      elseif strcmp(dimtok{j}, '{pos}')
        descr = [descr size(source.pos,1)];
      else
        descr = [descr numel(source.(dimtok{j}))];
      end
    end

    if prod(dim)~=prod(descr)
      ft_error('the dimensions of the source data are not consistent with the dimord (%s)', fn{i})
    end

    if iscell(element)
      element = reshape(element, [descr(1) 1]);
      for j=1:numel(element)
        element{j} = reshape(element{j}, [descr(2:end) 1]);
      end
    else
      element = reshape(element, [descr 1]);
    end

    % put the reshaped data element back into the structure
    source.(fn{i}) = element;

  elseif isfield(source, 'dimord')
    % check whether the general dimord applies to this element
    dimord = source.dimord;

    dimtok = tokenize(dimord, '_');
    descr  = [];
    for j=1:length(dimtok)
      if strcmp(dimtok{j}, 'pos')
        descr = [descr size(source.pos,1)];
      elseif strcmp(dimtok{j}, '{pos}')
        descr = [descr size(source.pos,1)];
      else
        descr = [descr numel(source.(dimtok{j}))];
      end
    end

    if prod(dim)==prod(descr)
      % the dimensions of the source data are consistent with the dimord
      if iscell(element)
        element = reshape(element, [descr(1) 1]);
        for j=1:numel(element)
          element{j} = reshape(element{j}, [descr(2:end) 1]);
        end
      else
        element = reshape(element, [descr 1]);
      end

      % create a dimord for this data element
      source.([fn{i} 'dimord']) = dimord;
      % put the reshaped data element back into the structure
      source.(fn{i}) = element;
    end

  else
    % try to create a dimord on the fly
    if     ~hasfreq && ~hastime && prod(dim)==npos
      if iscell(element)
        dimord = '{pos}';
      else
        dimord = 'pos';
      end
      descr  = npos;
    elseif  hasfreq && ~hastime && prod(dim)==prod([npos nfreq])
      if iscell(element)
        dimord = '{pos}_freq';
      else
        dimord = 'pos_freq';
      end
      descr  = [npos nfreq];
    elseif  hasfreq &&  hastime && prod(dim)==prod([npos nfreq ntime])
      if iscell(element)
        dimord = '{pos}_freq_time';
      else
        dimord = 'pos_freq_time';
      end
      descr  = [npos nfreq ntime];
    elseif ~hasfreq &&  hastime && prod(dim)==prod([npos ntime])
      if iscell(element)
        dimord = '{pos}_time';
      else
        dimord = 'pos_time';
      end
      descr  = [npos ntime];
    elseif prod(dim)==npos
      % although the source data has time and frequency, this element might be pos only
      if iscell(element)
        dimord = '{pos}';
      else
        dimord = 'pos';
      end
      descr  = npos;
    else
      dimord = '';
      descr  = [];
    end

    if ~isempty(dimord)
      % the dimensions of the source data are consistent with the dimord
      if iscell(element)
        element = reshape(element, [descr(1) 1]);
        for j=1:numel(element)
          element{j} = reshape(element{j}, [descr(2:end) 1]);
        end
      else
        element = reshape(element, [descr 1]);
      end

      % create a dimord for this data element
      source.([fn{i} 'dimord']) = dimord;
      % put the reshaped data element back into the structure
      source.(fn{i}) = element;
    end

  end % if isfield dimord
end % for each of the fields

if isfield(source, 'dimord')
  % remove the general dimord field, there should be one for each data element
  source = rmfield(source, 'dimord');
end

%if isvolume
%  % ensure consistent dimensions of the volumetric data
%  % reshape each of the volumes that is found into a 3D array
%  param = parameterselection('all', data);
%  dim   = data.dim;
%  for i=1:length(param)
%    tmp  = getsubfield(data, param{i});
%    tmp  = reshape(tmp, dim);
%    data = setsubfield(data, param{i}, tmp);
%  end
%end

% if issource || isvolume,
%   % the following section is to make a dimord-consistent representation of
%   % volume and source data, taking trials, time and frequency into account
%   if isequal(hasdimord, 'yes') && (~isfield(data, 'dimord') || ~strcmp(data.dimord, sourcedimord))
%
%     % determine the size of the data
%     if isfield(data, 'dimord'),
%       dimtok = tokenize(data.dimord, '_');
%       if ~isempty(strmatch('time', dimtok)), Ntime = length(data.time); else Ntime = 1; end
%       if ~isempty(strmatch('freq', dimtok)), Nfreq = length(data.freq); else Nfreq = 1; end
%     else
%       Nfreq = 1;
%       Ntime = 1;
%     end
%
%     % convert old style source representation into new style
%     if isfield(data, 'avg') && isfield(data.avg, 'mom') && isfield(data, 'freq') && strcmp(sourcedimord, 'rpt_pos'),
%       % frequency domain source representation convert to single trial power
%       Npos   = size(data.pos,1);
%       Nrpt   = length(data.cumtapcnt);
%       tmpmom = zeros(Npos, size(data.avg.mom{data.inside(1)},2));
%       tmpmom(data.inside,:) = cat(1,data.avg.mom{data.inside});
%       tmppow = zeros(Npos, Nrpt);
%       tapcnt = [0;cumsum(data.cumtapcnt)];
%       for k = 1:Nrpt
%         Ntap = tapcnt(k+1)-tapcnt(k);
%         tmppow(data.inside,k) = sum(abs(tmpmom(data.inside,(tapcnt(k)+1):tapcnt(k+1))).^2,2)./Ntap;
%       end
%       data.pow = tmppow';
%       data     = rmfield(data, 'avg');
%       if strcmp(inside, 'logical'),
%         data     = fixinside(data, 'logical');
%         data.inside = repmat(data.inside(:)',[Nrpt 1]);
%       end
%     elseif isfield(data, 'avg') && isfield(data.avg, 'mom') && isfield(data, 'freq') && strcmp(sourcedimord, 'rpttap_pos'),
%       % frequency domain source representation convert to single taper fourier coefficients
%       Npos   = size(data.pos,1);
%       Nrpt   = sum(data.cumtapcnt);
%       data.fourierspctrm = complex(zeros(Nrpt, Npos), zeros(Nrpt, Npos));
%       data.fourierspctrm(:, data.inside) = transpose(cat(1, data.avg.mom{data.inside}));
%       data   = rmfield(data, 'avg');
%     elseif isfield(data, 'avg') && isfield(data.avg, 'mom') && isfield(data, 'time') && strcmp(sourcedimord, 'pos_time'),
%       Npos   = size(data.pos,1);
%       Nrpt   = 1;
%       tmpmom = zeros(Npos, size(data.avg.mom{data.inside(1)},2));
%       tmpmom(data.inside,:) = cat(1,data.avg.mom{data.inside});
%       data.mom = tmpmom;
%       if isfield(data.avg, 'noise'),
%         tmpnoise = data.avg.noise(:);
%         data.noise = tmpnoise(:,ones(1,size(tmpmom,2)));
%       end
%       data = rmfield(data, 'avg');
%       Ntime = length(data.time);
%     elseif isfield(data, 'trial') && isfield(data.trial(1), 'mom') && isfield(data, 'time') && strcmp(sourcedimord, 'rpt_pos_time'),
%       Npos   = size(data.pos,1);
%       Nrpt   = length(data.trial);
%       Ntime  = length(data.time);
%       tmpmom = zeros(Nrpt, Npos, Ntime);
%       for k = 1:Nrpt
%         tmpmom(k,data.inside,:) = cat(1,data.trial(k).mom{data.inside});
%       end
%       data     = rmfield(data, 'trial');
%       data.mom = tmpmom;
%     elseif isfield(data, 'trial') && isstruct(data.trial)
%       Nrpt = length(data.trial);
%     else
%       Nrpt = 1;
%     end
%
%     % start with an initial specification of the dimord and dim
%     if (~isfield(data, 'dim') || ~isfield(data, 'dimord'))
%       if issource
%         % at least it should have a Nx3 pos
%         data.dim    = size(data.pos, 1);
%         data.dimord = 'pos';
%       elseif isvolume
%         % at least it should have a 1x3 dim
%         data.dim    = data.dim;
%         data.dimord = 'dim1_dim2_dim3';
%       end
%     end
%
%     % add the additional dimensions
%     if Nfreq>1
%       data.dimord = [data.dimord '_freq'];
%       data.dim    = [data.dim     Nfreq];
%     end
%     if Ntime>1
%       data.dimord = [data.dimord '_time'];
%       data.dim    = [data.dim     Ntime];
%     end
%     if Nrpt>1 && strcmp(sourcedimord, 'rpt_pos'),
%       data.dimord = ['rpt_' data.dimord];
%       data.dim    = [Nrpt   data.dim ];
%     elseif Nrpt>1 && strcmp(sourcedimord, 'rpttap_pos'),
%       data.dimord = ['rpttap_' data.dimord];
%       data.dim    = [Nrpt   data.dim ];
%     end
%
%     % the nested trial structure is not compatible with dimord
%     if isfield(data, 'trial') && isstruct(data.trial)
%       param = fieldnames(data.trial);
%       for i=1:length(param)
%         if isa(data.trial(1).(param{i}), 'cell')
%           concat = cell(data.dim(1), prod(data.dim(2:end)));
%         else
%           concat = zeros(data.dim(1), prod(data.dim(2:end)));
%         end
%         for j=1:length(data.trial)
%           tmp = data.trial(j).(param{i});
%           concat(j,:) = tmp(:);
%         end % for each trial
%         data.trial = rmfield(data.trial, param{i});
%         data.(param{i}) = reshape(concat, data.dim);
%       end % for each param
%       data = rmfield(data, 'trial');
%     end
%   end
%
%   % ensure consistent dimensions of the source reconstructed data
%   % reshape each of the source reconstructed parameters
%   if issource && isfield(data, 'dim') && prod(data.dim)==size(data.pos,1)
%     % volumetric represented source structure
%     dim = [prod(data.dim) 1];
%   elseif issource && any(~cellfun('isempty',strfind(fieldnames(data), 'dimord')))
%     % sparsely represented source structure new style
%     dim = [size(data.pos,1) 1];
%   elseif isfield(data, 'dim')
%     dim = [data.dim 1];
%   elseif isfield(data, 'dimord')
%     %HACK
%     dimtok = tokenize(data.dimord, '_');
%     for i=1:length(dimtok)
%       if strcmp(dimtok(i), 'pos')
%         dim(1,i) = size(getsubfield(data,dimtok{i}),1);
%       elseif strcmp(dimtok(i), 'rpt')
%         dim(1,i) = nan;
%       else
%         dim(1,i) = length(getsubfield(data,dimtok{i}));
%       end
%     end
%     i = find(isnan(dim));
%     if ~isempty(i)
%       n = fieldnames(data);
%       for ii=1:length(n)
%         numels(1,ii) = numel(getfield(data,n{ii}));
%       end
%       nrpt = numels./prod(dim(setdiff(1:length(dim),i)));
%       nrpt = nrpt(nrpt==round(nrpt));
%       dim(i) = max(nrpt);
%     end
%     if numel(dim)==1, dim(1,2) = 1; end
%   end
%
% end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cal] = eV2cal(eV)
% Convert energy or work from electron volts to calories.
% Chad A. Greene 2012
cal = eV*3.8267347377e-20;


text:
function dict_dist = find_dups

%% generate splits
[image_ids, class_ids, super_class_ids, path_list] = ...
    textread('/cvgl/group/Ebay_Dataset/Ebay_info.txt', '%d %d %d %s', ...
    'headerlines', 1);

% 1. Do one scan and find out how many classes each super class has
dict_dist = containers.Map;

for i = 1:length(image_ids)
    if mod(i,1000)==0
        fprintf('%d/%d\n', i, length(image_ids));
    end
    this_filename = path_list{i};
    this_image_idx  = image_ids(i);
    this_class_idx = class_ids(i);
    this_super_class_idx = super_class_ids(i);

    C = strsplit(this_filename, '/');
    this_name = C{2};

    if isKey(dict_dist, this_name)
        fprintf('[dup] %s %s\n', this_filename, dict_dist(this_name));
    else
        dict_dist(this_name) = this_filename;
    end
end

% train_images = {};
% val_images = {};
% dict = containers.Map('keytype', 'double', 'valuetype', 'any');
%
% for i = 1:length(image_ids)
%     this_filename = path_list{i};
%     this_image_idx  = image_ids{i};
%     this_class_idx = class_ids{i};
%     this_super_class_idx = super_class_ids{i};
%
%     fprintf('%d/%d, classid=%d, filename= %s\n', i, length(image_ids), ...
%         this_class_idx, this_filename);
%
%
%

author: rksltnl


text:
 function [smap, x, y, z] = ir_mri_sensemap_sim(varargin)
%function [smap, x. y, z] = ir_mri_sensemap_sim(varargin)
%|
%| Simulate 2D or 3D sensitivity maps for sensitivity-encoded MRI
%| based on grivich:00:tmf doi:10.1119/1.19461
%| This code makes maps for multiple coils, but does not model coupling
%| between coils so most likely it is an approximation at best.
%|
%| option
%|	nx, ny, nz		image size (default: [64 64 1])
%|	dx, dy, dz		pixel/voxel dimensions (default: [3 3 3])
%|	ncoil			# of coils total (default: 4)
%|	nring			# of rings of coils (default: 1)
%|	rcoil			coil radius (default: 100mm)
%|	dz_coil			ring spacing in z.  (def: nz*dz/nring)
%|				(3D geometry is a cylinder)
%|	coil_distance		distance of coil center from isocenter for
%|				central ring of coils as a multiple of FOVx,
%|				where FOVx=nx*dx (default: 1.2)
%|	orbit			default: 360
%|	scale			'' (default)
%|				'ssos_center' : make SSoS of center = 1
%|
%| out
%|	smap	[nx ny nz ncoil]	simulated sensitivity maps (complex!)
%|
%| all length parameters must have same units (e.g., mm or cm)
%|
%| Copyright 2005-6-20, Jeff Fessler and Amanda Funai, University of Michigan
%| 2014-08-19 JF more testing, verifying phase is correct
%| 2014-09-09 modified for 3D by Mai Le
%| 2016-05-03 JF fixes

if nargin < 1, ir_usage, end

if nargin == 1 && streq(varargin{1}, 'test', 4)
	ir_mri_sensemap_sim_test(varargin{1})
	return
end

arg.nx = 64;
arg.ny = [];
arg.nz = 1; % 2D
arg.dx = 3; % pixel size in mm
arg.dy = [];
arg.dz = [];
arg.ncoil = 4; % # of coils
arg.nring = 1;
arg.rcoil = 100; % coil radius in mm
arg.orbit = 360;
arg.orbit_start = 0; % can be [nring] to give each ring an offset [degrees]
arg.dz_coil = [];
arg.coil_distance = 1.2; % multiplies fov/2
arg.scale = '';
arg.chat = nargout == 0;

arg = vararg_pair(arg, varargin);

if isempty(arg.dy), arg.dy = arg.dx; end
if isempty(arg.dz), arg.dz = arg.dx; end
if isempty(arg.ny), arg.ny = arg.nx; end
if isempty(arg.rcoil), arg.rcoil = arg.dx * arg.nx / 2 * 0.50; end
if isempty(arg.dz_coil), arg.dz_coil = arg.dz * arg.nz / arg.nring; end

coils_per_ring = round(arg.ncoil / arg.nring);
if arg.nring * coils_per_ring ~= arg.ncoil
	fail('nring must be divisor of ncoil')
end

[ring_smap, x, y, z] = ir_mri_sensemap_sim_do(...
	arg.nx, arg.ny, arg.nz, ...
	arg.dx, arg.dy, arg.dz, ...
	arg.ncoil, coils_per_ring, arg.rcoil, arg.dz_coil, ...
	arg.orbit, arg.orbit_start, arg.coil_distance, arg.chat);

if arg.nz == 1
	smap = reshape(ring_smap, [arg.nx arg.ny arg.ncoil]);
	scale_center = 1 / sqrt(sum(abs(smap(end/2,end/2,:).^2)));
else
	smap = reshape(ring_smap, [arg.nx arg.ny arg.nz arg.ncoil]);
	scale_center = 1 / sqrt(sum(abs(smap(end/2,end/2,end/2,:).^2)));
end

switch arg.scale
case ''
case 'ssos_center'
	smap = smap * scale_center;
otherwise
	fail('unknown scale method "%s"', arg.scale)
end

% ir_mri_sensemap_sim_do()
function [smap x y z] = ir_mri_sensemap_sim_do(nx, ny, nz, ...
		dx, dy, dz, ncoil, ncoilpr, rcoil, dz_coil, ...
		orbit, orbit_start, coil_distance, chat)

nring = ncoil / ncoilpr;
rlist = rcoil * ones(ncoilpr,nring,'single'); % coil radii

plist = zeros(ncoilpr,nring,3,'single'); % position of coil center [x y z]
nlist = zeros(ncoilpr,nring,3,'single'); % normal vector (inward) from coil center
olist = zeros(ncoilpr,nring,3,'single'); % unit vector orthogonal to normal vector in x-y
ulist = zeros(ncoilpr,nring,3,'single'); % upward vector

if numel(orbit_start) == 1
	orbit_start = repmat(orbit_start, nring);
end

% cylindrical coil configuration, like abdominal coils
alist = deg2rad(orbit) * [0:(ncoilpr-1)] / ncoilpr; % coil angles [radians]
z_ring = ([1:nring]-(nring+1)/2) * dz_coil;
for ir = 1:nring
	for ic = 1:ncoilpr
		phi = alist(ic) + deg2rad(orbit_start(ir));
		Rad = max(nx/2 * dx, ny/2 * dy) * coil_distance;
		plist(ic,ir,:) = [Rad * [cos(phi) sin(phi)] z_ring(ir)];
		nlist(ic,ir,:) = -[cos(phi) sin(phi) 0*z_ring(ir)]; % cylinder
		olist(ic,ir,:) = [-sin(phi) cos(phi) 0];
		ulist(ic,ir,:) = [0 0 1];
	end
end

% object coordinates
x = ([1:nx] - (nx+1)/2) * dx;
y = ([1:ny] - (ny+1)/2) * dy;
z = ([1:nz] - (nz+1)/2) * dz;
[xx,yy,zz] = ndgrid(x,y,z);

smap = zeros(nx, ny, nz, ncoilpr, nring, 'single');
for ir = 1:nring
	for ic=1:ncoilpr
		% rotate coordinates to correspond to coil orientation
		zr =	(xx - plist(ic,ir,1)) .* nlist(ic,ir,1) + ...
			(yy - plist(ic,ir,2)) .* nlist(ic,ir,2) + ...
			(zz - plist(ic,ir,3)) .* nlist(ic,ir,3);
		xr =	xx .* nlist(ic,ir,2) - yy .* nlist(ic,ir,1);
		yr = zz - plist(ic,ir,3); % translate along object z axis

		if 0 % see coordinates
			im plc 1 2
			im(1, x, y, xr), xlabel x, ylabel y
			im(2, x, y, zr)
			keyboard
		end

		% compute sensitivity vectors in coil coordinates
		[sx,sy,sz] = ir_mri_smap1(xr, yr, zr, rlist(ic,ir));

		% coil response depends on tranverse magnetization only?
		% todo: unsure if this should depend on sy and ulist in 3D
		bx = sz * nlist(ic,ir,1) + sx * olist(ic,ir,1);
		by = sz * nlist(ic,ir,2) + sx * olist(ic,ir,2);
	%	bz = sz * nlist(ic,ir,3) + sx * olist(ic,ir,3);
		smap(:,:,:,ic,ir) = bx + 1i * by;

		if 0 && chat && nz == 1 && im % see field components
			im plc 2 2
			im(1, x, y, sx), cbar
			im(2, x, y, sy), cbar
			im(3, x, y, sz), cbar
			im subplot 4
			tmp = sqrt(sx.^2 + sz.^2);
			quiver(x, y, (sx./tmp)', (sz./tmp)', 0), axis square
			prompt
		end

		if 0 % see final field components vs phase
			im subplot 4
			bb = sqrt(bx.^2 + by.^2);
			quiver(x, y, (bx./bb)', (by./bb)', 0), axis square
			im(2, x, y, angle(smap(:,:,ic))), cbar
			keyboard
		end
	end
end
smap = smap * rlist(1) / (2*pi); % trick: scale so maximum is near unity

if chat && im % show smap and array geometry in z=0 plane
	if nz == 1
		ir_mri_sensemap_sim_show2(smap, x, y, dx, dy, nlist, plist, rlist)
	else
		ir_mri_sensemap_sim_show3(smap, x, y, z, dx, dy, dz, ...
		nlist, plist, rlist, olist, ulist, nring, ncoilpr, rcoil)
	end
end

% ir_mri_sensemap_sim_show3()
% shows coil geometry but not the 3D smap
function ir_mri_sensemap_sim_show3(smap, x, y, z, dx, dy, dz, ...
	nlist, plist, rlist, olist, ulist, nring, ncoilpr, rcoil)

pcolor = {'c', 'g', 'r'};
pcolor = @(i) pcolor{1+rem(i,3)};
clf, ir_plot3_cube(x,y,z)
xlabel x, ylabel y, zlabel z
hold on
plot3(plist(:,:,1), plist(:,:,2), plist(:,:,3), 'bo') % coil centers
if 1 % coil normals
	tmp1 = reshape(plist, [], 3);
	tmp2 = reshape(nlist, [], 3);
	quiver3(tmp1(:,1), tmp1(:,2), tmp1(:,3), ...
		tmp2(:,1), tmp2(:,2), tmp2(:,3), 0.2)
end
if 1 % coils
	for ir = 1:nring
	for ic = 1:ncoilpr
		tmp = linspace(0, 2*pi, 50)';
		tmp = cos(tmp) * squeeze(olist(ic,ir,:))' + ...
			sin(tmp) * squeeze(ulist(ic,ir,:))';
		tmp = repmat(squeeze(plist(ic,ir,:))', ...
			[nrow(tmp) 1]) + rcoil * tmp;
	%	plot3(tmp(:,1), tmp(:,2), tmp(:,3), 'g-')
		patch(tmp(:,1), tmp(:,2), tmp(:,3), pcolor(ir), ...
		'edgecolor', 'none', 'facealpha', 0.5)
	end
	end
end
hold off
axis equal
% end

function ir_plot3_cube(x,y,z)
x1 = x(1);
x2 = x(end);
y1 = y(1);
y2 = y(end);
z1 = z(1);
z2 = z(end);
x = [x1 x2 x2 x1 x1 x1 x2 x2 x1 x1];
y = [y1 y1 y2 y2 y1 y1 y1 y2 y2 y1];
z = [z1 z1 z1 z1 z1 z2 z2 z2 z2 z2];
plot3(x,y,z)

% ir_mri_sensemap_sim_show2()
function ir_mri_sensemap_sim_show2(smap, x, y, dx, dy, nlist, plist, rlist)
switch ndims(smap)
case 3
	[nx ny ncoil] = size(smap);
case 4
	[nx ny nz ncoil] = size(smap);
otherwise
	fail('unknown ndims(smap) = %d', ndims(smap))
end
nshow = min(max(ncoil,2), 4);
im('plc', 3, nshow)
for ii=1:min(ncoil,nshow)
	clim = [0 max(abs(smap(:)))];
	tmp = smap(:,:,ii);
	im(ii, x, y, abs(tmp), clim, 'Magnitude'), cbar
%	xmax = max(max(abs(x)), max(plist(:,1)));
%	ymax = max(max(abs(y)), max(plist(:,2)));
	xmax = max([max(abs(x)) max(abs(y)) max(col(plist(:,[1 2])))]);
%	axis([-xmax xmax -ymax ymax]*1.05)
	axis(xmax * [-1 1 -1 1] * 1.1)
	xtick([-1 0 1] * nx/2 * dx)
	ytick([-1 0 1] * ny/2 * dy)

	hold on
	plot(0,0,'.', plist(:,1), plist(:,2), 'bo')
	xdir = nlist(ii,2);
	ydir = nlist(ii,1);
	r = rlist(ii);
	plot(plist(ii,1)+r*xdir*[-1 1], plist(ii,2)+r*ydir*[1 -1], 'b-')
	hold off

%	ph = ir_unwrap(angle(tmp)); % trick: unwrap phase for pretty display
	ph = angle(tmp); % show raw phase (understandable with hsv colormap)
	im(ii+nshow, 'hsv', x, y, ph, [-pi pi], 'Phase'), cbar
	axis(xmax*[-1 1 -1 1]*1.1)
	xtick([-1 0 1] * nx/2 * dx)
	ytick([-1 0 1] * ny/2 * dy)
end

ssos = sqrt(sum(abs(smap).^2, ndims(smap)));
minmax(ssos)
ssos = ssos / ssos(end/2,end/2);
im(2*nshow+1, x, y, ssos, 'SSoS (normalized)'), cbar
xtick([-1 0 1] * nx/2 * dx)
ytick([-1 0 1] * ny/2 * dy)

if ncoil == 1
	subplot(122)
	bx = real(smap);
	by = imag(smap);
	quiver(x, y, bx', by'), title 'Field pattern in x-y plane'
	axis equal, axis tight
end

% clf, im(angle(smap(:,:,2))), cbar

% ir_mri_smap_r(r, z)
% function for testing near 0
function out = ir_mri_smap_r(r, z)
M = 4 * r ./ ((1 + r).^2 + z.^2); % = k^2, see ellipke
[K E] = ellipke(M);
out = 2 * z ./ r .* ((1 + r).^2 + z.^2).^(-0.5) .* ...
	((1 + r.^2 + z.^2) ./ ((1 - r).^2 + z.^2) .* E - K);

% ir_mri_smap1()
% based on grivich:00:tmf
% for a circular coil in "x-y plane" of radius "a"
% note that coil x-y plane is not same as object x-y plane!
% returns (i,j,k) components of B vector for each (x,y,z) location
function [smap_x smap_y smap_z] = ir_mri_smap1(x, y, z, a)
x = x ./ a; % normalized units
y = y ./ a;
z = z ./ a;
r = sqrt(x.^2 + y.^2);
M = 4 * r ./ ((1 + r).^2 + z.^2); % = k^2, see ellipke
[K E] = ellipke(M);
if ir_is_octave
	K = reshape(K, size(M)); % ellipke shape bug in octave
	E = reshape(E, size(M));
end

% the following is B_z in eqn (18) in grivich:00:tmf
% and same as eqn [10] in wang:00:dop to within constant scale factor
smap_z = 2 * ((1 + r).^2 + z.^2).^(-0.5) .* ...
	(K + (1 - r.^2 - z.^2) ./ ((1 - r).^2 + z.^2) .* E);
smap_z = smap_z / a;

if 0 && any(r(:) == 0) % test code to explore when r is near 0
	r0 = linspace(0,5e-7,101);
	z0 = 0.4;
	t0 = ir_mri_smap_r(r0, z0);
	slope = 3*pi * z0 / ((1+z0^2)^2.5);
	clf, plot(r0, t0, '-', r0, slope * r0, '--'); grid, prompt
end

% the following is B_r in eqn (17) in grivich:00:tmf
smap_r = 2 * z ./ r .* ((1+r).^2 + z.^2).^(-0.5) .* ...
	((1 + r.^2 + z.^2) ./ ((1-r).^2 + z.^2) .* E - K);
bad = abs(r) < 1e-6;
smap_r(bad) = 3 * pi * z(bad) ./ ((1 + z(bad).^2).^2.5) .* r(bad);
smap_r = smap_r / a;

if any(isnan(smap_r(:))) || any(isnan(smap_z(:)))
	keyboard
end

smap_x = smap_r .* div0(x, r);
smap_y = smap_r .* div0(y, r);
%smap_z = smap_r .* div0(z, r);

%phi = atan2(y, x);
%smap_x = smap_r .* cos(phi);
%smap_y = smap_r .* sin(phi);

% ir_mri_sensemap_sim_test0
% see ellipke
function ir_mri_sensemap_sim_test0
m = linspace(0,1,101);
[k e] = ellipke(m);
clf, plot(m, k, '-', m, e, '--'), legend('k', 'e')
yaxis_pi('0 p/2 p 3*p/2')

% ir_mri_sensemap_sim_test1
% test ir_mri_smap1 routine, cf Fig. 4 of grivich:00:tmf
function ir_mri_sensemap_sim_test1
a = 1;
x = linspace(-2,2,99);
y = linspace(-2,2,97);
zlist = [0.001 0.1 0.2 0.5 1.0];
zlist(1) = [];
[xx yy zz] = ndgrid(x, y, zlist);
[smap_x smap_y smap_z] = ir_mri_smap1(xx, yy, zz, a);
smap_b = sqrt(smap_x.^2 + smap_y.^2);
if im
	im('plc', 4, numel(zlist))
	ir_mri_sensemap_sim_test1_show(smap_x, x, y, 0, zlist, 'x')
	ir_mri_sensemap_sim_test1_show(smap_y, x, y, 1, zlist, 'y')
	ir_mri_sensemap_sim_test1_show(smap_z, x, y, 2, zlist, 'z')
	ir_mri_sensemap_sim_test1_show(smap_b, x, y, 3, zlist, 'b')
end

% ir_mri_sensemap_sim_test1_show()
function ir_mri_sensemap_sim_test1_show(map, x, y, offset, zlist, titl)
clim = [-20 20];
for iz = 1:numel(zlist)
	p = offset*numel(zlist) + iz;
	im(p, x, y, map(:,:,iz), titl, clim), cbar
	axis equal
	if iz == 1
		xtick([-2 2]), ytick([-2 2])
	else
		xtick off, ytick off
	end
	if zlist(iz) < 0.5
		blim = [7 12 19];
	else
		blim = [1 3 5];
	end
	hold on
	contour(x, y, abs(map(:,:,iz))', blim, 'b-')
	contour(x, y, abs(map(:,:,iz))', [0 0]+0.001, 'g-')
	hold off
	if streq(titl, 'b')
		xlabelf('z = %g', zlist(iz))
	end
end
drawnow

% ir_mri_sensemap_sim_test2
% 2D test case
function ir_mri_sensemap_sim_test2

[smap x y] = ir_mri_sensemap_sim('chat', 1, 'nx', 32, ...
	'rcoil', [], 'ncoil', 4, 'coil_distance', 1.2);

if 0 % check vs old version
	old = mri_sensemap_sim('chat', 1, 'nx', 32, ...
		'rcoil', [], 'ncoil', 4, 'coil_distance', 1.2);
	equivs(smap, old)
end

if 1 % check rotational symmetry in 4-coil case
	for ic=2:4
		tmp = rot90(smap(:,:,1), (ic-1));
		equivs(abs(tmp), abs(smap(:,:,ic)))
	%	im(8+ic, x, y, abs(tmp)), cbar
		p1 = angle(tmp) + (ic-1) * pi/2; % add pi/2 to rotated
		p2 = angle(smap(:,:,ic));
		tmp = exp(1i * (p2 - p1));
		equivs(tmp, ones(size(tmp))) % trick: equivs mod 2*pi
	end
end

% ir_mri_sensemap_sim_test3
% illustrate 3D sense maps with el
function ir_mri_sensemap_sim_test3
nring = 3;
ncoil = 4 * nring;
ig = image_geom('nx', 16, 'ny', 14, 'nz', 10, 'fov', 200, 'dz', 20); % 20cm fov
%ig = image_geom('nx', 72, 'ny', 48, 'nz', 12, 'fov', 22, 'zfov', 10); % michelle
%nring = 2; ncoil = 8; % michelle
ig.mask = ig.circ > 0;
smap = ir_mri_sensemap_sim('chat', 1, 'nx', ig.nx, 'ny', ig.ny, 'nz', ig.nz, ...
	'dx', ig.dx, ...
	'dz', ig.dz, ...
	'orbit_start', 1*[0 45 0], ...
...%	'orbit_start', 0*[0 0], ...
	'rcoil', 70, ...
...%	'rcoil', 3, ...
	'nring', nring, 'ncoil', ncoil, 'coil_distance', 1.2);
if im
	prompt
	tmp = smap .* repmat(ig.mask, [1 1 1 ncoil]);
	im clf, im('row', ncoil, reshape(abs(tmp), [ig.nx ig.ny ig.nz*ncoil]))
	prompt
	tmp = permute(tmp, [1 3 2 4]); % [nx nz ny ncoil] z cuts are smooth
	im('row', ncoil, reshape(abs(tmp), [ig.nx ig.nz ig.ny*ncoil]))
	prompt
	im('row', ncoil, abs(tmp(:,:,end/2,:)))
end

% ir_mri_sensemap_sim_test()
function ir_mri_sensemap_sim_test(arg)
switch(arg)
case 'test0'
	ir_mri_sensemap_sim_test0 % ellipk
case 'test1'
	ir_mri_sensemap_sim_test1 % basic test
case 'test2'
	ir_mri_sensemap_sim_test2 % 2D test
case 'test3'
	ir_mri_sensemap_sim_test3 % 3D test
case 'test'
	ir_mri_sensemap_sim_test1
	if im, prompt, end
	ir_mri_sensemap_sim_test2
	if im, prompt, end
	ir_mri_sensemap_sim_test3
otherwise
	fail('bad argument "%s"', arg)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function A = whtmtx(N)
%WHTMTX Generates sequency-ordered Walsh-Hadamard transformation matrix.
%   A = WHTMTX(N) returns the sequency-ordered transformation matrix for
%   a Walsh-Hadamard transform of dimension N, where N is a power of 2.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

bits = log2(N);
if rem(bits,1) > 0
   error('N must be a power of 2!');
end

% Scale MATLAB's naturally ordered Hadamard functions
H = hadamard(N) ./ sqrt(N);

% Create a vector of row numbers in binary. Subtract character '0' to make
% the vector of type logical rather than character.
rows = dec2bin(0:N-1,bits) - '0';

% Gray code the row numbers. Note that we are working from the most to
% least significant bit. That is, column 1 of gray contains the MSB.
gray = rows;
for i = 2:1:bits
   gray(:,i) = xor(rows(:,i),rows(:,i-1));
end

% Bit-reverse the gray coded row numbers and convert to new decimal row
% numbers. Then rearrange the rows.
rows = fliplr(gray) * pow2((bits-1:-1:0)');
A = H(rows+1,:);

author: dipum

repo: dipum-toolbox

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


text:
%% Simulation Parameters
name        ='SBBGK1d'; % Simulation Name
CFL         = 0.05;     % CFL condition
r_time      = 1/1000;  % Relaxation time
tEnd        = 0.1;      % End time
theta       = 0;        % {-1} BE, {0} MB, {1} FD.
quad        = 2;        % for NC = 1 , GH = 2
method      = 1;        % for TVD = 1, WENO3 = 2, WENO5 = 3
IC_case     = 1;        % IC: {1}Sod's, {2}LE, {3}RE, {4}DS, {5}SS, {6}Cavitation
plot_figs   = 1;        % 0: no, 1: yes please!
write_ans   = 1;        % 0: no, 1: yes please!
RRR=8.314;
% Using DG
P_deg       = 0;        % Polinomial Degree
Pp          = P_deg+1;  % Polinomials Points
% Using RK integration time step
RK_stages   = 4;        % Number of RK stages

gamma = 2.5; %1.4;        % Ratio of specific heats
 etpfix   = 0.90;	% {#} Harten's sonic entropy fix value, {0} no entropy fix
[MB_r,u,MB_p] = Euler_IC1d(x,1);
%% Microscopic Velocity Discretization (using Discrete Ordinate Method)
% that is to make coincide discrete values of microscopic velocities with
% values as the value points for using a quadrature method, so that we can
% integrate the velocity probability distribution to recover our
% macroscopics properties.
switch quad

    case{1} % Newton Cotes Quadrature:
    V  = [-20,20];  % range: a to b
    nv = 200;       % nodes desired (may not the actual value)
    [v,w,k] = cotes_xw(V(1),V(2),nv,5); % Using Netwon Cotes Degree 5

    case{2} % Gauss Hermite Quadrature:
    nv = 60;          % nodes desired (the actual value)
    [v,w] = GaussHermite(nv); % for integrating range: -inf to inf
    k = 1;            % quadrature constant.
    w = w.*exp(v.^2); % weighting function of the Gauss-Hermite quadrature

    otherwise
        error('Order must be between 1 and 2');
end

%% Velocity-Space Grid:
% The actual nv value will be computed using 'lenght' vector function:
nv = length(v);
% Using D.O.M.
    v = repmat(v,1,nx);     w = repmat(w,1,nx);
	p = zeros(nv,nx);
    %% Load Macroscopic Velocity, Temperature and Fugacity
    [rho_0,u0,t0] = BGK_IC1d(x,IC_case);
 rho_0 = repmat(rho_0,nv,1);                       ux = repmat(u0,nv,1);       t = repmat(t0,nv,1);
    f0 = f_equilibrium_1d(rho_0,ux,v,t,RRR);

    %% Marching Scheme
% First we need to define how big is our time step. Due to the discrete
% ordinate method the problem is similar to evolve the same problem for
% every mesoscopic velocity.
dt = dx*CFL/max(v(:,1));
dtdx = dt/dx;  % precomputed to save someflops

% Time domain discretization
time = 0:dt:tEnd;

repo: Aero-matlab