text:
function varargout=lpcar2rf(varargin)
%
% For calling details please see v_lpcar2rf.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcar2rf.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcar2rf.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcar2rf(varargin{:});
else
    v_lpcar2rf(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function w = fix_data_length(w, maxlen)
   %FIX_DATA_LENGTH adjust length of waveform data to allow batch processing
   %   waveform = fix_data_length(waveform)
   %       adjusts all waveforms to the length of the largest, while
   %       zero-padding all shorter waveforms
   %
   %   waveform = fix_data_length(waveform, maxlength)
   %       sets all waveform lengths to maxlength.  This use of the function
   %       has been superceeded by set(waveform,'samplelength',maxlength);
   %
   %  examples
   %       % let inWaves be a 1x2 waveform object
   %       % 3000 samples in inWaves(1)
   %       % 10025 samples in inWaves(2)
   %
   %       % set both waves' data to a length of to 10025 while padding the
   %       % smaller of the two with zeroes.
   %       outWaveforms = fix_data_length(inWaves)
   %
   %       % set both sample lengths to 500 truncating both of them...
   %       outWaveform = fix_data_length(inWaves, 500)
   %
   %       %The above example is nearly the same as
   %       outWaveform = set(inWaves,'sample_length',500)
   %
   %       Behaviorally, this differs from the set command because it
   %       automatically determines the maximum desired length, when not
   %       specified.
   %
   %
   % See also WAVEFORM/EXTRACT, WAVEFORM/DOUBLE, WAVEFORM/SET -- Sample_Length

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % Note: Glenn Thompson 2016 - this function could be simplified by using
   % waveform/get_time_range and waveform/pad
   Wcount = numel(w);

   if ~exist('maxlen','var')
      m(Wcount) = 0;
      for j = Wcount : -1 : 1
         m(j) = numel(w(j).data);
      end
      maxlen = max(m);
   end

   if length(maxlen) == 2
      st = maxlen(1);
      ed = maxlen(2);
   else
      st = 1;
      ed = maxlen(end);
   end

   for j = 1 : Wcount
      D = w(j).data;
      if ed > numel(D)
         D(ed) = 0;
      else
         D = D(st:ed);
      end
      w(j) = set(w(j),'data',D);
   end;
end

author: geoscience-community-codes

repo: GISMO


text:
function tf=moxunit_isa_test_skipped_exception(exception)
% return true if the input is a MOxUnit test-skipped exception
%
% tf=moxunit_isa_test_skipped_exception(exception)
%
% Input:
%   exception           an exception struct or object with field
%                       .identifier
%
% Output:
%   tf                  true if exception.identifer is equal to that
%                       thrown by moxunit_throw_test_skipped_exception
%
% Notes:
%   - if this exception is thrown when running a test through
%     a MOxUnitTestReport instance, this instance will record the test as
%     skipped, not as a failure
%
% NNO Jan 2014

    persistent cached_test_skipped_identifier

    if isempty(cached_test_skipped_identifier)
        try
            moxunit_throw_test_skipped_exception('error');
            assert(false,'should never get here');
        catch
            last_caught_error=lasterror();
        end

        cached_test_skipped_identifier=last_caught_error(1).identifier;
    end

    tf=strcmp(exception.identifier,...
                cached_test_skipped_identifier);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function f14=f14(x)

Bound=[-65.536 65.536];

if nargin==0
    f14 = Bound;
else
aij=[   -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;
    -32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];

x1=x(1,:)*ones(1,25);
x2=x(2,:)*ones(1,25);

j=cumsum(ones(1,25));

f14=1/(1/500+sum(1./(j+(x1-aij(1,:)).^6+(x2-aij(2,:)).^6)));
end

author: sfvsfv

repo: Mathematical-modeling


text:
function x = subsref( x, S, cheat )
error( 'Subscripts not allowed for this object.' );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [vol, sens, cfg] = prepare_headmodel(cfg, data)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that helps to prepare the electrodes/gradiometers and the volume
% this is used in sourceanalysis and dipolefitting
%
% This function will get the gradiometer/electrode definition from
%   cfg.channel
%   cfg.gradfile
%   cfg.elecfile
%   cfg.elec
%   cfg.grad
%   data.grad
%   data.elec
% and the volume conductor definition from
%   cfg.hdmfile
%   cfg.vol
%   data.vol
% Subsequently it will remove the gradiometers/electrodes that are not
% present in the data. Finally it with attach the gradiometers to a
% multi-sphere head model (if supplied) or attach the electrodes to
% a BEM head model.
%
% This function will return the electrodes/gradiometers in an order that is
% consistent with the order in cfg.channel, or in case that is empty in the
% order of the input electrode/gradiometer definition.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2004-2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: prepare_headmodel.m 952 2010-04-21 18:29:51Z roboos $

% set the defaults
if ~isfield(cfg, 'channel'), cfg.channel = 'all';   end
if ~isfield(cfg, 'order'),   cfg.order = 10;        end % order of expansion for Nolte method; 10 should be enough for real applications; in simulations it makes sense to go higher

if nargin<2
  data = [];
end

% get the volume conduction model
if isfield(cfg, 'hdmfile')
  fprintf('reading headmodel from file ''%s''\n', cfg.hdmfile);
  vol = ft_read_vol(cfg.hdmfile);
elseif isfield(cfg, 'vol')
  fprintf('using headmodel specified in the configuration\n');
  vol = cfg.vol;
elseif isfield(data, 'vol')
  fprintf('using headmodel specified in the data\n');
  vol = data.vol;
else
  error('no headmodel specified');
end

% get the gradiometer or electrode definition
if isfield(cfg, 'gradfile')
  fprintf('reading gradiometers from file ''%s''\n', cfg.gradfile);
  sens = read_sens(cfg.gradfile);
elseif isfield(cfg, 'grad')
  fprintf('using gradiometers specified in the configuration\n');
  sens = cfg.grad;
elseif isfield(data, 'grad')
  fprintf('using gradiometers specified in the data\n');
  sens = data.grad;
elseif isfield(cfg, 'elecfile')
  fprintf('reading electrodes from file ''%s''\n', cfg.elecfile);
  sens = read_sens(cfg.elecfile);
elseif isfield(cfg, 'elec')
  fprintf('using electrodes specified in the configuration\n');
  sens = cfg.elec;
elseif isfield(data, 'elec')
  fprintf('using electrodes specified in the data\n');
  sens = data.elec;
else
  error('no electrodes or gradiometers specified');
end

if isfield(data, 'topolabel')
  % the data reflects a componentanalysis, where the topographic and the
  % timecourse labels are different
  cfg.channel = ft_channelselection(cfg.channel, data.topolabel);
elseif isfield(data, 'label')
  % In the subsequent code, the matching channels in the sensor array and
  % in the configuration will be selected. To ensure that these channels
  % are also present in the data, update the configuration to match the data.
  cfg.channel = ft_channelselection(cfg.channel, data.label);
else
  % update the selected channels based on the electrode/gradiometer definition
  cfg.channel = ft_channelselection(cfg.channel, sens.label);
end

% ensure that these are a struct, which may be required in case configuration tracking is used
vol  = struct(vol);
sens = struct(sens);

% the prepare_vol_sens function from the forwinv module does most of the actual work
[vol, sens] = ft_prepare_vol_sens(vol, sens, 'channel', cfg.channel, 'order', cfg.order);

% update the selected channels in the configuration
cfg.channel = sens.label;

author: PatternRecognition

repo: OpenBMI


text:
function prediction = rmHrfSearchFit_oneGaussian(model, params, loopSlices, wProcess,varexp,allstimimages)
% rmHrfSearchFit_oneGaussian - make predictions without HRF
%
% prediction = rmHrfSearchFit_oneGaussian(model, params, loopSlices, wProcess)
%
% 10/2009: SD & WZ wrote it.

% intiate data stucture for our pRF values
tmp.x    = [];
tmp.y    = [];
tmp.s    = [];
tmp.ve   = [];

% Get the data required to build our pRF
for slice=loopSlices,
    % now we extract only the data from that slice
    s = rmSliceGet(model,slice);

    % if we have more than one slice, then check which voxels to process in
    % each slice
    if length(loopSlices) > 1
        wProcessSlice = wProcess{slice};
        varexpSlice  = varexp{slice};
    else
        wProcessSlice = wProcess;
        varexpSlice  = varexp;
    end

    % store
    tmp.x = [tmp.x s{1}.x0(wProcessSlice)];
    tmp.y = [tmp.y s{1}.y0(wProcessSlice)];
    tmp.s = [tmp.s s{1}.s(wProcessSlice)];
    tmp.ve = [tmp.ve varexpSlice(wProcessSlice)];
end

% build pRF and make our predictions
n = numel(tmp.x);
s = [[1:ceil(n./1000):n-2] n+1]; %#ok<NBRAK>
prediction = zeros(size(allstimimages,1),n);
fprintf(1,'[%s]:Making predictions without HRF for each voxel (%d):',mfilename,n);
drawnow;tic;
for n=1:numel(s)-1,
    % make rfs
    rf   = rfGaussian2d(params.analysis.X, params.analysis.Y,...
        tmp.s(s(n):s(n+1)-1), ...
        tmp.s(s(n):s(n+1)-1),...
        zeros(size(tmp.s(s(n):s(n+1)-1))), ...
        tmp.x(s(n):s(n+1)-1), ...
        tmp.y(s(n):s(n+1)-1));

    % convolve with stimulus
    pred = allstimimages*rf;

    % store
    prediction(:,s(n):s(n+1)-1) = pred;
    fprintf(1,'.');drawnow;
end;
clear n s rf pred;
fprintf(1, 'Done[%d min].\t(%s)\n', round(toc/60), datestr(now));
drawnow;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function yp = p06_fun ( neqn, t, y )

%*****************************************************************************80
%
%% P06_FUN evaluates the function for problem P06.
%
%  Discussion:
%
%    y1' = 2 y1 * ( 1 - y2 )
%    y2' = - y2 * ( 1 - y1 )
%    y1(0) = 1
%    y2(0) = 3
%
%    2 equations.
%    Enright and Pryce nonstiff problem #B1.
%    Autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wayne Enright, John Pryce,
%    Algorithm 648,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 1, pages 28-34.
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the derivative
%    function.
%
%    Output, real YP(NEQN), the value of the derivative function.
%
  yp = zeros ( neqn, 1 );

  yp(1) = 2.0 * y(1) * ( 1.0 - y(2) );
  yp(2) =     - y(2) * ( 1.0 - y(1) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Z_mean{i}: [num_components_for_subj_i x num_clusters]
%            contains posterior prob. of cluster membership
% cmp_keep:  [1 x num_subjects] list of preserved components for each
%           subject
% B_bar_mean{i,j}:  B-spline basis coefficients from cluster i to cluster j
%                   (note this is the transpose of my normal format)
% phi_t:        [num_b_splines x time] time-varying b-spline basis function
% Theta_diag:   [num_b_splines x num_PCs] FPCA weights for
%               self-connectivity
% Theta_offdiag: [num_b_splines x num_PCs] FPCA weights for
%               cross-connectivity
% NOTE: To obtain the 80x1 time-varying group PDC from cluster i to cluster
%       j -- F(i,j):
%       F(j,j) = phi_t'Theta_diag'B_bar_mean{j,j}
%       F(i,j) = phi_t'Theta_offdiag'B_bar_mean{i,j}  where j ~= i
%
% S_bar_mean:   [num_clusters x 3] MNI centroids for each cluster
%
%

%%
clear all; clc;

load real_data
load smoothing_outputs
load mcmc_outputs

%% Get posterior means spatial locations & connectivity coefs & indicators
S_bar_mean=zeros(M,3);
for i=1:size(S_BAR_array,1)
    S_bar_mean=S_bar_mean+S_BAR_array{i}/size(S_BAR_array,1);
end
B_bar_mean=cell(M,M);
for k1=1:M
    for k2=1:M
        B_bar_mean{k1,k2}=zeros(1,Q);
        for iter=1:size(B_BAR_array,1)
            B_bar_mean{k1,k2}=B_bar_mean{k1,k2}+B_BAR_array{iter}{k1,k2}/size(B_BAR_array,1);
        end
    end
end
Z_mean=cell(N,1);
for i=1:N
    Z_mean_i=Z_array{1}{i};
    for iter=2:size(Z_array,1)
        Z_mean_i=Z_mean_i+Z_array{iter}{i};
    end
    Z_mean_i=Z_mean_i/size(Z_array,1);
    Z_mean{i}=Z_mean_i;
end
%% Map posterior connectivities of group components
C_BAR_array=cell(size(B_BAR_array));
for iter=1:size(B_BAR_array,1)
    C_BAR_array{iter}=cell(M);
    for k1=1:M
        for k2=1:M
            if k1==k2
                C_BAR_array{iter}{k1,k2}=phi_t'*Theta_diag*B_BAR_array{iter}{k1,k2}';
            end
            if not(k1==k2)
                C_BAR_array{iter}{k1,k2}=phi_t'*Theta_offdiag*B_BAR_array{iter}{k1,k2}';
            end
        end
    end
end

labels = {'BA19/39 (+/-36,-70,23)' , 'BA31 (12,-41,-34)' , 'BA19 (-32,-70,22)' , 'BA32 (-16,13,34)' , 'BA32/ACC (1,30,29)' , 'BA24/MCC (1,-2,35)' , 'BA32 (20,7,32)' , 'BA6/preCG (52,-1,14)' , 'BA2/pCG (-33,-28,43)' , 'BA20 (-39,-9,-19)' , 'BA7/Cuneus (3,-73,32)' , 'Parietal WM (27,-35,37)'};

%% SELECT SIGNIFICANT COMPONENTS
num_subjects = length(subj_keep);
num_clusters = size(S_bar_mean,1);
T=80;
BASELINE = [-0.7 -0.5];
TAIL = 'both';
erWinCenterTimes = linspace(-0.75, 1.7188, T);

p_cmp = 0.5;  % probably threshold for cluster membership
p_sub = 0.5; % proportion of subjects that must be assigned (with probability > p_cmp) to a cluster in order to keep the cluster
p_M=zeros(num_clusters,num_subjects);  % p_M(k,i) probability that subject i has membership in cluster k
for k=1:num_clusters
    for i=1:num_subjects
        p_M(k,i)=max(Z_mean{i}(:,k));
    end
end
keep=zeros(num_clusters,1);
for k=1:num_clusters
    pr_k=sum(p_M(k,:)>p_cmp)/num_subjects;
    if pr_k>p_sub
        keep(k)=1;
    end
end
S_bar_mean=S_bar_mean(find(keep==1),:);
num_clusters = size(S_bar_mean,1);

%% PUT MNI COORDS IN DIPFIT STRUCTURE

%% insert symmetric dipole on opposite y-hemisphere
% indices of dual equivalent dipoles
dualEquivDipoles = 0;

for k=1:num_clusters
    centroids(k).momxyz = eps*ones(2,3);
    centroids(k).rv = 0;
    centroids(k).clustid = k;

    if ismember_bc(k,dualEquivDipoles)
        centroids(k).posxyz = [S_bar_mean(k,:); S_bar_mean(k,:).*[1 -1 1]];
        centroids(k).select = [1 2];
    else
        centroids(k).posxyz = [S_bar_mean(k,:); [0 0 0]];
        centroids(k).select = 1;
    end
end

%% do a dipplot of the cluster centroids

dipplot(centroids,'coordformat','Spherical','plot','on','projlines','on');

%% Create connectivity structure

M_keep=sum(keep);
ind_keep=find(keep==1);
alpha = 0.01;
if strcmpi(TAIL,'both')
    alpha=alpha/2;
end
% t_bl=10;    % time points 1:t_bl form baseline
q_l=alpha; % quantile for lower CI bound
q_u=1-alpha; % quantile for upper CI bound

baseidx = getindex(erWinCenterTimes,BASELINE);

C_bar_bl_mean=cell(size(C_BAR_array,1),1);
for iter=1:size(C_BAR_array,1)
    C_bar_bl_mean{iter}=zeros(M_keep);
    for k1=1:M_keep
        for k2=1:M_keep
            C_bar_bl_mean{iter}(k1,k2)=mean(C_BAR_array{iter}{ind_keep(k1),ind_keep(k2)}(baseidx(1):baseidx(2)));
        end
    end
end

C_bar_lower_keep=cell(M_keep);
C_bar_mean_keep=cell(M_keep);
C_bar_upper_keep=cell(M_keep);
C_bar_sgn_keep=cell(M_keep); % indicates which time points are significantly different from baseline mean
C_bar_keep=cell(M_keep);
for k1=1:M_keep
    for k2=1:M_keep
        C_bar_lower_keep{k1,k2}=zeros(T,1);
        C_bar_mean_keep{k1,k2}=zeros(T,1);
        C_bar_upper_keep{k1,k2}=zeros(T,1);
        C_bar_sgn_keep{k1,k2}=zeros(T,1);
        C_bar_keep{k1,k2}=zeros(T,1);
        for t=1:T
            tmp=zeros(size(C_BAR_array,1),1);
            for iter=1:size(C_BAR_array,1)
                tmp(iter)=C_BAR_array{iter}{ind_keep(k1),ind_keep(k2)}(t);
            end
            C_bar_lower_keep{k1,k2}(t)=prctile(tmp,100*q_l);
            C_bar_mean_keep{k1,k2}(t)=mean(tmp);
            C_bar_upper_keep{k1,k2}(t)=prctile(tmp,100*q_u);
            if strcmpi(TAIL,'both')
                if or((C_bar_lower_keep{k1,k2}(t)-C_bar_bl_mean{iter}(k1,k2))>0,(C_bar_upper_keep{k1,k2}(t)-C_bar_bl_mean{iter}(k1,k2))<0)
                    C_bar_sgn_keep{k1,k2}(t)=1; % significant
                end
            else
                if (C_bar_lower_keep{k1,k2}(t)-C_bar_bl_mean{iter}(k1,k2))>0
                    C_bar_sgn_keep{k1,k2}(t)=1; % significant
                end
            end
            C_bar_keep{k1,k2}(t)=exp(C_bar_mean_keep{k1,k2}(t))-1;
        end
    end
end

for i=1:M_keep
    for j=1:M_keep
        GroupConn.dDTF08(j,i,1,:) = C_bar_keep{i,j};
    end
end

GroupConn.erWinCenterTimes = erWinCenterTimes;
GroupConn.freqs = 5;  % 3-7 Hz

%% Create SIFT EEG structure
EEG_orig = pop_loadset('C:/aaa_wes/brains/multisubject_eeg/For_Wes/Data/eb79RespWrong_Connectivity.set');

SamplingRate = EEG_orig.srate;

EEG = eeg_emptyset;
EEG.data = zeros(EEG_orig.nbchan,T);
EEG.icaweights = eye(num_clusters,EEG_orig.nbchan);
EEG.icasphere =  EEG.icaweights';
EEG.icachansind = 1:num_clusters;
EEG.icaact = [];
EEG.srate = SamplingRate;
EEG.pnts = T;
EEG.trials = 1;
EEG.chanlocs = EEG_orig.chanlocs;
EEG.times = ((0:(EEG.pnts-1))/SamplingRate)*1000;   % ms
EEG.xmin = EEG.times(1);
EEG.xmax = EEG.times(end)/1000;  % sec
EEG.nbchan = EEG_orig.nbchan;
EEG.dipfit = EEG_orig.dipfit;
EEG.dipfit.model = centroids;
EEG.dipfit.chansel = 1:num_clusters;
EEG = eeg_checkset(EEG);
% pop_eegplot(EEG);

EEG.CAT.curComps = 1:num_clusters;
EEG.CAT.curComponentNames = labels(find(keep)); %strtrim(cellstr(num2str(EEG.CAT.curComps')))';
EEG.CAT.Conn = GroupConn;
EEG.CAT.MODEL = struct([]);
EEG.CAT.nbchan = num_clusters;
EEG.CAT.srcdata = EEG.data;

%% Make statistics
for i=1:size(EEG.CAT.Conn.dDTF08,1)
    for j=1:size(EEG.CAT.Conn.dDTF08,2)
        Stats.dDTF08.ci(1,j,i,1,:) = exp(C_bar_lower_keep{i,j})-1;  % lower bound
        Stats.dDTF08.ci(2,j,i,1,:) = exp(C_bar_upper_keep{i,j})-1;  % upper bound
        Stats.dDTF08.pval(j,i,1,:) = ~C_bar_sgn_keep{i,j}; %squeeze(Stats.dDTF08.ci(1,j,i,1,:))<=0 & squeeze(Stats.dDTF08.ci(2,j,i,1,:))>=0;  % insignificant if zero is contained within lower and upper bound
        Stats.alpha = 0.05;
    end
end

%%
pop_dipplot(EEG,1:num_clusters,'dipolelength',0,'spheres','on','projlines','on','coordformat','Spherical','num','on','dipnames',EEG.CAT.curComponentNames);

%% Visualize
[figureHandles tfgridcfg] = vis_TimeFreqGrid('ALLEEG',EEG,'Conn',EEG.CAT.Conn,'Stats',Stats,'MatrixLayout',{'Partial','UpperTriangle', 'dDTF08', 'LowerTriangle','dDTF08','Diagonal','none'},'TimesToPlot',[-0.7 1],'ColorLimits',[0 0.05],'Baseline',[],'Smooth2D',false,'Thresholding',{'Statistics','plotci',true,'ThresholdingMethod','pval','AlphaSignificance',0.05},'BackgroundColor',[0 0 0],'TextColor',[1 1 1],'LineColor',[0 0 0],'AxesFontSize',14);

%% threshold according to stats
load brainmovieconfig
EEGthresh = EEG;
EEGthresh.CAT.Conn.dDTF08(Stats.dDTF08.pval>=alpha) = 0;

cfg=pop_vis_causalBrainMovie3D(EEGthresh,cfg);

repo: REST


text:
% ***************************************************************************
% Parabolic Interpolation
% ***************************************************************************
% Author: Chaobin
% Email:  chaubyZou@163.com
% Date: October 2020
% ***************************************************************************
% Language: Matlab
% Also available in: Python
% Required library: None
% ***************************************************************************

classdef ParabolicInterpolation < handle
    properties
        name = 'parabolic interpolation';
        q_via = [];
        t_via = [];
    end

    methods
        % crate the objective
        % name: string
        % q_via: N x 3 array
        % t_via: N x 1 array
        function obj = ParabolicInterpolation(name, q_via, t_via)
            obj.name = name;
            obj.q_via = q_via;
            obj.t_via = t_via;
            if length(q_via(:,1)) ~= length(t_via)
                error('The q_via and t_via must have a same length');
            end
        end

        % parabolic interpolation with two data points
        % q0: the first data point
        % q1: the second data point
        % v0: the velocity of the first data point
        % v1: the velocity of the second data point
        % t0: the time of the first data point
        % t1: the time of the second data point
        % tf: the time of the flex point
        % qf: the position of the flex point
        % a0~a5: parameters
        function [a0, a1, a2, a3, a4, a5] = parabolic(obj, q0, q1, v0, v1, t0, t1, tf, qf)
            if abs(t0 - t1) < 1e-6
                error('t0 and t1 must be different');
            end

            if ((tf <= t0) || (tf >= t1))
                error('tf must satisfy t0 < tf < t1');
            end

            if ((qf <= min(q0, q1)) || (qf >= max(q0, q1)))
                error('qf must satisfy min(q0, q1) < qf < max(q0, q1)');
            end

            T = t1 - t0;
            h = q1 - q0;
            Ta = tf - t0;
            Td = t1 - tf;

            a0 = q0;
            a1 = v0;
            a2 = (2*h - v0*(T + Ta) - v1*Td)/(2*T*Ta);
            a3 = (2*q1*Ta + Td*(2*q0 + Ta*(v0 - v1)))/(2*T);
            a4 = (2*h - v0*Ta - v1*Td)/T;
            a5 = -(2*h - v0*Ta - v1*(T+Td))/(2*T*Td);
        end

        % parabolic interpolation for all data points
        % t: time
        % q: 1 x 3 array, output of the interpolation at time t
        function q = getPosition(obj, t)
            if (t < obj.t_via(1)) || (t > obj.t_via(end))
                error('The specific time error, time ranges error');
            end

            j = find(obj.t_via >= t, 1, 'first'); % find the index of t1
            if j == 1
                i = 1;
                j = 2;
            else
                i = j-1;
            end

            % get given position
            q0 = obj.q_via(i, 1);
            v0 = obj.q_via(i, 2);
            t0 = obj.t_via(i);

            q1 = obj.q_via(j, 1);
            v1 = obj.q_via(j, 2);
            t1 = obj.t_via(j);

            % symmetric acceleration
            tf = (t0 + t1)/2;
            qf = (q0 + q1)/2;

            % asymmetric acceleration, specify tf and qf by users
            % tf = ?
            % qf = ?

            [a0, a1, a2, a3, a4, a5] = obj.parabolic(q0, q1, v0, v1, t0, t1, tf, qf);

            if t <= tf
                q(1, 1) = a0 + a1*(t - t0) + a2*(t-t0)^2;
                q(1, 2) = a1 + 2*a2*(t - t0);
                q(1, 3) = 2*a2;
            else
                q(1, 1) = a3 + a4*(t - tf) + a5*(t-tf)^2;
                q(1, 2) = a4 + 2*a5*(t - tf);
                q(1, 3) = 2*a5;
            end
        end
    end
end

author: chauby

repo: PolynomialInterpolation


text:
function epo= proc_meanAcrossTime(epo, ival, clab)
%PROC_MEANACROSSTIME - Average signals across a specified time interval
%
%Synopsis:
%  EPO= proc_meanAcrossTime(EPO, <IVAL>, <CLAB>)
%
%Arguments:
%  DAT  - data structure of epoched data
%  IVAL - interval in which the average is to be calculated,
%         default [] which means the whole time range
%  CLAB - cell array of channels to be selected, default all
%
%Returns:
%  DAT  - updated data structure

misc_checkType(epo, 'STRUCT(x)');
misc_checkTypeIfExists('ival', 'DOUBLE[2]');
misc_checkTypeIfExists('clab', 'CHAR|CELL{CHAR}');

if nargin<2,
  ival= [];
end
if nargin<3,
  clab= {};
end

if isempty(ival),
  idx= 1:size(epo.x,1);
else
  idx= util_timeind(epo, ival);
end
epo.x= mean(epo.x(idx,:,:), 1);
if isfield(epo, 't') && length(idx)>1,
  epo.t= mean(epo.t(idx(1:end-1)));
end

if ~isempty(clab),
  epo= proc_selectChannels(epo, clab);
end

author: bbci


text:
function fe2d_r_fast_test ( )

%*****************************************************************************80
%
%% FE2D_R_FAST_TEST tests the FE2D_R_FAST code.
%
%  Discussion:
%
%    This function sets all parameter values and initial condition information
%    necessary to execute the "fast" version of the fe2d_r algorithm.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    28 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Reference:
%
%    Marcus R Garvie, John Burkardt, Jeff Morgan,
%    Simple Finite Element Methods for Approximating Predator-Prey Dynamics
%    in Two Dimensions using MATLAB,
%    Submitted to Bulletin of Mathematical Biology, 2014.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2D_R_FAST_TEST:\n' );
  fprintf ( 1, '  Test the FE2D_R_FAST function\n' );
  fprintf ( 1, '  which applies Robin boundary conditions as it\n' );
  fprintf ( 1, '  approximates a solution to a predator-prey system.\n' );
%
%  Set the parameters.
%
  alpha = 0.4;
  beta = 2.0;
  gamma = 0.6;
  delta = 1.0;
%
%  Use T=150.0 for normal run.
%  Use T=0.50 for a "quick" run that might take 15 minutes of computing.
%
% T = 150.0;
  T = 0.50;
  delt = 1.0 / 384.0;
  k1 = 0.01;
  k2 = 0.01;

  t = tic;
  fe2d_r_fast ( alpha, beta, gamma, delta, T, delt, @u0f, @v0f, k1, k2 );
  t = toc ( t );

  fprintf ( 1, '  Execution took %10.2g minutes \n', t / 60.0 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FE2D_R_FAST_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

function value = u0f ( x, y )

%*****************************************************************************80
%
%% U0F evaluates the initial condition for U.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for U at (X,Y).
%
  value = 6.0 / 35.0 - 2.0E-07 * ( x - 0.1 * y - 225.0 ) * ( x - 0.1 * y - 675.0 );

  return
end

function value = v0f ( x, y )

%*****************************************************************************80
%
%% V0F evaluates the initial condition for V.
%
%  Licensing:
%
%    Copyright (C) 2014 Marcus R. Garvie.
%    See 'mycopyright.txt' for details.
%
%  Modified:
%
%    26 April 2014
%
%  Author:
%
%    Marcus R. Garvie.
%
%  Parameters:
%
%    Input, real X, Y, a location in the region.
%
%    Output, real VALUE, the initial condition for V at (X,Y).
%
  value = 116.0 / 245.0 - 3.0E-05 * ( x - 450.0 ) - 1.2E-04 * ( y - 150.0 );

  return
end


text:
function [ grid_num, point_num ] = gl_levels_minmax_size ( dim_num, ...
  level_min, level_max )

%*****************************************************************************80
%
%% GL_LEVELS_MINMAX_SIZE returns sizes for GL_LEVELS_MINMAX.
%
%  Discussion:
%
%    This routine can be used to determine the necessary size to be
%    allocated to arrays GRID_ORDER and GRID_POINT in a call to
%    GL_GRIDS_MINMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MIN, LEVEL_MAX, the minimum and maximum values of
%    LEVEL, the sum of the levels in each spatial coordinate.
%
%    Output, integer GRID_NUM, the number of Gauss Legendre
%    grids whose LEVEL value is between LEVEL_MIN and LEVEL_MAX.
%
%    Output, integer POINT_NUM, the total number of points in the grids.
%

%
%  Determine the total number of points that will be generated
%  by "going through the motions".
%
  point_num = 0;
  grid_num = 0;

  for level = level_min : level_max

    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      order_1d = gl_level_to_order ( dim_num, level_1d );

      order_nd = prod ( order_1d(1:dim_num) );

      point_num = point_num + order_nd;

      grid_num = grid_num + 1;

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function tests = test_haar
  tests = functiontests(localfunctions);
end

function test_qmf(testCase)
    f = spx.wavelet.haar.quad_mirror_filter();
    verifyTrue(testCase, spx.norm.is_unit_norm_vec(f));
end

author: indigits

repo: sparse-plex


text:
%% housekeeping
clear
clc
close all

%% create data

tmp=importdata('ht2014data.txt');

vnames={'real_pce','cpixfe','ffr','m2','fsi'};

rawdb=ts('1988M12',tmp.data(:,2:end),vnames);

rawdb=pages2struct(rawdb);

%% plot the data

figure('name','raw data');
for ii=1:numel(vnames)
    v=vnames{ii};
    subplot(3,2,ii)
    plot(rawdb.(v))
    title(v)
end

%% transformed data
%
% we use levels of the federal funds rate and the stress index and growth
% rates of real personal consumption expenditures (PCE), money and prices.
db=struct();
db.R=rawdb.ffr;
db.S=rawdb.fsi;
db.C=100*log(rawdb.real_pce/rawdb.real_pce{-1});
db.M=100*log(rawdb.m2/rawdb.m2{-1});
db.P=100*log(rawdb.cpixfe/rawdb.cpixfe{-1});

descr=struct('C','Consumption growth',...
    'P','Inflation',...
    'R','Feds Funds rate',...
    'M','Money growth',...
    'S','Financial stress index');

varlist=fieldnames(descr);

figure('name','Transformed data');
for ii=1:numel(varlist)
    v=varlist{ii};
    subplot(3,2,ii)
    plot(db.(v))
    title(descr.(v))
    axis tight
end
%% Model 1: SVAR models with cholesky identification
clc

models={'1v1c','2v1c','3v1c','1v2c','2v2c','3v2c',...
    '3vS2c','3vSC2c','3vSCP2c','3vSRM2c','3vRM2c','3vRMC2c'};

nlags=2;

exog={};

constant=true;

panel=[];

priors=struct();
priors.var=svar.prior_template();
priors.var.type='sz';

date_range={db.C.start,db.C.finish};

sv0=svar.empty(1,0);

tic

for ii=1:numel(models)

    fprintf(1,' -------------------- %s -------------------- \n',models{ii});

    [markov_chains,switch_prior,restrictions]=build_model(models{ii},varlist);

    sv0(1,ii)=svar(varlist,exog,nlags,constant,panel,markov_chains);

    priors.nonvar=switch_prior;

    sv0(1,ii)=estimate(sv0(1,ii),db,date_range,priors,restrictions);

end

toc

% 2601.732325 seconds

%% estimates

for ii=1:numel(models)
    clc, close all

    fprintf(1,' -------------------- %s -------------------- \n',models{ii});

    pmode=posterior_mode(sv0(ii))

    pause
    % Printing estimates

    print_structural_form(sv0(ii))

    pause
    % Printing solution

    print_solution(sv0(ii))
    pause

    % historical probabilities
    plot_probabilities(sv0(ii))
    pause

    % plots probabilities against data
    close all
    plot_data_against_probabilities(sv0(ii),'regime')
    pause

end

%% A model with endogenous switching

clc

[mctvp,switch_prior,restrictions]=build_model('2v1c',varlist);

mctvp.endogenous_probabilities={
    'syncvol_tp_1_2=1/(1+exp(a12-b12*S))'
    'syncvol_tp_2_1=a21'
    };

mctvp.probability_parameters={'a12','b12','a21'};

mctvp.controlled_parameters={'s(3)'};%'s(1)','s(2)','s(4)',

switch_prior=rmfield(switch_prior,{'syncvol_tp_1_2','syncvol_tp_2_1'});

switch_prior.a12={0.1,0,5,'normal'};

switch_prior.b12={0.1,2,5,'gamma'};

switch_prior.a21={0.2,0.2,0.2,'beta'};

svtvp=svar(varlist,exog,nlags,constant,panel,mctvp);

priors.nonvar=switch_prior;

svtvp=estimate(svtvp,db,date_range,priors,restrictions);

author: jmaih


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

function CstMeshInitiator(mws)

%' switch on FD-TET setting for accurate farfields

FDSolver = invoke(mws,'FDSolver');
invoke(FDSolver,'ExtrudeOpenBC','True');

Mesh = invoke(mws,'Mesh');
invoke(Mesh,'FPBAAvoidNonRegUnite','True');
invoke(Mesh,'ConsiderSpaceForLowerMeshLimit','False');
invoke(Mesh,'MinimumStepNumber','5');

MeshSettings = invoke(mws,'MeshSettings');
invoke(MeshSettings, 'SetMeshType','Hex');
invoke(MeshSettings, 'Set','RatioLimitGeometry','20');
MeshSettings = invoke(mws,'MeshSettings');
invoke(MeshSettings, 'SetMeshType','HexTLM');
invoke(MeshSettings, 'Set','RatioLimitGeometry','20');

PostProcess1D = invoke(mws,'PostProcess1D');
invoke(PostProcess1D,'ActivateOperation','vswr','true');
invoke(PostProcess1D,'ActivateOperation','yz-matrices','true');
end

repo: CST-MATLAB-API


text:
function kegg2mol(cid,molfileDir,mets,takeMajorMS,pH,takeMajorTaut)
% Retreive molfiles from the KEGG Compound database.
%
% kegg2mol(cid,molfileDir,mets,takeMajorMS,pH,takeMajorTaut)
%
% INPUTS
% cid               m x 1 array of KEGG Compound identifiers.
%
% OPTIONAL INPUTS
% molfileDir        Path to directory where molfiles will be saved. Default
%                   is molfilesFromKegg.
% mets              m x 1 array of metabolite identifiers. Used to name
%                   molfiles. KEGG Compound ID will be used if mets is
%                   omitted.
% takeMajorMS       {(0),1}. If 0 (default), molfiles will be downloaded
%                   directly from KEGG. If 1, they will save as the major
%                   microspecies at the pH specified in input variable pH.
%                   Computing major microspecies requires ChemAxon's
%                   calculator plugins (cxcalc).
% pH                pH for computing major microspecies.
% takeMajorTaut     {0,(1)}. If 1 (default), molfiles will be saved as the
%                   major tautomers of the major microspecies at the
%                   specified pH.
%
% WRITTEN OUTPUTS
% One molfile for each metabolite with a KEGG Compound ID.

% If takeMajorMS == 1; approximately 2 seconds per kegg id
% If takeMajorMS == 0; approximately 0.3 seconds per kegg id

% Configure inputs
cid = reshape(cid,length(cid),1);
if ischar(cid)
    cid = strtrim(cellstr(cid));
end
if isnumeric(cid)
    eval(['cid = {' regexprep(sprintf('''C%05d''; ',cid),'(;\s)$','') '};']);
end

if ~exist('molfileDir','var')
    molfileDir = 'molfilesFromKegg';
elseif isempty(molfileDir)
    molfileDir = 'molfilesFromKegg';
end
[success,message] = mkdir(molfileDir);
if success ~= 1
    error(message)
end
if ~strcmp(molfileDir(end),filesep)
    molfileDir = [molfileDir filesep];
end

if ~exist('mets','var')
    mets = cid;
elseif isempty(mets)
    mets = cid;
end
mets = reshape(mets,length(mets),1);
if ischar(mets)
    mets = strtrim(cellstr(mets));
end
if isnumeric(mets)
    mets = strtrim(cellstr(num2str(mets)));
end
mets = regexprep(mets,'(\[\w\])$',''); % Remove compartment assignment

if ~exist('takeMajorMS','var')
    takeMajorMS = 0;
elseif isempty(takeMajorMS)
    takeMajorMS = 0;
end
if ~exist('pH','var')
    pH = 7;
elseif isempty(pH)
    pH = 7;
end
if ~exist('takeMajorMS','var')
    takeMajorTaut = 1;
elseif isempty(takeMajorMS)
    takeMajorTaut = 1;
end

% Only retreive molfiles for unique metabolites
bool = ~cellfun('isempty',cid);
umets = mets(bool);
ucid = cid(bool);
[umets,crossi] = unique(umets);
ucid = ucid(crossi);

% List all elements
elements = {'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca',...
    'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y',...
    'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce',...
    'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir',...
    'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm',...
    'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Uut', 'Fl', 'Uup', 'Lv', 'Uus', 'Uuo'};

% Retreive molfiles
nomol = {};
for i = 1:length(umets)
    met = umets{i};
    id = ucid{i};
    [mol,success] = urlread(sprintf('http://rest.kegg.jp/get/cpd:%s/mol', id)); % API

    if success == 0
        nomol = [nomol; {met}];
        continue;
    end

    mol = regexprep(mol,'\r',''); % Remove carriage returns
    mol = regexprep(mol,'[^\n]*\n',[id '\n'],'once'); % Replace top line with KEGG ID
    mol = regexprep(mol,'M  END.*','M  END'); % Remove all lines after end of molfile

    % Identify variable structures containing R groups, repeat units etc.
    % ChemAxon's calculator plugins cannot compute major microspecies for
    % such structures.
    atoms = {};
    if takeMajorMS == 1
        % Get atom list
        lines = regexp(mol,'\n','split');
        atomCount = str2double(lines{4}(1:3));
        if atomCount > 0
            atomBlock = lines(5:5 + atomCount - 1);
            pat = '[^a-z_A-z]+(?<atom>[a-z_A-Z]+)[^a-z_A-z]+';
            atoms = regexp(atomBlock,pat,'names');
            atoms = [atoms{:}]';
            atoms = {atoms.atom}';
        end

        % Check for repeat units
        if ~isempty(strmatch('M  STY',lines))
            atoms = {};
        end
    end

    if takeMajorMS == 0 || ~all(ismember(atoms,elements)) || isempty(atoms) || (length(atoms) == 2 && all(strcmp(atoms,'H'))) % Write raw molfile from KEGG directly to file
        fid = fopen([molfileDir met '.mol'],'w+');
        fprintf(fid,'%s',mol);
        fclose(fid);
    else % Get major microspecies and write to file
        fid = fopen('tmp.mol','w+'); % Write raw molfile from KEGG to temporary file
        fprintf(fid,'%s',mol);
        fclose(fid);

        if takeMajorTaut == 1
            majorTautOption = 'true';
        else
            majorTautOption = 'false';
        end

        status = system(['cxcalc -o ' molfileDir met '.mol majorms -H ' num2str(pH) ' -f mol -M ' majorTautOption ' tmp.mol']); % Call ChemAxon's calculator plugin (cxcalc) to compute major microspecies

        if status ~= 0
            nomol = [nomol; {met}];
        end
    end
end

delete('tmp.mol'); % Delete temporary file

if ~isempty(nomol)
    fprintf(['Could not retreive molfiles from KEGG for metabolites:\n' sprintf('%s\n', nomol{:})]);
end

author: opencobra

repo: cobratoolbox


text:
function x = reciprocal_cdf_inv ( cdf, a, b )

%*****************************************************************************80
%
%% RECIPROCAL_CDF_INV inverts the Reciprocal CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A <= B.
%
%    Output, real X, the corresponding argument of the CDF.
%
  if ( cdf <= 0.0 )
    x = 0.0;
  elseif ( 0.0 < cdf )
    x = b^cdf / a^( cdf - 1.0 );
  end

  return
end


text:
function [ vTransVectorMean, vTransVectorCov ] = UnscentedTransform( hTransformFunction, vVectorMean, mVectorCov )
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

vectorOrder     = size(vVectorMean, 1);
numSigmaPoints  = (2 * vectorOrder) + 1;

mVectorCovSqrt = chol(mVectorCov);
if(vectorOrder <= 3) %<! Used as k Factor in literature
    unscentedTransformScalingFactor = 3 - vectorOrder;
else
    unscentedTransformScalingFactor = 3;
end

vSigmaPoints                = zeros(vectorOrder, numSigmaPoints);
vTransformedSigmaPoints     = zeros(vectorOrder, numSigmaPoints);
vUnscentedTransformWeights  = zeros(numSigmaPoints, 1);

vTransVectorMean    = zeros(vectorOrder, 1);
vTransVectorCov     = zeros(vectorOrder, vectorOrder);

mVectorCovSqrtNormalized = sqrt(vectorOrder + unscentedTransformScalingFactor) .* mVectorCovSqrt;

vSigmaPoints(1) = 0;
vUnscentedTransformWeights(1) = unscentedTransformScalingFactor / (unscentedTransformScalingFactor + vectorOrder);
vUnscentedTransformWeights(2:end) = 1 / (2 * (unscentedTransformScalingFactor + vectorOrder));

for ii = 2:(vectorOrder + 1)
    vSigmaPoints(:, ii) = mVectorCovSqrtNormalized(:, (ii - 1));
end

vSigmaPoints(:, (vectorOrder + 2):end) = -vSigmaPoints(:, 2:(vectorOrder + 1));

for ii = 1:numSigmaPoints
    vSigmaPoints(:, ii) = vVectorMean + vSigmaPoints(:, ii);
    vTransformedSigmaPoints(:, ii) = hTransformFunction(vSigmaPoints(:, ii));
end

for ii = 1:numSigmaPoints
    vTransVectorMean = vTransVectorMean + (vUnscentedTransformWeights(ii) * vTransformedSigmaPoints(:, ii));
end

for ii = 1:numSigmaPoints
    vTransVectorCov = vTransVectorCov + ...
        (vUnscentedTransformWeights(ii) * (([vTransformedSigmaPoints(:, ii) - vTransVectorMean]) * ([vTransformedSigmaPoints(:, ii) - vTransVectorMean].')));
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function rgb_display ( rgb, header )

%*****************************************************************************80
%
%% RGB_DISPLAY displays an RGB dataset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Nick Berry,
%    A "Practical" Use for Genetic Programming,
%    http://www.datagenetics.com/blog.html
%
%  Parameters:
%
%    Input, uint8 RGB(256,256,3), the RGB information.
%
%    Input, string HEADER, an identifier used for the title and filename.
%
  imshow ( rgb )
  title ( header, 'FontSize', 16 )

  filename = sprintf ( '%s.png', header );
  print ( '-dpng', filename );

  return
end


text:
%%
 %  Copyright (c) 2014, Facebook, Inc.
 %  All rights reserved.
 %
 %  This source code is licensed under the BSD-style license found in the
 %  LICENSE file in the root directory of this source tree. An additional grant
 %  of patent rights can be found in the PATENTS file in the same directory.
 %
 %%
function [ap,level2_scores] = test_level2_cnn(level2_model,test_labels,test_features,test_gfeatures, config)
%test level2 cnn

is_score = (size(test_features{1}.features,2) == 1);
level2_features_test = generate_level2_features_cnn(test_features,test_gfeatures,config);
clear test_features;
for attr_id = 1:size(test_labels,2)
    %%%test%%%
    %ignore those whith labels = 0
    idx  = find(test_labels(:,attr_id) ~= 0);

    if is_score
        level2_scores{attr_id} = level2_model{attr_id}.svm_weights'*level2_features_test';
    else
        level2_scores{attr_id} = level2_model{attr_id}.svm_weights(1:end-1)*level2_features_test'+...
            level2_model{attr_id}.svm_weights(end);
    end
    ap(attr_id)= get_precision_recall(level2_scores{attr_id}(idx),test_labels(idx,attr_id));
end

for attr_id=1:length(ap)
   fprintf('%4.2f & ',ap(attr_id)*100);
end
fprintf(' mean AP=%4.3f\n',mean(ap)*100);

end

author: facebookarchive

repo: pose-aligned-deep-networks


text:
classdef TotalCoordinatesCalculatorFactory < handle
    methods (Access = public, Static)
        function obj = create(cParams)
            switch cParams.nodes.vert
                case 4
                    obj = IntersectionCoordComputer(cParams);
                case 6
                    obj = DiagonalCoordComputer(cParams);
            end
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function L = probitLikelihood(noise, mu, varsigma, y)

% PROBITLIKELIHOOD Likelihood of data under probit noise model.

% NOISE

D = size(y, 2);
for i = 1:D
  mu(:, i) = mu(:, i) + noise.bias(i);
end
L = cumGaussian((y.*mu)./(sqrt(noise.sigma2+varsigma)));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function inside = p07_inside ( m, n, point )

%*****************************************************************************80
%
%% P07_INSIDE reports if a point is inside the region in problem 07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, logical INSIDE(N), is TRUE if the point is in the region.
%
  inside(1:n) = 1;

  [ lo, hi ] = p07_box ( m );
%
%  Check whether points are in the bounding box.
%
  for j = 1 : n

    for i = 1 : m
      if ( point(i,j) < lo(i) | hi(i) < point(i,j) )
        inside(j) = 0;
        continue
      end

    end

  end
%
%  Check whether points in the bounding box are in the region.
%
  for j = 1 : n

    if ( ~inside(j) )
      continue
    end

    if ( cos ( point(1,j) ) < point(2,j) )
      inside(j) = 0;
      continue
    end

    if ( point(2,j) < -5.0 + 5.0 * point(1,j)^4 / ( 2.5 * pi )^4 )
      inside(j) = 0;
      continue
    end

  end

  return
end


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Torr Vision Group (TVG)
%  University of Oxford - UK
%
%  Qizhu Li <liqizhu@robots.ox.ac.uk>
%  August 2018
% ------------------------------------------------------------------------
% This file is part of the weakly-supervised training method presented in:
%    Qizhu Li*, Anurag Arnab*, Philip H.S. Torr,
%    "Weakly- and Semi-Supervised Panoptic Segmentation,"
%    European Conference on Computer Vision (ECCV) 2018.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
% This script demos generation of iterative ground truths for weakly-
% supervised experiments.
%
% It merges CAM predictions with MCG&Grabcut masks to produce GT for
% the first round of iterative training.
%  INPUT:
%  - results/cam/<split>/*.png: the CAMs obtained from a multi-class
%  classifier.
%  - results/mcg_and_grabcut/<split>/*.png: the combined cues from MCG
%  and Grabcut.
% ------------------------------------------------------------------------

clearvars;
addpath scripts
addpath utils
addpath visualisation

% Extract detections from Cityscapes ground truth file and save as .mat
% this only need to be done once, as they do not change over iterative
% training stages
demo_instanceTrainId_to_dets;

clearvars;
dataset = 'cityscapes';
split = 'train';
opts = get_opts(dataset, split);
opts.list_path = 'lists/demo_id.txt';
opts.pred_dir = fullfile('cam', split);
opts.sem_save_dir = 'pred_sem_cam_mandg_merged';
opts.ins_save_dir = 'pred_ins_cam_mandg_merged';
opts.run_score_thresh = false;
opts.visualise_results = true;
[opts, results] = run_sub(opts);

% visualise
visualise_results_cam_mandg(opts, results);

author: qizhuli

repo: Weakly-Supervised-Panoptic-Segmentation


text:
function pf = minus(pf1,pf2)
% implements pf1 - pf2
%
% overload the - operator, i.e. one can now write pf1 - pf2 in order to
% subtract two pole figures from each other
%
% See also
% PoleFigure/PoleFigure PoleFigure/plus PoleFigure/mtimes

pf = pf1 + (-pf2);

author: mtex-toolbox

repo: mtex

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


text:
% "Powel test function which has n decision variables, has no local optima.
%has one global optimum. all decision variables are bounded in [-4
% 5]. One global minimum at x=[3 -1 0, 1...,3 -1 0 1]  with f(x)=0

function sum=test10(x)
sum=0;

b=numel(x)/4;

for i=1:b
   sum=sum+(x(4*i-3)+10*x(4*i-2)).^2+5*(x(4*i-1)-x(4*i)).^2+(x(4*i-2)-x(4*i-1)).^4+10*(x(4*i-3)-x(4*i)).^4;
end

author: Sable

repo: mcbench-benchmarks


text:
function mrSESSION = mrvMakeEmptySession;
%
% mrSESSION = mrvMakeEmptySession: [no arguments]
%
% Create a mrSESSION struct without
% any functionals assigned, for pseudo-realtime.
%
%
% ras 04/05.

[p f] = fileparts(pwd);

mrSESSION.mrLoadRetVersion = '3.10'; % I figure it's updated at to least that :)
mrSESSION.sessionCode = f;
mrSESSION.description = 'mrvTurbo';
mrSESSION.subject = '';
mrSESSION.examNum = '';
mrSESSION.inplanes = [];
mrSESSION.functionals = [];

dataTYPES.name = 'Original';

vANATOMYPATH = '';

save mrSESSION mrSESSION dataTYPES

return


text:
function pass = test_chebpolyval(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% NOTE: Since CHEBPOLYVAL() is basically a wrapper to CHEBTECH/FEVAL(), this is
% just a simple test.

% Init:
seedRNG(42)
tol = 10*pref.chebfuneps;

% Use CHEBPOLYVAL():
n = 10;
c = rand(10, 2);
x = rand(3);
fx1 = chebpolyval(c, x);

% Use CHEBPOLY():
T = chebpoly(0:9);
Tc = T*flipud(c);
fx2 = feval(Tc, x);

% Error:
err = norm(fx1 - fx2, inf);
pass(1) = err < n*tol;

%%

x = chebfun('x');
f = chebpolyval(c, x);
pass(2) = norm(f - Tc, inf) < n*tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function stroud_test10 ( )

%*****************************************************************************80
%
%% TEST10 tests CIRCLE_CUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_2D_INDEX;

  r = 3.0;
  xc = 0.0;
  yc = 0.0;

  num = function_2d_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10 \n' );
  fprintf ( 1, '  CIRCLE_CUM approximates an integral over a circle. \n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use radius R = %f\n', r );
  fprintf ( 1, '  and center: \n' );
  fprintf ( 1, '  XC = %f\n', xc );
  fprintf ( 1, '  YC = %f\n', yc );
  fprintf ( 1, '\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '    Order:      2             4              8            16 \n' );
  fprintf ( 1, '  F(X) \n' );
  fprintf ( 1, '\n' );

  for i = 1 : num

    FUNC_2D_INDEX = i;

    for j = 1 : 4

      norder = 2^j;

      result(j) = circle_cum ( 'function_2d', xc, yc, r, norder );

    end

    fname = function_2d_name ( i );

    fprintf ( 1, '  %s  %12f  %12f  %12f  %12f\n', fname, result(1:4) );

  end

  return
end


text:
function gridlines_test ( )

%*****************************************************************************80
%
%% GRIDLINES_TEST tests the GRIDLINES library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GRIDLINES_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the GRIDLINES library.\n' );

  gridlines_test01 ( );
  gridlines_test02 ( );
  gridlines_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GRIDLINES_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function A = unfold(T, rowDims, colDims)
%UNFOLD   Unfold (flatten or matricize) a discrete tensor.
%   A = UNFOLD(T, ROWDIMS, COLDIMS) returns an unfolding of T, where
%   ROWDIMS and COLDIMS are vectors of integers denoting the dimensions of
%   T that are merged into the rows and columns of the matrix A,
%   respectively.
%
%   A = UNFOLD(T, ROWDIMS) chooses COLDIMS automatically to contain all
%   dimensions not contained in ROWDIMS. The entries of COLDIMS are
%   arranged in the ascending order.
%
%   Note that UNFOLD(T, 1:ndims(T)) is the same as T(:).
%
%   Example:
%   T = randn(2, 3, 4, 5);
%   A = unfold(T, [2 4]); % Returns a matrix of size 15 x 8.
%
% See also CHEBFUN3/FOLD.

% The structure of this code is similar to `matricize.m` from the HTUCKER
% toolbox of Tobler and Kressner.

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

if ( nargin == 3 )
    rowDims = rowDims(rowDims <= ndims(T));
    colDims = colDims(colDims <= ndims(T));
elseif ( nargin == 2 )
    % Put remaining modes onto columns.
    colDims = setdiff(1:ndims(T), rowDims);
else
    error('CHEBFUN:CHEBFUN3:unfold', ...
        'COL_DIMS must be a positive integer vector without duplicate entries.')
end

% Size of tensor T
sizeT = size(T);
if ( ndims(T) < 3 )
    % To work properly also with scalars, vectors and matrices.
    sizeT(3) = 1;
end

% Permute dimensions of T
T = permute(T, [rowDims, colDims]);

% Flatten the tensor T
A = reshape(T, prod(sizeT(rowDims)), prod(sizeT(colDims)) );

end


text:
function colors = solver(A,targetIndex)
   colors = [];
end


text:
function Ex_GFK()
% This shows how to use GFK in a 1-nearest neighbor classifier.

% ref: Geodesic Flow Kernel for Unsupervised Domain Adaptation.
% B. Gong, Y. Shi, F. Sha, and K. Grauman.
% Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Providence, RI, June 2012.

% Contact: Boqing Gong (boqinggo@usc.edu)

%-------------------------I. setup source/target domains----------------------
% Four domains: { Caltech10, amazon, webcam, dslr }
src = 'webcam';
tgt = 'Caltech10';

d = 20; % subspace dimension, the following dims are used in the paper:
% webcam-dslr: 10
% dslr-amazon: 20
% webcam-amazon: 10
% caltech-webcam: 20
% caltech-dslr: 10
% caltech-amazon: 20
% Note the dim from X to Y is the same as that from Y to X.

nPerClass = 20;
% 20 per class when Caltech/Amazon/Webcam is the source domain, and
% 8 when DSLR is the source domain.

%--------------------II. prepare data--------------------------------------
load(['data/' src '_SURF_L10.mat']);     % source domain
fts = fts ./ repmat(sum(fts,2),1,size(fts,2));
Xs = zscore(fts,1);    clear fts
Ys = labels;           clear labels
Ps = princomp(Xs);  % source subspace

load(['data/' tgt '_SURF_L10.mat']);     % target domain
fts = fts ./ repmat(sum(fts,2),1,size(fts,2));
Xt = zscore(fts,1);     clear fts
Yt = labels;            clear labels
% Pt = princomp(Xt);  % target subspace
Pt = pca(Xt);  % target subspace
% Pt = ftProc_pca_tr(Xt);  % target subspace

fprintf('\nsource (%s) --> target (%s):\n', src, tgt);
fprintf('round     accuracy\n');
%--------------------III. run experiments----------------------------------
round = 20; % 20 random trials
tot = 0;
for iter = 1 : round
    fprintf('%4d', iter);

    inds = split(Ys, nPerClass);
    Xr = Xs(inds,:);
    Yr = Ys(inds);

    %---------------III.A. PLS --------------------------------------------
    % Ps = PLS(Xr, OneOfKEncoding(Yr), 3*d);
    % PLS generally leads to better performance.
    % A nice implementation is publicaly available at http://www.utd.edu/~herve/

    G = GFK([Ps,null(Ps')], Pt(:,1:d));
    [~, accy] = my_kernel_knn(G, Xr, Yr, Xt, Yt);
    fprintf('\t\t%2.2f%%\n', accy*100);
    tot = tot + accy;
end
fprintf('mean accuracy: %2.2f%%\n\n', tot/round*100);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [prediction accuracy] = my_kernel_knn(M, Xr, Yr, Xt, Yt)
dist = repmat(diag(Xr*M*Xr'),1,length(Yt)) ...
    + repmat(diag(Xt*M*Xt')',length(Yr),1)...
    - 2*Xr*M*Xt';
[~, minIDX] = min(dist);
prediction = Yr(minIDX);
accuracy = sum( prediction==Yt ) / length(Yt);

function [idx1 idx2] = split(Y,nPerClass, ratio)
% [idx1 idx2] = split(X,Y,nPerClass)
idx1 = [];  idx2 = [];
for C = 1 : max(Y)
    idx = find(Y == C);
    rn = randperm(length(idx));
    if exist('ratio')
        nPerClass = floor(length(idx)*ratio);
    end
    idx1 = [idx1; idx( rn(1:min(nPerClass,length(idx))) ) ];
    idx2 = [idx2; idx( rn(min(nPerClass,length(idx))+1:end) ) ];
end

author: viggin

repo: domain-adaptation-toolbox


text:
function ker = get_kernel(ker_name, cu_filename)

global data;

switch data.gpu.float_t
case 'double',
    ker_name = [ker_name '_d'];
case 'single',
    ker_name = [ker_name '_f'];
end

if ~isfield(data.gpu.kernels, ker_name)
    cu_filename = fullfile('cuda', cu_filename);
    ptx_filename = regexprep(cu_filename, '\.cu$', '.ptx');
    ker = struct();
    ker.name = ker_name;
    ker.ker = parallel.gpu.CUDAKernel(ptx_filename, cu_filename, ker_name);
    ker.block_size = ker.ker.ThreadBlockSize;
    ker.grid_size = ker.ker.GridSize;
    data.gpu.kernels.(ker_name) = ker;
end

ker = data.gpu.kernels.(ker_name);

author: jimmy-ren


text:
function noise = ncnmNoiseExpandParam(noise, params)

% NCNMNOISEEXPANDPARAM Expand null category noise model's structure from param vector.
% FORMAT
% DESC returns a null category noise model structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG noise : the noise structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% noise structure.
% RETURN noise : noise structure with the given parameters in the
% relevant locations.
%
% SEEALSO : ncnmNoiseParamInit, ncnmNoiseExtractParam, noiseExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% NOISE

noise.bias = params(1:noise.numProcess);
noise.gamman = params(noise.numProcess+1);
if noise.gammaSplit
  noise.gammap = params(noise.numProcess+2);
else
  noise.gammap = noise.gamman;
end


text:
function dt_new = tree_getNewVersion(dt)

dt_new=dt;
cut_old = dt.cut;
dt_new.cut=cell(size(cut_old));
dt_new.cut= num2cell(dt.cut);
% dt.cut= num2cell(dt.cut);

ind = find(dt.var<0);
for k = 1:numel(ind)
    dt_new.cut{ind(k)} = dt_new.catsplit(cut_old(ind(k)), :);
end

author: zouchuhang

repo: LayoutNet


text:
function [eigVec, eigVal] = dtiSplitTensor(tensor)
% Derive eigenvector and eigenvalues from the volume of tensor data
%
% [eigVec, eigVal] = dtiSplitTensor(tensor)
%
% The input tensor array is in XxYxZx6xN format
% The order of the value in the 4th dimension is [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz])
%
% The routine returns a XxYxZx3x3xN volume of eigVec
% and                 a the XxYxZx3xN volume of eigVal arrays
%
% X,Y,Z are positions in the volume
% N is the number of subjects.
%
% SEE ALSO: dtiRebuildTensor
%
% HISTORY:
% 2003.12.08 ASH (armins@stanford.edu) Wrote it.
% 2004.02.03 DTM (merget@cs.stanford.edu) Sort eigVal from high to low
% 2004.02.17 ASH: added extra dimension for subjects
% 2005.01.06 ASH: truly added extra dimension for subjects
%
% The code is implemented in the mex file dtiSplitTensor.c. If you don't
% have a that file compiled for you system, then the (very slow) code below
% will be executed.  In one test, the following code ran in 2.85 minutes,
% and the compiled version ran in about 11 seconds.
%
% (c) Stanford VISTA Team 2003

disp('This function is mexified for speed- compile dtiSplitTensor.c.');

sz = size(tensor);
if (length(sz)<5),
    sz =[sz, 1];
end

% vec = zeros([sz(1:3), 3, 3, sz(5)]);
% val = zeros([sz(1:3), 3, sz(5)]);

h = mrvWaitbar(0, 'Computing tensors...');
for(x=1:sz(1))
    for(y=1:sz(2))
        for(z=1:sz(3))
            for(j=1:sz(5))
                D = [tensor(x, y, z, 1, j), tensor(x, y, z, 4, j), tensor(x, y, z, 5, j);
                     tensor(x, y, z, 4, j), tensor(x, y, z, 2, j), tensor(x, y, z, 6, j);
                     tensor(x, y, z, 5, j), tensor(x, y, z, 6, j), tensor(x, y, z, 3, j)];
                [vec, val] = eig(D);
                [val2, order] = sort(-diag(val));
                eigVec(x, y, z, :, :, j) = vec(:, order);
                eigVal(x, y, z, :, j) = -val2;
            end
        end
    end
    mrvWaitbar(x/sz(1),h);
end
close(h);

return;

author: vistalab

repo: vistasoft


text:
function [E_Img]   =  WNNM_DeNoising( N_Img, Par )

E_Img           = N_Img;                                                        % Estimated Image
[Height Width]  = size(E_Img);
TotalPatNum     = (Height-Par.patsize+1)*(Width-Par.patsize+1);                 %Total Patch Number in the image
Dim             = Par.patsize*Par.patsize;

[Neighbor_arr Num_arr Self_arr] =	NeighborIndex(N_Img, Par);                  % PreCompute the all the patch index in the searching window
            NL_mat              =   zeros(Par.patnum,length(Num_arr));          % NL Patch index matrix
            CurPat              =	zeros( Dim, TotalPatNum );
            Sigma_arr           =   zeros( 1, TotalPatNum);
            EPat                =   zeros( size(CurPat) );
            W                   =   zeros( size(CurPat) );

for iter = 1 : Par.Iter
    E_Img             	=	E_Img + Par.delta*(N_Img - E_Img);
    [CurPat Sigma_arr]	=	Im2Patch( E_Img, N_Img, Par );                      % image to patch and estimate local noise variance

    if (mod(iter-1,Par.Innerloop)==0)
        Par.patnum = Par.patnum-10;                                             % Lower Noise level, less NL patches
        NL_mat  =  Block_matching_WNNM(CurPat, Par, Neighbor_arr, Num_arr, Self_arr);% Caculate Non-local similar patches for each
        if(iter==1)
            Sigma_arr = Par.nSig * ones(size(Sigma_arr));                       % First Iteration use the input noise parameter
        end
    end

     [EPat, W]  =  PatEstimation( NL_mat, Self_arr, Sigma_arr, CurPat, Par );   % Estimate all the patches
     E_Img      =  Patch2Im( EPat, W, Par.patsize, Height, Width );
end
return;

author: lbasek

repo: image-denoising-benchmark


text:
function line_num = sphere_llq_grid_line_count ( lat_num, long_num )

%*****************************************************************************80
%
%% SPHERE_LLQ_GRID_LINE_COUNT counts lines for an LLQ grid on a sphere.
%
%  Discussion:
%
%    A SPHERE LLQ grid imposes a grid of quadrilaterals on a sphere,
%    using latitude and longitude lines.
%
%    The number returned is the number of pairs of points to be connected
%    to form all the line segments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LAT_NUM, LONG_NUM, the number of latitude and
%    longitude lines to draw.  The latitudes do not include the North and South
%    poles, which will be included automatically, so LAT_NUM = 5, for instance,
%    will result in points along 7 lines of latitude.
%
%    Output, integer LINE_NUM, the number of grid lines.
%
  line_num = long_num * ( lat_num + 1 ) ...
           + lat_num  * long_num;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [lin_fun, max_lin_fun] = RobertsonCRF(stack, stack_exposure, max_iterations, err_threshold, bNormalize)
%
%       [lin_fun, max_lin_fun] = RobertsonCRF(stack, stack_exposure, max_iterations, err_threshold, bNormalize)
%
%       This function computes camera response function using Mitsunaga and
%       Nayar method.
%
%        Input:
%           -stack: a stack of LDR images. If the stack is a single or
%           double values are assumed to be in [0,1].
%           -stack_exposure: an array containg the exposure time of each
%           image. Time is expressed in second (s)
%           -max_iterations: max number of iterations
%           -err_threshold: threshold after which the function can stop
%           -bNormalize: if 1 it enables function normalization
%
%        Output:
%           -lin_fun: tabled function
%           -max_lin_fun: maximum value of the inverse CRF
%
%     Copyright (C) 2016  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(isempty(stack))
    error('RobertsonCRF: a stack cannot be empty!');
end

if(isempty(stack_exposure))
    error('RobertsonCRF: a stack_exposure cannot be empty!');
end

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

if(~exist('bNormalize', 'var'))
    bNormalize = 0;
end

if(~exist('max_iterations', 'var'))
    max_iterations = 15;
end

if(isa(stack, 'double'))
    stack = uint8(ClampImg(round(stack * 255), 0, 255));
end

if(isa(stack, 'single'))
    stack = uint8(ClampImg(round(stack * 255), 0, 255));
end

if(isa(stack, 'uint16'))
    stack = uint8(stack / 255);
end

col = size(stack, 3);

lin_fun = zeros(256, col);

for i=1:col
    lin_fun(:, i) = (0:255) / 255;
end

global minM;
global maxM;

x = (0:255) / 255;

w = WeightFunction(x, 'Robertson', 0);

minM = 0;
for m=0:255
    if(w(m + 1) > 0.0)
        minM = m;
        break;
    end
end

maxM = 255;
for m=255:-1:0
    if(w(m + 1) > 0.0)
        maxM = m;
        break;
    end
end

global stack_min;
global stack_max;

[~, t_min] = min(stack_exposure);
[~, t_max] = max(stack_exposure);
stack_min = ClampImg(single(stack(:,:,:,t_min)) / 255.0, 0.0, 1.0);
stack_max = ClampImg(single(stack(:,:,:,t_max)) / 255.0, 0.0, 1.0);

global minSatTime;

minSatTime = 1e10 * ones(col, 1);

for j=1:col
    for i=1:length(stack_exposure)
        tmp = stack(:,:,j,i);

        if(~isempty(find(tmp == 255)))
            if(stack_exposure(i) < minSatTime(j))
                minSatTime(j) = stack_exposure(i);
            end
        end
    end
end

for i=1:max_iterations
    lin_fun_prev = lin_fun;

    %normalization CRF
    lin_fun = Normalization(lin_fun);

    %update X
    x_tilde = Update_X(stack, stack_exposure, lin_fun);

    %update the iCRF
    lin_fun = Update_lin_fun(x_tilde, stack, stack_exposure, lin_fun);

    %compute error
    delta = (lin_fun_prev - lin_fun).^2;
    err = mean(delta(:));

    %disp([i, err]);

    if(err < err_threshold)
        break;
    end
end

max_lin_fun = max(lin_fun(:));

if(bNormalize)
    for i=1:col
        lin_fun(:,i) = lin_fun(:,i) / max_lin_fun;
    end

    lin_fun = ClampImg(lin_fun, 0.0, 1.0);
end

% %poly-fit (rational)
% if(bPolyFit)
%     ft = fittype( 'rat33' );
%     opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
%     opts.Display = 'Off';
%     opts.StartPoint = ones(7, 1);
%
%     for i=1:col
%         [xData, yData] = prepareCurveData(x', lin_fun(:,i));
%         [fit_ret, ~] = fit( xData, yData, ft, opts );
%         lin_fun(:,i) = feval(fit_ret, x');
%     end
% end

end

function lin_fun = Normalization(lin_fun)
    col = size(lin_fun, 2);

    for j=1:col
        lf = lin_fun(:, j);

        index_min = 1;
        for minIndx=1:256
            if(lf(minIndx) > 0.0)
                index_min = minIndx;
                break;
            end
        end

        index_max = 256;
        for maxIndx=256:-1:1
            if(lf(maxIndx) > 0.0)
                index_max = maxIndx;
                break;
            end
        end

        index = index_min + round((index_max - index_min) / 2);

        mid = lf(index);

        if(mid == 0.0)
            while(index < 256 && lf(index) == 0.0)
                index = index + 1;
            end

            mid = lf(index);
        end

        if(mid > 0.0)
            lin_fun(:,j) = lin_fun(:,j) / mid;
        end
    end
end

function imgOut = Update_X(stack, stack_exposure, lin_fun)
    global stack_min;
    global stack_max;
    global minM;
    global maxM;

    [r, c, col, n] = size(stack);

    %for each LDR image...
    imgOut    = zeros(r, c, col, 'single');
    totWeight = zeros(r, c, col, 'single');

    max_t = -ones(r, c, col);
    min_t = 1e10 * ones(r, c, col);

    for i=1:n
        t = stack_exposure(i);
        m = stack(:,:,:,i);

        indx = find(m > maxM);
        min_t(indx) = min(min_t(indx), t);

        indx = find(m < minM);
        max_t(indx) = max(max_t(indx), t);

        %normalizing m values in [0,1]
        tmpStack = ClampImg(single(m) / 255.0, 0.0, 1.0);

        %computing the weight function
        weight  = WeightFunction(tmpStack, 'Robertson', 0);

        tmpStack = tabledFunction(tmpStack, lin_fun);

        %summing things up...
        indx = find(tmpStack > stack_min & tmpStack < stack_max);

        if(t > 0.0 && ~isempty(indx))
            imgOut(indx) = imgOut(indx) + (weight(indx) .* tmpStack(indx)) * t;
            totWeight(indx) = totWeight(indx) + weight(indx) * t * t;
        end
    end

    imgOut = imgOut ./ totWeight;

    %taking care of saturated pixels
    saturation = 1e-4;
    imgOut(totWeight < saturation & totWeight > 0.0) = -1.0;

    for i=1:col
        io = imgOut(:,:,i);
        tw = totWeight(:,:,i);
        mxt = max_t(:,:,i);
        mnt = min_t(:,:,i);

        indx = find(tw == 0.0 & mxt > -1.0);
        io(indx) = lin_fun(minM, i) ./ mxt(indx);

        indx = find(tw == 0.0 & mnt < 1e10);
        io(indx) = lin_fun(maxM, i) ./ mnt(indx);

        imgOut(:,:,i) = io;
    end

end

function f_out = Update_lin_fun(x_tilde, stack, stack_exposure, lin_fun)
    col = size(x_tilde, 3);

    n = length(stack_exposure);
    f_out = zeros(size(lin_fun));

    global minSatTime;

    for i=1:col
        x_tilde_i = x_tilde(:,:,i);

        cardEm = zeros(256, 1);
        sumEm = zeros(256, 1);

        %for m values in [0,254]
        for m=0:254
            mp = m + 1;

            for k=1:n
                t = stack_exposure(k);

                tmp = stack(:,:,i,k);

                x = x_tilde_i((tmp == m) & (x_tilde_i > 0.0));

                sumEm(mp) = sumEm(mp) + t * sum(x(:));

                cardEm(mp) = cardEm(mp) + length(x(:));
            end
        end

        %for m = 255
        cardEm(256) = 1;

        m = 255;
        sumEm(256) = 1e10;
        for k=1:n
        	t = stack_exposure(k);

            tmp = stack(:,:,i,k);

            x = x_tilde_i((tmp == m) & (x_tilde_i > 0.0));

            if(~isempty(x) & (t == minSatTime))
               sumEm(256) = min([sumEm(256); t * x(:)]);
            end
        end

        %computing f_out + filling
        prev = 0.0;
        for m=1:256
            if(cardEm(m) > 0.0)
                f_out(m,i) = sumEm(m) / cardEm(m);
                prev = f_out(m,i);
            else
                f_out(m,i) = prev;
            end
        end
    end
end

author: banterle


text:
% SP_REFINE: construct a refined space from a given one. The function only
%                refines the space, the mesh must be refined separately.
%
%     [sp_fine, Proj] = sp_refine (space, msh, nsub, [degree], [regularity]);
%
% The same number of subdivisions, degree and regularity is applied to every patch
%
% INPUTS:
%
%     space:      the coarse space, an object of the sp_multipatch class (see sp_multipatch)
%     msh:        an object of the msh_multipatch class, usuallly the refined mesh (see msh_multipatch)
%     nsub:       number of uniform subdivisions to apply on each knot span, and for each direction
%     degree:     degree of the fine space, and for each direction
%     regularity: regularity for the new space, and for each direction
%
% OUTPUT:
%
%     sp_fine: the refined space, an object of the class sp_multipatch (see sp_multipatch)
%     Proj:    the coefficients relating 1D splines of the coarse and the fine spaces for each patch.
%                A cell-array of size 1 x npatch, each entry containing the
%                coefficients for the patch (either for scalar or vector-valued spaces).
%
% Copyright (C) 2015, 2016, 2018 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [sp_fine, Proj] = sp_refine (space, msh, nsub, varargin)

  if (space.ndof == sum (space.ndof_per_patch))
    gluing = false;
  else
    gluing = true;
  end

  sp_fine = cell (1, space.npatch);
  Proj = cell (1, space.npatch);
  for iptc = 1:space.npatch
    if (nargout == 2)
      [sp_fine{iptc}, Proj{iptc}] = sp_refine (space.sp_patch{iptc}, msh.msh_patch{iptc}, nsub, varargin{:});
    else
      sp_fine{iptc} = sp_refine (space.sp_patch{iptc}, msh.msh_patch{iptc}, nsub, varargin{:});
    end
  end
  if (~isempty (space.boundary))
    sp_fine = sp_multipatch (sp_fine, msh, space.interfaces, space.boundary.interfaces, gluing);
  else
    sp_fine = sp_multipatch (sp_fine, msh, space.interfaces, [], gluing);
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function dip = fixdipole(dip)

% FIXDIPOLE ensures that the dipole position and moment are
% consistently represented throughout FieldTrip functions.

% Copyright (C) 2009, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

[m, n] = size(dip.pos);

if n==3
  % the input representation is Nx3, which is what we want
elseif m==3
  % it is possible to translate it into a Nx3 unambiguously
  ft_warning('input dipole positions should be specified as Nx3 matrix');
  dip.pos = dip.pos';
elseif m==1
  % it is possible to translate it into a Nx3 unambiguously
  ft_warning('input dipole positions should be specified as Nx3 matrix');
  dip.pos = reshape(dip.pos, 3, n/3)';
else
  % it is not clear how to convert to a Nx3 matrix
  ft_error('input dipole positions should be specified as Nx3 matrix');
end

if isfield(dip, 'mom')
  ndip = size(dip.pos,1);
  if numel(dip.mom)==ndip*3
    ntime = 1;
  else
    ntime = numel(dip.mom)/(ndip*3);
  end
  dip.mom = reshape(dip.mom, ndip*3, ntime);
else
  dip.mom = [];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Standard ODFs
% predefined ODFs
%%
%


text:
function tc = tc_visualizeGlm_sparklines(tc);
%
% tc = tc_visualizeGlm_sparklines(tc);
%
% Visualize results of a general linear model on a time
% course, using the sparkline method. (See Edward Tufte,
% "Beautiful Evidence".)
%
% This is a test to see whether this would work.
%
% ras. 01/2007.
if notDefined('tc')
    tc = get(gcf, 'UserData');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apply a GLM if needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(tc, 'glm')
    tc = tc_applyGlm(tc);
end
X = tc.glm.designMatrix;
Y = tc.wholeTc(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up existing objects in figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
otherAxes = findobj('Type', 'axes','Parent', tc.ui.fig);
delete(otherAxes);
otherUiControls = findobj('Type', 'uicontrol', 'Parent', tc.ui.fig);
delete(otherUiControls);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% params
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nConds = sum(tc.trials.condNums>0);
nRows = nConds + 1;
rowSize = min(.15, 1/(nRows+1)); % normalized height of rows
t = [1:length(tc.wholeTc)] .* tc.TR; % time poitns for wholeTc

lo = min(tc.wholeTc); % min and max time series values
hi = max(tc.wholeTc);

if isnumeric(tc.params.glmHRF)
    opts = {sprintf('Mean trial for conditions \n %s',num2str(tc.params.snrConds)), ...
        'Boynton gamma function', 'SPM difference-of-gammas' ...
        'Dale & Buckner ''97'};
    hrfName = opts{tc.params.glmHRF};
else
    hrfName = tc.params.glmHRF; hrfName(hrfName=='_') = ' ';
end
title({'HRF function used: ' hrfName}, 'FontWeight', 'bold')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for row = 1:nRows
   subplot('Position', [.3 1-row*rowSize .65 rowSize]);

   if row==1        % plot time series data
       plot(t, tc.wholeTc, 'k');
   else             % plot predictor
       plot(t, tc.glm.designMatrix(:,row-1) .* tc.glm.betas(:,row-1), ...
            'Color', tc.trials.condColors{row});
   end
   axis([t(1) t(end) lo hi]);
   axis off

   % label condition
   if row==1
       text(-5, [lo + .3*(hi-lo)], 'Data');
   else
       text(-5, [lo + .3*(hi-lo)], tc.trials.condNames{row});
   end
end

return


text:
%
% Input
%  I: binary image (black = true)
%
function plot_image_only(I)
    assert(size(I,1)==105);
    assert(size(I,2)==105);
    newI = 1-I;
    newI(isinf(I)) = 0.5;
    image([1 105],[1 105],repmat(newI,[1 1 3]));
    set(gca,'YDir','reverse','XTick',[],'YTick',[]);
    xlim([1 105]);
    ylim([1 105]);
end

author: brendenlake

repo: omniglot


text:
function area = quad_area2_2d ( q )

%*****************************************************************************80
%
%% QUAD_AREA2_2D computes the area of a quadrilateral in 2D.
%
%  Discussion:
%
%    A quadrilateral is a polygon defined by 4 vertices.
%
%    This algorithm computes the area of the related
%    Varignon parallelogram first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(2,4), the vertices, specified in
%    counter clockwise order.
%
%    Output, real AREA, the area of the quadrilateral.
%

%
%  Define a parallelogram by averaging consecutive vertices.
%
  p(1:2,1:3) = ( q(1:2,1:3) + q(1:2,2:4) ) / 2.0;
  p(1:2,  4) = ( q(1:2,  4) + q(1:2,1  ) ) / 2.0;
%
%  Compute the area.
%
  area = parallelogram_area_2d ( p );
%
%  The quadrilateral's area is twice that of the parallelogram.
%
  area = 2.0 * area;

  return
end


text:
function [serverResponse] = buildFluxDistLayout( minerva, model, solution, identifier, hexColour, maxThickness, content)
% Builds a layout for MINERVA from a flux distribution. If a dictionary
% of identifiers is not provided it is assumed that the map and the COBRA
% model's nomenclature is coherent. Sends the layout to the remote MINERVA
% instance
%
% USAGE:
%
%    [serverResponse] = buildFluxDistLayout( minerva, model, solution, identifier, hexColour, thickness, content)
%
% INPUTS:
%    minerva:           Struct with the information of minerva instance:
%                       address, login, password and model (map)
%    model:             COBRA model structure
%    solution.v:        optimizeCb solution structure with a flux vector
%    identifier:        Name for the layout in MINERVA
%
% OPTIONAL INPUT:
%    hexColour          colour of overlay (hex color format)
%                       e.g. '#009933' corresponds to http://www.color-hex.com/color/009933
%                       If you want to make a color gradient, you can input
%                       an array of 2 or 3 colors like ["#ff0000", "#6617B5", "#0000ff"]
%                       note that they should be declared with (") rather
%                       than with (')
%    maxThickness:      maximum thickness
%    content:           character array with the following format for each
%                       reaction to be displayed. Bypasses the use of solution.v to set the format.
%                       'name%09reactionIdentifier%09lineWidth%09color%0D'
%
% OUTPUT:
%    serverResponse:          Response of the MINERVA
%
% .. Author: - Alberto Noronha Jan/2016
%            - Ines Thiele April/2020, fixed issue with using ReconMap-3 as target map.

if ~exist('thickness', 'var')
    maxThickness = 10;
end

useThickness = true; % flag to change thickness according to the flux

if ~exist('hexColour','var')
    defaultColor = '#57c657';
else
    hexColour = convertStringsToChars(hexColour);
    if ischar(hexColour)
        defaultColor = hexColour;
    elseif length(hexColour) == 1
        defaultColor = hexColour{1};
    else
        useThickness = false;

        if length(hexColour) >= 2
            cmap = makeColorGradient(hexColour{2}, hexColour{1}, maxThickness + 1);
        end

        if length(hexColour) >= 3
            ncmap = makeColorGradient(hexColour{3}, hexColour{2}, maxThickness + 1);
        end
    end
end

%nRxn=length(solution.v);
%normalizedFluxes = min(ones(nRxn,1),normalizeFluxes(abs(solution.v))-8);

% build input data for minerva
if ~exist('content','var')
    normalizedFluxes = normalizeFluxes(abs(solution.v), maxThickness);
    content = 'name%09reactionIdentifier%09lineWidth%09color%0D';

    for i=1:length(solution.v)
        mapReactionId = model.rxns{i};

        % if not ReconMap 2.01 use new reaction notation
        if ~strcmp(minerva.map, 'ReconMap-2.01')
            mapReactionId = strcat('R_', mapReactionId);
        end

        if solution.v(i) ~= 0

            if useThickness
                thickness = normalizedFluxes(i);
                color = defaultColor;
            else
                thickness = 1;
                if solution.v(i) < 0 && exist('ncmap','var')
                    color = ncmap{round(normalizedFluxes(i)) + 1};
                else
                    color = cmap{round(normalizedFluxes(i)) + 1};
                end
            end

            line = strcat('%09', mapReactionId, '%09', num2str(thickness), '%09', color, '%0D');
            content = strcat(content, line);
        end

        if contains(mapReactionId,'%') || contains(mapReactionId,' ')
            error('ReactionID cannot contain delimiting characters, such as: %')
        end
    end
end

%   get all the parameters
login = minerva.login;
password = minerva.password;
googleLicenseContent = minerva.googleLicenseConsent;
map = minerva.map;
%     have to turn it into string
%     disp(content);
% content = sprintf(content);
serverResponse = postMINERVArequest(login, password, map, googleLicenseContent, identifier, content);
if isempty(serverResponse)
    warning('Minerva server did not respond')
end

%end of function
end

%% Normalize a flux into a range of 1 to 10
function [ normalized_value ] = normalizeFluxes(fluxDistribution, thickness)

if ~exist('thickness','var') || nargin < 2
    thickness = 8;
end

m = min(fluxDistribution);
range = max(fluxDistribution) - m;
fluxDistribution = (fluxDistribution - m) / range;
range2 = - thickness;
normalized_value = (fluxDistribution*range2) + thickness;

end

author: opencobra

repo: cobratoolbox


text:
function [x, funVal, ValueL]=tree_mtLeastR(A, y, z, opts)
%
%%
% Function mcLeastR:
%      Least Squares Loss for Multi-task Learning
%             with the tree structured group Lasso Regularization
%
%% Problem
%
%  min  1/2 sum_j || A_j x_j - y_j||^2 + z * sum_i sum_j w_j ||x^ij_{G_ji}||
%
%  x^i denotes the i-th row of x
%  x_j denotes the j-th column of x
%  y_j denotes the j-th column of y
%
%  G_i's are nodes with tree structure
%
%  We assume that the tasks are of a tree structure
%
%  The tree structured group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%
%  opts.ind_t(1,:) contains the starting index
%  opts.ind_t(2,:) contains the ending index
%  opts.ind_t(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x^j is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind_t(1,1)=-1; opts.ind_t(2,1)=-1; and opts.ind_t(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x^j( opts.ind_t(1,j):opts.ind_t(2,j) ) denotes x^j_{G_j}. In this case,
%  the entries in opts.ind_t(1:2,:) are within 1 and k (the number of tasks).
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x^j( opts.G ( opts.ind_t(1,j):opts.ind_t(2,j) ) ) denotes x^j_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind_t(1:2,:) are within 1 and length(opts.G).
%
% The following example shows how G and ind works:
%
% G={ {1, 2}, {4, 5}, {3, 6}, {7, 8},
%     {1, 2, 3, 6}, {4, 5, 7, 8},
%     {1, 2, 3, 4, 5, 6, 7, 8} }.
%
% ind_t={ [1, 2, 100]', [3, 4, 100]', [5, 6, 100]', [7, 8, 100]',
%       [9, 12, 100]', [13, 16, 100]', [17, 24, 100]' },
%
% where each node has a weight of 100.
%
%  3) we use opts.ind to denote the starting and ending indices for the
%  samples of different tasks.
%  Samples for the the first task is in
%    A ( (opts.ind(1)+1):opts.ind(2), :  )
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mxk)
%  z -        Regularization parameter (z >=0)
%  opts-      optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         Solution (of size n x k)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on October 3, 2010.
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts (.ind, .ind_t) should be specified!\n');
end

[m,n]=size(A);

if (size(y,1) ~=m)
    error('\n Check the length of y!\n');
end

if (z<0)
    error('\n z should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% Detailed initialization
%%
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize ind
if ~isfield(opts,'ind')
    error('\n In tree_mtLeastR, .ind should be specified');
else
    ind=opts.ind;
    k=length(ind)-1;

    if ind(k+1)~=m
        error('\n Check opts.ind');
    end
end

% Initialize ind
if (~isfield(opts,'ind_t'))
    error('\n In tree_mtLeastR, the field .ind_t should be specified');
else
    ind_t=opts.ind_t;

    if (size(ind_t,1)~=3)
        error('\n Check opts.ind_t');
    end
end

GFlag=0;
% if GFlag=1, we shall apply general_altra
if (isfield(opts,'G'))
    GFlag=1;

    G=opts.G;
    if (max(G) >k || max(G) <1)
        error('\n The input G is incorrect. It should be within %d and %d',1,k);
    end
end

%% Starting point initialization

ATy=zeros(n, k);
% compute AT y
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

    if (opts.nFlag==0)
        tt =A(ind_i,:)'*y(ind_i,1);
    elseif (opts.nFlag==1)
        tt= A(ind_i,:)'*y(ind_i,1) - sum(y(ind_i,1)) * mu';
        tt=tt./nu(ind_i,1);
    else
        invNu=y(ind_i,1)./nu(ind_i,1);
        tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
    end

    ATy(:,i)= tt;
end

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
%     if (z<0 || z>1)
%         error('\n opts.rFlag=1, and z should be in [0,1]');
%     end

    computedFlag=0;
    if (isfield(opts,'lambdaMax'))
        if (opts.lambdaMax~=-1)
            lambda=z*opts.lambdaMax;
            computedFlag=1;
        end
    end

    if (~computedFlag)
        if (GFlag==0)
            lambda_max=findLambdaMax_mt(ATy, n, k, ind_t, size(ind_t,2));
        else
            lambda_max=general_findLambdaMax_mt(ATy, n, k, G, ind_t, size(ind_t,2));
        end

        % As .rFlag=1, we set lambda as a ratio of lambda_max
        lambda=z*lambda_max;
    end
end

% The following is for computing lambdaMax
% we use opts.lambdaMax=-1 to show that we need the computation.
%
% One can use this for setting up opts.lambdaMax
if (isfield(opts,'lambdaMax'))

    if (opts.lambdaMax==-1)
        if (GFlag==0)
            lambda_max=findLambdaMax_mt(ATy, n, k, ind_t, size(ind_t,2));
        else
            lambda_max=general_findLambdaMax_mt(ATy, n, k, G, ind_t, size(ind_t,2));
        end

        x=lambda_max;
        funVal=lambda_max;
        ValueL=lambda_max;

        return;
    end
end

% initialize a starting point
if opts.init==2
    x=zeros(n,k);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n || size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

Ax=zeros(m,1);
% compute Ax: Ax_i= A_i * x_i
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
    m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

    if (opts.nFlag==0)
        Ax(ind_i,1)=A(ind_i,:)* x(:,i);
    elseif (opts.nFlag==1)
        invNu=x(:,i)./nu; mu_invNu=mu * invNu;
        Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
    else
        Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
        Ax(ind_i,1)=Ax./nu(ind_i,1);
    end
end

if (opts.init==0)
    % ------  This function is not available
    %
    % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail
    %

    % Here, we only support starting from zero, due to the complex tree
    % structure

    x=zeros(n,k);
end

%% The main program

%% The Armijo Goldstein line search schemes + accelearted gradient descent

if (opts.mFlag==0 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute ATAs : n x k
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

            if (opts.nFlag==0)
                tt =A(ind_i,:)'*As(ind_i,1);
            elseif (opts.nFlag==1)
                tt= A(ind_i,:)'*As(ind_i,1) - sum(As(ind_i,1)) * mu';
                tt=tt./nu(ind_i,1);
            else
                invNu=As(ind_i,1)./nu(ind_i,1);
                tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
            end

            ATAs(:,i)= tt;
        end

        % obtain the gradient g
        g=ATAs-ATy;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L;

            % tree overlapping group Lasso projection
            ind_work(1:2,:)=ind_t(1:2,:);
            ind_work(3,:)=ind_t(3,:) * (lambda / L);

            if (GFlag==0)
                x=altra_mt(v, n, k, ind_work, size(ind_work,2));
            else
                x=general_altra_mt(v, n, k, G, ind_work, size(ind_work,2));
            end

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax: Ax_i= A_i * x_i
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

                if (opts.nFlag==0)
                    Ax(ind_i,1)=A(ind_i,:)* x(:,i);
                elseif (opts.nFlag==1)
                    invNu=x(:,i)./nu; mu_invNu=mu * invNu;
                    Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
                else
                    Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
                    Ax(ind_i,1)=Ax./nu(ind_i,1);
                end
            end

            Av=Ax -As;
            r_sum=norm(v,'fro')^2; l_sum=Av'*Av;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=Axy'*Axy/2;

        for i=1:n
            xRow=x(i,:);

            if (GFlag==0)
                tree_norm=treeNorm(xRow, k, ind_t, size(ind_t,2));
            else
                tree_norm=general_treeNorm(xRow, k, G, ind_t, size(ind_t,2));
            end

            funVal(iterStep)=funVal(iterStep)+ lambda*tree_norm;
        end

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end

        % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp=x; Axp=Ax; xxp=zeros(n,k); L =L/2;
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
end

author: thomaskuestner


text:
function numdiff = ndiff(A)
%
% Accounts for the number of different elements that compose the matrix A

B = A(1);
for i=2:numel(A)
    if isempty(find(A(i) == B))
        B = [B A(i)];
    end
end
numdiff = length(B);