text:
clear; clc; close all;

for repeat = 1:5
close all;

dataset = 'movielens10m';
load(strcat('../../noncvx-lowrank(ICDM 2015)/data/recsys/', dataset, '.mat'));
% load('matlab.mat');

clear item user Movie;

[row, col, val] = find(data);
idx = randperm(length(val));

val = val - mean(val);
val = val/std(val);

traIdx = idx(1:floor(length(val)*0.5));
tstIdx = idx(ceil(length(val)*0.5): end);

clear idx;

traData = sparse(row(traIdx), col(traIdx), val(traIdx), size(data,1), size(data,2));

para.maxIter = 100000;
para.tol = 1e-5;

para.test.row  = row(tstIdx);
para.test.col  = col(tstIdx);
para.test.data = val(tstIdx);
para.test.m = size(traData, 1);
para.test.n = size(traData, 2);

% ---------------------------------------------------------------
switch (dataset)
    case 'movielens100k'
        lambdaMax = 30;
        gridLambda = lambdaMax*(0.9).^(0:9);
    case 'movielens1m'
        lambdaMax = 35;
        gridLambda = lambdaMax*(0.95).^(0:9);
    case 'movielens10m'
        gridLambda = 60;
end

clear lambdaMax;

gridRMSE = zeros(2, size(gridLambda,2 ));
gridRank = zeros(1, size(gridLambda,2 ));
for g = 1:length(gridLambda)
    lambda = gridLambda(g);

    [U, S, V] = SoftImpute(traData, lambda, para );
    [U, S, V, out] = AISImpute(traData, lambda, para );
    gridRank(g) = nnz(S);

    gridRMSE(1, g) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

    [ U, S, V ] = PostProcess(traData, U, V, S);

    gridRMSE(2, g) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

    if(g > 1 && gridRMSE(2, g) > gridRMSE(2, g - 1))
        break;
    end
end

gridRMSE = gridRMSE(2,1:g);
[~, lambda] = min(gridRMSE);
gndRank = gridRank(lambda);

switch (dataset)
    case 'movielens100k'
        para.maxR = ceil(gndRank*1.2);
    case 'movielens1m'
        para.maxR = ceil(gndRank*1.2);
    case 'movielens10m'
        para.maxR = ceil(gndRank*1.5);
end

lambda = gridLambda(lambda);

clear gridRMSE gridRank g X U S V gridLambda out data;

% active --------------------------------------------------------
method = 1;
t = tic;
[U, S, V, out{method}] = ActiveSubspace(traData, lambda, para );
Time(method) = toc(t);
RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

[ U, S, V ] = PostProcess(traData, U, V, S);
RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
clear U S V t;

out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));

figure(1);
plot(out{method}.Time, out{method}.RMSE);
hold on;
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
hold on;

% boost ---------------------------------------------------------
method = 2;
t = tic;
[U, S, V, out{method}] = Boost( traData, lambda, para);
Time(method) = toc(t);
Time(method) = toc(t);
RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

[ U, S, V ] = PostProcess(traData, U, V, S);
RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
clear U S V t;

out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

% % TR ------------------------------------------------------------
% method = 3;
% t = tic;
% [U, S, V, out{method}] = MMBS( traData, lambda, para );
% Time(method) = toc(t);
% Time(method) = toc(t);
% RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
%
% [ U, S, V ] = PostProcess(traData, U, V, S);
% RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
% clear U S V t;
%
% out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));
%
% figure(1);
% plot(out{method}.Time, out{method}.RMSE);
% figure(2);
% semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

% ALT-Impute ----------------------------------------------------
method = 4;
t = tic;
[U, S, V, out{method}] = SoftImputeALS( traData, lambda, para.maxR, para );
Time(method) = toc(t);
RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

[ U, S, V ] = PostProcess(traData, U, V, S);
RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
clear U S V t;

out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

% % SSGD ----------------------------------------------------------
% method = 5;
% t = tic;
% [U, S, V, out{method}] = SSGD( traData, lambda, gndRank, para );
% Time(method) = toc(t);
% RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
%
% [ U, S, V ] = PostProcess(traData, U, V, S);
% RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
% clear U S V t;
%
% out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));
%
% figure(1);
% plot(out{method}.Time, out{method}.RMSE);
% figure(2);
% semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

% LMaFit --------------------------------------------------------
% gridRMSE = zeros(1, 10);
% for g = 1:10
%     [U, S, V] = FixedRank( traData, g, para );
%
%     gridRMSE(g) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
%
% %     if(g > 1 && gridRMSE(g) > gridRMSE(g - 1))
% %         break;
% %     end
% end
%
% gridRMSE = gridRMSE(1:g);
% [~, rnk] = min(gridRMSE);
%
% clear gridRMSE X g U S V;

switch (dataset)
    case 'movielens100k'
        rnk = 3;
    case 'movielens1m'
        rnk = 7;
    case 'movielens10m'
        rnk = 12;
end

method = 6;
t = tic;
[U, S, V, out{method}] = FixedRank( traData, rnk, para );
Time(method) = toc(t);

RMSE(2, method) = out{method}.RMSE(end);

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

% %% APG -----------------------------------------------------------
% method = 7;
% t = tic;
% [U, S, V, out{method}] = APGMatComp( traData, lambda, para );
% Time(method) = toc(t);
% RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
%
% [ U, S, V ] = PostProcess(traData, U, V, S);
% RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
% clear U S V t;
%
% out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));
%
% figure(1);
% plot(out{method}.Time, out{method}.RMSE);
% figure(2);
% semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

%% Soft-Impute ---------------------------------------------------
method = 8;
t = tic;
[U, S, V, out{method}] = SoftImpute( traData, lambda, para );
Time(method) = toc(t);
RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

[ U, S, V ] = PostProcess(traData, U, V, S);
RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
clear U S V t;

out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

%% AISImpute -----------------------------------------------------
method = 9;
t = tic;
[U, S, V, out{method}] = AISImpute( traData, lambda, para );
Time(method) = toc(t);
RMSE(1, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));

[ U, S, V ] = PostProcess(traData, U, V, S);
RMSE(2, method) = MatCompRMSE(U, V, S, row(tstIdx), col(tstIdx), val(tstIdx));
clear U S V t;

out{method}.RMSE = out{method}.RMSE - (RMSE(1, method) - RMSE(2, method));

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));

% clear col data gndRank lambda method para rnk row traData traIdx tstIdx val Age movies;

figure(1);
xlabel('time (sec)');
ylabel('testing RMSE');
legend('active', 'boost', 'TR', 'ALT-Impute', 'SSGD', 'LaMFit', 'APG', 'Soft-Impute', 'AIS-Impute')

figure(2);
xlabel('time (sec)');
ylabel('relative error(obj)');
legend('active', 'boost', 'TR', 'ALT-Impute', 'SSGD', 'LaMFit', 'APG', 'Soft-Impute', 'AIS-Impute')

clear traData lambda gndRank method row S t traIdx tstIdx U V val col rnk para;

save(strcat(dataset, '-', num2str(repeat),'.mat'));

end

author: HKUST-KnowComp

repo: FMG

save_path: github-repos/MATLAB/HKUST-KnowComp-FMG


text:
function xf = compact_fourier_coeff(xf)

% Creates a compact fourier series representation by removing the strict right
% half plane.

if iscell(xf)
    xf = cellfun(@(xf) xf(:,1:(size(xf,2)+1)/2,:), xf, 'uniformoutput', false);
else
    xf = xf(:,1:(size(xf,2)+1)/2,:);
end

author: martin-danelljan

repo: ECO


text:
%% script demoing how to run JLabel

clear classes;
%% set up path

[~,computername] = system('hostname');
computername = strtrim(computername);

switch computername,

  case 'bransonk-lw1',

    JCtrax_path = 'E:\Code\JCtrax';
    %FlyBowlAnalysis_path = 'E:\Code\FlyBowlAnalysis';
    %rootdatadir = 'E:\Code\Jdetect\larva\fly_data\TrainingData';
    %settingsdir = 'E:\Code\FlyBowlAnalysis\settings';
    configfilename = 'params\JLabelParams.xml';
    npool = 4;

  case 'bransonk-lw2',

    JCtrax_path = 'C:\Code\JCtrax';
    %FlyBowlAnalysis_path = 'C:\Code\FlyBowlAnalysis';
    %rootdatadir = 'C:\Code\Jdetect\larva\fly_data\TrainingData';
    %settingsdir = 'C:\Code\FlyBowlAnalysis\settings';
    configfilename = 'params\JLabelParamsKristinChase.xml';
    npool = 4;

  case 'bransonk-desktop',

    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/larva/fly_data/TrainingData/';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = 'params/JLabelParams_bransonk-desktop.xml';
    npool = 8;

  case 'kabram-ws.janelia.priv',

    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/larva/fly_data/TrainingData/';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = 'params/JLabelParamsMayank.xml';
    npool = 8;

  case 'robiea-ws'
    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/sciserv/flyolympiad/Olympiad_Screen/fly_bowl/bowl_data';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = '/groups/branson/home/robiea/Projects_data/JLabel/JLabelParams.xml';
    npool = 8;

  otherwise

    warning('Unknown computer %s. Paths may not be setup correctly',computername);
    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/larva/fly_data/TrainingData/';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = 'params/JLabelParams_bransonk-desktop.xml';
    npool = 1;

end

addpath(fullfile(JCtrax_path,'misc'));
addpath(fullfile(JCtrax_path,'filehandling'));
jlabelpath = fileparts(which('JLabel'));
addpath(fullfile(jlabelpath,'compute_perframe_features'));
%addpath(FlyBowlAnalysis_path);

%%
try
  if npool > 1,
    matlabpool(npool);
  end
catch %#ok<CTCH>
end
%% start JLabel

hfig = JLabel('configfilename',configfilename); handles = guidata(hfig);

author: kristinbranson

repo: JAABA


text:
function [MPa] = hPa2MPa(hPa)
% Convert pressure from hectopascals to megapascals.
% Chad Greene 2012
MPa = hPa*0.000100000;

author: Sable

repo: mcbench-benchmarks


text:
% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix
function make(opt)
fprintf('=> Building liblinear.\n');
if nargin < 1
    try
        % This part is for OCTAVE
        if(exist('OCTAVE_VERSION', 'builtin'))
            % Use -std=c++11 for newer versions of Octave
            if ispc
                setenv('CFLAGS','-std=gnu99 -O3')
                setenv('CC','gcc')
            else
                setenv('CFLAGS','-O3 -fstack-protector-strong -Wformat -Werror=format-security')
            end
            %mex libsvmread.c
            %mex libsvmwrite.c
            mex -I.. -O3 svmtrain.cpp linear_model_matlab.cpp ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
            mex -I.. -O3 svmpredict.cpp linear_model_matlab.cpp ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
            delete *.o
            % This part is for MATLAB
            % Add -largeArrayDims on 64-bit machines of MATLAB
        else
            if ispc
                %mex COMPFLAGS="\$COMPFLAGS -std=c99 -O3" -largeArrayDims libsvmread.c
                %mex COMPFLAGS="\$COMPFLAGS -std=c99 -O3" -largeArrayDims libsvmwrite.c
                mex COMPFLAGS="\$COMPFLAGS -O3" -I.. -largeArrayDims svmtrain.cpp linear_model_matlab.cpp ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
                mex COMPFLAGS="\$COMPFLAGS -O3" -I.. -largeArrayDims svmpredict.cpp linear_model_matlab.cpp ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
            else
                %mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c
                %mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c
                mex CFLAGS="\$CFLAGS" -I.. -largeArrayDims svmtrain.cpp linear_model_matlab.cpp ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
                mex CFLAGS="\$CFLAGS" -I.. -largeArrayDims svmpredict.cpp linear_model_matlab.cpp ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
            end
        end
    catch err
        fprintf('Error: %s failed (line %d)\n', err.stack(1).file, err.stack(1).line);
        disp(err.message);
        fprintf('=> Please check README for detailed instructions.\n');
    end
elseif nargin == 1
    switch lower(opt)
        case 'clean'
            delete *.o
        case 'cleanall'
            delete *.o
            delete *.mexa64
        otherwise
            error('make option "%s" not recognized', opt)
    end
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
%% Pei&Lin Normalization
% This program demonstrates Pei-Lin Normalization.
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.3.1/modules/ximgproc/samples/peilin.cpp>
%

%%
% Source images
fname1 = fullfile(mexopencv.root(), 'test', 'peilin_plane.png');
fname2 = fullfile(mexopencv.root(), 'test', 'peilin_shape.png');
if exist(fname1, 'file') ~= 2
    disp('Downloading Image...')
    url = 'https://cdn.rawgit.com/opencv/opencv_contrib/3.3.1/modules/ximgproc/samples/peilin_plane.png';
    urlwrite(url, fname1);
end
if exist(fname2, 'file') ~= 2
    disp('Downloading Image...')
    url = 'https://cdn.rawgit.com/opencv/opencv_contrib/3.3.1/modules/ximgproc/samples/peilin_shape.png';
    urlwrite(url, fname2);
end

%%
% Load images
I = cv.imread(fname1, 'Grayscale',true);
J = cv.imread(fname2, 'Grayscale',true);

%%
% Apply normalization
N = cv.warpAffine(I, cv.PeiLinNormalization(I));
D = cv.warpAffine(I, cv.PeiLinNormalization(J), 'WarpInverse',true);

%%
% Show results
subplot(221), imshow(I), title('I')
subplot(222), imshow(N), title('N')
subplot(223), imshow(J), title('J')
subplot(224), imshow(D), title('D')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function V = ReferenceVectorAdaptation(PopObj,V)
% Reference vector adaption strategy

%------------------------------- 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".
%--------------------------------------------------------------------------

    V = V.*repmat(max(PopObj,[],1)-min(PopObj,[],1),size(V,1),1);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function geodemo_4b
%% GEODEMO_4b
% Compare horizontal slices from two different CF compliant structured
% grid models (EFDC and CBOFS) at a particular time step and depth, using
% geosubset to subset the data

%% EFDC
url{2}='http://comt.sura.org/thredds/dodsC/data/comt_1_archive/estuarine_hypoxia/VIMS_EFDC/2004_DO3d';
var{2}='salt';
titl{2}='EFDC';

%% ROMS
url{1}='http://comt.sura.org/thredds/dodsC/data/comt_1_archive/estuarine_hypoxia/VIMS_CBOFS/2004-2005';
var{1}='salt';
titl{1}='CBOFS';

%% Create geosubset object
dat=[2004 4 10 6 0 0];  % Apr 4, 2004 06:00 UTC
depth=-5;  % horizontal slice 5 m from surface
ax=[ -76.5220  -75.7105   36.8248   37.7850]; %lon/lat range
cax=[0 33];  %color range
lat_mid=38; % for scaling plots

s.time=datenum(dat);
s.lon=ax(1:2);
s.lat=ax(3:4);

%% Perform analysis without using dataset dependant code
% Access datasets, subset data, interpolate data to a constant z, plot results at z depth

figure;
for i=1:length(url);
  nc{i}=ncgeodataset(url{i});
  % create a salinity geovariable object.  No data read yet.
  svar{i}=geovariable(nc{i},var{i});
  disp(['reading data from ' titl{i} '...'])
  % using geosubset here, which allows for subsetting and striding, reading
  % multiple time steps, only certain zlevels and more.
  sub{i}=svar{i}.geosubset(s);
  sz{i}=zsliceg(squeeze(sub{i}.data),squeeze(sub{i}.grid.z),depth);
  a{i}=subplot(1,length(url),i);
  pcolorjw(sub{i}.grid.lon,sub{i}.grid.lat,double(sz{i}));colorbar
  axis(ax);
  caxis(cax);
  title(sprintf('%s, depth=%f: %s',titl{i},depth,datestr(sub{1}.grid.time)));
  set (a{i}, 'DataAspectRatio', [1 cos(lat_mid*pi/180) 1000] );
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, varargin)
% Plot: data panel
%
% case 'save', dosave = 1;
% case {'savedir', 'mysavedir'}, mysavedir = varargin{i+1};
% case {'figs', 'dotfigs'}, savedotfigs = 1;
%
% Doc not complete yet.  Please update me!
%
% See parcel_clusters.m
% See parcel_cl_nmds.m
%
% load Parcellation_info/parcellation.mat
%
% :Examples:
% ::
%
%    parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, 'save')
%    parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, 'save', 'savedir', 'Parcellation_info')
%    parcel_cl_nmds_plots(parcel_cl_avgs, NMDS, 'save', 'savedir', 'Parcellation_info_tor_mask_try1', 'savedotfig')

    % ..
    %    SETUP
    % ..

    dosave = 0;
    mysavedir = [];
    savedotfigs = 0;

    for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}

                % functional commands
                case 'save', dosave = 1;

                case {'savedir', 'mysavedir'}, mysavedir = varargin{i+1};

                case {'figs', 'dotfigs'}, savedotfigs = 1;

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
    end

    N = length(parcel_cl_avgs(1).timeseries);

    if ~isfield(NMDS, 'basecolors')
        NMDS.basecolors = {[1 0 0] [0 1 0] [0 0 1] [1 1 0] [0 1 1] [1 0 1] ...
            [1 .5 0] [.5 1 0] [.5 0 1] [1 0 .5] [0 1 .5] [0 .5 1]};
    end

    basecolors = NMDS.basecolors;

    % NMDS plot
    % --------------------------------------------------
    % Make figure
    disp('Visualizing results')
    create_figure('nmdsfig');

    nmdsfig(NMDS.stats_mds.GroupSpace,'classes',NMDS.stats_mds.ClusterSolution.classes, ...
        'names', [],'sig',NMDS.stats.fdrsig, 'fill', 'colors', basecolors);
    % --------------------------------------------------

    % --------------------------------------------------
    % Get class clusters (convenient format)
    % --------------------------------------------------

    % re-define class clusters
    clear class_clusters
    for i = 1:max(NMDS.stats_mds.ClusterSolution.classes)
        wh = find(NMDS.stats_mds.ClusterSolution.classes == i);
        class_clusters{i} = parcel_cl_avgs(wh);

        % refine class (network) membership
        [parcel_cl_avgs(wh).from_class] = deal(i);
    end

    nclasses = length(class_clusters); % redefine based on results

    % --------------------------------------------------
    % get index of which parcels (clusters) are in which classes
    % --------------------------------------------------

    parcelindx = cat(1, parcel_cl_avgs.from_class);
    [parcelindx, ivals] = sort(parcelindx);
    cl = parcel_cl_avgs(ivals);

    % --------------------------------------------------
    % Get data
    % --------------------------------------------------

    data = [];
    for i = 1:length(cl)

        % zscore within, or at least center to avoid individual diffs in
        % baseline driving correlations
        for s = 1:N
            cl(i).timeseries{s} = scale(cl(i).timeseries{s});
        end

        data(:, i) = cat(1, cl(i).timeseries{:});
    end

    for s = 1:N-1, len(s) = length(cl(1).timeseries{s}); end
    len = cumsum([1 len]);

    % --------------------------------------------------
    % Panel TS figure
    % --------------------------------------------------

    create_figure('Panel_timeseries', nclasses, 1);
    for i = 1:nclasses
        subplot(nclasses, 1, i);
        hold on;
        mycolor = basecolors{i};

        [lh, fh] = plot_error(data(:, parcelindx == i)');
        mymean = mean(data(:, parcelindx == i)');

        set(lh, 'Color', mycolor);
        set(fh, 'FaceColor', mycolor);
        set(gca, 'YLim', [min(mymean) max(mymean)]);
        axis tight
        hold on
        plot_vertical_line(len);

        if i == 1
            for j = 1:length(len), text(len(j), max(mymean), ['Subj' num2str(j)]); end
        end

        axis off
        drawnow
    end

    % --------------------------------------------------
    % Correlation figure
    % --------------------------------------------------

    % get within-parcel correlations
    %-----------------------------------
    for i = 1:length(cl)
        dd = cat(1, cl(i).all_data{:});
        rr = corrcoef(dd);
        rr = rr .* (1 - eye(size(rr)));
        within_corr(i) = mean(squareform(rr));
    end

    create_figure('Correlations');
    [r, p] = corrcoef(data);

    % put withins on diagonal
    r = r .* (1 - eye(size(r)));
    r = r + diag(within_corr);
    %r(p > .05) = NaN;

    % put spaces in to mark off classes
    %-----------------------------------
    %wh_spaces = [0; diff(parcelindx)];
    rnew = [];
    pinew = [];
    for i = 1:nclasses
        rnew = [rnew; r(parcelindx == i, :)];
        rnew(end + 1, :) = NaN;

        pinew = [pinew; parcelindx(parcelindx == i)];
        pinew(end+1) = NaN;

    end

    rnew2 = [];
    for i = 1:nclasses
        rnew2 = [rnew2 rnew(:, parcelindx == i)];
        rnew2(:, end + 1) = NaN;
    end
    r = rnew2;

    % Image the correlation matrix
    %-----------------------------------
    imagesc(r); colorbar
    set(gca,'YDir', 'Reverse')
    axis tight
    newcm = colormap_tor([0 0 .7], [1 1 0], [1 .5 0]);
    colormap(newcm)
    title('Correlations among parcels');
    xlabel('Parcel index number');

    % Color bars for class ID
    %-----------------------------------
    axpos = get(gca, 'Position');
    axh = axes('Position', [.05 axpos(2) .04 axpos(4)]);
    set(axh, 'YDir', 'Reverse', 'YLim', [1 length(pinew)]);
    hold on;
    for i = 1:nclasses

        yy = find(pinew == i);
        plot(.7 * ones(size(yy)), yy, '-', 'Color', basecolors{i}, 'LineWidth', 10);

        text(-2, round(median(yy)), num2str(i), 'FontSize', 18, 'Color', 'k');

    end
    axis off
    drawnow

    % --------------------------------------------------
    % Orthviews of parcels
    % --------------------------------------------------

    nclasses = length(class_clusters); % redefine based on results

    while length(basecolors) < nclasses
        basecolors{end+1} = basecolors{end - 11} + basecolors{end - 6} ./ 2;
    end

    parcelindx = cat(1, parcel_cl_avgs.from_class);

    colors = nmdsfig_tools('class_colors',parcelindx,basecolors);

    addstr = 'noadd';

    for i = 1:length(parcel_cl_avgs)

        cluster_orthviews(parcel_cl_avgs(i), colors(i), addstr, 'solid');

        addstr = 'add';
    end

    % --------------------------------------------------
    % Montages of parcels
    % --------------------------------------------------

    disp(['Saving montages in ' mysavedir])

    for i = 1:length(class_clusters)

        mycolor = basecolors(i);

        montage_clusters([], class_clusters{i}, mycolor);
        drawnow

        savename = fullfile(mysavedir, ['Montage_class' num2str(i)]);
        saveas(gcf, savename,'png')

    end

    if dosave
        save_figures(1, savedotfigs, mysavedir);
    end

    % --------------------------------------------------
    % Glass brains
    % --------------------------------------------------
    for i = 1:length(class_clusters)
        create_figure('glass');

        mysig = NMDS.stats.fdrsig;
        mysig = mysig(:, parcelindx == i);
        mysig = mysig(parcelindx == i, :);

        cluster_nmdsfig_glassbrain( ...
            class_clusters{i}, ...
            ones(length(class_clusters{i}), 1), ...
            NMDS.basecolors(i), ...
            mysig, [], 'blobs');

        if dosave
           saveas(gcf, fullfile(mysavedir, ['glass_network_' num2str(i)]), 'png');
           if savedotfigs
               saveas(gcf, fullfile(mysavedir, ['glass_network_' num2str(i)]), 'fig');
           end
        end

    end

end

% -------------------------------------------------------------------------
% Save all figures
% -------------------------------------------------------------------------
function save_figures(verbose, savedotfigs, mysavedir)

    if nargin == 0, verbose = 1; end
    if nargin == 1, savedotfigs = 0; end
    if nargin == 2 || isempty(mysavedir), mysavedir = '.'; end

    fignames = {'nmdsfig' 'Panel_timeseries' 'Correlations'};

    for i = 1:length(fignames)
        name = fignames{i};

        h = findobj('Tag', name);
        if ~isempty(h)

            if length(h) > 1
                if verbose, disp(['Warning: More than one plot with tag ' name '. Saving highest-numbered figure.']); end
                h = max(h);
            end

            figure(h)
            scn_export_papersetup(500);

            if savedotfigs
                if verbose
                    fprintf('Saving: %s%s\n', fullfile(mysavedir, name), '.fig');
                end

                saveas(h, fullfile(mysavedir, name), 'fig');
            end

            if verbose
                fprintf('Saving: %s%s\n', fullfile(mysavedir, name), '.png');
            end
            saveas(h, fullfile(mysavedir, name), 'png');
        end
    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function F = cumprod(F)
%CUMPROD   Indefinite product integral.
%   CUMPROD(F) is the indefinite product integral of the CHEBFUN F, which is
%   defined as exp(cumsum(log(F))).
%
% See also PROD.

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

% Loop over the columns:
for k = 1:numel(F)
    %[TODO]: Is this right? What about row CHEBFUNS?
    F(k) = exp(cumsum(log(F(k))));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function I_d = Akashi2016(I)
%Akashi2016 I_d = Akashi2016(I)
%  The core mechanism is herein implemented fully in accordance with their
%  paper. I have not implemented the repeated score trials. It is pending.
%  This is because this method is somewhat slow.
%
%  Convergence criterion is based on Yamamoto and Nakazawa (exp(-15)), you
%  can change it:
%    edit Akashi2016.m
%  Search 'eps' and change it.
%
%  Also, because of random initialization, I believe averaging the
%  estimated diffuse component over many applications of the method will
%  greatly improve the quality.
%
%  See also SIHR.

assert(isa(I, 'float'), 'SIHR:I:notTypeSingleNorDouble', ...
    'Input I is not type single nor double.')
assert(min(I(:)) >= 0 && max(I(:)) <= 1, 'SIHR:I:notWithinRange', ...
    'Input I is not within [0, 1] range.')
[n_row, n_col, n_ch] = size(I);
assert(n_row > 1 && n_col > 1, 'SIHR:I:singletonDimension', ...
    'Input I has a singleton dimension.')
assert(n_ch == 3, 'SIHR:I:notRGB', ...
    'Input I is not a RGB image.')

I = 255 * I;
M = 3; % = number of color channels
N = n_row * n_col; % = number of pixels
R = 7; % = number of color clusters - 1
I = reshape(I, [N, M])';
i_s = ones(3, 1, class(I)) / sqrt(3);
H = 254 * rand(R, N, class(I)) + 1;
W_d = 254 * rand(3, R-1, class(I)) + 1;
W_d = my_normc(W_d); % W_d ./ vecnorm(W_d, 2, 1); % normalize(W_d,1,'norm');
W = [i_s, W_d];
A = ones(M, class(I));
lambda = 3;
eps = exp(-15);
F_t_1 = Inf;
iter = uint16(0);
max_iter = uint16(10e3); % change for later convergence (takes way longer)
% tic
while true
    W_bar = my_normc(W);
    H = H .* ((W_bar') * I) ./ ...
        ((W_bar') * W_bar * H + lambda);
    H_d = H(2:end, :);
    Vl = max(0, I-i_s*H(1, :));
    % W_d = W(:,2:end); % not sure
    W_d_bar = W_bar(:, 2:end);
    W_d = W_d_bar .* (Vl * (H_d') + W_d_bar .* (A * W_d_bar * H_d * (H_d'))) ./ ...
        (W_d_bar * H_d * (H_d') + W_d_bar .* (A * Vl * (H_d')));
    W = [i_s, W_d];
    F_t = 0.5 * norm((I - W * H), 'fro') + lambda * sum(H(:));
    if abs(F_t-F_t_1) < eps * abs(F_t) || iter >= max_iter
        break
    end
    F_t_1 = F_t;
    iter = iter + 1;
end
% toc
W_d = W(:, 2:end);
% h_s = H(1, :);
H_d = H(2:end, :);
% I_s = i_s * h_s;
I_d = W_d * H_d;
% I_s = reshape(full(I_s)', [n_row, n_col, n_ch]) / 255;
I_d = reshape(full(I_d)', [n_row, n_col, n_ch]) / 255;

% figure(1), imshow(I_s)
% figure(2), imshow(I_d)

end

author: vitorsr

repo: SIHR

save_path: github-repos/MATLAB/vitorsr-SIHR


text:
function [sel1, sel2] = match_str(a, b, fullout)

% MATCH_STR looks for matching labels in two lists of strings
% and returns the indices into both the 1st and 2nd list of the matches.
% They will be ordered according to the first input argument.
%
% Use as
%   [sel1, sel2] = match_str(strlist1, strlist2)
%
% The strings can be stored as a char matrix or as an vertical array of
% cells, the matching is done for each row.
%
% When including a 1 as the third input argument, the output lists of
% indices will be expanded to the size of the largest input argument.
% Entries that occur only in one of the two inputs will correspond to a 0
% in the output, in this case. This can be convenient in rare cases if the
% size of the input lists is meaningful.

% Copyright (C) 2000-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% ensure that both are cell-arrays
if isempty(a)
  a = {};
elseif ~iscell(a)
  a = cellstr(a);
end
if isempty(b)
  b = {};
elseif ~iscell(b)
  b = cellstr(b);
end

% regardless of what optimizations are implemented, the code should remain
% functionally compatible to the original, which is
% sel1 = [];
% sel2 = [];
% for i=1:length(a)
%   for j=1:length(b)
%     if strcmp(a(i),b(j))
%        sel1 = [sel1; i];
%        sel2 = [sel2; j];
%      end
%    end
% end

% ensure that both are column vectors
a = a(:);
b = b(:);
Na = numel(a);
Nb = numel(b);

% According to the original implementation empty numeric elements are
% allowed, but are not returned as match. This is different to empty string
% elements, which are returned as match.
% See also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1808
empty_a = cellfun(@isnumeric, a) & cellfun(@isempty, a);
empty_b = cellfun(@isnumeric, b) & cellfun(@isempty, b);
% the following allows the unique function to operate normally
a(empty_a) = {''};
b(empty_b) = {''};

% replace all unique strings by a unique number and use the fact that
% numeric comparisons are much faster than string comparisons
[dum1, dum2, c] = unique([a; b]);
a = c(1:Na);
b = c((Na+1):end);

% empty numeric elements should never be returned as a match
a(empty_a) = nan;
b(empty_b) = nan;

if nargin < 3 || ~fullout
  sel1 = [];
  sel2 = [];
  for i=1:length(a)
    % s = find(strcmp(a(i), b));  % for string comparison
    s = find(a(i)==b);            % for numeric comparison
    sel2 = [sel2; s];
    s(:) = i;
    sel1 = [sel1; s];
  end
else
  sel1 = zeros(max(Na,Nb),1);
  sel2 = zeros(max(Na,Nb),1);
  for i=1:length(a)
    s = find(a(i)==b);
    sel2(s) = s;
    sel1(s) = i;
  end
end

author: spm


text:
function [xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk,dxpdalpha] = project_points2(X,f,c,k,alpha)

%project_points2.m
%
%[xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk] = project_points2(X,om,T,f,c,k,alpha)
%
%Projects a 3D structure onto the image plane.
%
%INPUT: X: 3D structure in the world coordinate frame (3xN matrix for N points)
%       (om,T): Rigid motion parameters between world coordinate frame and camera reference frame
%               om: rotation vector (3x1 vector); T: translation vector (3x1 vector)
%       f: camera focal length in units of horizontal and vertical pixel units (2x1 vector)
%       c: principal point location in pixel units (2x1 vector)
%       k: Distortion coefficients (radial and tangential) (4x1 vector)
%       alpha: Skew coefficient between x and y pixel (alpha = 0 <=> square pixels)
%
%OUTPUT: xp: Projected pixel coordinates (2xN matrix for N points)
%        dxpdom: Derivative of xp with respect to om ((2N)x3 matrix)
%        dxpdT: Derivative of xp with respect to T ((2N)x3 matrix)
%        dxpdf: Derivative of xp with respect to f ((2N)x2 matrix if f is 2x1, or (2N)x1 matrix is f is a scalar)
%        dxpdc: Derivative of xp with respect to c ((2N)x2 matrix)
%        dxpdk: Derivative of xp with respect to k ((2N)x4 matrix)
%
%Definitions:
%Let P be a point in 3D of coordinates X in the world reference frame (stored in the matrix X)
%The coordinate vector of P in the camera reference frame is: Xc = R*X + T
%where R is the rotation matrix corresponding to the rotation vector om: R = rodrigues(om);
%call x, y and z the 3 coordinates of Xc: x = Xc(1); y = Xc(2); z = Xc(3);
%The pinehole projection coordinates of P is [a;b] where a=x/z and b=y/z.
%call r^2 = a^2 + b^2.
%The distorted point coordinates are: xd = [xx;yy] where:
%
%xx = a * (1 + kc(1)*r^2 + kc(2)*r^4 + kc(5)*r^6)      +      2*kc(3)*a*b + kc(4)*(r^2 + 2*a^2);
%yy = b * (1 + kc(1)*r^2 + kc(2)*r^4 + kc(5)*r^6)      +      kc(3)*(r^2 + 2*b^2) + 2*kc(4)*a*b;
%
%The left terms correspond to radial distortion (6th degree), the right terms correspond to tangential distortion
%
%Finally, convertion into pixel coordinates: The final pixel coordinates vector xp=[xxp;yyp] where:
%
%xxp = f(1)*(xx + alpha*yy) + c(1)
%yyp = f(2)*yy + c(2)
%
%
%NOTE: About 90 percent of the code takes care fo computing the Jacobian matrices
%
%
%Important function called within that program:
%
%rodrigues.m: Computes the rotation matrix corresponding to a rotation vector
%
%rigid_motion.m: Computes the rigid motion transformation of a given structure

[m,n] = size(X);

Y = X;

inv_Z = 1./Y(3,:);

x = (Y(1:2,:) .* (ones(2,1) * inv_Z)) ;

bb = (-x(1,:) .* inv_Z)'*ones(1,3);
cc = (-x(2,:) .* inv_Z)'*ones(1,3);

if nargout > 1,
    dxdom = zeros(2*n,3);
    dxdom(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(1:3:end,:) + bb .* dYdom(3:3:end,:);
    dxdom(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(2:3:end,:) + cc .* dYdom(3:3:end,:);

    dxdT = zeros(2*n,3);
    dxdT(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(1:3:end,:) + bb .* dYdT(3:3:end,:);
    dxdT(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(2:3:end,:) + cc .* dYdT(3:3:end,:);
end;

% Add distortion:

r2 = x(1,:).^2 + x(2,:).^2;

if nargout > 1,
    dr2dom = 2*((x(1,:)')*ones(1,3)) .* dxdom(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdom(2:2:end,:);
    dr2dT = 2*((x(1,:)')*ones(1,3)) .* dxdT(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdT(2:2:end,:);
end;

r4 = r2.^2;

if nargout > 1,
    dr4dom = 2*((r2')*ones(1,3)) .* dr2dom;
    dr4dT = 2*((r2')*ones(1,3)) .* dr2dT;
end

r6 = r2.^3;

if nargout > 1,
    dr6dom = 3*((r2'.^2)*ones(1,3)) .* dr2dom;
    dr6dT = 3*((r2'.^2)*ones(1,3)) .* dr2dT;
end;

% Radial distortion:

cdist = 1 + k(1) * r2 + k(2) * r4 + k(5) * r6;

if nargout > 1,
    dcdistdom = k(1) * dr2dom + k(2) * dr4dom + k(5) * dr6dom;
    dcdistdT = k(1) * dr2dT + k(2) * dr4dT + k(5) * dr6dT;
    dcdistdk = [ r2' r4' zeros(n,2) r6'];
end;

xd1 = x .* (ones(2,1)*cdist);

if nargout > 1,
    dxd1dom = zeros(2*n,3);
    dxd1dom(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdom;
    dxd1dom(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdom;
    coeff = (reshape([cdist;cdist],2*n,1)*ones(1,3));
    dxd1dom = dxd1dom + coeff.* dxdom;

    dxd1dT = zeros(2*n,3);
    dxd1dT(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdT;
    dxd1dT(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdT;
    dxd1dT = dxd1dT + coeff.* dxdT;

    dxd1dk = zeros(2*n,5);
    dxd1dk(1:2:end,:) = (x(1,:)'*ones(1,5)) .* dcdistdk;
    dxd1dk(2:2:end,:) = (x(2,:)'*ones(1,5)) .* dcdistdk;
end;

% tangential distortion:

a1 = 2.*x(1,:).*x(2,:);
a2 = r2 + 2*x(1,:).^2;
a3 = r2 + 2*x(2,:).^2;

delta_x = [k(3)*a1 + k(4)*a2 ;
    k(3) * a3 + k(4)*a1];

%ddelta_xdx = zeros(2*n,2*n);
aa = (2*k(3)*x(2,:)+6*k(4)*x(1,:))'*ones(1,3);
bb = (2*k(3)*x(1,:)+2*k(4)*x(2,:))'*ones(1,3);
cc = (6*k(3)*x(2,:)+2*k(4)*x(1,:))'*ones(1,3);

if nargout > 1,
    ddelta_xdom = zeros(2*n,3);
    ddelta_xdom(1:2:end,:) = aa .* dxdom(1:2:end,:) + bb .* dxdom(2:2:end,:);
    ddelta_xdom(2:2:end,:) = bb .* dxdom(1:2:end,:) + cc .* dxdom(2:2:end,:);

    ddelta_xdT = zeros(2*n,3);
    ddelta_xdT(1:2:end,:) = aa .* dxdT(1:2:end,:) + bb .* dxdT(2:2:end,:);
    ddelta_xdT(2:2:end,:) = bb .* dxdT(1:2:end,:) + cc .* dxdT(2:2:end,:);

    ddelta_xdk = zeros(2*n,5);
    ddelta_xdk(1:2:end,3) = a1';
    ddelta_xdk(1:2:end,4) = a2';
    ddelta_xdk(2:2:end,3) = a3';
    ddelta_xdk(2:2:end,4) = a1';
end;

xd2 = xd1 + delta_x;

if nargout > 1,
    dxd2dom = dxd1dom + ddelta_xdom ;
    dxd2dT = dxd1dT + ddelta_xdT;
    dxd2dk = dxd1dk + ddelta_xdk ;
end;

% Add Skew:

xd3 = [xd2(1,:) + alpha*xd2(2,:);xd2(2,:)];

% Compute: dxd3dom, dxd3dT, dxd3dk, dxd3dalpha
if nargout > 1,
    dxd3dom = zeros(2*n,3);
    dxd3dom(1:2:2*n,:) = dxd2dom(1:2:2*n,:) + alpha*dxd2dom(2:2:2*n,:);
    dxd3dom(2:2:2*n,:) = dxd2dom(2:2:2*n,:);
    dxd3dT = zeros(2*n,3);
    dxd3dT(1:2:2*n,:) = dxd2dT(1:2:2*n,:) + alpha*dxd2dT(2:2:2*n,:);
    dxd3dT(2:2:2*n,:) = dxd2dT(2:2:2*n,:);
    dxd3dk = zeros(2*n,5);
    dxd3dk(1:2:2*n,:) = dxd2dk(1:2:2*n,:) + alpha*dxd2dk(2:2:2*n,:);
    dxd3dk(2:2:2*n,:) = dxd2dk(2:2:2*n,:);
    dxd3dalpha = zeros(2*n,1);
    dxd3dalpha(1:2:2*n,:) = xd2(2,:)';
end;

% Pixel coordinates:
if length(f)>1,
    xp = xd3 .* (f(:) * ones(1,n))  +  c(:)*ones(1,n);
    if nargout > 1,
        coeff = reshape(f(:)*ones(1,n),2*n,1);
        dxpdom = (coeff*ones(1,3)) .* dxd3dom;
        dxpdT = (coeff*ones(1,3)) .* dxd3dT;
        dxpdk = (coeff*ones(1,5)) .* dxd3dk;
        dxpdalpha = (coeff) .* dxd3dalpha;
        dxpdf = zeros(2*n,2);
        dxpdf(1:2:end,1) = xd3(1,:)';
        dxpdf(2:2:end,2) = xd3(2,:)';
    end;
else
    xp = f * xd3 + c*ones(1,n);
    if nargout > 1,
        dxpdom = f  * dxd3dom;
        dxpdT = f * dxd3dT;
        dxpdk = f  * dxd3dk;
        dxpdalpha = f .* dxd3dalpha;
        dxpdf = xd3(:);
    end;
end;

if nargout > 1,
    dxpdc = zeros(2*n,2);
    dxpdc(1:2:end,1) = ones(n,1);
    dxpdc(2:2:end,2) = ones(n,1);
end;

return;

% Test of the Jacobians:

n = 10;

X = 10*randn(3,n);
om = randn(3,1);
T = [10*randn(2,1);40];
f = 1000*rand(2,1);
c = 1000*randn(2,1);
k = 0.5*randn(5,1);
alpha = 0.01*randn(1,1);

[x,dxdom,dxdT,dxdf,dxdc,dxdk,dxdalpha] = project_points2(X,om,T,f,c,k,alpha);

% Test on om: OK

dom = 0.000000001 * norm(om)*randn(3,1);
om2 = om + dom;

[x2] = project_points2(X,om2,T,f,c,k,alpha);

x_pred = x + reshape(dxdom * dom,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on T: OK!!

dT = 0.0001 * norm(T)*randn(3,1);
T2 = T + dT;

[x2] = project_points2(X,om,T2,f,c,k,alpha);

x_pred = x + reshape(dxdT * dT,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on f: OK!!

df = 0.001 * norm(f)*randn(2,1);
f2 = f + df;

[x2] = project_points2(X,om,T,f2,c,k,alpha);

x_pred = x + reshape(dxdf * df,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on c: OK!!

dc = 0.01 * norm(c)*randn(2,1);
c2 = c + dc;

[x2] = project_points2(X,om,T,f,c2,k,alpha);

x_pred = x + reshape(dxdc * dc,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on k: OK!!

dk = 0.001 * norm(k)*randn(5,1);
k2 = k + dk;

[x2] = project_points2(X,om,T,f,c,k2,alpha);

x_pred = x + reshape(dxdk * dk,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on alpha: OK!!

dalpha = 0.001 * norm(k)*randn(1,1);
alpha2 = alpha + dalpha;

[x2] = project_points2(X,om,T,f,c,k,alpha2);

x_pred = x + reshape(dxdalpha * dalpha,2,n);

norm(x2-x)/norm(x2 - x_pred)

author: jianxiongxiao

repo: ProfXkit


text:
classdef MrGlm < MrCopyData
    % Class providing General Linear Model of fMRI data
    % (for mass-univariate analysis, i.e. per-voxel regression)
    %
    %
    % EXAMPLE
    %   MrGlm
    %
    %   See also

    % Author:   Saskia Klein & Lars Kasper
    % Created:  2014-07-08
    % Copyright (C) 2014 Institute for Biomedical Engineering
    %                    University of Zurich and ETH Zurich
    %
    % This file is part of the TAPAS UniQC 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/>.
    %

    properties

        % Multiple regressors, i.e. confounds that will be regressed directly,
        % without convolution with the hemodynamic response function(s).
        % fields:
        %      realign     Realignment parameters
        regressors = struct( ...
            'realign', [], ...
            'physio', [], ...
            'other', [] ...
            );

        % Multiple conditions, i.e. behavioral/neuronal regressors that will be
        % convolved with the hemodynamic response function(s).
        conditions = struct( ...
            'names', [], ...
            'onsets', [], ...
            'durations', [] ...
            );

        % The final design matrix used for the voxel-wise regression
        designMatrix = [];

        % timing parameters
        timingUnits = '';
        repetitionTime = '';

        % HRF derivatives
        hrfDerivatives = '';

        % SPM Directory
        parameters = struct( ...
            'save', struct( ...
            'path', '', ... % path where the SPM file is stored, e.g. the MrSeries path
            'spmDirectory', '')... % name of the SPM Directory for the SPM file
            );

        % masking threshold, defined as proportion of globals
        maskingThreshold = 0.8;

        % explicit mask for the analysis, e.g. the segmentation results for a
        % whithin brain mask
        explicitMasking = '';

        % serial correlations, AR(1) or FAST
        serialCorrelations = 'AR(1)';

        % estimation method (classical or Bayesian)
        estimationMethod = 'classical'

        % AR model order (for Bayesian estimation only)
        ARModelOrderBayes = 3;

        % contrasts need to be specified before hand as well
        gcon = struct(...
            'name', {'main_positive', 'main negative'}, ...
            'convec', {1, -1});

    end % properties

    methods

        % Constructor of class
        function this = MrGlm()
        end

        % NOTE: Most of the methods are saved in separate function.m-files in this folder;
        %       except: constructor, delete, set/get methods for properties.

    end % methods

end

author: translationalneuromodeling

repo: tapas


text:
function h=vline(varargin)
%VLINE draws vertical lines
%   VLINE(X) draws one vertical line for each element in vector X
%   VLINE(AX,X) draws the lines to the axes specified in AX
%   VLINE(X,...) accepts HG param/value pairs for line objects
%   H=VLINE(...) returns a handle to each line created

%   Copyright 2007-2008 The MathWorks, Inc.
%   $Revision: 1.1.6.3 $  $Date: 2008/05/23 15:39:00 $

if isempty(varargin) || (ishandle(varargin{1}) && length(varargin)==1)
    error('stats:hline:NotEnoughArgs','Not enough arguments');
end

if ishandle(varargin{1})
    ax=varargin{1};
    varargin=varargin(2:end);
else
    ax=gca;
end

x = varargin{1};
varargin=varargin(2:end);

hh=graph2d.constantline(x,'parent',ax,varargin{:});
hh = changedependvar(hh,'x');

if nargout>0
    h=hh;
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function hfssThickenSheet(fid,ObjectList,Thickness, Units)
% ----------------------------------------------------------------------------
% function hfssThickenSheet(fid, ObjectList, Thickness, Units)
%
% Description :
% -------------
% Thickens a given set of sheets.
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% ObjectList - a cell array of objects that need get thickened.
% Thickness - a float that defines the thickness value
% Units   - units for the thickness (use 'in', 'mm', 'meter' or
%           anything else defined in HFSS).
%
% Note :
% ------
%
% Example :
% ---------
% hfssThickenSheet(fid, {'Rectangle1', 'Rectangle2'}, 0.4, 'mm');
% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
% Written by Franz Schroeder
% franz.dude@gmail.com / franz.schroeder@tu-dresden.de
% 16 September 2013
% -------------------------------------------------------------------------
% ---
% ----------------------------------------------------------------------------

nObjects = length(ObjectList);

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'oEditor.ThickenSheet _\n');
fprintf(fid, 'Array("NAME:Selections", _\n');

% Object Selections.
fprintf(fid, '"Selections:=", "');
for iObj = 1:nObjects
	fprintf(fid, '%s', ObjectList{iObj});
	if (iObj ~= nObjects)
		fprintf(fid, ',');
	end
end
fprintf(fid, '", _\n');

%set parameters
fprintf(fid, '"NewPartsModelFlag:=", _\n');
fprintf(fid, '"Model"), _\n');
fprintf(fid, 'Array("NAME:SheetThickenParameters", _\n');
fprintf(fid, '"Thickness:=", "%f%s", _\n', Thickness, Units);
fprintf(fid, '"BothSides:=", false)\n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function chosenROI = chooseROIwithMouseMontage(view)
%
% chosenROI = chooseROIwithMouseMontage(view)
%
% Use mouse click to pick an ROI.  The mouse click must be within
% a small fraction of a pixel of one of the ROIs.
%
% djh, 1/10/98
% ras, updated for montage view, 09/04

% First put up some instructions
disp('Click on desired ROI to choose it.');

[y,x,b] = ginput(1);
newCoords = montage2Coords(view, [y x]');
x = newCoords(1);
y = newCoords(2);
z = newCoords(3);

chosenROI=0;
for r=1:length(view.ROIs)
    % Check if clicked on one of the coords.
    coords = view.ROIs(r).coords;
    % Get min distance between selected pixel and ROI
    if ~isempty(coords)
        xdiff = coords(1,:)-x;
        ydiff = coords(2,:)-y;
        zdiff = coords(3,:)-z;
        diff = xdiff.^2 + ydiff.^2 + zdiff.^2;
        if find (diff < 1e-2)
            chosenROI=r;
        end
    end
end

if chosenROI
    disp(['You have chosen: ',view.ROIs(chosenROI).name]);
else
    myErrorDlg('You must click ON an ROI to choose it.');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function scale_sample = extract_scale_sample(im, pos, base_target_sz, scaleFactors, scale_model_sz, use_mexResize)

% Get scale filter sample.

if nargin < 6
    use_mexResize = true;
end

nScales = length(scaleFactors);

% Downsample factor.
df = floor(min(scaleFactors));
if df > 1
    im = im(1:df:end,1:df:end,:);
    pos = (pos - 1) / df + 1;
    scaleFactors = scaleFactors / df;
end

for s = 1:nScales
    patch_sz = floor(base_target_sz * scaleFactors(s));

    xs = floor(pos(2)) + (1:patch_sz(2)) - floor(patch_sz(2)/2);
    ys = floor(pos(1)) + (1:patch_sz(1)) - floor(patch_sz(1)/2);

    %check for out-of-bounds coordinates, and set them to the values at
    %the borders
    xs(xs < 1) = 1;
    ys(ys < 1) = 1;
    xs(xs > size(im,2)) = size(im,2);
    ys(ys > size(im,1)) = size(im,1);

    %extract image
    im_patch = im(ys, xs, :);

    % resize image to model size
    if use_mexResize
        im_patch_resized = mexResize(im_patch, scale_model_sz, 'auto');
    else
        im_patch_resized = imresize(im_patch, scale_model_sz, 'bilinear', 'Antialiasing',false);
    end

    % extract scale features
    temp_hog = fhog(single(im_patch_resized), 4);

    if s == 1
        dim_scale = size(temp_hog,1)*size(temp_hog,2)*31;
        scale_sample = zeros(dim_scale, nScales, 'single');
    end

    scale_sample(:,s) = reshape(temp_hog(:,:,1:31), dim_scale, 1);
end


text:
function pass = test_fevalm( pref )
% Test spherefun/fevalm

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

tol = 100*pref.cheb2Prefs.chebfun2eps;
rng(2016);

% Check empty spherefun:
f = spherefun;
s = pi*(2*rand(5,1) - 1);
t = pi/2*rand(5,1);
B = fevalm(f, s, t);
pass(1) = isempty( B );

% Check rank 1 spherefun:
f = chebfun2(@(lam,th) cos(lam).*sin(th));
s = pi*(2*rand(5,1) - 1);
t = pi/2*rand(5,1);
[ss, tt] = meshgrid( s, t);
A = feval(f, ss, tt);
B = fevalm(f, s, t);
pass(2) = norm( A - B ) < tol;

% Check essentially one dimensional function:
f = spherefun(@(lam,th) exp(-(cos(th)-1).^2));
s = pi*(2*rand(5,1) - 1);
t = pi/2*rand(5,1);
[ss, tt] = meshgrid( s, t);
A = feval(f, ss, tt);
B = fevalm(f, s, t);
pass(3) = norm( A - B ) < tol;

end


text:
function [ t_stop, y_stop ] = p28_stop ( neqn )

%*****************************************************************************80
%
%% P28_STOP returns the stopping point for problem p28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
   -3.726957553088175e-01; ...
   -6.230137949234190e-01 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function M = multinomialdoublystochasticfactory(n)
% Manifold of n-by-n doubly-stochastic matrices with positive entries.
%
% function M = multinomialdoublystochasticfactory(n)
%
% M is a Manopt manifold structure to optimize over the set of n-by-n
% matrices with (strictly) positive entries and such that the entries of
% each column and each row sum to one.
%
% Points on the manifold and tangent vectors are represented naturally as
% matrices of size n. The Riemannian metric imposed on the manifold is the
% Fisher metric, that is, if X is a point on the manifold and U, V are two
% tangent vectors:
%
%     M.inner(X, U, V) = <U, V>_X = sum(sum(U.*V./X)).
%
% The retraction here provided is only first order. Consequently, the
% slope test in the checkhessian tool is only valid at points X where the
% gradient is zero. Furthermore, if some entries of X are very close to
% zero, this may cause numerical difficulties that can also lead to a
% failed slope test. More generally, it is important that the solution of
% the optimization problem should have positive entries, sufficiently far
% away from zero to avoid numerical issues.
%
% The file is based on developments in the research paper
% A. Douik and B. Hassibi, "Manifold Optimization Over the Set
% of Doubly Stochastic Matrices: A Second-Order Geometry"
% ArXiv:1802.02628, 2018.
%
% Link to the paper: https://arxiv.org/abs/1802.02628.
%
% Please cite the Manopt paper as well as the research paper:
% @Techreport{Douik2018Manifold,
%   Title   = {Manifold Optimization Over the Set of Doubly Stochastic
%              Matrices: {A} Second-Order Geometry},
%   Author  = {Douik, A. and Hassibi, B.},
%   Journal = {Arxiv preprint ArXiv:1802.02628},
%   Year    = {2018}
% }
%
% See also: multinomialsymmetricfactory multinomialfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Ahmed Douik, March 6, 2018.
% Contributors: Nicolas Boumal
% Change log:
%
%    Apr. 24, 2018 (AD):
%        Changed pinv() to a particular solution to the equation.
%
%    July 24, 2018 (AD):
%        A bugfix related to the pinv() change, with effects in many places.
%
%    Sep.  6, 2018 (NB):
%        Removed M.exp() as it was not implemented.
%
%    Aug. 19, 2019 (AD, BM, NB):
%        Fixed typos in comments; replaced some linear solves with pcg
%        for efficiency when n is large. By default, pcg is used:
%        comments in the code indicate other possibilities and how they
%        differ. Added maxDSiters to doubly_stochastic argument.
%        The main change has been to factor out these linear solves.

    e = ones(n, 1);

    maxDSiters = 100 + 2*n;

    function [alpha, beta] = mylinearsolve(X, b)
        % zeta = sparse(A)\b; % sparse might not be better perf.-wise.
        % where A = [eye(n) X ; X' eye(n)];
        %
        % For large n use the pcg solver instead of \.
        % [zeta, ~, ~, iter] = pcg(sparse(A), b, 1e-6, 100);
        %
        % Even faster is to create a function handle
        % computing A*x (x is a given vector).
        % Make sure that A is not created, and X is only
        % passed with mylinearsolve and not A.
        [zeta, ~, ~, iter] = pcg(@mycompute, b, 1e-6, 100);

        function Ax = mycompute(x)
            xtop = x(1:n,1);
            xbottom = x(n+1:end,1);
            Axtop = xtop + X*xbottom;
            Axbottom = X'*xtop + xbottom;
            Ax = [Axtop; Axbottom];
        end

        alpha = zeta(1:n, 1);
        beta = zeta(n+1:end, 1);
    end

    M.name = @() sprintf('%dx%d doubly-stochastic matrices with positive entries', n, n);

    M.dim = @() (n-1)^2;

    % Fisher metric
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        ip = sum((eta(:).*zeta(:))./X(:));
    end

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(X, Y) error('multinomialdoublystochasticfactory.dist not implemented yet.');

    % The manifold is not compact as a result of the choice of the metric,
    % thus any choice here is arbitrary. This is notably used to pick
    % default values of initial and maximal trust-region radius in the
    % trustregions solver.
    M.typicaldist = @() n;

    % Pick a random point on the manifold
    M.rand = @random;
    function X = random()
        Z = abs(randn(n, n));     % Random point in the ambient space
        X = doubly_stochastic(Z, maxDSiters); % Projection on the Manifold
    end

    % Pick a random vector in the tangent space at X.
    M.randvec = @randomvec;
    function eta = randomvec(X) % A random vector in the tangent space
        % A random vector in the ambient space
        Z = randn(n, n);
        % Projection of the vector onto the tangent space
        b = [sum(Z, 2) ; sum(Z, 1)'];
        [alpha, beta] = mylinearsolve(X, b);
        eta = Z - (alpha*e' + e*beta').*X;
        % Normalizing the vector
        nrm = M.norm(X, eta);
        eta = eta / nrm;
    end

    % Projection of vector eta in the ambient space to the tangent space.
    M.proj = @projection;
    function etaproj = projection(X, eta) % Projection of the vector eta in the ambeint space onto the tangent space
        b = [sum(eta, 2) ; sum(eta, 1)'];
        [alpha, beta] = mylinearsolve(X, b);
        etaproj = eta - (alpha*e' + e*beta').*X;
    end

    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;

    % Conversion of Euclidean to Riemannian gradient
    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad) % projection of the euclidean gradient
        mu = (X.*egrad);
        b = [sum(mu, 2) ; sum(mu, 1)'];
        [alpha, beta] = mylinearsolve(X, b);
        rgrad = mu - (alpha*e' + e*beta').*X;
    end

    % First-order retraction
    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y = X.*exp(t*(eta./X));
        Y = doubly_stochastic(Y, maxDSiters);
        Y = max(Y, eps);
    end

    % Conversion of Euclidean to Riemannian Hessian
    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, eta)

        % computing the directional derivative of the Riemannian
        % gradient
        gamma = egrad.*X;
        gammadot = ehess.*X + egrad.*eta;

        b = [sum(gamma, 2) ; sum(gamma, 1)'];
        bdot = [sum(gammadot, 2) ; sum(gammadot, 1)'];
        [alpha, beta] = mylinearsolve(X, b);
        [alphadot, betadot] = mylinearsolve(X, bdot- [eta*beta; eta'*alpha]);

        S = (alpha*e' + e*beta');
        deltadot = gammadot - (alphadot*e' + e*betadot').*X- S.*eta;

        % projecting gamma
        delta = gamma - S.*X;

        % computing and projecting nabla
        nabla = deltadot - 0.5*(delta.*eta)./X;
        rhess = projection(X, nabla);
    end

    % Miscellaneous manifold functions
    M.hash = @(X) ['z' hashmd5(X(:))];
    M.lincomb = @matrixlincomb;
    M.zerovec = @(X) zeros(n, n);
    M.transp = @(X1, X2, d) projection(X2, d);
    M.vec = @(X, U) U(:);
    M.mat = @(X, u) reshape(u, n, n);
    M.vecmatareisometries = @() false;

end

author: NicolasBoumal

repo: manopt


text:
function [net,stats] = train_rand(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% 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).
addpath(fullfile(vl_rootnn, 'examples'));

%%%%%%%%%%%%% new fields %%%%%%%%%%%%%
opts.saveInterval = 2;
%%%%%%%%%%%%% new fields %%%%%%%%%%%%%

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.epochSize = inf;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = [] ;  % Empty array means use the default SGD solver
[opts, varargin] = vl_argparse(opts, varargin) ;
if ~isempty(opts.solver)
    assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
        'Invalid solver; expected a function handle with two outputs.') ;
    % Call without input arguments, to get default options
    opts.solverOpts = opts.solver() ;
end

opts.momentum = 0.9 ;
opts.saveSolverState = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isscalar(opts.train) && isnumeric(opts.train) && isnan(opts.train)
    opts.train = [] ;
end
if isscalar(opts.val) && isnumeric(opts.val) && isnan(opts.val)
    opts.val = [] ;
end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
    if isempty(opts.derOutputs)
        error('DEROUTPUTS must be specified when training.\n') ;
    end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
    fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
    % [KH] make sure to use the input opts, saved maybe outdated
    opt = net.layers(end).block.opt;
    [net, state, stats] = loadState(modelPath(start)) ;
    net.layers(end).block.opt = opt;
else
    state = [] ;
end

for epoch=start+1:opts.numEpochs

    % Set the random seed based on the epoch and opts.randomSeed.
    % This is important for reproducibility, including when training
    % is restarted from a checkpoint.

    rng(epoch + opts.randomSeed) ;
    prepareGPUs(opts, epoch == start+1) ;

    % Train for one epoch.
    params = opts ;
    params.epoch = epoch ;
    params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
    params.train = opts.train; %(randperm(numel(opts.train))) ; % shuffle
    params.train = params.train; %(1:min(opts.epochSize, numel(opts.train)));
    params.val = opts.val; %(randperm(numel(opts.val))) ;
    params.imdb = imdb ;
    params.getBatch = getBatch ;

    if numel(opts.gpus) <= 1
        [net, state] = processEpoch(net, state, params, 'train') ;
        [net, state] = processEpoch(net, state, params, 'val') ;
        if ~evaluateMode && ~mod(epoch, opts.saveInterval)
            saveState(modelPath(epoch), net, state) ;
        end
        lastStats = state.stats ;
    else
        spmd
            [net, state] = processEpoch(net, state, params, 'train') ;
            [net, state] = processEpoch(net, state, params, 'val') ;
            if labindex == 1 && ~evaluateMode && ~mod(epoch, opts.saveInterval)
                saveState(modelPath(epoch), net, state) ;
            end
            lastStats = state.stats ;
        end
        lastStats = accumulateStats(lastStats) ;
    end

    stats.train(epoch) = lastStats.train ;
    stats.val(epoch) = lastStats.val ;
    clear lastStats ;
    saveStats(modelPath(epoch), stats) ;

    if opts.plotStatistics
        switchFigure(1) ; clf ;
        plots = setdiff(...
            cat(2,...
            fieldnames(stats.train)', ...
            fieldnames(stats.val)'), {'num', 'time'}) ;
        for p = plots
            p = char(p) ;
            values = zeros(0, epoch) ;
            leg = {} ;
            for f = {'train', 'val'}
                f = char(f) ;
                if isfield(stats.(f), p)
                    tmp = [stats.(f).(p)] ;
                    values(end+1,:) = tmp(1,:)' ;
                    leg{end+1} = f ;
                end
            end
            subplot(1,numel(plots),find(strcmp(p,plots))) ;
            plot(1:epoch, values','o-') ;
            xlabel('epoch') ;
            title(p) ;
            legend(leg{:}) ;
            grid on ;
        end
        drawnow ;
        print(1, modelFigPath, '-dpdf') ;
    end

    if ~isempty(opts.postEpochFn)
        if nargout(opts.postEpochFn) == 0
            opts.postEpochFn(net, params, state) ;
        else
            [lr, params] = opts.postEpochFn(net, params, state) ;
            if ~isempty(lr), opts.learningRate = lr; end
            if opts.learningRate == 0, break; end
            opts.train = params.train;
        end
    end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
    state.solverState = cell(1, numel(net.params)) ;
    state.solverState(:) = {0} ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
    net.move('gpu') ;
    for i = 1:numel(state.solverState)
        s = state.solverState{i} ;
        if isnumeric(s)
            state.solverState{i} = gpuArray(s) ;
        elseif isstruct(s)
            state.solverState{i} = structfun(@gpuArray, s, 'UniformOutput', false) ;
        end
    end
end
if numGpus > 1
    parserv = ParameterServer(params.parameterServer) ;
    net.setParameterServer(parserv) ;
else
    parserv = [] ;
end

% profile
if params.profile
    if numGpus <= 1
        profile clear ;
        profile on ;
    else
        mpiprofile reset ;
        mpiprofile on ;
    end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

if numel(subset) > 0
logInfo('Epoch %d/%d: %s LR=%d', epoch, params.numEpochs, ...
    func2str(params.solver), params.learningRate);
start = tic;  a = tic;
objs = [];
for t=1:params.batchSize:numel(subset)
    batchSize = min(params.batchSize, numel(subset) - t + 1) ;

    for s=1:params.numSubBatches
        % get this image batch and prefetch the next
        batchStart = t + (labindex-1) + (s-1) * numlabs ;
        batchEnd = min(t+params.batchSize-1, numel(subset)) ;
        batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
        num = num + numel(batch) ;
        if numel(batch) == 0, continue ; end

        inputs = params.getBatch(params.imdb, batch) ;

        if params.prefetch
            if s == params.numSubBatches
                batchStart = t + (labindex-1) + params.batchSize ;
                batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
            else
                batchStart = batchStart + numlabs ;
            end
            nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
            params.getBatch(params.imdb, nextBatch) ;
        end

        if strcmp(mode, 'train')
            net.mode = 'normal' ;
            net.conserveMemory = true;
            net.accumulateParamDers = (s ~= 1) ;
            net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
        else
            net.mode = 'test' ;
            net.eval(inputs) ;
        end
    end

    % Accumulate gradient.
    if strcmp(mode, 'train')
        if ~isempty(parserv), parserv.sync() ; end
        state = accumulateGradients(net, state, params, batchSize, parserv) ;
    end

    % Get statistics.
    time = toc(start) + adjustTime ;
    batchTime = time - stats.time ;
    stats.num = num ;
    stats.time = time ;
    stats = params.extractStatsFn(stats,net) ;
    currentSpeed = batchSize / batchTime ;
    averageSpeed = (t + batchSize - 1) / time ;
    if t == 3*params.batchSize + 1
        % compensate for the first three iterations, which are outliers
        adjustTime = 4*batchTime - time ;
        stats.time = time + adjustTime ;
    end

    if toc(a) > 30
        fprintf('%s-%s ep%02d: %3d/%3d (%d):', ...
            net.layers(end).block.opt.dataset, mode, epoch, ...
            fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize), ...
            params.batchSize) ;
        %fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
        fprintf(' %.1fHz', averageSpeed) ;
        for f = setdiff(fieldnames(stats)', {'num', 'time'})
            f = char(f) ;
            fprintf(' %s: %.3f', f, stats.(f)) ;
        end
        fprintf('\n') ;
        a = tic;
    end
    objs = [objs, stats.objective];
end
logInfo('Epoch %d, Avg %s obj = %g', epoch, mode, mean(objs(~isnan(objs))));
end  % if N>0

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
    if numGpus <= 1
        state.prof.(mode) = profile('info') ;
        profile off ;
    else
        state.prof.(mode) = mpiprofile('info');
        mpiprofile off ;
    end
end
if ~params.saveSolverState
    state.solverState = [] ;
else
    for i = 1:numel(state.solverState)
        s = state.solverState{i} ;
        if isnumeric(s)
            state.solverState{i} = gather(s) ;
        elseif isstruct(s)
            state.solverState{i} = structfun(@gather, s, 'UniformOutput', false) ;
        end
    end
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

    if ~isempty(parserv)
        parDer = parserv.pullWithIndex(p) ;
    else
        parDer = net.params(p).der ;
    end

    switch net.params(p).trainMethod
    case 'average' % mainly for batch normalization
        thisLR = net.params(p).learningRate ;
        net.params(p).value = vl_taccum(...
            1 - thisLR, net.params(p).value, ...
            (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
        thisDecay = params.weightDecay * net.params(p).weightDecay ;
        thisLR = params.learningRate * net.params(p).learningRate ;

        if thisLR>0 || thisDecay>0
            % Normalize gradient and incorporate weight decay.
            try
            parDer = vl_taccum(1/batchSize, parDer, ...
                thisDecay, net.params(p).value) ;
        catch, net.params(p), keyboard, end

            if isempty(params.solver)
                % Default solver is the optimised SGD.
                % Update momentum.
                state.solverState{p} = vl_taccum(...
                    params.momentum, state.solverState{p}, ...
                    -1, parDer) ;

                % Nesterov update (aka one step ahead).
                if params.nesterovUpdate
                    delta = params.momentum * state.solverState{p} - parDer ;
                else
                    delta = state.solverState{p} ;
                end

                % Update parameters.
                net.params(p).value = vl_taccum(...
                    1,  net.params(p).value, thisLR, delta) ;

            else
                % call solver function to update weights
                [net.params(p).value, state.solverState{p}] = ...
                    params.solver(net.params(p).value, state.solverState{p}, ...
                    parDer, params.solverOpts, thisLR) ;
            end
        end
    otherwise
        error('Unknown training method ''%s'' for parameter ''%s''.', ...
            net.params(p).trainMethod, ...
            net.params(p).name) ;
    end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------
for s = {'train', 'val'}
    s = char(s) ;
    total = 0 ;

    % initialize stats stucture with same fields and same order as
    % stats_{1}
    stats__ = stats_{1} ;
    names = fieldnames(stats__.(s))' ;
    values = zeros(1, numel(names)) ;
    fields = cat(1, names, num2cell(values)) ;
    stats.(s) = struct(fields{:}) ;

    for g = 1:numel(stats_)
        stats__ = stats_{g} ;
        num__ = stats__.(s).num ;
        total = total + num__ ;

        for f = setdiff(fieldnames(stats__.(s))', 'num')
            f = char(f) ;
            stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

            if g == numel(stats_)
                stats.(s).(f) = stats.(s).(f) / total ;
            end
        end
    end
    stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
    if net.layers(sel(i)).block.ignoreAverage, continue; end;
    stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;
logInfo('saved: %s', fileName);

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
    try, save(fileName, 'stats', '-append') ; end
else
    try, save(fileName, 'stats') ; end
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
    warning('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
    try
        set(0,'CurrentFigure',n) ;
    catch
        figure(n) ;
    end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
    % check parallel pool integrity as it could have timed out
    pool = gcp('nocreate') ;
    if ~isempty(pool) && pool.NumWorkers ~= numGpus
        delete(pool) ;
    end
    pool = gcp('nocreate') ;
    if isempty(pool)
        parpool('local', numGpus) ;
        cold = true ;
    end

end
if numGpus >= 1 && cold
    fprintf('%s: resetting GPU\n', mfilename)
    clearMex() ;
    if numGpus == 1
        gpuDevice(opts.gpus)
    else
        spmd
            clearMex() ;
            gpuDevice(opts.gpus(labindex))
        end
    end
end

author: kunhe

repo: FastAP-metric-learning


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

%*****************************************************************************80
%
%% HALF_NORMAL_CDF evaluates the Half Normal CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= a )
    cdf = 0.0;
  else
    cdf2 = normal_cdf ( x, a, b );
    cdf = 2.0 * cdf2 - 1.0;
  end

  return
end


text:
classdef nme_gizmo_acp < mp.nme_gizmo_ac & mp.form_acp

%   MATPOWER
%   Copyright (c) 2019-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function y = negLogLogitTransform(x, transform)

% NEGLOGLOGITTRANSFORM Constrains a parameter to be positive.
% FORMAT
% DESC contains commands to constrain parameters to be positive via
% log(1+exp(x)).
% ARG x : input argument.
% ARG y : return argument.
% ARG transform : type of transform, 'atox' maps a value into
% the transformed space (i.e. makes it positive). 'xtoa' maps the
% parameter back from transformed space to the original
% space. 'gradfact' gives the factor needed to correct gradients
% with respect to the transformed parameter.
%
% SEEALSO : expTransform, sigmoidTransform
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2007

% OPTIMI

y = zeros(size(x));
limVal = 36;
switch transform
 case 'atox'
  index = find(x<-limVal);
  y(index) = eps;
  x(index) = NaN;
  index = find(x<limVal);
  y(index) = log(1+exp(x(index)));
  x(index) = NaN;
  index = find(~isnan(x));
  y(index) = x(index);
 case 'xtoa'
  index = find(x<limVal);
  y(index) = log(exp(x(index))-1);
  index = find(x>=limVal);
  y(index) = x(index);
 case 'gradfact'
  index = find(x>limVal);
  y(index) = 1;
  index = find(x<=limVal);
  y(index) = 1-exp(-x(index));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = ellipsoidMesh(elli, varargin)
%ELLIPSOIDMESH Convert a 3D ellipsoid to face-vertex mesh representation.
%
%   [V, F] = ellipsoidMesh(ELLI)
%   ELLI is given by:
%   [XC YC ZC  A B C  PHI THETA PSI],
%   where (XC, YC, ZC) is the ellipsoid center, A, B and C are the half
%   lengths of the ellipsoid main axes, and PHI THETA PSI are Euler angles
%   representing ellipsoid orientation, in degrees.
%
%
%   See also
%   meshes3d, drawEllipsoid, sphereMesh, inertiaEllipsoid
%

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

%% Default values

% number of meridians
nPhi    = 32;

% number of parallels
nTheta  = 16;

%% Extract input arguments

% Parse the input (try to extract center coordinates and radius)
if nargin == 0
    % no input: assumes ellipsoid with default shape
    elli = [0 0 0 5 4 3 0 0 0];
end

% default set of options for drawing meshes
options = {'FaceColor', 'g', 'linestyle', 'none'};

while length(varargin) > 1
    switch lower(varargin{1})
        case 'nphi'
            nPhi = varargin{2};

        case 'ntheta'
            nTheta = varargin{2};

        otherwise
            % assumes this is drawing option
            options = [options varargin(1:2)]; %#ok<AGROW>
    end

    varargin(1:2) = [];
end

%% Parse numerical inputs

% Extract ellipsoid parameters
xc  = elli(:,1);
yc  = elli(:,2);
zc  = elli(:,3);
a   = elli(:,4);
b   = elli(:,5);
c   = elli(:,6);
k   = pi / 180;
ellPhi   = elli(:,7) * k;
ellTheta = elli(:,8) * k;
ellPsi   = elli(:,9) * k;

%% Coordinates computation

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

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

%% parametrisation of ellipsoid

% spherical coordinates
theta   = linspace(0, pi, nTheta+1);
phi     = linspace(0, 2*pi, nPhi+1);

% convert to cartesian coordinates
sintheta = sin(theta);
x = cos(phi') * sintheta;
y = sin(phi') * sintheta;
z = ones(length(phi),1) * cos(theta);

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

% convert to FV mesh
[vertices, faces] = surfToMesh(x, y, z, 'xPeriodic', false, 'yPeriodic', true);

% format output
varargout = formatMeshOutput(nargout, vertices, faces);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%%WTF??????? I wrote this and I thought I saved periodically, but it looks
%%%like there's no trace of it after the @#*(%&(*& xwindows crashed.
%
% To restore, copy code from Bob's dtiFiberUI dialog (e.g., xform/translate
% ROIs. Pass in as cells a long list of the parameters users can edit, plus
% default values from ctrCreate. Don't forget to grab and output what the
% user does. Also, fix the Q#(*%&(*#%& ctrInit.fig\
%
% ALSO: we pass in params and pass out params and change a field here as a
% flag to indicate that the ctrInitParamsFile should use the values already
% in PARAMS rather than reset them.


text:
function chi = RoadInclination_Arkustangens(z_I_Ground_m,rx_I_IAji_m,ry_I_IAji_m,phiz_IK_rad)
%% Computation of Road Angles from Street Profile in Inertial Frame
% Input parameters:
% phiz_IK_rad     [rad]    Euler-Angle around z-Axis from Inertial Reference Frame I to Vehicle Fixed Reference Frame K
% rx_I_IAji_m     [m]      Track x-Position [x1 x2 x3 x4] [1x4]
% ry_I_IAji_m     [m]      Track y-Position [y1 y2 y3 y4] [1x4]
% z_I_Ground_m    [m]      Input Vector with z-Profile of Street in Inertial Coordinate System at the Axle Locations [z1 z2 z3 z4] [3x4]
% Output parameters:
% chi_rad         [rad]    Approximated Cardan Angles from Inertial Reference Frame to Road Fixed Reference Frame [chi_x chi_y Psi] [3x1]

%% Computation of Road Angles

% Distances [dxl dxr dyf dyr] between suspension hardpoints
dxl=sqrt( (rx_I_IAji_m(1)-rx_I_IAji_m(3)).^2 + (ry_I_IAji_m(1)-ry_I_IAji_m(3)).^2 );
dxr=sqrt( (rx_I_IAji_m(2)-rx_I_IAji_m(4)).^2 + (ry_I_IAji_m(2)-ry_I_IAji_m(4)).^2 );
dyf=sqrt( (rx_I_IAji_m(1)-rx_I_IAji_m(2)).^2 + (ry_I_IAji_m(1)-ry_I_IAji_m(2)).^2 );
dyr=sqrt( (rx_I_IAji_m(3)-rx_I_IAji_m(4)).^2 + (ry_I_IAji_m(3)-ry_I_IAji_m(4)).^2 );

% Inclination Angle in x-direction for negative z-road-profile
chi_y = -atan( 0.5*((z_I_Ground_m(1)+z_I_Ground_m(2))-(z_I_Ground_m(3)+z_I_Ground_m(4)))/(0.5*(dxl+dxr)) );
% Inclination Angle in y-direction for negative z-road-profile
chi_x = -atan( 0.5*((z_I_Ground_m(1)+z_I_Ground_m(3))-(z_I_Ground_m(2)+z_I_Ground_m(4)))/(0.5*(dyf+dyr)) );

%% Output Parameters
chi = [chi_x;chi_y;phiz_IK_rad];

end

author: TUMFTM


text:
function model = setFoodConstraints(model, foodMenu, compensateDiet)
% This function sets diet/food constraints for uptake reactions of
% a whole-body metabolic model.
%
% function model = setFoodConstraints(model, Diet, factor)
%
% INPUT
%   model         model structure
%   foodMenu      {column vector of food items, column vector of serving mass (g))
%Optional Input:
%   compensateDiet   boolean 0 and 1; 1 indicates to compensate for gaps in nutrition components in the EUAverageDiet. default is 1

% OUTPUT
%   model         updated model
%
% AUTHORS
%   Bronson R. Weston 2021-2022

%Add any food reactions if necessary
foodMenu=makeDietFormatConsistant(model,foodMenu);
[model,~] = addFoodRxns2HM(model,foodMenu(:,1));

%Clear any current dietary/food reaction fluxes
fdRxns=find(contains(model.rxns,'Food_EX_'));
metRxns=find(contains(model.rxns,'Diet_EX_'));
model.lb(fdRxns)=0;
model.ub(fdRxns)=0;
model.lb(metRxns)=0;
model.ub(metRxns)=0;

%Add diet components
for i=1:length(foodMenu(:,1))
    %     foodMenu{i,1}
    %     find(ismember(model.rxns,foodMenu{i,1}))
    if ischar(foodMenu{i,2})
        val= -1*str2num(foodMenu{i,2});
    else
        val= -1*foodMenu{i,2};
    end
    model.ub(find(ismember(model.rxns,foodMenu{i,1}))) = val;
    model.lb(find(ismember(model.rxns,foodMenu{i,1}))) = val;
    if isempty(find(ismember(model.rxns,foodMenu{i,1})))
        error(['"' foodMenu{i,1} '" is not a valid food item/ metabolite'])
    end
end

if ~exist('compensateDiet','var')
    compensateDiet=1;
end

% load 'fdTable.mat'
% load 'fdCategoriesTable.mat'
metFlux = getMetaboliteFlux(foodMenu);

if compensateDiet~=1
    model.SetupInfo.DietComposition = metFlux;
    return
end

% If compensating diet, make changes in accordance with setDietConstraints
% and EUAverageDietNew

EUAverageDietNew
Diet = getMetaboliteFlux(Diet);

for i=1:length(Diet(:,1))
    [missingEUmets,ind]=setdiff(Diet(:,1),metFlux(:,1));
    metFlux=[metFlux; Diet(ind,:)];
end

microbiotaEnabling =1; % consider essential microbial metabolites to be added to diet

% load essential metabolite list for AGORA models
AGORAEssentialMetabolites;
AGORAessential = regexprep(AGORAessential,'EX_','Diet_EX_');
AGORAessential = regexprep(AGORAessential,'\[u\]','\[d\]');

% ensure uptake of metabolites required for microbiota (may not all be
% needed)
if microbiotaEnabling == 1
    MissingUptakes = setdiff(AGORAessential,metFlux(:,1));
    % open uptake for those reactions
    % set constraints to a default
    model.lb(ismember(model.rxns,MissingUptakes))=-0.1;
end

% these compounds are in the diet and needed for the proper function of the
% model but are not reported in our diet database
if 1
    PotentialMissingDietCompounds={'Diet_EX_asn_L[d]'; 'Diet_EX_gln_L[d]';'Diet_EX_chol[d]';'Diet_EX_crn[d]';'Diet_EX_elaid[d]';...
        'Diet_EX_hdcea[d]';'Diet_EX_dlnlcg[d]';'Diet_EX_adrn[d]';'Diet_EX_hco3[d]';...
        'Diet_EX_sprm[d]'; 'Diet_EX_carn[d]';'Diet_EX_7thf[d]';...
        'Diet_EX_Lcystin[d]';%??
        'Diet_EX_hista[d]';'Diet_EX_orn[d]';...
        'Diet_EX_ptrc[d]';'Diet_EX_creat[d]';
        'Diet_EX_cytd[d]'; 'Diet_EX_so4[d]'
        };
    MissingDietCompounds=setdiff(PotentialMissingDietCompounds,metFlux(:,1));
    model.lb(ismember(model.rxns,MissingDietCompounds))=-50;
end

% micronutrient - defined to have mole/day/person rate below 1e-6 mol/day/person
% lower bounds will be relaxed by factor 10
micronutrients ={%'Diet_EX_adpcbl[d]'
    % I changed the exchange ID since it was generated based on the
    % metabolite (ID adocbl). AH 16/12/01
    'Diet_EX_adocbl[d]'
    'Diet_EX_vitd2[d]'
    'Diet_EX_vitd3[d]'
    'Diet_EX_psyl[d]'
    'Diet_EX_gum[d]'
    'Diet_EX_bglc[d]'
    'Diet_EX_phyQ[d]'
    'Diet_EX_fol[d]'
    'Diet_EX_5mthf[d]'
    'Diet_EX_q10[d]'
    'Diet_EX_retinol_9_cis[d]'
    'Diet_EX_pydxn[d]'
    'Diet_EX_pydam[d]'
    'Diet_EX_pydx[d]'
    'Diet_EX_pheme[d]'
    'Diet_EX_ribflv[d]'
    'Diet_EX_thm[d]'
    % added 24.08.2016
    'Diet_EX_avite1[d]'
    'Diet_EX_pnto_R[d]'
    };
%no uptake enforced but max uptake rate
ions={    'Diet_EX_na1[d]'	%'0.056546816'
    'Diet_EX_cl[d]'	%'0.056412715'
    'Diet_EX_k[d]'	%'0.12020983'
    'Diet_EX_pi[d]'	%'0.007143207'
    'Diet_EX_zn2[d]'
    'Diet_EX_cu2[d]'
    };
so4={
    %02.06.2017
    'Diet_EX_so4[d]' % see text above about addition of met for brain demand
    };

for i = 1:  size(metFlux(:,1),1)
    R = metFlux{i,1};
    index=find(ismember(model.rxns,R));
    % exception for micronutrients to avoid numberical issues
    if ~isempty(find(ismember(micronutrients,R)))&& metFlux{i,2}<=10
        model.lb(index) = -10;
        model.ub(index) = -10;
%     elseif ~isempty(find(ismember(micronutrients,R))) &&  metFlux{i,2}<0.1
%         model.lb(index) = -1*metFlux{i,2}*99 +model.lb(index);
%         model.ub(index) = model.lb(index);
    elseif ~isempty(find(ismember(ions,R)))
        model.lb(index) = -1* metFlux{i,2}*99+model.lb(index);
        model.ub(index) = model.lb(index);
    elseif ~isempty(find(ismember(so4,R)))
        model.lb(index) = -1000;
        model.ub(index) = -1000;
    elseif ~isempty(find(ismember(missingEUmets,R)))
        model.lb(index) = -1* metFlux{i,2}+model.lb(index);
        model.ub(index) = model.lb(index);
    end
end

model.SetupInfo.DietComposition = metFlux;

author: opencobra

repo: cobratoolbox


text:
function L=GetLevel(i,Ms,n)
%return the level of i-th coefficient
col = floor(i/n);
row = rem(i,n);

dlevels = length(Ms);

p=n/(2^(dlevels-1));

for j=1:dlevels
    if (col<=p) && (row<=p)
        L=j;
        break;
    else
        p=2*p;
    end
end

author: thomaskuestner


text:
function test_pull1663

% MEM 14gb
% WALLTIME 00:60:00
% DEPENDENCY ft_prepare_sourcemodel headsurface ft_prepare_leadfield ft_freqanalysis ft_sourceanalysis

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

% This function creates a set of input-structures to be used for testing
% the DUNEuro MEG forward solution.
% The structure of this script is more or less
% 1. create input data (dull segmentation, hex and tet meshes, dipoles, sensors)
% 2. compute and compare MEG leadfield for hex and tet meshes
%   a. create volume conductor (ft_prepare_headmodel)
%   b. create source grid (ft_prepare_sourcemodel)
%   c. compute leadfield (ft_prepare_leadfield)
%   d. compare hex leadfield with tet leadfield

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. create input data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% create mesh

% from test_pull1427.m
% dull segmentation
segprob = [];
segprob.brain = false(10,10,10); segprob.brain(4:7,4:7,4:7) = true;
segprob.skull = false(10,10,10); segprob.skull(3:8,3:8,3:8) = true;
segprob.scalp = false(10,10,10); segprob.scalp(2:9,2:9,2:9) = true;
%segprob.air = true(10,10,10);
%segprob.air(2:9,2:9,2:9) = false;
%segprob.brain = false(11,11,11); segprob.brain(4:8,4:8,4:8) = true;
%segprob.skull = false(11,11,11); segprob.skull(3:9,3:9,3:9) = true;
%segprob.scalp = false(11,11,11); segprob.scalp(2:10,2:10,2:10) = true;
segprob.dim = size(segprob.brain);
segprob.unit = 'cm';
segprob.coordsys = 'ctf';
segprob.transform = eye(4);
segprob.transform(1,4) = -0.5;
segprob.transform(2,4) = -0.5;
segprob.transform(3,4) = -0.5;

% visualize the segmentation
% it is more difficult to visualize a probabilistic segmentation than an indexed one
segindx = ft_datatype_segmentation(segprob, 'segmentationstyle', 'indexed');

cfg = [];
cfg.funparameter = 'tissue';
cfg.method = 'ortho';
cfg.location = [5 5 5]; % this is the center of the volume, in this plot it will be rounded off to the nearest voxel
ft_sourceplot(cfg, segindx);

% hexa mesh
cfg = [];
cfg.shift = 0.3;
cfg.method = 'hexahedral';
mesh_vol_hex = ft_prepare_mesh(cfg, segprob);

figure
ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
hold on
ft_plot_mesh(mesh_vol_hex, 'surfaceonly', false, 'facecolor', 'none', 'edgecolor', 'm');
view(120, 30)

% tetra mesh
cfg = [];
cfg.method = 'tetrahedral';
mesh_vol_tet = ft_prepare_mesh(cfg, segprob);

figure
ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
hold on
ft_plot_mesh(mesh_vol_tet, 'surfaceonly', false, 'facecolor', 'none', 'edgecolor', 'm');
view(120, 30)

%% define sensors

% i manually passed 5 coils and fixed projections
% or maybe take some from a ctf file? something like this (from test_pull1377.m)

% for MEG data + sensor info
%
% load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf151.mat'), 'data');
% datameg = data;
% clear data

coils = [5 5 12; 5 12 5; -2 5 5; 12 5 5; 5 -2 5];
projections = [0 0 1 ; 0 1 0; -1 0 0; 1 0 0; 0 -1 0 ];

figure
ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
hold on
quiver3(coils(:,1),coils(:,2),coils(:,3),projections(:,1),projections(:,2),projections(:,3),'bo')

meg_labels = cellstr(strings(1,size(coils,1)));
for i=1:size(coils,1)
  meg_labels(i) = cellstr(strcat('meg',num2str(i)));
end

sens = [];
sens.coilpos = coils;
sens.coilori = projections;
sens.chanpos = coils;
sens.chanori = projections;
sens.label = meg_labels;
sens.type = 'meg';
sens.unit = 'mm';
sens.tra = eye(5);
sens = ft_convert_units(sens,'mm');

%% define dipoles
dip_pos = [5.5 6.5 6.5; 5.5 5.5 3.5; 3.5 5.5 5.5];
dip_mom = [0 1 0; 0 0 -1; -1 0 0 ];

figure
ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
hold on
quiver3(dip_pos(:,1),dip_pos(:,2),dip_pos(:,3),dip_mom(:,1),dip_mom(:,2),dip_mom(:,3),'bo')
view(-200, 15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. compute the leadfield
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% prepare headmodel

cfg              = [];
cfg.method       = 'duneuro';
cfg.conductivity = [0.33, 0.43, 0.53];  % vector, conductivity values for tissues: check the order here
vol_duneuro_hex  = ft_prepare_headmodel(cfg, mesh_vol_hex);
vol_duneuro_tet  = ft_prepare_headmodel(cfg, mesh_vol_tet);

%% prepare sourcemodel

% hex
cfg = [];
cfg.sourcemodel.pos = dip_pos';
cfg.sourcemodel.inside = ones(size(dip_pos,2),1);
cfg.grad = sens;
cfg.headmodel = vol_duneuro_hex;
sm_duneuro_hex = ft_prepare_sourcemodel(cfg);

% tet
cfg = [];
cfg.sourcemodel.pos = dip_pos';
cfg.sourcemodel.inside = ones(size(dip_pos,2),1);
cfg.grad = sens;
cfg.headmodel = vol_duneuro_tet;
sm_duneuro_tet = ft_prepare_sourcemodel(cfg);

%% prepare leadfield

% hex
cfg                 = [];
cfg.sourcemodel     = sm_duneuro_hex;
cfg.sourcemodel.mom = dip_mom;
cfg.headmodel       = vol_duneuro_hex;
cfg.grad            = sens;
cfg.reducerank      = 3;
out_hex = ft_prepare_leadfield(cfg);
lf_hex = cell2mat(out_hex.leadfield);

% tet
cfg                 = [];
cfg.sourcemodel     = sm_duneuro_tet;
cfg.sourcemodel.mom = dip_mom;
cfg.headmodel       = vol_duneuro_tet;
cfg.grad            = sens;
cfg.reducerank      = 3;
out_tet = ft_prepare_leadfield(cfg);
lf_tet = cell2mat(out_tet.leadfield);

%% compare leadfields

figure, plot(lf_hex,'r'), hold on, plot(lf_tet,'b')
% set a limit for an error
rel_err_perc = 100*(abs(lf_hex-lf_tet)./norm(lf_hex));
assert(max(max(rel_err_perc))<10)

end % main function

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pass = test_biharm(pref)
% Test BIHARM command.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e6*pref.cheb3Prefs.chebfun3eps;

% Function to be used:
ff = @(x,y,z) x.^2.*y.^2 + x.^2.*z.^2 + y.^2.*z.^2;

% Bihamrmonic operator applied to ff:
f = chebfun3(ff);
fB = biharm(f);

% Exact solution:
gg = @(x,y,z) 24;
g = chebfun3(gg);

% Compare:
pass = norm(fB - g) < tol;

end

author: chebfun


text:
function [ detectedPts, scale ] = dog(  vertex, faces, num_octaves, params )
%vertex, faces
%num_octaves - number of octaves to use (20 works fine)
%params.sigscale - dog scaling method. default is good
%params.ExcludeBoundery - exclude boundary vertices from detection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('params','var')
    params = [];
end

if ~isfield(params,'sigscale')
    params.sigscale = 0;
end
if ~isfield(params,'ExcludeBoundery')
    params.ExcludeBoundery = 1;
end

%dog find geometry difference of gaussians features

    ss = size(vertex);
    if (ss(1) == 3)
        vertex = vertex';
    end

%create ocataves
    %W = my_euclidean_distance(triangulation2adjacency(faces),vertex);
    W = triangulation2adjacency(faces);
    W = W + speye(length(vertex));

    W = spdiags(1./sum(W,2),0,length(vertex),length(vertex))*W;
    vertex2=vertex;
    for ind = 1:num_octaves
        vertex2 = (W*vertex2);
        octaves {ind} = vertex2;
        if (ind > 1)
            diff{ind} =  sum((octaves{ind - 1} - vertex2)'.^2)';
        else
            diff{ind} =  sum((vertex - vertex2)'.^2)';
        end
        ds(ind) = sum(diff{ind});
        if (params.sigscale == 1)
            diff{ind} = diff{ind}*ind;
        elseif (ds(ind) > 0)
            diff{ind} = diff{ind}/ds(ind);
        end
    end

    detectedPts = find_local_maxima(faces,vertex2,diff);

    D = my_euclidean_distance(triangulation2adjacency(faces),vertex);
    d = sum(D);
    w = sum(triangulation2adjacency(faces));
   scale = full((d./w)');
   if (params.ExcludeBoundery)
       boundary = find_boundary_vertex(faces);
       for k = 1:length(detectedPts)
           detectedPts{k} = setdiff(detectedPts{k}, boundary);
       end
   end
  %detectedPts{1} = {};
   %detectedPts{2} = {};
%    a = sum(abs(vertex - repmat([5.92 -63.36 -596.2],length(vertex),1)),2)<0.1;
% b = sum(abs(vertex - repmat([-111.4 -35.49 -580.7],length(vertex),1)),2)<0.1;
% detectedPts = {};
% detectedPts{1}  = [find(a) find(b)];
 end

function W = my_euclidean_distance(A,vertex)

if size(vertex,1)<size(vertex,2)
    vertex = vertex';
end

[i,j,s] = find(sparse(A));
d = sqrt(sum( (vertex(i,:) - vertex(j,:)).^2, 2));
W = sparse(i,j,d);

end

author: Sable

repo: mcbench-benchmarks