text:
%% Example 1 Overview
%
% GLMsingle is new tool that provides efficient, scalable, and accurate
% single-trial fMRI response estimates.
%
% The purpose of this Example 1 notebook is to guide the user through basic
% calls to GLMsingle, using a representative, small-scale test dataset (in
% this case, an example session from a rapid event-related visual fMRI
% dataset - the Natural Scenes Dataset core experiment).
%
% The goal is to examine the effect of GLMsingle on the reliability of
% single-trial fMRI response estimates. By default, the tool implements a
% set of optimizations that improve upon generic GLM approaches by: (1)
% identifying an optimal hemodynamic response function (HRF) at each voxel,
% (2) deriving a set of useful GLM nuisance regressors via "GLMdenoise" and
% picking an optimal number to include in the final GLM, and (3) applying a
% custom amount of ridge regularization at each voxel using an efficient
% technique called "fracridge". The output of GLMsingle are GLM betas
% reflecting the estimated percent signal change in each voxel in response
% to each experimental stimulus or condition being modeled.
%
% Beyond directly improving the reliability of neural responses to repeated
% stimuli, these optimized techniques for signal estimation can have a
% range of desirable downstream effects such as: improving cross-subject
% representational similarity within and between datasets; improving the
% single-image decodability of evoked neural patterns via MVPA; and,
% decreasing the correlation in spatial patterns observed at neighboring
% timepoints in analysis of fMRI GLM outputs. See our video presentation at
% V-VSS 2020 for a summary of these phenomena as observed in recent
% massive-scale fMRI datasets (the Natural Scenes Dataset and BOLD5000):
% https://www.youtube.com/watch?v=yb3Nn7Han8o
%
% Example 1 contains a full walkthrough of the process of loading an
% example dataset and design matrix, estimating neural responses using
% GLMsingle, estimating the reliability of responses at each voxel, and
% comparing those achieved via GLMsingle to those achieved using a baseline
% GLM. After loading and visualizing formatted fMRI time-series and their
% corresponding design matrices, we will describe the default behavior of
% GLMsingle and show how to modify hyperparameters if the user desires.
% Throughout the notebook we will highlight important metrics and outputs
% using figures, print statements, and comments.
%
% Users encountering bugs, unexpected outputs, or other issues regarding
% GLMsingle shouldn't hesitate to raise an issue on GitHub:
% https://github.com/kendrickkay/GLMsingle/issues

%% Add dependencies and download the example dataset

% Start fresh
clear
clc
close all

this_dir    = fileparts(which('example1.m'));

% Add path to GLMsingle
run(fullfile(this_dir, '..', '..', 'setup.m'));

% Name of directory to which outputs will be saved
outputdir = fullfile(this_dir, 'example1outputs');

% Download files to data directory
input_dir = fullfile(this_dir, 'data');
if ~exist(input_dir, 'dir')
    mkdir('data')
end

input_file = fullfile(input_dir, 'nsdcoreexampledataset.mat');
URL = 'https://osf.io/k89b2/download';

download_data(URL, input_file);

load(input_file)

% Data comes from the NSD dataset (subj01, nsd01 scan session).
% https://www.biorxiv.org/content/10.1101/2021.02.22.432340v1.full.pdf

%% Data overview
clc
whos

% data -> consists of several runs of 4D volume files (x,y,z,t)  where
% (t)ime is the 4th dimention. In this example, data consists of only a
% single slice and has been prepared with a TR = 1s

% ROI -> manually defined region in the occipital cortex. It is a binary
% matrix where (x,y,z) = 1 corresponds to the cortical area that responded
% to visual stimuli used in the NSD project.

fprintf('There are %d runs in total.\n',length(design));
fprintf('The dimensions of the data for the first run are %s.\n',mat2str(size(data{1})));
fprintf('The stimulus duration is %.3f seconds.\n',stimdur);
fprintf('The sampling rate (TR) is %.3f seconds.\n',tr);

%%

figure(1);clf

%Show example design matrix.
for d = 1
    imagesc(design{d}); colormap gray; drawnow
    xlabel('Conditions')
    ylabel('TRs')
    title(sprintf('Design matrix for run%i',d))
    axis image
end

xticks(0:53:length(design{d}))
set(gcf,'Position',[418   412   782   605])

%%
% design -> Each run has a corresponding design matrix where each column
% describes a single condition (conditions are repeated across runs). Each
% design matrix is binary with 1 specfing the time (TR) when the stimulus
% is presented on the screen.
%
% In this NSD scan session, there are a total of 750 trials, in which a
% total of 583 distinct images are shown. (Thus, some images were presented
% more than once.) In the design matrix shown, there are 583 predictor
% columns/conditions, one per distinct image. Notice that white rectangles
% are pseudo randomized and they indicate when the presentation of each
% image occurs. Note that in some runs not all images are shown; if a
% column does not have a white rectangle it means that this image is shown
% in a different run within the session.

%%

% Show an example slice of the first fMRI volume.
figure(2);clf
imagesc(data{1}(:,:,:,1));
colormap(gray);
axis equal tight;
c=colorbar;
title('fMRI data (first volume)');
set(gcf,'Position',[418   412   782   605])
axis off
c.Label.String = 'T2*w intensity';
set(gca,'FontSize',15)

%% Call GLMestimatesingletrial with default parameters

% Outputs and figures will be stored in a folder (you can specify its name
% as the 5th input to GLMestimatesingletrial). Model estimates can be also
% saved to the 'results' variable which is the only output of
% GLMestimatesingletrial.

% Optional parameters below can be assigned to a structure, i.e., opt =
% struct('wantlibrary',1,'wantglmdenoise',1); Options are the 6th input to
% GLMestimatesingletrial.

% There are many options that can be specified; here, we comment on the
% main options that one might want to modify/set. Defaults for the options
% are indicated below.

% wantlibrary = 1 -> Fit HRF to each voxel
% wantglmdenoise = 1 -> Use GLMdenoise
% wantfracridge = 1  -> Use ridge regression to improve beta estimates
% chunknum = 50000 -> is the number of voxels that we will
%     process at the same time. For setups with lower memory, you may need to
%     decrease this number.

% wantmemoryoutputs is a logical vector [A B C D] indicating which of the
%     four model types to return in the output <results>. The user must be
%     careful with this, as large datasets can require a lot of RAM. If you
%     do not request the various model types, they will be cleared from
%     memory (but still potentially saved to disk). Default: [0 0 0 1]
%     which means return only the final type-D model.

% wantfileoutputs is a logical vector [A B C D] indicating which of the
%     four model types to save to disk (assuming that they are computed). A
%     = 0/1 for saving the results of the ONOFF model, B = 0/1 for saving
%     the results of the FITHRF model, C = 0/1 for saving the results of the
%     FITHRF_GLMdenoise model, D = 0/1 for saving the results of the
%     FITHRF_GLMdenoise_RR model. Default: [1 1 1 1] which means save all
%     computed results to disk.

% numpcstotry (optional) is a non-negative integer indicating the maximum
%     number of GLMdenoise PCs to enter into the model. Default: 10.

% fracs (optional) is a vector of fractions that are greater than 0
%     and less than or equal to 1. We automatically sort in descending
%     order and ensure the fractions are unique. These fractions indicate
%     the regularization levels to evaluate using fractional ridge
%     regression (fracridge) and cross-validation. Default:
%     fliplr(.05:.05:1). A special case is when <fracs> is specified as a
%     single scalar value. In this case, cross-validation is NOT performed
%     for the type-D model, and we instead blindly use the supplied
%     fractional value for the type-D model.

% For the purpose of this example, we will keep all outputs in the memory.
opt = struct('wantmemoryoutputs',[1 1 1 1]);

% This example saves output .mat files to the folder
% "example1outputs/GLMsingle". If these outputs don't already exist, we
% will perform the time-consuming call to GLMestimatesingletrial.m;
% otherwise, we will just load from disk.
if ~exist(fullfile(outputdir, 'GLMsingle', 'TYPEB_FITHRF.mat'),'file') || ...
   ~exist(fullfile(outputdir, 'GLMsingle', 'TYPEC_FITHRF_GLMDENOISE.mat'),'file') || ...
   ~exist(fullfile(outputdir, 'GLMsingle', 'TYPED_FITHRF_GLMDENOISE_RR.mat'),'file')

    [results designSINGLE] = GLMestimatesingletrial(design,data,stimdur,tr,[outputdir '/GLMsingle'],opt);

    % We assign outputs of GLMestimatesingletrial to "models" structure.
    % Note that results{1} contains GLM estimates from an ONOFF model,
    % where all images are treated as the same condition. These estimates
    % could be potentially used to find cortical areas that respond to
    % visual stimuli. We want to compare beta weights between conditions
    % therefore we are not going to store the ONOFF GLM results.

    clear models;
    models.FIT_HRF = results{2};
    models.FIT_HRF_GLMdenoise = results{3};
    models.FIT_HRF_GLMdenoise_RR = results{4};

else
    % Load existing file outputs if they exist
    results = load([outputdir '/GLMsingle/TYPEB_FITHRF.mat']);
    models.FIT_HRF = results;
    results = load([outputdir '/GLMsingle/TYPEC_FITHRF_GLMDENOISE.mat']);
    models.FIT_HRF_GLMdenoise = results;
    results = load([outputdir '/GLMsingle/TYPED_FITHRF_GLMDENOISE_RR.mat']);
    models.FIT_HRF_GLMdenoise_RR = results;

end

%% Summary of important outputs

% The outputs of GLMestimatesingletrial.m are formally documented in its
% header. Here, we highlight a few of the more important outputs:
%
% R2 -> is model accuracy expressed in terms of R^2 (percentage).
%
% modelmd -> is the full set of single-trial beta weights (X x Y x Z x
% TRIALS). Beta weights are arranged in chronological order.
%
% HRFindex -> is the 1-index of the best fit HRF. HRFs can be recovered
% with getcanonicalHRFlibrary(stimdur,tr)
%
% FRACvalue -> is the fractional ridge regression regularization level
% chosen for each voxel. Values closer to 1 mean less regularization.

%% Plot a slice of brain showing GLMsingle outputs

% We are going to plot several outputs from the FIT_HRF_GLMdenoise_RR GLM,
% which contains the full set of GLMsingle optimizations.

slice = 1;

% we will plot betas, R2, optimal HRF indices, and the voxel frac values
val2plot = {'modelmd';'R2';'HRFindex';'FRACvalue'};
cmaps = {cmapsign2;hot;jet;copper};

figure(3);clf

for v = 1 : length(val2plot)

    f=subplot(2,2,v);

    if contains('modelmd',val2plot{v})
        % When plotting betas, for simplicity just average across all image
        % presentations This will yield a summary of whether voxels tend to
        % increase or decrease their activity in response to the
        % experimental stimuli (similar to outputs from an ONOFF GLM)
        imagesc(nanmean(models.FIT_HRF_GLMdenoise_RR.(val2plot{v})(:,:,slice),4),[-5 5]); axis off image;
        title('Average GLM betas (750 stimuli)')

    else
        % Plot all other voxel-wise metrics as outputted from GLMsingle
        imagesc(models.FIT_HRF_GLMdenoise_RR.(val2plot{v})(:,:,slice)); axis off image;
        title(val2plot{v})

    end

    colormap(f,cmaps{v})
    colorbar
    set(gca,'FontSize',15)
end

set(gcf,'Position',[418   412   782   605])

%% Run a baseline GLM to compare with GLMsingle

% Additionally, for comparison purposes we are going to run a standard GLM
% without HRF fitting, GLMdenoise, or ridge regression regularization. We
% will change the default settings by using the "opt" structure.
opt.wantlibrary = 0; % switch off HRF fitting
opt.wantglmdenoise = 0; % switch off GLMdenoise
opt.wantfracridge = 0; % switch off ridge regression
opt.wantfileoutputs = [0 1 0 0];
opt.wantmemoryoutputs = [0 1 0 0];

% If these outputs don't already exist, we will perform the call to
% GLMestimatesingletrial.m; otherwise, we will just load from disk.
if ~exist(fullfile(outputdir, 'GLMbaseline', 'TYPEB_FITHRF.mat'),'file')

    [ASSUME_HRF] = GLMestimatesingletrial(design,data,stimdur,tr,[outputdir '/GLMbaseline'],opt);
    models.ASSUME_HRF = ASSUME_HRF{2};

else

    % Note that even though we are loading TYPEB_FITHRF betas, HRF fitting
    % has been turned off and this struct field will thus contain the
    % outputs of a GLM fit using the canonical HRF.
    results = load([outputdir '/GLMbaseline/TYPEB_FITHRF.mat']);
    models.ASSUME_HRF = results;

end

% We assign outputs from GLMestimatesingletrial to "models" structure.
% Again, results{1} contains GLM estimates from an ONOFF model so we are
% not going to extract it.

%%

% Now, "models" variable holds solutions for 4 GLM models
disp(fieldnames(models))

%% Get indices of repeated conditions to use for reliability calculations

% To compare the results of different GLMs we are going to calculate the
% voxel-wise split-half reliablity for each model. Reliability values
% reflect a correlation between beta weights for repeated presentations of
% the same conditions. In short, we are going to check how
% reliable/reproducible are the single trial responses to repeated
% conditions estimated with each GLM type.

% This NSD scan session has a large number of images that are just shown
% once during the session, some images that are shown twice, and a few that
% are shown three times. In the code below, we are attempting to locate the
% indices in the beta weight GLMsingle outputs modelmd(x,y,z,trials) that
% correspond to repeated images. Here we only consider stimuli that have
% been presented at least twice. For the purpose of the example we ignore
% the 3rd repetition of the stimulus.

% Consolidate design matrices
designALL = cat(1,design{:});

% Construct a vector containing 1-indexed condition numbers in
% chronological order.

corder = [];
for p=1:size(designALL,1)
    if any(designALL(p,:))
        corder = [corder find(designALL(p,:))];
    end
end

%%

% Let's take a look at the first few entries
corder(1:3)

% Note that [375 497 8] means that the first stimulus trial involved
% presentation of the 375th condition, the second stimulus trial involved
% presentation of the 497th condition, and so on.

%%

% In order to compute split-half reliability, we have to do some indexing.
% we want to find images with least two repetitions and then prepare a
% useful matrix of indices that refer to when these occur.
repindices = [];  % 2 x images containing stimulus trial indices.

% The first row refers to the first presentation; the second row refers to
% the second presentation.
for p=1:size(designALL,2)  % loop over every condition
    temp = find(corder==p);
    if length(temp) >= 2
        repindices = cat(2,repindices,[temp(1); temp(2)]);  % note that for images with 3 presentations, we are simply ignoring the third trial
    end
end

% Let's take a look at a few entries
repindices(:,1:3)

% Notice that the first condition is presented on the 217th stimulus trial
% and the 486th stimulus trial, the second condition is presented on the
% 218th and 621st stimulus trials, and so on.

fprintf('There are %i repeated images in the experiment \n',length(repindices))

% Now, for each voxel we are going to correlate beta weights describing the
% response to images presented for the first time with beta weights
% describing the response from the repetition of the same image. With 136
% repeated conditions, the correlation for each voxel will reflect the
% relationship between two vectors with 136 beta weights each.

%% Compute median split-half reliability for each GLM version

% Finally, let's compute split-half reliability. We are going to loop
% through our 4 models and calculate split-half reliability for each of
% them.

% We first arrange models from least to most sophisticated (for
% visualization purposes)
model_names = fieldnames(models);
model_names = model_names([4 1 2 3]);

% Create output variable for reliability values
vox_reliabilities = cell(1,length(models));

% For each GLM...
for m = 1 : length(model_names)

    % Get the GLM betas
    betas = models.(model_names{m}).modelmd(:,:,:,repindices);  % use indexing to pull out the trials we want
    betas_reshaped = reshape(betas,size(betas,1),size(betas,2),size(betas,3),2,[]);  % reshape to X x Y x Z x 2 x CONDITIONS

    % compute reliabilities using an efficient (vectorized) utility
    % function
    vox_reliabilities{m} = calccorrelation(betas_reshaped(:,:,:,1,:),betas_reshaped(:,:,:,2,:),5);

    % Note that calccorrelation.m is a utility function that computes
    % correlations in a vectorized fashion (for optimal speed).

end

%% Compare visual voxel reliabilities between beta versions

figure(4);clf
subplot(1,2,1);
cmap = [0.2314    0.6039    0.6980
    0.8615    0.7890    0.2457
    0.8824    0.6863         0
    0.9490    0.1020         0];

% For each GLM type we calculate median reliability for voxels within the
% visual ROI and plot it as a bar plot.
for m = 1 : 4
    bar(m,nanmedian(vox_reliabilities{m}(ROI==1)),'FaceColor','None','Linewidth',3,'EdgeColor',cmap(m,:)); hold on
end
ylabel('Median reliability')
legend(model_names,'Interpreter','None','Location','NorthWest')
set(gca,'Fontsize',16)
set(gca,'TickLabelInterpreter','none')
xtickangle(0)
xticks([])
ylim([0.1 0.2])
set(gcf,'Position',[418   412   782   605])
title('Median voxel split-half reliability of GLM models')

subplot(1,2,1);

% Comparison is the final output (FIT_HRF_GLMDENOISE_RR) vs. the baseline
% GLM (ASSUME_HRF)
vox_reliability = vox_reliabilities{4} - vox_reliabilities{1};
underlay = data{1}(:,:,:,1);
ROI(ROI~=1) = NaN;
overlay = vox_reliability;

underlay_im = cmaplookup(underlay,min(underlay(:)),max(underlay(:)),[],gray(256));
overlay_im = cmaplookup(overlay,-0.3,0.3,[],cmapsign2);

mask = ROI==1;

subplot(1,2,2);
hold on
imagesc(underlay_im);
imagesc(overlay_im, 'AlphaData', mask);
hold off
axis image
colormap(cmapsign2)
c = colorbar;
c.Ticks = [0 0.5 1];
c.TickLabels = {'-0.3';'0';'0.3'};
title('change in nsdgeneral voxel reliability** due to GLMsingle (\Delta{\itr})')
set(gca,'Fontsize',16)
xlabel('**plotting (FITHRF_GLMDENOISE_RR - ASSUMEHRF) reliabilities','Interpreter','none','FontSize',12);
xticks([])
yticks([])

set(gcf,'Position',[36 343 1116 674])

% Notice that there is systematic increase in reliability moving from the
% first to the second to the third to the final fourth version of the GLM
% results. These increases reflect, respectively, the addition of HRF
% fitting, the derivation and use of data-driven nuisance regressors, and
% the use of ridge regression as a way to regularize the instability of
% closely spaced experimental trials. Depending on one's experimental
% goals, it is possible with setting of option flags to activate a subset
% of these analysis features.
%
% Also, keep in mind that in the above figure, we are simply showing the
% median as a metric of the central tendency (you may want to peruse
% individual voxels in scatter plots, for example).

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function view = blurDialog(view);

if isempty(view)
   myErrorDlg('Volume window must be open to blur phase.');
   return;
end
if ~strcmp('Gray',view.viewType)
   myErrorDlg('Blur phase only in gray view.');
   return;
end
if (isempty(view.ph) | isempty(view.co))
   myErrorDlg('Correlation matrices not loaded.  Load them from the File menu');
   return;
end

% Make ui.blur structure
if ~isfield(view.ui,'blur');
    view.ui.blur = [];
end
if ~isfield(view.ui.blur,'Lambda');
   view.ui.blur.Lambda = 0;
end
if ~isfield(view.ui.blur,'weightFlag');
   view.ui.blur.weightFlag = 0;
end

windowPos = [150 200 300 50];

blurwin = figure('MenuBar','none');

set(blurwin,'Position', windowPos);
set(blurwin, 'Name', 'Blurring');

view.ui.blur.WindowHandle = gcf;

radioPos = [0.05 0.5 0.3 0.4];
radioOff = [0 0.4 0 0];
view.ui.blur.Radio(1) = uicontrol('Style','radiobutton',...
    'String','Blur',...
    'Units','normalized',...
    'Position',radioPos,...
    'Callback', [view.name,'.ui.blur.weightFlag=0; ',...
    'selectButton(',view.name,'.ui.blur.Radio,1)']);
view.ui.blur.Radio(2) = uicontrol('Style','radiobutton',...
    'String','Weighted blur',...
    'Units','normalized',...
    'Position',radioPos-radioOff,...
    'Callback', [view.name,'.ui.blur.weightFlag=1; ',...
    'selectButton(',view.name,'.ui.blur.Radio,2)']);

 % Select the appropriate radio button
 if view.ui.blur.weightFlag
    selectButton(view.ui.blur.Radio,2);
 else
    selectButton(view.ui.blur.Radio,1);
 end

view.ui.blur.SliderPos = [0.4 0.5 0.4 0.4];
view.ui.blur.LabelOffset = [0 0.4 0 0];

range = [-1 1];
sliderCallback = [view.name '.ui.blur.Lambda = get(' view.name '.ui.blur.Slider,''Value''); ',...
      'uicontrol(''Style'',''text'',''Units'',''normalized'',''Position'',',...
       view.name '.ui.blur.SliderPos - ',...
      view.name '.ui.blur.LabelOffset, ''String'',[''Blurring: '' num2str(10^(',...
      view.name '.ui.blur.Lambda))])'];
view.ui.blur.Slider = uicontrol('Style','slider',...
   'Units','normalized',...
   'Position',view.ui.blur.SliderPos,...
   'Min',range(1), 'Max',range(2),...
   'Value',view.ui.blur.Lambda,...
	'Callback',sliderCallback);

uicontrol('Style','text',...
   'Units','normalized',...
   'Position',view.ui.blur.SliderPos - view.ui.blur.LabelOffset,...
   'String',['Blurring: ' num2str(10^(view.ui.blur.Lambda))]);
OKposition = [0.85 0.1 0.1 0.8];
uicontrol('String','Blur',...
   'Units','normalized',...
   'Position', OKposition,...
   'Callback', ['close(',num2str(view.ui.blur.WindowHandle),'); ',...
      view.name,' = blurPh(',view.name,...
      ',',view.name,'.ui.blur.weightFlag',...
      ',10^(',view.name,'.ui.blur.Lambda)); ',...
      view.name,'=refreshScreen(',view.name,',1);']);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef MaOEAIT < ALGORITHM
% <multi/many> <real/integer> <constrained/none>
% Many-objective evolutionary algorithms based on an independent two-stage
% approach
% Evaluation1 --- 20000 --- Number of evaluations for NDWA
% Evaluation1 ---  6000 --- Number of evaluations for reference lines mapping
% epsilon     --- 0.999 --- Threshold in PCA

%------------------------------- Reference --------------------------------
% Y. Sun, B. Xue, M. Zhang, G. G. Yen, A new two-stage evolutionary
% algorithm for many-objective optimization, IEEE Transactions on
% Evolutionary Computation, 2019, 23(5): 748-761.
%------------------------------- 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
        function main(Algorithm,Problem)
            %% Parameter setting
            [Evaluation1,Evaluation2,epsilon] = Algorithm.ParameterSet(20000,6000,0.999);

            %% Generate random population
            Population = Problem.Initialization();
            Archive    = [];

            %% Generate the weight vectors
            [W,W_num] = UniformPoint(Problem.N,Problem.M);
            UW        = [W;W(end:-1:1,:)];
            num_W     = size(UW,1);

            %% The objective value of each solution on single-objective in different stage
            SF = zeros(Problem.N+floor(Problem.N/2)*2,1);

            %% The objective value of each solution on single-objective in phase 1
            SF(1:Problem.N) = sum(Population.objs.*UW(1,:),2);

            %% Non-dominated dynamic weight aggregation (NDWA)
            while Algorithm.NotTerminated(Population) && Problem.FE < Evaluation1
                W_index = mod(ceil(Problem.FE/Problem.N)+1,num_W);
                if W_index == 0
                    W_index = num_W;
                end
                MatingPool          = TournamentSelection(2,Problem.N,SF(1:Problem.N));
                Offspring           = OperatorGA(Problem,Population(MatingPool),{0.9,20,1,20});
                SF(Problem.N+1:end) = sum(Offspring.objs.*UW(W_index,:),2);
                [Population,SF]     = NDWA_EnvironmentalSelection([Population,Offspring],SF,Problem.N);
                Archive             = [Archive,Population];
                Archive             = Archive(NDSort(Archive.objs,1)==1);
            end

            %% Pareto-optimal subspace learning
            learned_data = FindSubspace(Archive.decs,epsilon);
            mean_val     = mean(learned_data,1);
            u_limit      = ones(size(Problem.upper));
            l_limit      = zeros(size(Problem.lower));
            for i = 1:numel(mean_val)
                if abs(learned_data(1,i)-mean_val(i)) < 0.1
                    l_limit(i) = round(mean_val(i),1);
                    u_limit(i) = round(mean_val(i),1);
                else
                    l_limit(i) = Problem.lower(i);
                    u_limit(i) = Problem.upper(i);
                end
            end

            %% Reference lines mapping
            % Optimize Problem.M single-objective optimization problems
            rf_list       = diag(ones(1,Problem.M));
            SO_maxeval    = floor(Evaluation2/Problem.M);
            Extreme_point = repmat(Population(1),1,Problem.M);
            for i = 1 : Problem.M
                Curr_evl        = Problem.FE;
                PopDec          = rand(Problem.N, Problem.D).*repmat(u_limit-l_limit,Problem.N,1) + repmat(l_limit,Problem.N,1);
                Population      = Problem.Evaluation(PopDec);
                SF(1:Problem.N) = cos_v_func(Population.objs,rf_list(i,:));
                while Algorithm.NotTerminated(Population) && Problem.FE < (Curr_evl+SO_maxeval)
                    MatingPool          = TournamentSelection(2,Problem.N,SF(1:Problem.N));
                    Offspring           = Operator(Problem,Population(MatingPool).decs,l_limit,u_limit);
                    SF(Problem.N+1:end) = cos_v_func(Offspring.objs,rf_list(i,:));
                    [~,Rank]            = sort(SF,1);
                    Population          = [Population,Offspring];
                    Population          = Population(Rank(1:Problem.N));
                    SF(1:Problem.N)     = SF(Rank(1:Problem.N));
                end
                Extreme_point(i) = Population(1);
            end
            Ideal_point = min(Extreme_point.objs,[],1);
            Nadir_point = max(Extreme_point.objs,[],1);
            RefPoint    = W.*repmat(Nadir_point-Ideal_point,size(W,1),1) + repmat(Ideal_point,size(W,1),1);

            %% Diversity maintaining
            % Optimize W_num single-objective optimization problems
            SO_maxeval = floor((Problem.maxFE-Evaluation1-Evaluation2)/W_num);
            Result     = repmat(Population(1),1,size(RefPoint,1));
            for i = 1 : W_num
                Curr_evl        = Problem.FE;
                PopDec          = rand(Problem.N,Problem.D).*repmat(u_limit-l_limit,Problem.N, 1) + repmat(l_limit,Problem.N,1);
                Population      = Problem.Evaluation(PopDec);
                SF(1:Problem.N) = cos_v_func(Population.objs,RefPoint(i,:));
                while Algorithm.NotTerminated(Population) && Problem.FE < (Curr_evl+SO_maxeval)
                    MatingPool          = TournamentSelection(2,Problem.N,SF(1:Problem.N));
                    Offspring           = Operator(Problem,Population(MatingPool).decs,l_limit,u_limit);
                    SF(Problem.N+1:end) = cos_v_func(Offspring.objs,RefPoint(i,:));
                    [~,Rank]            = sort(SF,1);
                    Population          = [Population,Offspring];
                    Population          = Population(Rank(1:Problem.N));
                    SF(1:Problem.N)     = SF(Rank(1:Problem.N));
                    if i == W_num && Problem.FE == (Curr_evl+SO_maxeval)
                        Result(W_num) = Population(1);
                        Population    = Result;
                    end
                end
                Result(i) = Population(1);
            end
        end
    end
end

function [Population,SF] = NDWA_EnvironmentalSelection(Population,SF,N)
% Environmental selection

    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;
    Last = find(FrontNo==MaxFNo);
    Next(Last(1:N-sum(Next))) = true;
    Population = Population(Next);
    SF(1:N)    = SF(Next);
end

function SF = cos_v_func(PopObj,v)
% Calculate the objective value of each solution on each single-objective
% optimization problem

    r1 = sum(PopObj.*v,2);
    r2 = sqrt(sum(PopObj.^2,2));
    SF = -r1./r2;
end

function Offspring = Operator(Problem,ParentDec,lower,upper)
% Simulated binary crossover and polynomial mutation

    [proC,disC,proM,disM] = deal(1,20,1,20);
    Parent1 = ParentDec(1:floor(end/2),:);
    Parent2 = ParentDec(floor(end/2)+1:floor(end/2)*2,:);
    [N,D]   = size(Parent1);

    beta = zeros(N,D);
    mu   = rand(N,D);
    beta(mu<=0.5) = (2*mu(mu<=0.5)).^(1/(disC+1));
    beta(mu>0.5)  = (2-2*mu(mu>0.5)).^(-1/(disC+1));
    beta = beta.*(-1).^randi([0,1],N,D);
    beta(rand(N,D)<0.5) = 1;
    beta(repmat(rand(N,1)>proC,1,D)) = 1;
    Offspring = [(Parent1+Parent2)/2+beta.*(Parent1-Parent2)/2
                 (Parent1+Parent2)/2-beta.*(Parent1-Parent2)/2];

    Lower = repmat(lower,2*N,1);
    Upper = repmat(upper,2*N,1);
    Site  = rand(2*N,D) < proM/D;
    mu    = rand(2*N,D);
    temp  = Site & mu<=0.5;
    Offspring       = min(max(Offspring,Lower),Upper);
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                      (1-(Offspring(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp = Site & mu>0.5;
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                      (1-(Upper(temp)-Offspring(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));

    Offspring = Problem.Evaluation(Offspring);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Copyright (C) 2007 Peter Carbonetto. All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% Author: Peter Carbonetto
%         Dept. of Computer Science
%         University of British Columbia
%         May 19, 2007

function [dR, dS] = computeJGGradient (qR, qS, auxdata)
  [K C f Rv Rf Sv Sf NS d] = deal(auxdata{:});

  nr = length(Rv);   % The number of large regions.
  ns = length(Sv);   % The number of small regions.

  % Reshape the input vectors.
  qR = reshapemarginals(qR,Rv,K);
  qS = reshapemarginals(qS,Sv,K);

  % Compute the gradient terms for the large regions. Repeat for each
  % large region.
  dR = cell(1,nr);
  for r = 1:nr
    is    = Rv{r};  % The variable nodes in the large region.
    table = 1 ./ qR{r};

    % Multiply the table by all the factors in the large region.
    for j = Rf{r}
      table = multiplyfactors(table,is,f{j},C{j});
    end

    dR{r} = -log(table) + 1;
  end

  % Compute the gradient terms for the small regions. Repeat for each
  % small region.
  dS = cell(1,ns);
  for s = 1:ns
    is    = Sv{s};  % The variable nodes in the small region.
    table = 1 ./ qS{s};

    % Multiply the table by all the factors in the small region.
    for j = Sf{s}
      table = multiplyfactors(table,is,f{j},C{j});
    end

    dS{s} = (1 - d(s)) * (-log(table) + 1);
  end

  % Convert the computed gradients into vectors.
  dR = vectorize(dR);
  dS = vectorize(dS);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function result = polygon_x_2d ( n, v )

%*****************************************************************************80
%
%% POLYGON_X_2D integrates the function X over a polygon in 2D.
%
%  Discussion:
%
%    The polygon is bounded by the points (X(1:N), Y(1:N)).
%
%    INTEGRAL = (1/6) * sum ( 1 <= I <= N )
%      ( X(I)^2 + X(I) * X(I-1) + X(I-1)^2 ) * ( Y(I) - Y(I-1) )
%
%    where X(0) and Y(0) should be replaced by X(N) and Y(N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    S F Bockman,
%    Generalizing the Formula for Areas of Polygons to Moments,
%    American Mathematical Society Monthly,
%    1989, pages 131-132.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%    N should be at least 3 for a nonzero result.
%
%    Input, real V(2,N), the coordinates of the vertices
%    of the polygon.  These vertices should be given in counter-clockwise order.
%
%    Output, real RESULT, the value of the integral.
%
  result = 0.0;

  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_X_2D - Warning!\n' );
    fprintf ( 1, '  The number of vertices must be at least 3.\n' );
    fprintf ( 1, '  The input value of N = %d\n', n );
    return
  end

  for i = 1 : n

    if ( i == 1 )
      im1 = n;
    else
      im1 = i - 1;
    end

    result = result + ( v(1,i).^2 + v(1,i) * v(1,im1) + v(1,im1).^2 ) ...
      * ( v(2,i) - v(2,im1) );

  end

  result = result / 6.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [pvec, pstruct] = tapas_bayes_optimal_categorical_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% This is a dummy function

pvec    = [];
pstruct = struct;

return;

author: translationalneuromodeling

repo: tapas


text:
function value = ball_volume_nd ( n, r )

%*****************************************************************************80
%
%% BALL_VOLUME_ND computes the volume of a ball in ND.
%
%  Integration region:
%
%    Points X(1:N) such that
%
%      Sum ( X(1:N)**2 ) <= R**2
%
%  Discussion:
%
%    N  Volume
%
%    2             PI    * R**2
%    3  (4/3)    * PI    * R**3
%    4  (1/2)    * PI**2 * R**4
%    5  (8/15)   * PI**2 * R**5
%    6  (1/6)    * PI**3 * R**6
%    7  (16/105) * PI**3 * R**7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the space.
%
%    Input, real R, the radius of the ball.
%
%    Output, real BALL_VOLUME_ND, the volume of the ball.
%
  value = ball_unit_volume_nd ( n ) * r^n;

  return
end


text:
function [g] = spm_gx_SHC(x,v,P)
% maps to state of a SCH to a 2-D position in the world
% FORMAT [g] = spm_gx_SHC(x,v,P)
%
% x    - vector of hidden sates
% P.g  - state-space location associated with each hidden states
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_gx_SHC.m 3113 2009-05-11 15:25:13Z karl $

% expected position (coordinates in P.g)
%--------------------------------------------------------------------------
px     = exp(x(:));
px     = px/sum(px);
g      = P.g*px;

author: spm


text:
 clc;
 clear all;
I1=imread('E:\sun\ͼƬ\clockA.bmp');
I2=imread('E:\sun\ͼƬ\clockB.bmp');
%  I1=colorspace('yuv<-rgb',I1);
% % %I2=colorspace('yuv<-rgb',I2);
%  I1=I1(:,:,1);
% %Y2=I1(:,:,2);
%
% % %I1=imread('D:\MATLAB701\work\Lena.jpg');
% % I1=imread('E:\sun\ͼƬ\clockA.bmp');
% % I2=imread('E:\sun\ͼƬ\clockB.bmp');
I1=double(I1);
I2=double(I2);
I=(I1+I2)/2;        %??
n=1;
[Y1,h1] = dtwavexfm2(I1,n,'near_sym_b','qshift_b');
[Y2,h2] = dtwavexfm2(I2,n,'near_sym_b','qshift_b');
[Y,h] = dtwavexfm2(I,n,'near_sym_b','qshift_b');

[r,c]=size(h1{:,:,1});
for i=1:6
   h_seg(:,:,i)=seg(h{1,1}(:,:,i));%??
   F=zeros(r,c);
  L2=bwlabel(h_seg(:,:,i));
  ma=max(max(L2));
      for k=1:ma
    I=L2==k;
     index1=find(h_seg(:,:,i)==1);
     F(index1)=W(i,k)*I1(index)+(1-W(i,k))*I2(index1);
      end
  Fu(:,:,i)=F;
end
Fuse= dtwaveifm2(Yl,Fu,'near_sym_b','qshift_b');
h=Entrophy(Fuse);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function piece_num = slice ( dim_num, slice_num )

%*****************************************************************************80
%
%% SLICE: maximum number of pieces created by a given number of slices.
%
%  Discussion:
%
%    If we imagine slicing a pizza, each slice produce more pieces.
%    The position of the slice affects the number of pieces created, but there
%    is a maximum.
%
%    This function determines the maximum number of pieces created by a given
%    number of slices, applied to a space of a given dimension.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Banks,
%    Slicing Pizzas, Racing Turtles, and Further Adventures in
%    Applied Mathematics,
%    Princeton, 1999,
%    ISBN13: 9780691059471,
%    LC: QA93.B358.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer SLICE_NUM, the number of slices.
%
%    Input, integer PIECE_NUM, the maximum number of pieces that can
%    be created by the given number of slices applied in the given dimension.
%
  piece_num = 0;
  for j = 0 : min ( dim_num, slice_num )
    piece_num = piece_num + i4_choose ( slice_num, j );
  end

  return
end


text:
% VL_VERSION  Obtain VLFeat version information
%   VL_VERSION() prints VLFeat version number. VL_VERSION('VERBOSE')
%   prints additional information.
%
%   STR = VL_VERSION(...) returns the ouptut in a string rather than
%   printing it.
%
%   See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [invA]=invltod(A,n)

% computes the inverse of a n*n lower triangular matrix with ones on the main diagonal
% uses result (XXX) from Appendix (XXX)

% first obtain the B matrix from A
B=sparse(tril(A,-1));

% compute the summation term
summ=sparse(n,n);
prodt=speye(n);
for ii=1:n-1
prodt=prodt*B;
summ=summ+(-1)^ii*prodt;
end

% compute the inverse of A
invA=full(speye(n)+summ);

author: european-central-bank

repo: BEAR-toolbox


text:
function data = ReadChoreSpineFile(filename)

nspinepts = 11;
nparams = 2*nspinepts + 2;

data = [];

fid = fopen(filename,'r');
if fid < 1,
  error('Could not open file %s for reading',filename);
end

ids = [];
expname = '';
while true,
  l = fgetl(fid);
  if ~ischar(l),
    break;
  end
  if isempty(l),
    continue;
  end

  % first word is the experiment name
  [expnamecurr,count,~,nextindex] = sscanf(l,'%s',1);
  if count < 1,
    warning('Could not read experiment name from line >%s<',l);
    continue;
  end
  if ~isempty(expname) && ~strcmp(expname,expnamecurr),
    warning('Experiment name %s does not match previous experiment name %s',expnamecurr,expname);
  end

  % read in the rest of the data
  l = l(nextindex:end);
  [datacurr,count] = sscanf(l,'%f',nparams);
  if count < nparams,
    warning('Could not parse data from line >%s<',l);
    continue;
  end

  id = datacurr(1);
  timestamps = datacurr(2);
  value = datacurr(3:2+nspinepts*2);
  x = value(1:2:end-1);
  y = value(2:2:end);

  idi = find(id == ids,1);
  if isempty(idi),
    idi = numel(ids)+1;
    newdata = struct('id',id,'timestamps',timestamps,'xspine',x(:),'yspine',y(:));
    if isempty(ids),
      data = newdata;
    else
      data(idi) = newdata;
    end
    ids(idi) = id;
  else
    data(idi).timestamps(end+1) = timestamps;
    data(idi).xspine(:,end+1) = x;
    data(idi).yspine(:,end+1) = y;
  end

end

% sort by ids
[ids,order] = sort(ids); %#ok<ASGLU>
data = data(order);

% sort by timestamps
for i = 1:numel(data),
  [data(i).timestamps,order] = sort(data(i).timestamps);
  data(i).xspine = data(i).xspine(:,order);
  data(i).yspine = data(i).yspine(:,order);
end

fclose(fid);

author: kristinbranson

repo: JAABA


text:
function r = atan(p)
%atan: this function calculates the atan of a x, y, and z values of
% the structures that I use in the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

r = p;
r.value = atan(p.value);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
% Local Regression and Likelihood, Figure 11.4.
% Author: Catherine Loader
%
% Local Adaptive Smooth - Dopler example.
%
% NEED: Improve plot - curve is hard to see.

x = (0:2047)'/2047;
m = 20*sqrt(x.*(1-x)) .* sin(2*pi*1.05 ./ (x+0.05));
y = m + normrnd(0,1,2048,1);
fit = locfit(x,y,'pen',4,'acri','cp','maxk',500);
figure('Name','fig11_4: local adaptive smooth: dopler example');
lfplot(fit,'red');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function jed = ymdf_to_jed_zoroastrian ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_ZOROASTRIAN converts a Zoroastrian YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%
  jed_epoch = epoch_to_jed_zoroastrian ( );

  jed = jed_epoch + ( d - 1 ) + 30 * ( m - 1 ) + 365 * ( y - 1 ) + f;

  return
end


text:
% qpwls_std_test
% examine FFT-based prediction of QPWLS estimate pixel variance
% particularly w.r.t. effect of "mask"

if ~isvar('G2')
	down = 1;
	ig1 = image_geom('nx', 512, 'ny', 504, 'fov', 500, 'down', down);
	ig2 = ig1;
	ig2.mask = ig2.circ(220) > 0;
	sg = sino_geom('ge1', 'down', down);

	% now a double sized version
	ig3 = image_geom('nx', ig1.nx*2, 'ny', ig1.nx*2, 'fov', ig1.fov*2);

	G1 = Gtomo_dd(sg, ig1);
	G2 = Gtomo_dd(sg, ig2);
	G3 = Gtomo_dd(sg, ig3);
	im pl 2 1
	im(1, ig1.mask, 'mask1') % square
	im(2, ig2.mask, 'mask1') % circle
prompt
end

if ~isvar('std2')
	l2bs = [5:14];

	for ii=1:length(l2bs)
		l2b = l2bs(ii);
		R1 = Robject(ig1.mask);
		R2 = Robject(ig2.mask);
		R3 = Robject(ig3.mask);

		[psf var fw1(ii)] = qpwls_psf(G1, R1, 2^l2b, ig1.mask);
		std1(ii) = sqrt(var);
		[psf var fw2(ii)] = qpwls_psf(G2, R2, 2^l2b, ig2.mask);
		std2(ii) = sqrt(var);
		[psf var fw3(ii)] = qpwls_psf(G3, R3, 2^l2b, ig3.mask);
		std3(ii) = sqrt(var);
	end
end

if im
	clf
	plot(l2bs, std2, 'o', l2bs, std1, 's', l2bs, std3, '+')
	xlabel 'log_2(\beta)'
	ylabel 'relative standard deviation'
	legend('circle mask', 'full mask', 'double size')
	ir_savefig fig_qpwls_std_fan
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef(Abstract) AbstractFirstOrderIntegrator < AbstractIntegrator
    %AbstractFirstOrderIntegrator Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        [t,y,te,ye,ie] = integrate(obj, odefun, tspan, y0, evtsFunc, odeOutputFun)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [M, moral_edges] = moralize(G)
% MORALIZE Ensure that for every child, all its parents are married, and drop directionality of edges.
% [M, moral_edges] = moralize(G)

M = G;
n = length(M);
for i=1:n
  fam = family(G,i);
  M(fam,fam)=1;
end
M = setdiag(M,0);
moral_edges = sparse(triu(max(0,M-G),1));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef ImageProcessingTestsSuite < handle

    methods

        function obj = ImageProcessingTestsSuite()
            path = './tests/Source/ImageProcessingTests/ImageProcessingTests.m';
            suite = matlab.unittest.TestSuite.fromFile(path, 'Tag','ImageProcessing');
            results = suite.run;
            table(results)
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function vw = contrastGUI(vw, flag)
%
% vw = contrastGUI([vw=cur view])
%
% Dialog to set parameters for a contrast map.
%
% Inputs:
%  stim: stim ('stim') structure from er_concatParfiles.
%  initParams: initial params structure specifying initial values for the
%       relevant fields (same format as output p struct below).
%
% Outputs: parameters structure p with fields:
%   active: initial list of active conditions, or weights of active conditions.
%   control:  initial list of control conditions, or weights.
%   name: contrast name.
%   options: [to be added]
%
%
% ras, 01/2007.

% first, check if this is a callback from one of the uicontrols --
% if it is, it will have a flag as the third argument, describing what to
% do next:
if exist('flag','var')
    s = get(gcf,'UserData'); % this struct will be useful for the callbacks

    switch flag
        case 1,     num = get(gcbo,'UserData');  setControlCond(s, num);
        case 2,     num = get(gcbo,'UserData');  setActiveCond(s, num);
        case 3,     setSaveNameUI(s);
        case 4,     vw = callComputeContrastMap(s, vw);
        otherwise,  error('Illegal flag entered as third param.')
    end

    return
end

if notDefined('vw'),  vw = getCurView;  end

%% open the figure
s = openFig(vw);

%% append an 'Advanced Options' params
advancedOptionsPanel(s);

return
% /---------------------------------------------------------------------/ %

%% /---------------------------------------------------------------------/ %
function s = openFig(vw)
% s = openFig(stim);
% creates the interface window and returns the GUI structure s.

stim = er_concatParfiles(vw);
nConds = length(stim.condNums);
winHeight = (30 * nConds) + 180; % do this in pixels, 30 pix per condition row
winWidth = 450;
bgColor = [.8 .8 .8];
fgColor = [0 0 0];
s.h(1) = figure('Name', 'Statistical Contrast Map', 'Color', bgColor, ...
    'Position', [300 300 winWidth winHeight],...
    'NumberTitle', 'off', 'MenuBar','none');

% add title text fields
tittxt = sprintf('Compute a Contrast Map');
uicontrol('Style','text','String',tittxt,...
    'FontName','Helvetica','FontSize',18,'FontWeight','bold',...
    'HorizontalAlignment','center',...
    'ForegroundColor', fgColor, 'BackgroundColor', bgColor, ...
    'Units','Normalized', ...
    'Position',[10/winWidth (winHeight-30)/winHeight 430/winWidth 30/winHeight]);

uicontrol('Style', 'text', 'String', viewGet(vw, 'annotation'),...
    'FontName','Helvetica','FontSize',16,'FontWeight','bold',...
    'HorizontalAlignment','center',...
    'ForegroundColor', fgColor, 'BackgroundColor', bgColor,...
    'Units','Normalized','Position',[10/winWidth (winHeight-60)/winHeight 430/winWidth 30/winHeight]);

% make a top row of labels
uicontrol('Style', 'text', 'String', 'Cond #', ...
    'ForegroundColor', fgColor, 'BackgroundColor', bgColor, ...
    'FontName', 'Helvetica', 'FontSize', 12, 'FontWeight', 'bold', ...
    'Units', 'Normalized', ...
    'Position', [30/winWidth (winHeight-90)/winHeight 60/winWidth 30/winHeight]);

uicontrol('Style','text','String','Control?', ...
    'FontName', 'Helvetica', 'FontSize', 12, 'FontWeight', 'bold', ...
    'ForegroundColor', 'w', 'BackgroundColor', [1 0 0],...
    'Units','Normalized','Position',[360/winWidth (winHeight-90)/winHeight 60/winWidth 30/winHeight]);

uicontrol('Style','text','String','Condition Name',...
    'FontName', 'Helvetica', 'FontSize', 12, 'FontWeight', 'bold', ...
    'ForegroundColor', fgColor, 'BackgroundColor', bgColor, ...
    'Units','Normalized','Position',[210/winWidth (winHeight-90)/winHeight 120/winWidth 30/winHeight]);

uicontrol('Style','text','String','Active?', ...
    'FontName', 'Helvetica', 'FontSize', 12, 'FontWeight', 'bold', ...
    'ForegroundColor', 'w', 'BackgroundColor', [0 0 1],...
    'Units','Normalized','Position',[120/winWidth (winHeight-90)/winHeight 60/winWidth 30/winHeight]);

% make labels for each cond num
for i = 1:nConds
    ypos = (winHeight - 90 -  30*i)/winHeight;
    s.condNumHandles(i) = uicontrol('Style', 'text', ...
        'String', num2str(stim.condNums(i)), ...
        'FontName', 'Helvetica', 'FontSize', 12, ...
        'ForegroundColor', fgColor, 'BackgroundColor', bgColor,...
        'Units', 'Normalized', ...
        'Position', [30/winWidth ypos 60/winWidth 30/winHeight]);
end

% make control checkboxes for each cond
cbstr = sprintf('contrastGUI(%s, 1);', vw.name);
for i = 1:nConds
    ypos = (winHeight - 90 -  30*i)/winHeight;
    s.controlHandles(i) = uicontrol('Style','checkbox','Value',0,...
        'FontName','Helvetica','FontSize',12,...
        'ForegroundColor',fgColor,'BackgroundColor',bgColor,...
        'UserData',i,'Callback',cbstr,...
        'Units','Normalized','Position',[360/winWidth ypos 60/winWidth 30/winHeight]);

    % add a weights control
    s.controlWeightEdits(i) = uicontrol('Style', 'edit', 'String','',...
        'FontName', 'Helvetica', 'FontSize', 8,...
        'ForegroundColor', [0 0 0], 'BackgroundColor', [1 1 1],...
        'UserData', i, 'Callback', cbstr, 'Enable', 'off', ...
        'Visible', 'off', 'Units', 'Normalized', ...
        'Position', [(380/winWidth) (ypos + 8/winWidth) 40/winWidth 20/winHeight]);

end

% make text fields for each cond name
for i = 1:nConds
    ypos = (winHeight - 90 -  30*i)/winHeight;
    s.condNameHandles(i) = uicontrol('Style', 'text', ...
        'String', stim.condNames{i},...
        'FontName', 'Helvetica', 'FontSize', 12,...
        'ForegroundColor', fgColor, 'BackgroundColor', [1 1 1],...
        'Units', 'Normalized', ...
        'Position',[210/winWidth ypos 120/winWidth 30/winHeight]);
end

% make active checkboxes for each cond
cbstr = sprintf('contrastGUI(%s, 2);', vw.name);
for i = 1:nConds
    ypos = (winHeight - 90 -  30*i)/winHeight;
    s.activeHandles(i) = uicontrol('Style','check','Value',0,...
        'FontName','Helvetica','FontSize',12,'FontWeight','bold',...
        'ForegroundColor', fgColor, 'BackgroundColor',bgColor,...
        'UserData',i,'Callback',cbstr,...
        'Units', 'Normalized', ...
        'Position', [120/winWidth ypos 60/winWidth 30/winHeight]);

    s.activeWeightEdits(i) = uicontrol('Style', 'edit', 'String', '', ...
        'FontName', 'Helvetica', 'FontSize', 8, ...
        'ForegroundColor', fgColor, 'BackgroundColor', [1 1 1],...
        'UserData', i, 'Callback', cbstr, 'Enable', 'off', ...
        'Visible', 'off', 'Units', 'Normalized', ...
        'Position', [(140/winWidth) (ypos + 8/winWidth) 40/winWidth 20/winHeight]);

end

% make button to set save file path
cbstr = sprintf('contrastGUI(%s, 3);', vw.name);
s.setPathHandle = uicontrol('Style', 'pushbutton', 'String', 'Name:',...
    'FontName', 'Helvetica', 'FontSize',12,...
    'ForegroundColor', 'k', 'BackgroundColor', [.8 .8 .8],...
    'Callback', cbstr, 'Enable', 'off', ...
    'Units', 'Normalized', ...
    'Position', [120/winWidth 50/winHeight 60/winWidth 30/winHeight]);

% make field with save filename
s.saveHandle = uicontrol('Style', 'edit', 'String', 'V ',...
    'FontName', 'Helvetica', 'FontSize',12, 'FontWeight', 'bold', ...
    'ForegroundColor', [0 0 0], 'BackgroundColor', [1 1 1], ...
    'Enable','off',...
    'Units', 'Normalized', ...
    'Position', [120/winWidth 20/winHeight 210/winWidth 30/winHeight]);

% make 'GO' button
cbstr = sprintf('%s = contrastGUI(%s, 4);',vw.name, vw.name);
s.goHandle = uicontrol('Style','pushbutton', 'String', 'GO',...
    'FontName', 'Helvetica', 'FontSize', 24, 'FontWeight', 'bold',...
    'ForegroundColor',[1 1 1],'BackgroundColor',[.1 .7 .1],...
    'UserData', stim.condNums,...
    'Callback', cbstr, 'Enable', 'off', ...
    'Units', 'Normalized', ...
    'Position', [360/winWidth 20/winHeight 60/winWidth 60/winHeight]);

% put the s struct, containing important info, in the figure UserData
s.condNames = stim.condNames;
s.condNums = stim.condNums;
set(gcf, 'UserData', s);

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function s = advancedOptionsPanel(s)
% Add a uipanel containing advanced controls for a contrast map,
% such as setting the type of statistical test (F or T), the
% units of the resulting map (-log(p), T, F, p, ces, etc), and
% whether to activate the weights edit controls.

s.optsPanel = mrvPanel('below', .25);

% popup menu to set test type
s.testPopup = uicontrol('Parent', s.optsPanel, 'Style', 'popup', ...
    'Units', 'normalized', 'Position', [.1 .4 .2 .15], ...
    'FontName', 'Arial', 'FontSize', 12, ...
	'BackgroundColor', get(gcf, 'Color'), ...
    'String', {'T test' 'F test'}, 'Value', 1);

uicontrol('Parent', s.optsPanel, 'Style', 'text', ...
    'Units', 'normalized', 'Position', [.1 .6 .2 .1], ...
    'String', 'Test Type', 'FontName', 'Helvetica', ...
	'BackgroundColor', get(gcf, 'Color'), ...
    'FontSize', 10, 'FontWeight', 'bold', 'FontAngle', 'italic');

% popup menu to set units of map
vals = {'-log(p)' 'p' 'T' 'F' 'Contrast Effect Size'};
s.unitPopup = uicontrol('Parent', s.optsPanel, 'Style', 'popup', ...
    'Units', 'normalized', 'Position', [.4 .4 .2 .15], ...
    'FontName', 'Arial', 'FontSize', 12, ...
	'BackgroundColor', get(gcf, 'Color'), ...
    'String', vals, 'Value', 1);

uicontrol('Parent', s.optsPanel, 'Style', 'text', ...
    'Units', 'normalized', 'Position', [.4 .6 .2 .1], ...
    'String', 'Units For Map?', 'FontName', 'Helvetica', ...
	'BackgroundColor', get(gcf, 'Color'), ...
    'FontSize', 10, 'FontWeight', 'bold', 'FontAngle', 'italic');

% checkbox to toggle visibility of edit fields for condition weights
cb = ['TMP = get(gcf, ''UserData''); ' ...
      'vis = get(gcbo, ''Value'') + 1; ' ...
      'onoff = {''off'' ''on''}; ' ...
      'set(TMP.activeWeightEdits, ''Visible'', onoff{vis}); ' ...
      'set(TMP.controlWeightEdits, ''Visible'', onoff{vis}); ' ...
      'clear TMP vis onoff '];
uicontrol('Parent', s.optsPanel, 'Style', 'checkbox', ...
    'Units', 'normalized', 'Position', [.1 .1 .5 .15], ...
    'String', 'Set Condition Weights', 'FontName', 'Helvetica', ...
	'BackgroundColor', get(gcf, 'Color'), ...
    'FontSize', 10, 'Callback', cb);

% initialize panel to invisible
mrvPanelToggle(s.optsPanel, 'off');

% add a button to the main GUI to toggle this panel
cb = ['TMP = get(gcf, ''UserData''); ' ...
      'mrvPanelToggle(TMP.optsPanel); ' ...
      'clear TMP '];
s.panelToggle = uicontrol('Style', 'pushbutton', ...
    'String', 'More...', 'Units', 'norm', 'Position', [.1 .03 .15 .06], ...
    'BackgroundColor', [.9 .9 .9], 'ForegroundColor', 'k', ...
    'Callback', cb);

% update userdata of GUI to include this info
set(gcf, 'UserData', s);

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function setControlCond(s,num)
% set the specified num condition as a control condition,
% highlighting the selected condition number and checking if ready to
% proceed

[active control w] = loadActiveControl(s);

bgColor = get(gcf, 'Color');

if get(s.controlHandles(num),'Value')==1
	set(s.controlHandles(num),'BackgroundColor',[1 0 0]);
	set(s.activeHandles(num),'Value',0,'BackgroundColor',bgColor);
    set(s.condNameHandles(num),'ForegroundColor','w','BackgroundColor',[1 0 0]);
    set(s.condNumHandles(num),'ForegroundColor','w','BackgroundColor',[1 0 0]);
    set(s.controlWeightEdits(num),'Enable','On');
    set(s.activeWeightEdits(num),'Enable','Off','String','');

    if isempty(get(s.controlWeightEdits(num),'String'))
        set(s.controlWeightEdits(num),'String','1');
    end

else
	set(s.controlHandles(num),'BackgroundColor',bgColor);
    set(s.condNameHandles(num),'ForegroundColor','k','BackgroundColor',[1 1 1]);
    set(s.condNumHandles(num),'ForegroundColor','k','BackgroundColor',bgColor);
    set(s.controlWeightEdits(num),'Enable','Off','String','');

end

setEvenWeights(s);

checkIfReady(s);

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function setActiveCond(s,num)
% set the specified num condition as an active condition,
% highlighting the selected condition number and checking if ready to
% proceed

[active control w] = loadActiveControl(s);

bgColor = get(gcf, 'Color');

if get(s.activeHandles(num),'Value')==1
	set(s.activeHandles(num),'BackgroundColor',[0 0 1]);
    set(s.controlHandles(num),'Value',0,'BackgroundColor',bgColor);
    set(s.condNameHandles(num),'ForegroundColor','w','BackgroundColor',[0 0 1]);
    set(s.condNumHandles(num),'ForegroundColor','w','BackgroundColor',[0 0 1]);
    set(s.activeWeightEdits(num),'Enable','On');
    set(s.controlWeightEdits(num),'Enable','Off','String','');

    if isempty(get(s.activeWeightEdits(num), 'String'))
        set(s.activeWeightEdits(num), 'String', '1');
    end

else
	set(s.activeHandles(num),'BackgroundColor',bgColor);
    set(s.condNameHandles(num),'ForegroundColor','k','BackgroundColor',[1 1 1]);
    set(s.condNumHandles(num),'ForegroundColor','k','BackgroundColor',bgColor);
    set(s.activeWeightEdits(num),'Enable','Off','String','');

end

setEvenWeights(s);

checkIfReady(s);

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function checkIfReady(s)
% update the name of the save file based on the highlighted conditions, and
% check if at least one control and one active condition are selected -- if
% so, enable the GO button.
atLeast1C = 0; atLeast1A = 0;

for i = 1:length(s.condNameHandles)
    s.condNames{i} = get(s.condNameHandles(i),'String');
end

saveName = '';
for i = 1:length(s.activeHandles)
    if get(s.activeHandles(i),'Value')==1
        saveName = [saveName rmBlanks(s.condNames{i})];
        atLeast1C = 1;
    end
end

saveName = [saveName 'V'];

for i = 1:length(s.controlHandles)
    if get(s.controlHandles(i),'Value')==1
        saveName = [saveName rmBlanks(s.condNames{i})];
        atLeast1A = 1;
    end
end

set(s.saveHandle,'String',saveName);

if atLeast1C && atLeast1A
    set(s.setPathHandle,'Enable','on');
    set(s.saveHandle,'Enable','on');
    set(s.goHandle,'Enable','on');
else
    set(s.setPathHandle,'Enable','off');
    set(s.saveHandle,'Enable','off');
    set(s.goHandle,'Enable','off');
end

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function str = rmBlanks(str)
% removes all the blanks from str.
ok = ones(size(str));
loc = findstr(' ',str);
ok(loc) = 0;
str = str(find(ok));
return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function setSaveNameUI(s,num)
[fname,pth] = myUIPutFile('Inplane/Original/','*.mat','Select a save file name / path: ' );
saveName = fullfile(pth,fname);
set(s.saveHandle,'String',saveName);
return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function setEvenWeights(s)
% make the weights for all selected active and control conditions
% equal. (sum of active weights = 1; sum of control weights = -1).

[active control w] = loadActiveControl(s);

if sum(w) ~= 0
    if ~isempty(active)
        w(active) = 1 / length(active);
    end

    if ~isempty(control)
        w(control) = -1 / length(control);
    end
end

for i = 1:length(s.activeHandles)
    if get(s.activeHandles(i),'Value')==1
        set(s.activeWeightEdits(i), 'String', num2str(w(i)));

    elseif get(s.controlHandles(i),'Value')==1
        set(s.controlWeightEdits(i), 'String', num2str(w(i)));

    end
end

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function [active, control, w] = loadActiveControl(s)
% [active, control, w] = loadActiveControl(s)
% get the active and control conditions based on the GUI settings.
% active: list of active conditions.
% control: list of control conditions.
% w: weights vector, where graded weights can be used for
% each condition. (+ weight = active, -weight = control).

active = zeros(1,length(s.activeHandles));
control = zeros(1,length(s.activeHandles));
w = zeros(1, length(s.activeHandles));

for i = 1:length(s.activeHandles)
    if get(s.activeHandles(i),'Value')==1
        active(i) = 1;

        if ~isempty( get(s.activeWeightEdits(i), 'String') )
            w(i) = str2num(get(s.activeWeightEdits(i),'String'));
        end
    elseif get(s.controlHandles(i),'Value')==1
        control(i) = 1;

        if ~isempty( get(s.controlWeightEdits(i), 'String') )
            w(i) = -str2double(get(s.controlWeightEdits(i),'String'));
        end
    end
end

active = find(active);
control = find(control);

return
% /---------------------------------------------------------------------/ %

% /---------------------------------------------------------------------/ %
function vw = callComputeContrastMap(s, vw)
% callback for the big 'GO' button; this takes the UI info
% contained in the s struct, and calls the computeContrastMap
% function:
hgui = gcf;

[active control w] = loadActiveControl(s);

% account for null condition: -1 offset to get to condNums
active = active - 1;
control = control - 1;

% get advanced options (test type, units)
testList = {'T' 'F'};
test = testList{ get(s.testPopup, 'Value') };
unitList = {'log10p' 'p' 't' 'f' 'ces'};
units = unitList{ get(s.unitPopup, 'Value') };

saveName = get(s.saveHandle,'String');
if isempty(findstr(filesep,saveName)) % check that absolute path not specified
    %     savePath = fullfile(dataDir(vw),['contrastMap_' saveName]);
    savePath = fullfile(dataDir(vw),saveName);
else
    savePath = saveName;
end

condNums = get(gcbo,'UserData');

% to do: use weights vector
vw = computeContrastMap2(vw, active, control, saveName, ...
                           'test', test, 'mapUnits', units);

close(hgui);

return
% /---------------------------------------------------------------------/ %


text:
%DEMO_AUDIOSHRINK  Decomposition into tonal and transient parts
%
%   This demos shows how to do audio coding and "tonal + transient"
%   decomposition using group lasso shrinkage of two |wmdct| transforms
%   with different time-frequency resolutions.
%
%   The signal is transformed using two orthonormal |wmdct| bases.
%   Then group lasso shrinkage is applied to the two transforms
%   in order to:
%
%     * select fixed frequency lines of large |wmdct| coefficients on the
%       wide window |wmdct| transform
%
%     * select fixed time lines of large |wmdct| coefficients on the
%       narrow window |wmdct| transform
%
%   The corresponding approximated signals are computed with the
%   corresponding inverse, |iwmdct|.
%
%   .. figure::
%
%      Plots and time-frequency images
%
%      The upper plots in the figure show the tonal parts of the signal, the
%      lower plots show the transients. The TF-plots on the left are the
%      corresponding wmdct coefficients found by appropriate group lasso
%      shrinkage
%
%   Corresponding reconstructed tonal and transient sounds may be
%   listened from arrays rec1 and rec2 (sampling rate: 44.1 kHz)
%

% Load audio signal and add noise
% -------------------------------
% Use the 'glockenspiel' signal.
sig=gspi;
fs=44100;

% Shorten signal
siglen = 2^16;
sig = sig(1:siglen);

% Add Gaussian white noise
nsig = sig + 0.01*randn(size(sig));

% Tonal layer
% -----------

% Create a WMDCT basis with 256 channels
F1=frametight(frame('wmdct','gauss',256));

% Group lasso and invert
c1 = franagrouplasso(F1,nsig,0.8,'soft','freq');
rec1 = frsyn(F1,c1);

% Transient layer
% ---------------

% Create a WMDCT basis with 32 channels
F2=frametight(frame('wmdct','gauss',32));

c2 = franagrouplasso(F2,nsig,0.5,'soft','time');
rec2 = frsyn(F2,c2);

% Plots
% -----

% Dynamic range for plotting
dr=50;
xplot=(0:siglen-1)/fs;

figure(1);
subplot(2,2,1);
plot(xplot,rec1);
xlabel('Time (s)');
axis tight;

subplot(2,2,2);
plotframe(F1,c1,fs,dr);

subplot(2,2,3);
plot(xplot,rec2);
xlabel('Time (s)');
axis tight;

subplot(2,2,4);
plotframe(F2,c2,fs,dr);

% Count the number of non-zero coefficients
N1=sum(abs(c1)>0);
N2=sum(abs(c2)>0);

p1 = 100*N1/siglen;
p2 = 100*N2/siglen;
p=p1+p2;

fprintf('Percentage of retained coefficients: %f + %f = %f\n',p1,p2,p);

disp('To play the original, type "soundsc(sig,fs)"');
disp('To play the tonal part, type "soundsc(rec1,fs)"');
disp('To play the transient part, type "soundsc(rec2,fs)"');

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function asSetAllTitlesToImageString
global asObjs
for i = 1 : length(asObjs)
%     asObjs(i).roi.setPosition(pos);
    asObjs(i).toggleTitleAsImageText;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function combo_test07 ( )

%*****************************************************************************80
%
%% COMBO_TEST07 tests I4_FACTORIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMBO_TEST07:\n' );
  fprintf ( 1, '  I4_FACTORIAL evaluates the factorial function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     X       Exact F       FACTORIAL(X)\n' );
  fprintf ( 1, '\n' );

  n = 0;

  while ( 1 )

    [ n, x, fx ] = i4_factorial_values ( n );

    if ( n == 0 )
      break
    end

    if ( x <= 0.0 )
      continue
    end

    fx2 = i4_factorial ( x );

    fprintf ( 1, '  %4d  %12d  %12d\n', x, fx, fx2 );

  end

  return
end


text:
function axisset(mode)
% AXISSET User axis setting.

% Author: Craig Borghesani <cborg@terasoft.com>
% Date: 1/19/98 10:36PM

% obtain handle information
if mode == 1,
 f = gcf;

else
 f2 = gcf;
 f = get(f2,'userdata');

end

ui_data = get(f,'userdata');
ui_han = ui_data{1};
last_lims = ui_han(36);
ui_han2 = get(ui_han(54),'userdata');

cur_axs = get(ui_han(32),'userdata');

axs_data = get(cur_axs,'userdata');
axs_typ = axs_data(1);

if mode == 1, % axis limit entering

 win_str = get(ui_han(48),'userdata');

% define various shades of grey
 grey = get(0,'defaultuicontrolbackground');
 ltgrey = [0.5,0.5,0.5]*1.5;
 dkgrey = [0.5,0.5,0.5]*0.5;

 f2 = findobj(0,'type','figure','tag','AxisLimits');

 if isempty(f2),
  scrn_size = get(0,'screensize');
  fig_h = 80;
  fig_w = 360;
  fig_lft = (scrn_size(3) - fig_w)/2;
  fig_btm = (scrn_size(4) - (fig_h+20))/2;
  f2 = figure('pos',[fig_lft,fig_btm,fig_w,fig_h],'menubar','none','numbertitle','off',...
              'color',grey,'name',['Axis Limits',win_str],...
              'userdata',f,'defaultuicontrolbackgroundcolor',grey,...
              'tag','AxisLimits');

  lft = 10; btm = fig_h - 25;
  ui_han2(1) = uicontrol('style','text',...
                        'pos',[lft,btm,80,17],...
                        'string','Current Axis:',...
                        'horiz','left');
  ui_han2(2) = uicontrol('style','text',...
                        'pos',[lft+85,btm,150,17],...
                        'string','',...
                        'horiz','left');

  btm = btm - 22;
  ui_han2(3) = uicontrol('style','text',...
                        'pos',[lft,btm,50,17],...
                        'string','X Min',...
                        'horiz','left');
  ui_han2(4) = uicontrol('style','edit',...
                        'pos',[lft+55,btm,60,20],...
                        'string','',...
                        'horiz','right',...
                        'back','w');
  ui_han2(5) = uicontrol('style','text',...
                        'pos',[lft+130,btm,50,17],...
                        'string','X Max',...
                        'horiz','left');
  ui_han2(6) = uicontrol('style','edit',...
                        'pos',[lft+185,btm,60,20],...
                        'string','',...
                        'horiz','right',...
                        'back','w');
  btm = btm - 22;
  ui_han2(7) = uicontrol('style','text',...
                        'pos',[lft,btm,50,17],...
                        'string','Y Min',...
                        'horiz','left');
  ui_han2(8) = uicontrol('style','edit',...
                        'pos',[lft+55,btm,60,20],...
                        'string','',...
                        'horiz','right',...
                        'back','w');
  ui_han2(9) = uicontrol('style','text',...
                        'pos',[lft+130,btm,50,17],...
                        'string','Y Max',...
                        'horiz','left');
  ui_han2(10) = uicontrol('style','edit',...
                        'pos',[lft+185,btm,60,20],...
                        'string','',...
                        'horiz','right',...
                        'back','w');

  btm = fig_h - 25;
  lft = fig_w - 70;
  ui_han2(11) = uicontrol('style','push',...
                         'pos',[lft,btm,60,20],...
                         'string','OK',...
                         'callback','axisset(2)');

  btm = btm - 25;
  ui_han2(12) = uicontrol('style','push',...
                         'pos',[lft,btm,60,20],...
                         'string','Cancel',...
                         'callback','axisset(3)');

  set(ui_han(54),'userdata',ui_han2);
  set(f2,'vis','on');

 end

 figure(f2);

 if axs_data(1) == 1, set(ui_han2(2),'string','Nyquist');
 elseif axs_data(1) == 2, set(ui_han2(2),'string','Bode (Magnitude)');
 elseif axs_data(1) == 3, set(ui_han2(2),'string','Bode (Phase)');
 elseif axs_data(1) == 4, set(ui_han2(2),'string','Nichols');
 elseif axs_data(1) == 5, set(ui_han2(2),'string','Root Locus');
 elseif axs_data(1) == 6, set(ui_han2(2),'string','Gain (Magnitude)');
 elseif axs_data(1) == 7, set(ui_han2(2),'string','Gain (Phase)');
 elseif axs_data(1) == 8, set(ui_han2(2),'string','Time Response'); end

 Xlimits = get(cur_axs,'xlim');
 Ylimits = get(cur_axs,'ylim');
 set(ui_han2(4),'string',Xlimits(1));
 set(ui_han2(6),'string',Xlimits(2));
 set(ui_han2(8),'string',Ylimits(1));
 set(ui_han2(10),'string',Ylimits(2));

elseif mode == 2, % axis limits acceptance

 Xlimits = get(cur_axs,'xlim');
 Ylimits = get(cur_axs,'ylim');
 set(last_lims,'userdata',[cur_axs,Xlimits,Ylimits]);

 Xmin = str2num(get(ui_han2(4),'string'));
 Xmax = str2num(get(ui_han2(6),'string'));
 Ymin = str2num(get(ui_han2(8),'string'));
 Ymax = str2num(get(ui_han2(10),'string'));

 if any(axs_typ == [1, 5]),
  set(cur_axs,{'xlim','ylim'},ctsqaxes(cur_axs,[Xmin,Xmax],[Ymin,Ymax]));
 else
  set(cur_axs,'xlim',[Xmin,Xmax],'ylim',[Ymin,Ymax]);
 end

 set(f2,'vis','off');

elseif mode == 3, % Cancel

  set(f2,'vis','off');

end

author: Sable

repo: mcbench-benchmarks


text:
function [Knots_n] = refine_cubic_grid_2d(Knots, old_spacing, old_vsz, new_spacing, new_vsz, varargin)
    %[Knots_n] = refine_linear_grid_3d(Knots, old_spacing, ds, volsz,
    %  {upsampling_type}
    Nd = size(Knots, 4);

    upsampling_type = 'sample_std';
    if nargin >= 6
        upsampling_type = varargin{1};
    end

    ksz_new = ceil(new_vsz(1:Nd) ./ new_spacing(1:Nd)) + 3;
    k = (old_vsz(1:Nd)) ./ (new_vsz(1:Nd));
    d = 2 * (1-k);
    interp_type = 1;

    tmp = cat(3, volresize_kd(squeeze(Knots(:,:,1)), k,d, interp_type), ...
                 volresize_kd(squeeze(Knots(:,:,2)), k,d, interp_type));
    if any(size(tmp) < [ksz_new,2])
        warning('Knots upsampling unexpected behavior');
        tmp2 = zeros([ksz_new, 2]);
        szs = min([size(tmp,1), size(tmp,2)], ksz_new);
        tmp2(1:szs(1), 1:szs(2),  :) = tmp(1:szs(1), 1:szs(2),  :);
        tmp = tmp2;
    end
    Kn = tmp(1:ksz_new(1), 1:ksz_new(2), :);

    Knots_n = Kn;
    if strcmp(upsampling_type, 'variation')
        Tmin = cubic_disp_2d(squeeze(Knots), old_vsz, old_spacing);
        Tmin_u = cat(3, imresize_my(Tmin(:,:, 1), new_vsz, 1), ...
                        imresize_my(Tmin(:,:, 2), new_vsz, 1));
        objf = @(x) align_knots(x, Tmin_u, new_vsz, [ksz_new, 2], new_spacing);
        uopt = []; uopt.method = 'cg';
        uopt.MaxIter = 10;
        uopt.Corr = 15;
        uopt.Display = 'off';
        uopt.DerivativeCheck = 'off';
        uopt.LS_type = 0; uopt.LS_init = 8;
        K2 = minFunc(objf, Kn(:), uopt);
        Knots_n = reshape(K2, [ksz_new, 2]);
    end
end

function [f, gr] = align_knots(K, T, szv, szk, grid_spacing)
    K = reshape(K, szk);
    Kx = cubic_disp_2d(K, szv, grid_spacing);
    df = Kx - T;
    f = sum(df(:).^2)/2;
    [gr1, gr2] = cubic_partial_conv_2d(df(:,:, 1), df(:,:, 2), size(K), grid_spacing);
    gr = [gr1(:); gr2(:)];
end

repo: pTVreg


text:
function f21=f21(x)
Bound=[0 10];

if nargin==0
    f21 = Bound;
else
aij=[4 4 4 4;
     1 1 1 1;
     8 8 8 8;
     6 6 6 6;
     3 7 3 7];
ci=[.1 .2 .2 .4 .4];

f21=0;

for i=1:5
    f21=f21+(norm(x-aij(i,:)').^2+ci(:,i))^-1;
end

f21=-f21;

end

author: sfvsfv

repo: Mathematical-modeling


text:
function dn = jed_to_datenum ( jed )

%*****************************************************************************80
%
%% JED_TO_DATENUM converts a JED to a MATLAB date number.
%
%  Discussion:
%
%    The MATLAB "datenum" function accepts a string defining
%    a date and returns a datenumber:
%
%      dn = datenum ( 'Aug 17 1939' )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, real DN, a MATLAB date number.
%
  dn = jed - 1721058.5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef MivFrameworkAppDef < handle
    % MivFrameworkAppDef. Defines application-dependent behaviour for the
    % Framework
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        Config
        ContextDef
        Directories
        ClassFactory
    end

    methods
        function obj = MivFrameworkAppDef
            obj.Config = MimConfig();
            obj.ContextDef = MivContextDef();
            obj.Directories = MimDirectories(obj.GetApplicationParentDirectory, obj.Config);
            obj.ClassFactory = MimClassFactory();
        end

        function context_def = GetContextDef(obj)
            context_def = obj.ContextDef;
        end

        function class_factory = GetClassFactory(obj)
            class_factory = obj.ClassFactory;
        end

        function parent_directory = GetApplicationParentDirectory(obj)
            parent_directory = CoreDiskUtilities.GetUserDirectory;
        end

        function bin_directory = GetBinDirectory(obj)
            bin_directory = fullfile(CoreDiskUtilities.GetUserDirectory, 'mim', 'bin');
        end

        function files_to_compile = GetFilesToCompile(obj, reporting)
            files_to_compile = MivGetMexFilesToCompile(reporting);
        end

        function config = GetFrameworkConfig(obj)
            config = obj.Config;
        end

        function directories = GetFrameworkDirectories(obj)
            directories = obj.Directories;
        end

        function debug_mode = IsDebugMode(obj)
            debug_mode = false;
        end

        function time_functions = TimeFunctions(obj)
            time_functions = obj.Config.TimeFunctions;
        end

        function NewDatasetLoaded(obj, dataset_uid, dataset, reporting)
        end

        function plugins_path = GetUserPath(~)
            full_path = mfilename('fullpath');
            [path_root, ~, ~] = fileparts(full_path);
            plugins_path = fullfile(path_root, '..', 'User');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [results] = niftiCheckIsEqual(nifti1,nifti2,quiet)
%
% results = niftiCheckIsEqual([nifti1],[nifti2],[quiet])
%
%
% DESCRIPTION:
%   Run 'isequal' on a pair of nifti files, field by field. If ~ quiet print
%   a summary of equal and ~equal values do the command window.
%
% USAGE:
%   results = niftiCheckIsEqual(nifti1,nifti2)
%
% INPUT:
%   nifti1  - The filename/path or structure of a niftifile
%   nifti2  - The filename/path or structure of a niftifile
%   quiet   - Supress output to command window. Only field names are
%             returned.
%
% OUTPUT:
%   results - A structure containing the following fields
%           - .allequal >> boolean where 1 == all equal and 0 == not.
%           - .equal    >> fieldnames of equal-fields
%           - .notequal >> filednames of unequal-fields
%           - .nonfield >> fieldnames that are absent in at least one
%                          nifti.
%
%
% (C) Stanford University, VISTA 2014 [ lmperry@stanford.edu ]
%
%
%#ok<*AGROW>

%% HANDLE INPUTS

results = {};

if notDefined('nifti1')
    nifti1 = mrvSelectFile('r','*.nii*','Select first nifti file',pwd);
    if isnumeric(nifti1) || isempty(nifti1)
        return
    end
end

if notDefined('nifti2')
    nifti2 = mrvSelectFile('r','*.nii*','Select second nifti file',pwd);
    if isnumeric(nifti2) || isempty(nifti2)
        return
    end
end

% Read the nifti files if not passed in as structs
if ~isstruct(nifti1)
    n1 = niftiRead(nifti1);
else
    n1 = nifti1;
end

if ~isstruct(nifti2)
    n2 = niftiRead(nifti2);
else
    n2 = nifti2;
end

% Quiet param will cause window output to be supressed
if notDefined('quiet')
    quiet = 0;
end

if ~quiet
    fprintf('Running %s ... \n\tNifti 1 = %s\n\tNifti 2 = %s\n',mfilename,n1.fname,n2.fname);
end

%% COMPARE NIFTI FIELDS

fields1 = fieldnames(n1);
fields2 = fieldnames(n2);

if ~isequal(fields1,fields2)
    warning('Nifti files do not have the same fields.');
    allfields = unique(sort(horzcat(fields1,fields2)));
else
    allfields = fields1;
end

% Initialize the return structure fields
notequal = {};
equal    = {};
nonfield = {};

% Go over the fields and compare
for ii = 1:numel(allfields)
    if isfield(n1, allfields{ii}) && isfield(n2, allfields{ii})
        if isequal(n1.(allfields{ii}),n2.(allfields{ii}))
            equal{end+1} = allfields{ii};
        else
            notequal{end+1} = allfields{ii};
        end
    else
       nonfield{end+1} = allfields{ii};
    end
end

%% DISPLAY THE RESULTS

if ~quiet
    if ~isempty(nonfield)
        fprintf('\nNON-EXISTING FIELDS:\n');
        for jj = 1:numel(nonfield)
            fprintf('\t%s \n',nonfield{jj});
        end
    end

    % Display results for equal fields
    if ~isempty(equal)
        fprintf('\nEQUAL FIELDS:\n');
        for jj = 1:numel(equal)
            fprintf('\t%s \n',equal{jj});
        end

        fprintf('\n:::: EQUAL DETAILED OUTPUT ::::\n')
        for jj = 1:numel(equal)
            disp(['>>> ' equal{jj} ' <<<']);

            % Do not display the data field
            if ~strcmpi(equal{jj},'data')
                disp(n1.(equal{jj}));
                disp('----');
            else disp(' Data fields are equal.');
                disp(' ----')
            end
        end
    end

    % Display results for non-equal fields
    if ~isempty(notequal)
        fprintf('\nNON-EQUAL FIELDS:\n');
        for jj = 1:numel(notequal)
            fprintf('\t%s \n',notequal{jj});
        end

        fprintf('\n:::: NON-EQUAL DETAILED OUTPUT ::::\n')
        for kk = 1:numel(notequal)
            disp(['>>> ' notequal{kk} ' <<<']);

            % Do not display the data field
            if ~strcmpi(notequal{kk},'data')
                disp(n1.(notequal{kk}));
                disp(n2.(notequal{kk}));
                disp('----');
            else disp(' Data field is not equal. Generating difference image.');
                disp(' ----')
                showMontage(n1.data-n2.data);
                set(gcf,'Name','Nifti DATA Difference Image (nifti1.data - nifti2.data)');
            end
        end
    else
        disp('ALL FIELDS ARE EQUAL');
    end
end

%% RETURN THE RESULTS

results.allequal = isequal(n1,n2);
results.equal    = equal;
results.notequal = notequal;
results.nonfield = nonfield;

if ~quiet
    fprintf('\nSUMMARY:\n\t%s fields were equal\n\t%s fields were not equal \n\t%s  fields were unique to one nifti\n\n',num2str(numel(equal)),num2str(numel(notequal)),num2str(numel(nonfield)));
end

return

author: vistalab

repo: vistasoft


text:
function f = hsvargplvmObjective(params, model)

% HSVARGPLVMOBJECTIVE Wrapper function for hierarchical var-GP-LVM objective.
% FORMAT
% DESC provides a wrapper function for the varihierarchical var-GP-LVM, it
% takes the negative of the log likelihood, feeding the parameters
% correctly to the model.
% ARG params : the parameters of the variational GP-LVM model.
% ARG model : the model structure in which the parameters are to be
% placed.
% RETURN f : the negative of the log likelihood of the model.
%
% SEEALSO : hsvargplvmCreate, hsvargplvmLogLikelihood, hsvargplvmExpandParam
%
% COPYRIGHT : Andreas C. Damianou 2012

% HSVARGPLVM

model = modelExpandParam(model, params);
f = - modelLogLikelihood(model);

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function net = concat_2net(net1,net2)
net1 = net1.saveobj();
net2 = net2.saveobj();
net.vars = [net1.vars,net2.vars];
net.layers = [net1.layers,net2.layers];
net.params = [net1.params,net2.params];
net.meta = net1.meta;
net = dagnn.DagNN.loadobj(net);
end

author: layumi

repo: Image-Text-Embedding


text:
function ta = scatter3_2(X,str)
% SCATTER3_2 - scatter plot 3 dimensional data
% (slight modification of MATLAB's scatter3.m)
%
% First argument - input is an nx3 matrix. output is a scatter plot graph, where each point
% comes from one of the rows of the input matrix. The matrix must have 3
% columns.
%
% Optional second argument - string to determine plotting style, using same
% options as in MATLAB's plot function
%
% EXAMPLE:  >> A = rand(20,3);
%           >> scatter3_2(A);
%

% Motivation: MATLAB's scatter3 requires 3 column vectors, eg, scatter3(x,y,z), where
% each 'row' of (x,y,z) coordinates is plotted. We would instead like to be
% able to just input one nx3 matrix. It will have n points to plot, and
% each point is found in each row. The columns do not need to be separated
% into three different arguments.

if nargin<2
    str='b';
end

scatter3(X(:,1),X(:,2),X(:,3),str);


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters all_theta
%   p = PREDICT(all_theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(all_theta'*x) >= 0.5, predict 1)

m = size(X, 1);
num_labels = size(all_theta, 2);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

C = sigmoid(X*all_theta');
[M , p] = max(C , [] , 2);

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
% HEIGHT = maxPyrHt(IMSIZE, FILTSIZE)
%
% Compute maximum pyramid height for given image and filter sizes.
% Specifically: the number of corrDn operations that can be sequentially
% performed when subsampling by a factor of 2.

% Eero Simoncelli, 6/96.

function height = maxPyrHt(imsz, filtsz)

imsz = imsz(:);
filtsz = filtsz(:);

if any(imsz == 1) % 1D image
  imsz = prod(imsz);
  filtsz = prod(filtsz);
elseif any(filtsz == 1)              % 2D image, 1D filter
  filtsz = [filtsz(1); filtsz(1)];
end

if any(imsz < filtsz)
  height = 0;
else
  height = 1 + maxPyrHt( floor(imsz/2), filtsz );
end

repo: SelfExSR


text:
function [Pts, varargout] = Glo2Loc( pcData, varargin )
[ Dim, DataLen] = size(pcData);
if Dim ~= 2 && Dim ~= 3
    error('the input pcData must be 2 * N or 3 * N format!');
end
if (nargin - 1) == 2
    R = varargin{1};
    T = varargin{2};
end
if (nargin - 1) == 1
    Pose = varargin{1};
    [R T] = ExtractRT(Pose);
end
Pts = R * ( bsxfun( @minus, pcData, T ) );
if (nargout-1) == 1
    varargout{1} = [inv(R) -R*T];
end
if (nargout-1) == 2
    varargout{1} = inv(R);
    varargout{2} = -R*T;
end
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function [curInputs, curLabel] = PrepareInputFeatures(curImgsLDR, curExpo, curLabel, isTest)

if(~exist('isTest', 'var'))
    isTest = false;
end

numImgs = length(curImgsLDR);

curInLDR = curImgsLDR;

%%% handling the boundary regions where flow is not valid
nanInds1 = isnan(curInLDR{1});
curInLDR{1}(nanInds1) = LDRtoLDR(curInLDR{2}(nanInds1), curExpo(2), curExpo(1));

nanInds2 = isnan(curInLDR{3});
curInLDR{3}(nanInds2) = LDRtoLDR(curInLDR{2}(nanInds2), curExpo(2), curExpo(3));

if(~isTest)
    darkRef = curInLDR{2} < 0.5;
    badRef = (darkRef & nanInds2) | (~darkRef & nanInds1);
    curLabel(badRef) = (LDRtoHDR(curInLDR{2}(badRef), curExpo(2)));
end

%%% concatenating inputs in the LDR and HDR domains
curLabel = (curLabel);
curInputs = (curInLDR{1});
for k = 2 : numImgs
    curInputs = cat(3, curInputs, (curInLDR{k}));
end

for k = 1 : numImgs
    curInputs = cat(3, curInputs, (LDRtoHDR(curInLDR{k}, curExpo(k))));
end

author: qingsenyangit

repo: AHDRNet