text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [result] = weiss(evaluation_and_training_days, setup, fid)

    % load parameters of algorithm
    dataset = setup.dataset;
    household = setup.household;
    granularity = setup.granularity;
    filteringMethod = setup.filtering;
    filtLength = setup.filtLength;
    plevelMinLength = setup.plevelMinLength;
    maxEventDuration = setup.maxEventDuration;
    eventPowerStepThreshold = setup.eventThreshold;
    r = setup.r;
    osc = setup.osc;
    evaluation_days = evaluation_and_training_days{1};
    training_days = evaluation_and_training_days{2};
    appliance = setup.appliance;

    % set variables
    edgeThreshold = 2;
    num_measurements = size(evaluation_days,1)*86400;

    % build signature database
    global caching;
    if caching == 1
        if exist('cache_sigdat.mat') == 2
            load('cache_sigdat');
        else
            signature_database = buildSignatureDatabase(setup, training_days);
            save('cache_sigdat', 'signature_database');
        end
    else
        signature_database = buildSignatureDatabase(setup, training_days);
    end

    signatures = signature_database.signatures;
    names_of_signatures = signature_database.names;
    phases_of_signatures = signature_database.phases;
    % calculate length (2-norm) of each signature
    numOfSignatures = size(signatures,1);
    signatureLength = zeros(numOfSignatures,1);
    for j = 1:numOfSignatures
        signatureLength(j,1) = norm(signatures(j,:));
    end
    signature_database.signatureLength = signatureLength;

    % write signatures to text file
    fprintf(fid,'%20s %13s %16s %9s\n', 'appliance:', 'true power:', 'reactive power:', 'phase:');
    for i = 1:size(signatures,1)
        fprintf(fid,'%20s %13.2f %16.2f %9s\n', cell2mat(names_of_signatures(i)), signatures(i,1),...
             signatures(i,2), num2str(phases_of_signatures(i,1)));
    end
    fprintf(fid, '\n\n');

    result = struct;
    result.events = [];
    result.appliance_names = {};

    input_params = struct;
    input_params.dataset = dataset;
    input_params.household = household;
    input_params.evaluation_days = evaluation_days;
    input_params.granularity = granularity;
    input_params.filtering_method = filteringMethod;
    input_params.filtLength = filtLength;
    input_params.edgeThreshold = filtLength;
    input_params.plevelMinLength = plevelMinLength;
    input_params.eventPowerStepThreshold = eventPowerStepThreshold;
    input_params.maxEventDuration = maxEventDuration;

    % Disaggregate all appliances or individual appliance?
    if strcmp(appliance, 'Stove') == 1
        [events, times] = get_events_of_multi_phase_appliance(input_params);
        matching_ids = find(strcmp(signature_database.names, 'Stove'));
        stove_sig = signatures(matching_ids,:);
        result = infer_events_stove(result, events, times, min(stove_sig(stove_sig > 0))-200);
        [result.usage, result.usage_times_start] = infer_usage(appliance, result.events(:,3), result.events(:,1), num_measurements, setup.usage_duration);
        result.consumption = infer_consumption_stove(result.events(:,3), result.events(:,1), num_measurements, setup.usage_duration);

    elseif strcmp(appliance, 'Dishwasher') == 1
        applianceID = getApplianceID(appliance);
        phase = getPhase(household, applianceID, dataset);
        [events, times] = get_events_of_single_phase_appliance(phase, input_params);
        result = infer_events(result, events, times, signature_database, phase, setup, appliance);
        [result.usage, result.usage_times_start] = infer_usage(appliance, result.events(:,3), result.events(:,1), num_measurements, setup.usage_duration);
        result.consumption = infer_consumption(result, evaluation_and_training_days, setup.usage_duration, appliance, setup);

    elseif strcmp(appliance, 'Water kettle') == 1 || strcmp(appliance, 'TV') == 1 || strcmp(appliance, 'Stereo') == 1
            applianceID = getApplianceID(appliance);
            phase = getPhase(household, applianceID, dataset);
            [events, times] = get_events_of_single_phase_appliance(phase, input_params);
            result = infer_events(result, events, times, signature_database, phase, setup, appliance);
            [result.usage, result.usage_times_start] = infer_usage(appliance, result.events(:,3), result.events(:,1), num_measurements, setup.usage_duration);
            result.consumption = infer_consumption(result, evaluation_and_training_days, setup.usage_duration, appliance, setup);
    elseif strcmp(appliance, 'Fridge') == 1 || strcmp(appliance, 'Freezer') == 1
            applianceID = getApplianceID(appliance);
            phase = getPhase(household, applianceID, dataset);
            [events, times] = get_events_of_single_phase_appliance(phase, input_params);
            result = infer_events(result, events, times, signature_database, phase, setup, appliance);
            result.consumption = infer_consumption(result, evaluation_and_training_days, setup.usage_duration, appliance, setup);
    elseif strcmp(appliance, 'Laptop') == 1
            applianceID = getApplianceID(appliance);
            phase = getPhase(household, applianceID, dataset);
            [events, times] = get_events_of_single_phase_appliance(phase, input_params);
            result = infer_events(result, events, times, signature_database, phase, setup, appliance);
%             result.consumption = infer_consumption(result, evaluation_and_traini
    elseif strcmp(appliance, 'All')

        %% "Old": all appliances
        signatures_in_previous_phases = 0;
        for phase = 1:3
            if caching == 1
                error('Does not work with caching - run without caching');
            end
           [event_vecs, timeOfEvents] = get_events_of_single_phase_appliance(phase, input_params);
            % assign each event to its best match in the signature database
            [signatureIDs, dist] = knnsearch(signatures(phases_of_signatures == phase, 1:2), event_vecs(:,1:2));
            if ~isempty(signatureIDs)
                dist_threshold = r*signatureLength(signatureIDs,1) + event_vecs(:,4);
                matching_valid = dist < dist_threshold;
                %% should te next line not be earlier???
                signatureIDs = signatureIDs + signatures_in_previous_phases;
                result.events = [result.events; timeOfEvents(matching_valid), signatureIDs(matching_valid), event_vecs(matching_valid, 1:3)];
                signatures_in_previous_phases = signatures_in_previous_phases + nnz(phases_of_signatures == phase);
            end
        end

    else
        error('Appliance not supported');
    end

    % store labeled events
    for i = 1:length(names_of_signatures)
        if ismember(cell2mat(names_of_signatures(i)), result.appliance_names)
            result.events(result.events(:,2) == i,2) = find(ismember(result.appliance_names, names_of_signatures{i}));
        else
            % do nothing if signature is not in signature database
            % result.appliance_names{end+1} = cell2mat(names_of_signatures(i));
            % result.events(result.events(:,2) == i,2) = length(result.appliance_names);
        end
    end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function g = s2g(s, epsilon);

% G = s2g(S, EPSILON)
%
% Scattering to Hybrid-G transformation
% G and S are matrices of size [2,2,F]
% where F is the number of frequencies
% (the number of ports is always 2)
%
% EPSILON is a limit used in finding correspondent Hybrid-G matrices in the
% vicinity of singularities; by default 1e-12, should be enough for most
% realistic problems; could be increased for a gain in speed
%
% written by tudor dima, tudima@zahoo.com, change the z into y

if nargin < 2, epsilon = 1e-12; end;

d = (1+s(1,1,:)).*(1-s(2,2,:)) + s(1,2,:).*s(2,1,:);
[n,i] = min(abs(d));
exact_g = 1;
while n <= epsilon
    exact_g = 0;
    p1 = 1+round(rand); p2 = 1+round(rand);
    s(p1,p2,i) = s(p1,p2,i)+(rand-0.5)*epsilon;
    d = (1+s(1,1,:)).*(1-s(2,2,:)) + s(1,2,:).*s(2,1,:);
    [n,i] = min(abs(d));
end;

if exact_g == 0
    fprintf(1,'%s\n%s\n', 's2g: correspondent G-hybrid matrix non-existent', ...
        'an approximation is produced');
end;

g(1,1,:) = ((1 - s(1,1,:)) .* (1 - s(2,2,:)) - s(1,2,:).*s(2,1,:))./d;
g(1,2,:) = -2*s(1,2,:)./d;
g(2,1,:) = 2*s(2,1,:)./d;
g(2,2,:) = ((1 + s(1,1,:)) .* (1 + s(2,2,:)) - s(1,2,:).*s(2,1,:))./d;

author: Sable

repo: mcbench-benchmarks


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Problem 5
%                                    0.9^n, 0<=n<=3
%  z-Transform computation of f[n]=  2^-n , 4<=n<=6
%                                     1,    7<=n<=10

n1=0:3;
f1=0.9.^n1;
n2=4:6;
f2=2.^(-n2);
n3=7:10;
f3=ones(size(n3));
n=[n1 n2 n3];
f=[f1 f2 f3];
stem(n,f);
axis([-.5 10.5 -.1 1.1]);

figure
syms n z
n1=4;
n2=7;
n3=11;
f1=0.9^n;
f2=2^(-n);
f3=1;
f=f1*(heaviside(n)-heaviside(n-n1)) +f2*(heaviside(n-n1)-heaviside(n-n2)) +f3*(heaviside(n-n2)-heaviside(n-n3));
n_s=0:10;
f_s=subs(f,n_s);
stem(n_s,f_s);
axis([-.5 10.5 -.1 1.1]);

% a)z-Transform of f[n]
F1=ztrans(f,z)

% b)z-Transform of f[n]
n1=0:3;
f1=0.9.^n1;
n2=4:6;
f2=2.^(-n2);
n3=7:10;
f3=ones(size(n3));
n=[n1 n2 n3];
f=[f1 f2 f3];
F2=sum(f.*(z.^-n))

% c) confirmation
figure
syms n
ftest=iztrans(F1,n)
n=0:10;
ftest1=subs(ftest,n);
stem(n,ftest1);
axis([-.5 10.5 -.1 1.1]);


text:
function lambda = minij_eigenvalues ( n )

%*****************************************************************************80
%
%% MINIJ_EIGENVALUES returns the eigenvalues of the MINIJ matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  for i = 1 : n
    angle = ( 2 * i - 1 ) * pi / ( 2 * n + 1 );
    lambda(i,1) = 0.5 / ( 1.0 - cos ( angle ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function K = covMaternard(d, hyp, x, z, i)

% Matern covariance function with nu = d/2 and Automatic Relevance Detemination
% (ARD) distance measure. For % d=1 the function is also known as the exponential
% covariance function or the Ornstein-Uhlenbeck covariance. The covariance
% function is:
%
%   k(x^p,x^q) = s2f * f( sqrt(d)*r ) * exp(-sqrt(d)*r)
%
% with f(t)=1 for d=1, f(t)=1+t for d=3 and f(t)=1+t.*(1+t/3) for d=5.
% Here r is the distance sqrt((x^p-x^q)'*inv(P)*(x^p-x^q)), P is ell times
% the unit matrix and sf2 is the signal variance. The hyperparameters are:
%
% hyp = [ log(ell_1)
%         log(ell_2)
%          .
%         log(ell_D)
%         log(sqrt(sf2)) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-09-10.
%                  Ali Abusnina, 2014
%
% See also COVFUNCTIONS.M.

if nargin<3, K = '(D+1)'; return; end              % report number of parameters
if nargin<4, z = []; end                                   % make sure, z exists
xeqz = numel(z)==0; dg = strcmp(z,'diag') && numel(z)>0;        % determine mode

[n,D] = size(x);
ell = exp(hyp(1:D));                               % characteristic length scale
sf2 = exp(2*hyp(D+1));

switch d
  case 1, f = @(t) 1;               df = @(t) 1;
  case 3, f = @(t) 1 + t;           df = @(t) t;
  case 5, f = @(t) 1 + t.*(1+t/3);  df = @(t) t.*(1+t)/3;
end
          m = @(t,f) f(t).*exp(-t); dm = @(t,f) df(t).*t.*exp(-t);

% precompute distances
if dg                                                               % vector kxx
  K = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
      K = sqrt( sq_dist(sqrt(d)* (bsxfun(@rdivide, x', ell))    ));

  else                                                   % cross covariances Kxz

    K = sqrt( sq_dist(sqrt(d)* (bsxfun(@rdivide, x', ell)) , sqrt(d) * (bsxfun(@rdivide, z', ell))) );

  end
end

if nargin<5                                                        % covariances
  K = sf2*m(K,f);
else                                                               % derivatives
  if i<=D
    K = sf2*dm(K,f);
  elseif i<=D+1
    K = 2*sf2*m(K,f);
  else
    error('Unknown hyperparameter')
  end
end

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function [final_sur ] = detect_sym( obblist, Rmat, Dmat, RELATIONS, LABELLEN )

% detect the surround and symmetry relations in the obblist
% The symmetry relation is disabled in this version

% checked_label = [];
% load(pathparam.label2vec_filename);
% box_vec = label2vec_map('Box');
% box_label = binary2dec(box_vec{1});

% search for the candidate_list
% each candidate is a list of objects with similar sizes which may be
% involved in a symmetry or surround relation
objnum = size(obblist,1);
sizelist = obblist(:,10:12);
size_thre = 0.03;
matched = [];
candidate_list = {};
for i = 1:objnum
    idx = find(matched==i);
    if(length(idx)==0)
        delta = abs(sizelist - repmat(sizelist(i,:),objnum,1));
        delta = delta - size_thre;
        delta = delta(:,1)+delta(:,2)+delta(:,3);
        idx = find(delta<=0);
        matched = union(matched,idx);
        if(length(idx)>1)
            candidate_list{length(candidate_list)+1} = idx;
        end
    end
end

%% extract the objects with similar heights in each candidate
ncandidate_list = {};
up = [0;1;0];
height_thre = 0.1;
for i = 1:length(candidate_list)
    clist = candidate_list{i};
    heights = obblist(clist,1:3)*up;
    matchednums = zeros(length(clist),1);
    for j = 1:length(clist)
        tmplist = abs(heights-heights(j));
        ind = find(tmplist<height_thre);
        matchednums(j) = length(ind);
    end
    [~,I] = max(matchednums);
    tmplist = abs(heights-heights(j));
	ind = find(tmplist<height_thre);
  	if(length(ind)>1)
        ncandidate_list{length(ncandidate_list)+1} = clist(ind);
    end
end
candidate_list = ncandidate_list;

%% search for all sets without central obj
% [ sym ] = search_symmetryset( candidate_list, obblist, Rmat, Dmat, LABELLEN );
%% search for all sets surrounding one obj
[ cen_sym, cen_sur ] = search_surroundset( candidate_list, obblist, Rmat, Dmat, RELATIONS, LABELLEN );
sur = {};
for i = 1:size(cen_sym,1)
    sur{size(sur,1)+1,1} = cen_sym{i,1};
    sur{size(sur,1),2} = cen_sym{i,2};
end
for i = 1:size(cen_sur,1)
    sur{size(sur,1)+1,1} = cen_sur{i,1};
    sur{size(sur,1),2} = cen_sur{i,2};
end
% sur = [cen_sym,cen_sur];

% sort the sur pairs based on the distance
dlist = zeros(size(sur,1),1);
for i = 1:size(sur,1)
    list = sur{i,1};
    c = sur{i,2};
    dlist(i) = Dmat(list(1),c);
end
[~,I] = sort(dlist);
sur = sur(I,:);

% except for the central object, each object can only be involved in one surround relation
isinvolved = zeros(size(obblist,1),1);
final_sur = {};
for i = 1:size(sur,1)
    list = sur{i,1};
    if(sum(isinvolved(list))==0)
        % add this surround pair
        s = [sur{i,2};list(:)];
        final_sur{length(final_sur)+1} = s;
    end
end

end

repo: GRAINS


text:
function matlabVal = java2matlab(javaVal)

  if (strcmp(class(javaVal), 'java.util.Vector') == 1)
    if (javaVal.size() ~= 2)
      error('Incorrect Vector size (~= 2)');
    end
    fieldType=javaVal.get(0);
    value=javaVal.get(1);

    switch fieldType
     case 1 % function_handle
      matlabVal=str2func(value(2:end));

     case 2 % struct
      if (strcmp(class(value), 'java.util.Vector') == 0)
        error('Struct array value is not a java.util.Vector');
      end
      if (value.size() < 1)
        error('Struct array is empty.');
      end
      dim = value.get(0);
      matlabVal=struct('x', cell(dim'));
      matlabVal=rmfield(matlabVal, 'x');

      for k=1:numel(matlabVal)
        struc = value.get(k);
        keys= struc.keySet();
        iterator=keys.iterator();

        while iterator.hasNext()
          key=iterator.next();
          val=struc.get(key);
          matlabVal(k).(key)=java2matlab(val);
        end
      end

     case 3 % cell
      if (strcmp(class(value), 'java.util.Vector') == 0)
        error('Cell array value is not a java.util.Vector');
      end
      if (value.size() < 1)
        error('Java cell-array vector is empty.');
      end
      dim=value.get(0);
      matlabVal=cell(dim');
      if (value.size() < numel(matlabVal)+1)
        error('Java cell-array vector is too small.');
      end
      for i=1:numel(matlabVal)
        matlabVal{i} = java2matlab(value.get(i));
      end

     case 4 % matlab matrix
      if (strcmp(class(value), 'java.util.Vector') == 0)
        error('Matrix value is not a java.util.Vector');
      end
      if (value.size() < 1)
        error('Java matrix vector is empty.');
      end
      dim=value.get(0);
      matlabVal=reshape(value.get(1), reshape(dim,1,length(dim)));

     case 5 % java
      matlabVal = value;

     case 6 % slaveRef
      matlabVal = slaveRef([], value);

    end
  else
    matlabVal = javaVal;
  end

author: huashiyiqike

repo: LSTM-MATLAB


text:
clear all; close all; clc

n=100; L=4;
x=linspace(0,L,n);
f=(x.^2).';  % parabola with 100 data points

M=21;  % polynomial degree
for j=1:M
    phi(:,j)=(x.').^(j-1);  % build matrix A
end

trials=[2 10 100];
for j=1:3
    for jj=1:trials(j)
       f=(x.^2+0.2*randn(1,n)).';
       a1=pinv(phi)*f; f1=phi*a1; E1(jj)=norm(f-f1)/norm(f);
       a2=phi\f; f2=phi*a2; E2(jj)=norm(f-f2)/norm(f);
       [a3,stats]=lasso(phi,f,'Lambda',0.1); f3=phi*a3; E3(jj)=norm(f-f3)/norm(f);
       A1(:,jj)=a1; A2(:,jj)=a2; A3(:,jj)=a3;
    end
    A1m=mean(A1.'); A2m=mean(A2.'); A3m=mean(A3.');
    Err=[E1; E2; E3];

    subplot(3,3,j), bar(A1m), axis([0 21 -1 1.2])
    subplot(3,3,3+j), bar(A2m), axis([0 21 -1 1.2])
    subplot(3,3,6+j), bar(A3m), axis([0 21 -1 1.2])
end

%subplot(3,3,j), bar(A1m,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis([0 21 -1 1.2])
%subplot(3,3,3+j), bar(A2m,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis([0 21 -1 1.2])
%subplot(3,3,6+j), bar(A3m,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis([0 21 -1 1.2])

%%
figure(2)

Atot=[A1m; A2m; A3m];  % average loadings of three methods
Atot2=(Atot>0.2).*Atot; % threshold
Atot3=[Atot; Atot2];    % combine both thresholded and not

figure(3), bar3(Atot.'), axis([0 4 0 20 0 1]), view(-70,38), set(gca,'Xtick',[],'Fontsize',[18])
figure(4), bar3(Atot2.'), axis([0 4 0 20 0 1]), view(-70,38), set(gca,'Xtick',[],'Fontsize',[18])

n=200; L=8;
x=linspace(0,L,n);
x1=x(1:100);   % train (interpolation)
x2=x(101:200); % test (extrapolation)

ftrain=(x1.^2).';  % interpolated parabola x=[0,4]
ftest=(x2.^2).';   % extrapolated parbola x=[4,5]

for j=1:M
    phi_i(:,j)=(x1.').^(j-1); % interpolation key
    phi_e(:,j)=(x2.').^(j-1); % extrapolation key
end

for jj=1:6 % compute inter/extra-polation scores
    ani=Atot3(jj,:).';
    fnai=phi_i*ani;
    Eni(jj)=norm(ftrain-fnai)/norm(ftrain);
    fnae=phi_e*ani;
    Ene(jj)=norm(ftest-fnae)/norm(ftest);
end

figure(5)
subplot(3,2,1), bar(Eni)
subplot(3,2,2), bar(Ene)

subplot(3,2,3), bar(Eni), axis([0 7 0 0.01])
subplot(3,2,4), bar(Ene), axis([0 7 0 0.1])

figure(3)
legend('','','','Location','NorthEast')
legend boxoff

%% Dendrograms

author: dynamicslab


text:
function [mu ul ll] = circ_mean(alpha, w, dim)
%
% mu = circ_mean(alpha, w)
%   Computes the mean direction for circular data.
%
%   Input:
%     alpha	sample of angles in radians
%     [w		weightings in case of binned angle data]
%     [dim  compute along this dimension, default is 1]
%
%     If dim argument is specified, all other optional arguments can be
%     left empty: circ_mean(alpha, [], dim)
%
%   Output:
%     mu		mean direction
%     ul    upper 95% confidence limit
%     ll    lower 95% confidence limit
%
% PHB 7/6/2008
%
% References:
%   Statistical analysis of circular data, N. I. Fisher
%   Topics in circular statistics, S. R. Jammalamadaka et al.
%   Biostatistical Analysis, J. H. Zar
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html

if nargin < 3
  dim = 1;
end

if nargin < 2 || isempty(w)
  % if no specific weighting has been specified
  % assume no binning has taken place
	w = ones(size(alpha));
else
  if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
    error('Input dimensions do not match');
  end
end

% compute weighted sum of cos and sin of angles
r = sum(w.*exp(1i*alpha),dim);

% obtain mean by
mu = angle(r);

% confidence limits if desired
if nargout > 1
  t = circ_confmean(alpha,0.05,w,[],dim);
  ul = mu + t;
  ll = mu - t;
end

author: brainstorm-tools


text:
function qr = quatrot(w);
% QUATROT - returns the quaternion representing the rotation given by
%           the 3D vector w (rotation of ||w|| around the direction given by w)
%
%   qr = quatrot(w);
%
%  ON - 3/99
%

th = norm(w);

qr = [cos(th/2) sin(th/2)*w(:)'/th];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ J, grad, u ] = optimal_control ( x, e_conn, q, y1, y2, Md, Nd, L )

%*****************************************************************************80
%
%% OPTIMAL_CONTROL is a script to solve the optimal control problem.
%
%  Discussion:
%
%    The differential equation being solved has the form:
%
%      ( q(x) u_x(x) )_x = f(x)
%
%    The solution is approxiamted using 1D linear finite elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Parameters:
%
%    Input, X
%
%    Input, E_CONN
%
%    Input, Q
%
%    Input, Y1
%
%    Input, Y2
%
%    Input, MD
%
%    Input, ND
%
%    Input, L
%
%    Output, J
%
%    Output, GRAD
%
%    Output, U
%
  alpha = 0.000003;

  [n_nodes   , n_dimensions] = size(x     );
  [n_elements, nel_dof     ] = size(e_conn);

  n_equations      = n_nodes-2;
  ide(1)           = -1;
  ide(2:n_nodes-1) = 1:n_equations;
  ide(n_nodes)     = -2;

  dir              = zeros(2,1);

  n_gauss          = 3;  % number of points used in Gaussian integration

% q = x./(1-3*x.^2);
% q = x.^3;
% q = x.*(1-x);
%---------------------------------------------------------------------
%  Build the finite element matrices
%---------------------------------------------------------------------

  M = sparse(n_nodes    ,n_nodes    );
  K = sparse(n_equations,n_equations);
  b = zeros (n_equations,1          );
  u = zeros (n_nodes    ,1          );

  [r,w] = oned_gauss(n_gauss);

  o_g   = ones(size(r));

  for n_el=1:n_elements

    % compute value of each test function and their spatial derivaties
    % at the integration points
    nodes_local       = e_conn(n_el,:);
    x_local           = x(nodes_local,:);
    [x_g,w_g,phi,p_x] = oned_shape(x_local,r,w);

    q_local           = q(nodes_local);
    q_g               = phi*q_local;

    % compute the value of functions at the Gauss points
    f_g   = f_function(x_g, y1, y2, Md, Nd, L);

    %--------------------------------------------------------------------
    %  Integrate the weak form of the equations (element contributions)
    %--------------------------------------------------------------------
    M_loc = oned_bilinear(o_g, phi, phi, w_g);
    K_loc =-oned_bilinear(q_g, p_x, p_x, w_g);
    b_loc = oned_f_int   (f_g, phi,      w_g);

    %-----------------------------------------------------------------
    % Assemble contributions into the global system matrix
    %-----------------------------------------------------------------

    M(nodes_local,nodes_local) = M(nodes_local,nodes_local) + M_loc;

    for n_t=1:nel_dof
      n_test = ide(nodes_local(n_t));
      if (n_test > 0)  % this is an unknown, fill the row
        for n_u=1:nel_dof
          n_unk = ide(nodes_local(n_u));
          if (n_unk > 0)
            K(n_test,n_unk) = K(n_test,n_unk) + K_loc(n_t,n_u);
          end
        end
        b(n_test) = b(n_test) + b_loc(n_t);
      end
    end

  end

  %-----------------------------------------------------------------------
  %  Perform Implicit Solve
  %-----------------------------------------------------------------------
%  figure(5); spy(A);  % view the sparsity pattern in A

  du = K\b;

  %-----------------------------------------------------------------------
  %  Construct the Solution (here, we simply apply the Dirichlet bc's)
  %-----------------------------------------------------------------------
  for n=1:n_nodes
    i = ide(n);
    if (i>0)
      %  get the nodal value out of the linear system solve
      u(n) = du(i);
    else
      %  get the nodal value from the specified Dirichlet bc
      u(n) = dir(-i);
    end
  end

%  figure(12)
%  plot(x,u)
%  hold on
%  up = u_hat(x, y);
%  plot(x,up,'r+')

  lam = zeros(n_nodes  ,1);
  c = zeros(n_equations,1);
  %   Compute the right hand side of the adjoint equation
  for n_el=1:n_elements

    nodes_local       = e_conn(n_el,:);
    x_local           = x(nodes_local,:);
    [x_g,w_g,phi,p_x] = oned_shape(x_local,r,w);

    u_local           = u(nodes_local);
    u_g               = phi*u_local;

    % compute the value of functions at the Gauss points
    uhat_g   = u_hat(x_g, y1, Md, L);

    %--------------------------------------------------------------------
    %  Integrate the weak form of the equations (element contributions)
    %--------------------------------------------------------------------
    c_loc = oned_f_int   ( uhat_g-u_g, phi,      w_g);

    %-----------------------------------------------------------------
    % Assemble contributions into the global system matrix
    %-----------------------------------------------------------------
    for n_t=1:nel_dof
      n_test = ide(nodes_local(n_t));
      if (n_test > 0)  % this is an unknown, fill the row
        c(n_test) = c(n_test) + c_loc(n_t);
      end
    end

  end

  J = 0.5*sqrt(c'*M(2:end-1,2:end-1)*c);
  dlam = K'\c;

  for n=1:n_nodes
    i = ide(n);
    if (i>0)
      lam(n) = dlam(i);
    else
      lam(n) = dir(-i);
    end
  end

%  figure
%  plot(x,lam)

  grad = zeros(n_nodes,1);

  % Compute the gradient
  for n_el=1:n_elements

    nodes_local       = e_conn(n_el,:);
    x_local           = x(nodes_local,:);
    [x_g,w_g,phi,p_x] = oned_shape(x_local,r,w);

    u_local           = u(nodes_local);
    u_g               = phi*u_local;
    ux_g              = p_x*u_local;

    q_local           = q(nodes_local);
    q_g               = phi*q_local;
    qx_g              = p_x*q_local;

    lam_local         = lam(nodes_local);
    lam_g             = phi*lam_local;
    lamx_g            = p_x*lam_local;

%     grad_loc = alpha*oned_f_int( q_g, phi, w_g) ...
%               - oned_f_int( lamx_g.*ux_g, phi, w_g);
    grad_loc = alpha*oned_f_int( qx_g, p_x, w_g) ...
              - oned_f_int( lamx_g.*ux_g, phi, w_g);

    for n_t=1:nel_dof
      n_test = nodes_local(n_t);
      grad(n_test) = grad(n_test) + grad_loc(n_t);
    end

  end

%  figure
%  plot(x,grad)
%  title('gradient')

  return
end


text:
function [ error_per_image, err_pp, err_pp_dim, frontal_ids ] = compute_error_small( ground_truth_all, detected_points_all )
%compute_error
%   compute the average point-to-point Euclidean error normalized by the
%   inter-ocular distance (measured as the Euclidean distance between the
%   outer corners of the eyes)
%
%   Inputs:
%          grounth_truth_all, size: num_of_points x 2 x num_of_images
%          detected_points_all, size: num_of_points x 2 x num_of_images
%   Output:
%          error_per_image, size: num_of_images x 1

num_of_images = size(ground_truth_all,3);
ground_truth_all_orig = ground_truth_all;

ground_truth_all = ground_truth_all([1:60,62:64,66:end],:,:);
ground_truth_all = ground_truth_all(18:end,:,:);

num_of_points = size(ground_truth_all,1);
frontal_ids = true(num_of_images,1);

if(size(detected_points_all,1) == 68)
    detected_points_all = detected_points_all([1:60,62:64,66:end],:,:);
    detected_points_all = detected_points_all(18:end,:,:);
end

error_per_image = zeros(num_of_images,1);
err_pp = zeros(num_of_images, num_of_points);
err_pp_dim = zeros(num_of_images, num_of_points, 2);

for i =1:num_of_images
    detected_points      = detected_points_all(:,:,i);
    ground_truth_points  = ground_truth_all(:,:,i);
    visible = ground_truth_points(:,1) > 0;
    ground_truth_points_orig = ground_truth_all_orig(:,:,i);
    visible_orig = ground_truth_points_orig(:,1) > 0;

    if(sum(visible_orig) < 55)
        frontal_ids(i) = false;
    end

    normalization_x = max(ground_truth_points_orig(visible_orig,1)) - min(ground_truth_points_orig(visible_orig,1));
    normalization_y = max(ground_truth_points_orig(visible_orig,2)) - min(ground_truth_points_orig(visible_orig,2));
    normalization = (normalization_x + normalization_y)/2;

    sum_c=0;
    for j=1:num_of_points
        if(visible(j))
            sum_c = sum_c+norm(detected_points(j,:)-ground_truth_points(j,:));
            err_pp(i,j) = norm(detected_points(j,:)-ground_truth_points(j,:));
            err_pp_dim(i,j,1) = detected_points(j,1)-ground_truth_points(j,1);
            err_pp_dim(i,j,2) = detected_points(j,2)-ground_truth_points(j,2);
        end
    end
    error_per_image(i) = sum_c/(sum(visible)*normalization);
    err_pp(i,:) = err_pp(i,:) ./ normalization;
    err_pp_dim(i,:) = err_pp_dim(i,:) ./ normalization;
end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
% Find out how big the cliques are in an HHMM as a function of depth
% (This is how we get the complexity bound of O(D K^{1.5D}).)

if 0
Qsize = [];
Fsize = [];
Nclqs = [];
end

ds = 1:15;

for d = ds
  allQ = 1;
  [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(d, allQ);

  N = length(intra);
  ns = 2*ones(1,N);

  bnet = mk_dbn(intra, inter, ns);
  for i=1:N
    bnet.CPD{i} = tabular_CPD(bnet, i);
  end

  if 0
    T = 5;
    dag = unroll_dbn_topology(intra, inter, T);
    engine = jtree_unrolled_dbn_inf_engine(bnet, T, 'constrained', 1);
    S = struct(engine);
    S1 = struct(S.sub_engine);
  end

  engine = jtree_dbn_inf_engine(bnet);
  S = struct(engine);
  J = S.jtree_struct;

  ss = 2*d+1;
  Qnodes2 = Qnodes + ss;
  QQnodes = [Qnodes Qnodes2];

  % find out how many Q nodes in each clique, and how many F nodes
  C = length(J.cliques);
  Nclqs(d) = 0;
  for c=1:C
    Qsize(c,d) = length(myintersect(J.cliques{c}, QQnodes));
    Fsize(c,d) = length(myintersect(J.cliques{c}, Fnodes));
    if length(J.cliques{c}) > 1 % exclude observed leaves
      Nclqs(d) = Nclqs(d) + 1;
    end
  end
  %pred_max_Qsize(d) = ceil(d+(d+1)/2);
  pred_max_Qsize(d) = ceil(1.5*d);

  fprintf('d=%d\n', d);
  %fprintf('D=%d, max F = %d. max Q = %d, pred max Q = %d\n', ...
	%  D, max(Fsize), max(Qsize), ceil(D+(D+1)/2));

  %histc(Qsize,1:max(Qsize)) % how many of each size?
end % next d

Q = 2;
pred_mass = ds.*(Q.^ds) + Q.^(ceil(1.5 * ds))
pred_mass2 = Q.^(ceil(1.5 * ds))

for d=ds
  mass(d) = 0;
  for c=1:C
    mass(d) = mass(d) + Q^Qsize(c,d);
  end
end

if 0
%plot(ds, max(Qsize), 'o-',  ds, pred_max_Qsize, '*--');
%plot(ds, max(Qsize), 'o-',  ds, 1.5*ds, '*--');
%plot(ds, mass, 'o-',  ds, pred_mass, '*--');
D = 15;
%plot(ds(1:D), mass(1:D), 'bo-',  ds(1:D), pred_mass(1:D), 'g*--', ds(1:D), pred_mass2(1:D), 'k+-.');
plot(ds(1:D), log(mass(1:D)), 'bo-',  ds(1:D), log(pred_mass(1:D)), 'g*--', ds(1:D), log(pred_mass2(1:D)), 'k+-.');

grid on
xlabel('depth of hierarchy')
title('max num Q nodes in any clique vs. depth')
legend('actual', 'predicted')

%previewfig(gcf, 'width', 3, 'height', 1.5, 'color', 'bw');
%exportfig(gcf, '/home/cs/murphyk/WP/ConferencePapers/HHMM/clqsize2.eps', ...
%          'width', 3, 'height', 1.5, 'color', 'bw');

end

if 0
for d=ds
  effnumclqs(d) = length(find(Qsize(:,d)>0));
end
ds = 1:10;
Qs = 2:10;
maxC = size(Qsize, 1);
cost = [];
cost_bound = [];
for qi=1:length(Qs)
  Q = Qs(qi);
  for d=ds
    cost(d,qi) = 0;
    for c=1:maxC
      if length(Qsize(c,d) > 0) % this clique contains Q nodes
	cost(d,qi) = cost(d,qi) + Q^Qsize(c,d)*2^Fsize(c,d);
      end
    end
    %cost_bound(d,qi) = effnumclqs(d) * 8 * Q^(max(Qsize(:,d)));
    cost_bound(d,qi) = (effnumclqs(d)*8) + Q^(max(Qsize(:,d)));
  end
end

qi=2; plot(ds, cost(:,qi), 'o-',  ds, cost_bound(:,qi), '*--');
end

if 0
% convert numbers in cliques into names
for d=1:D
  Fdecode(Fnodes(d)) = d;
end
for c=8:15
  clqs = J.cliques{c};
  fprintf('clique %d: ', c);
  for k=clqs
    if myismember(k, Qnodes)
      fprintf('Q%d ', k)
    elseif myismember(k, Fnodes)
      fprintf('F%d ', Fdecode(k))
    elseif isequal(k, Onode)
      fprintf('O ')
    elseif myismember(k, Qnodes2)
      fprintf('Q%d* ', k-ss)
    else
      error(['unrecognized node ' k])
    end
  end
  fprintf('\n');
end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [annotation, img] = LMimscale(annotation, img, scaling, method)
%
% Scales one image (as in imresize) and the associated annotation.
% The scale factor is given by 'scaling'.
%
% When scaling<1, the image is downsampled.
%
% [annotation, img] = LMimscale(annotation, img, scaling);
%
% or, if you are using the LabelMe struct D:
%
% [D(i).annotation, img] = LMimscale(D(i).annotation, img, scaling);
%
% Also, the scaling can be given in absolute coordinates:
%
% [annotation, img] = LMimscale(annotation, img, [NROWS NCOLS]);

% Store the original numerical format of the image and turn it into 'double'.

imgtype = whos('img');
img = single(img);

if nargin < 4
    method = 'bilinear';
end
if length(scaling) > 1
  scalingx = scaling(2)/size(img,2);
  scalingy = scaling(1)/size(img,1);
  if isnan(scalingx)
      scalingx = scalingy;
  end
  if isnan(scalingy)
      scalingy = scalingx;
  end
else
  scalingx = scaling;
  scalingy = scaling;
end

if scaling ~= 1
    if nargout > 1
      % Image resampling:
      %img = imresizefast(img, scaling, method);
      img = imresize(img, scaling, method);
    end

    % Change the size of the polygon coordinates
    if isfield(annotation, 'object')
        Nobjects = length(annotation.object); n=0;
        for i = 1:Nobjects
            [x,y] = getLMpolygon(annotation.object(i).polygon);
            x = round(x*scalingx);
            y = round(y*scalingy);
            annotation.object(i).polygon = setLMpolygon(x,y);

%             Npoints = length(annotation.object(i).polygon.pt);
%             for j = 1:Npoints
%                 % Scale each point:
%                 x=str2num(annotation.object(i).polygon.pt(j).x);
%                 y=str2num(annotation.object(i).polygon.pt(j).y);
%
%                 x = round(x*scalingx);
%                 y = round(y*scalingy);
%
%                 annotation.object(i).polygon.pt(j).x = num2str(x);
%                 annotation.object(i).polygon.pt(j).y = num2str(y);
%             end
        end
    end
end

% add/modify image size field
annotation.imagesize.nrows = size(img,1);
annotation.imagesize.ncols = size(img,2);

if nargout > 1
  % return the image in its original numeric format
  img = feval(imgtype.class, img);
end

function img = imresizefast(img, scaling, method, init)

if nargin<4
    init = 0;
end

if max(scaling) > .5
    img = imresize(img, scaling, method);
else
    c = size(img,3);
    for n = 1:c
        img(:,:,n) = conv2(img(:,:,n), [1 2 1; 2 4 2; 1 2 1]/16, 'same');
    end
    img = img(init+1:2:end, init+1:2:end, :);
    %img = convn(img, [1 2 1]/4, 'same');
    %img = img(:,init+1:2:end,:);
    img = imresizefast(img, 2*scaling, method, 1-init);
end

author: CSAILVision

repo: LabelMeToolbox


text:
function nvar = p18_nvar ( option )

%*****************************************************************************80
%
%% P18_NVAR sets the number of variables for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, integer NVAR, the number of variables.
%
  nvar = 12;

  return
end


text:
function [NewTessFile, iSurface, I, J] = tess_downsize( TessFile, newNbVertices, Method )
% TESS_DOWNSIZE: Reduces the number of vertices in a surface file.
%
% USAGE:  [NewTessFile, iSurface, I, J] = tess_downsize(TessFile, newNbVertices=[ask], Method=[ask]);
%
% INPUT:
%    - TessFile      : Full path to surface file to decimate
%    - newNbVertices : Desired number of vertices
%    - Method        : {'reducepatch', 'reducepatch_subdiv', 'iso2mesh', 'iso2mesh_project'}
% OUTPUT:
%    - NewTessFile : Filename of the newly created file
%    - iSurface    : Index of the new surface file
%    - I,J         : Indices of the vertices that were kept (see intersect function)

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

%% ===== PARSE INPUTS =====
if (nargin < 3) || isempty(Method)
    Method = [];
end
if (nargin < 2) || isempty(newNbVertices)
    newNbVertices = [];
end
% File name: string or cell array of strings
MultipleFiles = [];
if iscell(TessFile)
    if (length(TessFile) > 1)
        MultipleFiles = TessFile;
    end
    TessFile = TessFile{1};
end
% Save current modifications
panel_scout('SaveModifications');
% Initialize returned values
NewTessFile = '';
iSurface = [];
I = [];
J = [];

%% ===== ASK FOR MISSING OPTIONS =====
% Get the number of vertices
VarInfo = whos('-file',file_fullpath(TessFile),'Vertices');
oldNbVertices = VarInfo.size(1);
% If new number of vertices was not provided: ask user
if isempty(newNbVertices)
    % Ask user the new number of vertices
    newNbVertices = java_dialog('input', 'New number of vertices:', ...
                                         'Resample surface', [], num2str(oldNbVertices));
    if isempty(newNbVertices)
        return
    end
    % Read user input
    newNbVertices = str2double(newNbVertices);
end
% Check if new number of vertices is valid
if isempty(newNbVertices) || isnan(newNbVertices)
    error('Invalid vertices number');
end
if (newNbVertices >= oldNbVertices)
    NewTessFile = TessFile;
    disp(sprintf('TESS> Surface has %d vertices, cannot downsample to %d vertices.', oldNbVertices, newNbVertices));
    return;
end

% Ask for resampling method
if isempty(Method)
    % Ask method
    ind = java_dialog('radio', 'Select the resampling method:', 'Resample surface', [], ...
                      {['Matlab''s reducepatch:
' ...
                        '   | - Inhomogeneous mesh: large faces at the top of the gyri
' ...
                        '   | - Keeps the atlases and the subjects co-registration'], ...
                       ['Matlab''s reducepatch + subdivide large faces:
' ...
                        '   | - The large faces at the top of the gyri are subdivided in three
' ...
                        '   | - Deletes the atlases and the subjects co-registration'], ...
                       ['iso2mesh/CGAL library:
' ...
                        '   | - Homogeneous mesh: all the faces have similar sizes
' ...
                        '   | - Deletes the atlases and the subject co-registration
' ...
                        '   | - If the downsample looks dark, right-click > Swap faces']}, 1);
%                        ['iso2mesh/CGAL + project on the original surface:
' ...
%                         '   | - Homogeneous mesh but possible topological problems
' ...
%                         '   | - Damages the atlases and the subject co-registration']}, 1);
    if isempty(ind)
        return
    end
    % Select corresponding method name
    switch (ind)
        case 1,  Method = 'reducepatch';
        case 2,  Method = 'reducepatch_subdiv';
        case 3,  Method = 'iso2mesh';
        case 4,  Method = 'iso2mesh_project';
    end
end

%% ===== PROCESS MULTIPLE FILES =====
if ~isempty(MultipleFiles)
    for i = 1:length(MultipleFiles)
        [NewTessFile, iSurface, I, J] = tess_downsize( MultipleFiles{i}, newNbVertices, Method );
    end
    return;
end

%% ===== LOAD FILE =====
% Progress bar
bst_progress('start', 'Resample surface', 'Loading file...');
% Load file
TessMat = in_tess_bst(TessFile);
% Prepare variables
TessMat.Faces    = double(TessMat.Faces);
TessMat.Vertices = double(TessMat.Vertices);
dsFactor = newNbVertices / size(TessMat.Vertices, 1);

%% ===== RESAMPLE =====
bst_progress('start', 'Resample surface', ['Resampling surface: ' TessMat.Comment '...']);
% Resampling methods
switch (Method)
    % ===== REDUCEPATCH =====
    % Matlab's reducepatch
    case 'reducepatch'
        % Reduce number of vertices
        [NewTessMat.Faces, NewTessMat.Vertices] = reducepatch(TessMat.Faces, TessMat.Vertices, dsFactor);
        % Find the vertices that were kept by reducepatch
        [tmp, I, J] = intersect(TessMat.Vertices, NewTessMat.Vertices, 'rows');
        % Re-order the vertices so that they are in the same order in the output surface
        [I, iSort] = sort(I);
        NewTessMat.Vertices = TessMat.Vertices(I,:);
        J = J(iSort);
        % Re-order the vertices in the faces
        iSortFaces(J) = 1:length(J);
        NewTessMat.Faces = iSortFaces(NewTessMat.Faces);
        MethodTag = '';
        % Set the
        J = (1:length(J))';

    % ===== REDUCEPATCH + SUBDIV =====
    % Reducepatch + subdivide the large faces into smaller triangles
    case 'reducepatch_subdiv'
        % Reduce number of vertices
        [NewTessMat.Faces, NewTessMat.Vertices] = reducepatch(TessMat.Faces, TessMat.Vertices, dsFactor * 0.94);
        % Find the vertices that were kept by reducepatch
        [tmp, I, J] = intersect(TessMat.Vertices, NewTessMat.Vertices, 'rows');

        % Progress bar
        bst_progress('start', 'Resample surface', 'Analyzing surface...');
        % Calulate face areas and perimeter
        FaceArea  = tess_area(NewTessMat.Vertices, NewTessMat.Faces);
        % Vertex connectivity, normals, Curvature
        VertConn    = tess_vertconn(NewTessMat.Vertices, NewTessMat.Faces);
        [VertNormals, FaceNormals] = tess_normals(NewTessMat.Vertices, NewTessMat.Faces, VertConn);
        Curvature   = tess_curvature(NewTessMat.Vertices, VertConn, VertNormals);
        % Get center of each face
        FaceCenter = (NewTessMat.Vertices(NewTessMat.Faces(:,1),:) + NewTessMat.Vertices(NewTessMat.Faces(:,2),:) + NewTessMat.Vertices(NewTessMat.Faces(:,3),:)) ./ 3;
        % Get center of mass of the vertices
        SurfCenter = mean(NewTessMat.Vertices, 1);
        % Get large faces to subdivide (perimeter or area)
        iBigFaces = find((sum(FaceNormals .* bst_bsxfun(@minus, FaceCenter, [0 0 SurfCenter(3)]), 2) > 0.04) & ...  % Faces pointing outwards (normal in the same direction as position vector)
                         (sum(Curvature(NewTessMat.Faces) > 0, 2) >= 2) & ...       % Curvature has to be > 0
                         (FaceArea > mean(FaceArea) + 1*std(FaceArea)));            % Face area threshold
        % If there are not enough points to add (white matter): perform search on all the surface
        if (length(iBigFaces) < .75 * (newNbVertices - size(NewTessMat.Vertices,1)))
            iBigFaces = find(FaceArea > mean(FaceArea) + 2.5 * std(FaceArea));
        end
        % Display message
        disp(sprintf('BST> Subdividing %d faces from the %d faces generated by reducepatch.', length(iBigFaces), length(NewTessMat.Faces)));

% figure;
        % Loop over each face
        iRmFaces = [];
        bst_progress('start', 'Resample surface', 'Subdividing large faces...', 1, length(iBigFaces));
        for i = 1:length(iBigFaces)
            bst_progress('inc', 1);
            % Get the face and, the positions of its vertices, and the center of the face
            f = NewTessMat.Faces(iBigFaces(i),:);
            v = NewTessMat.Vertices(f,:);
            c = mean(v,1);

            % === BOUNDING BOX ===
            % Get maximum distance to consider around the face
            dmax = 1.2 * max(sqrt(sum(bst_bsxfun(@minus, v, c) .^ 2, 2)));
            % Select the vertices of the high-res surface that in a small sphere around the center of the face
            iVertBox = find(sum(bst_bsxfun(@minus, TessMat.Vertices, c) .^ 2, 2) < dmax.^2);

%             % Display selection
%             cla; hold on;
%             plot3(TessMat.Vertices(iVertBox,1), TessMat.Vertices(iVertBox,2), TessMat.Vertices(iVertBox,3), '.', 'tag', 'ptri');
%             plot3(c(1), c(2), c(3), '*y');
%             patch('Vertices', v, 'Faces', [1 2 3], 'FaceColor', 'r');
%             axis vis3d equal; drawnow; rotate3d on;

            % Get the vertices for the target face in the hi-resolution surface
            s1 = find(I(J == f(1)) == iVertBox);
            s2 = find(I(J == f(2)) == iVertBox);
            s3 = find(I(J == f(3)) == iVertBox);
            % Error?
            if isempty(s1) || isempty(s2) || isempty(s3)
                disp(sprintf('BST> Cannot subdivide big face #%d (box too small), skipping...', i));
                continue;
            end
            % Get a subset of the vertex connectivity matrix
            boxVertConn = TessMat.VertConn(iVertBox, iVertBox);

            % === FIND PATH TO COMMON NODES ===
            % Expand areas around all the vertices until they all overlap
            iter_max = 20;
            iter = 1;
            sx = [];
            while isempty(sx) && (iter < iter_max)
                s1 = union(s1, find(any(boxVertConn(s1,:),1)));
                s2 = union(s2, find(any(boxVertConn(s2,:),1)));
                s3 = union(s3, find(any(boxVertConn(s3,:),1)));
                sx = intersect(intersect(s1, s2), s3);
                iter = iter + 1;
            end

            % Expand areas around all the vertices until they all overlap
            iter_max = 50;
            iter = 1;
            istop = 0;
            d1 = 0;
            d2 = 0;
            d3 = 0;
            while (istop < 2) && (iter <= iter_max)
                % Grow from vertex #1
                i1 = find(any(boxVertConn(s1,:),1));
                s1 = [s1, setdiff(i1,s1)];
                d1 = [d1, i1*0+iter];
                % Grow from vertex #2
                i2 = find(any(boxVertConn(s2,:),1));
                s2 = [s2, setdiff(i2,s2)];
                d2 = [d2, i2*0+iter];
                % Grow from vertex #1
                i3 = find(any(boxVertConn(s3,:),1));
                s3 = [s3, setdiff(i3,s3)];
                d3 = [d3, i1*0+iter];
                % If all the vertices are in the region: stop immediately
                if (length(s1) == length(iVertBox)) && (length(s2) == length(iVertBox)) && (length(s3) == length(iVertBox))
                    istop = 10;
                % Do one more iterations after all the vertices are identified
                elseif (istop > 0) || (all(ismember([s2 s3],s1)) && all(ismember([s1 s3],s2)) && all(ismember([s1 s2],s3)))
                    istop = istop + 1;
                else
                    iter = iter + 1;
                end
            end
            % If an error occured: skip face
            if (iter > iter_max)
                disp(sprintf('BST> Cannot subdivide big face #%d (more than %d nodes distance), skipping...', i, iter_max));
                continue;
            end
            % Take intersection of the three regions
            [sx,ix,jx] = intersect(s1, s2);
            d1 = d1(ix);
            d2 = d2(jx);
            [sx,ix,jx] = intersect(sx, s3);
            d1 = d1(ix);
            d2 = d2(ix);
            d3 = d3(jx);
            dx = d1 + d2 + d3;

%             delete(findobj(0, 'tag', 'ptri')); plot3(TessMat.Vertices(iVertBox(s1),1), TessMat.Vertices(iVertBox(s1),2), TessMat.Vertices(iVertBox(s1),3), '.g', 'tag', 'ptri');
%             delete(findobj(0, 'tag', 'ptri')); plot3(TessMat.Vertices(iVertBox(s2),1), TessMat.Vertices(iVertBox(s2),2), TessMat.Vertices(iVertBox(s2),3), '.b', 'tag', 'ptri');
%             delete(findobj(0, 'tag', 'ptri')); plot3(TessMat.Vertices(iVertBox(s3),1), TessMat.Vertices(iVertBox(s3),2), TessMat.Vertices(iVertBox(s3),3), '.y', 'tag', 'ptri');

            % === SELECT VERTICES INSIDE THE FACE ===
            % Convert sx back to full indices list
            sx = iVertBox(sx);
            % Keep only the ones that project on the face INSIDE the triangle
            isInside = bst_intriangle(v(1,:), v(2,:), v(3,:), TessMat.Vertices(sx,:));
            sx = sx(isInside);
            dx = dx(isInside);
            if isempty(sx)
                disp(sprintf('BST> Cannot subdivide big face #%d (no candidate in the triangle), skipping...', i));
                continue;
            end
            % plot3(TessMat.Vertices(sx,1), TessMat.Vertices(sx,2), TessMat.Vertices(sx,3), '.y');

            % === SELECT VERTICES INSIDE THE FACE ===
            % Keep the closest path to all the nodes
            pathLength = min(dx);
            sx = sx(dx <= pathLength + 2);
%             plot3(TessMat.Vertices(sx,1), TessMat.Vertices(sx,2), TessMat.Vertices(sx,3), '.y');

            % === KEEP THE MOST CENTRAL LOCATION ===
            % Find the closest to the face center
            [d,imin] = min(sqrt(sum(bst_bsxfun(@minus, TessMat.Vertices(sx,:), c) .^ 2, 2)));
            sx = sx(imin);
            % Make sure it is not already in the destination surface
            if ismember(sx, I)
                disp(sprintf('BST> Cannot subdivide big face #%d (vertex already selected), skipping...', i));
                continue;
            end
%             plot3(TessMat.Vertices(sx,1), TessMat.Vertices(sx,2), TessMat.Vertices(sx,3), 'og');

            % === ADD VERTEX ===
            % Add the vertex to the list of vertices
            NewTessMat.Vertices = [NewTessMat.Vertices; c];
            iVertNew = size(NewTessMat.Vertices,1);
            I(end+1) = sx;
            J(end+1) = iVertNew;
            % Add the three new faces to the new surface
            NewTessMat.Faces = [NewTessMat.Faces; ...
                                f(1), f(2), iVertNew; ...
                                f(1), iVertNew, f(3); ...
                                iVertNew, f(2), f(3)];
            iRmFaces(end+1) = iBigFaces(i);
        end
        % Verify the unicity of the vertex selection
        if (length(I) ~= length(unique(I)))
            disp('BST> Error: The same vertex was selected multiple times in the high-resolution brain.');
            disp('BST> Using the basic reducepatch results instead...');
            % Call reducepatch only
            [NewTessFile, iSurface, I, J] = tess_downsize( TessFile, newNbVertices, 'reducepatch');
            return;
        end
        % Remove the deleted faces
        NewTessMat.Faces(iRmFaces,:) = [];
        MethodTag = '_subdiv';
        I = [];
        J = [];

    % ===== ISO2MESH =====
    % Using iso2mesh toolbox: good surfaces with equal triangle sizes, but no
    % correspondence of vertices in the original surface, and impossible to reconstruct the info
    case 'iso2mesh'
        % Reduce number of vertices
        NewTessMat = iso2mesh_resample(TessMat, dsFactor);
        if isempty(NewTessMat)
            return;
        end
        % Do not return any correspondence with the original vertices
        MethodTag = '_iso2mesh';
        I = [];
        J = [];

    % ===== ISO2MESH + PROJECT =====
    % Using iso2mesh toolbox: good surfaces with equal triangle sizes, but no
    % correspondence of vertices in the original surface, and impossible to reconstruct the info
    case 'iso2mesh_project'
        % Reduce number of vertices
        NewTessMat = iso2mesh_resample(TessMat, dsFactor);
        if isempty(NewTessMat)
            return;
        end
        Vertices = NewTessMat.Vertices;
        Faces    = NewTessMat.Faces;
        % Progress bar
        bst_progress('start', 'Resample surface', 'Analyzing surface...');
        % Calculate new normals
        newVertNormals = tess_normals(Vertices, Faces);
        % Remove duplicate vertices for Delaunay tesselation
        [delaunayVert, iSrc, iDest] = unique(TessMat.Vertices, 'rows');
        % Get the nearest neighbors
        I = bst_nearest(delaunayVert, Vertices);
        % Convert back to initial indices
        I = iSrc(I);
        % Compute the scalar product of the norms between the nearest neighbors and all the original vertices
        nv = size(TessMat.Vertices, 1);
        P = newVertNormals(:,1) .* TessMat.VertNormals(I) + ...
            newVertNormals(:,2) .* TessMat.VertNormals(nv + I) + ...
            newVertNormals(:,3) .* TessMat.VertNormals(2*nv + I);
        % Values with negative P are most likely on the other side of the sulcus
        iErr = find(P < -0.6);
        iKeep = find(P >= -0.6);

        % Loop on those points to fix them one by one with a different neighbor
        bst_progress('start', 'Resample surface', 'Fixing surface...', 1, length(iErr));
        for i = 1:length(iErr)
            bst_progress('inc', 1);
            % Calculate the scalar product of the normals of the current vertex with all the original vertices
            prodNorm = sum(bst_bsxfun(@times, newVertNormals(iErr(i),:), TessMat.VertNormals), 2);
            % Get indices for which the scalar product is positive
            iProdOk = find(prodNorm > -0.6);
            % Remove the vertices that are already in the mesh
            iProdOk = intersect(iProdOk, iKeep);
            % Find the nearest neighbor
            [m,iFix] = min(sum(bst_bsxfun(@minus, Vertices(iErr(i),:), TessMat.Vertices(iProdOk,:)) .^ 2, 2));
            I(iErr(i)) = iProdOk(iFix);
            % Add the corrected vertex to the list of valid vertices
            iKeep(end+1) = iProdOk(iFix);
        end

        % Find repeated vertices
        if (length(unique(I)) ~= length(I))
            disp('BST> ERROR: Found some duplicated vertices. Surface topology is incorrect...');
        end
        % Replace vertices with their nearest neighbor in the original surface
        Vertices = TessMat.Vertices(I,:);
        % Output structure
        NewTessMat.Faces    = Faces;
        NewTessMat.Vertices = Vertices;
        MethodTag = '_iso2mesh_proj';
end

%% ===== REMOVE FOLDED FACES =====
% Find equal faces
tmpFaces = sort(NewTessMat.Faces, 2);
[tmpFaces, iFaces] = unique(tmpFaces, 'rows');
% If there are some folded faces: delete them
if (length(iFaces) ~= size(NewTessMat.Faces,1))
    iRmFaces = setdiff(1:size(NewTessMat.Faces,1), iFaces);
    NewTessMat.Faces(iRmFaces,:) = [];
end

%% ===== CREATE NEW SURFACE STRUCTURE =====
% Build new filename and Comment
[filepath, filebase, fileext] = bst_fileparts(file_fullpath(TessFile));
NewComment = TessMat.Comment;
% Remove previous '_nbvertV' tags from Comment field
if (NewComment(end) == 'V')
    iUnderscore = strfind(NewComment, '_');
    if isempty(iUnderscore)
        iUnderscore = strfind(NewComment, ' ');
    end
    if ~isempty(~iUnderscore)
        NewComment = NewComment(1:iUnderscore(end)-1);
    end
end
% Remove previous '_nbvertV' tags from filename
if (filebase(end) == 'V')
    iUnderscore = strfind(filebase, '_');
    filebase = filebase(1:iUnderscore(end)-1);
end
% Add a '_nbvertV' tag
NewTessFile = file_unique(bst_fullfile(filepath, sprintf('%s_%dV%s', filebase, newNbVertices, fileext)));
NewComment  = sprintf('%s%s_%dV', NewComment, MethodTag, size(NewTessMat.Vertices,1));

% As per MMII convention - there should be one tessellation file per envelope
% A downsized version of e.g. a cortex is considered as a different envelope
% ans is therefore saved in a separate tessellation file than the original.
NewTessMat.Comment  = NewComment;
% Copy history field
if isfield(TessMat, 'History')
    NewTessMat.History = TessMat.History;
end
% History: Downsample surface
NewTessMat = bst_history('add', NewTessMat, 'downsample', sprintf('Downsample surface: %d -> %d vertices', oldNbVertices, newNbVertices));

%% ===== DOWNSAMPLE SCOUTS =====
% Existing atlases
if isfield(TessMat, 'Atlas') && ~isempty(TessMat.Atlas) && ~isempty(I)
    % Copy scout structure
    NewTessMat.Atlas = TessMat.Atlas;
    % Loop on all the scouts, and keep only those vertices
    for iAtlas = 1:length(NewTessMat.Atlas)
        iRmScout = [];
        for iScout = 1:length(NewTessMat.Atlas(iAtlas).Scouts)
            % Replace the old vertices index with the new ones
            [a,b,c] = intersect(NewTessMat.Atlas(iAtlas).Scouts(iScout).Vertices, I);
            NewTessMat.Atlas(iAtlas).Scouts(iScout).Vertices = reshape(sort(J(c)), 1, []);
            % If scout has no vertex left: tag for deletion
            if isempty(NewTessMat.Atlas(iAtlas).Scouts(iScout).Vertices)
                iRmScout(end+1) = iScout;
            end
        end
        % Remove empty scouts
        if ~isempty(iRmScout)
            NewTessMat.Atlas(iAtlas).Scouts(iRmScout) = [];
        end
        % Set scouts seeds
        NewTessMat.Atlas(iAtlas).Scouts = panel_scout('SetScoutsSeed', NewTessMat.Atlas(iAtlas).Scouts, NewTessMat.Vertices);
    end
end
% Selected atlas
if isfield(TessMat, 'iAtlas') && ~isempty(TessMat.iAtlas)
    NewTessMat.iAtlas = TessMat.iAtlas;
end

%% ===== DOWNSAMPLE REGISTRATION MAPS =====
% FreeSurfer spheres
if isfield(TessMat, 'Reg') && isfield(TessMat.Reg, 'Sphere') && isfield(TessMat.Reg.Sphere, 'Vertices') && ~isempty(TessMat.Vertices) && (length(TessMat.Reg.Sphere.Vertices) == length(TessMat.Vertices))
    % Keep only the selected indices
    if ~isempty(I)
        newSphVert = TessMat.Reg.Sphere.Vertices(I,:);
        NewTessMat.Reg.Sphere.Vertices = newSphVert;
    else
        NewTessMat.Reg.Sphere = [];
    end
end

if isfield(TessMat, 'Reg') && isfield(TessMat.Reg, 'SphereLR') && isfield(TessMat.Reg.SphereLR, 'Vertices') && ~isempty(TessMat.Vertices) && (length(TessMat.Reg.SphereLR.Vertices) == length(TessMat.Vertices))
    % Keep only the selected indices
    if ~isempty(I)
        newSphVert = TessMat.Reg.SphereLR.Vertices(I,:);
        NewTessMat.Reg.SphereLR.Vertices = newSphVert;
    else
        NewTessMat.Reg.SphereLR = [];
    end
end

% BrainSuite squares
if isfield(TessMat, 'Reg') && isfield(TessMat.Reg, 'Square') && isfield(TessMat.Reg.Square, 'Vertices') && ~isempty(TessMat.Reg.Square.Vertices) && (length(TessMat.Reg.Square.Vertices) == length(TessMat.Vertices))
    % Keep only the selected indices
    if ~isempty(I)
        newSqVert = TessMat.Reg.Square.Vertices(I,:);
        NewTessMat.Reg.Square.Vertices = newSqVert;
    else
        NewTessMat.Reg.Square = [];
    end
    NewTessMat.Reg.AtlasSquare=TessMat.Reg.AtlasSquare;
end

%% ===== UPDATE DATABASE =====
% Save downsized surface file
bst_save(NewTessFile, NewTessMat, 'v7');
% Make output filename relative
NewTessFile = file_short(NewTessFile);
% Get subject
[sSubject, iSubject] = bst_get('SurfaceFile', TessFile);
% Register this file in Brainstorm database
iSurface = db_add_surface(iSubject, NewTessFile, NewComment);

% Close progress bar
bst_progress('stop');

end

%% ===== iso2mesh_resample =====
% Resample a surface using iso2mesh/CGAL library
% Author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
function NewTessMat = iso2mesh_resample(TessMat, dsFactor)
    % Check if iso2mesh is installed
    if ~exist('meshresample', 'file')
        bst_error(['Please install iso2mesh on your computer:

' ...
               '  1) Visit the website: http://iso2mesh.sourceforge.net
' ...
               '  2) Download the iso2mesh package for your operating system.
' ...
               '  3) Unzip it on your local hard drive.
' ...
               '  4) Add the iso2mesh folder to your Matlab path.
' ...
               '  5) Try again downsampling the surface.

'], 'Install iso2mesh', 0);
        web('http://sourceforge.net/projects/iso2mesh/files/iso2mesh/1.5.0%20%28iso2mesh%202013%29/', '-browser');
        NewTessMat = [];
        return;
    end
    % Running iso2mesh routine
    [Vertices,Faces] = meshresample(TessMat.Vertices, TessMat.Faces, dsFactor);
    % Error handling
    if isempty(Vertices) || isempty(Faces)
        error(['Iso2mesh failed downsampling this surface:' 10 'See Matlab command window for more information.']);
    end
    % Report results
    NewTessMat.Vertices = Vertices;
    NewTessMat.Faces    = Faces;
    % Swap faces
    % NewTessMat.Faces = NewTessMat.Faces(:,[2 1 3]);
end

% OLD VERSION: FUNCTION cgalsimp2 WAS INCLUDED IN BRAINSTORM DISTRIBUTION
% %% ===== iso2mesh_resample =====
% % Resample a surface using iso2mesh/CGAL library
% % Author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% function NewTessMat = iso2mesh_resample(TessMat, dsFactor, nCall)
%     % First call
%     if (nargin < 3) || isempty(nCall)
%         nCall = 1;
%     end
%     % Get the executable name
%     switch(bst_get('OsType'))
%         case {'linux32', 'linux64'}, exePath = 'cgalsimp2.mexglx';
%         case 'mac32',                exePath = 'cgalsimp2.mexmaci';
%         case 'mac64',                exePath = 'cgalsimp2.mexmaci64';
%         case {'win32', 'win64'},     exePath = 'cgalsimp2.exe';
%         otherwise, error('CGAL executable is not available on your OS.');
%     end
%     % Add the full path
%     exePath = bst_fullfile(bst_get('BrainstormHomeDir'), 'external', 'iso2mesh', exePath);
%     % Get temporary mesh file
%     fin  = file_unique(bst_fullfile(bst_get('BrainstormTmpDir'), 'mesh_in.off'));
%     fout = file_unique(bst_fullfile(bst_get('BrainstormTmpDir'), 'mesh_out.off'));
%     % Save input file
%     out_tess_off(TessMat, fin);
%     % Execute cgalsimp2 with a system call
%     [status, result] = system(['"' exePath '" "' fin '" ' num2str(dsFactor) ' "' fout '"']);
%     if status
%         file_delete(fin, 1);
%         error(['CGAL failed downsampling this surface:' 10 result]);
%     end
%     % Read results
%     NewTessMat = in_tess_off(fout);
%     % Delete mesh files
%     file_delete({fin, fout}, 1);
%     % If no results are produced: Fix the surface and call it again
%     if isempty(NewTessMat.Vertices)
%         % If it is already the second call
%         if (nCall > 2)
%             error('CGAL failed downsampling this surface.');
%         end
%         % Remove duplicate faces
%         TessMat.Faces = removedupelem(TessMat.Faces);
%         % Remove isolated nodes
%         [TessMat.Vertices, TessMat.Faces] = removeisolatednode(TessMat.Vertices, TessMat.Faces);
%         % Run again the function
%         NewTessMat = iso2mesh_resample(TessMat, dsFactor, nCall + 1);
%     end
% end


text:
function SaveNet(itr, useHourCost, sigGraph, learner, LOG, para)

if isfield(para, 'saveModelEveryXIter') && mod(itr,para.saveModelEveryXIter)~=0
    return;
end

modelfile = [para.output sprintf('.itr%d.hour%d.LR%s', itr, round(LOG.nHourSeen), FormatFloat4Name(learner.learningRate))];
if useHourCost
    cost_cv = LOG.cost_cv_hour.taskCost(end);
else
    cost_cv = LOG.cost_cv.taskCost(end);
end
if cost_cv>0.1
    modelfile = sprintf('%s.CV%2.3f.mat', modelfile, cost_cv);
else
    modelfile = sprintf('%s.CV%s.mat', modelfile, FormatFloat4Name(cost_cv));
end

sigGraph.useGPU = 0;
sigGraph = sigGraph.setPrecision();
sigGraph = sigGraph.cleanUp();
save(modelfile, 'sigGraph', 'para', 'learner', 'LOG');

end

author: singaxiong

repo: SignalGraph


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [yc,y1,y2,y3] = grid2grid(yc,m,in,out)
%
% transfers input grid (c,s,n) to output grid (c,s,n) using inter- and extrapolation
% central operation are reduce/expand, acting on each dimension individually
%
% Input:
%   yc        input grid
%   m         number of discretization points
%   in        type of input grid
%   out       type of output grid
% Output:
%   yc        output grid
%   y1,y2,y3  components of yc
%
%==============================================================================

function [yc,y1,y2,y3] = grid2grid(yc,m,in,out)

if nargin == 0, % help and minimal example
  runMinmalExample;
  return;
end;

dim = length(m);  y{1} = []; y{3} = []; y{3} = [];
if ~exist('out','var'), out = in;  end;

switch in,
  case {'centered','cell-centered'}, n = prod(m);
    % decompose Yin
    for j=1:dim, y{j} = reshape(yc((j-1)*n+(1:n)),m); end;
    switch out,
      case {'centered','cell-centered'},
      case 'staggered', for j=1:dim, y{j} = extend(y{j},j); end;
      case 'nodal',
        for j=1:dim,
          for k=1:dim, y{j} = extend(y{j},k); end;
        end;
    end;

  case 'staggered',
    % decompose Yin, get dimensions right
    E = eye(dim) + m'*ones(1,dim); n = cumsum([0,prod(E)]);
    for j=1:dim, y{j} = reshape(yc(n(j)+1:n(j+1)),E(:,j)'); end;
    switch out,
      case {'centered','cell-centered'}, for j=1:dim, y{j} = reduce(y{j},j); end;
      case 'staggered',
      case 'nodal',
        for j=1:dim,
          for k=setdiff(1:dim,j), y{j} = extend(y{j},k); end;
        end;
    end;

  case 'nodal', n = prod(m+1);
    % decompose Yin
    for j=1:dim, y{j} = reshape(yc((j-1)*n+(1:n)),m+1); end;
    switch out,
      case {'centered','cell-centered'},
        for j=1:dim, for k=1:dim, y{j} = reduce(y{j},k); end; end;
      case 'staggered',
        for j=1:dim, for k=setdiff(1:dim,j) y{j} = reduce(y{j},k); end; end;
      case 'nodal',
    end;

  otherwise,
    error('can not interpret flags')
end;
y1 = y{1}; y2 = y{2}; y3 = y{3}; yc = [y1(:);y2(:);y3(:)];

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

% the following operators act on x=y{k} in the j-th direction,
% permuting x makes j the first direction, averaging (extend creates two
% additional outside points based on linear BC) and permute back

function x = reduce(x,j)
m = size(x); J = [j,setdiff(1:length(size(x)),j)];
x = permute(x,J); x = (x(1:end-1,:,:)+x(2:end,:,:))/2; x = ipermute(x,J);

function x = extend(x,j)
m = size(x); J = [j,setdiff(1:length(size(x)),j)]; x = permute(x,J);
x = [1.5*x(1,:,:)-0.5*x(2,:,:);(x(1:end-1,:,:)+x(2:end,:,:))/2;...
  1.5*x(end,:,:)-0.5*x(end-1,:,:)];
x = ipermute(x,J);

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

function runMinmalExample
help(mfilename);
omega = [0,6,0,4]; m = [4 3];

h   = (omega(2:2:end)-omega(1:2:end))./m;
eps = 0.2*min(h);

yC  = getCellCenteredGrid(omega,m); %yC = yC + eps*randn(size(yC));
yN  = reshape(getNodalGrid(omega,m),[],2);
yS  = getStaggeredGrid(omega,m);

yCN = grid2grid(yC,m,'centered','nodal');
yCS = grid2grid(yC,m,'centered','staggered');
ySN = grid2grid(yS,m,'staggered','nodal');
ySC = grid2grid(yS,m,'staggered','centered');
yNC = grid2grid(yN,m,'nodal','centered');
yNS = grid2grid(yN,m,'nodal','staggered');

yC  = reshape(yC, [],2);
yN  = reshape(yN, [],2);
yCN = reshape(yCN,[],2);
ySN = reshape(ySN,[],2);
ySC = reshape(ySC,[],2);
yNC = reshape(yNC,[],2);

FAIRfigure(1); clf;

subplot(3,2,1);
plotGrid(yN,omega,m,'color','g'); hold on;
plot(yC(:,1),yC(:,2),'bs');
plot(yCN(:,1),yCN(:,2),'m*');
title(sprintf('%s: example cell-centered to nodel',mfilename));

subplot(3,2,2);
plotGrid(yN,omega,m,'color','g'); hold on;
plot(yC(:,1),yC(:,2),'bs');
[y11,y12,y21,y22] = splittStaggeredGrid(yCS,omega,m);
plot(y11,y12,'m>',y21,y22,'m^');
title(sprintf('%s: example cell-centered to staggered',mfilename));

[y11,y12,y21,y22] = splittStaggeredGrid(yS,omega,m);
subplot(3,2,3);
plotGrid(yN,omega,m,'color','g'); hold on;
plot(y11,y12,'b>',y21,y22,'b^');
plot(ySC(:,1),ySC(:,2),'ms');
title(sprintf('%s: example staggered to cell-centered',mfilename));

subplot(3,2,4);
plotGrid(yN,omega,m,'color','g'); hold on;
plot(y11,y12,'b>',y21,y22,'b^');
plot(ySN(:,1),ySN(:,2),'m*');
title(sprintf('%s: example staggered to nodal',mfilename));

subplot(3,2,5);
plotGrid(yN,omega,m,'color','g'); hold on;
plot(yN(:,1),yN(:,2),'b*');
plot(yNC(:,1),yNC(:,2),'ms');
title(sprintf('%s: example nodal to cell-centered',mfilename));

subplot(3,2,6);
plotGrid(yN,omega,m,'color','g'); hold on;
plot(yN(:,1),yN(:,2),'b*');
[y11,y12,y21,y22] = splittStaggeredGrid(yNS,omega,m);
plot(y11,y12,'m>',y21,y22,'m^');

title(sprintf('%s: example nodal to staggered',mfilename));

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

function [y11,y12,y21,y22] = splittStaggeredGrid(yS,omega,m);
dim = length(omega)/2;
e  = @(i) (1:dim == i); % i-th unit vector
ns = cumsum([0;prod(ones(dim,1)*m+eye(dim),2)]);
yN = reshape(grid2grid(yS,m,'staggered','nodal'),[m+1,2]);
y11 = reshape(yS(ns(1)+1:ns(2)),m+e(1));
y12 = (yN(:,1:end-1,2)+yN(:,2:end,2))/2;
y21 = (yN(1:end-1,:,1)+yN(2:end,:,1))/2;
y22 = reshape(yS(ns(2)+1:ns(3)),m+e(2));

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


text:
function s = s_cap ( s )

%*****************************************************************************80
%
%% S_CAP replaces any lowercase letters by uppercase ones in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character S(*), the string to be transformed.
%
%    Output, character S(*), the capitalized string.
%
  nchar = s_len_trim ( s );

  for i = 1 : nchar

    s(i) = ch_cap ( s(i) );

  end

  return
end


text:
function prob = ml_gaussPDF(Data, Mu, Sigma)
%ML_GAUSSPDF
% This function computes the Probability Density Function (PDF) of a
% multivariate Gaussian represented by means and covariance matrix.
%
% Author:	Sylvain Calinon, 2009
%			http://programming-by-demonstration.org
%
% Inputs -----------------------------------------------------------------
%   o Data:  D x N array representing N datapoints of D dimensions.
%   o Mu:    D x 1 array representing the centers of the K GMM components.
%   o Sigma: D x D x 1 array representing the covariance matrices of the
%            K GMM components.
% Outputs ----------------------------------------------------------------
%   o prob:  1 x N array representing the probabilities for the
%            N datapoints.

[nbVar,nbData] = size(Data);

%      (D x N) - repmat((D x 1),1,N)
%      (D x N) - (D x N)
Data = (Data - repmat(Mu,1,nbData))';

% Data = (N x D)

% (N x 1)
prob = sum((Data*inv(Sigma)).*Data, 2);
prob = exp(-0.5*prob) / sqrt((2*pi)^nbVar * (abs(det(Sigma))+realmin));

author: epfl-lasa


text:
%
% [fh, varfh] = gppred(X, f, Covf, Xh, kX, logtheta, ...)
%
% p(y|X) = N(y; f(X), inv(V))
function [f, varf] = spgpr(logtheta, covfunc, X, y, V, varargin)
%function [fh, varfh] = gppred(X, f, Covf, Xh, kX, logtheta, varargin)

% If you need more than just variances for PCA, you could give some sparse
% matrix to indicate which covariances should be evaluated. But again, this
% covariances should be zero for the latent function values fh a
% priori.. Things would be much easier if I just factored the principal
% components too.. :)

opts = struct('cholkx', [], 'khx', [], 'kh', []);
%opts = struct('cholkx', [], 'khx', [], 'kh', []);

[ opts, errmsg, wrnmsg ] = argschk( opts, varargin{:} );
if ~isempty(errmsg), error( errmsg ), end
if ~isempty(wrnmsg), warning( wrnmsg ), end

nh = cols(Xh);

if isempty(opts.cholkx)
  Kpp = feval(kX, X, X, logtheta);
  Kpp = Kpp + 1e-6*eye(size(Kpp));
  Lp = chol(Kpp, 'lower');
else
  Lp = opts.cholkx;
end

if isempty(opts.khx)
  Kxp = feval(kX, Xh, X, logtheta);
else
  Kxp = opts.khx;
end

if isempty(opts.kh)
  Kx = feval(kX, Xh, [], logtheta);
else
  Kx = opts.kh;
end

fh = Kxp * solve_triu(Lp', solve_tril(Lp, f));

if nargout >= 2
  varfh = zeros(nh,1);
  for i=1:nh
    r = solve_tril(Lp, Kxp(i,:)');
    s = solve_triu(Lp', r);
    varfh(i) = Kx(i) - r'*r + s'*Covf*s;
  end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function E=TeagerEnergy(x,dim)
%%TEAGERENERGY Compute the Teager energy of a discrete signal as defined in
%              [1]. The Teager energy is a type of instantaneous "energy"
%              signal that produces a result related to the rate of change
%              of the amplitude and frequency of the original signal. It
%              has found use in a number of applications, primarily related
%              to acoustic processing, such as in the wavelet-based speech
%              enhancement of [2].
%
%INPUTS: x A vector signal or a matrix of signals whose Teager energy is
%          desired.
%      dim The dimension of x over which the Teager energy is to be found.
%          For example, for a NXnumSigs set of signals, one would use
%          dims=1 as each column is a separate signal.
%
%OUTPUTS: E The Teager energy. This has the same size as x.
%
%EXAMPLE 1:
%Here we compute the Teager energy signal of an up-chirp followeb by a down
%chirp.
% fStart=10;
% fEnd=100;
% T=1;
%
% %We want the chirp to be sampled at 100 times the Nyquist frequency.
% T0=1/(100*2*fEnd);%The sampling period.
% %We want to sample just the duration of the signal. Thus, time goes from
% %T=0 to T=1. This means that we need N=T/T0 samples.
% N=T/T0;
%
% %Make a chirp up and then a chirp down
% [x1,t]=LFMChirp(T,fStart,fEnd,T0);
% x2=LFMChirp(T,fEnd,fStart,T0);
% x=real([x1(:);x2(:)]);
% t=[t(:);t(end)+t(2)-t(1)+t(:)];
%
% E=TeagerEnergy(x);
%
% figure(1)
% clf
% subplot(2,1,1)
% hold on
% plot(t,x,'-r','linewidth',2)
% title('Original Signal')
% subplot(2,1,2)
% plot(t,E,'-r','linewidth',2)
% title('Transformed Signal')
%
%EXAMPLE 2:
%Here we compute the Teager energy of a sum of sinusoids. It essentially
%just follows the beat pattern.
% n=(0:100).';
% x=sin(pi/6*n)+sin(pi/4*n);
%
% E=TeagerEnergy(x);
%
% figure(2)
% clf
% subplot(2,1,1)
% hold on
% plot(n,x,'-r','linewidth',2)
% title('Original Signal')
% subplot(2,1,2)
% plot(n,E,'-r','linewidth',2)
% title('Transformed Signal')
%
%REFERENCES:
%[1] J. F. Kaiser, "On a simple algorithm to calculate the 'energy' of a
%    signal," in Proceedings of the International Conference on Acoustics,
%    Speech and Signal Processing, Albuquerque, NM, 3-6 Apr. 1990, pp.
%    381-384.
%[2] M. Bahoura and J. Rouat, "Wavelet speech enhancement based on the
%    Teager energy operator," IEEE Signal Processing Letters, vol. 8, no.
%    1, pp. 10-12, Jan. 2001.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(dim))
    E=x.^2-circshift(x,-1).*circshift(x,1);
else
    E=x.^2-circshift(x,-1,dim).*circshift(x,1,dim);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function b = InitBatch(bfile)

%InitBatch - Initialize a new batch job.
%
%  USAGE
%
%    b = InitBatch(bfile)
%
%    bfile          batch file listing the parameters for each iteration
%
%  SEE
%
%    See also StartBatch, ShowBatch
%

% Copyright (C) 2007-2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Check number of parameters
if nargin < 1,
	error('Incorrect number of parameters (type ''help InitBatch'' for details).');
end

% Make sure batch file is valid
s = dbstack;
if ~isstring_FMAT(bfile) || ~exist(bfile,'file'),
	error(['Batch file not found (type ''help ' s.name ''' for details).']);
end

% Open batch file
f = fopen(bfile,'r');
if f == -1, error(['Could not open file ''' bfile '''.']); end

item = 1;
l = 1;
nFields = 0;
% Read batch file
while ~feof(f),
	field = 0;
	% Read next line
	line = fgetl(f);
	while ~isempty(line),
		% Get next field
		[token,line] = strtok(line);
		% Skip rest of line if this is a comment
		if isempty(token) | token(1) == '%', break; end
		field = field + 1;
		% Determine if this is a number, vector or string
		n = str2num(token);
		if ~isempty(n),
			% It is a number, convert it to numerical format
			b.field{item,field} = n;
		elseif token(1) ~= '[',
			% It is a string, keep it as it is
			b.field{item,field} = token;
		else
			% It is a vector (or matrix)
			[token,line] = strtok([token line],'[]');
			n = str2num(token);
			if isempty(n),
				error(['Incorrect vector or matrix in batch file (line ' int2str(l) ', element ' int2str(field) ').']);
			end
			b.field{item,field} = n;
			line = line(2:end);
		end
	end
	if field > 0,
		% Make sure all lines have the same number of elements
		if item == 1,
			nFields = field;
		elseif field ~= nFields,
			if l == 2,
				error(['Incoherent numbers of elements in batch file (' int2str(nFields)  ' elements on line 1, but ' int2str(field) ' elements on line 2).']);
			else
				error(['Incoherent numbers of elements in batch file (' int2str(nFields)  ' elements on lines 1-' int2str(l-1) ', but ' int2str(field) ' elements on line ' int2str(l) ').']);
			end
		end
		% Next item
		b.line(item) = l;
		item = item + 1;
	end
	% Update line counter (for error messages)
	l = l + 1;
end

% Close file
fclose(f);

% Reset iterator
b.currentItem = 0;
b.currentField = 0;

% Empty fields
b.log = -1;
b.mfile = [];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function dataset = voc2012_test(dataset, usage, use_flip)
% Pascal voc 2012 test set
% set opts.imdb_train opts.roidb_train
% or set opts.imdb_test opts.roidb_train

% change to point to your devkit install
devkit                      = voc2012_devkit();

switch usage
    case {'train'}
        dataset.imdb_train    = {  imdb_from_voc(devkit, 'test', '2012', use_flip) };
        dataset.roidb_train   = cellfun(@(x) x.roidb_func(x), dataset.imdb_train, 'UniformOutput', false);
    case {'test'}
        dataset.imdb_test     = imdb_from_voc(devkit, 'test', '2012', use_flip) ;
        dataset.roidb_test    = dataset.imdb_test.roidb_func(dataset.imdb_test);
    otherwise
        error('usage = ''train'' or ''test''');
end

end

author: ShaoqingRen


text:
Network vgg16 {
Layer CONV1 {
Type: CONV
Dimensions { K 64,C 3,R 3,S 3,Y 224,X 224 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
}
}

Layer CONV2 {
Type: CONV
Dimensions { K 64,C 64,R 3,S 3,Y 224,X 224 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
}
}

Layer CONV3 {
Type: CONV
Dimensions { K 128,C 64,R 3,S 3,Y 112,X 112 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
}
}

Layer CONV4 {
Type: CONV
Dimensions { K 128,C 128,R 3,S 3,Y 112,X 112 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV5 {
Type: CONV
Dimensions { K 256,C 128,R 3,S 3,Y 56,X 56 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV6 {
Type: CONV
Dimensions { K 256,C 256,R 3,S 3,Y 56,X 56 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV7 {
Type: CONV
Dimensions { K 256,C 256,R 3,S 3,Y 56,X 56 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV8 {
Type: CONV
Dimensions { K 512,C 256,R 3,S 3,Y 28,X 28 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV9 {
Type: CONV
Dimensions { K 512,C 512,R 3,S 3,Y 28,X 28 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV10 {
Type: CONV
Dimensions { K 512,C 512,R 3,S 3,Y 28,X 28 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV11 {
Type: CONV
Dimensions { K 512,C 512,R 3,S 3,Y 14,X 14 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV12 {
Type: CONV
Dimensions { K 512,C 512,R 3,S 3,Y 14,X 14 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

Layer CONV13 {
Type: CONV
Dimensions { K 512,C 512,R 3,S 3,Y 14,X 14 }
Dataflow {
			SpatialMap(1,1) C;
			TemporalMap(64,64) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
			TemporalMap(1,1) Y';
			TemporalMap(1,1) X';
			Cluster(64, P);
			SpatialMap(1,1) K;
			TemporalMap(3,3) R;
			TemporalMap(3,3) S;
		}
}

}

author: maestro-project

repo: maestro


text:
%CURLABLIST Get current label list
%
% [LABLISTNUMBER,LABLISTNAME,T0,T1] = CURLABLIST(A)
%
% INPUT
%   A             - Dataset
%
% OUTPUT
%   LABLISTNUMBER - Number of current lablist
%   LABLISTNAME   - Name of current label list
%		T0            - Start column current targets
%   T1            - End column current targets
%
% DESCRIPTION
% The number and name of the current label list are returned.
% If T1<T0: no tergets (or soft labels) set for current label list.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MULTI_LABELING, ADDLABLIST, CHANGELABLIST

author: marianux

repo: ecg-kit

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


text:
d = 500;
n = 20;
C = 100;
k = 10;
q = 50;
q2 = 50;

ODLSI = C*k*(k*d+d*n + q*k*n) + C*q*k*d^3
EDLSI = C*k*(k*d+d*n + q*k*n) + C*d^3 + C*q*d*k*(q*k + d)
OFDDL = C^2*d*k*(n+C*k+C*n)+C*k^2*q*(d + C^2*n)
EFDDL = C^2*k*((q+1)*k*(d+C*n) + 2*d*n)

OCOPAR = C^3*k^2*(2*d + C*k + q*n) + C*q*k*d^3
ECOPAR = C^3*k^2*(2*d + C*k + q*n) + C*d^3 + C*q*d*k*(q*k + d)

LRSDL = C^2*k*((q+1)*k*(d + C*n) + 2*d*n) + C^2*d*k*n + (q + q2)*d*k^2

%%
ODLSID = C*q*k*d^3
EDLSID = C*d^3 + C*q*d*k*(q*k + k)
OFDDLX = C^2*k*(d*n + q*C*k*n + C*d*k)
EFDDLX = C^2*k*(d*n + q*C *n*k + d*k)
OFDDLD = C*d*k*(q*k + C^2*n)
EFDDLD = C*d*k*(C*n + C*q*k) + C^3*k^2 *n

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%PlotBrowser
% Creates GUI window which allows the user to quickly view plots of
% simulation. User can either select plots from list box, or scroll through
% plots using the previous and next page buttons.
%
%Uses meta class to find names of all public static  plot methods in
%analysis classes. Their signature should be:
%   function plot*(sim, axes, time, compartments)
%   end
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 8/6/2011
classdef PlotBrowser < handle
    properties (SetAccess = protected)
        simulation

        plotNames
        plotFunctions

        figure
        buttonsPanel
        compartmentListBoxTitle
        compartmentListBox
        plotListBoxTitle
        plotListBox
        plotButton
        prevButton
        nextButton
        saveButton
        printButton
        plotPanel
    end

    properties (Constant = true)
        plotClasses = {
            'Chromosome'
            'MetabolicReaction'
            'Transcription'
            'Translation'
            }
    end

    methods
        function this = PlotBrowser(sim, selectedPlotNames, selectedCompartmentNames)
            this = this@handle;
            this.simulation = sim;

            %find plots
            this.discoverPlots();

            %open window
            this.open();

            %layout window fullscreen
            screenSize = get(0, 'ScreenSize');
            screenSize(2) = 30;
            screenSize(4) = screenSize(4) - screenSize(2);
            screenSize([3 4]) = max(1, screenSize([3 4]));
            set(this.figure, 'OuterPosition', screenSize);

            %plot
            if nargin < 2
                selectedPlotNames = this.plotNames{1};
            end
            if nargin < 3
                selectedCompartmentNames = 'Sum';
            end
            this.selectPlots(selectedPlotNames, selectedCompartmentNames);
        end
    end

    methods
        function selectPlots(this, selectedPlotNames, selectedCompartmentNames)
            %select plots
            [~, selectedPlots] = ismember(selectedPlotNames, this.plotNames);
            set(this.plotListBox, 'Value', selectedPlots);

            %select compartments
            if nargin < 3
                selectedCompartmentNames = 'Sum';
            end
            [~, selectedCompartments] = ismember(selectedCompartmentNames, get(this.compartmentListBox, 'String'));
            set(this.compartmentListBox, 'Value', selectedCompartments);

            %plot
            this.graphSelectedPlots();
        end

        function open(this)
            %figure
            this.figure = figure('Name', 'Whole Cell Simulation :: Plot Browser', ...
                'NumberTitle', 'off', 'ResizeFcn', @this.layoutWindow); %#ok<CPROP,PROP>

            %list box panel
            this.buttonsPanel = uipanel('Title', 'Configure Plots');

            %select compartment
            this.compartmentListBoxTitle = uicontrol(this.buttonsPanel, 'Style', 'text',...
                'String', 'Select Compartments', 'HorizontalAlignment', 'Center');
            this.compartmentListBox = uicontrol(this.buttonsPanel, 'Style', 'listbox',...
                'String', [{'Sum'; 'Mean'; 'Median'; 'Std Dev'}; this.simulation.compartment.names],...
                'Min', 1, 'Max', 4 + this.simulation.compartment.count, 'Value', 1);

            %plot list box
            this.plotListBoxTitle = uicontrol(this.buttonsPanel,'Style', 'text',...
                'String','Select Plots', 'HorizontalAlignment', 'Center');
            this.plotListBox = uicontrol(this.buttonsPanel, 'Style', 'listbox',...
                'String', this.plotNames, 'Min', 1, 'Max', length(this.plotNames), 'Value', 1);

            %plot button
            this.plotButton = uicontrol(this.buttonsPanel, 'Style', 'pushbutton',...
                'String', 'Plot', 'Callback', @this.graphSelectedPlots);

            %previous page button
            this.prevButton = uicontrol(this.buttonsPanel, 'Style', 'pushbutton',...
                'String', 'Prev', 'Callback', @this.graphPrevPlots);

            %next page button
            this.nextButton = uicontrol(this.buttonsPanel, 'Style', 'pushbutton',...
                'String', 'Next', 'Callback', @this.graphNextPlots);

            %save button
            this.saveButton = uicontrol(this.buttonsPanel, 'Style', 'pushbutton',...
                'String', 'Save', 'Callback', @this.save);

            %print button
            this.printButton = uicontrol(this.buttonsPanel, 'Style', 'pushbutton',...
                'String', 'Print', 'Callback', @this.print);

            %plot panel
            this.plotPanel = uipanel('Title','Plots');
        end

        function close(this)
            close(this.figure);

            this.figure = [];
            this.buttonsPanel = [];
            this.compartmentListBoxTitle = [];
            this.compartmentListBox = [];
            this.plotListBoxTitle = [];
            this.plotListBox = [];
            this.plotButton = [];
            this.prevButton = [];
            this.nextButton = [];
            this.saveButton = [];
            this.printButton = [];
            this.plotPanel = [];
        end
    end

    %helper methods
    methods (Access = protected)
        %Return list of names of plotting functions
        function discoverPlots(this)
            this.plotNames = cell(0, 1);
            this.plotFunctions = cell(0, 1);

            metaPackage = meta.package.fromName('edu.stanford.covert.cell.sim.analysis');
            for i = 1:numel(metaPackage.Classes)
                metaClass = metaPackage.Classes{i};

                %ignore test clases, handle clases
                if any(ismember(cellfun(@(class) class.Name, metaClass.SuperClasses, 'UniformOutput', false), {'TestCase', 'handle'}))
                    continue;
                end

                %ignore unsupported classes
                className = metaClass.Name(find(metaClass.Name == '.', 1, 'last') + 1:end);
                if ~ismember(className, this.plotClasses)
                    continue;
                end

                %find plotting methods
                for j = 1:numel(metaClass.Methods)
                    metaMethod = metaClass.Methods{j};

                    %filter out methods to match required signature
                    if ...
                            length(metaMethod.Name) < 5 || ...
                            ~strcmp(metaMethod.Name(1:4), 'plot') || ...
                            ~metaMethod.Static || ...
                            ~strcmp(metaMethod.Access, 'public') || ...
                            length(metaMethod.InputNames) < 4
                        continue;
                    end

                    this.plotNames = [this.plotNames; className ':' strrep(metaMethod.Name(5:end), '_', ' ')];
                    this.plotFunctions = [this.plotFunctions; metaClass.Name '.' metaMethod.Name];
                end

                %sort plots by class, name
                [~, order] = sort(this.plotNames);
                this.plotNames = this.plotNames(order);
                this.plotFunctions = this.plotFunctions(order);
            end
        end

        function graphSelectedPlots(this, ~, ~)
            sim = this.simulation;
            selectedPlots = sort(get(this.plotListBox, 'Value'));
            time = permute(sim.state('Time').values, [3 1 2]);
            compartments = sort(get(this.compartmentListBox, 'Value') - 4);

            rows = ceil(sqrt(length(selectedPlots)));
            cols = ceil(length(selectedPlots) / rows);
            for i = 1:length(selectedPlots)
                axes = subplot(rows, cols, i, 'Parent', this.plotPanel);
                cla(axes);
                legend(axes, 'off');
                set(axes, 'Box', 'on');

                feval(this.plotFunctions{selectedPlots(i)}, sim, axes, time, compartments);

                set(axes, 'fontSize', 6);
                title(axes, this.plotNames{selectedPlots(i)}, 'fontSize', 10);
            end
        end
    end

    %Callbacks
    methods (Access = protected)
        function layoutWindow(this, ~, ~)
            %figure size
            figureSize = get(this.figure, 'Position');
            figureWidth = figureSize(3);
            figureHeight = figureSize(4);

            %list box panel
            panelWidth = 300;
            set(this.buttonsPanel, 'Position', [0 0 panelWidth / figureWidth 1]);

            %select compartment
            compartmentListBoxWidth = panelWidth - 20;
            compartmentListBoxHeight = 100;
            compartmentListBoxSize = [10 120 compartmentListBoxWidth compartmentListBoxHeight];
            set(this.compartmentListBoxTitle, 'Position', [10 120+compartmentListBoxHeight compartmentListBoxWidth 20]);
            set(this.compartmentListBox, 'Position', compartmentListBoxSize);

            %plot list box
            plotListBoxWidth = compartmentListBoxWidth;
            plotListBoxHeight = figureHeight-140-compartmentListBoxHeight-35-10-20;
            plotListBoxSize = [10 120+compartmentListBoxHeight+10+35 plotListBoxWidth plotListBoxHeight];
            set(this.plotListBoxTitle, 'Position', [10 120+compartmentListBoxHeight+10+35+plotListBoxHeight compartmentListBoxWidth 20]);
            set(this.plotListBox, 'Position', plotListBoxSize);

            %plot button
            buttonWidth = panelWidth-20;
            buttonHeight = 20;
            buttonSize = [10 70 buttonWidth buttonHeight];
            set(this.plotButton, 'Position', buttonSize);

            %previous page button
            buttonWidth = (panelWidth-20-10)/2;
            buttonSize = [10 40 buttonWidth buttonHeight];
            set(this.prevButton, 'Position', buttonSize);

            %next page button
            buttonSize = [10+buttonWidth+10 40 buttonWidth buttonHeight];
            set(this.nextButton, 'Position', buttonSize);

            %save button
            buttonSize = [10 10 buttonWidth buttonHeight];
            set(this.saveButton, 'Position', buttonSize);

            %print button
            buttonSize = [10 + buttonWidth + 10 10 buttonWidth buttonHeight];
            set(this.printButton, 'Position', buttonSize);

            %plot panel
            set(this.plotPanel, 'Position', [panelWidth / figureWidth 0 1 - panelWidth/figureWidth 1]);
        end

        function graphPrevPlots(this, ~, ~)
            selectedPlots = get(this.plotListBox, 'Value');

            unselectedPlots = setdiff(1:length(this.plotNames), selectedPlots);

            %select next plots
            if ...
                    max(selectedPlots) - min(selectedPlots) + 1 == length(selectedPlots) || ...
                    max(unselectedPlots) - min(unselectedPlots) + 1 == length(unselectedPlots)
                selectedPlots = mod(selectedPlots - length(selectedPlots) - 1, length(this.plotNames)) + 1;
            else
                selectedPlots = (1:length(selectedPlots)) + min(selectedPlots) - 1;
            end
            set(this.plotListBox ,'Value', sort(selectedPlots));

            %graph
            this.graphSelectedPlots();
        end

        function graphNextPlots(this, ~, ~)
            selectedPlots = get(this.plotListBox, 'Value');

            unselectedPlots = setdiff(1:length(this.plotNames), selectedPlots);

            %select next plots
            if ...
                    max(selectedPlots) - min(selectedPlots) + 1 == length(selectedPlots) || ...
                    max(unselectedPlots) - min(unselectedPlots) + 1 == length(unselectedPlots)
                selectedPlots = mod(selectedPlots+length(selectedPlots)-1,length(this.plotNames))+1;
            else
                selectedPlots = (1:length(selectedPlots)) + min(selectedPlots) - 1;
            end
            set(this.plotListBox, 'Value', sort(selectedPlots));

            %graph
            this.graphSelectedPlots();
        end

        function save(this, ~, ~)
            filemenufcn(this.figure, 'FileSaveAs')
        end

        function print(this, ~, ~)
            printdlg(this.figure);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function r = tenRank( A, tol )
%TENRANK Calculates the multilinear rank.
%   R = TENRANK( A, TOL ) calculates the multilinear rank R of the
%   given tensor A. All singular values below the specified tolerance
%   TOL are considered to be zero.
%
%   R = TENRANK( A ) calculates the multilinear rank R of the given
%   tensor A. The default tolerance 1e-5 is chosen for the cutoff.
%

%   GeomCG Tensor Completion. Copyright 2013 by
%   Michael Steinlechner
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

    if ~exist('tol','var')
        tol = 1e-5;
    end

    r = zeros( 1, ndims( A ) );

    for i = 1:ndims( A )
        r(i) = rank( double( tenmat( A, i ) ), tol);
    end
end

author: andrewssobral