text:
function []=ps_select(reest_flag,plot_flag)
% PS_SELECT select PS based on gamma and D_A
%   PS_SELECT(REEST_FLAG,PLOT_FLAG)
%   REEST_FLAG = 1 skip reestimation completely (default 0)
%   REEST_FLAG = 2 reuse previously reestimated gamma (default 0)
%   REEST_FLAG = 3 reestimate gamma for all candidate pixels (default 0)
%   PLOT_FLAG = 1 for plots (default 0)
%
%   Andy Hooper, June 2006
%
%   ======================================================================
%   07/2006 AH: Use specific bperp for re-estimation of gamma
%   08/2006 AH: Load workspaces into structures
%   09/2006 AH: Check if no gamma threshold meets criteria
%   09/2006 AH: add small baselines
%   09/2009 AH: add option to select on density instead of percentage
%   01/2010 KS: fixed badly conditioned polyfit errors by scaling and
%               centering
%   02/2010 AH: Only bin by D_A if enough candidate pixels
%   02/2010 AH: Leave out ifgs in drop_ifg_index from noise calculation
%   09/2010 JJS+MA: To make use oversampling dataset
%   12/2010 AH: Fix error message for density selection
%   02/2011 DB: Fix error with keep_ix
%   05/2012 AH: subtract only pixel being tested, not zero whole grid cell
%   01/2013 AH: Set default threshold if not enough random pixels
%   04/2015 DB: Give a warning to remove patch from patch list when no PS are left
%   09/2015 DB: Store information on nr of PS left. Support auto procesing
%   01/2016 DB: Replace save with stamps_save
%   02/2016 DB: Identified bug when patch size is smaller than filter size.
%               For now drop this patch. This needs to be fixed better.
%   08/2017 AH: Ensure coh_thresh not negative.
%   ======================================================================
logit;
logit('Selecting stable-phase pixels...')

if nargin<1
    reest_flag=0;
end

if nargin<2
    plot_flag=0;
end

psver=getpsver;
if psver>1
   setpsver(1)
end

slc_osf=getparm('slc_osf',1);  % [MA]
clap_alpha=getparm('clap_alpha',1);
clap_beta=getparm('clap_beta',1);
n_win=getparm('clap_win',1);
select_method=getparm('select_method',1);
if strcmpi(select_method,'PERCENT')
    max_percent_rand=getparm('percent_rand',1);
else
    max_density_rand=getparm('density_rand',1);
end
gamma_stdev_reject=getparm('gamma_stdev_reject',1);
small_baseline_flag=getparm('small_baseline_flag',1);
drop_ifg_index=getparm('drop_ifg_index',1);

if strcmpi(small_baseline_flag,'y')
    low_coh_thresh=15; % equivalent to coh of 15/100
else
    low_coh_thresh=31; % equivalent to coh of 31/100
end

load psver

psname=['ps',num2str(psver)];
phname=['ph',num2str(psver)];
pmname=['pm',num2str(psver)];
selectname=['select',num2str(psver)];
daname=['da',num2str(psver)];
bpname=['bp',num2str(psver)];

ps=load(psname);

ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);

if exist([phname,'.mat'],'file')
    phin=load(phname);
    ph=phin.ph;
    clear phin
else
    ph=ps.ph;
end

bperp=ps.bperp;
n_ifg=ps.n_ifg;
if ~strcmpi(small_baseline_flag,'y')
    master_ix=sum(ps.master_day>ps.day)+1;
    no_master_ix=setdiff([1:ps.n_ifg],ps.master_ix);
    ifg_index=setdiff(ifg_index,ps.master_ix);
    ifg_index(ifg_index>master_ix)=ifg_index(ifg_index>master_ix)-1;
    ph=ph(:,no_master_ix);
    bperp=bperp(no_master_ix);
    n_ifg=length(no_master_ix);
end
n_ps=ps.n_ps;
xy=ps.xy;

pm=load(pmname);
if exist([daname,'.mat'],'file')
    da=load(daname);
    D_A=da.D_A;
    clear da
else
    D_A=[];
end

if ~isempty(D_A) & size(D_A,1)>=10000
    % chunk up PSC
    D_A_sort=sort(D_A);
    if size(D_A,1)>=50000
      bin_size=10000;
    else
      bin_size=2000;
    end
    D_A_max=[0;D_A_sort(bin_size:bin_size:end-bin_size);D_A_sort(end)];

else
    D_A_max=[0;1];
    D_A=ones(size(pm.coh_ps));
end

if ~strcmpi(select_method,'PERCENT')
    patch_area=prod(max(xy(:,2:3),[],1)-min(xy(:,2:3),[],1))/1e6; % in km
    max_percent_rand=max_density_rand*patch_area/(length(D_A_max)-1);
end

min_coh=zeros(length(D_A_max)-1,1);
D_A_mean=zeros(size(D_A_max,1)-1,1);
Nr_dist=pm.Nr;

if reest_flag==3 % reestimate for all candidate pixels
    coh_thresh=0;
    coh_thresh_coeffs=[];
else
  for i=1:length(D_A_max)-1
    coh_chunk=pm.coh_ps(D_A>D_A_max(i) & D_A<=D_A_max(i+1));
    D_A_mean(i)=mean(D_A(D_A>D_A_max(i) & D_A<=D_A_max(i+1)));
    coh_chunk=coh_chunk(coh_chunk~=0);  % discard PSC for which coherence was not calculated
    Na=hist(coh_chunk,pm.coh_bins);
    Nr=Nr_dist*sum(Na(1:low_coh_thresh))/sum(Nr_dist(1:low_coh_thresh));
    if i==length(D_A_max)-1 & plot_flag==1
        figure
        plot(pm.coh_bins,Na,'g')
        hold on
        plot(pm.coh_bins,Nr,'r')
        legend('data','random')
        title('Before Gamma Reestimation')
    end
    Na(Na==0)=1; % avoid divide by zero
    if strcmpi(select_method,'PERCENT')
        percent_rand=fliplr(cumsum(fliplr(Nr))./cumsum(fliplr(Na))*100);
    else
        percent_rand=fliplr(cumsum(fliplr(Nr))); % absolute number
    end
    ok_ix=find(percent_rand<max_percent_rand);
    if isempty(ok_ix)
        min_coh(i)=1; % no threshold meets criteria
    else
        min_fit_ix=min(ok_ix)-3;
        if min_fit_ix<=0
            %min_coh(i)=low_coh_thresh/100;
            min_coh(i)=nan;
        else
            max_fit_ix=min(ok_ix)+2;
            max_fit_ix(max_fit_ix>100)=100; % make sure not higher than length of percent_rand
            [p,S,mu]=polyfit(percent_rand(min_fit_ix:max_fit_ix),[min_fit_ix*0.01:0.01:max_fit_ix*0.01],3); % KS
            min_coh(i)=polyval(p,max_percent_rand,[],mu); % KS
%             p=polyfit(percent_rand(min_fit_ix:max_fit_ix),[min_fit_ix*0.01:0.01:max_fit_ix*0.01],3);
%             min_coh(i)=polyval(p,max_percent_rand);
        end
    end
  end

  nonnanix=~isnan(min_coh);
  if sum(nonnanix)<1
      warning('Not enough random phase pixels to set gamma threshold - using default threshold of 0.3')
      coh_thresh=0.3;
      coh_thresh_coeffs=[];
  else
  min_coh=min_coh(nonnanix);
  D_A_mean=D_A_mean(nonnanix);
  if size(min_coh,1)>1
    coh_thresh_coeffs=polyfit(D_A_mean,min_coh,1);  % fit polynomial to the curve
    if coh_thresh_coeffs(1)>0 % positive slope (as expected)
        coh_thresh=polyval(coh_thresh_coeffs,D_A);
    else % unable to ascertain correct slope
        coh_thresh=polyval(coh_thresh_coeffs,0.35); % set an average threshold for all D_A
        coh_thresh_coeffs=[];
    end
  else
    coh_thresh=min_coh;
    coh_thresh_coeffs=[];
  end
 end
end

coh_thresh(coh_thresh<0)=0; % to ensures pixels with coh=0 are rejected

logit(sprintf('Initial gamma threshold: %.3f at D_A=%.2f to %.3f at D_A=%.2f',min(coh_thresh),min(D_A),max(coh_thresh),max(D_A)))

if  plot_flag==1
    figure
    plot(D_A_mean,min_coh,'*')
    hold on
    if ~isempty(coh_thresh_coeffs)
        plot(D_A_mean,polyval(coh_thresh_coeffs,D_A_mean),'r')
    end
    ylabel('\gamma_{thresh}')
    xlabel('D_A')
end

ix=find(pm.coh_ps>coh_thresh); % select those below threshold
n_ps=length(ix);
logit(sprintf('%d PS selected initially',n_ps))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% reject part-time PS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if gamma_stdev_reject>0
    ph_res_cpx=exp(j*pm.ph_res(:,ifg_index));
    coh_std=zeros(length(ix),1);
    for i=1:length(ix)
        coh_std(i)=std(bootstrp(100,@(ph) abs(sum(ph))/length(ph), ph_res_cpx(ix(i),ifg_index)));
    end
    clear ph_res_cpx
    ix=ix(coh_std<gamma_stdev_reject);
    n_ps=length(ix);
    logit(sfprintf('%d PS left after pps rejection',n_ps))
end

if reest_flag~=1

  if reest_flag~=2 % reestimate coh with the PS removed from filtered patch
    for i=1:length(drop_ifg_index);
      if strcmpi(small_baseline_flag,'y')
        logit(sprintf('%s-%s is dropped from noise re-estimation',datestr(ps.ifgday(drop_ifg_index(i),1)),datestr(ps.ifgday(drop_ifg_index(i),2))))
      else
        logit(sprintf('%s is dropped from noise re-estimation',datestr(ps.day(drop_ifg_index(i)))))
      end
    end
    pm=rmfield(pm,{'ph_res'});
    pm=rmfield(pm,{'ph_patch'});
    ph_patch2=zeros(n_ps,n_ifg,'single');
    ph_res2=zeros(n_ps,n_ifg,'single');
    ph=ph(ix,:);

    if length(coh_thresh)>1
        coh_thresh=coh_thresh(ix);
    end

    n_i=max(pm.grid_ij(:,1));
    n_j=max(pm.grid_ij(:,2));
    K_ps2=zeros(n_ps,1);
    C_ps2=zeros(n_ps,1);
    coh_ps2=zeros(n_ps,1);
    ph_filt=zeros(n_win,n_win,n_ifg);

    for i=1:n_ps
        ps_ij=pm.grid_ij(ix(i),:);
        i_min=max(ps_ij(1)-n_win/2,1);
        i_max=i_min+n_win-1;
        if i_max>n_i
            i_min=i_min-i_max+n_i;
            i_max=n_i;
        end
        j_min=max(ps_ij(2)-n_win/2,1);
        j_max=j_min+n_win-1;
        if j_max>n_j
            j_min=j_min-j_max+n_j;
            j_max=n_j;
        end

        % it could occur that your patch size is smaller than the filter size
        % crude bug fix is to drop this patch. It needs fixing in future...
        if j_min<1 || i_min<1
            % THIS NEEDS TO BECOME AN ACTUAL FIX, but not sure how...
            ph_patch2(i,:) =0;
        else

            % remove the pixel for which the smoothign is computed
            ps_bit_i=ps_ij(1)-i_min+1;
            ps_bit_j=ps_ij(2)-j_min+1;
            ph_bit=pm.ph_grid(i_min:i_max,j_min:j_max,:);
            ph_bit(ps_bit_i,ps_bit_j,:)=0;

            %ph_bit(ps_bit_i,ps_bit_j,:)=ph_bit(ps_bit_i,ps_bit_j,:)-shiftdim(pm.ph_weight(i,:),-1);

            % JJS oversample update for PS removal + [MA] general usage update
            ix_i=ps_bit_i-(slc_osf-1):ps_bit_i+(slc_osf-1);
            ix_i=ix_i(ix_i>0&ix_i<=size(ph_bit,1));
            ix_j=ps_bit_j-(slc_osf-1):ps_bit_j+(slc_osf-1);
            ix_j=ix_j(ix_j>0&ix_j<=size(ph_bit,2));
            ph_bit(ix_i,ix_j)=0;

            for i_ifg=1:n_ifg
                ph_filt(:,:,i_ifg)=clap_filt_patch(ph_bit(:,:,i_ifg),clap_alpha,clap_beta,pm.low_pass);
            end

            ph_patch2(i,:)=squeeze(ph_filt(ps_bit_i,ps_bit_j,:));
        end
        if i/10000==floor(i/10000)
            logit(sprintf('%d patches re-estimated',i))
        end

    end

    pm=rmfield(pm,{'ph_grid'});
    bp=load(bpname);
    bperp_mat=bp.bperp_mat(ix,:);
    clear bp

    for i=1:n_ps
        psdph=ph(i,:).*conj(ph_patch2(i,:));
        if sum(psdph==0)==0  % insist on a non-null value in every ifg
            psdph=psdph./abs(psdph);
            [Kopt,Copt,cohopt,ph_residual]=ps_topofit(psdph(ifg_index),bperp_mat(i,ifg_index)',pm.n_trial_wraps,'n');
            K_ps2(i)=Kopt(1);
            C_ps2(i)=Copt(1);
            coh_ps2(i)=cohopt(1);
            ph_res2(i,ifg_index)=angle(ph_residual);
        else
            K_ps2(i)=nan;
            coh_ps2(i)=nan;
        end
        if i/10000==floor(i/10000)
            logit(sprintf('%d coherences re-estimated',i))
        end
    end
  else % reest_flag==2, use previously recalculated coh
    sl=load(selectname);
    ix=sl.ix;
    coh_ps2=sl.coh_ps2;
    K_ps2=sl.K_ps2;
    C_ps2=sl.C_ps2;
    ph_res2=sl.ph_res2;
    ph_patch2=sl.ph_patch2;
  end
    pm.coh_ps(ix)=coh_ps2;

    % calculate threshold again based on recalculated coh
    for i=1:length(D_A_max)-1
        coh_chunk=pm.coh_ps(D_A>D_A_max(i) & D_A<=D_A_max(i+1));
        D_A_mean(i)=mean(D_A(D_A>D_A_max(i) & D_A<=D_A_max(i+1)));
        coh_chunk=coh_chunk(coh_chunk~=0);  % discard PSC for which coherence was not calculated
        Na=hist(coh_chunk,pm.coh_bins);
        Nr=Nr_dist*sum(Na(1:low_coh_thresh))/sum(Nr_dist(1:low_coh_thresh));
        if i==length(D_A_max)-1 & plot_flag==1
            figure
            plot(pm.coh_bins,Na,'g')
            hold on
            plot(pm.coh_bins,Nr,'r')
            legend('data','random')
            title('After Gamma Reestimation')
        end
        Na(Na==0)=1; % avoid divide by zero
        if strcmpi(select_method,'PERCENT')
            percent_rand=fliplr(cumsum(fliplr(Nr))./cumsum(fliplr(Na))*100);
        else
            percent_rand=fliplr(cumsum(fliplr(Nr))); % despite the name, percent_rand here is the absolute number
        end
        ok_ix=find(percent_rand<max_percent_rand);
        if isempty(ok_ix)
            min_coh(i)=1;
        else
            min_fit_ix=min(ok_ix)-3;
            if min_fit_ix<=0
                %min_coh(i)=low_coh_thresh/100;
                min_coh(i)=nan;
            else
                max_fit_ix=min(ok_ix)+2;
                max_fit_ix(max_fit_ix>100)=100; % make sure not higher than length of percent_rand
                [p,S,mu]=polyfit(percent_rand(min_fit_ix:max_fit_ix),[min_fit_ix*0.01:0.01:max_fit_ix*0.01],3); % KS
                min_coh(i)=polyval(p,max_percent_rand,[],mu); % KS
%                 p=polyfit(percent_rand(min_fit_ix:max_fit_ix),[min_fit_ix*0.01:0.01:max_fit_ix*0.01],3);
%                 min_coh(i)=polyval(p,max_percent_rand);
            end
        end
    end

    nonnanix=~isnan(min_coh);
    if sum(nonnanix)<1
        coh_thresh=0.3;
        coh_thresh_coeffs=[];
    else
        min_coh=min_coh(nonnanix);
        D_A_mean=D_A_mean(nonnanix);
        if length(min_coh)>1
        %keyboard
            coh_thresh_coeffs=polyfit(D_A_mean,min_coh,1);  % fit polynomial to the curve
            if coh_thresh_coeffs(1)>0 % positive slope (as expected)
                coh_thresh=polyval(coh_thresh_coeffs,D_A(ix));
            else % unable to ascertain correct slope
                coh_thresh=polyval(coh_thresh_coeffs,0.35); % set an average threshold for all D_A
                coh_thresh_coeffs=[];
            end
        else
            coh_thresh=min_coh;
            coh_thresh_coeffs=[];
        end
    end

    coh_thresh(coh_thresh<0)=0; % to ensures pixels with coh=0 are rejected

    logit(sprintf('Reestimation gamma threshold: %.3f at D_A=%.2f to %.3f at D_A=%.2f',min(coh_thresh),min(D_A),max(coh_thresh),max(D_A)))

    bperp_range=max(bperp)-min(bperp);
    keep_ix=coh_ps2>coh_thresh & abs(pm.K_ps(ix)-K_ps2)<2*pi/bperp_range;
    clear pm
    logit(sprintf('%d ps selected after re-estimation of coherence',sum(keep_ix)))
    clear pm

else % reest_flag==1, skip re-estimation
    pm=rmfield(pm,{'ph_grid'});
    ph_patch2=pm.ph_patch(ix,:);
    ph_res2=pm.ph_res(ix,:);
    K_ps2=pm.K_ps(ix);
    C_ps2=pm.C_ps(ix);
    coh_ps2=pm.coh_ps(ix);
    keep_ix=true(size(ix));
end % end-if reest_flag

%%% Keep information about number of PS left.
if exist('no_ps_info.mat','file')~=2
   stamps_step_no_ps = zeros([5 1 ]);       % keep for the first 5 steps only
else
   load('no_ps_info.mat');
   % reset as we are currently re-processing
   stamps_step_no_ps(3:end)=0;
end
if sum(keep_ix)==0
   fprintf('***No PS points left. Updating the stamps log for this****\n')
   % update the flag indicating no PS left in step 3
   stamps_step_no_ps(3)=1;
end
save('no_ps_info.mat','stamps_step_no_ps')

if  plot_flag==1
    figure
    plot(D_A_mean,min_coh,'*')
    hold on
    if ~isempty(coh_thresh_coeffs)
        plot(D_A_mean,polyval(coh_thresh_coeffs,D_A_mean),'r')
    end
    ylabel('\gamma_{thresh}')
    xlabel('D_A')
end

% save(selectname,'ix','keep_ix','ph_patch2','ph_res2','K_ps2','C_ps2','coh_ps2','coh_thresh','coh_thresh_coeffs','clap_alpha','clap_beta','n_win','max_percent_rand','gamma_stdev_reject','small_baseline_flag','ifg_index');
stamps_save(selectname,ix,keep_ix,ph_patch2,ph_res2,K_ps2,C_ps2,coh_ps2,coh_thresh,coh_thresh_coeffs,clap_alpha,clap_beta,n_win,max_percent_rand,gamma_stdev_reject,small_baseline_flag,ifg_index);

logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
clear, clc;

% This is an example for running the function LeastR
%
%  min  1/2 || A x - y||^2 + 1/2 * rsL2 * ||x||_2^2 + rho * ||x||_1
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
% [3]  Jun Liu, Jianhui Chen, and Jieping Ye,
%      Large-Scale Sparse Logistic Regression, KDD, 2009.
%
%% ------------   History --------------------
%
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1;

m=1000;  n=1000;    % The data matrix is of size m x n

% for reproducibility
randNum=1;

% ---------------------- Generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);       % the data matrix

randn('state',(randNum-1)*3+2);
xOrin=randn(n,1);

randn('state',(randNum-1)*3+3);
noise=randn(m,1);
y=A*xOrin +...
    noise*0.01;     % the response
weight=0.001:0.001:1;%%*********************************************add weight to rho!!!!!!!!!!1
rho=0.5*weight';            % the regularization parameter
                    % it is a ratio between (0,1), if .rFlag=1

%----------------------- Set optional items ------------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% termination criterion
opts.tFlag=0;       % run .maxIter iterations
%opts.maxIter=100;   % maximum number of iterations

% normalization
opts.nFlag=0;       % without normalization

% regularization
opts.rFlag=1;       % the input parameter 'rho' is a ratio in (0, 1)
%opts.rsL2=0.01;     % the squared two norm term

%----------------------- Run the code LeastR -----------------------
fprintf('\n mFlag=0, lFlag=0 \n');
opts.mFlag=0;       % treating it as compositive function
opts.lFlag=0;       % Nemirovski's line search
tic;
[x1, funVal1, ValueL1]= LeastRweight(A, y, 0.5, opts);
[x11, funVal1, ValueL1]= LeastRweight(A, y, rho*2, opts);
[x111, funVal1, ValueL1]= LeastRweight(A, y, flipud(rho*2), opts);
toc;

opts.maxIter=1000;

fprintf('\n mFlag=1, lFlag=0 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=0;       % Nemirovski's line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x2, funVal2, ValueL2]= LeastRweight(A, y, rho, opts);
toc;

fprintf('\n mFlag=1, lFlag=1 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=1;       % adaptive line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x3, funVal3, ValueL3]= LeastRweight(A, y, rho, opts);
toc;

figure;
plot(funVal1,'-r');
hold on;
plot(funVal2,'--b');
hold on;
plot(funVal3,':g');
legend('mFlag=0, lFlag=0', 'mFlag=1, lFlag=0', 'mFlag=1, lFlag=1');
xlabel('Iteration (i)');
ylabel('The objective function value');

% % --------------------- compute the pathwise solutions ----------------
opts.fName='LeastR';      % set the function name to 'LeastR'
Z=[0.5, 0.2, 0.1, 0.01];  % set the parameters

% run the function pathSolutionLeast
fprintf('\n Compute the pathwise solutions, please wait...');
X=pathSolutionLeast(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
function s = getAlgorithmsHeaderSizeWordAligned(this)

% Assume a word is 4 bytes
s = getAlgorithmsHeaderSize(this);
s = ceil(s/4)*4;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [scales, scalesK] = hsvargplvmRetainedScales(model, thresh, displPlots, displ)

if nargin < 4 || isempty(displ)
    displ = false;
end

if nargin < 3 || isempty(displPlots)
    displPlots = true;
end

if nargin < 2 || isempty(thresh)
    thresh = 0.01;
end

for h=1:model.H
    scalesK{h} = zeros(model.layer{h}.M, model.layer{h}.q);
    scalesAll{h} = zeros(model.layer{h}.M, model.layer{h}.q);
    scalesAllTmp{h} = svargplvmScales('get', model.layer{h});

    for m=1:model.layer{h}.M
        scales{h}{m} = vargplvmRetainedScales(model.layer{h}.comp{m},thresh);
        scalesK{h}(m,scales{h}{m}) = 1;
        scalesAll{h}(m,:) = scalesAllTmp{h}{m};
        scalesAll{h}(m,:) = scalesAll{h}(m,:) / max(scalesAll{h}(m,:)); %% Scale so that 1 is max
    end
end

if displPlots
    for h=1:model.H
        if model.layer{h}.M > 1
            figure
            imagesc(scalesAll{h}); title(['All scales for layer ' num2str(h)]);     set(gca,'XGrid','on')
            figure
            imagesc(scalesK{h}); title(['Binarized scales for layer ' num2str(h)]);     set(gca,'XGrid','on')
        end
    end
end

if displ
    for h=1:model.H
        fprintf('# Layer %d\n', h)
        fprintf(' q  | Scales\n')
        fprintf('------------------\n')
        for q=1:model.layer{h}.M
            fprintf(' %d | %s \n',q,num2str(scales{h}{q}));
        end
        fprintf('\n');
    end
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function [h, t] = vl_plotbox(boxes, varargin)
% PLOTBOX  Plot boxes
%   PLOTBOX(BOX) plots a box specified by the 4-dimensional vector
%   BOXES = [XMIN YMIN XMAX YMAX]'. If BOXES is a 4 x N matrix, a box
%   for each of the N columns is plotted.
%
%   H = PLOTBOX(BOXES) returns an handle to the line drawings
%   representing the boxes. For multiple boxes, H is a row vector with
%   one handle per box.
%
%   PLOTBOX(BOXES, 'LABEL', LABEL) annotates the box with the string
%   LABEL. If BOXES contains multiple boxes, then LABEL can be a cell
%   array with one entry for each box. H is then a 2 x N array with
%   handles to boxes and labels.
%
%   PLOTBOX(BOXES, ...) passes any extra argument to the underlying
%   plotting function. The first optional argument can be a line
%   specification string such as the one used by PLOT().
%
%   See also:: VL_PLOTFRAME().

% Author:: Andrea Vedaldi

% Copyright (C) 2008-13 Andrea Vedaldi
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.label = {} ;

% if the first optional argument is a linespec expand it
if length(varargin) > 0
  lineprop = vl_linespec2prop(varargin{1}) ;
  varargin = {lineprop{:}, varargin{2:end}} ;
end

% parse optional arguments
[opts, varargin] = vl_argparse(opts, varargin) ;
if ischar(opts.label)
  opts.label = {opts.label} ;
end

if size(boxes,2) == 0
  h = [] ;
  return ;
end

if size(boxes,1) ~= 4
  error('BOXES must be a 4 x N matrix') ;
end

Lx = [1 0 0 0 ;
      0 0 1 0 ;
      0 0 1 0 ;
      1 0 0 0 ;
      1 0 0 0 ] ;

Ly = [0 1 0 0 ;
      0 1 0 0 ;
      0 0 0 1 ;
      0 0 0 1 ;
      0 1 0 0 ] ;

fig = newplot ;
h = line(Lx * boxes, Ly * boxes, varargin{:}) ;
h = h' ;

if ~isempty(opts.label)
  ish = ishold ;
  hold on ;
  t = zeros(1,length(opts.label)) ;
  for r = 1:size(boxes,2)
    cl = get(h(1,r), 'Color') ;
    q = mod(r - 1, length(opts.label)) + 1 ;
    h(2,r) = text(mean(boxes([1 3],r)), boxes(2,r), opts.label{q}, ...
                  'Background', cl, ...
                  'VerticalAlignment', 'bottom', ...
                  'HorizontalAlignment', 'center') ;
  end
  if ~ish, hold off ; end
end

if nargout == 0
  clear h ;
end

author: jianxiongxiao

repo: ProfXkit


text:
function angles = tetrahedron_face_angles_3d ( tetra )

%*****************************************************************************80
%
%% TETRAHEDRON_FACE_ANGLES_3D returns the 12 face angles of a tetrahedron 3D.
%
%  Discussion:
%
%    The tetrahedron has 4 triangular faces.  This routine computes the
%    3 planar angles associated with each face.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real TETRA(3,4) the tetrahedron vertices.
%
%    Output, real ANGLES(3,4), the face angles.
%

%
%  Face 123
%
  tri(1:3,1:3) = tetra(1:3,1:3);
  angles(1:3,1) = triangle_angles_3d ( tri );
%
%  Face 124
%
  tri(1:3,1:2) = tetra(1:3,1:2);
  tri(1:3,3) = tetra(1:3,4);
  angles(1:3,2) = triangle_angles_3d ( tri );
%
%  Face 134
%
  tri(1:3,1) = tetra(1:3,1);
  tri(1:3,2:3) = tetra(1:3,3:4);
  angles(1:3,3) = triangle_angles_3d ( tri );
%
%  Face 234
%
  tri(1:3,1:3) = tetra(1:3,2:4);
  angles(1:3,4) = triangle_angles_3d ( tri );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ unique_num, a1, a2 ] = r8vec2_sorted_unique ( n, a1, a2 )

%*****************************************************************************80
%
%% R8VEC2_SORTED_UNIQUE keeps unique elements in a sorted R8VEC2.
%
%  Discussion:
%
%    An R8VEC2 is two R8VEC's.
%
%    An R8VEC is a vector of R8 values.
%
%    Item I is stored as the pair A1(I), A2(I).
%
%    The items must have been sorted, or at least it must be the
%    case that equal items are stored in adjacent vector locations.
%
%    If the items were not sorted, then this routine will only
%    replace a string of equal values by a single representative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of items.
%
%    Input, real A1(N), A2(N), the array of N items.
%
%    Output, integer UNIQUE_NUM, the number of unique items.
%
%    Output, real A1(UNIQUE_NUM), A2(UNIQUE_NUM), the array of unique items.
%
  if ( n <= 0 )
    unique_num = 0;
    return
  end

  unique_num = 1;

  for itest = 2 : n

    if ( a1(itest) ~= a1(unique_num) || a2(itest) ~= a2(unique_num) )

      unique_num = unique_num + 1;

      a1(unique_num) = a1(itest);
      a2(unique_num) = a2(itest);

    end

  end

  a1 = a1(1:unique_num);
  a2 = a2(1:unique_num);

  return
end


text:
function hypersphere_properties_test01 ( )

%*****************************************************************************80
%
%% HYPERSPHERE_PROPERTIES_TEST01 tests the coordinate conversion routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERSPHERE_PROPERTIES_TEST01\n' );
  fprintf ( 1, '  Test the coordinate conversion routines:\n' );
  fprintf ( 1, '  CARTESIAN_TO_HYPERSPHERE: X       -> R,Theta\n' );
  fprintf ( 1, '  HYPERSPHERE_TO_CARTESIAN: R,Theta -> X.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pick a random X, and compute X2 by converting X\n' );
  fprintf ( 1, '  to hypersphere and back.  Consider norm of difference.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  M    || X - X2 ||\n' );
  seed = 123456789;
  n = 1;
  for m = 1 : 5
    fprintf ( 1, '\n' );
    for test = 1 : 5
      [ x, seed ] = r8mat_uniform_01 ( m, n, seed );
      [ c, seed ] = r8vec_uniform_01 ( m, seed );
      [ r, theta ] = cartesian_to_hypersphere ( m, n, c, x );
      x2 = hypersphere_to_cartesian ( m, n, c, r, theta );
      err = norm ( x - x2 );
      fprintf ( 1, '  %d  %g\n', m, err );
    end
  end

  return
end


text:
%%
% Test for discrete 2D sampling.

N = 512;

rep = '../results/sampling/';
if not(exist(rep))
    mkdir(rep);
end

name = 'ring';
name = 'peaks';

% helpers
rho = .8;
t = linspace(0,1,256);
cm = @(c0,c1)(t').^(rho)*c0(:)' + (1-(t').^(rho))*c1(:)';

gaussian = @(m,x,y) exp( -(x-y).^2 );

t = (0:N-1)'/N;
[Y,X] = meshgrid(t,t);
switch name
    case 'peaks'
        a = peaks(N);
        a = abs(a);
        col = [0 0 1];
    case 'ring'
        r = 1*sqrt((X-.5).^2+(Y-.5).^2);
        a = abs(sin(r))^2./(r+eps);
        a = cos(10*r).^2.* exp(-14*r.^2.2);
        a = min(a,.65);
%        surf(a); shading interp;
        col = [0 0 1];
end
a = a/sum(a(:));

% discrete samples
K = 1000;

for K=[100 500 1000 4000 10000]
    I = randsample(N^2,K,true,a(:));
    x = X(I); y = Y(I);

    clf; hold on;
    %imagesc(t,t,a);
    %colormap(cm(col,[1 1 1]));
    ms = 15;
    plot(x,y, '.', 'MarkerSize', ms, 'color', col);
    axis equal;
    axis([0 1 0 1]);  box on;
    set(gca, 'XTick', [], 'YTick', [], 'box', 'on');
    saveas(gcf, [rep name '-disc-' num2str(K) '.eps'], 'epsc');
end

col_cont = [1/2 0 1/2];
clf; hold on;
imagesc(t,t,a);
colormap(cm(col_cont,[1 1 1]));
q = 8;
contour(t,t,a, linspace(0,max(a(:)),q), 'color', col_cont);
caxis([0 max(a(:))]);
% axis off;
set(gca, 'XTick', [], 'YTick', [], 'box', 'on');
axis equal;
saveas(gcf, [rep name '-dens.png'], 'png');

%% pixelize
col_euler = [1 0 0];
b = a/max(a(:));
for n = [8 16 32 64]
P = N/n;
A = zeros(N,N,3);
for i=1:n
    for j=1:n
        seli = (i-1)*P+1:i*P;
        selj = (j-1)*P+1:j*P;
        u = mean(mean(b(seli,selj)));
        for s=1:3
            A(seli,selj,s) = col_euler(s)*u + (1-u);
        end
    end
end
A = A/max(A(:));
imwrite(A, [rep name '-pix-n' num2str(n) '.png'], 'png');
end

author: optimaltransport


text:
%FASTHOUGHTRANSFORM  Calculates 2D Fast Hough transform of an image
%
%     dst = cv.FastHoughTransform(src)
%     dst = cv.FastHoughTransform(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ The source (input) image.
%
% ## Output
% * __dst__ The destination image, result of transformation.
%
% ## Options
% * __DDepth__ The depth of destination image. Default `int32`.
% * __Op__ The operation to be applied. This specifies binary operations, that
%   is such ones which involve two operands. Formally, a binary operation `f`
%   on a set `S` is a binary relation that maps elements of the Cartesian
%   product `SxS` to `S`: `f: SxS -> S`. Default 'Addition'. One of
%   * __Minimum__ Binary minimum operation. The constant specifies the binary
%     minimum operation `f` that is defined as follows: `f(x, y) = min(x, y)`.
%   * __Maximum__ Binary maximum operation. The constant specifies the binary
%     maximum operation `f` that is defined as follows: `f(x, y) = max(x, y)`.
%   * __Addition__ Binary addition operation. The constant specifies the binary
%     addition operation `f` that is defined as follows: `f(x, y) = x + y`.
%   * __Average__ Binary average operation. The constant specifies the binary
%     average operation `f` that is defined as follows: `f(x, y) = (x + y)/2`.
% * __AngleRange__ The part of Hough space to calculate. Each member specifies
%   primarily direction of lines (horizontal or vertical) and the direction of
%   angle changes. Direction of angle changes is from multiples of 90 to odd
%   multiples of 45. The image considered to be written top-down and
%   left-to-right. Angles are started from vertical line and go clockwise.
%   Separate quarters and halves are written in orientation they should be in
%   full Hough space. Default `ARO_315_135`. One of:
%   * **ARO_0_45** Vertical primarily direction and clockwise angle changes.
%   * **ARO_45_90** Horizontal primarily direction and counterclockwise angle
%     changes.
%   * **ARO_90_135** Horizontal primarily direction and clockwise angle
%     changes.
%   * **ARO_315_0** Vertical primarily direction and counterclockwise angle
%     changes.
%   * **ARO_315_45** Vertical primarily direction.
%   * **ARO_45_135** Horizontal primarily direction.
%   * **ARO_315_135** Full set of directions.
%   * **ARO_CTR_HOR** `90 +/- atan(0.5)`, interval approximately from `64.5`
%     to `116.5` degrees. It is used for calculating Fast Hough Transform for
%     images skewed by `atan(0.5)`.
%   * **ARO_CTR_VER** `0 +/- atan(0.5)`, interval approximately from `333.5`
%     (`-26.5`) to `26.5` degrees. It is used for calculating Fast Hough
%     Transform for images skewed by `atan(0.5)`.
% * __MakeSkew__ Specifies to do or not to do skewing of Hough transform
%   image. The enum specifies to do or not to do skewing of Hough transform
%   image so it would be no cycling in Hough transform image through borders
%   of image. Default 'Deskew'. One of:
%   * __Raw__ Use raw cyclic image.
%   * __Deskew__ Prepare deskewed image.
%
% The function calculates the fast Hough transform for full, half or quarter
% range of angles.
%
% See also: cv.HoughPoint2Line, cv.HoughLines, hough, houghlines, houghpeaks
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [xc_incoherent_collapsed_pow xc_incoherent_collapsed_frq n_comb_xc n_comb_sp xc_incoherent_single xc_incoherent sp_incoherent sp] = ...
    xcorr_pss(capbuf,f_search_set,ds_comb_arm,fc,sampling_carrier_twist,k_factor,xc)

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU Affero 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 Affero General Public License for more details.
%
% You should have received a copy of the GNU Affero General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Perform the main correlations necessary to detect the PSS

error(nargchk(7,7,nargin));
error(chk_param(capbuf,'capbuf','vector','horizontal'));
error(chk_param(f_search_set,'f_search_set','vector','real'));
error(chk_param(ds_comb_arm,'ds_comb_arm','scalar','real','integer','>=',0));
error(chk_param(fc,'fc','scalar','real','>',0));

% % Create the time domain pss signals
% persistent pss_td
% persistent pss_td_populated
% if (isempty(pss_td_populated))
%   pss_td=NaN(3,128+9);
%   for t=0:2
%     temp=pss(t);
%     %temp=temp(6:67);
%     temp=[0 temp(32:end) zeros(1,65) temp(1:31)];
%     temp_td=idft(temp)*sqrt(128/62);
%     pss_td(t+1,:)=[temp_td(end-8:end) temp_td];
%   end
%   pss_td_populated=1;
% end

n_cap=length(capbuf);
n_f=length(f_search_set);

% % Correlate against the PSS
% xc=NaN(3,n_cap-136,n_f);
% for foi=1:n_f
% %   f_off=f_search_set(foi);
%   for t=1:3
% %     temp=conj(fshift(pss_td(t,:),f_off,fs_lte/16)/137);
% %     for k=1:n_cap-136
% %       xc(t,k,foi)=sum(temp.*capbuf(k:k+136));
% %     end
%     col_idx = (t-1)*n_f + foi;
%     xc(t,:,foi)=corr_store(1:(n_cap-136),col_idx);
%   end
% end

% Calculate the received power in the vicinity of the possible PSS
sp=NaN(1,n_cap-136-137);
capbufx2=absx2(capbuf);
sp(1)=sum(capbufx2(1:137*2));
for k=2:n_cap-136-137
  sp(k)=sp(k-1)-capbufx2(k-1)+capbufx2(k+137*2-1);
end
sp=sp/(137*2);
assert(any(isnan(sp))==0);

% Perform incoherent combining
xc_incoherent=NaN(3,9600,n_f);
xc_incoherent_single=NaN(3,9600,n_f);
n_comb_xc=floor((size(xc,2)-100)/9600);
for foi=1:n_f

  if sampling_carrier_twist == 1
      f_off=f_search_set(foi);
      % fc*k_factor is the receiver's actual RX center frequency.
      k_factor=(fc-f_off)/fc;
%   else
%       k_factor = 1; % because it is already corrected outside
  end

  for t=1:3
    %xc_incoherent_single(t,:,foi)=sum(transpose(reshape(absx2(xc(t,1:n_comb_xc*9600,foi)),9600,n_comb_xc)),1)/n_comb_xc;
    %xc_incoherent_single(t,:,foi)=absx2(xc(t,1:9600,foi));
    xc_incoherent_single(t,:,foi)=zeros(1,9600);
    for m=1:n_comb_xc
      actual_time_offset=(m-1)*.005*k_factor;
      actual_start_index=round(actual_time_offset*fs_lte/16)+1;
%       xc_incoherent_single(t,:,foi)=xc_incoherent_single(t,:,foi)+absx2(xc(t,actual_start_index:actual_start_index+9599,foi));
      xc_incoherent_single(t,:,foi)=xc_incoherent_single(t,:,foi)+xc(t,actual_start_index:actual_start_index+9599,foi);
    end
  end
  xc_incoherent_single(:,:,foi)=xc_incoherent_single(:,:,foi)/n_comb_xc;
  % Combine adjacent samples that might be different taps from the same channel.
  xc_incoherent(:,:,foi)=xc_incoherent_single(:,:,foi);
  for t=1:ds_comb_arm
    for k=1:3
      xc_incoherent(k,:,foi)=xc_incoherent(k,:,foi)+tshift(xc_incoherent_single(k,:,foi),-t)+tshift(xc_incoherent_single(k,:,foi),t);
    end
  end
  xc_incoherent(:,:,foi)=xc_incoherent(:,:,foi)/(2*ds_comb_arm+1);
end
sp_incoherent=NaN(1,9600);
n_comb_sp=floor(length(sp)/9600);
sp_incoherent=sum(transpose(reshape(sp(1:n_comb_sp*9600),9600,n_comb_sp)),1)/n_comb_sp;

% Align the correlations and the signal power measurements.
% xc(:,1) represents the correlation from capture buffer offset 1 to 137.
% sp(1) represents the signal power for samples 1 through 274.
% Suppose capbuf(1) contains the first sample of an SSS and thus the PSS is
% located from samples 138 to 274. The pss correlation peak will be located
% at time sample 138. This PSS correlation peak must be compared against
% the signal power for samples 1 through 274. Thus, the signal power must
% be cyclically shifted right by 137 samples.
sp_incoherent=tshift(sp_incoherent,137);

% For each time offset and PSS index, find the peak correlation
% among all the frequencies.
xc_incoherent_collapsed_pow=NaN(3,9600);
xc_incoherent_collapsed_frq=NaN(3,9600);
for t=1:3
  [pow frq]=max(transpose(shiftdim(xc_incoherent(t,:,:),1)),[],1);
  xc_incoherent_collapsed_pow(t,:)=pow;
  xc_incoherent_collapsed_frq(t,:)=frq;
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function plot_Bell_Contraction()

offset = 0.125;
tP = 0.25;
total_time = 3*tP + 2*offset;
dt = tP/125;

%Coefficients for Polynomial Phase-Interpolation
a = 2.739726027397260;  % y1(t) = at^2
b = 2.739726027397260;  % y3(t) = -b(t-1)^2+1
c = -2.029426686960933; % y2(t) = ct^3 + dt^2 + gt + h
d = 3.044140030441400;
g = -0.015220700152207;
h = 0.000253678335870;

t1 = 0.05*tP;
t2 = 0.95*tP;
tt1 = 0:dt:tP;
for i=1:length(tt1)
    t = tt1(i);
    if (t <t1)
        r1(i) = a*(t/tP)^2;                             %a*pow((t/tP),2);
    elseif ( (t >=t1 ) && ( t<t2) )
        r1(i) = c*(t/tP)^3 + d*(t/tP)^2 + g*(t/tP) + h; %c*pow((t/tP),3) + d*pow((t/tP),2) + g*(t/tP) + h;
    elseif (t>=t2)
        r1(i) = -b*( (t/tP) - 1.0 )^2 + 1;                  %-b*pow( ((t/tP) - 1),2) + 1;
    end
end

tt2 = tP:dt:tP+2*offset;
tt4 = (2*tP+2*offset):dt:(3*tP+2*offset);
r2 = ones(1,length(tt2));
r4 = zeros(1,length(tt4));
%%RHS Expansion Phase
tt3 = (tP+2*offset):dt:(2*tP+2*offset);
r3 = r1 - 0.5;
r3 = -r3 + 0.5;

%plot(tt1,r1,'r.','MarkerSize',6); hold on;
%plot(tt2,r2,'r.','MarkerSize',6); hold on;
%plot(tt3,r3,'r.','MarkerSize',6); hold on;
%plot(tt4,r4,'r.','MarkerSize',6); hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LHS
%
%%%%%%%%%%%%%%%%%%%%%%%%%

tprev = offset;
t1 = 0.05*tP+offset;
t2 = 0.95*tP+offset;
frac= 0.9;
ttt1 = offset:dt:(tP+offset);
for i=1:length(ttt1)
    t = ttt1(i);
    if (t <t1)
        Lgeo(i) =  a*( (t-tprev) /tP)^2;                             %a*pow((t/tP),2);
    elseif ( (t >=t1 ) && ( t<t2) )
        Lgeo(i) = c*( (t-tprev)/tP)^3 + d*( (t-tprev)/tP)^2 + g*( (t-tprev)/tP) + h; %c*pow((t/tP),3) + d*pow((t/tP),2) + g*(t/tP) + h;
    elseif (t>=t2)
        Lgeo(i) = -b*( ((t-tprev)/tP) - 1.0 )^2 + 1;                  %-b*pow( ((t/tP) - 1),2) + 1;
    end
end

tt0 = 0:dt:offset;
L0 = zeros(1,length(tt0));

%Contraction
L1 = frac*Lgeo;

%Phase Two -> Expansion
ttt2 = (tP+offset):dt:(2*tP+offset);
L2 = Lgeo - 0.5;
L2 = -1.2*L2 + 0.4 - (1-frac);

%Phase Three -> Contraction to equilibrium
ttt3 = (2*tP+offset):dt:(3*tP+offset);
L3 = Lgeo - 1.0;
L3 = (0.3)*L3;

YAXISy = -0.75:dt:1.6;
YAXISx = zeros(1,length(YAXISy));

XAXISx = -0.25:dt:1.5;
XAXISy = zeros(1,length(XAXISx));

plot(tt1,r1,'r.','MarkerSize',8); hold on;
plot(tt0,L0,'.','MarkerSize',8); hold on;

plot(XAXISx,XAXISy,'k-','LineWidth',5);
plot(YAXISx,YAXISy,'k-','LineWidth',5);

plot(tt1,r1,'r.','MarkerSize',8); hold on;
plot(tt0,L0,'.','MarkerSize',8); hold on;

plot(tt2,r2,'r.','MarkerSize',8); hold on;
plot(tt3,r3,'r.','MarkerSize',8); hold on;
plot(tt4,r4,'r.','MarkerSize',8); hold on;

plot(ttt1,L1,'.','MarkerSize',8); hold on;
plot(ttt2,L2,'.','MarkerSize',8); hold on;
plot(ttt3,L3,'.','MarkerSize',8); hold on;
axis([-offset 3*tP+3*offset -0.75 1.5]);

vL1y = -0.75:0.25:0;
vL1x = offset*ones(1,length(vL1y));

vL2y = -0.75:0.25:0.9;
vL2x = (tP+offset)*ones(1,length(vL2y));

vL3y = -0.75:0.25:-0.21;
vL3x = (2*tP+offset)*ones(1,length(vL3y));

vL4y = -0.75:0.25:0.0;
vL4x = (3*tP+offset)*ones(1,length(vL4y));

vL5y = -0.75:0.25:0.0;
vL5x = (3*tP+2*offset)*ones(1,length(vL4y));

vR1y = 1:0.2:1.6;
vR1x = tP*ones(1,length(vR1y));

vR2y = 1:0.2:1.6;
vR2x = (tP+2*offset)*ones(1,length(vR2y));

vR3y = 0:0.2:1.6;
vR3x = (2*tP+2*offset)*ones(1,length(vR3y));

plot(vR1x,vR1y,'k--','LineWidth',1); hold on;
plot(vR2x,vR2y,'k--','LineWidth',1); hold on;
plot(vR3x,vR3y,'k--','LineWidth',1); hold on;

plot(vL1x,vL1y,'k--','LineWidth',1); hold on;
plot(vL2x,vL2y,'k--','LineWidth',1); hold on;
plot(vL3x,vL3y,'k--','LineWidth',1); hold on;
plot(vL4x,vL4y,'k--','LineWidth',1); hold on;
plot(vL5x,vL5y,'k--','LineWidth',1); hold on;
axis([-2*offset 3*tP+4*offset -1.0 1.75]);

legend('RIGHT','LEFT');

author: nickabattista


text:
function gaussian_test ( )

%*****************************************************************************80
%
%% GAUSSIAN_TEST examines the gaussian correlation.
%
%  Discussion:
%
%    Note that I can't use the original value EIGEN_NUM = 20.  The
%    EIGS function fails to converge.  In fact, it still fails at 12.
%
%    This code is based substantially on a document by Toby Driscoll.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Toby Driscoll,
%    Mercer's theorem and the Karhunen-Loeve expansion,
%    Oxford University Mathematical Institute,
%    http://www2.maths.ox.ac.uk/chebfun/examples/stats/pdf/MercerKarhunenLoeve.pdf
%
  rmpath ( '/usr/local/matlab/toolbox/datafeed/datafeed' )
  addpath ( '../chebfun' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSSIAN_TEST:\n' );
  fprintf ( 1, '  Demonstrate Mercer''s theorem and the KL expansion\n' );
  fprintf ( 1, '  for the gaussian kernel.\n' );
%
%  Set the interval.
%
  a = 0.0;
  b = 10.0;
  fprintf ( 1, '  Using interval [%g,%g]\n', a, b );
  s_num = 21;
  s_vec = linspace ( a, b, s_num );
%
%  FRED is a function from the CHEBFUN library.
%
%  It constructs a "chebop" representing the Fredholm integral operator
%  with kernel K for functions in domain [A,B].
%
  F = fred ( @gaussian_correlation, domain ( [ a, b ] ) );
%
%  EIGS has been extended to be able to compute the eigenvalues Lambda
%  and eigenfunctions Psi of the Fredholm integral operator represented by F.
%
%  The "LM" switch requests that we return the eigenvalues of largest magnitude.
%
%  Each Psi is a CHEBFUN, that is, it takes a real number argument.
%
  eigen_num = 20;
  [ Psi, Lambda ] = eigs ( F, eigen_num, 'lm' );

  eigen_found = length ( Lambda );
  fprintf ( 1, '  Requested %d eigenmodes, computed %d\n', eigen_num, eigen_found );
  eigen_num = min ( eigen_num, eigen_found );
%
%  Print the eigenvalues.
%
  lambda_vec = diag ( Lambda );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I       Lambda(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %2d  %10g\n', i, lambda_vec(i) );
  end
%
%  Plot the eigenvalues.
%
  figure ( 1 )
  clf
  plot ( lambda_vec, 'Linewidth', 2 );
  hold on
  plot ( lambda_vec, 'b.', 'Markersize', 20 );
  title ( 'gaussian: Mercer eigenvalues' );
  xlabel ( '<--- N --->')
  grid on
  print -dpng 'gaussian_figure1.png'
%
%  Plot selected eigenfunctions.
%
  figure ( 2 )
  subplot ( 4, 1, 1 )
  plot ( Psi(:,1), 'Linewidth', 2 );
  title ( 'gaussian: Mercer eigenfunction PSI(1)' )
  grid on
  subplot ( 4, 1, 2 )
  plot ( Psi(:,2), 'Linewidth', 2 );
  title ( 'gaussian: Mercer eigenfunction PSI(2)' )
  grid on
  subplot ( 4, 1, 3 )
  plot ( Psi(:,5), 'Linewidth', 2 );
  title ( 'gaussian: Mercer eigenfunction PSI(5)' )
  grid on
  subplot ( 4, 1, 4 )
  plot ( Psi(:,10), 'Linewidth', 2 );
  title ( 'gaussian: Mercer eigenfunction PSI(10)' )
  grid on
  print -dpng 'gaussian_figure2.png'
%
%  Orthonormality check.
%
  ptp = Psi' * Psi;
  error_frobenius = r8mat_is_identity ( eigen_num, ptp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of I - Psi'' * Psi = %g\n', error_frobenius );
%
%  K(S,S) should be exactly 1.
%  Because we are using a truncated representation of K, our estimate of K(S,S)
%  will be smaller than 1.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Truncated estimate of K(s,s) = 1 for S in the interval.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        S        K(s,s) estimate\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 21
    s = s_vec(i);
    ptlp = Psi(s,:) * Lambda * Psi(s,:)';
    fprintf ( 1, '  %10g  %14g\n', s, ptlp );
  end
%
%  Look at eigenvalue decay.
%
  x = diff ( log ( ( 1:eigen_num ) ) );
  y = diff ( log ( lambda_vec ) )';
  c = y ./ x;
  figure ( 3 );
  clf
  plot ( c, 'Linewidth', 2 );
  grid on
  hold on
  plot ( c, 'b.', 'Markersize', 20 );
  xlabel ( '<-- N -->' )
  title ( 'gaussian: Eigenvalue decay rate' );
  print -dpng 'gaussian_figure3.png'
%
%  Look at eigenvalue sum.
%
%  The trace of K(s,s) over [a,b] should be the integral of 1 over [a,b],
%  that is, b - a.
%
%  We compare the trace to the partial sums of the lambda's to see how much
%  of the variance of the process we have captured.
%
  trace_K = b - a;
  lambda_cum = cumsum ( lambda_vec ) / trace_K;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index    Cumulative Eigenvalue sum\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %5d  %10g\n', i, lambda_cum(i) );
  end
%
%  We decide to use the first 10 eigenfunctions.
%
  eigen_use = 10;
%
%  Find 400 realizations of the process by selecting, for each realization,
%  10 random parameters Z in the truncated KL expansion.
%
  Z = randn ( eigen_use, 400 );
  X = Psi(:,1:eigen_use) * ( sqrt ( Lambda(1:eigen_use,1:eigen_use) ) * Z );
%
%  Plot 40 of the realizations;
%  Plot their mean, computed from all 400.
%
  figure ( 4 )
  clf
  plot ( X(:,1:40) )
  mu = sum ( X, 2 ) / 400;
  hold on;
  plot ( mu, 'k', 'linewidth', 3 );
  title ( 'gaussian: 40 Random Realizations X(t,omega), and their Mean.' )
  print -dpng 'gaussian_figure4.png'
%
%  Estimate the covariance from the data.
%
  figure ( 5 )
  clf
  [ S, T ] = meshgrid ( s_vec, s_vec );
  C = cov ( X(s_vec,:)' );
  mesh ( S, T, C );
  hold on;
  D = gaussian_correlation ( S, T );
  plot3 ( S, T, D, 'k.', 'markersize', 10 )
  title ( 'gaussian: Covariance K(S,T) (dots), and Estimate from Realizations (mesh)' )
  print -dpng 'gaussian_figure5.png'
%
%  Using just 10 functions in the expansion,
%  reduce the correlation length,
%  and examine the sum of the lambda's.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use a fixed number of eigenfunctions = %d\n', eigen_use );
  fprintf ( 1, '  but vary the correlation length RHOBAR.\n' );
  fprintf ( 1, '  (We used RHOBAR = 1 above.)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The sum of the eigenvalues, divided by (B-A),\n' );
  fprintf ( 1, '  discloses the relative amount of the total variation\n' );
  fprintf ( 1, '  that is captured by the truncated expansion.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      RHOBAR      VARSUM\n' );
  fprintf ( 1, '\n' );
  rhobar = 4.0;
  for i = 1 : 10
    K = @ ( s, t ) exp ( - ( ( s - t ) / rhobar ).^2 );
    F = fred ( K, domain ( [ a, b ] ) );
    lambda_vec = eigs ( F, eigen_use, 'lm' );
    varsum = sum ( lambda_vec(1:eigen_use) ) / ( b - a );
    fprintf ( 1, '  %10g  %10g\n', rhobar, varsum );
    rhobar = rhobar / 2.0;
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAUSSIAN_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../chebfun' )

  return
end


text:
% Compile sharedmatrix mex file.
%
% Copyright (c) 2010,2011 Joshua V Dillon
% All rights reserved. (See file header for details.)

fprintf('Compiling sharedmatrix mex files...\n');

numbits=computer;numbits=str2num(numbits(end-1:end))

if numbits==64
	% For 64-bit machines:
	mex -largeArrayDims -O -v sharedmatrix.c
else
	% For 32-bit machines:
	mex -O -v sharedmatrix.c
end

fprintf('Done.\n');

% --- Miscellaneous Comments ------------------------------------------------

% Valgrind doesn't seem to work here.
% matlab -nodesktop -nojvm -nosplash -r "x=rand(3,4);sharedmatrix('clone',12345,x);quit" -D"valgrind --error-limit=no --tool=memcheck -v --log-file=valgrind.log --leak-check=full --track-origins=yes"

% mex -largeArrayDims -O -v CFLAGS='$CFLAGS -Wall -Werror -ansi -DDEBUG -D_XOPEN_SOURCE -fPIC' sharedmatrix.c

% --- Possibly useful notes. ------------------------------------------------

% --- Note 1 ---

% If you have trouble compiling in Matlab (GLIB error) try prefixing the matlab
% command with a LD_PRELOAD, ie,
% $ LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6 matlab

% --- Note 2 ---

% http://morganbye.net/blog/2011/05/matlab-ubuntu-1104
%
% MATLAB with Ubuntu 11.04
%
% Under the new Ubuntu 11.04 (Natty Narwhal), when you try and run MATLAB
% you'll get the following error
%
%     /matlab/bin/util/oscheck.sh: 605: /lib64/libc.so.6: not found
%
% To resolve this problem, go to a terminal window (Alt + F2 from anywhere) and
% type:
%
% For 64 bit:
%     sudo ln -s /lib64/x86_64-linux-gnu/libc-2.13.so /lib64/libc.so.6
%
% For 32 bit:
%     sudo ln -s /lib/i386-linux-gnu/libc-2.13.so /lib/libc.so.6
%
% This should restore the missing library that was uninstalled during the
% Ubuntu update process.

author: Sable

repo: mcbench-benchmarks


text:
function geometry_test015 ( )

%*****************************************************************************80
%
%% TEST015 tests CIRCLE_EXP2IMP_2D, TRIANGLE_CIRCUMCIRCLE_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST015\n' );
  fprintf ( 1, '  CIRCLE_EXP2IMP_2D computes the radius and \n' );
  fprintf ( 1, '  center of the circle through three points.\n' );
  fprintf ( 1, '  TRIANGLE_CIRCUMCIRCLE_2D computes the radius and \n' );
  fprintf ( 1, '  center of the circle through the vertices of\n' );
  fprintf ( 1, '  a triangle.\n' );

  p1test(1:2,1:n) = [ ...
    4.0, 2.0; ...
    1.0, 5.0; ...
   -2.0, 2.0 ]';

  p2test(1:2,1:n) = [ ...
    4.0, 2.0; ...
    5.0, 4.0; ...
    6.0, 6.0 ]';

  p3test(1:2,1:n) = [ ...
    4.0, 2.0; ...
    1.0, 5.0; ...
    4.0, 2.0 ]';

  for i = 1 : n

    p1(1:2,1) = p1test(1:2,i);
    p2(1:2,1) = p2test(1:2,i);
    p3(1:2,1) = p3test(1:2,i);

    r8vec_print ( 2, p1, '  P1:' );
    r8vec_print ( 2, p2, '  P2:' );
    r8vec_print ( 2, p3, '  P3:' );

    [ r, center ] = circle_exp2imp_2d ( p1, p2, p3 );

    circle_imp_print_2d ( r, center, '  The implicit circle:' )

    t(1:2,1:3) = [ p1(1:2,1)'; p2(1:2,1)'; p3(1:2,1)' ]';

    [ r, center ] = triangle_circumcircle_2d ( t );

    circle_imp_print_2d ( r, center, '  The triangle''s circumcircle:' )

  end

  return
end


text:
function [depVarValue, depVarUnit, taskStr, refBodyInfo] = lvd_getDepVarValueUnit(i, subLog, taskStr, refBodyId, celBodyData, onlyReturnTaskStr, inFrame)
    %lvd_getDepVarValueUnit Summary of this function goes here
    %   Detailed explanation goes here

    if(~isempty(refBodyId))
        refBodyInfo = getBodyInfoByNumber(refBodyId, celBodyData);
    else
        refBodyInfo = [];
    end

    if(onlyReturnTaskStr == true)
        depVarValue = NaN;
        depVarUnit = NaN;

        return;
    end

    switch taskStr
        case 'Yaw Angle'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'yaw', inFrame);
            depVarUnit = 'deg';
        case 'Pitch Angle'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'pitch', inFrame);
            depVarUnit = 'deg';
        case 'Roll Angle'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'roll', inFrame);
            depVarUnit = 'deg';
        case 'Bank Angle'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'bank', inFrame);
            depVarUnit = 'deg';
        case 'Angle of Attack'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'angleOfAttack', inFrame);
            depVarUnit = 'deg';
        case 'SideSlip Angle'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'sideslip', inFrame);
            depVarUnit = 'deg';
        case 'Total Angle of Attack'
            depVarValue = lvd_SteeringAngleTask(subLog(i), 'totalAngleOfAttack', inFrame);
            depVarUnit = 'deg';
        case 'Throttle'
            depVarValue = lvd_ThrottleTask(subLog(i), 'throttle', inFrame);
            depVarUnit = 'Percent';
        case 'Thrust to Weight Ratio'
            depVarValue = lvd_ThrottleTask(subLog(i), 't2w', inFrame);
            depVarUnit = ' ';
        case 'Total Thrust'
            depVarValue = lvd_ThrottleTask(subLog(i), 'totalthrust', inFrame);
            depVarUnit = 'kN';
        case 'Two-Body Time To Impact'
            depVarValue = lvd_TwoBodyImpactPointTasks(subLog(i), 'timeToImpact');
            depVarUnit = 'sec';
        case 'Two-Body Impact Latitude'
            depVarValue = lvd_TwoBodyImpactPointTasks(subLog(i), 'latitude');
            depVarUnit = 'degN';
        case 'Two-Body Impact Longitude'
            depVarValue = lvd_TwoBodyImpactPointTasks(subLog(i), 'longitude');
            depVarUnit = 'degE';
        case 'Drag Coefficient'
            depVarValue = lvd_AeroTasks(subLog(i), 'dragCoeff', []);
            depVarUnit = 'm^2';
        case 'Event Number'
            depVarValue = lvd_EventNumTask(subLog(i), 'eventNum');
            depVarUnit = '';
        case 'Total Effective Isp'
            depVarValue = lvd_PropulsionTasks(subLog(i), 'totalEffIsp');
            depVarUnit = 'sec';
        case 'Total Thrust Vector X Component'
            depVarValue = lvd_ThrottleTask(subLog(i), 'thrust_x', inFrame);
            depVarUnit = 'kN';
        case 'Total Thrust Vector Y Component'
            depVarValue = lvd_ThrottleTask(subLog(i), 'thrust_y', inFrame);
            depVarUnit = 'kN';
        case 'Total Thrust Vector Z Component'
            depVarValue = lvd_ThrottleTask(subLog(i), 'thrust_z', inFrame);
            depVarUnit = 'kN';
        case 'Power Net Charge Rate'
            [depVarValue, depVarUnit] = lvd_ElectricalPowerGlobalTasks(subLog(i), 'netChargeRate');
        case 'Power Cumulative Storage State of Charge'
            [depVarValue, depVarUnit] = lvd_ElectricalPowerGlobalTasks(subLog(i), 'cumStorageSoC');
        case 'Power Maximum Available Storage'
            [depVarValue, depVarUnit] = lvd_ElectricalPowerGlobalTasks(subLog(i), 'maxAvailableStorage');

        case 'Position Vector (X)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'rVectX', inFrame);
            depVarUnit = 'km';
        case 'Position Vector (Y)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'rVectY', inFrame);
            depVarUnit = 'km';
        case 'Position Vector (Z)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'rVectZ', inFrame);
            depVarUnit = 'km';
        case 'Velocity Vector (X)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'vVectX', inFrame);
            depVarUnit = 'km/s';
        case 'Velocity Vector (Y)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'vVectY', inFrame);
            depVarUnit = 'km/s';
        case 'Velocity Vector (Z)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'vVectZ', inFrame);
            depVarUnit = 'km/s';
        case 'Position Vector Magnitude'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'radius', inFrame);
            depVarUnit = 'km';
        case 'Velocity Vector Magnitude'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'velocity', inFrame);
            depVarUnit = 'km/s';
        case 'Surface Velocity'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'horzVel', inFrame);
            depVarUnit = 'km/s';
        case 'Vertical Velocity'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'vertVel', inFrame);
            depVarUnit = 'km/s';

        case 'Semi-major Axis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'sma', inFrame);
            depVarUnit = 'km';
        case 'Eccentricity'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'ecc', inFrame);
            depVarUnit = '';
        case 'Inclination'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'inc', inFrame);
            depVarUnit = 'deg';
        case 'Right Asc. of the Asc. Node'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'raan', inFrame);
            depVarUnit = 'deg';
        case 'Argument of Periapsis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'arg', inFrame);
            depVarUnit = 'deg';
        case 'True Anomaly'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'tru', inFrame);
            depVarUnit = 'deg';
        case 'Mean Anomaly'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'mean', inFrame);
            depVarUnit = 'deg';
        case 'Orbital Period'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'period', inFrame);
            depVarUnit = 'sec';
        case 'Radius of Periapsis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'rPe', inFrame);
            depVarUnit = 'km';
        case 'Radius of Apoapsis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'rApo', inFrame);
            depVarUnit = 'km';
        case 'Altitude of Apoapsis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'altApo', inFrame);
            depVarUnit = 'km';
        case 'Altitude of Periapsis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'altPeri', inFrame);
            depVarUnit = 'km';
        case 'Equinoctial H1'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'H1', inFrame);
            depVarUnit = '';
        case 'Equinoctial K1'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'K1', inFrame);
            depVarUnit = '';
        case 'Equinoctial H2'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'H2', inFrame);
            depVarUnit = '';
        case 'Equinoctial K2'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'K2', inFrame);
            depVarUnit = '';
        case 'Flight Path Angle'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'FPA', inFrame);
            depVarUnit = 'deg';

        case 'Altitude'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'alt', inFrame);
            depVarUnit = 'km';
        case 'Longitude (East)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'long', inFrame);
            depVarUnit = 'deg';
        case 'Latitude (North)'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'lat', inFrame);
            depVarUnit = 'deg';
        case 'Velocity Azimuth'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'velAz', inFrame);
            depVarUnit = 'deg';
        case 'Velocity Elevation'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'velEl', inFrame);
            depVarUnit = 'deg';
        case 'Longitudinal Drift Rate'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'longDriftRate', inFrame);
            depVarUnit = 'deg/hr';

        case 'C3 Energy'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'c3', inFrame);
            depVarUnit = 'km^2/s^2';
        case 'Seconds Past Periapsis'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'tau', inFrame);
            depVarUnit = 'sec';

        case 'Central Body ID'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'centralBodyId', inFrame);
            depVarUnit = '';

        case  'Hyperbolic Velocity Unit Vector X'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'hyperVelUnitVectX', inFrame);
            depVarUnit = '';
        case  'Hyperbolic Velocity Unit Vector Y'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'hyperVelUnitVectY', inFrame);
            depVarUnit = '';
        case  'Hyperbolic Velocity Unit Vector Z'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'hyperVelUnitVectZ', inFrame);
            depVarUnit = '';
        case  'Hyperbolic Velocity Vector Right Ascension'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'hyperVelUnitVectRA', inFrame);
            depVarUnit = 'deg';
        case  'Hyperbolic Velocity Vector Declination'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'hyperVelUnitVectDec', inFrame);
            depVarUnit = 'deg';
        case  'Hyperbolic Velocity Magnitude'
            depVarValue = lvd_KinematicStateTasks(subLog(i), 'hyperVelMag', inFrame);
            depVarUnit = 'km/s';
        case 'Drag Force'
            depVarValue = lvd_AeroTasks(subLog(i), 'dragForce', inFrame);
            depVarUnit = 'kN';

        otherwise %is a programmatically generated string that we'll handle here
            fluidTypeMassPattern = '^Fluid Type (\d+?) Mass - ".*"';

            tankMassPattern = '^Tank (\d+?) Mass - ".*"';
            tankMassDotPattern = '^Tank (\d+?) Mass Flow Rate - ".*"';

            stageDryMassPattern = '^Stage (\d+?) Dry Mass - ".*"';
            stageActivePattern = '^Stage (\d+?) Active State - ".*"';

            engineActivePattern = '^Engine (\d+?) Active State - ".*"';

            stopwatchValuePattern = '^Stopwatch (\d+?) Value - ".*"';

            extremaValuePattern = '^Extrema (\d+?) Value - ".*"';

            grdObjAzValuePattern = '^Ground Object (\d+?) Azimuth to S/C - ".*"';
            grdObjElValuePattern = '^Ground Object (\d+?) Elevation to S/C - ".*"';
            grdObjRngValuePattern = '^Ground Object (\d+?) Range to S/C - ".*"';
            grdObjLoSValuePattern = '^Ground Object (\d+?) Line of Sight to S/C - ".*"';

            calcObjValuePattern = '^Calculus (\d+?) Value - ".*"';

            pwrStorageActivePattern = '^Power Storage (\d+?) Active State - ".*"';
            pwrStorageSoCPattern = '^Power Storage (\d+?) State of Charge - ".*"';
            pwrSinkActivePattern = '^Power Sink (\d+?) Active State - ".*"';
            pwrSinkDischargeRatePattern = '^Power Sink (\d+?) Discharge Rate - ".*"';
            pwrSrcActivePattern = '^Power Source (\d+?) Active State - ".*"';
            pwrSrcChargeRatePattern = '^Power Source (\d+?) Charge Rate - ".*"';

            geoPtPosXPattern = '^Point (\d+?) Position \(X\) - ".*"';
            geoPtPosYPattern = '^Point (\d+?) Position \(Y\) - ".*"';
            geoPtPosZPattern = '^Point (\d+?) Position \(Z\) - ".*"';
            geoPtVelXPattern = '^Point (\d+?) Velocity \(X\) - ".*"';
            geoPtVelYPattern = '^Point (\d+?) Velocity \(Y\) - ".*"';
            geoPtVelZPattern = '^Point (\d+?) Velocity \(Z\) - ".*"';

            geoVectXPattern = '^Vector (\d+?) X Component - ".*"';
            geoVectYPattern = '^Vector (\d+?) Y Component - ".*"';
            geoVectZPattern = '^Vector (\d+?) Z Component - ".*"';
            geoVectMagPattern = '^Vector (\d+?) Magnitude - ".*"';
            geoVectOriginXPattern = '^Vector (\d+?) Origin Position \(X\) - ".*"';
            geoVectOriginYPattern = '^Vector (\d+?) Origin Position \(Y\) - ".*"';
            geoVectOriginZPattern = '^Vector (\d+?) Origin Position \(Z\) - ".*"';

            geoAngleMagPattern = '^Angle (\d+?) Magnitude - ".*"';

            geoPlaneNormXPattern = '^Plane (\d+?) Normal Vector X Component - ".*"';
            geoPlaneNormYPattern = '^Plane (\d+?) Normal Vector Y Component - ".*"';
            geoPlaneNormZPattern = '^Plane (\d+?) Normal Vector Z Component - ".*"';

            pluginValuePattern = '^Plugin (\d+?) Value - ".*"';

            if(not(isempty(regexpi(taskStr, fluidTypeMassPattern)))) %#ok<*RGXPI>
                tokens = regexpi(taskStr, fluidTypeMassPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                fluidTypeInd = str2double(tokens);

                [~,fluidTypes] = subLog(i).launchVehicle.tankTypes.getFluidTypesGraphAnalysisTaskStrs();
                fluidType = fluidTypes(fluidTypeInd);

                depVarValue = lvd_FluidTypeMassTasks(subLog(i), 'totalFluidTypeMass', fluidType);
                depVarUnit = 'mT';

            elseif(not(isempty(regexpi(taskStr, tankMassPattern))))
                tokens = regexpi(taskStr, tankMassPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                tankInd = str2double(tokens);

                [~,tanks] = subLog(i).launchVehicle.getTanksGraphAnalysisTaskStrs();
                tank = tanks(tankInd);

                depVarValue = lvd_TankMassTasks(subLog(i), 'tankMass', tank);
                depVarUnit = 'mT';

            elseif(not(isempty(regexpi(taskStr, tankMassDotPattern))))
                tokens = regexpi(taskStr, tankMassDotPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                tankInd = str2double(tokens);

                [~,tanks] = subLog(i).launchVehicle.getTanksGraphAnalysisTaskStrs();
                tank = tanks(tankInd);

                depVarValue = lvd_TankMassTasks(subLog(i), 'tankMDot', tank);
                depVarUnit = 'mT/s';

            elseif(not(isempty(regexpi(taskStr, stageDryMassPattern))))
                tokens = regexpi(taskStr, stageDryMassPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                stageInd = str2double(tokens);

                [~,stages] = subLog(i).launchVehicle.getStagesGraphAnalysisTaskStrs();
                stage = stages(stageInd);

                depVarValue = lvd_StageTasks(subLog(i), 'dryMass', stage);
                depVarUnit = 'mT';

            elseif(not(isempty(regexpi(taskStr, stageActivePattern))))
                tokens = regexpi(taskStr, stageActivePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                stageInd = str2double(tokens);

                [~,stages] = subLog(i).launchVehicle.getStagesGraphAnalysisTaskStrs();
                stage = stages(stageInd);

                depVarValue = lvd_StageTasks(subLog(i), 'active', stage);
                depVarUnit = '';

            elseif(not(isempty(regexpi(taskStr, engineActivePattern))))
                tokens = regexpi(taskStr, engineActivePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                engineInd = str2double(tokens);

                [~,engines] = subLog(i).launchVehicle.getEnginesGraphAnalysisTaskStrs();
                engine = engines(engineInd);

                depVarValue = lvd_EngineTasks(subLog(i), 'active', engine);
                depVarUnit = '';

            elseif(not(isempty(regexpi(taskStr, stopwatchValuePattern))))
                tokens = regexpi(taskStr, stopwatchValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                stopwatchInd = str2double(tokens);

                [~,stopwatches] = subLog(i).launchVehicle.getStopwatchGraphAnalysisTaskStrs();
                stopwatch = stopwatches(stopwatchInd);

                depVarValue = lvd_StopwatchMassTasks(subLog(i), 'swValue', stopwatch);
                depVarUnit = 'sec';

            elseif(not(isempty(regexpi(taskStr, extremaValuePattern))))
                tokens = regexpi(taskStr, extremaValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                calcObjInd = str2double(tokens);

                [~,calcObjs] = subLog(i).launchVehicle.getExtremaGraphAnalysisTaskStrs();
                calcObj = calcObjs(calcObjInd);

                [depVarValue, depVarUnit] = lvd_ExtremaTasks(subLog(i), 'extremumValue', calcObj);

            elseif(not(isempty(regexpi(taskStr, grdObjAzValuePattern))))
                tokens = regexpi(taskStr, grdObjAzValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                ind = str2double(tokens);

                [~, grdObjs] = subLog(i).lvdData.groundObjs.getGrdObjAzGraphAnalysisTaskStrs();
                grdObj = grdObjs(ind);

                [depVarValue, depVarUnit] = lvd_GrdObjTasks(subLog(i), 'azimuth', grdObj, inFrame);

            elseif(not(isempty(regexpi(taskStr, grdObjElValuePattern))))
                tokens = regexpi(taskStr, grdObjElValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                ind = str2double(tokens);

                [~, grdObjs] = subLog(i).lvdData.groundObjs.getGrdObjAzGraphAnalysisTaskStrs();
                grdObj = grdObjs(ind);

                [depVarValue, depVarUnit] = lvd_GrdObjTasks(subLog(i), 'elevation', grdObj, inFrame);

            elseif(not(isempty(regexpi(taskStr, grdObjRngValuePattern))))
                tokens = regexpi(taskStr, grdObjRngValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                ind = str2double(tokens);

                [~, grdObjs] = subLog(i).lvdData.groundObjs.getGrdObjAzGraphAnalysisTaskStrs();
                grdObj = grdObjs(ind);

                [depVarValue, depVarUnit] = lvd_GrdObjTasks(subLog(i), 'range', grdObj, inFrame);

            elseif(not(isempty(regexpi(taskStr, grdObjLoSValuePattern))))
                tokens = regexpi(taskStr, grdObjLoSValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                ind = str2double(tokens);

                [~, grdObjs] = subLog(i).lvdData.groundObjs.getGrdObjAzGraphAnalysisTaskStrs();
                grdObj = grdObjs(ind);

                [depVarValue, depVarUnit] = lvd_GrdObjTasks(subLog(i), 'LoS', grdObj, inFrame);

            elseif(not(isempty(regexpi(taskStr, calcObjValuePattern))))
                tokens = regexpi(taskStr, calcObjValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                calcObjInd = str2double(tokens);

                [~,calcObjs] = subLog(i).launchVehicle.getCalculusCalcObjGraphAnalysisTaskStrs();
                calcObj = calcObjs(calcObjInd);

                [depVarValue, depVarUnit] = lvd_CalculusCalculationTasks(subLog(i), 'calcObjValue', calcObj);

            elseif(not(isempty(regexpi(taskStr, pwrStorageActivePattern))))
                tokens = regexpi(taskStr, pwrStorageActivePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pwrStorageObjInd = str2double(tokens);

                [~, pwrStorages] = subLog(i).launchVehicle.getPowerStorageGraphAnalysisTaskStrs();
                pwrStorage = pwrStorages(pwrStorageObjInd);

                [depVarValue, depVarUnit] = lvd_ElectricalPowerStorageTasks(subLog(i), 'active', pwrStorage);

            elseif(not(isempty(regexpi(taskStr, pwrStorageSoCPattern))))
                tokens = regexpi(taskStr, pwrStorageSoCPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pwrStorageObjInd = str2double(tokens);

                [~, pwrStorages] = subLog(i).launchVehicle.getPowerStorageGraphAnalysisTaskStrs();
                pwrStorage = pwrStorages(pwrStorageObjInd);

                [depVarValue, depVarUnit] = lvd_ElectricalPowerStorageTasks(subLog(i), 'stateOfCharge', pwrStorage);

            elseif(not(isempty(regexpi(taskStr, pwrSinkActivePattern))))
                tokens = regexpi(taskStr, pwrSinkActivePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pwrSinkObjInd = str2double(tokens);

                [~, pwrSinks] = subLog(i).launchVehicle.getPowerSinksGraphAnalysisTaskStrs();
                pwrSink = pwrSinks(pwrSinkObjInd);

                [depVarValue, depVarUnit] = lvd_ElectricalPowerSinkTasks(subLog(i), 'active', pwrSink);

            elseif(not(isempty(regexpi(taskStr, pwrSinkDischargeRatePattern))))
                tokens = regexpi(taskStr, pwrSinkDischargeRatePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pwrSinkObjInd = str2double(tokens);

                [~, pwrSinks] = subLog(i).launchVehicle.getPowerSinksGraphAnalysisTaskStrs();
                pwrSink = pwrSinks(pwrSinkObjInd);

                [depVarValue, depVarUnit] = lvd_ElectricalPowerSinkTasks(subLog(i), 'dischargeRate', pwrSink);

            elseif(not(isempty(regexpi(taskStr, pwrSrcActivePattern))))
                tokens = regexpi(taskStr, pwrSrcActivePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pwrSrcObjInd = str2double(tokens);

                [~, powerSrcs] = subLog(i).launchVehicle.getPowerSrcsGraphAnalysisTaskStrs();
                powerSrc = powerSrcs(pwrSrcObjInd);

                [depVarValue, depVarUnit] = lvd_ElectricalPowerSrcTasks(subLog(i), 'active', powerSrc);

            elseif(not(isempty(regexpi(taskStr, pwrSrcChargeRatePattern))))
                tokens = regexpi(taskStr, pwrSrcChargeRatePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pwrSrcObjInd = str2double(tokens);

                [~, powerSrcs] = subLog(i).launchVehicle.getPowerSrcsGraphAnalysisTaskStrs();
                powerSrc = powerSrcs(pwrSrcObjInd);

                [depVarValue, depVarUnit] = lvd_ElectricalPowerSrcTasks(subLog(i), 'chargeRate', powerSrc);

            elseif(not(isempty(regexpi(taskStr, geoPtPosXPattern))))
                tokens = regexpi(taskStr, geoPtPosXPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pointInd = str2double(tokens);

                [~, points] = subLog(i).lvdData.geometry.points.getPointPositionXGraphAnalysisTaskStrs();
                point = points(pointInd);

                [depVarValue, depVarUnit] = lvd_GeometricPointTasks(subLog(i), 'PosX', point, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPtPosYPattern))))
                tokens = regexpi(taskStr, geoPtPosYPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pointInd = str2double(tokens);

                [~, points] = subLog(i).lvdData.geometry.points.getPointPositionXGraphAnalysisTaskStrs();
                point = points(pointInd);

                [depVarValue, depVarUnit] = lvd_GeometricPointTasks(subLog(i), 'PosY', point, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPtPosZPattern))))
                tokens = regexpi(taskStr, geoPtPosZPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pointInd = str2double(tokens);

                [~, points] = subLog(i).lvdData.geometry.points.getPointPositionXGraphAnalysisTaskStrs();
                point = points(pointInd);

                [depVarValue, depVarUnit] = lvd_GeometricPointTasks(subLog(i), 'PosZ', point, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPtVelXPattern))))
                tokens = regexpi(taskStr, geoPtVelXPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pointInd = str2double(tokens);

                [~, points] = subLog(i).lvdData.geometry.points.getPointPositionXGraphAnalysisTaskStrs();
                point = points(pointInd);

                [depVarValue, depVarUnit] = lvd_GeometricPointTasks(subLog(i), 'VelX', point, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPtVelYPattern))))
                tokens = regexpi(taskStr, geoPtVelYPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pointInd = str2double(tokens);

                [~, points] = subLog(i).lvdData.geometry.points.getPointPositionXGraphAnalysisTaskStrs();
                point = points(pointInd);

                [depVarValue, depVarUnit] = lvd_GeometricPointTasks(subLog(i), 'VelY', point, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPtVelZPattern))))
                tokens = regexpi(taskStr, geoPtVelZPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pointInd = str2double(tokens);

                [~, points] = subLog(i).lvdData.geometry.points.getPointPositionXGraphAnalysisTaskStrs();
                point = points(pointInd);

                [depVarValue, depVarUnit] = lvd_GeometricPointTasks(subLog(i), 'VelZ', point, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectXPattern))))
                tokens = regexpi(taskStr, geoVectXPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getVectorXComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'VectorX', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectYPattern))))
                tokens = regexpi(taskStr, geoVectYPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getVectorYComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'VectorY', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectZPattern))))
                tokens = regexpi(taskStr, geoVectZPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getVectorZComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'VectorZ', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectMagPattern))))
                tokens = regexpi(taskStr, geoVectMagPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getVectorMagComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'VectMag', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectOriginXPattern))))
                tokens = regexpi(taskStr, geoVectOriginXPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getOriginPosXComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'OriginX', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectOriginYPattern))))
                tokens = regexpi(taskStr, geoVectOriginYPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getOriginPosYComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'OriginY', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoVectOriginZPattern))))
                tokens = regexpi(taskStr, geoVectOriginZPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                vectInd = str2double(tokens);

                [~, vectors] = subLog(i).lvdData.geometry.vectors.getOriginPosZComponentGraphAnalysisTaskStrs();
                vector = vectors(vectInd);

                [depVarValue, depVarUnit] = lvd_GeometricVectorTasks(subLog(i), 'OriginZ', vector, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoAngleMagPattern))))
                tokens = regexpi(taskStr, geoAngleMagPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                angleInd = str2double(tokens);

                [~, angles] = subLog(i).lvdData.geometry.angles.getAngleMagGraphAnalysisTaskStrs();
                angle = angles(angleInd);

                [depVarValue, depVarUnit] = lvd_GeometricAngleTasks(subLog(i), 'Mag', angle, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPlaneNormXPattern))))
                tokens = regexpi(taskStr, geoPlaneNormXPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                planeInd = str2double(tokens);

                [~, planes] = subLog(i).lvdData.geometry.planes.getPlaneNormalXComponentGraphAnalysisTaskStrs();
                plane = planes(planeInd);

                [depVarValue, depVarUnit] = lvd_GeometricPlaneTasks(subLog(i), 'NormVectorX', plane, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPlaneNormYPattern))))
                tokens = regexpi(taskStr, geoPlaneNormYPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                planeInd = str2double(tokens);

                [~, planes] = subLog(i).lvdData.geometry.planes.getPlaneNormalYComponentGraphAnalysisTaskStrs();
                plane = planes(planeInd);

                [depVarValue, depVarUnit] = lvd_GeometricPlaneTasks(subLog(i), 'NormVectorY', plane, inFrame);

            elseif(not(isempty(regexpi(taskStr, geoPlaneNormZPattern))))
                tokens = regexpi(taskStr, geoPlaneNormZPattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                planeInd = str2double(tokens);

                [~, planes] = subLog(i).lvdData.geometry.planes.getPlaneNormalZComponentGraphAnalysisTaskStrs();
                plane = planes(planeInd);

                [depVarValue, depVarUnit] = lvd_GeometricPlaneTasks(subLog(i), 'NormVectorZ', plane, inFrame);

            elseif(not(isempty(regexpi(taskStr, pluginValuePattern))))
                tokens = regexpi(taskStr, pluginValuePattern, 'tokens');
                tokens = tokens{1};
                tokens = tokens{1};
                pluginInd = str2double(tokens);

                lvdData = subLog(i).lvdData;
                [~, plugins] = lvdData.plugins.getAllPluginGraphAnalysisTaskStrs();
                plugin = plugins(pluginInd);

                depVarValue = lvd_PluginValueTask(subLog(i), 'plugin_value', plugin, lvdData, inFrame);
                depVarUnit = '';
            else
                error('Unknown LVD task string: "%s"', taskStr);
            end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [f] = myrosenbrock(x)
f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;

repo: BUPT-ICS-Courseware


text:
% erplab_eegscanner
%
% PURPOSE: checks EEG dataset according to set of criteria
% Format:
% serror = erplab_eegscanner(EEG, funcname, chckmultieeg, chckemptyeeg, chckepocheeg, chckeventlist)
%
% Inputs:
%
% EEG           - EEGLAB dataset
% funcname      - function name that calls erplab_eegscanner.
%
% Criteria:
% chckmultieeg  = check for multiple dataset (md).    0= do not accept md;    1=accept only md;  2=do not care
% chckemptyeeg  = check for empty dataset (ed).       0= do not accept ed;    1=accept only ed;  2=do not care
% chckepocheeg  = check for epoched dataset (epd).    0= do not accept epd;   1=accept only epd; 2=do not care
% chcknoevents  = check for no event codes (ev).      0= do not accept if no event   1=accept if no events;
% chckeventlist = check for EVENTLIST structure (ES). 0= do not accept ES;    1=accept only ES;  2=do not care
%
% Option (parameter-value)
%
% 'ErrorMsg'   - 'popup': display popup window with error message when criteria are not met
%              - 'error': display error message (command window) when criteria are not met
%              - 'off'  ; % do not display error msg (just get outcome through "serror")
%
%
% Output:
% serror        - 0 means Ok.no problem found. 1 means criteria were not met.
%
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

function serror = erplab_eegscanner(EEG, funcname, chckmultieeg, chckemptyeeg, chckepocheeg, chcknoevents, chckeventlist, varargin)

serror = 0; % no problem by default
if nargin<1
        help erplab_eegscanner
        return
end

%
% Parsing inputs
%
p = inputParser;
p.FunctionName  = mfilename;
p.CaseSensitive = false;
p.addRequired('EEG');
p.addRequired('funcname',@ischar);
p.addRequired('chckmultieeg', @isnumeric);
p.addRequired('chckemptyeeg', @isnumeric);
p.addRequired('chckepocheeg', @isnumeric);
p.addRequired('chcknoevents', @isnumeric);
p.addRequired('chckeventlist',@isnumeric);
% Parameters
p.addParamValue('ErrorMsg', 'popup', @ischar);
p.parse(EEG, funcname, chckmultieeg, chckemptyeeg, chckepocheeg, chcknoevents, chckeventlist, varargin{:});

if strcmpi(p.Results.ErrorMsg, 'popup') || strcmpi(p.Results.ErrorMsg, 'on')
        errormsg = 1; % popup window with error message
elseif strcmpi(p.Results.ErrorMsg, 'error')
        errormsg = 2;       % error message
else % off
        errormsg = 0; % do not display error msg
end

title = sprintf('ERPLAB: Error at %s', funcname);
N = length(EEG);
if chckmultieeg==0 && N>1
        serror = 1;
        msgboxText =  'Unfortunately, %s does not work with multiple datasets';
        if errormsg==1
                errorfound(sprintf(msgboxText, funcname), title);
                return
        elseif errormsg==2
                error('prog:input', msgboxText, funcname)
        end
        return
end
if chckmultieeg==1 && N<2
        serror = 1;
        msgboxText =  'Unfortunately, %s only works with multiple datasets';
        if errormsg==1
                errorfound(sprintf(msgboxText, funcname), title);
                return
        elseif errormsg==2
                error('prog:input', msgboxText, funcname)
        end
        return
end

for k=1:N
        if chckemptyeeg==0 && isempty(EEG(k).data)
                msgboxText =  '%s cannot read an empty dataset!';
                serror = 1;
                break
        end
        if chckemptyeeg==1 && ~isempty(EEG(k).data)
                msgboxText =  '%s only works with empty dataset(s)';
                serror = 1;
                break
        end
        if chckepocheeg==0 && ~isempty(EEG(k).epoch)
                msgboxText =  '%s only works with continuous data';
                serror = 1;
                break
        end
        if chckepocheeg==1 && isempty(EEG(k).epoch)
                msgboxText =  '%s only works with epoched data';
                serror = 1;
                break
        end
        if chcknoevents==0 && isempty(EEG(k).event)
                msgboxText =  '%s : There is not event codes in this dataset!';
                serror = 1;
                break
        end
        if chcknoevents==1 && ~isempty(EEG(k).event)
                msgboxText =  '%s : This dataset has event codes already!';
                serror = 1;
                break
        end
        if chckeventlist~=2
                if isfield(EEG(k), 'EVENTLIST')
                        if chckeventlist==1 % accept only
                                if isempty(EEG(k).EVENTLIST)
                                        msgboxText = '%s : EEG.EVENTLIST structure is empty!';
                                        serror = 1;
                                        break
                                end
                                if isfield(EEG(k).EVENTLIST, 'eventinfo')
                                        if chckeventlist==1 && isempty(EEG(k).EVENTLIST.eventinfo)
                                                msgboxText = '%s : EVENTLIST.eventinfo structure is empty!';
                                                serror = 1;
                                                break
                                        end
                                else
                                        msgboxText =  '%s : EVENTLIST.eventinfo structure was not found!';
                                        serror = 1;
                                        break
                                end
                        else % do not accept EVENTLIST
                                if ~isempty(EEG(k).EVENTLIST)
                                        msgboxText = '%s does not work with EEG.EVENTLIST structure.';
                                        serror = 1;
                                        break
                                end
                        end
                else
                        if chckeventlist==1 % accept only
                                msgboxText =  '%s : EVENTLIST structure was not found!';
                                serror = 1;
                                break
                        end
                end
        end
end
if serror>0
        if errormsg==1
                errorfound(sprintf(msgboxText, funcname), title);
                return
        elseif errormsg==2
                error('prog:input', ['ERPLAB says: ' msgboxText], funcname)
        end
end
return

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [y,deriv] = log_distance_mv2df(w,input_data,new_dim)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% The function projects each column of input_data to a subspace and then
% computes log distance from a centroid. The input_data is fixed, but
% the projection and centroid parameters are variable.
%
%   W = reshape(w);
%   y.' = log sum((W(:,1:end-1).'*input_data - W(:,end)).^2,1)
%
%      W is the augmented matrix [M c] where M maps an input  vector
%      to a lower dimensional space and c is the centroid in
%      the lower dimensional space.
%
% Parameters:
%     w: the vectorized version of the W matrix
%     input_data: is an K-by-T matrix of input vectors of length K, for
%     each of T trials.
%     new_dim: the dimension of vectors in the lower dimensional space.
%

if nargin==0
    test_this();
    return;
end

if isempty(w)
    [dim, num_trials] = size(input_data);
    map = @(w) map_this(w,input_data,dim,new_dim);
    transmap = @(w) transmap_this(w,input_data,num_trials,new_dim);
    delta = linTrans(w,map,transmap);
    y = logsumsquares_fh(new_dim,1,delta);
    return;
end

if isa(w,'function_handle')
    f = log_distance_mv2df([],input_data,new_dim);
    y = compose_mv(f,w,[]);
    return;
end

f = log_distance_mv2df([],input_data,new_dim);
if nargout==1
    y = f(w);
else
    [y,deriv] = f(w);
end

function y = map_this(w,input_data,dim,new_dim)
% V = [input_data; ones(1,num_trials)];
W = reshape(w,new_dim,dim+1);
y = bsxfun(@minus,W(:,1:end-1)*input_data,W(:,end));
y = y(:);

function dx = transmap_this(dy,input_data,num_trials,new_dim)
dY = reshape(dy,new_dim,num_trials);
% V = [input_data; ones(1,num_trials)];
% Vt = V.';
% dX = dY*Vt;
dYt = dY.';
dYtSum = sum(dYt,1);
dX = [input_data*dYt;-dYtSum].';
dx = dX(:);

function test_this()
K = 5;
N = 10;
P = 3;
M = randn(P,N);
c = randn(P,1);
W = [M c];
w = W(:);
input_data = randn(N,K);

f = log_distance_mv2df([],input_data,P);
test_MV2DF(f,w);

author: nesl


text:
function G =  getGroupOverlapColor(row, col)

    N = row*col;

    g = sparse(zeros(N,1));
    g = diag(g);

    %% build overlapping group

    % top let corner group

    g(1, 1 ) = 1;
    g(1, 2 ) = 1;
    g(1, 1+row ) = 1;
    g(1, 2+row ) = 1;

     % bottom let corner group

    g(row, row-1 ) = 1;
    g(row, row ) = 1;
    g(row, row-1+row ) = 1;
    g(row, row+row ) = 1;

    % top right corner group

    g((col-1)*row+1, (col-1)*row+1 ) = 1;
    g((col-1)*row+1, (col-2)*row+1 ) = 1;
    g((col-1)*row+1, (col-1)*row+2 ) = 1;
    g((col-1)*row+1, (col-2)*row+2 ) = 1;

    % bottom right corner group

    g((col-1)*row+row , (col-1)*row+row-1 ) = 1;
    g((col-1)*row+row , (col-2)*row+row-1 ) = 1;
    g((col-1)*row+row , (col-1)*row+row ) = 1;
    g((col-1)*row+row , (col-2)*row+row ) = 1;

    % boundary group

    for i=2:col-1
        % top row rgoup

        j = 1;

        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;

        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;

        g( (i-1)*row+j  , i*row+j   ) = 1;
        g( (i-1)*row+j  , i*row+j+1 ) = 1;

        % bottom row group

        j = row;

        g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;

        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;

        g( (i-1)*row+j  , i*row+j-1 ) = 1;
        g( (i-1)*row+j  , i*row+j   ) = 1;

    end

    for j=2:row-1
          % left column rgoup

        i=1;

        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;
        g( (i-1)*row+j  , i*row+j-1 ) = 1;
        g( (i-1)*row+j  , i*row+j   ) = 1;
        g( (i-1)*row+j  , i*row+j+1 ) = 1;

        % right  column group

        i=col;

        g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;

    end

    for i=2:col-1
        for j=2:row-1

            g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
            g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
            g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;
            g( (i-1)*row+j  , i*row+j-1 ) = 1;
            g( (i-1)*row+j  , i*row+j   ) = 1;
            g( (i-1)*row+j  , i*row+j+1 ) = 1;

        end
    end

    G = sparse(3*N, 3*N);

    G(1:N, 1:N) = g;
    G(N+1:2*N, N+1:2*N) =g;
    G(2*N+1:3*N, 2*N+1:3*N) =g;

    % G =[g; g; g];

end

author: andrewssobral

repo: lrslibrary