text:
%EIMM_DEMO2 Bearings Only Tracking of a Manouvering Target demonstration
%
% Simple demonstration for non-linear IMM using the following models:
%  1. Standard Wiener process velocity model
%  2. Coordinated turn model
%
% The measurement model is non-linear bearings only model.
% See BOT-demo or the documentation for more details.
%
% Copyright (C) 2007-2008 Jouni Hartikainen
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

print_figures = 1;
save_figures = 1;

% Dimensionality of the state space
dims = 5;

% The number of models in use
nmodels = 2;

% Step size
dt = 0.1;

% Space for function handles and parameters
a_func = {};
ia_func = {};
a_param = {};
h_func = {};
h_param = {};
dh_dx_func = {};

a_func{1} = [];
a_func{2} = @f_turn;
ia_func{1} = [];
ia_func{2} = @f_turn_inv;
a_param{1} = [];
a_param{2} = {dt};

% Space for model parameters
ind = cell(1,nmodels);
F   = cell(1,nmodels);
L   = cell(1,nmodels);
Qc  = cell(1,nmodels);
A   = cell(1,nmodels);
Q   = cell(1,nmodels);
H   = cell(1,nmodels);
R   = cell(1,nmodels);

% Index vector of model 1
ind{1} = [1 2 3 4]';

% Transition matrix for the continous-time velocity model.
F{1} = [0 0 1 0;
        0 0 0 1;
        0 0 0 0;
        0 0 0 0];

% Noise effect matrix for the continous-time system.
L{1} = [0 0;
        0 0;
        1 0;
        0 1];

% Process noise variance
q1 = 0.01;
Qc{1} = diag([q1 q1]);

% Discretization of the continous-time system.
[A{1},Q{1}] = lti_disc(F{1},L{1},Qc{1},dt);

% Process noise variance for model 1 using EKF
EKF_q1 = .05;
EKF_Qc1 = diag([EKF_q1 EKF_q1]);

% Discretization of the continous-time system.
[EKF_A1,EKF_Q1] = lti_disc(F{1},L{1},EKF_Qc1,dt);

%%% Specification of the turning model

% System components. 5th parameter is the turning rate
ind{2} = [1 2 3 4 5]';

% Derivative of the dynamic function
A{2} = @f_turn_dx;
% Process noise for the turning rate
Qc{2} = 0.15;

% Noise effect matrix
L{2} = [0 0 0 0 1]';

% Process noise covariance
Q{2} = L{2}*Qc{2}*L{2}'*dt;

hdims = 2;

%mu_ip = [0.90 0.05 0.05];
mu_ip = [0.95 0.05];
mu_0j = mu_ip;
%p_ij = [0.65 0.35;
%        0.10 0.90];

p_ij = [0.90 0.10;
        0.10 0.90];

% Number of data points
n = 200;

% Space for real states and modes
X_r = zeros(dims,n);
mstate = zeros(1,n);

%%%%%%% Creation of trajectory %%%%%%%

% Start with constant velocity 1 toward right
mstate(1:40) = 1;
X_r(:,1) = [0 0 1 0 0]';

% At 4s make a turn left with rate 1
mstate(41:90) = 2;
X_r(5,40) = 1;

% At 9s move straight for 2 seconds
mstate(91:110) = 1;

% At 11s commence another turn right with rate -1
mstate(111:160) = 2;
X_r(5,110) = -1;

% At 16s move straight for 4 seconds
mstate(161:200) = 1;

% Generate object state values
for i = 2:n
   st = mstate(i);
   if isstr(a_func{st}) | strcmp(class(a_func{st}),'function_handle')
       X_r(ind{st},i) = feval(a_func{st},X_r(ind{st},i-1),a_param{st});
   else
       X_r(ind{st},i) = A{st}*X_r(ind{st},i-1);
   end
end

% Positions of sensors
% $$$ S1 = [-0.5; 3];
% $$$ S2 = [-0.5;-3];
% $$$ S3 = [   7;-3];
% $$$ S4 = [   7; 3];

S1 = [-0.5; 3.5];
S2 = [-0.5;-3.5];
S3 = [   7;-3.5];
S4 = [   7; 3.5];

s = [S1 S2 S3 S4];

% Handles to measurement models
H{1} = @bot_dh_dx;
H{2} = @bot_dh_dx;

h_func{1} = @bot_h;
h_func{2} = @bot_h;

h_param{1} = s;
h_param{2} = s;

% Standard deviation
sd = 0.1*ones(1,size(s,2));
R = {};
R{1} = diag(sd.^2);
R{2} = R{1};
% Generate measurements
Y = bot_h(X_r,s);
% Add noise
for i = 1:size(Y,2)
    st = mstate(i);
    for j = 1:size(Y,1)
        Y(j,i) = Y(j,i) + sqrt(R{st}(j,j)) * randn;
    end
end

% Print the trajectory
if print_figures
    h = plot(X_r(1,:),X_r(2,:),'-g',...
             s(1,:),s(2,:),'k^',...
             X_r(1,1),X_r(2,1),'ro','MarkerSize',12);
    legend('Real trajectory',...
           'Positions of sensors',...
           'Starting position', 'Location', 'North');
    set(h,'markersize',5);
    set(h,'linewidth',0.5);
    set(gca,'FontSize',8);
    xlim([-1 7.5])
    ylim([-4 4])
    if save_figures
        print('-depsc','eimm2_trajectory.eps');
    end
    pause
end

m = [0 0 0 -1 0]';
P = diag([10.1 10.1 1.1 1.1 1]);

%% Space for the estimates.

% EKF with model 1
EKF_MM = zeros(size(A{1},1), size(Y,2));
EKF_PP = zeros(size(A{1},1), size(A{1},1), size(Y,2));

% UKF with model 1
UKF_MM = zeros(size(A{1},1), size(Y,2));
UKF_PP = zeros(size(A{1},1), size(A{1},1), size(Y,2));

% EKF based IMM
EIMM_MM = zeros(size(m,1), size(Y,2));
EIMM_PP = zeros(size(m,1), size(m,1), size(Y,2));
EIMM_MM_i = cell(nmodels,n);
EIMM_PP_i = cell(nmodels,n);
EIMM_MU = zeros(nmodels,size(Y,2));

% UKF based IMM
UIMM_MM = zeros(size(m,1), size(Y,2));
UIMM_PP = zeros(size(m,1), size(m,1), size(Y,2));
UIMM_MM_i = cell(nmodels,n);
UIMM_PP_i = cell(nmodels,n);
UIMM_MU = zeros(nmodels,size(Y,2));

%%% Initial estimates %%%

% EKF with model 1
EKF_M = [0 0 0 -1]';
EKF_P = diag([1.1 1.1 0.1 0.1]);

% UKF with model 1
UKF_M = [0 0 0 -1]';
UKF_P = diag([1.1 1.1 0.1 0.1]);

% EKF based IMM
x_ip1{1} = [0 0 1 0]';
x_ip1{2} = [0 0 1 0 0]';
mu_ip1 = mu_ip;

P_ip1{1} = diag([0.1 0.1 0.1 0.1]);
P_ip1{2} = diag([0.1 0.1 0.1 0.1 1]);

% UKF based IMM
x_ip2{1} = [0 0 1 0]';
x_ip2{2} = [0 0 1 0 0]';
mu_ip2 = mu_ip;

P_ip2{1} = diag([0.1 0.1 0.1 0.1]);
P_ip2{2} = diag([0.1 0.1 0.1 0.1 1]);

% Filtering steps.
for i = 1:size(Y,2)
    % EKF with model 1
    [EKF_M,EKF_P] = kf_predict(EKF_M,EKF_P,EKF_A1,EKF_Q1);
    [EKF_M,EKF_P] = ekf_update1(EKF_M,EKF_P,Y(:,i),H{1},R{1},h_func{1},[],h_param{1});

    EKF_MM(:,i)   = EKF_M;
    EKF_PP(:,:,i) = EKF_P;

    % UKF with model 2
    [UKF_M,UKF_P] = kf_predict(UKF_M,UKF_P,EKF_A1,EKF_Q1);
    [UKF_M,UKF_P] = ukf_update1(UKF_M,UKF_P,Y(:,i),h_func{1},R{1},h_param{1});

    UKF_MM(:,i)   = UKF_M;
    UKF_PP(:,:,i) = UKF_P;

    % EKF based IMM
    [x_p1,P_p1,c_j1] = eimm_predict(x_ip1,P_ip1,mu_ip1,p_ij,ind,dims,A,a_func,a_param,Q);
    [x_ip1,P_ip1,mu_ip1,m1,P1] = eimm_update(x_p1,P_p1,c_j1,ind,dims,Y(:,i),H,h_func,R,h_param);
    EIMM_MM(:,i)   = m1;
    EIMM_PP(:,:,i) = P1;
    EIMM_MU(:,i)   = mu_ip1';
    EIMM_MM_i(:,i) = x_ip1';
    EIMM_PP_i(:,i) = P_ip1';

    % UKF based IMM
    [x_p2,P_p2,c_j2] = uimm_predict(x_ip2,P_ip2,mu_ip2,p_ij,ind,dims,A,a_func,a_param,Q);
    [x_ip2,P_ip2,mu_ip2,m2,P2] = uimm_update(x_p2,P_p2,c_j2,ind,dims,Y(:,i),H,h_func,R,h_param);

    UIMM_MM(:,i)   = m2;
    UIMM_PP(:,:,i) = P2;
    UIMM_MU(:,i)   = mu_ip2';
    UIMM_MM_i(:,i) = x_ip2';
    UIMM_PP_i(:,i) = P_ip2';

    % Plot the estimates so far
    if print_figures
        plot(EKF_MM(1,1:i),EKF_MM(2,1:i),'y-',...
             EIMM_MM(1,1:i),EIMM_MM(2,1:i),'r-',...
             UIMM_MM(1,1:i),UIMM_MM(2,1:i),'b-',...
             X_r(1,1:i),X_r(2,1:i),'g-');

        % Measurement directions
        hold on
        for k = 1:size(s,2)
            len = sqrt(sum((X_r(1:2,i)-s(:,k)).^2,1));
            dx = len*cos(Y(k,i));
            dy = len*sin(Y(k,i));

            plot([s(1,k);s(1,k)+dx], [s(2,k);s(2,k)+dy], 'k--')
        end
        plot(s(1,:),s(2,:),'k^')
        xlim([-1 7.5])
        ylim([-4 4])
        hold off
        drawnow
    end
end

% Smooth with EKF based IMM smoother
[SMI_1,SPI_1,SMI_i_1,SPI_i_1,MU_S_1] = eimm_smooth(EIMM_MM,EIMM_PP,EIMM_MM_i,EIMM_PP_i,EIMM_MU,p_ij,mu_0j,ind,dims,A,ia_func,a_param,Q,R,H,h_func,h_param,Y);

% Smooth with UKF based IMM smoother
[SMI_2,SPI_2,SMI_i_2,SPI_i_2,MU_S_2] = uimm_smooth(UIMM_MM,UIMM_PP,UIMM_MM_i,UIMM_PP_i,UIMM_MU,p_ij,mu_0j,ind,dims,A,ia_func,a_param,Q,R,H,h_func,h_param,Y);

% Smooth the EKF estimates with RTS smoother using model 1
[SM1,SP1] = rts_smooth(EKF_MM,EKF_PP,EKF_A1,EKF_Q1);

% Smooth the UKF estimates with RTS smoother using model 1
[SM2,SP2] = rts_smooth(UKF_MM,UKF_PP,EKF_A1,EKF_Q1);

% Calculate the MSEs
MSE_EKF1_1 = mean((X_r(1,:)-EKF_MM(1,:)).^2);
MSE_EKF1_2 = mean((X_r(2,:)-EKF_MM(2,:)).^2);
MSE_EKF1 = 1/2*(MSE_EKF1_1 + MSE_EKF1_2);

MSE_EKS1_1 = mean((X_r(1,:)-SM1(1,:)).^2);
MSE_EKS1_2 = mean((X_r(2,:)-SM1(2,:)).^2);
MSE_EKS1 = 1/2*(MSE_EKS1_1 + MSE_EKS1_2);

% Calculate the MSEs
MSE_UKF1_1 = mean((X_r(1,:)-UKF_MM(1,:)).^2);
MSE_UKF1_2 = mean((X_r(2,:)-UKF_MM(2,:)).^2);
MSE_UKF1 = 1/2*(MSE_UKF1_1 + MSE_UKF1_2);

MSE_UKS1_1 = mean((X_r(1,:)-SM2(1,:)).^2);
MSE_UKS1_2 = mean((X_r(2,:)-SM2(2,:)).^2);
MSE_UKS1 = 1/2*(MSE_UKS1_1 + MSE_UKS1_2);

MSE_EIMM1 = mean((X_r(1,:)-EIMM_MM(1,:)).^2);
MSE_EIMM2 = mean((X_r(2,:)-EIMM_MM(2,:)).^2);
MSE_EIMM = 1/2*(MSE_EIMM1 + MSE_EIMM2);

MSE_EIMMS1 = mean((X_r(1,:)-SMI_1(1,:)).^2);
MSE_EIMMS2 = mean((X_r(2,:)-SMI_1(2,:)).^2);
MSE_EIMMS = 1/2*(MSE_EIMMS1 + MSE_EIMMS2);

MSE_UIMM1 = mean((X_r(1,:)-UIMM_MM(1,:)).^2);
MSE_UIMM2 = mean((X_r(2,:)-UIMM_MM(2,:)).^2);
MSE_UIMM = 1/2*(MSE_UIMM1 + MSE_UIMM2);

MSE_UIMMS1 = mean((X_r(1,:)-SMI_2(1,:)).^2);
MSE_UIMMS2 = mean((X_r(2,:)-SMI_2(2,:)).^2);
MSE_UIMMS = 1/2*(MSE_UIMMS1 + MSE_UIMMS2);

fprintf('Mean square errors of position estimates:\n');
fprintf('EKF1-RMSE = %.4f\n',MSE_EKF1);
fprintf('EKS1-RMSE = %.4f\n',MSE_EKS1);
fprintf('UKF1-RMSE = %.4f\n',MSE_UKF1);
fprintf('UKS1-RMSE = %.4f\n',MSE_UKS1);
fprintf('EIMM-RMSE = %.4f\n',MSE_EIMM);
fprintf('EIMMS-RMSE = %.4f\n',MSE_EIMMS);
fprintf('UIMM-RMSE = %.4f\n',MSE_UIMM);
fprintf('UIMMS-RMSE = %.4f\n',MSE_UIMMS);

% Plot the final filtering and smoothing results
if print_figures
    h = plot(X_r(1,:),X_r(2,:),'g-',...
             EKF_MM(1,:),EKF_MM(2,:),'-k',...
             SM1(1,:),SM1(2,:),'--k',...
             EIMM_MM(1,:),EIMM_MM(2,:),'-r',...
             SMI_1(1,:),SMI_1(2,:),'--r',...
             UIMM_MM(1,:),UIMM_MM(2,:),'-b',...
             SMI_2(1,:),SMI_2(2,:),'--b');
    legend('True trajectory',...
           'EKF',...
           'ERTS',...
           'IMM-EKF',...
           'IMM-EKS',...
           'IMM-UKF',...
           'IMM-UKS');
    %title('Estimates produced by IMM-filter.')
    set(h,'markersize',2);
    set(h,'linewidth',0.5);
    set(gca,'FontSize',8);
    xlim([-1 7.5])
    ylim([-4 4])
    if save_figures
        print('-depsc','eimm2_1.eps');
    end
    pause

    % Determine the real model probabilities
    p_models = zeros(nmodels,n);
    I1 = find(mstate == 1);
    p_models(1,I1) = 1;
    I2 = find(mstate == 2);
    p_models(2,I2) = 1;

    % Plot model 2 probability for filters
    h = plot(1:n, p_models(2,:),'g--',...
             1:n,EIMM_MU(2,:)','-r',...
             1:n,MU_S_1(2,:)','--r',...
             1:n,UIMM_MU(2,:)','-b',...
             1:n,MU_S_2(2,:)','--b');
    %1:n,MU_S_2(2,:)','b-');
    legend('True',...
           'IMM-EKF',...
           'IMM-EKS',...
           'IMM-UKF',...
           'IMM-UKS');
    %title('Probability of model 2');
    ylim([-0.1,1.1]);
    set(h,'markersize',2);
    set(h,'linewidth',0.5);
    set(gca,'FontSize',8);
    if save_figures
        print('-depsc','eimm2_2.eps');
    end
    pause

    % Collect the turn rate estimates
    EMM_t = zeros(5,n);
    EMM_s = zeros(5,n);
    UMM_t = zeros(5,n);
    UMM_s = zeros(5,n);
    for i = 1:n
        EMM_t(:,i) = EIMM_MM_i{2,i};
        EMM_s(:,i) = SMI_i_1{2,i};
        UMM_t(:,i) = UIMM_MM_i{2,i};
        UMM_s(:,i) = SMI_i_2{2,i};
    end

    % Plot the filtered turn rates
    h = plot(1:n, X_r(5,:),'-g',...
             1:n,EMM_t(5,:),'-r',...
             1:n,EMM_s(5,:),'--r',...
             1:n,UMM_t(5,:),'-b',...
             1:n,UMM_s(5,:),'--b');
    %title('Turn rate estimates')
    legend('True',...
           'IMM-EKF',...
           'IMM-EKS',...
           'IMM-UKF',...
           'IMM-UKS');
    set(h,'markersize',2);
    set(h,'linewidth',0.5);
    set(gca,'FontSize',8);
    if save_figures
        print('-depsc','eimm2_3.eps');
    end

end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function plot_results(varargin)
%PLOT_RESULTS Plots all results for a given dataset.
%   Name-value pairs are supported, listed below.
%   E.g.: plot_results('dataset', 'imagenet', 'set', 'val')

  opts.dataset = 'cifar';  % choose dataset
  opts.set = 'train';  % train or val
  opts.var = 'err';  % err or obj (error/objective)
  opts.save = false;  % if true, saves as a PDF file in the base directory

  opts.baseDir = [vl_rootnn() '/data/curveball/<dataset>'];  % results location
  opts.subdir = '';  % optional subdirectory

  opts = vl_argparse(opts, varargin);

  % compose directory where experiments are located
  base = strrep(opts.baseDir, '<dataset>', opts.dataset);
  base = [base '/' opts.subdir];

  % list experiments (folders)
  folders = dir(base);

  figure(10);
  clf();
  hold('on');
  ax = gca();
  set(ax, 'lineStyleOrder', {'-','--','-.',':','*','+'});

  for i = 1:numel(folders)
    name = folders(i).name;

    if folders(i).isdir && ~strncmp(name, '.', 1)
      % try to load final results first
      filename = [base '/' name '/results.mat'];
      if ~exist(filename, 'file')
        % if it doesn't exist, look for the most recent checkpoint
        filename = models.checkpoint([base '/' name '/epoch-*.mat']);
      end

      if ~isempty(filename)
        % load result statistics
        data = load(filename, 'stats');

        % get the relevant curve
        values = data.stats.values(opts.set, opts.var);

        % display it
        plot(values, 'DisplayName', name, 'LineWidth',2);
      end
    end
  end

  % make plots look nice
  legend('show');
  grid('on');

  xlabel('Epochs');
  ylabel([opts.set ' ' opts.var]);

  % save to PDF
  if opts.save
    filename = [base '/results.pdf'];
    print(10, filename, '-dpdf');
  end
end

repo: curveball


text:
function Reaction_time_task_GUI(varargin)
fh = figure('Visible','off', ...
            'HandleVisibility','callback', ...
            'NumberTitle','off', ...  % turns off MATLAB GUI window title
            'Name','Reaction_time_task_GUI', ...    % now, define my own
            'Units','normalized', ...
            'Position',[1.000000e-001 1.000000e-001 6.000000e-001 6.000000e-001], ...
            'Color','white', ...   % match bg of MODE_new.jpg
            'Resize','off');
set(fh,'MenuBar','none');  % removes MATLAB default menu bar
x0 = 0.1; y0 = 0.85;
x0run = 0.825; y0run = 0.075; % lower-left corner coord. of the run button
dxst1 = 0.075; dyst1 = 0.05;  % for static text 1
dxet1 = 0.100; dyet1 = 0.05;  % for edit text 1
dxst2 = 0.400; dyst2 = 0.10;  % for static text 2
dxpb1 = 0.050; dypb1 = 0.05;  % for PDF logo pushbuttons
dxrun = 0.100; dyrun = 0.075; % for "Run" pushbutton
% First, create a Heading for the secondary GUI with a static text box
dx = 0.6; dy = 0.125;
st0h = uicontrol(fh, ...
            'style','text', ...
            'String',['Lateral intraparietal area' char(10) '(Eq. 28, p. 1370: Grossberg & Pilly (2008))'], ...
            'Units','normalized', ...
            'Position',[x0,y0,dx,dy], ...
            'FontSize',12, ...
            'FontWeight','bold');

% Create static box, edit box, static box, logo (that points to a doc; optional)
y0 = y0 - dy*1.05; space = dxst1*0.20;
data.str{1} = 'A9';
data.default{1} = 1.000000;
  st1h(1) = uicontrol(fh,'style','text', ...
      'String','A9', ...
      'Units','normalized', ...
      'Position',[x0,y0,dxst1,dyst1], ...
      'FontSize',10, ...
      'FontWeight','bold');
  et1h(1) = uicontrol(fh,'Style','edit', ...
      'String','1', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+space,y0,dxet1,dyet1], ...
      'FontSize',10, ...
      'FontWeight','normal', ...
      'callback',{@edit1_Callback 'A9'});
  st2h(1) = uicontrol(fh,'style','text', ...
      'String','Parameter that ...', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...
      'FontSize',10, ...
      'FontWeight','bold');
  pb1h(1) = uicontrol(fh,'Style','pushbutton', ...
      'String','PDF', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...
      'FontSize',10, ...
      'FontWeight','bold', ...
      'callback',{@opendoc1_Callback 'A9_heur.pdf'});
  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;
data.str{2} = 'g_f';
data.default{2} = 1.000000;
  st1h(2) = uicontrol(fh,'style','text', ...
      'String','g_f', ...
      'Units','normalized', ...
      'Position',[x0,y0,dxst1,dyst1], ...
      'FontSize',10, ...
      'FontWeight','bold');
  et1h(2) = uicontrol(fh,'Style','edit', ...
      'String','1', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+space,y0,dxet1,dyet1], ...
      'FontSize',10, ...
      'FontWeight','normal', ...
      'callback',{@edit1_Callback 'g_f'});
  st2h(2) = uicontrol(fh,'style','text', ...
      'String','Gain of recurrent self-excitation ...', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...
      'FontSize',10, ...
      'FontWeight','bold');
  pb1h(2) = uicontrol(fh,'Style','pushbutton', ...
      'String','PDF', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...
      'FontSize',10, ...
      'FontWeight','bold', ...
      'callback',{@opendoc1_Callback 'g_f_heur.pdf'});
  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;
data.str{3} = 'g_h';
data.default{3} = 5.000000;
  st1h(3) = uicontrol(fh,'style','text', ...
      'String','g_h', ...
      'Units','normalized', ...
      'Position',[x0,y0,dxst1,dyst1], ...
      'FontSize',10, ...
      'FontWeight','bold');
  et1h(3) = uicontrol(fh,'Style','edit', ...
      'String','5', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+space,y0,dxet1,dyet1], ...
      'FontSize',10, ...
      'FontWeight','normal', ...
      'callback',{@edit1_Callback 'g_h'});
  st2h(3) = uicontrol(fh,'style','text', ...
      'String','Gain of recurrent inhibition ...', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...
      'FontSize',10, ...
      'FontWeight','bold');
  pb1h(3) = uicontrol(fh,'Style','pushbutton', ...
      'String','PDF', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...
      'FontSize',10, ...
      'FontWeight','bold', ...
      'callback',{@opendoc1_Callback 'g_h_heur.pdf'});
  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;
data.str{4} = 'lambda';
data.default{4} = 5.000000;
  st1h(4) = uicontrol(fh,'style','text', ...
      'String','lambda', ...
      'Units','normalized', ...
      'Position',[x0,y0,dxst1,dyst1], ...
      'FontSize',10, ...
      'FontWeight','bold');
  et1h(4) = uicontrol(fh,'Style','edit', ...
      'String','5', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+space,y0,dxet1,dyet1], ...
      'FontSize',10, ...
      'FontWeight','normal', ...
      'callback',{@edit1_Callback 'lambda'});
  st2h(4) = uicontrol(fh,'style','text', ...
      'String','Gain of bottom-up excitation ...', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...
      'FontSize',10, ...
      'FontWeight','bold');
  pb1h(4) = uicontrol(fh,'Style','pushbutton', ...
      'String','PDF', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...
      'FontSize',10, ...
      'FontWeight','bold', ...
      'callback',{@opendoc1_Callback 'lambda_heur.pdf'});
  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;
data.str{5} = 'sigLIP';
data.default{5} = 5.000000;
  st1h(5) = uicontrol(fh,'style','text', ...
      'String','sigLIP', ...
      'Units','normalized', ...
      'Position',[x0,y0,dxst1,dyst1], ...
      'FontSize',10, ...
      'FontWeight','bold');
  et1h(5) = uicontrol(fh,'Style','edit', ...
      'String','5', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+space,y0,dxet1,dyet1], ...
      'FontSize',10, ...
      'FontWeight','normal', ...
      'callback',{@edit1_Callback 'sigLIP'});
  st2h(5) = uicontrol(fh,'style','text', ...
      'String','Gain of recurrent ...', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...
      'FontSize',10, ...
      'FontWeight','bold');
  pb1h(5) = uicontrol(fh,'Style','pushbutton', ...
      'String','PDF', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...
      'FontSize',10, ...
      'FontWeight','bold', ...
      'callback',{@opendoc1_Callback 'sigLIP_heur.pdf'});
  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;
data.str{6} = 'Tdec';
data.default{6} = 55.000000;
  st1h(6) = uicontrol(fh,'style','text', ...
      'String','Tdec', ...
      'Units','normalized', ...
      'Position',[x0,y0,dxst1,dyst1], ...
      'FontSize',10, ...
      'FontWeight','bold');
  et1h(6) = uicontrol(fh,'Style','edit', ...
      'String','55', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+space,y0,dxet1,dyet1], ...
      'FontSize',10, ...
      'FontWeight','normal', ...
      'callback',{@edit1_Callback 'Tdec'});
  st2h(6) = uicontrol(fh,'style','text', ...
      'String','Threshold LIP ...', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+space*2,y0-0.02,dxst2,dyst2], ...
      'FontSize',10, ...
      'FontWeight','bold');
  pb1h(6) = uicontrol(fh,'Style','pushbutton', ...
      'String','PDF', ...
      'Units','normalized', ...
      'Position',[x0+dxst1+dxet1+dxst2+space*3,y0,dxpb1,dypb1], ...
      'FontSize',10, ...
      'FontWeight','bold', ...
      'callback',{@opendoc1_Callback 'Tdec_heur.pdf'});
  y0 = y0 - max([dyst1 dyst2 dyet1 dypb1])*1.1;
% save the edit text handles to data
data.et1h = et1h(1:end);
for i=1:size(data.str,1)
   eval(['data.' data.str{i} '= data.default{i};']);
end
data.guiName = 'Reaction_time_task_GUI';
setappdata(fh, 'myData', data);
%align(st0h,'HorizontalAlignment','center');
st3h = uicontrol(fh,'style','text', ...
      'String','After changing each parameter, press "Enter".', ...
      'Units','normalized', ...
      'Position',[0.05,0.08,0.60,0.05], ...
      'HorizontalAlignment','Left', ...
      'foregroundColor',[1 0 0], ...
      'FontSize',10, ...
      'FontWeight','bold');
st4h = uicontrol(fh,'style','text', ...
      'String','To see each parameter heuristic, please click on the PDF logo.', ...
      'Units','normalized', ...
      'Position',[0.05,0.02,0.60,0.05], ...
      'HorizontalAlignment','Left', ...
      'foregroundColor',[1 0 0], ...
      'FontSize',10, ...
      'FontWeight','bold');

runh = uicontrol(fh,'Style','pushbutton', ...
      'String','RUN', ...
      'Units','normalized', ...
      'Position',[x0run,y0run,dxrun,dyrun], ...
      'FontSize',12, ...
      'FontWeight','bold', ...
      'callback',{@run1_Callback 'RTtasknew(A9,g_f,g_h,lambda,sigLIP,Tdec)'});
PDFbuttons(pb1h);

etbgcolor = get(et1h(1),'BackgroundColor');
set(fh,'Color',etbgcolor);
set(fh,'Visible','on');

function opendoc1_Callback(hObject,eventdata,varargin)
% hObject    handle to Exit_menu_item (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
open(varargin{1});

function edit1_Callback(hObject, eventdata, varargin)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB

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

fh = get(hObject,'parent');
data = getappdata(fh, 'myData');
value = str2num(get(hObject,'String'));
eval(['data.' varargin{1} '= value;']);
setappdata(fh,'myData',data);

function run1_Callback(hObject, eventdata, varargin)
% hObject    handle to run1_Callback
% eventdata  reserved - to be defined in a future version of MATLAB
% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double
fh = get(hObject,'parent');
data = getappdata(fh, 'myData');
for i=1:length(data.str)
   eval([data.str{i} '= data.' data.str{i} ';']);  % e.g., A9 = data.A9;
end
% Prevent user from changing the contents of edit boxes
editBoxEnable2(data,'Off');
old_color=get(hObject,'BackgroundColor'); % save original button color
set(hObject, 'BackgroundColor', [1 0 0]); % Change button color
set(hObject, 'String', 'running ...');
% open('RTcaption.pdf');
drawnow    % to make sure that the "Run" button change color right away
eval(varargin{1})   % run the m-file specified by varargin{1}
set(hObject, 'String', 'Run');
set(hObject, 'BackgroundColor', old_color);
editBoxEnable2(data,'On');

author: Sable

repo: mcbench-benchmarks


text:
function u = tora_control(x)

load('tora_data', 'K', 'U', 'domain');

p = [x(3) x(4)];
u = tpcontroller(p, x, K, U, domain);


text:
classdef nnoffset < nntest
  methods (Test)
    function basic(test)
      param = [.34, .5] ;
      x = test.randn(4,5,10,3) ;
      y = vl_nnnoffset(x,param) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nnnoffset(x,param,dzdy) ;
      test.der(@(x) vl_nnnoffset(x,param), x, dzdy, dzdx, 1e-3*test.range) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [ w, xy ] = triangle_unit_o06 ( )

%*****************************************************************************80
%
%% TRIANGLE_UNIT_O06 returns a 6 point quadrature rule for the unit triangle.
%
%  Discussion:
%
%    This rule is precise for monomials through degree 4.
%
%    The integration region is:
%
%      0 <= X
%      0 <= Y
%      X + Y <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(6), the weights.
%
%    Output, real XY(2,6), the abscissas.
%
  w(1:6,1) = [ ...
    0.22338158967801146570, ...
    0.22338158967801146570, ...
    0.22338158967801146570, ...
    0.10995174365532186764, ...
    0.10995174365532186764, ...
    0.10995174365532186764 ];

  xy(1:2,1:6) = [ ...
    0.10810301816807022736,  0.44594849091596488632; ...
    0.44594849091596488632,  0.10810301816807022736; ...
    0.44594849091596488632,  0.44594849091596488632; ...
    0.81684757298045851308,  0.091576213509770743460; ...
    0.091576213509770743460,  0.81684757298045851308; ...
    0.091576213509770743460,  0.091576213509770743460 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s=get_name(a)
s=[get_name(a.algorithm)];
%if nargin==2
s=[s '(' get_name(a.child) ')'];
s=[s sprintf(' bags=%d m=%d',a.bags,a.m)];
%end
eval_name

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef LaunchVehicleStopwatchState < matlab.mixin.SetGet & matlab.mixin.Copyable
    %LaunchVehicleStopwatchState Summary of this class goes here
    %   Detailed explanation goes here

    properties
        stopwatch LaunchVehicleStopwatch
        value(1,1) double = 0;
        running(1,1) StopwatchRunningEnum = StopwatchRunningEnum.NotRunning;
    end

    methods
        function obj = LaunchVehicleStopwatchState(stopwatch)
            obj.stopwatch = stopwatch;
        end
    end

    methods(Access=protected)
        function copyObj = copyElement(obj)
            copyObj = copyElement@matlab.mixin.Copyable(obj);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ad,omega] = calcAngleDistribution(oR,varargin)
% compute the angle distribution of a uniform ODF for a crystal symmetry
%
% Syntax
%   [ad,omega] = calcAngleDistribution(oR)
%   [ad,omega] = calcAngleDistribution(oR,omega)
%
% Input
%  oR    - @orientationRegion
%  omega - angle
%
% Output
%  ad - angle distribution
%  omega - angles
%

if isempty(varargin)
  omega = linspace(0,oR.maxAngle,300);
else
  % restrict omega
  omega = varargin{1};
  omega = omega(omega < oR.maxAngle + 1e-8);
end

ad = zeros(size(omega));
sR = oR.axisSector;
S2G = equispacedS2Grid(sR,'resolution',0.05*degree);

for i = 1:length(omega)
  sR = oR.axisSector(omega(i));
  ad(i) = 2 * volume(sR,S2G) * sin(omega(i)/2)^2;
end

ad = ad ./ mean(ad);

author: mtex-toolbox

repo: mtex

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


text:
% KRON   Kronecker tensor product.
%    KRON(X,Y) is the Kronecker tensor product of X and Y.
%    The result is a large matrix formed by taking all possible
%    products between the elements of X and those of Y. For
%    example, if X is 2 by 3, then KRON(X,Y) is
%
%       [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y
%         X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]
%
%    If either X or Y is sparse, only nonzero elements are multiplied
%    in the computation, and the result is sparse.
%
%    Class support for inputs X,Y:
%       float: double, single
%       integers: uint8, int8, uint16, int16, uint32, int32, uint64, int64
%
%    Reference page in Doc Center
%       doc kron
%
%    Other functions named kron
%
%       splanar/kron    sym/kron
%

author: jmaih


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function cmap_proba = apply_sigmoid(cmap, thr, fq)
 if nargin<3, fq=1; end,
 if nargin<2, thr =0.5; end
 cmap_proba=cmap;
 cmap_proba(:) = 1./(1+exp(-fq*(cmap(:)-thr)));

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function pg = msTest(imsegs, data, maps, labelclassifier, segclassifier, normalize)
% [vacc, hacc, vcm, hcm] = testMultipleSegmentationsCV2(imsegs, labdata, segdata, maps, vclassifier, hclassifier, sclassifier, ncv)

if ~exist('normalize', 'var') || isempty(normalize)
    normalize = 1;
end

pg = cell(numel(imsegs), 1);
allcount = 0;

nclasses = size(labelclassifier.wcs,2);

for f = 1:numel(imsegs)

    %disp(num2str(f))

    nsp = imsegs(f).nseg;

    pg{f} = zeros(nsp, nclasses);

    segs = cell(nsp, 1);

    count = zeros(nsp, 1);

    if ~normalize
        sum_s = zeros(nsp, 1);
    end

    for k = 1:size(data, 2)
        if numel(data{f,k})>0

        isvalid = false(size(data{f, k}, 1), 1);
        for s = 1:size(data{f, k}, 1)
            [segs, ind] = checksegs(segs, maps{f}(:, k), s);
            isvalid(s) = ~isempty(ind);
        end

        yconf_all = test_boosted_dt_mc(labelclassifier, data{f, k}(isvalid, :));
        if exist('segclassifier','var')
        sconf_all = test_boosted_dt_mc(segclassifier, data{f, k}(isvalid, :));
        end

        isvalid = find(isvalid);
        for s = 1:numel(isvalid)

            ind = find(maps{f}(:, k)==isvalid(s));

            count(ind) = count(ind) + 1;

            yconf = yconf_all(s, :);
            yconf = 1 ./ (1+exp(-yconf));

            if normalize
                yconf = yconf / sum(yconf);
            end

            if exist('segclassifier','var')
            sconf = sconf_all(s, :);
            sconf = 1 ./ (1+exp(-sconf));
            pgs = yconf*sconf;
            else
              pgs = yconf;
            end

            if ~normalize
                sum_s(ind) = sum_s(ind) + sconf;
            end
            pg{f}(ind, :) = pg{f}(ind, :) + repmat(pgs, numel(ind), 1);

        end

        end
    end

    if normalize
        pg{f} = pg{f} ./ max(repmat(sum(pg{f}, 2), 1, size(pg{f}, 2)), 0.00001);
    else
        pg{f} = pg{f} ./ repmat(sum_s, [1 size(pg{f}, 2)]);
    end

    allcount = allcount + mean(count);

end

allcount = allcount / numel(imsegs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [segs, ind] = checksegs(segs, map, s)
% Checks whether this segment has been seen before

ind = find(map==s);

if 0

if numel(ind)==1 % already accounted for by superpixels
    ind = [];
    return;
end

end

oldsegs = segs{ind(1)};

for k = 1:numel(oldsegs)
    if (numel(oldsegs{k})==numel(ind)) && oldsegs{k}(1)==ind(1) && all(oldsegs{k}==ind)
        ind = [];
        return;
    end
end

segs{ind(1)}{end+1} = ind;

author: zouchuhang

repo: LayoutNet


text:
function varargout = drawCuboid(cuboid, varargin)
%DRAWCUBOID Draw a 3D cuboid, eventually rotated.
%
%   drawCuboid(CUBOID)
%   Displays a 3D cuboid on current axis. CUBOID is given by:
%   [XC YC ZC L W D YAW PITCH ROLL],
%   where (XC, YC, ZC) is the cuboid center, L, W and H are the lengths of
%   the cuboid main axes, and YAW PITCH ROLL are Euler angles representing
%   the cuboid orientation, in degrees.
%
%   If cuboid is axis-aligned, it can be specified using only center and
%   side lengths:
%   CUBOID = [XC YC ZC L W H]
%
%   Example
%   % Draw a basic rotated cuboid
%     figure; hold on;
%     drawCuboid([10 20 30   90 40 10   10 20 30], 'FaceColor', 'g');
%     axis equal;
%     view(3);
%
%     % Draw three "borromean" cuboids
%     figure; hold on;
%     drawCuboid([10 20 30 90 50 10], 'FaceColor', 'r');
%     drawCuboid([10 20 30 50 10 90], 'FaceColor', 'g');
%     drawCuboid([10 20 30 10 90 50], 'FaceColor', 'b');
%     view(3); axis equal;
%     set(gcf, 'renderer', 'opengl')
%
%   See also
%   meshes3d, polyhedra, createCube, drawEllipsoid, drawCube
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

phi   = 0;
theta = 0;
psi   = 0;

%% Parses the input
if nargin == 0
    % no input: assumes cuboid with default shape
    xc = 0;	yc = 0; zc = 0;
    a = 5; b = 4; c = 3;

else
    % one argument: parses elements
    xc  = cuboid(:,1);
    yc  = cuboid(:,2);
    zc  = cuboid(:,3);
    a   = cuboid(:,4);
    b   = cuboid(:,5);
    c   = cuboid(:,6);
    if size(cuboid, 2) >= 9
        k   = pi / 180;
        phi   = cuboid(:,7) * k;
        theta = cuboid(:,8) * k;
        psi   = cuboid(:,9) * k;
    end
end

%% Compute cuboid coordinates

% create unit centered cube
[v, f] = createCube;
v = bsxfun(@minus, v, mean(v, 1));

% convert unit basis to ellipsoid basis
sca     = createScaling3d(a, b, c);
rotZ    = createRotationOz(phi);
rotY    = createRotationOy(theta);
rotX    = createRotationOx(psi);
tra     = createTranslation3d([xc yc zc]);

% concatenate transforms
trans   = tra * rotZ * rotY * rotX * sca;

% transform mesh vertices
[x, y, z] = transformPoint3d(v, trans);

%% Process output
if nargout == 0
    % no output: draw the cuboid
    drawMesh([x y z], f, varargin{:});

elseif nargout == 1
    % one output: draw the cuboid and return handle
    varargout{1} = drawMesh([x y z], f, varargin{:});

elseif nargout == 3
    % 3 outputs: return computed coordinates
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end


text:
function G = grad_intrinsic(l,F)
  % GRAD_INTRINSIC Construct an intrinsic gradient operator.
  %
  % Inputs:
  %  l  #F by 3 list of edge lengths
  %  F  #F by 3 list of triangle indices into some vertex list V
  % Outputs:
  %  G  #F*2 by #V gradient matrix: G=[Gx;Gy] where x runs along the 23 edge and
  %    y runs in the counter-clockwise 90° rotation.
  %

  x = (l(:,2).^2-l(:,1).^2-l(:,3).^2)./(-2.*l(:,1));
  y = sqrt(l(:,3).^2 - x.^2);
  n = max(F(:));
  m = size(F,1);
  Z = zeros(m,1);
  V2 = [x y;Z Z;l(:,1) Z];
  F2 = reshape(1:3*m,m,3);
  G2 = grad(V2,F2);
  P = sparse(F2,F,1,m*3,n);
  G = G2*P;

end

author: alecjacobson

repo: gptoolbox


text:
function points3d = mytriangualation(matchedPoints1,matchedPoints2,cam1, cam2)
points2d(:,:,1)=matchedPoints1;
points2d(:,:,2)=matchedPoints2;
cam(:,:,1)=cam1;
cam(:,:,2)=cam2;
nPoints = size(points2d, 1);
points3d = zeros(nPoints, 3, 'like', points2d);

for i = 1:nPoints
    pairs=squeeze(points2d(i, :, :))';
    A = zeros(4, 4);
    for j = 1:2
        P = cam(:,:,j)';
        A(2*j-1,:)=pairs(j, 1)*P(3,:)-P(1,:);
        A(2*j,:)=pairs(j, 2)*P(3,:)-P(2,:);
    end
    [~,~,V] = svd(A);
    X = V(:, end);
    X = X/X(end);
    points3d(i, :) = X(1:3)';
end

author: yihui-he


text:
function [y]=col(x, k, varargin)
%Computes the columns of a block TT-tensor
%   [Y]=COL(X,K)
%   It is useful for block algorithms when r(d+1) is not equal to 1
%   For analogous thing on r(1) see ROW(X,K)
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

d=x.d;
r=x.r;

if nargin > 2
    for j = 1 : nargin-2
        k = [k, varargin{j}];
    end
end

BlockSize = r(d+1);
if max(k) > BlockSize
    error('TT-tensor.col:Index exceeds dimensions in TT block.');
end

CutMatr = eye(BlockSize);
CutMatr = CutMatr(:, k);
y = x*CutMatr;

return
end

author: oseledets

repo: TT-Toolbox


text:
function [mGrid]=calc_SynRateGrid(nNeqM4,fMinLon,fMaxLon,fMinLat,fMaxLat,dfLon,dfLat,nMode,mCat)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: [mGrid] = calc_SynRateGrid(12,-118.2,-114.5,32.1,36.1,0.25,0.25,2,a);
%
% This function (calc_SynRateGrid) calculates the rates on the grid for the
% background rate used for calc_SynCat (program package ETESProject by karen
% Felzer was used.)
%
% Author: van Stiphout, Thomas
% Email: vanstiphout@sed.ethz.ch
% Created: 12. Nov 2007
% Changed: -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Variables
% Input:
% nNeqM4         No. of Events with greater M>=4 eqs/year
% fMinLon    Min Longitude of rate grid
% fMaxLon    Max Longitude of rate grid
% fMinLat    Min Latitude of rate grid
% fMaxLat    Max Latitiude of rate grid
% dfLon      discretization for Longitude grid
% dfLat      discretization for Latitude grid
% nMode      rate grid mode 1: uniform rate, 2: rate according to real data
%              3: probability of an earthquake to occure in this polygon
%
% Output:
% SynCat Simulated ETES catalog [yr lon lat month day mag depth hrs min sec]
% vMain            Vector with assignments if main (=1) or aftershock (=0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vLon=[fMinLon:dfLon:fMaxLon]';
vLat=[fMinLat:dfLat:fMaxLat]';

for i=1:size(vLon,1)-1
    for j=1:size(vLat,1)-1
        mGrid((i-1)*(size(vLat,1)-1)+j,1:4)=[vLon(i),vLon(i+1), vLat(j),vLat(j+1)];
    end
end

% switch btw different rategrids (uniform (nMode = 1)  according to real
% data (n Mode = 2)
switch nMode
    case 1
        mGrid(:,5)=nNeqM4/size(mGrid,1);
    case 2
        nNeqM4=sum(mCat(:,6)>=4)/(max(mCat(:,3))-min(mCat(:,3)))
        for i=1:size(mGrid)
            mGrid_sum(i)= sum(mCat(:,1)>=mGrid(i,1) & mCat(:,1)<mGrid(i,2) & ...
                mCat(:,2)>=mGrid(i,3) & mCat(:,2)<mGrid(i,4));
        end
        % normalize values in mGrid_sum
        mGridNorm=mGrid_sum/sum(mGrid_sum);
        mGrid(:,5)=nNeqM4*mGridNorm;
    case 3
%         nNeqM4=sum(mCat(:,6)>=4)/(max(mCat(:,3))-min(mCat(:,3)))
        for i=1:size(mGrid)
            mGrid_sum(i)= sum(mCat(:,1)>=mGrid(i,1) & mCat(:,1)<mGrid(i,2) & ...
                mCat(:,2)>=mGrid(i,3) & mCat(:,2)<mGrid(i,4));
        end
        % normalize values in mGrid_sum
        mGridNorm=mGrid_sum/max(mGrid_sum);
        mGrid(:,5)=mGridNorm;
end

author: CelsoReyes


text:
function F=msst(X,y,N)
%+++ Multi-Scale significance test
if nargin<3;N=1000;end

[n,p]=size(X);
r0=15/p;
ratio=[0.5 0.6 0.7 0.8 0.9 0.95];
for i=1:N
  for j=1:length(ratio)

    [Xcal,ycal]=traintestselect(X,y,ratio(j));
    kp=find(ycal==1);
    kn=find(ycal~=1);
    for k=1:p
%      [ptemp,h] = ranksum(Xcal(kp,k),Xcal(kn,k));
     [h,ptemp] = ttest2(Xcal(kp,k),Xcal(kn,k));
     P(i,j,k)=-log10(ptemp);
    end

  end
  fprintf('The %ith sampling finished.\n',i);
end
%+++ Output
F.P=P;

author: viggin

repo: domain-adaptation-toolbox


text:
function [z,output] = nlsb_gndl(F,dF,lb,ub,z0,options)
%NLSB_GNDL Bound-constrained NLS by projected Gauss-Newton dogleg TR.
%   [z,output] = nlsb_gndl(F,dF,lb,ub,z0) starts at z0 and attempts to find
%   a local minimizer of the real-valued function f(z), which is the
%   nonlinear least squares objective function f(z) := 0.5*(F(z)'*F(z))
%   subject to the constraints real(lb) <= real(z) <= real(ub) and
%   imag(lb) <= imag(z) <= imag(ub). The input variable z (and lb, ub) may
%   be a scalar, vector, matrix, tensor or even a (nested) cell array of
%   tensors and its contents may be real or complex. This method may be
%   applied in the following ways:
%
%   1. F is function of both z and conj(z).
%
%      Method 1: general medium-scale problems.
%      nlsb_gndl(F,dF,lb,ub,z0) where F(z) returns a column vector of
%      complex. Set dF equal to the string 'Jacobian-C' for automatic
%      numerical approximation of the complex Jacobian, or supply the
%      complex Jacobian manually with a structure dF containing:
%
%         dF.dzc     - The function dF.dzc(zk) should return the complex
%                      Jacobian [dF(zk)/d(z^T) dF(zk)/d(conj(z)^T)], which
%                      is defined as the matrix in which the m-th row is
%                      equal to [(dFm(zk)/dz); (dFm(zk)/d(conj(z)))]^T,
%                      where Fm is the m-th component of F.
%
%      Method 2: general large-scale problems.
%      nlsb_gndl(F,dF,lb,ub,z0) where F(z) returns a column vector of
%      complex residuals and dF is a structure containing:
%
%         dF.dzx     - The function dF.dzx(zk,x,'notransp') should return
%                      the matrix-vector product [dF(zk)/d(z^T)]*x and
%                      dF.dzx(zk,x,'transp') should return the
%                      matrix-vector product [dF(zk)/d(z^T)]'*x.
%         dF.dconjzx - The function dF.dconjzx(zk,x,'notransp') should
%                      return the matrix-vector product
%                      [dF(zk)/d(conj(z)^T)]*x and
%                      dF.dconjzx(zk,x,'transp') should return the matrix-
%                      vector product [dF(zk)/d(conj(z)^T)]'*x.
%
%   2. F is function only of z.
%
%      Method 1: analytic medium-scale problems.
%      nlsb_gndl(F,dF,lb,ub,z0) where F(z) returns a column vector of
%      complex residuals. Set dF equal to the string 'Jacobian' for
%      automatic numerical approximation of the Jacobian, respectively. Or,
%      supply the Jacobian manually with a structure dF containing:
%
%         dF.dz      - The function dF.dz(zk) should return the Jacobian
%                      dF(zk)/d(z^T), which is defined as the matrix in
%                      which the m-th row is equal to (dFm(zk)/dz)^T, where
%                      Fm is the m-th component of F.
%
%      Method 2: analytic large-scale problems.
%      nlsb_gndl(F,dF,lb,ub,z0) where F(z) returns a column vector of
%      complex residuals and dF is a structure containing:
%
%         dF.dzx     - The function dF.dzx(zk,x,'notransp') should return
%                      the matrix-vector product [dF(zk)/d(z^T)]*x and
%                      dF.dzx(zk,x,'transp') should return the matrix-
%                      vector product [dF(zk)/d(z^T)]'*x.
%
%      Method 3: analytic problems in a modest number of variables z and
%                large number of residuals F(z).
%      nlsb_gndl(f,dF,lb,ub,z0) where f(z) := 0.5*(F(z)'*F(z)) and dF is a
%      structure containing:
%
%         dF.JHF     - The function dF.JHF(zk) should return
%                      [dF(zk)/d(z^T)]'*F(zk), which is also equal to
%                      2*df(zk)/d(conj(z)) = 2*conj(df(zk)/d(z)) if z is
%                      complex, or equal to df(xk)/dx if it is real.
%         dF.JHJ     - The function dF.JHF(zk) should return the Gramian
%                      [dF(zk)/d(z^T)]'*[dF(zk)/d(z^T)].
%
%      Method 4: analytic problems in a large number of variables z and
%                large number of residuals F(z).
%      nlsb_gndl(f,dF,lb,ub,z0) where f(z) := 0.5*(F(z)'*F(z)) and dF is a
%      structure containing:
%
%         dF.JHF     - The function dF.JHF(zk) should return
%                      [dF(zk)/d(z^T)]'*F(zk), which is also equal to
%                      2*df(zk)/d(conj(z)) = 2*conj(df(zk)/d(z)) if z is
%                      complex, or equal to df(xk)/dx if it is real.
%         dF.JHJx    - The function dF.JHF(zk,x) should return the matrix-
%                      vector product ([dF(zk)/d(z^T)]'*[dF(zk)/d(z^T)])*x.
%
%   The structure output returns additional information:
%
%      output.alpha        - The plane search step lengths in every
%                            iteration, if a plane search is selected.
%      output.cgiterations - The number of CG/LSQR iterations to compute
%                            the Gauss-Newton step in every iteration.
%                            (large-scale methods only).
%      output.cgrelres     - The relative residual norm of the computed
%                            Gauss-Newton step (large-scale methods only).
%      output.delta        - The trust region radius at every step attempt.
%      output.fval         - The value of the objective function f in every
%                            iteration.
%      output.info         - The circumstances under which the procedure
%                            terminated:
%                               1: Objective function tolerance reached.
%                               2: Step size tolerance reached.
%                               3: Maximum number of iterations reached.
%      output.infops       - The circumstances under which the plane search
%                            terminated in every iteration.
%      output.iterations   - The number of iterations.
%      output.relfval      - The difference in objective function value
%                            between every two successive iterates,
%                            relativeto its initial value.
%      output.relstep      - The step size relative to the norm of the
%                            current iterate in every iteration.
%      output.rho          - The trustworthiness at every step attempt.
%
%   nlsb_gndl(F,dF,lb,ub,z0,options) may be used to set the following
%   options:
%
%      options.CGMaxIter = 15     - The maximum number of CG/LSQR
%                                   iterations for computing the
%                                   Gauss-Newton step (large-scale methods
%                                   only).
%      options.CGTol = 1e-6       - The tolerance for the CG/LSQR method to
%                                   compute the Gauss-Newton step
%                                   (large-scale methods only).
%      options.Delta =            - The initial trust region radius. If
%      0.3*max(1,norm(z0))          equal to NaN, the initial radius will
%                                   be equal to length of the first
%                                   Gauss-Newton step.
%      options.Display = 1        - Displays the objective function value,
%                                   its difference with the previous
%                                   iterate relative to the first iterate
%                                   and the relative step size each
%                                   options.Display iterations. Set to 0 to
%                                   disable.
%      options.JHasFullRank       - If set to true, the Gauss-Newton step
%      = false                      is computed as a least squares
%                                   solution, if possible. Otherwise, it is
%                                   computed using a more expensive
%                                   pseudo-inverse.
%      options.MaxIter = 200      - The maximum number of iterations.
%      options.PlaneSearch        - The plane search used to minimize the
%      = false                      objective function in the plane spanned
%                                   by the steepest descent direction and
%                                   the Gauss-Newton step. Disables dogleg
%                                   trust region strategy. The method
%                                   should have the function signature
%                                   options.PlaneSearch(F,dF,z,p1,p2, ...
%                                   state,options.PlaneSearchOptions).
%      options.PlaneSearchOptions - The options structure passed to the
%                                   plane search search routine.
%      options.TolFun = 1e-12     - The tolerance for output.relfval. Note
%                                   that because the objective function is
%                                   a squared norm, TolFun can be as small
%                                   as eps^2.
%      options.TolX = 1e-6        - The tolerance for output.relstep.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables", SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.
%   [2] C.T. Kelley, "Iterative methods for optimization," SIAM Frontiers
%       in Applied Mathematics, No. 18, 1999.

% Check the objective function f, derivative dF and first iterate z0.
if ~isa(F,'function_handle')
    error('nlsb_gndl:F','The first argument must be a function.');
end
if ischar(dF)
    type = dF;
    if strcmp(type,'Jacobian-C'), fld = 'dzc'; else fld = 'dz'; end
    dF = struct(fld,@derivjac);
end
if ~isstruct(dF)
    error('nlsb_gndl:dF','Second argument not valid.');
else
    if isfield(dF,'dzc')
        method = 'F+dFdzc';
    elseif isfield(dF,'dzx') && isfield(dF,'dconjzx')
        method = 'F+dFdzx+dFdconjzx';
    elseif isfield(dF,'dz')
        method = 'F+dFdz';
    elseif isfield(dF,'dzx')
        method = 'F+dFdzx';
    elseif isfield(dF,'JHJ')  && isfield(dF,'JHF')
        method = 'f+JHJ+JHF';
        f = F;
    elseif isfield(dF,'JHJx') && isfield(dF,'JHF')
        method = 'f+JHJx+JHF';
        f = F;
    else
        error('nlsb_gndl:dF', ...
             ['The structure dF should supply [dF.dzc] or ' ...
              '[dF.dzx and dF.dconjzx] or [dF.dz] or [dF.dzx] or ' ...
              '[dF.JHJ and dF.JHF] or [dF.JHJx and dF.JHF].']);
    end
end

% Define projection and reduction operators.
function z = proj(z)
    z = median([real(lb) real(z) real(ub)],2)+ ...
        median([imag(lb) imag(z) imag(ub)],2)*1i;
end
function H = red(H,A)
    H = bsxfun(@times,H,~A(1:size(H,1)));
    H = bsxfun(@times,H,~A(1:size(H,1)).');
    H(1:size(H,1)+1:end) = H(1:size(H,1)+1:end)+A(1:size(H,1));
end

% Evaluate the function value at z0.
lb = serialize(lb);
ub = serialize(ub);
dim = structure(z0);
z0 = proj(serialize(z0));
A = [real(z0)<=real(lb) | real(z0)>=real(ub) ...
     imag(z0)<=imag(lb) | imag(z0)>=imag(ub)];
z = deserialize(z0,dim);
switch method
    case {'F+dFdzc','F+dFdzx+dFdconjzx','F+dFdz','F+dFdzx'}
        Fval = F(z); Fval = Fval(:);
        fval = 0.5*sum(Fval'*Fval);
    case {'f+JHJ+JHF','f+JHJx+JHF'}
        fval = f(z);
end

% Numerical approximaton of complex derivatives.
function J = derivjac(zk)
    J = deriv(F,zk,Fval,type);
end

% In the case 'F+dFdzx+dFdconjzx', compute a reduced version of J'*(J*x).
function y = JH_Jx(x)
    x1 = ~A(:,1).*x(1:end/2)+~A(:,2).*x(end/2+1:end)*1i;
    dFdzx = dF.dzx(z,x1,'notransp');
    dFdconjzconjx = dF.dconjzx(z,conj(x1),'notransp');
    y = real(dFdzx)+real(dFdconjzconjx)+ ...
        (imag(dFdzx)+imag(dFdconjzconjx))*1i;
	dFdzx = dF.dzx(z,y,'transp');
    dFdconjzconjx = dF.dconjzx(z,y,'transp');
    y = [~A(:,1).*(real(dFdzx)+real(dFdconjzconjx)); ...
         ~A(:,2).*(imag(dFdzx)-imag(dFdconjzconjx))];
    y(A(:)) = x(A(:));
end

% In the case 'F+dFdzx', compute a reduced version of dFdz'*(dFdz*x).
function y = dFdzH_dFdzx(x)
    x1 = ~A(:,1).*x(1:end/2)+~A(:,2).*x(end/2+1:end)*1i;
    y = dF.dzx(z,x1,'notransp');
    y = dF.dzx(z,y,'transp');
    y = [~A(:,1).*real(y);~A(:,2).*imag(y)];
    y(A(:)) = x(A(:));
end

% In the case 'f+JHJx+JHF', compute a reduced version of JHJ*x.
function y = JHJx(x)
    x1 = ~A(:,1).*x(1:end/2)+~A(:,2).*x(end/2+1:end)*1i;
    y = dF.JHJx(z,x1);
    y = [~A(:,1).*real(y);~A(:,2).*imag(y)];
    y(A(:)) = x(A(:));
end

% Modify the preconditioner, if available.
if isfield(dF,'M') && ~isempty(dF.M), dF.PC = @PC; else dF.PC = []; end
function x = PC(b)
    x = dF.M(z,b(1:end/2)+b(end/2+1:end)*1i);
    x = [real(x);imag(x)];
end

% Check the options structure.
if nargin < 6, options = struct; end
if ~isfield(options,'CGMaxIter'), options.CGMaxIter = 15; end
if ~isfield(options,'CGTol'), options.CGTol = 1e-6; end
if ~isfield(options,'Delta'), options.Delta = 0.3*max(1,norm(z0)); end
if ~isfield(options,'Display'), options.Display = 1; end
if ~isfield(options,'JHasFullRank'), options.JHasFullRank = false; end
if ~isfield(options,'MaxIter'), options.MaxIter = 200; end
if ~isfield(options,'PlaneSearch'), options.PlaneSearch = false; end
if ~isfield(options,'PlaneSearchOptions')
    options.PlaneSearchOptions = struct;
end
if ~isfield(options,'TolFun'), options.TolFun = 1e-12; end
if ~isfield(options,'TolX'), options.TolX = 1e-6; end

% Gauss-Newton with dogleg trust region.
output.alpha = [];
output.cgiterations = [];
output.cgrelres = [];
output.delta = options.Delta;
output.fval = fval;
output.info = false;
output.infops = [];
output.iterations = 0;
output.relfval = [];
output.relstep = [];
output.rho = [];
while ~output.info

    % Compute the (in)exact Gauss-Newton step pgn.
    switch method
        case 'F+dFdzc'
            % Compute the Gauss-Newton step pgn.
            dFdzc = dF.dzc(z);
            dFdz = dFdzc(:,1:end/2);
            dFdconjz = dFdzc(:,end/2+1:end);
            J = [real(dFdz)+real(dFdconjz),imag(dFdconjz)-imag(dFdz); ...
                 imag(dFdz)+imag(dFdconjz),real(dFdz)-real(dFdconjz)];
            grad = dFdz'*Fval+dFdconjz.'*conj(Fval);
            if options.JHasFullRank || issparse(J)
                pgn = red(J'*J,A)\[-real(grad);-imag(grad)];
            else
                pgn = pinv(red(J'*J,A))*[-real(grad);-imag(grad)];
            end
            pgn = pgn(1:end/2)+pgn(end/2+1:end)*1i;
            % Compute the Cauchy point pcp = -alpha*grad.
            gg = grad'*grad;
            gBg = dFdz*grad+dFdconjz*conj(grad);
            gBg = gBg'*gBg;
            alpha = gg/gBg;
        case 'F+dFdzx+dFdconjzx'
            % Compute the Cauchy point pcp = -alpha*grad.
            grad = dF.dzx(z,Fval,'transp')+ ...
                   conj(dF.dconjzx(z,Fval,'transp'));
            gg = grad'*grad;
            gBg = dF.dzx(z,grad,'notransp')+ ...
                  dF.dconjzx(z,conj(grad),'notransp');
            gBg = gBg'*gBg;
            alpha = gg/gBg;
            if ~isfinite(alpha), alpha = 1; end;
            % Compute the Gauss-Newton step pgn.
            [pgn,~,output.cgrelres(end+1),output.cgiterations(end+1)] = ...
                mpcg(@JH_Jx,[-real(grad);-imag(grad)], ...
                     options.CGTol,options.CGMaxIter,dF.PC,[], ...
                     -alpha*[real(grad);imag(grad)]);
            pgn = pgn(1:end/2)+pgn(end/2+1:end)*1i;
        case 'F+dFdz'
            % Compute the Gauss-Newton step pgn.
            dFdz = dF.dz(z);
            grad = dFdz'*Fval;
            JHJr = dFdz'*dFdz;
            gradr = -grad;
            JHJIsReal = isreal(JHJr);
            if ~JHJIsReal
                JHJr = [real(JHJr) -imag(JHJr); imag(JHJr) real(JHJr)];
                gradr = [-real(grad);-imag(grad)];
            end
            if options.JHasFullRank || issparse(dFdz)
                pgn = red(JHJr,A)\gradr;
            else
                pgn = pinv(red(JHJr,A))*gradr;
            end
            if ~JHJIsReal, pgn = pgn(1:end/2)+pgn(end/2+1:end)*1i; end
            % Compute the Cauchy point pcp = -alpha*grad.
            gg = grad'*grad;
            gBg = dFdz*grad;
            gBg = gBg'*gBg;
            alpha = gg/gBg;
        case 'F+dFdzx'
            % Compute the Cauchy point pcp = -alpha*grad.
            grad = dF.dzx(z,Fval,'transp');
            gg = grad'*grad;
            gBg = dF.dzx(z,grad,'notransp');
            gBg = gBg'*gBg;
            alpha = gg/gBg;
            if ~isfinite(alpha), alpha = 1; end;
            % Compute the Gauss-Newton step pgn.
            [pgn,~,output.cgrelres(end+1),output.cgiterations(end+1)] = ...
                mpcg(@dFdzH_dFdzx,-[real(grad);imag(grad)], ...
                     options.CGTol,options.CGMaxIter,dF.PC,[], ...
                     -alpha*[real(grad);imag(grad)]);
            pgn = pgn(1:end/2)+pgn(end/2+1:end)*1i;
        case 'f+JHJ+JHF'
            % Compute the Gauss-Newton step pgn.
            grad = serialize(dF.JHF(z));
            JHJ = dF.JHJ(z);
            JHJr = JHJ;
            gradr = -grad;
            JHJIsReal = isreal(JHJr);
            if ~JHJIsReal
                JHJr = [real(JHJr) -imag(JHJr); imag(JHJr) real(JHJr)];
                gradr = [-real(grad);-imag(grad)];
            end
            if options.JHasFullRank || issparse(JHJ)
                pgn = red(JHJr,A)\gradr;
            else
                pgn = pinv(red(JHJr,A))*gradr;
            end
            if ~JHJIsReal, pgn = pgn(1:end/2)+pgn(end/2+1:end)*1i; end
            % Compute the Cauchy point pcp = -alpha*grad.
            gg = grad'*grad;
            gBg = real(grad'*JHJ*grad);
            alpha = gg/gBg;
        case 'f+JHJx+JHF'
            % Compute the Cauchy point pcp = -alpha*grad.
            grad = serialize(dF.JHF(z));
            gg = grad'*grad;
            gBg = real(grad'*dF.JHJx(z,grad));
            alpha = gg/gBg;
            if ~isfinite(alpha), alpha = 1; end;
            % Compute the Gauss-Newton step pgn.
            [pgn,~,output.cgrelres(end+1),output.cgiterations(end+1)] = ...
                mpcg(@JHJx,-[real(grad);imag(grad)], ...
                     options.CGTol,options.CGMaxIter,dF.PC,[], ...
                     -alpha*[real(grad);imag(grad)]);
            pgn = pgn(1:end/2)+pgn(end/2+1:end)*1i;
    end

    % Project pgn and pcp.
    if ~all(isfinite(pgn)), pgn = -alpha*grad; end
    pgn = proj(z0+pgn)-z0;
    pcp = proj(z0-alpha*grad)-z0;
    gg = pcp'*pcp;

    % Plane search in the plane spanned by {pgn,pcp}.
    if ~all(isfinite(pgn)), pgn = -alpha*grad; end
    if isa(options.PlaneSearch,'function_handle')
        state = output; state.grad = grad;
        [alpha,outputps] = options.PlaneSearch( ...
            F,dF,z,deserialize(pgn,dim),deserialize(pcp,dim), ...
            state,options.PlaneSearchOptions);
        output.alpha(:,end+1) = alpha;
        if length(alpha) < 3, alpha(3) = 1; end
        p = alpha(1)*pgn+alpha(2)*pcp;
        z1 = deserialize(alpha(3)*(z0+p),dim);
        relstep = norm(p)/norm(z0); if isnan(relstep), relstep = 0; end
        if isfield(outputps,'fval'), fval = outputps.fval;
        else
            switch method
                case {'F+dFdzc','F+dFdzx+dFdconjzx','F+dFdz','F+dFdzx'}
                    Fval = F(z); fval = 0.5*sum(Fval(:)'*Fval(:));
                case {'f+JHJ+JHF','f+JHJx+JHF'}
                    fval = f(z);
            end
        end
        if isfield(outputps,'info')
            output.infops(end+1) = outputps.info;
        end
        rho = 1;
    else
        rho = -inf;
    end

    % Dogleg trust region.
    normpgn = norm(pgn);
    if isnan(output.delta(end)), output.delta(end) = max(1,normpgn); end
    while rho <= 0

        % Compute the dogleg step p.
        % Assume the projection did not alter the pgn or pcp too much,
        % computing dfval's would otherwise be quite expensive.
        delta = output.delta(end);
        if normpgn <= delta
            p = pgn;
        elseif sqrt(gg) >= delta
            p = delta/sqrt(gg)*pcp;
        else
            bma = pgn-pcp; bmabma = bma'*bma;
            c = real(pcp'*bma);
            if c <= 0
                beta = (-c+sqrt(c^2+bmabma*(delta^2-gg)))/bmabma;
            else
                beta = (delta^2-gg)/(c+sqrt(c^2+bmabma*(delta^2-gg)));
            end
            p = pcp+beta*bma;
        end

        % Estimate objective function improvement.
        % Because of projection, we can not avoid a matrix-vector product.
        switch method
            case 'F+dFdzc'
                dfval = dFdz*p+dFdconjz*conj(p);
                dfval = -real(p'*grad)-0.5*(dfval'*dfval);
            case 'F+dFdzx+dFdconjzx'
                dfval = dF.dzx(z,p,'notransp')+ ...
                        dF.dconjzx(z,conj(p),'notransp');
                dfval = -real(p'*grad)-0.5*(dfval'*dfval);
            case 'F+dFdz'
                dfval = dFdz*p;
                dfval = -real(p'*grad)-0.5*(dfval'*dfval);
            case 'F+dFdzx'
                dfval = dF.dzx(z,p,'notransp');
                dfval = -real(p'*grad)-0.5*(dfval'*dfval);
            case 'f+JHJ+JHF'
                dfval = -real(p'*grad)-0.5*real(p'*JHJ*p);
            case 'f+JHJx+JHF'
                dfval = -real(p'*grad)-0.5*real(p'*dF.JHJx(z,p));
        end

        % Compute the trustworthiness rho.
        if dfval > 0
            z1 = deserialize(z0+p,dim);
            switch method
                case {'F+dFdzc','F+dFdzx+dFdconjzx','F+dFdz','F+dFdzx'}
                    Fval = F(z1); Fval = Fval(:);
                    fval = 0.5*sum(Fval'*Fval);
                case {'f+JHJ+JHF','f+JHJx+JHF'}
                    fval = f(z1);
            end
            rho = (output.fval(end)-fval)/dfval;
            if isnan(rho), rho = -inf; end
            output.rho(end+1) = rho;
        end

        % Update trust region radius delta.
        if rho > 0.5
            output.delta(end+1) = max(delta,2*norm(p));
        else
            sigma = (1-0.25)/(1+exp(-14*(rho-0.25)))+0.25;
            if normpgn < sigma*delta && rho < 0
                e = ceil(log2(normpgn/delta)/log2(sigma));
                output.delta(end+1) = sigma^e*delta;
            else
                output.delta(end+1) = sigma*delta;
            end
        end

        % Check for convergence.
        relstep = norm(p)/norm(z0); if isnan(relstep), relstep = 0; end
        if rho <= 0 && relstep <= options.TolX
            output.rho(end+1) = rho;
            fval = output.fval(end);
            break;
        end

    end

    % Save current state.
    if rho > 0
        z = z1;
        z0 = serialize(z);
        A = [real(z0)<=real(lb) | real(z0)>=real(ub) ...
             imag(z0)<=imag(lb) | imag(z0)>=imag(ub)];
    end

    % Update the output structure.
    output.fval(end+1) = fval;
    output.iterations = output.iterations+1;
    output.relfval(end+1) = ...
        abs(diff(output.fval(end:-1:end-1)))/abs(output.fval(1));
    output.relstep(end+1) = relstep;
    if output.relfval(end) <= options.TolFun, output.info = 1; end
    if output.relstep(end) <= options.TolX, output.info = 2; end
    if output.iterations >= options.MaxIter, output.info = 3; end

    % Display progress.
    if options.Display > 0 && (output.iterations == 1 || output.info || ...
       mod(output.iterations,options.Display) == 0)
        if output.iterations == 1
            bold = '%s';
            [~,~,~,~,v] = regexp(version('-release'),'([0-9]+)([ab])');
            if usejava('Desktop') && str2double(v{1}{1}) > 2011 || ...
               (str2double(v{1}{1}) == 2011 && strcmpi(v{1}{2},'b'))
                bold = '<strong>%s</strong>';
            end
        end
        if output.iterations == 1 || ...
           mod(output.iterations,15*options.Display) == 0
            fprintf('\n%7s%s','',sprintf(bold,'fval'));
            fprintf('%13s%s','',sprintf(bold,'relfval'));
            fprintf('%10s%s','',sprintf(bold,'relstep'));
            if isa(options.PlaneSearch,'function_handle')
                fprintf('%10s%s','',sprintf(bold,'alpha'));
            else
                fprintf('%10s%s','',sprintf(bold,'delta'));
                fprintf('%8s%s','',sprintf(bold,'rho'));
            end
            fprintf('\n%21s%9s = %4.e %6s = %4.e\n\n','=1/2*norm(F)^2', ...
                    'TolFun',options.TolFun,'TolX',options.TolX);
        end
        if output.iterations == 1
            fprintf('%4i: % 14.8e |\n',0,output.fval(1));
        end
        if isa(options.PlaneSearch,'function_handle')
            stralpha = [repmat('%10.4e ',1,size(output.alpha,1)) '\n'];
            fprintf(['%4i: % 14.8e | %14.8e | %14.8e | ' stralpha], ...
                    output.iterations,output.fval(end), ...
                    output.relfval(end),output.relstep(end), ...
                    abs(output.alpha(:,end)));
        else
            fprintf(['%4i: % 14.8e | %14.8e | %14.8e | ' ...
                     '%10.4e | %10.4e\n'],...
                    output.iterations,output.fval(end), ...
                    output.relfval(end),output.relstep(end), ...
                    output.delta(end),output.rho(end));
        end
    end

end

% Display termination message.
if options.Display > 0
    ahref = '\n%s\n\n';
    x = round(linspace(0,output.iterations,min(500,output.iterations)));
    if length(bold) > 2
        ahref = sprintf(['\n<a href="matlab:semilogy(%s,%s);' ...
            'xlabel(''iteration'');legend(''fval'',' ...
            '''relfval'',''relstep'')">%%s</a>\n\n'],mat2str(x'), ...
            mat2str([output.fval(x+1)' [nan output.relfval(x(2:end))]' ...
                    [nan output.relstep(x(2:end))]'],3));
    end
    switch output.info
        case 1, fprintf(ahref,'Objective function tolerance reached.');
        case 2, fprintf(ahref,'Step size tolerance reached.');
        case 3, fprintf(ahref,'Maximum number of iterations reached.');
    end
end

end

function [z,offset] = deserialize(z,dim,offset)
    if iscell(dim)
        v = z;
        z = cell(size(dim));
        if nargin < 3, offset = 0; end
        for i = 1:numel(z)
            if iscell(dim{i})
                [z{i},offset] = deserialize(v,dim{i},offset);
            else
                n = prod(dim{i}(:));
                z{i} = reshape(v(offset+(1:n)),dim{i});
                offset = offset+n;
            end
        end
    elseif ~isempty(dim)
        z = reshape(z,dim);
    end
end

function z = serialize(z)
    if iscell(z)
        for i = find(cellfun(@iscell,z(:).'))
            z{i} = serialize(z{i});
        end
        s = cellfun(@numel,z(:)); o = [0; cumsum(s)];
        c = z; z = zeros(o(end),1);
        for i = 1:length(s), z(o(i)+(1:s(i))) = c{i}(:); end
    else
        z = z(:);
    end
end

function dim = structure(z)
    if iscell(z)
        dim = cellfun(@size,z,'UniformOutput',false);
        for i = find(cellfun(@iscell,z(:).'))
            dim{i} = structure(z{i});
        end
    else
        dim = size(z);
        if numel(z) == dim(1), dim = []; end
    end
end

author: andrewssobral

repo: lrslibrary


text:
classdef FEMInputWriter < handle

    properties (Access = private)
        fileName
        xmax
        ymax
        Nx
        Ny
        P
        xmesh
        ymesh
        zmesh
        mesh
        nDirichlet
        nNeumann
        DoF
        dispPrescribed
        pointLoads
        edgeNodes
    end

    methods (Access = public)
        function obj = FEMInputWriter(cParams)
            obj.init(cParams);
        end

        function createTest(obj)
            obj.computeMeshGrid();
            obj.createMesh();
            obj.computeBoundaryConditions();
            obj.computePrescribedDisplacementsMatrix();
            obj.computePrescribedPointLoads();
            obj.computeEdgeNodes();
            obj.writeFile();
        end
    end

    methods (Access = private)
        function init(obj,cParams)
            obj.fileName = cParams.testName;
            obj.xmax     = cParams.x1;
            obj.ymax     = cParams.y1;
            obj.Nx       = cParams.N;
            obj.Ny       = cParams.M;
            obj.P        = cParams.P;
            obj.DoF      = cParams.DoF;
        end

        function computeMeshGrid(obj)
            x1        = linspace(0,obj.xmax,obj.Nx);
            x2        = linspace(0,obj.ymax,obj.Ny);
            [X,Y]     = meshgrid(x1,x2);
            Z         = zeros(size(X));
            obj.xmesh = X;
            obj.ymesh = Y;
            obj.zmesh = Z;
        end

        function createMesh(obj)
            [F,V]    = mesh2tri(obj.xmesh,obj.ymesh,obj.zmesh,'f');
            s.coord  = V(:,1:2);
            s.connec = F;
            obj.mesh = Mesh(s);
            obj.mesh.plot;
        end

        function computeBoundaryConditions(obj)
            t              = 0.3*obj.ymax;
            m              = obj.mesh;
            root           = m.coord(:,1) == 0;
            tipLength      = m.coord(:,1) == obj.xmax;
            tipWidth       = m.coord(:,2) < obj.ymax-t & m.coord(:,2) > t;
            tip            = tipLength & tipWidth;
            obj.nDirichlet = find(root);
            obj.nNeumann   = find(tip);
        end

        function computePrescribedDisplacementsMatrix(obj)
            obj.dispPrescribed = obj.computeBoundaryConditionMatrix(obj.DoF,obj.nDirichlet);
        end

        function computePrescribedPointLoads(obj)
            Fmat  = obj.computeBoundaryConditionMatrix(obj.DoF,obj.nNeumann);
            nnode = size(Fmat,1)/obj.DoF;
            Pnod  = obj.P/nnode;
            for i = 2:2:size(Fmat,1)
                Fmat(i,3) = Pnod;
            end
            obj.pointLoads = Fmat;
        end

        function computeEdgeNodes(obj)
            COOR          = obj.mesh.coord;
            x             = COOR(:,1);
            y             = COOR(:,2);
            edge          = find(x==0 | x==max(x) | y==0 | y==max(y));
            obj.edgeNodes = sort(edge);
        end

        function writeFile(obj)
            fileID = fopen([obj.fileName],'w');
            obj.writeProblemData(fileID);
            obj.writeCoordinates(fileID);
            obj.writeConnectivities(fileID);
            obj.writeDirichletData(fileID);
            obj.writeNeumannData(fileID);
            obj.writeNodesSolid(fileID);
            obj.writeExternalBorderNodes(fileID);
            fclose(fileID);
        end

        function writeCoordinates(obj,fileID)
            COOR        = zeros(size(obj.mesh.coord,1),4);
            COOR(:,2:3) = obj.mesh.coord;
            COOR(:,1)   = (1:1:size(COOR,1))';
            fprintf(fileID,'%%%% Coordinates\n%% Node\n');
            fprintf(fileID,'gidcoord = [\n');
            for i=1:size(COOR,1)
                fprintf(fileID,'%d %f %f %f;\n',COOR(i,:));
            end
            fprintf(fileID,'];\n');
        end

        function writeConnectivities(obj,fileID)
            Tnod        = zeros(size(obj.mesh.connec,1),5);
            Tnod(:,2:4) = obj.mesh.connec;
            Tnod(:,1)   = (1:1:size(Tnod,1))';
            fprintf(fileID,'%%%% Connectivities\n%% Node\n');
            fprintf(fileID,'gidlnods = [\n');
            for i=1:size(Tnod,1)
                fprintf(fileID,'%d %d %d %d %d;\n',Tnod(i,:));
            end
            fprintf(fileID,'];\n');
        end

        function writeDirichletData(obj,fileID)
            fprintf(fileID,'%%%% Variable prescribed\n%% Node\n');
            fprintf(fileID,'lnodes = [\n');
            for i=1:size(obj.dispPrescribed,1)
                fprintf(fileID,'%d %d %f;\n',obj.dispPrescribed(i,:));
            end
            fprintf(fileID,'];\n');
        end

        function writeNeumannData(obj,fileID)
            fprintf(fileID,'%%%% Point loads\n%% Node\n');
            fprintf(fileID,'pointload_complete = [\n');
            for i=1:size(obj.pointLoads,1)
                fprintf(fileID,'%d %d %f;\n',obj.pointLoads(i,:));
            end
            fprintf(fileID,'];\n');
        end

        function writeExternalBorderNodes(obj,fileID)
            fprintf(fileID,'%%%% External Border Nodes\n%% Node\n');
            fprintf(fileID,'External_border_nodes = [\n');
            for i=1:size(obj.edgeNodes,1)
                fprintf(fileID,'%d;\n',obj.edgeNodes(i));
            end
            fprintf(fileID,'];\n');
        end
    end

    methods (Access = private, Static)
        function bc = computeBoundaryConditionMatrix(DoF,n)
            bc = zeros(DoF*length(n),3);
            for i = 1:length(n)
                bc(2*i-1,1:2) = [n(i),1];
                bc(2*i,1:2)   = [n(i),2];
            end
        end
        function writeProblemData(fileID)
            fprintf(fileID,'%%%% Data\nData_prb = {\n');
            fprintf(fileID,'''TRIANGLE'';\n''SI'';\n''2D'';\n''Plane_Stress'';\n');
            fprintf(fileID,'''ELASTIC'';\n''MACRO'';\n};\n');
        end
        function writeNodesSolid(fileID)
            fprintf(fileID,'%%%% Nodes solid\n%% Node\n');
            fprintf(fileID,'nodesolid = unique(pointload_complete(:,1));\n');
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function out = setfisx(fis, arg1, arg2, arg3, arg4, arg5, arg6)
%SETFISX Set fuzzy inference system properties.
%
%   See setfis for syntax and explanation.
%
%   It is almost clone of MATALAB's setfis
%   with some modifications, so it's compatible
%   with extendent fuzzy rule structure.
%
%   Compare also lines 90-97 of this function
%   with lines 97-106 of the original setfis.

%   Per Konstantin A. Sidelnikov, 2009.

numInputs=length(fis.input);
numOutputs=length(fis.output);

switch nargin
case 1
   indent=32*ones(1,8);
   fprintf('      Name\n');
   fprintf('      Type\n');
   fprintf('      AndMethod\n');
   fprintf('      OrMethod\n');
   fprintf('      ImpMethod\n');
   fprintf('      AggMethod\n');
   fprintf('      DefuzzMethod\n');
   fprintf('      InMFParams\n');
   fprintf('      OutMFParams\n');
   fprintf('      RuleList\n');

case 3
   propName=lower(arg1);
   newVal=arg2;

   switch propName
   case 'name'
      fis.name=newVal;
      out=fis;
   case 'type'
      fis.type=newVal;
      out=fis;
   case 'andmethod'
      fis.andMethod=newVal;
      out=fis;
   case 'ormethod'
      fis.orMethod=newVal;
      out=fis;
   case 'impmethod'
      fis.impMethod=newVal;
      out=fis;
   case 'aggmethod'
      fis.aggMethod=newVal;
      out=fis;
   case 'defuzzmethod'
      fis.defuzzMethod=newVal;
      out=fis;

   case {'inlabels','outlabels','inmflabels','outmflabels','inrange','outrange','inmftypes','outmftypes'}
      error('You may not set this property directly');

   case 'inmfparams'
      for i=1:numInputs
         numInputMFs(i)=length(fis.input(i).mf);
      end
      totalInputMFs=sum(numInputMFs);
      k=1;
      for i=1:numInputs
         for j=1:numInputMFs(i)
            fis.input(i).mf(j).params=newVal(k,:);
            k=k+1;
         end
      end
      out=fis;

   case 'outmfparams'
      for i=1:length(fis.output)
         numOutputMFs(i)=length(fis.output(i).mf);
      end
      totalOutputMFs=sum(numOutputMFs);
      k=1;
      for i=1:numOutputs
         for j=1:numOutputMFs(i)
            fis.output(i).mf(j).params=newVal(k,:);
            k=k+1;
         end
      end
      out=fis;

   case 'rulelist'
       newVal = newVal';
       fis.rule = struct( ...
           'antecedent', newVal(1, :), ...
           'consequent', newVal(2, :), ...
           'weight', newVal(3, :), ...
           'connection', newVal(4, :));
       out = fis;

   otherwise
      error(['There is no FIS system property called ', propName]);
   end

case 5
   % Name assignment
   % ===========================================
   varType=lower(arg1);
   varIndex=arg2;
   varProp=lower(arg3);
   newVal=arg4;

   switch varType
   case 'input'

      if varIndex>length(fis.input),
         error(['There are not that many input variables.']);
      end

      switch varProp
      case 'name'
         fis.input(varIndex).name=newVal;
         out=fis;
      case 'range'
         fis.input(varIndex).range=newVal;
         out=fis;
      case 'nummfs'
         error('You may not set this property directly');
      case 'mflist'
         error('You may not set this property directly');
      end

   case 'output'
      % Range checking
      if varIndex>length(fis.output),
         error(['There are not that many output variables.']);
      end

      switch varProp
      case 'name'
         fis.output(varIndex).name=newVal;
         out=fis;
      case 'range'
         fis.output(varIndex).range=newVal;
         out=fis;
      case 'nummfs'
         error('You may not set this property directly');
      case 'mflist'
         error('You may not set this property directly');
      end

   otherwise
      disp(['Variable type must be either "input" or "output"']);

   end

   % ===============================================
   % Handle MEMBERSHIP FUNCTIONS
   % ===============================================
case 7
   % Name assignment
   % ===========================================
   varType=lower(arg1);
   varIndex=arg2;
   MFIndex=arg4;
   MFProp=lower(arg5);
   newVal=arg6;

   % New value preparation
   % ===========================================
   switch varType
   case 'input'

      % Range checking
      % =======================================
      if varIndex>length(fis.input)
         error(['There are not that many input variables.']);
      end

      if MFIndex>length(fis.input(varIndex).mf),
         errStr=['There are only ',int2str(length(fis.input(varIndex).mf)), ...
               ' MFs associated with that variable'];
         error(errStr)
      end

      switch MFProp
      case 'name'
         fis.input(varIndex).mf(MFIndex).name=newVal;
         out=fis;
      case 'type'
         fis.input(varIndex).mf(MFIndex).type=newVal;
         out=fis;
      case 'params'
         fis.input(varIndex).mf(MFIndex).params=newVal;
         out=fis;
      end

   case 'output'
      % Range checking
      % =======================================
      if MFIndex>length(fis.output(varIndex).mf),
         errStr=['There are only ',int2str(length(fis.output(varIndex).mf)), ...
               ' MFs associated with that variable'];
         error(errStr)
      end

      switch MFProp
      case 'name'
         fis.output(varIndex).mf(MFIndex).name=newVal;
         out=fis;
      case 'type'
         fis.output(varIndex).mf(MFIndex).type=newVal;
         out=fis;
      case 'params'
         fis.output(varIndex).mf(MFIndex).params=newVal;
         out=fis;
      end
   end

end


text:
function c=ref_rdgt(f,g,a,M)
%REF_RDGT  Reference Real DGT
%   Usage:  c=ref_rdgt(f,g,a,M);
%
%   Linear algebra version of the algorithm. Create big matrix
%   containing all the basis functions and multiply with the transpose.

L=size(f,1);

b=L/M;
N=L/a;

Mhalf=ceil(M/2);

F=zeros(L,M*N);

l=(0:L-1).';

for n=0:N-1

  % Do the unmodulated coefficient.
  F(:,M*n+1)=circshift(g,n*a);

  for m=1:Mhalf-1
    F(:,M*n+2*m)=sqrt(2)*cos(2*pi*m*l/M).*circshift(g,n*a);

    F(:,M*n+2*m+1)=sqrt(2)*sin(2*pi*m*l/M).*circshift(g,n*a);

  end;

  if mod(M,2)==0
    F(:,M*(n+1))=cos(pi*l).*circshift(g,n*a);
  end;

end;

% dot-transpose will work because F is real.
c=F.'*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function y = gvar(x,n)
%GVAR Variance of a grouped sample.
%  In some scientific works, once the data have been gathered from a
%  population of interest, it is often difficult to get a sense of what
%  the data indicate when they are presented in an unorganized fashion.
%  Assembling the raw data into a meaningful form, such as a frequency
%  distribution, makes the data easier to understand and interpret. It is
%  in the context of frequency distributions that the importance of
%  conveying in a succinct way numerical information contained in the data
%  is encountered.
%  So, grouped data is data that has been organized into groups known as
%  classes.  The raw dataset can be organized by constructing a table
%  showing the frequency distribution of the variable (whose values are
%  given in the raw dataset). Such a frequency table is often referred to
%  as grouped data.
%  Here, we developed a m-code to calculate the variance of a grouped data.
%  One can input the returns or modified vectors n and xout containing the
%  frequency counts and the bin locations of the hist m-function, in a
%  column form matrix.
%  Normalizes Y by (N-1), where N is the sample size.  This is an unbiased
%  estimator of the variance of the population from which X is drawn.
%  Y = VAR(X,1) normalizes by N and produces the second moment of the
%  sample about its mean. VAR(X,0) is the same as VAR(X).
%
%  Variance calculation uses the formula,
%
%
%              V =  Sum(F*(MC - M)^2)/D
%
%  where:
%  F  = class frequency
%  MC = mark class
%  M = grouped mean
%  D = N - 1 or N, whether normalizes by n-1 or by n
%  N = sample size [=sum(F)]
%
%  Syntax: function y = gvar(x)
%
%  Inputs:
%       x - data matrix (Size of matrix must be n-by-2; absolut frequency=
%           column 1, class mark=column 2)
%       n - [normalized by n-1] = 0 (default), [normalized by n] = 1
%  Outputs:
%       y  - variance of the values in x
%
%  Example: Suppose we have the next frequency table:
%
%                     ----------------
%                       MC         F
%                     ----------------
%                        5         6
%                       15        16
%                       25        24
%                       35        25
%                       45        17
%                     ----------------
%
%  Taken from:
%  http://mlsc.lboro.ac.uk/resources/statistics/var_stand_deviat_group.pdf
%
%  Where we are interested to get the variance value normalized by n.
%
%  Data input:
%  x=[5 6;15 16;25 24;35 25;45 17];
%
%  Calling on Matlab the function:
%          y = gvar(x,1)
%
%  Answer is:
%
%  y = 138.7268
%
%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.edu.mx
%
%  Copyright (C)  September 18, 2012.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2012).
%     gvar:Varience of a grouped sample. [WWW document].
%     URL http://www.mathworks.com/matlabcentral/fileexchange/
%     38281-gvar
%
%  Reference:
%  Jayaraman, K. (1999), A Statistical Manual for Foresty Research. Foresty
%             Research Support Programme for Asia and the Pacific. FAO-
%             Corporate Document Repository. Forestry Statistics and Data
%             Collection.
%             URL http://www.fao.org/DOCREP/003/X6831E/X6831E00.HTM
%             PDF ftp://ftp.fao.org/docrep/fao/003/X6831E/X6831E00.pdf
%

if  nargin < 2,
    n = 0; %default
end

c = size(x,2);

if c ~= 2
    error('stats:gvar:BadData','X must have two colums.');
end

mc = x(:,1); %class mark
f = x(:,2); %absolut frequency
s = sum(f.*mc);
m = s/sum(f);
ss = sum(f.*(mc - m).^2);

if n == 0; %default
    d = sum(f) - 1;
else n = 1;
    d = sum(f);
end

y = ss/d;

return,

author: Sable


text:
% StackExchange Signal Processing Q62024
% https://dsp.stackexchange.com/questions/62024
% Proximal Gradient Method (PGM) for a Function Model with More than 2 Functions (Sum of Functions)
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     25/11/2019
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = OFF;

%% Simulation Parameters

numElements = 40;
paramLambda1 = 0.5; %<! L1 Norm
paramLambda2 = 0.75; %<! TV Norm

numIterations = 1000; %<! For the ADMM

%% Generate Data

vY = 10 * randn(numElements, 1);

% Generate the Diff Operator (1D Gradient) by Finite Differences
mD = spdiags([-ones(numElements, 1), ones(numElements, 1)], [0, 1], numElements - 1, numElements);

hSolveProxTv = @(vY, paramLambda) SolveProxTvAdmm(vY, mD, paramLambda, numIterations);

% Objective Function
hObjFun = @(vX) (0.5 * sum( (vX - vY) .^ 2)) + (paramLambda1 * sum(abs(vX))) + (paramLambda2 * sum(abs(mD * vX)));

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable vX(numElements);
    minimize( (0.5 * pow_pos(norm(vX - vY, 2), 2)) + (paramLambda1 * norm(vX, 1)) + (paramLambda2 * norm(mD * vX, 1)));
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp([' ']);

%% Solution by Analytical Solution
%{
Solving $ \arg \min_x \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} +
{\lambda}_{1} {\left\| x \right\|}_{1} + {\lambda}_{2} {\left\| D x \right\|}_{1} $
%}

vX = SolveProxL1(hSolveProxTv(vY, paramLambda2), paramLambda1);

disp([' ']);
disp(['Analytical Direct Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp([' ']);

%% Solution by Analytical Solution
%{
Solving $ \arg \min_x \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} +
{\lambda}_{1} {\left\| x \right\|}_{1} + {\lambda}_{2} {\left\| D x \right\|}_{1} $
%}

vX = hSolveProxTv(SolveProxL1(vY, paramLambda1), paramLambda2);

disp([' ']);
disp(['Analytical Direct Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
global data Ysiz d1 d2 numFrame;

%% select file
if ~exist('nam', 'var') || isempty(nam)
    % choose files manually
    try
        load .dir.mat; %load previous path
    catch
        dir_nm = [cd(), filesep]; %use the current path
    end
    [file_nm, dir_nm] = uigetfile(fullfile(dir_nm, '*.tif;*.mat;*.h5'));
    nam = [dir_nm, file_nm];  % full name of the data file
    [dir_nm, file_nm, file_type] = fileparts(nam);
else
    % use pre-specified file
    if exist(nam, 'file')
        [dir_nm, file_nm, file_type] = fileparts(nam);
    else
        dir_nm = 0;
    end
end
if dir_nm~=0
    save .dir.mat dir_nm;
else
    fprintf('no file was selected. STOP!\n');
    return;
end

%% convert the data to mat file
nam_mat = [dir_nm, filesep, file_nm, '.mat'];
if strcmpi(file_type, '.mat')
    fprintf('The selected file is *.mat file\n');
elseif  exist(nam_mat, 'file')
    % the selected file has been converted to *.mat file already
    fprintf('The selected file has been replaced with its *.mat version.\n');
elseif or(strcmpi(file_type, '.tif'), strcmpi(file_type, '.tiff'))
    % convert
    tic;
    fprintf('converting the selected file to *.mat version...\n');
    nam_mat = tif2mat(nam);
    fprintf('Time cost in converting data to *.mat file:     %.2f seconds\n', toc);
elseif or(strcmpi(file_type, '.h5'), strcmpi(file_type, '.hdf5'))
    fprintf('the selected file is hdf5 file\n');
    temp = h5info(nam);
    dataset_nam = ['/', temp.Datasets.Name];
    dataset_info = h5info(nam, dataset_nam);
    dims = dataset_info.Dataspace.Size;
    ndims = length(dims);
    d1 = dims(2);
    d2 = dims(3);
    numFrame = dims(end);
    Ysiz = [d1, d2, numFrame];
    fprintf('\nThe data has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1, d2, numFrame, prod(Ysiz)*8/(2^30));
    return;
else
    fprintf('The selected file type was not supported yet! email me to get support (zhoupc1988@gmail.com)\n');
    return;
end

%% information of the data
data = matfile(nam_mat);
Ysiz = data.Ysiz;
d1 = Ysiz(1);   %height
d2 = Ysiz(2);   %width
numFrame = Ysiz(3);    %total number of frames

fprintf('\nThe data has been mapped to RAM. It has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1, d2, numFrame, prod(Ysiz)*8/(2^30));

author: zhoupc


text:
function add_display_menu(version)
    % adds a display menu to the current figure
    %
    % appently mostly used by the view_  functions

    % when called, hzma is the current axes;
    ZG=ZmapGlobal.Data;

    op2e = uimenu('Label','Display');
    switch version
        case 1
            circlefun=@plotci2;
            fixscalefun=@(data)fix_caxis(data,'horiz');
            uimenu(op2e,'Label','Plot Map in lambert projection','MenuSelectedFcn',@(~,~)plotmap);
            %overlayfun=@overlay;
        case 2
            circlefun=@plotci2;
            fixscalefun=@(data)fix_caxis(data,'');
            uimenu(op2e,'Label','Plot Map in lambert projection','MenuSelectedFcn',@(~,~)plotmap);
            %overlayfun=@overlay;
        case 3
            circlefun=@plotci3;
            fixscalefun=@(data)fix_caxis(data,'horiz');
            %overlayfun=@()zmap_update_displays();
        case 4
            circlefun=@plotci2;
            fixscalefun=@(data)fix_caxis(data,'horiz');
            uimenu(op2e,'Label','Plot Map in lambert projection','MenuSelectedFcn',@(~,~)plotmap)
            uimenu(op2e,'Label','Plot map on top of topography (white background)',...
                'MenuSelectedFcn', @(~,~) dramap_z('dramap2_z','w', valueMap)); % this is different from case
            uimenu(op2e,'Label','Plot map on top of topography (black background)',...
                'MenuSelectedFcn', @(~,~)dramap_z('dramap2_z','k', valueMap)); % this is different from case #1
            %overlayfun=@overlay;
        case 5
            circlefun=@plotci2;
            fixscalefun=@(data)fix_caxis(data,'horiz');
            add_colormap_section(op2e);
            add_shading_section(op2e);
            add_brighten_section(op2e);
            %overlayfun=@overlay;
    end

    uimenu(op2e,'Label','Fix color (z) scale','MenuSelectedFcn',@(~,~)fixscalefun(ZGvalueMap));
    uimenu(op2e,'Label','Show Grid','MenuSelectedFcn',@cb_showgrid);
    uimenu(op2e,'Label','Show Circles','MenuSelectedFcn',@(~,~)circlefun);
    add_colormap_section(op2e);
    add_shading_section(op2e);
    add_brighten_section(op2e);
    uimenu(op2e,'Label','Redraw Overlay',...
        'MenuSelectedFcn', "set(gca,'NextPlot','add');zmap_update_displays();"); % this is different from case #1

    function cb_shader(style)
        % set default shading style and apply to current axes
        axes(gca);
        ZG.shading_style=style;
        shading(ZG.shading_style);
    end

    function cb_showgrid(src,~)
        set(gca,'NextPlot','add');
        plot(newgri(:,1),newgri(:,2),'+k')
    end
    function add_brighten_section(parent)
        uimenu(parent,'Label','Brighten +0.4','MenuSelectedFcn', @(~,~)brighten(0.4));
        uimenu(parent,'Label','Brighten -0.4','MenuSelectedFcn', @(~,~)brighten(-0.4))
    end
    function add_colormap_section(parent)
        uimenu(parent,'Label','Colormap InvertGray',...
            'MenuSelectedFcn',@(~,~)flip_and_brighten(@gray,0.4));
        uimenu(parent,'Label','Colormap Invertjet',...
            'MenuSelectedFcn',@(~,~)flip_and_brighten(@jet,0));

        function flip_and_brighten(colorfn,brightenamount)
            colormap(flipud(colorfn(64)));
            if brightenamount; brighten(brightenamount); end
        end
    end

    function add_shading_section(parent)
        %TODO make this 1 option, simple inputdlg box, or flip the names
        uimenu(parent,'Label','shading flat',...
            'MenuSelectedFcn',@(~,~)cb_shader('flat'))
        uimenu(parent,'Label','shading interpolated',...
            'MenuSelectedFcn',@(~,~)cb_shader('interp'))
    end

end


text:
function root_num = p15_root_num ( )

%*****************************************************************************80
%
%% P15_ROOT_NUM returns the number of known roots for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer ROOT_NUM, the number of known roots.
%
  root_num = 1;

  return
end

author: johannesgerer


text:
classdef IDTLZ1 < PROBLEM
% <multi/many> <real> <large/none> <expensive/none>
% Inverted DTLZ1

%------------------------------- Reference --------------------------------
% H. Jain and K. Deb, An evolutionary many-objective optimization algorithm
% using reference-point based non-dominated sorting approach, part II:
% Handling constraints and extending to an adaptive approach, IEEE
% Transactions on Evolutionary Computation, 2014, 18(4): 602-622.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M+4; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            g      = 100*(obj.D-obj.M+1+sum((PopDec(:,obj.M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,obj.M:end)-0.5)),2));
            PopObj = (1+repmat(g,1,obj.M))/2 - 0.5*repmat(1+g,1,obj.M).*fliplr(cumprod([ones(size(PopDec,1),1),PopDec(:,1:obj.M-1)],2)).*[ones(size(PopDec,1),1),1-PopDec(:,obj.M-1:-1:1)];
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = (1-UniformPoint(N,obj.M))/2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,1,10)';
                R = {(1-a*a')/2,(1-a*(1-a'))/2,(1-(1-a)*ones(size(a')))/2};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%
% my_append(file, lines)
%
% Similate the function of "double right arrow file writing" function in Unix
% The inputs are
%       file: the destination
%       lines: the text that will be written
%
% Author: Xiong Xiao,
%         School of Computer Engineering, Nanyang Technological University
% Date Created: 02 Feb, 2009
% Last Modified: 02 Feb, 2009
%
function my_append(file, lines, start, stop)

FID = fopen(file, 'a');
if nargin < 3
    start = 1;
    stop = length(lines);
end

if iscell(lines)
    for i=start:min(stop,length(lines))
        if length(lines{i})==0
            continue;
        elseif strcmp(lines{i}(end), '\n') == 1
            fprintf(FID, '%s', lines{i});
        else
            fprintf(FID, '%s\n', lines{i});
        end
    end
else
    fprintf(FID, '%s\n', lines);
end

fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function [EB] = bit2EB(bit)
% Convert computery things from bits to exabytes.
% Chad A. Greene 2012
EB = bit*2^-63;

repo: mcbench-benchmarks