text:
function cvx_optval = matrix_frac( x,Y )

%MATRIX_FRAC   Internal cvx version.

narginchk(2,2);
if ndims( Y ) > 2 || size( Y, 1 ) ~= size( Y, 2 ), %#ok

    error( 'Second argument must be a square matrix.' );

elseif ndims( x ) > 2 || size( x, 2 ) > 1, %#ok

    error( 'First argument must be a column vector.' );

elseif size( x, 1 ) ~= size( Y, 1 ),

    error( 'Size of first argument (vector) must match size of second argument (matrix).' );

elseif cvx_isconstant( x ) && cvx_isconstant( Y ),

    cvx_optval = cvx( matrix_frac( cvx_constant( x ), cvx_constant(Y) ) );

elseif cvx_isaffine( x ) && cvx_isaffine( Y ),

    n = size( x, 1 );
    z = [];
    cvx_begin
        epigraph variable z
        [Y x; x' z] == semidefinite( n+1 ); %#ok
    cvx_end

else

    error( 'Disciplined convex programming error:\n    MATRIX_FRAC is convex and nonmonotonic, so its input must be affine.' );

end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [aw,tt1, tt2, tmc, mag_zone]=bvalue(this, mcType, method)
    %BVALUE evaluate b-value, a-value and magnitude of completeness
    % of an earthquake catalog stored in a Catalog object.
    %
    % BVALUE(COBJ, MCTYPE) produces a Gutenberg-Richter type plot
    %    with the best fit line and display of b-,a-values and Mc
    %    for the catalog object COBJ. MCTYPE is a number from 1-5
    %    to select the algorithm used for calculation of the
    %    magnitude of completeness. Options are:
    %
    %    1: Maximum curvature
    %    2: Fixed Mc = minimum magnitude (Mmin)
    %    3: Mc90 (90% probability)
    %    4: Mc95 (95% probability)
    %    5: Best combination (Mc95 - Mc90 - maximum curvature)

    % Liberally adapted from original code in ZMAP.
    % Author: Silvio De Angelis, 27/07/2012 00:00:00
    % Modified and included in Catalog by Glenn Thompson,
    % 14/06/2014

    % This program is free software; you can redistribute it and/or modify
    % it under the terms of the GNU General Public License cobj.magas published by
    % the Free Software Foundation; either version 2 of the License, or
    % (at your option) any later version.
    %
    % This program is distributed in the hope that it will be useful,
    % but WITHOUT ANY WARRANTY; without even the implied warranty of
    % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    % GNU General Public License for more details.
    %
    % You should have received a copy of the GNU General Pucobj.magblic License
    % along with this program; if not, write to the
    % Free Software Foundation, Inc.,
    % 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

    if nargin < 2
        disp('--------------------------------------------------------')
        disp('ERROR: Usage is: bvalue(cobj, mcType). mcType not specified')
        disp('--------------------------------------------------------')
        disp('mcType can be:')
        disp('1: Maximum curvature')
        disp('2: Fixed Mc = minimum magnitude (Mmin)')
        disp('3: Mc90 (90% probability)')
        disp('4: Mc95 (95% probability)')
        disp('5: Best combination (Mc95 - Mc90 - maximum curvature)')
        return
    end

    % form magnitude vector - removing any NaN values with find
    good_magnitude_indices = find(this.data > 0.0);
    if strcmp(method, 'power')
        mag = log10(this.data(good_magnitude_indices));
    elseif strcmp(method, 'exponential')
        mag = this.data(good_magnitude_indices);
    end

    %MIN AND MAX MAGNITUDE IN CATALOG
    minimum_mag = min(mag);
    maximum_mag = max(mag);

    %COUNT EVENTS IN EACH MAGNITUDE BIN
    if strcmp(method, 'power')
        magrange = minimum_mag:0.1:maximum_mag;
    elseif strcmp(method, 'exponential')
        magrange = 10.^(log10(minimum_mag):0.1:log10(maximum_mag));
    end
    [bval, xt2] = hist(mag, magrange);

    %CUMULATIVE NUMBER OF EVENTS IN EACH MAGNITUDE BIN
    bvalsum = cumsum(bval);

    %NUMBER OF EVENTS IN EACH BIN IN REVERSE ORDER
    bval2 = bval(length(bval):-1:1);

    %NUMBER OF EVENTS IN EACH MAGNITUDE BIN IN REVERSE ORDER
    bvalsum3 = cumsum(bval(length(bval):-1:1));

    %BINS IN REVERSE ORDER
    xt3 = fliplr(magrange);
    backg_ab = log10(bvalsum3);

    %CREATE FIGURE WINDOW AND MAKE FREQUENCY-MAGNITUDE PLOT
    figure('Color','w','Position',[0 0 600 600])

    pl = semilogy(xt3,bvalsum3,'sb');

    set(pl, 'LineWidth', [1.0],'MarkerSize', [10],'MarkerFaceColor','r','MarkerEdgeColor','k');
    axis square
    hold on

    %pl1 = semilogy(xt3,bval2,'^b');
    %set(pl1, 'LineWidth',[1.0],'MarkerSize',[10],'MarkerFaceColor','w','MarkerEdgeColor','k');
    if strcmp(method, 'power')
        %xlabel('Log_1_0(Amplitude)','Fontsize', 12)
        xlabel('Magnitude','Fontsize', 12)
    elseif strcmp(method, 'exponential')
        xlabel('Amplitude','Fontsize', 12)
    end

    ylabel('Cumulative Minutes','Fontsize',12)
    set(gca,'visible','on','FontSize',12,'FontWeight','normal',...
        'FontWeight','bold','LineWidth',[1.0],'TickDir','in','Ticklength',[0.01 0.01],...
        'Box','on','Tag','cufi','color','w')

    %ESTIMATE B-VALUE (MAX LIKELIHOOD ESTIMATE)
    Nmin = 10;
    fMccorr = 0;
    fBinning = 0.1;

    if length(mag) >= Nmin

        %GOODNESS-OF-FIT TO POWER LAW
        %%%%%%%%%%%%%%%%%% mcperc_ca3.m start %%%%%%%%%%%%%%%%%%%%
        % This is a completeness determination test

        if strcmp(method, 'power')
            [bval,xt2] = hist(mag,-2:0.1:6);
        elseif strcmp(method, 'exponential')
            [bval,xt2] = hist(log10(mag),-2:0.1:6);
        end
        l = max(find(bval == max(bval)));
        magco0 =  xt2(l)

        dat = [];

        %for i = magco0-0.6:0.1:magco0+0.2
        for i = magco0-0.5:0.1:magco0+0.7
            if strcmp(method, 'power')
                l = mag >= i - 0.0499;
            elseif strcmp(method, 'exponential')
                l = mag >= 10^(i - 0.0499);
            end
            nu = length(mag(l));
            if length(mag(l)) >= 25;
                %[bv magco stan av] =  bvalca3(catZmap(l,:),2,2);
                if strcmp(method, 'power')
                    [mw bv2 stan2 av] =  bvalue_lib.bmemag(mag(l));
                elseif strcmp(method, 'exponential')
                    [mw bv2 stan2 av] =  bvalue_lib.bmemag(log10(mag(l)));
                end
                bvalue_lib.synthb_aut;
                dat = [ dat ; i res2];
            else
                dat = [ dat ; i nan];
            end

        end

        j =  min(find(dat(:,2) < 10 ));
        if isempty(j) == 1; Mc90 = nan ;
        else;
            Mc90 = dat(j,1);
        end

        j =  min(find(dat(:,2) < 5 ));
        if isempty(j) == 1; Mc95 = nan ;
        else;
            Mc95 = dat(j,1);
        end

        j =  min(find(dat(:,2) < 10 ));
        if isempty(j) == 1; j =  min(find(dat(:,2) < 15 )); end
        if isempty(j) == 1; j =  min(find(dat(:,2) < 20 )); end
        if isempty(j) == 1; j =  min(find(dat(:,2) < 25 )); end
        j2 =  min(find(dat(:,2) == min(dat(:,2)) ));
        %j = min([j j2]);

        Mc = dat(j,1);
        magco = Mc;
        prf = 100 - dat(j2,2);
        if isempty(magco) == 1; magco = nan; prf = 100 -min(dat(:,2)); end
        %display(['Completeness Mc: ' num2str(Mc) ]);
        %%%%%%%%%%%%%%%%%% mcperc_ca3.m end %%%%%%%%%%%%%%%%%%%%%%

        %CALCULATE MC
        [fMc] = bvalue_lib.calc_Mc(mag, mcType, fBinning, fMccorr);
        l = mag >= fMc-(fBinning/2);
        if length(mag(l)) >= Nmin
            [fMeanMag, fBValue, fStd_B, fAValue] =  bvalue_lib.calc_bmemag(mag(l), fBinning);
        else
            [fMc, fBValue, fStd_B, fAValue] = deal(NaN);
        end

        %STANDARD DEV OF a-value SET TO NAN;
        [fStd_A, fStd_Mc] = deal(NaN);

    else
        [fMc, fStd_Mc, fBValue, fStd_B, fAValue, fStd_A, ...
            fStdDevB, fStdDevMc] = deal(NaN);
    end

    magco = fMc; % magnitude of completeness?
    index_low=find(xt3 < magco+.05 & xt3 > magco-.05);
    mag_hi = xt3(1);
    index_hi = 1;
    mz = xt3 <= mag_hi & xt3 >= magco-.0001;
    mag_zone=xt3(mz);
    y = backg_ab(mz);

    %PLOT MC IN FIGURE
    Mc = semilogy(xt3(index_low),bvalsum3(index_low)*1.5,'vk');
    set(Mc,'LineWidth',[1.0],'MarkerSize',7)
    Mc = text(xt3(index_low)+0.2,bvalsum3(index_low)*1.5,'Mc');
    set(Mc,'FontWeight','normal','FontSize',12,'Color','k')

    %CREATE AND PLOT FIT LINE
    sol_type = 'Maximum Likelihood Solution';
    bw=fBValue;
    aw=fAValue;
    ew=fStd_B;
    p = [ -1*bw aw];
    f = polyval(p,mag_zone);
    f = 10.^f;

    hold on
    ttm= semilogy(mag_zone,f,'k');
    set(ttm,'LineWidth',[2.0])
    std_backg = ew;

    %ERROR CALCULATIONS
    %b = mag;
    bv = [];
    si = [];

    set(gca,'XLim',[min(mag)-0.5  max(mag+0.5)])
    %set(gca,'YLim',[0.9 length(mag+30)*2.5]);

    p=-p(1,1);
    p=fix(100*p)/100;
    tt1=num2str(bw,3);
    tt2=num2str(std_backg,1);
    tt4=num2str(bv,3);
    tt5=num2str(si,2);
    tmc=num2str(magco,2);
    rect=[0 0 1 1];
    h2=axes('position',rect);
    set(h2,'visible','off');
    a0 = aw-log10((max(this.dnum)-min(this.dnum))/365);

    text(.53,.88, ['b-value = ',tt1,' +/- ',tt2,',  a value = ',num2str(aw,3)],'FontSize',12);
    %text(.53,.85,sol_type,'FontSize',12 );
    text(.53,.82,['Magnitude of Completeness = ',tmc],'FontSize',12);

    % Glenn 20150111 add R^2 value
    thiscorr = corrcoef(mag_zone, f)
    r2 = thiscorr(1,2);
    thiscorr2 = corrcoef(mag_zone, log10(f))
    r22 = thiscorr2(1,2);
    %text(.53,.76,['R^2 = ',num2str(r2)],'FontSize',12);
    %text(.53,.70,['R^2 = ',num2str(r22)],'FontSize',12);

end

author: geoscience-community-codes

repo: GISMO


text:
%Combine forward and backward solutions

%The input arguments is highly dependent on the implementation. In the future I
%will probably change script with some other implementation independent
%version.

%dzBck=PbckInv*dx_bckwrd
function [P_smt, Cen_smt, h_smt, Vn_smt, qbn_smt, imu_smt, dxDif]=combine_2filt(fwdnav, fwdcov, Cen_bck, h_bck, Vn_bck, qbn_bck, dzBck, PbckInv)

%Forward filter results
Cen_fwd=reshape(fwdnav(1:9),[3,3]);
h_fwd=fwdnav(10); %elipsoid height
Vn_fwd=fwdnav(11:13);
qbn_fwd=fwdnav(14:17);
imu_fwd=fwdnav(18:end);

Pfwd=mat2vec_v000(fwdcov);
%%IMPORTATNT::In a real implementation you must replace inv operation with cholesky based inverse. DO NOT EVER compute inverse with "inv"
X=inv(eye(size(Pfwd))+Pfwd*PbckInv);
Wt=X*Pfwd;      %Note that (Maybeck, v2p10) uses the transpose of W. I did not understand why he prefers the transpose. (It seems it complicates the results. Could be because of some numeric properties?)

%Compute the difference between the forward solution and the nominal backward
%solution (we are going to estimate the errors on nominal forward solution
%using backward. We choose this approach because PfwdInv*dxFwd=0 as dxFwd=0)
dxDif=zeros(size(Pfwd,1),1);

%pos difference
Rappx=6378137/(sqrt(1.0-0.00669437999014*Cen_fwd(3,3)^2))+h_fwd;
vr_a=Rappx*Cen_fwd(1:2,:)*Cen_bck(3,:)';
dxDif(1:3)=[vr_a;h_bck-h_fwd];

%Vel dif
dxDif(4:6)=Vn_fwd-Vn_bck;

%Att dif
qnb_bck=[qbn_bck(1);-qbn_bck(2:4)];
vr_a=quatmult_v000(qbn_fwd, qnb_bck);
mx_a=quat2dcm_v000(vr_a);
dxDif(7)=mx_a(2,3);
dxDif(8)=mx_a(3,1);
dxDif(9)=mx_a(1,2);

%imu errors
dxDif(10:end)=imu_fwd;

%Smoothed covariance
%P_smt=inv(PfwdInv+PbckInv);
P_smt=X*Pfwd*X'+Wt*PbckInv*Wt';

%Smoothed error estimate on the forward results
dxSmt=P_smt*(PbckInv*dxDif+dzBck);

%Smoothed results
[qbn_smt, Vn_smt, Cen_smt, h_smt]=correctnav_Cen_v000(qbn_fwd, Vn_fwd, Cen_fwd, h_fwd, dxSmt(1:9), 2, 2);
imu_smt=imu_fwd-dxSmt(10:end);

repo: ekfmonoslam


text:
function L = log_prior(CPD)
% LOG_PRIOR Return log P(theta) for a generic CPD  - we return 0
% L = log_prior(CPD)

L = 0;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function sqi = rsqi(qrs,fs,CI,debug)
% compute smoothness of hrv or of hr time series given a confidence interval (CI).
% The underlying assumption for using this function is that the more smooth
% the qrs time series is the most likely it is to be a meaningful qrs time series.
%SMI = length(find(abs(hrv_CI)>30));

% inputs
%   qrs:    qrs fiducials (required, in data points number)
%   secDer: use second derivative? (i.e. HRV instead of HR to
%           compute SMI, default: 1)
%   fs:     sampling frequency (default: 1kHz)
%   CI:     confidence interval (default: 0.96)
%   segL:   length of the ecg segment in seconds (default: 60sec)
%
% outputs
%   sqi     percentage of intervals inside CI
%
% References:
% [1] Johnson, A. E. W., Behar, J., Andreotti, F., Clifford, G. D. and Oster, J. (2015).
% Multimodal heart beat detection using signal quality indices, Physiological Measurement
% 36 (2015): 1665-1677.
%
%
%
% Multimodal peak detection using ECG, ABP, PPG or SV
% Johnson, A. E. W., Behar, J., Andreotti, F., Clifford, G. D. and Oster, J.
%
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

if size(qrs,2) > size(qrs,1); qrs = qrs';end

if length(qrs) < 3, sqi = 0; return; end;

% == manage inputs
if nargin<1; error('rsqi: wrong number of input arguments \n'); end;
if nargin<2; fs = 1000; end;
if nargin<3; CI = 1; end;
if nargin<4; debug = 0; end;
if length(qrs) < 5;
    sqi = 0;
    disp('rSQI: too few qrs detection points')
    return
end

% == core function
%try
% == compute variability
hr = 60./(diff(qrs)/fs);

% rather than looking at the distribution of hr this option is looking at
% the distribution of hrv. This makes more sense because this way
% we are looking into high variation in deltaHR from a measure to the
% following one rather than the variability of absolute value of HR (which
% might be high if the foetus HR is changing significantly)

% now taking the derivative of smoothed hear rate. This will give the
% hrv
hrv = sort(diff(hr));
hrv_N = length(hrv);
% plot(hr); hold on, plot(yi,'r');
% we tolerate some mistakes using a confidence interval
if CI~=1
    CI_sup = ceil(hrv_N*(CI+(1-CI)/2));
    CI_inf = ceil(hrv_N*((1-CI)/2));
    hrv_CI = hrv(CI_inf:CI_sup);
else
    hrv_CI = hrv;
end
% output the std of the hrv
% SMI = std(hrv_CI); % OLD version

% new version (25-08-2013)
SMI = length(find(abs(hrv_CI)>20));
% this looks at the absolute number of outliers in hrv_CI with >
% 30bpm drop or increase from a point to the next.

sqi = 1 - SMI/hrv_N;
if sqi <0
    disp('What')
end

% == plots
if debug
    hist(hrv_CI,40); xlabel('hr or hrv histogram');
    title(['assess regularity in term of hr or hrv, REGULARITY:' SMI]);
    set(findall(gcf,'type','text'),'fontSize',14,'fontWeight','bold');
end

end

author: fernandoandreotti


text:
classdef DCF < dagnn.ElementWise
%DCF  layer
%   Discriminant Correlation Filters(DCF)
%
%   QiangWang, 2016
% -------------------------------------------------------------------------------------------------------------------------
    properties
        win_size = [3,3];
        sigma = 1;
    end
    properties (Transient)
        yf = [];
        lambda = 1e-4;
    end
    methods
        function outputs = forward(obj, inputs, params)

            xf = fft2(inputs{1});% target region
            zf = fft2(inputs{2});% search region
            xf_conj = conj(xf);
            [h,w,c,~] = size(xf);
            hwc = h*w*c;

            useGPU = isa(xf, 'gpuArray');
            if isempty(obj.yf)
                obj.initYF(useGPU);
            end

            kxxf = sum(xf .* xf_conj, 3) ./ hwc;
            alphaf = bsxfun(@rdivide, obj.yf, (kxxf + obj.lambda));
            kzxf = sum(zf .* xf_conj, 3) ./ hwc;
            outputs{1} = real(ifft2(alphaf .* kzxf));
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)

            dldrf = fft2(derOutputs{1});
            xf = fft2(inputs{1});% target region
            zf = fft2(inputs{2});% search region
            xf_conj = conj(xf);

            [h,w,c,~] = size(xf);
            hwc = h*w*c;

            kxxf = sum(xf .* xf_conj, 3) ./ hwc +obj.lambda;

            alphaf = bsxfun(@rdivide,obj.yf, kxxf);
            dldz = real(ifft2(bsxfun(@times,dldrf.*conj(alphaf),xf)))/hwc;
            kzxf = sum(zf .* xf_conj, 3) ./ hwc;
            dldx = real(ifft2(bsxfun(@times,conj(dldrf).*alphaf,zf)-...
            2*bsxfun(@times,xf,real(dldrf.*conj(alphaf.*kzxf)./kxxf))))/hwc;

            derInputs{1} = dldx;
            derInputs{2} = dldz;
            derParams = {};
        end

        function initYF(obj, useGPU)
            yf_ = single(fft2(gaussian_shaped_labels(obj.sigma, obj.win_size)));
            lambda_ = gather(obj.lambda);
            if useGPU
                obj.yf = gpuArray(yf_);
                obj.lambda = gpuArray(lambda_);
            else
                obj.yf = yf_;
                obj.lambda = lambda_;
            end
        end

        function obj = reset(obj)
            obj.yf = [] ;
            obj.lambda = 1e-4;
        end

        function obj = DCF(varargin)
            obj.load(varargin);
            obj.win_size = obj.win_size;
            obj.sigma = obj.sigma ;
        end
    end
end

function labels = gaussian_shaped_labels(sigma, sz)%kcf
[rs, cs] = ndgrid((1:sz(1)) - floor(sz(1)/2), (1:sz(2)) - floor(sz(2)/2));
labels = exp(-0.5 / sigma^2 * (rs.^2 + cs.^2));
labels = circshift(labels, -floor(sz(1:2) / 2) + 1);
assert(labels(1,1) == 1)
end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
%compute central speed in the central-head direction
function [data,units]=compute_velcentralch(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velcentralch=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    % this is just slightly faster
    velcentralch{i} = trx(larva).dxcentral_mm.*cos(trx(larva).centralheadang(1:end-1)) + ...
      trx(larva).dycentral_mm.*sin(trx(larva).centralheadang(1:end-1));
    %velcentralch1{1,i}=trx(larva).velmagcentral.*(cos(trx(larva).velangcentral).*cos(trx(larva).centralheadang(1,1:end-1))+sin(trx(larva).velangcentral).*sin(trx(larva).centralheadang(1,1:end-1)));
end

units=parseunits('mm/s');
data=velcentralch;

author: kristinbranson

repo: JAABA


text:
function sF = dot(sVF, a, varargin)
%
% Syntax
%   sVF = dot(sVF,sVF2)
%   sVF = dot(sVF,v)
%
% Input
%  sVF  - @S2VectorFieldTri
%  sVF2 - @S2VectorField
%  v    - @vector3d
%
% Output
%  sF - @S2FunTri
%

% first should be S2VectorFieldTri
if ~isa(sVF,'S2VectorFieldTri'), [sVF,a] = deal(a,sVF); end

if isa(a,'vector3d')
  sF = S2FunTri(sVF.tri, dot(sVF.values, a, varargin{:}));
elseif isa(a,'S2AxisFieldTri') || isa(a,'S2VectorFieldTri')
  sF = S2FunTri(sVF.tri, dot(sVF.values, b.values, varargin{:}));
else
   sF = S2FunTri(sVF.tri, dot(sVF.values, b.eval(sVF.vertices), varargin{:}));
end

end

author: mtex-toolbox

repo: mtex

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


text:
% eeglab_options() - handle EEGLAB options. This script (not function)
%                    set the various options in the eeg_options() file.
%
% Usage:
%   eeglab_options;
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2006-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, 2006-
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% load local file
% ---------------
homefolder = '';
try,
    %clear eeg_options; % note: we instead clear this function handle in pop_editoptions()

    if iseeglabdeployed
        com1 = readtxtfile(fullfile(eeglabexefolder, 'eeg_optionsbackup.txt'));
        com2 = readtxtfile(fullfile(eeglabexefolder, 'eeg_options.txt'));
        eval( com1 );
        eval( com2 );
    else
        eeg_optionsbackup;
        icadefs;

        % folder for eeg_options file (also update the pop_editoptions)
        if ~isempty(EEGOPTION_PATH)
             homefolder = EEGOPTION_PATH;
        elseif ispc
             if ~exist('evalc'), eval('evalc = @(x)(eval(x));'); end;
             homefolder = deblank(evalc('!echo %USERPROFILE%'));
        else homefolder = '~';
        end;

        option_file = fullfile(homefolder, 'eeg_options.m');
        oldp = pwd;
        try
            if ~isempty(dir(option_file))
                cd(homefolder);
            else
                tmpp2 = fileparts(which('eeglab_options.m'));
                cd(tmpp2);
            end;
        catch, end;
        eeg_options; % default one with EEGLAB
        cd(oldp);
    end;
    option_savematlab = ~option_savetwofiles;

catch
    lasterr
    disp('Warning: could not access the local eeg_options file');
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function sys=relaxdouble(X)
%RELAXDOUBLE Return numerial value treating nonlinear variables as independent

solution = yalmip('getsolution');
lmi_variables = X.lmi_variables;
opt_variables = solution.variables;

values = zeros(1+length(lmi_variables),1);
values(1)=1;
for i=1:length(lmi_variables)
    opt_index = find(lmi_variables(i)==opt_variables);
    if isempty(opt_index)
        values(i+1,1)=NaN;
    else
        values(i+1,1)=solution.optvar(opt_index);
    end
end

sys = X.basis*values;
sys = full(reshape(sys,X.dim(1),X.dim(2)));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [gx,dgdx,dgdp] = g_demo_susceptibility(Xt,P,ut,in)

%- initializations
idX1 = 1:in.n5(end);
nreg=length(in.n5);
nresp=length(in.r);

gxr=zeros(nresp,1);
dgdxr=zeros(length(Xt),nresp);
dgdpr=zeros(length(P),nresp);

%== compute gradients
[gxn,dgdxn,dgdpn] = g_HRF3(Xt(idX1),P,ut,in) ;
if nresp>0
[gxr,dgdxr(in.r,:),dgdpr] = g_dummy(Xt(in.r),P,ut,in);
else
    gxr=[];
    dgdxr=[];
    dgdpr=[];
end
%== concatenate BOLD and behavioral responses
%- state
gx = [gxn ; gxr] ;

%- Jacobian
dgdx = zeros(length(Xt),length(gx));
dgdx(idX1,1:nreg) = dgdxn ;
dgdx(:,nreg+(1:nresp)) = dgdxr ;

%- gradient wrt parameters
dgdp = zeros(length(P),nreg+nresp);
dgdp(1:size(dgdpn,1),1:nreg) = dgdpn ;
dgdp(:,nreg+(1:nresp)) = dgdpr ;

end

function [gxr,dgdxr,dgdpr]=g_dummy(Xt,P,ut,in)
gxr = Xt * 10;
dgdxr =  10;
dgdpr = 0;
end

author: MBB-team

repo: VBA-toolbox


text:
function tests = test_readstack(varargin)
%TEST_READSTACK  Test case for the file readstack
%
%   Test case for the file readstack

%   Example
%   test_readstack
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-09-27,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function test_TifList(testCase)
% Test call of function without argument

fName = fullfile('ratBrainMriSlices', 'BRNOR39e5p1-000.tif');
img = readstack(fName);

exp = [96 96 96];
assertEqual(testCase, exp, size(img));

function test_TifListVerbose(testCase)
% Test call of function without argument

fName = fullfile('ratBrainMriSlices', 'BRNOR39e5p1-000.tif');
img = readstack(fName, 'verbose');

exp = [96 96 96];
assertEqual(testCase, exp, size(img));

function test_TifListQuiet(testCase)
% Test call of function without argument

fName = fullfile('ratBrainMriSlices', 'BRNOR39e5p1-000.tif');
img = readstack(fName, 'quiet');

exp = [96 96 96];
assertEqual(testCase, exp, size(img));

function test_TifListRange(testCase)
% Test call of function without argument

fName = fullfile('ratBrainMriSlices', 'BRNOR39e5p1-000.tif');
range = 21:80;
img = readstack(fName, range);

exp = [96 96 length(range)];
assertEqual(testCase, exp, size(img));

function test_TifStack(testCase)
% Test call of function without argument

fName = fullfile('ratBrainMriSlices', 'BRNOR39e5p1stack.tif');
img = readstack(fName);

exp = [96 96 96];
assertEqual(testCase, exp, size(img));

function test_TifSlicesRGB(testCase)
% Test call of function without argument

fName = 'mriOverlay/mriOvr-000.tif';
img = readstack(fName);

exp = [128 128 3 27];
assertEqual(testCase, exp, size(img));

function test_TifStackRGB(testCase)
% Test call of function without argument

fName = 'mriOverlay.tif';
img = readstack(fName);

exp = [128 128 3 27];
assertEqual(testCase, exp, size(img));

function test_rawImage(testCase)

fileName = 'img_Int16_rampXYZ_6x5x4_msb.raw';
img = readstack(fullfile('mhd_images', fileName), [6 5 4], 'int16');

exp = [5 6 4];
assertEqual(testCase, exp, size(img));

function test_LSM(testCase)

fName = fullfile('images', '16c_col0_607.lsm');
img = readstack(fName);
exp = [512 512 89];
assertEqual(testCase, exp, size(img));

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function varargout=subsref(obj,subs)
% SUBSREF Subscripted reference
% An overloaded function...
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: subsref.m 4136 2010-12-09 22:22:28Z guillaume $

if isempty(subs), return; end

switch subs(1).type
    case '{}'
        error('Cell contents reference from a non-cell array object.');
    case '.'
        varargout = access_fields(obj,subs);
        return;
end

if numel(subs)~=1, error('Expression too complicated'); end;

dim  = [size(obj) ones(1,16)];
nd   = find(dim>1,1,'last')-1;
sobj = struct(obj);

if ~numel(subs.subs)
    [subs.subs{1:nd+1}] = deal(':');
elseif length(subs.subs) < nd
    l   = length(subs.subs);
    dim = [dim(1:(l-1)) prod(dim(l:end))];
    if numel(sobj) ~= 1
        error('Can only reshape simple file_array objects.');
    else
        if numel(sobj.scl_slope)>1 || numel(sobj.scl_inter)>1
            error('Can not reshape file_array objects with multiple slopes and intercepts.');
        end
        sobj.dim = dim;
    end
end

di   = ones(16,1);
args = cell(1,length(subs.subs));
for i=1:length(subs.subs)
    if ischar(subs.subs{i})
        if ~strcmp(subs.subs{i},':'), error('This shouldn''t happen....'); end
        if length(subs.subs) == 1
            args{i} = 1:prod(dim); % possible overflow when int32()
            k = 0;
            for j=1:length(sobj)
                sobj(j).dim = [prod(sobj(j).dim) 1];
                sobj(j).pos = [k+1 1];
                k           = k + sobj(j).dim(1);
            end
        else
            args{i} = 1:dim(i);
        end
    else
        args{i} = subs.subs{i};
    end
    di(i) = length(args{i});
end

if length(sobj)==1
    t = subfun(sobj,args{:});
else
    dt  = datatypes;
    dt  = dt([dt.code]==sobj(1).dtype); % assuming identical datatypes
    t = zeros(di',func2str(dt.conv));
    for j=1:length(sobj)
        ps = [sobj(j).pos ones(1,length(args))];
        dm = [sobj(j).dim ones(1,length(args))];
        for i=1:length(args)
            msk      = find(args{i}>=ps(i) & args{i}<(ps(i)+dm(i)));
            args2{i} = msk;
            args3{i} = double(args{i}(msk))-ps(i)+1;
        end

        t  = subsasgn(t,struct('type','()','subs',{args2}),subfun(sobj(j),args3{:}));
    end
end
varargout = {t};

%==========================================================================
% function t = subfun(sobj,varargin)
%==========================================================================
function t = subfun(sobj,varargin)

%sobj.dim = [sobj.dim ones(1,16)];
try
    args = cell(size(varargin));
    for i=1:length(varargin)
        args{i} = int32(varargin{i});
    end
    t = file2mat(sobj,args{:});
catch
    t = multifile2mat(sobj,varargin{:});
end
if ~isempty(sobj.scl_slope) || ~isempty(sobj.scl_inter)
    slope = 1;
    inter = 0;
    if ~isempty(sobj.scl_slope), slope = sobj.scl_slope; end
    if ~isempty(sobj.scl_inter), inter = sobj.scl_inter; end
    if numel(slope)>1
        slope = resize_scales(slope,sobj.dim,varargin);
        t     = double(t).*slope;
    else
        t     = double(t)*slope;
    end
    if numel(inter)>1
        inter = resize_scales(inter,sobj.dim,varargin);
    end;
    t = t + inter;
end

%==========================================================================
% function c = access_fields(obj,subs)
%==========================================================================
function c = access_fields(obj,subs)

sobj = struct(obj);
c    = cell(1,numel(sobj));
for i=1:numel(sobj)
    %obj = class(sobj(i),'file_array');
    obj = sobj(i);
    switch(subs(1).subs)
        case 'fname',      t = fname(obj);
        case 'dtype',      t = dtype(obj);
        case 'offset',     t = offset(obj);
        case 'dim',        t = dim(obj);
        case 'scl_slope',  t = scl_slope(obj);
        case 'scl_inter',  t = scl_inter(obj);
        case 'permission', t = permission(obj);
        otherwise
            error(['Reference to non-existent field "' subs(1).subs '".']);
    end
    if numel(subs)>1
        t = subsref(t,subs(2:end));
    end
    c{i} = t;
end

%==========================================================================
% function val = multifile2mat(sobj,varargin)
%==========================================================================
function val = multifile2mat(sobj,varargin)

% Convert subscripts into linear index
[indx2{1:length(varargin)}] = ndgrid(varargin{:},1);
ind = sub2ind(sobj.dim,indx2{:});

% Work out the partition
dt  = datatypes;
dt  = dt([dt.code]==sobj.dtype);
sz  = dt.size;
try
    mem = spm('Memory'); % in bytes, has to be a multiple of 16 (max([dt.size]))
catch
    mem = 200 * 1024 * 1024;
end
s   = ceil(prod(sobj.dim) * sz / mem);

% Assign indices to partitions
[x,y] = ind2sub([mem/sz s],ind(:));
c     = histc(y,1:s);
cc    = [0 reshape(cumsum(c),1,[])];

% Read data in relevant partitions
obj = sobj;
val = zeros(length(x),1,func2str(dt.conv));
for i=reshape(find(c),1,[])
    obj.offset = sobj.offset + mem*(i-1);
    obj.dim = [1 min(mem/sz, prod(sobj.dim)-(i-1)*mem/sz)];
    val(cc(i)+1:cc(i+1)) = file2mat(obj,int32(1),int32(x(y==i)));
end
r   = cellfun('length',varargin);
if numel(r) == 1, r = [r 1]; end
val = reshape(val,r);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
 function inv2 = de_ftab_inv2(fit, sl, varargin)
%function inv2 = de_ftab_inv2(fit, sl, [options])
%|
%| Build object that does polynomial inverse of BH function.
%| (To map log data into corresponding to line-integrals of material density.)
%|
%| in
%|	fit	strum	initialized by de_ftab_fit()
%|	sl	cell	{s1, ..., sL} see sls.sl from de_ftab_sls.m
%|
%| option
%|	'itype'	char	'T' or 'fixlinear' or 'fit' todo
%|				'fixlinear' is default which means
%|				use ftab.T for linear terms and then
%|				fit higher-order polynomial for rest
%|	'T'	[L M]	ftab.T - linearization at 0, if desired
%|	'show'	0|1	plot it?
%|
%| out
%|	inv2	strum
%|	methods:
%|		inv2.fun(hf)		map log values into corrected values
%|		inv2.plot(fit, sl)	show fits
%|
%| Copyright 2008-09-28, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(fit, 'test'), fail('run de_ftab test'), return, end

% option defaults
arg.emax = 4;
arg.T = [];
arg.esum = 6;
%arg.itype = 'T';
%arg.itype = 'fit'; % general fit
arg.itype = 'fixlinear'; % fix linear coefficients based on T, fit rest
%arg.wt = sll{2}.^3;	% weighting
%arg.wt = 1 ./ (0*sll{2} + 1);	% weighting
arg.show = false;

arg = vararg_pair(arg, varargin);

inv2 = de_ftab_inv2_setup(fit, sl, arg);

%
% de_ftab_inv2_setup()
% add methods to ftab related to polynomial approximation to its inverse
%
function inv2 = de_ftab_inv2_setup(fit, sl, arg)

LL = fit.LL;
MM = fit.MM;

if fit.MM < fit.LL
	warn('not making inv2() for MM=%d < LL=%d', fit.MM, fit.LL)
	inv2 = [];
return
end

switch arg.itype
case 'T' % 1st-order poly based on T
	if isempty(arg.T), fail 'need T', end
	coef = arg.T'; % [K L], K=M
	expo = eye(MM); % [K M], K=M

case {'fit', 'fixlinear'}
	sll = ndgrid_jf('mat', sl{:}); % [(Ns) L]
	sll = reshapee(sll, [], LL); % [*Ns L]
	ftmp = fit.fmfun(sll); % [*Ns M]

	expo = de_poly_expo(MM, 'emax', arg.emax, 'esum', arg.esum);
	expo(1,:) = []; % remove DC term
	basis = de_poly_eval(ftmp, [], expo, 'basis', 1); % [*Ns K]

% todo: consider weighting?
	if streq(arg.itype, 'fixlinear')
		if isempty(arg.T), fail 'need T', end
		jf_equal(expo(1:MM,1:MM), eye(MM)) % linear terms
		KK = size(expo,1); % # of coefficients
		coef = zeros(KK, LL);
		coef(1:MM,:) = arg.T'; % [M L]
		lin = de_poly_eval(ftmp, coef(1:MM,:), expo(1:MM,:)); % [*Ns L]
		basis = basis(:,(MM+1):end); % [*Ns K-M]
		coef((MM+1):end,:) = basis \ (sll - lin); % [K-M L]
	else
		coef = basis \ sll; % [K L]
	end

otherwise
	fail('bad itype %s', arg.itype)
end

arg.coef = coef;
arg.expo = expo;

if 1 % display polynomial fit errors
	tmp = de_poly_eval(ftmp, coef, expo);
	for ll=1:LL
		printm('worst inverse poly error l=%d: %g of %g', ll, ...
			max(abs(tmp(:,ll) - sll(:,ll))), max(sll(:,ll)))
	end
end

%
% old way: (todo: remove?)
% Build polynomial approximation to shat = (T F)^{-1}(z).
% Use weighted fitting so that soft tissue is fit well!
%

if 0	% this polynomial way seems insufficiently accurate!?
	ftab.inv.basis_func = ir_poly2_fun(3);
	tmp = ftab.inv.basis_func(z1(:), z2(:));
	ftab.inv.nbasis = length(ftab.inv.basis_func(0,0));
	ftab.inv.coef = zeros(ftab.inv.nbasis,2);
	for ll=1:ftab.LL
		ftab.inv.coef(:,ll) = (diag_sp(wt(:)) * tmp) \ ...
			(wt(:) .* sll{ll}(:));
	end
end

%ftab.inv.eval = @(ftab,fha) de_ftab_invert(ftab,fhat);

fun = @(inv2, fhat) de_ftab_inv2_eval(inv2, fhat);
meth = {'fun', fun, '(fhat [() M]) -> [() L]';
	'plot', @de_ftab_inv2_plot, '(fit, sl)'};
inv2 = strum(arg, meth);

%
% de_ftab_inv2_eval()
% todo
%
function shat = de_ftab_inv2_eval(inv2, fhat)
shat = de_poly_eval(fhat, inv2.coef, inv2.expo);

%
% de_ftab_inv2_plot()
% see how well a grid of s values is recovered
%
function out = de_ftab_inv2_plot(inv2, fit, sl)
LL = fit.LL;
if LL == 2
	sll = ndgrid_jf('mat', sl{:});
	ftmp = fit.fmfun(sll);
	shat = inv2.fun(ftmp);
	if im
		clf
	%	plot(sl{1}, sl{2}, shat, '.')
		plot(shat(:,:,1), shat(:,:,2), '.')
		axis tight
		xlabel 's1', ylabel 's2'
		title 's -> (fit) -> f -> (inv2) -> s'
	prompt
	end
else
	warn 'inv2_plot only for L=2'
end

if nargout, out = []; end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ fe, next, ierr ] = chfev ( x1, x2, f1, f2, d1, d2, ne, xe )

%*****************************************************************************80
%
%% CHFEV evaluates a cubic polynomial given in Hermite form.
%
%  Discussion:
%
%    This routine evaluates a cubic polynomial given in Hermite form at an
%    array of points.  While designed for use by SPLINE_PCHIP_VAL, it may
%    be useful directly as an evaluator for a piecewise cubic
%    Hermite function in applications, such as graphing, where
%    the interval is known in advance.
%
%    The cubic polynomial is determined by function values
%    F1, F2 and derivatives D1, D2 on the interval [X1,X2].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2005
%
%  Author:
%
%    Original FORTRAN77 version by Fred Fritsch.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Fred Fritsch, Ralph Carlson,
%    Monotone Piecewise Cubic Interpolation,
%    SIAM Journal on Numerical Analysis,
%    Volume 17, Number 2, April 1980, pages 238-246.
%
%    David Kahaner, Cleve Moler, Steven Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1988.
%
%  Parameters:
%
%    Input, real X1, X2, the endpoints of the interval of
%    definition of the cubic.  X1 and X2 must be distinct.
%
%    Input, real F1, F2, the values of the function at X1 and
%    X2, respectively.
%
%    Input, real D1, D2, the derivative values at X1 and
%    X2, respectively.
%
%    Input, integer NE, the number of evaluation points.
%
%    Input, real XE(NE), the points at which the function is to
%    be evaluated.  If any of the XE are outside the interval
%    [X1,X2], a warning error is returned in NEXT.
%
%    Output, real FE(NE), the value of the cubic function
%    at the points XE.
%
%    Output, integer NEXT(2), indicates the number of extrapolation points:
%    NEXT(1) = number of evaluation points to the left of interval.
%    NEXT(2) = number of evaluation points to the right of interval.
%
%    Output, integer IERR, error flag.
%    0, no errors.
%    -1, NE < 1.
%    -2, X1 == X2.
%
  if ( ne < 1 )
    ierr = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CHFEV - Fatal error!\n' );
    fprintf ( 1, '  Number of evaluation points is less than 1.\n' );
    fprintf ( 1, '  NE = %d\n', ne );
    error ( 'CHFEV - Fatal error!' )
  end

  h = x2 - x1;

  if ( h == 0.0 )
    ierr = -2;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CHFEV - Fatal error!\n' );
    fprintf ( 1, '  The interval [X1,X2] is of zero length.\n' );
    error ( 'CHFEV - Fatal error!' )
  end
%
%  Initialize.
%
  ierr = 0;
  next(1) = 0;
  next(2) = 0;
  xmi = min ( 0.0, h );
  xma = max ( 0.0, h );
%
%  Compute cubic coefficients expanded about X1.
%
  delta = ( f2 - f1 ) / h;
  del1 = ( d1 - delta ) / h;
  del2 = ( d2 - delta ) / h;
  c2 = -( del1 + del1 + del2 );
  c3 = ( del1 + del2 ) / h;
%
%  Evaluation loop.
%
  for i = 1 : ne

    x = xe(i) - x1;
    fe(i) = f1 + x * ( d1 + x * ( c2 + x * c3 ) );
%
%  Count the extrapolation points.
%
    if ( x < xmi )
      next(1) = next(1) + 1
    end

    if ( xma < x )
      next(2) = next(2) + 1
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = oselect(ob, y)
% recall that imask/omask and idim/odim are swapped on ctranspose
out = fatrix2_select(ob.omask, ob.odim, y);


text:
% reformatting frame_turn
figure;plot(frame_turn(:,17))
ylim([-10,200])
hold on;plot(frame_turn(:,18),'r')
%%
st2 = frame_turn(:,18);
IX_start = find(diff([1;st2])>10)
IX_stop = find(diff([1;st2])<-10)
% IX_stop = [IX_stop,8248];
%%
for i = 1:2,
    IX = IX_start(i):IX_stop(i)-1;
    frame_turn(IX,17) = frame_turn(IX,17)+40;
end

author: xiuyechen

repo: FishExplorer


text:
function gen_hermite_poly_test ( )

%******************************hermite****************************************80
%
%% GEN_HERMITE_POLY_TEST tests GEN_HERMITE_POLY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2015
%
%  Author:
%
%    John Burkardt
%
  i_test = 6;
  n = 10;

  mu_test = [ 0.0, 0.0, 0.1, 0.1, 0.5, 1.0 ];
  x_test = [ 0.0, 1.0, 0.0, 0.5, 0.5, 0.5 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEN_HERMITE_POLY_TEST\n' );
  fprintf ( 1, '  GEN_HERMITE_POLY evaluates the generalized Hermite \n' );
  fprintf ( 1, '  polynomials.\n' );

  for i = 1 : i_test

    x = x_test(i);
    mu = mu_test(i);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Table of L(N,MU)(X) for\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '    N(max) = %d\n', n );
    fprintf ( 1, '    MU =     %f\n', mu );
    fprintf ( 1, '    X =      %f\n', x );
    fprintf ( 1, '\n' );

    c = gen_hermite_poly ( n, x, mu );

    for j = 0 : n
      fprintf ( 1, '  %4d  %12f\n', j, c(j+1) );
    end

  end

  return
end


text:
function handles = changeFontsize(handles, change)

currentFontsize = handles.fontsizePanels;
newFontsize = currentFontsize + change;

% Allow max fontsize 22 and min fontsize 6:
maxFontsize = 22;
minFontsize = 6;

if ( (newFontsize > maxFontsize) || (newFontsize < minFontsize) )
    % Maximum or minimum fontsize reached -- do nothing.
    return
end

% Get the field names in CHEBGUI:
names = fieldnames(handles);

% Find all the elements whos font size we want to change:
textLocs = strfind(names, 'text_');
inputLocs = strfind(names, 'input_');
buttonLocs = strfind(names, 'button_');
panelLocs = strfind(names, 'panel_');
toggleLocs = strfind(names, 'toggle_');
popupLocs = strfind(names, 'popupmenu_');
iterLocs = strfind(names, 'iter_');
editLocs = strfind(names, 'edit_');

% Combine all the locations of elements whose font we wish to change
allLocs = strcat(textLocs, inputLocs, buttonLocs, panelLocs, toggleLocs, ...
    popupLocs, iterLocs, editLocs);

% Deal with variable maximum for input and noninput-type handles
for fieldCounter = 1:length(textLocs)
    if ( ~isempty(allLocs{fieldCounter}) )
        % Access the field values dynamically using the .( ) call
        currFS = get(handles.(names{fieldCounter}), 'FontSize');
        newFS = currFS + change;
        set(handles.(names{fieldCounter}), 'FontSize', newFS)
    end
end

% Change font sizes of the plots:
set(handles.fig_sol,  'fontsize', newFontsize)
set(handles.fig_norm, 'fontsize', newFontsize)

% Store the new font size:
handles.fontsizePanels = newFontsize;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [vMags, vClusTime]=cluslengthcomp(vmain, vcluster, mCatalog)
% function cluslengthcomp(vmain, vcluster, mCatalog);
% -------------------------------------------------------------
% Compare actual cluster length with applied windowing technique
% using the ECOS example
%
% Incoming variables:
% vmain: Vector of mainshocks
% vcluster : Vector with cluster numbers
% mCatalog : EQ catalog in ZMAP format
%
% Outgoing variable:
% vClusTime : Vector of length of Cluster
% E.G. tmp=load('/home/jowoe/PhD/Decluster/Data/ecos_all_cluster');
% J. Woessner
% last update: 30.07.02

mTmpCat = [vmain vcluster mCatalog];
%mTmpCat = mTmpCat(1:1000,:);
% Select cluster
for nCevent = 1: max(mTmpCat(:,1))
    vSelClus = (mTmpCat(:,2) == nCevent);
    vSelMain = (mTmpCat(:,1) == nCevent);
    mTmpCat2 = mTmpCat(vSelClus,:);
%     if isempty(mTmpCat2)
%         disp('oo');
%     end
    mTmpCat3 = mTmpCat(vSelMain,:);
    vMinClusTime(nCevent) = min(mTmpCat2(:,5));
    vMaxClusTime(nCevent) = max(mTmpCat2(:,5));
    vMainTime(nCevent) = mTmpCat3(:,5); % If this does not work, then there is more than one mainshock, which doesn't make sense
    if vMainTime(nCevent) < vMinClusTime(nCevent)
        vMinClusTime(nCevent) = vMainTime(nCevent);
    end % END of if-check for minimum date of cluster
    if vMainTime(nCevent) > vMaxClusTime(nCevent)  % This should actually not happen
        vMaxClusTime(nCevent) = vMainTime(nCevent);
    end % END of if-check for maximum date of cluster
    vClusTime(nCevent) = (vMaxClusTime(nCevent)-vMinClusTime(nCevent))*365;
    vMags(nCevent) = mTmpCat3(:,8);
end % End of FOR over nCevent
figure;
%[i]=find(vClusTime > 10.^(0.5409*vMags-0.547)) % Gardner & Knopoff
%[i]=find(vClusTime > 10.^(0.5055*vMags-0.1329)) % Gruenthal, 1985
%[i]=find(vClusTime > exp(-3.95+sqrt(0.62+17.32*vMags))) % Gruenthal, pers

%vMags(i) = NaN ;
semilogy(vMags,vClusTime,'*');
hold on;
vMagnitude = (0:0.1:10);
vMagnitudea = (0:0.1:6.5);
vMagnitudeb = (6.5:0.1:10);

% Gardner & Knopoff
vTimeGaKn74 = 10.^(0.5409*vMagnitudea-0.547);
vTimeGaKn74b = 10.^(0.032*vMagnitudeb+2.7389); % M>=6.5
semilogy(vMagnitudea,vTimeGaKn74,'Color',[1 0 0],'Linewidth', 2);
semilogy(vMagnitudeb,vTimeGaKn74b,'Color',[1 0 0],'Linewidth', 2);

% % Gruenthal, 1985
% vSpaceGr85 = 10.^(0.1060*vMagnitude+1.0982);
% vTimeGr85 = 10.^(0.5055*vMagnitude-0.1329);
% semilogy(vMagnitude,vTimeGr85,'Color',[1 0 0],'Linewidth', 2);

% Gruenthal, pers. communication
% vSpaceGr = exp(1.77+sqrt(0.037+1.02*vMagnitude));
% vTimeGra = exp(-3.95+sqrt(0.62+17.32*vMagnitudea));
% vTimeGrb = 10.^(2.8+0.024*vMagnitudeb); % M >= 6.5
% semilogy(vMagnitudea,vTimeGra,'Color',[0 0.8 0],'Linewidth', 2);
% semilogy(vMagnitudeb,vTimeGrb,'Color',[0 0.8 0],'Linewidth', 2);

set(gca,'Xlim', [0 8]);
xlabel('Magnitude');
ylabel('Time / [days]');

author: CelsoReyes


text:
function s=str(obj)
% Return string representation of MoxUnitTestCase instance
%
% s=str(obj)
%
% Input:
%   obj             MoxUnitTestCase instance.
%
% Output:
%   s               string representation of obj.
%

    s=sprintf('%s(%s,%s)>', class(obj), getName(obj), getLocation(obj));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [proj_lg, BHCalib] = BHCorrection(datafolder, geo, ScanXML, proj_lg,gpuids)
% Entry Function For BH Correction
%   Detailed explanation goes here

disp('BH correction is very memory intensive. The effect is to be justified. ');
typein = input('Continue BH correction? Y or N (recommended): ', 's');
if(contains(typein, 'n', 'IgnoreCase', true))
    BHCalib = NaN;
    disp('BH correction is skipped.');
    return;
end

disp('Beam Hardening Correction is on-going: be patient... ');

% Key calibration information
BHCalib = BHCalibFromXML(datafolder, ScanXML);

% Precompute filter attenuated spectrum: debug pass
BHCalib = BH_SpectrumFilter(BHCalib);

% Precompute bowtie attenuated spectra LUT
BHCalib = BH_SpectrumBowtieLUT(geo, BHCalib);

% Build reference object (water) attanuation LUT
BHCalib = BH_ObjectCalibLUT(BHCalib);

% BH correction via reference object (water)
BHCalib = BH_RemappingFunc(BHCalib);

proj_lg = BH_ObjectRemapping(BHCalib, proj_lg, gpuids);

disp('BH correction is done.')

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function [ip, deltaGPn, iter] = point_scene_to_image(S, sicd_meta, varargin)
%POINT_SCENE_TO_IMAGE Transforms a 3D ECF point to pixel row, col
% This function implements the SICD Image Projections Description Document:
% http://www.gwg.nga.mil/ntb/baseline/docs/SICD/index.html
%
% [ip, deltaGPn, iter] = point_scene_to_image(S, sicd_meta, 'PropertyName', PropertyValue, ...)
%
% Inputs:
%    S            - [3xN] N Ground Points in ECF coordinates
%    sicd_meta    - SICD meta data structure
%
%       Property name     Description
%       delta_gp_max      Ground plane displacement tol (m), default = quarter pixel
%       delta_arp         ARP position adjustable parameter (ECF, m).  Default 0.
%       delta_varp        VARP position adjustable parameter (ECF, m/s).  Default 0.
%       range_bias        Range bias adjustable parameter (m).  Default 0.
%       adj_params_frame  Coordinate frame used for expressing delta_arp
%                         and delta_varp adjustable parameters.  Allowed
%                         values: 'ECF', 'RIC_ECF', 'RIC_ECI'. Default ECF.
%
% Outputs:
%    ip           - [2xN] (row; column) coordinates of N points in image
%                   (or subimage if FirstRow/FirstCol are nonzero).
%                   Zero-based, following SICD convention (rather than
%                   MATLAB convention, which is one-based); that is,
%                   upper-left pixel is [0;0].
%    deltaGPn     - Residual ground plane displacement (m)
%    iter         - # iterations required
%
% Authors: Thomas McDowall, Harris Corporation
%          Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Extract the relevant SICD info
% Convert ground / scene point to image coordinates
SCP_Row = double(sicd_meta.ImageData.SCPPixel.Row);
SCP_Col = double(sicd_meta.ImageData.SCPPixel.Col);

First_Row = double(sicd_meta.ImageData.FirstRow);
First_Col = double(sicd_meta.ImageData.FirstCol);
% We override this so that point_image_to_ground will work in coordinate
% space of full image. We will add it back in at the end.
sicd_meta.ImageData.FirstRow = 0;
sicd_meta.ImageData.FirstCol = 0;

Row_SS = sicd_meta.Grid.Row.SS;
Col_SS = sicd_meta.Grid.Col.SS;

% Parse input parameters
p = inputParser;
p.addParamValue('delta_gp_max',0.25 * sqrt(Row_SS^2 + Col_SS^2), @isscalar); % meters
% Adjustable parameters
p.addParamValue('delta_arp',[0 0 0], @(x) numel(x)==3); % ECF? (meters)
p.addParamValue('delta_varp',[0 0 0], @(x) numel(x)==3); % ECF? (meters/s)
p.addParamValue('range_bias',0, @isscalar); % meters
p.addParamValue('adj_params_frame','ECF', @(x) any(strcmpi(x,{'ECF','RIC_ECF','RIC_ECI'})));
p.FunctionName = mfilename;
p.parse(varargin{:});

uRow = [sicd_meta.Grid.Row.UVectECF.X;...
    sicd_meta.Grid.Row.UVectECF.Y;...
    sicd_meta.Grid.Row.UVectECF.Z];

uCol = [sicd_meta.Grid.Col.UVectECF.X;...
    sicd_meta.Grid.Col.UVectECF.Y;...
    sicd_meta.Grid.Col.UVectECF.Z];

%% 3.1 SCP Projection Equations
SCP = [sicd_meta.GeoData.SCP.ECF.X; sicd_meta.GeoData.SCP.ECF.Y; sicd_meta.GeoData.SCP.ECF.Z];

ARP_SCP_COA = [sicd_meta.SCPCOA.ARPPos.X; sicd_meta.SCPCOA.ARPPos.Y; sicd_meta.SCPCOA.ARPPos.Z];
VARP_SCP_COA = [sicd_meta.SCPCOA.ARPVel.X; sicd_meta.SCPCOA.ARPVel.Y; sicd_meta.SCPCOA.ARPVel.Z];

ARP_SCP_COAminusSCP = ARP_SCP_COA - SCP;
%R_SCP_COA = norm(ARP_SCP_COAminusSCP);
%Rdot_SCP_COA = (1/R_SCP_COA) * dot(VARP_SCP_COA, ARP_SCP_COAminusSCP);

% Normal to instantaneous slant plane that contains SCP at SCP COA is
% tangent to R/Rdot contour at SCP. Points away from center of Earth. Use
% LOOK to establish sign.
if strcmp(sicd_meta.SCPCOA.SideOfTrack, 'L')
    LOOK = 1;
else
    LOOK = -1;
end
SPN_SCP_COA = LOOK * cross(VARP_SCP_COA, -ARP_SCP_COAminusSCP);
uSPN = SPN_SCP_COA./repmat(sqrt(sum(SPN_SCP_COA.^2)),3,1);

%% 6.1 Scene To Image: Single Scene Point

% Spherical earth ground plane normal (exact orientation of plane is not
% critical, as long as it contains S)
if isvector(S), S = S(:); end; % Assure orientation
uGPN = S./repmat(sqrt(sum(S.^2)),3,1);

% Ground plane points are projected along straight lines to the image plane. The
% GP to IP direction is along the SCP COA slant plane normal. Also, compute
% image plane unit normal, uIPN. Compute projection scale factor SF.
uPROJ = uSPN;
IPN = cross(uRow, uCol); % should match SICD.PFA.IPN for PFA data
uIPN = IPN./repmat(sqrt(sum(IPN.^2)),3,1); % pointing away from center of earth
SF = dot(uPROJ, uIPN);

% Initialize
Gn = S;

% 2.4 - Image Plane parameters
% The following section is for ground to image with non-orthogonal axes.
% theta col is angle between uRow and uCol if not 0.
cos_theta_col = dot(uRow, uCol);
sin_theta_col = sqrt(1 - cos_theta_col.^2);
GI = (1 / sin_theta_col^2) * [1 -cos_theta_col; -cos_theta_col 1];

done = false;

% Iterate the ground to image transform
iter = ones(1,size(S,2));
to_iter = true(1,size(S,2));
ip = zeros(2,size(S,2));
deltaPn = zeros(3,size(S,2));
deltaGPn = zeros(1,size(S,2));
while(~done && all(iter < 6))
    % Project ground plane point Gn to image plane point In. The projection
    % distance is DISTn. Compute image coordinates xrow and ycol.
    DISTn = (1/SF) * dot((repmat(SCP,1,size(Gn,2)) - Gn), repmat(uIPN,1,size(Gn,2)));
    In = Gn + repmat(DISTn,3,1) .* repmat(uPROJ,1,size(DISTn,2));

    % For a point at IPP, corresponding SCP pixel-centered coord is
    deltaIPP = In - repmat(SCP,1,size(In,2));
    ip_iter = GI * ...
        [dot(deltaIPP, repmat(uRow,1,size(In,2))); ...
        dot(deltaIPP, repmat(uCol,1,size(In,2)))];

    xrow = ip_iter(1,:);
    ycol = ip_iter(2,:);

    irow = xrow / Row_SS;
    icol = ycol / Col_SS;

    row = irow + SCP_Row;
    col = icol + SCP_Col;

    ip(:,to_iter) = [row; col];

    % Transform to ground plane containing the scene point(s) S
    Pn = point_image_to_ground(ip(:,to_iter), sicd_meta, ...
        'projection_type', 'plane', 'gref', S(:,to_iter), 'ugpn', uGPN(:,to_iter), ...
        ... % Pass through adjustable parameters
        'delta_arp', p.Results.delta_arp, ...
        'delta_varp', p.Results.delta_varp, ...
        'range_bias', p.Results.range_bias, ...
        'adj_params_frame', p.Results.adj_params_frame);

    %   gndPt = ecfToGeodetic(Pn');
    %   fprintf('Ground Point: %f %f %f\n', gndPt(1), gndPt(2), gndPt(3));

    % Compute displacement between Pn and S
    deltaPn(:,to_iter) = S(:,to_iter) - Pn;
    deltaGPn(to_iter) = sqrt(sum(deltaPn(:,to_iter).^2));

    old_iter = to_iter;
    to_iter = deltaGPn > p.Results.delta_gp_max; % Need to iterate further on these points
    if any(to_iter)
        Gn = Gn(:,to_iter(old_iter)) + deltaPn(:,to_iter);
        iter(to_iter) = iter(to_iter) + 1;
    else
        done = true;
    end

end
if iter > 5
    warning('point_scene_to_image:SolutionDidNotConverge',...
        'Solution for image coordinate did not converge.');
end

ip(1,:) = ip(1,:) - First_Row;
ip(2,:) = ip(2,:) - First_Col;

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [haus, meanAbsSurfDist] = calc_HausdorffMetric(structNum1,structNum2,planC)
%function [haus, meanAbsSurfDist] = calc_HausdorffMetric(structNum1,structNum2,planC)
%
%This function computes Hausdorff distance between structNum1 and structNum2
%
%APA, 11/20/2014

if ~exist('planC','var')
    global planC
end

[~, x1V, y1V, z1V, planC] = getStructSurface(structNum1,planC);
[~, x2V, y2V, z2V, planC] = getStructSurface(structNum2,planC);

haus = NaN;
meanAbsSurfDist = NaN;
if ~isempty(x1V) && ~isempty(x1V)
    [haus,meanAbsSurfDist] = hausdorff([x1V(:) y1V(:) z1V(:)],[x2V(:) y2V(:) z2V(:)]);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
[X,Y,Z] = meshgrid(-2:.07:2,-2:.07:2,-2:.07:2);
V = X.*exp(-X.^2-Y.^2-Z.^2);

author: Sable

repo: mcbench-benchmarks


text:
%GENTRUNK Generation of Trunk's classification problem of 2 Gaussian classes
%
%   A = GENTRUNK(N,K)
%
% INPUT
%   N       Dataset size, or 2-element array of class sizes (default: [50 50]).
%   K       Dimensionality of the dataset to be generated (default: 10).
%
% OUTPUT
%   A       Dataset.
%
% DESCRIPTION
% Generation of a K-dimensional 2-class dataset A of N objects. Both classes
% are Gaussian distributed with the idenity matrix as covariance matrix.
% The means of the first class are defined by ua(j) = 1/sqrt(j). The means
% for the second class are ub = -ua. These means are such that the Nearest
% Mean Classifier always shows peaking for a finite training set.
%
% REFERENCES
% 1. G.V. Trunk, A Problem of Dimensionality: A Simple Example, IEEE Trans.
% Pattern Analysis and Machine Intelligence, vol. 1, pp. 306-307, 1979
% 2. A.K. Jain, R.P.W. Duin, and J. Mao, Statistical Pattern Recognition:
% A Review, IEEE Transactions on Pattern Analysis and Machine Intelligence,
% vol. 22, pp. 4-37, 2000.
%
% EXAMPLE
% a = gentrunk([1000 1000],200);
% e = clevalf(a,nmc,[1:9 10:5:25 50:25:200],[5 5],25);
% plote(e)
%
% SEE ALSO (PRTools Guide)
% DATASETS, PRDATASETS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: gentrunk.m,v 1.2 2009/01/27 13:01:42 duin Exp $

function A = gendats (N,k)

		if (nargin < 1), N = [50 50]; end
	if (nargin < 2), k = 10;	end

	% Set equal priors and generate random class sizes according to these.
	p = [0.5 0.5]; N = genclass(N,p);

	% Unit covariance matrices
	% GA = eye(k); GB = eye(k);
  % Trunk means
  ma = 1./sqrt(1:k);
  mb = -ma;

  A = prdataset([randn(N(1),k) + repmat(ma,N(1),1); randn(N(2),k) + repmat(mb,N(2),1)]);
  A = prdataset(A,genlab(N));
  A = setprior(A,p);

% 	U = prdataset([ma;mb],[1 2]');
% 	U = setprior(U,p);
%
% 	% Create dataset.
% 	A = gendatgauss(N,U,cat(3,GA,GB));
	A = setname(A,'Trunk''s Problem');

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function out = EN_MS_LZcomplexity(y,n,preProc)
% EN_MS_LZcomplexity Lempel-Ziv complexity of a n-bit encoding of a time series
%
%---INPUTS:
% y, the input time series
% n, the (integer) number of bits to encode the data into
% preProc [opt], first apply a given preProcessing to the time series. For now,
%               just 'diff' is implemented, which zscores incremental
%               differences and then applies the complexity method.
%
%---OUTPUT: the normalized Lempel-Ziv complexity: i.e., the number of distinct
%           symbol sequences in the time series divided by the expected number
%           of distinct symbols for a noise sequence.

% Uses Michael Small's code: 'complexity' (renamed MS_complexity here).
%
% cf. M. Small, Applied Nonlinear Time Series Analysis: Applications in Physics,
% Physiology, and Finance (book) World Scientific, Nonlinear Science Series A,
% Vol. 52 (2005)
% Code is available at http://small.eie.polyu.edu.hk/matlab/
%
% The code is a wrapper for Michael Small's original code and uses the
% associated mex file compiled from complexitybs.c (renamed MS_complexitybs.c
% here).
% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2 || isempty(n)
    n = 2; % n-bit encoding
end
if nargin < 3
    preProc = []; % no preprocessing
end

% Apply some pre-processing to the time series before performing the analysis
if ischar(preProc)
    switch preProc
    case 'diff'
        y = zscore(diff(y));
    otherwise
        error('Unknown preprocessing setting ''%s''', preProc);
    end
end

% Run Michael Small's (mexed) code for calcaulting the Lempel-Ziv complexity:
out = MS_complexity(y,n);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ vX ] = ProjectL1BallDual( vY, ballRadius, vLowerBound, vUpperBound )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = ProjectL1Ball( vY, ballRadius, vLowerBound, vUpperBound )
%   Solving the Orthogonal Porjection Problem of the input vector onto the
%   L1 Ball with Box Constraints using Dual Function.
% Input:
%   - vY            -   Input Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - ballRadius    -   Ball Radius.
%                       Sets the Radius of the L1 Ball. For Unit L1 Ball
%                       set to 1.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
%   - vLowerBound   -   Lower Bound Vector.
%                       Sets the lower bound values of the solution
%                       (Element wise).
%                       Structure: Vector.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - vUpperBound   -   Upper Bound Vector.
%                       Sets the upper bound values of the solution
%                       (Element wise).
%                       Structure: Vector.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - vX            -   Output Vector.
%                       The projection of the Input Vector onto the L1
%                       Ball.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  https://math.stackexchange.com/a/2830242/33.
% Remarks:
%   1.  S
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     24/06/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

DEBUG_MODE = OFF;

paramLambda     = 0; %<! Initialization value

% Check feasibility of the problem
minSum = 0;
for ii = 1:size(vY, 1)
    if(sign(vLowerBound(ii)) == sign(vUpperBound(ii)))
        % Zero isn't within the boundary (Or both are zero)
        minSum = minSum + min(abs(vLowerBound(ii)), abs(vUpperBound(ii)));
    end
end

if(minSum > ballRadius)
    % The problem is infeasible
    vX = mean([vLowerBound, vUpperBound], 2);
    return;
end

% The dual objective function which should be maximized to find the optimal
% 'paramLambda'.
% The functios is negated as we want to maximize while using
% MATLAB's minimization function
hObjFun = @(paramLambda) -ObjectiveDualFunction(vY, ballRadius, vLowerBound, vUpperBound, paramLambda); %<! Objective function

if(DEBUG_MODE == ON)
    numSamples = 1000;

    vParamLambda    = linspace(0, 5, numSamples);
    vObjVal         = zeros([numSamples, 1]);

    for ii = 1:numSamples
        vObjVal(ii) = hObjFun(vParamLambda(ii));
    end

    figure();
    plot(vParamLambda, vObjVal);
end

paramLambda = fminsearch(hObjFun, paramLambda); %<! Objective Function isn't smooth hence can't use Newton Method
paramLambda = max(paramLambda, 0);

% Solution as derived for the updated forumaltion of the problem.
[~, vX] = ObjectiveDualFunction(vY, ballRadius, vLowerBound, vUpperBound, paramLambda);

end

function [ valObj, vX ] = ObjectiveDualFunction( vY, ballRadius, vL, vU, paramLambda )

vX      = zeros([size(vY, 1), 1]);

for ii = 1:size(vY, 1)

    if(sign(vL(ii)) == sign(vU(ii)))
        valX = ProjBoxFunction(vY(ii) - (paramLambda * sign(vL(ii))), vL(ii), vU(ii));
    elseif(sign(vY(ii)) == sign(vL(ii)))
        % Implictily vL(ii) <= 0, vU(ii) >= 0 and vY(ii) <= 0
        valX = ProjBoxFunction(vY(ii) + paramLambda, vL(ii), 0); %<! Making sure sign of valX is non positive
    elseif(sign(vY(ii)) == sign(vU(ii)))
        % Implictily vU(ii) >= 0, vL(ii) <= 0 and vY(ii) >= 0
        valX = ProjBoxFunction(vY(ii) - paramLambda, 0, vU(ii)); %<! Making sure sign of valX is non negative
    end

    vX(ii) = valX;

end

valObj = (0.5 * sum((vX - vY) .^ 2)) + (paramLambda * (sum(abs(vX)) - ballRadius));

end

function [ outVal ] = ProjBoxFunction( valIn, valL, valU )

outVal = min(max(valIn, valL), valU);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [out] = analyseThermoConstrainedModel(model, cumNormProbCutoff, printLevel, resultsBaseFileName)
% USAGE:
%
%    [out] = analyseThermoConstrainedModel(model, cumNormProbCutoff, printLevel, resultsBaseFileName)
%
% INPUTS:
%    model:                  structure with fields:
%
%                              * .DfGt0 - `m x 1` array of estimated standard transformed Gibbs
%                                energies of formation.
%                              * .DrGt0 - `n x 1` array of estimated standard transformed
%                                reaction Gibbs energies.
%                              * .DfGtMin - `m x 1` array of estimated lower bounds on
%                                transformed Gibbs energies of formation.
%                              * .DfGtMax - `m x 1` array of estimated upper bounds on
%                                transformed Gibbs energies of formation.
%                              * .DrGtMin - `n x 1` array of estimated lower bounds on
%                                transformed reaction Gibbs energies.
%                              * .DrGtMax - `n x 1` array of estimated upper bounds on
%                                transformed reaction Gibbs energies.
%                              * .quantDir - `n x 1` array indicating quantitatively assigned
%                                reaction directionality. 1 for reactions that are
%                                irreversible in the forward direction, -1 for
%                                reactions that are irreversible in the reverse
%                                direction, and 0 for reversible reactions.
%    cumNormProbCutoff:      default = 0.2
%    printLevel:             verbose level, default = 1
%    resultsBaseFileName:
%
% OUTPUT:
%    out:

out=[];

if ~exist('printLevel','var')
    printLevel=1;
end
if ~exist('cumNormProbCutoff','var')
    cumNormProbCutoff=0.2;
else
    if cumNormProbCutoff<0
        error('cumNormProbCutoff cannot be less than zero')
    end
end

if 1
    %statistics of directions
    fprintf('%s\n','directionalityStats...');
    model=directionalityStats(model,cumNormProbCutoff,printLevel);
end
%OUTPUT
% model.directions    a structue of boolean vectors with different directionality
%               assignments where some vectors contain subsets of others
%
% qualtiative -> quantiative changed reaction directions
%   .directions.forwardForward
%   .directions.forwardReverse
%   .directions.forwardReversible
%   .directions.reversibleFwd
%   .directions.reversibleRev
%   .directions.reversibleReversible
%   .directions.tightened
%
% subsets of qualtiatively forward  -> quantiatively reversible
%   .directions.forwardReversible_bydGt0
%   .directions.forwardReversible_bydGt0LHS
%   .directions.forwardReversible_bydGt0Mid
%   .directions.forwardReversible_bydGt0RHS
%
%   .directions.forwardReversible_byConc_zero_fixed_DrG0
%   .directions.forwardReversible_byConc_negative_fixed_DrG0
%   .directions.forwardReversible_byConc_positive_fixed_DrG0
%   .directions.forwardReversible_byConc_negative_uncertain_DrG0
%   .directions.forwardReversible_byConc_positive_uncertain_DrG0

if 1
    %report on directionality changes
    fprintf('%s\n','directionalityChangeReport...');
    directionalityChangeReport(model,cumNormProbCutoff,printLevel,resultsBaseFileName)
end

% if 0
%     fprintf('\n%s\n','...standardGibbsFormationEnergyStats');
%     [nKeq,nGC,nNone]=standardGibbsFormationEnergyStats(model,figures);
% end

if 0
    %pie charts with proportions of reaction directionalities and changes in
    %directionality
    fprintf('%s\n','directionalityStatFigures...');
    directionalityStatsFigures(model.directions,resultsBaseFileName)
end

%figures of reaction directionality changes
%qualitatively forward now quantiatiavely reversible
if any(model.directions.forward2Reversible)
    fprintf('%s\n','forwardReversibleFigures...');
    if 0
        forwardReversibleFigures(model)
    else
        thorStandard=0;
        forwardReversibleFiguresCC(model,model.directions,thorStandard)
    end
end

return

%create a vertical errorbar figure of the qualitatively forward transport reactions
%that are quantitatively reversible, whether from group contribution or
%Keq, that then need to be assigned to be forward, to limit the growth rate
% i.e. abc transporters or reactions involving protons
if any(model.directions.forwardReversible & model.transportRxnBool)
     forwardTransportQuantReverseFigures(model)
end

%vertical errorbar of qualitatively forward but quantitatively reverse
%forward, probably reverse
%   .directions.forwardReversible_bydGt0RHS
%   .directions.forwardReversible_byConc_positive_fixed_DrG0
%   .directions.forwardReversible_byConc_positive_uncertain_DrG0
fwdProbReverse=model.directions.forwardReversible_bydGt0RHS | ...
    model.directions.forwardReversible_byConc_positive_uncertain_DrG0| ...
    model.directions.forwardReversible_byConc_positive_fixed_DrG0;
if any(fwdProbReverse) && 0
    forwardProbablyReverseFigures(model,directions)
end

author: opencobra

repo: cobratoolbox


text:
function circle2 = transformCircle3d(circle, tfm)
%TRANSFORMCIRCLE3D Transform a 3D circle with a 3D affine transformation.
%
%   CIRCLE2 = transformPlane3d(CIRCLE, TFM)
%
%   Example
%     circle = [1 1 1 2 45 45 0];
%     tfm = createRotationOz(pi);
%     circle2 = transformCircle3d(circle, tfm);
%     figure('color','w'); hold on; axis equal tight; view(-10,25);
%     xlabel('x'); ylabel('y'); zlabel('z');
%     drawCircle3d(circle,'r'); drawPoint3d(circle(1:3),'r+')
%     drawCircle3d(circle2,'g'); drawPoint3d(circle2(1:3),'g+')
%
%   See also
%     transforms3d, transformPoint3d, transformVector3d, transformLine3d,
%     transformPlane3d
%

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2022-12-03, using MATLAB 9.13.0.2080170 (R2022b) Update 1
% Copyright 2022

parser = inputParser;
addRequired(parser, 'circle', @(x) validateattributes(x, {'numeric'},...
    {'ncols',7,'real','finite','nonnan'}));
addRequired(parser, 'tfm', @isTransform3d);
parse(parser, circle, tfm);
circle = parser.Results.circle;
tfm = parser.Results.tfm;

% Compute transformation from local basis to world basis
initialTfm = localToGlobal3d(circle(1:3), circle(5), circle(6), circle(7));
% Add the additional transformation
newTfm = tfm*initialTfm;

% Convert to Euler angles
[phi, theta, psi] = rotation3dToEulerAngles(newTfm, 'ZYZ');

% Create transformed circle
circle2 = [transformPoint3d(circle(1:3), tfm), circle(4), theta, phi, psi];

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [b, mpe] = back1_mpe(engine, bfuture, f, ev1, t)

if t ~= 1
  error('mixed up time stamps')
end
bnet = bnet_from_engine(engine);
ss = bnet.nnodes_per_slice;
maximize = 1;

int = engine.interface;
D = engine.in_clq; % from J2
C = engine.int_clq1; % from J1
phiD = marginalize_pot(bfuture.clpot{D}, int, maximize);
phiC = marginalize_pot(f.clpot{C}, int, maximize);
ratio = divide_by_pot(phiD, phiC);
f.clpot{C} = multiply_by_pot(f.clpot{C}, ratio);

[mpe, b.clpot] = find_max_config(engine.jtree_engine1, f.clpot, f.seppot, ev1);
for c=1:length(b.clpot)
  [b.clpot{c}, ll(c)] = normalize_pot(b.clpot{c});
end
b.t = t;

author: bayesnet

repo: bnt


text:
classdef Bicubic < dagnn.ElementWise
%%% Bicubic resizing layer %%%
%
% Performs bicubic resizing on inputs{1} and stores it in outputs{1}.
% outputs{1} =imresize(input, obj.scale,'bicubic');
% *Back-propagation (backward function) not implemented*

    properties
        scale = 2;
    end

    methods
        function outputs = forward(obj, inputs, params)
            outputs{1} = imresize(inputs{1}, obj.scale, 'bicubic');
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            derInputs{1} = 0*inputs{1} ;
            derParams = {} ;
        end

        function obj = Bicubic(varargin)
            obj.load(varargin);
            obj.scale = obj.scale;
        end
    end
end

author: sooyekim

repo: Deep-SR-ITM


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the logistic
%  regression exercise. You will need to complete the following functions
%  in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the exam scores and the third column
%  contains the label.

data = load('ex2data1.txt');
X = data(:, [1, 2]); y = data(:, 3);

%% ==================== Part 1: Plotting ====================
%  We start the exercise by first plotting the data to understand the
%  the problem we are working with.

fprintf(['Plotting data with + indicating (y = 1) examples and o ' ...
         'indicating (y = 0) examples.\n']);

plotData(X, y);

% Put some labels
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============ Part 2: Compute Cost and Gradient ============
%  In this part of the exercise, you will implement the cost and gradient
%  for logistic regression. You neeed to complete the code in
%  costFunction.m

%  Setup the data matrix appropriately, and add ones for the intercept term
[m, n] = size(X);

% Add intercept term to x and X_test
X = [ones(m, 1) X];

% Initialize fitting parameters
initial_theta = zeros(n + 1, 1);

% Compute and display initial cost and gradient
[cost, grad] = costFunction(initial_theta, X, y);

fprintf('Cost at initial theta (zeros): %f\n', cost);
fprintf('Gradient at initial theta (zeros): \n');
fprintf(' %f \n', grad);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 3: Optimizing using fminunc  =============
%  In this exercise, you will use a built-in function (fminunc) to find the
%  optimal parameters theta.

%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);

%  Run fminunc to obtain the optimal theta
%  This function will return theta and the cost
[theta, cost] = ...
	fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);

% Print theta to screen
fprintf('Cost at theta found by fminunc: %f\n', cost);
fprintf('theta: \n');
fprintf(' %f \n', theta);

% Plot Boundary
plotDecisionBoundary(theta, X, y);

% Put some labels
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============== Part 4: Predict and Accuracies ==============
%  After learning the parameters, you'll like to use it to predict the outcomes
%  on unseen data. In this part, you will use the logistic regression model
%  to predict the probability that a student with score 45 on exam 1 and
%  score 85 on exam 2 will be admitted.
%
%  Furthermore, you will compute the training and test set accuracies of
%  our model.
%
%  Your task is to complete the code in predict.m

%  Predict probability for a student with score 45 on exam 1
%  and score 85 on exam 2

prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n\n'], prob);

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

author: sfvsfv

repo: Mathematical-modeling


text:
function ori = byAxisAngle(axis,angle,varargin)
% define orientations by rotational axis and rotational angle
%
% Syntax
%   ori = orientation.byAxisAngle(v,omega,CS) % an orientation
%   ori = orientation.byAxisAngle(h,omega) % a misorientation
%
% Input
%  v     - rotational axis @vector3d
%  m     - rotational axis @Miller
%  omega - rotation angle
%  CS    - @crystalSymmetry
%
% Output
%  ori - @orientation
%
% See also
% orientation/orientation orientation/byEuler orientation/byMatrix orientation/map

% TODO: better implementation required, call byAxisAngle@rotattion

ori = orientation(axis2quat(axis,angle),varargin{:});

% copy crystal symmetry if possible
if isa(axis,'Miller'), ori.CS = axis.CS; ori.SS = axis.CS; end

author: mtex-toolbox

repo: mtex


text:
function I = shift_reconstruct(Y,shifts,diffphase,us_fac,Nr,Nc,Np,method,add_value)

% applies 3-d sub-pixel shifts to an input image
% INPUTS:
% Y:            input image (double 3d tensor) in space (real) or frequency (complex) domain
% shifts:       shifts
% diffphase:    phase difference
% us_fac:       upsampling factor for subpixel shifts
% method:       method for treating boundaries
% add_value:    value to add when zero-ing out boundaries

% OUTPUT:
% I:        output image

% Written by Eftychios A. Pnevmatikakis, Simons Foundation, 2016

if isreal(Y);
    buf2ft = fftn(Y);
else
    buf2ft = Y;
end

if nargin < 9 || isempty(add_value); add_value = 0; end
if nargin < 8 || isempty(method); method = 'zero'; end
if nargin < 4 || isempty(us_fac); us_fac = 50; end

[nr,nc,np]=size(Y);
if any(shifts)
    if nargin < 5 || isempty(Nr); Nr = ifftshift(-fix(nr/2):ceil(nr/2)-1); end
    if nargin < 6 || isempty(Nc); Nc = ifftshift(-fix(nc/2):ceil(nc/2)-1); end
    if nargin < 7 || isempty(Np); Np = ifftshift(-fix(np/2):ceil(np/2)-1); end

    %shifts = shiftdim(shifts,3);

    row_shift = shifts(1);
    col_shift = shifts(2);
    if ismatrix(Y); pln_shift = 0; else pln_shift = shifts(3); end
    shifts = [row_shift,col_shift,pln_shift];

    if us_fac > 0
        if isvector(Nc); [Nc,Nr,Np] = meshgrid(Nc,Nr,Np); end
        Greg = buf2ft.*exp(1i*2*pi*(-row_shift*Nr/nr-col_shift*Nc/nc-pln_shift*Np/np));
    elseif us_fac == 0
        Greg = buf2ft;
    end
    Greg = Greg*exp(1i*diffphase);
    I = real(ifftn(Greg));
    I = remove_boundaries(I,shifts,method,add_value);
else
    if isreal(Y)
        I = Y;
    else
        I = real(ifftn(Y));
    end
end

end

author: flatironinstitute

repo: NoRMCorre


text:
function mortrage()
clc;clear all; close all;
% definitions
% P = principle amount
% J = Monthly intrest rate
% R = Total monthly payment
% R = J*P+ amount applied to principle
% new principle amount
% P+j*P-R => P(1+J)-R => P*m-R: m=1+J

% inputs
pri=input('Principle amount ');
Yr=input(' No of Year ' );
rate=input('Annual percentage rate ' );
peryear=1/12;
percent=1/100;
ttl_mnt=Yr*12;
syms m J P R A N;

% the principal after n payments can be written as
% P = A?mn ? R? (mn ? 1)/(m? 1).
solve(A*m^N - R*(m^N - 1)/(m - 1), R);
R = subs(ans, m, J + 1);
format bank;
 disp( '     Interest Rate      Payment')
%for rate= 6:0.2:12
    disp([rate, double(subs(R, [A, N, J], [pri, ttl_mnt, rate*percent*peryear]))])
%end
ttl_payment=subs(R, [A, N, J], [pri, ttl_mnt, rate*percent*peryear])*ttl_mnt;
disp('total amunt paid =')
disp(ttl_payment)


text:
function fx = p02_fun ( n, x )

%*****************************************************************************80
%
%% P02_FUN evaluates the integrand for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gwynne Evans,
%    Practical Numerical Integration,
%    Wiley, 1993,
%    ISBN: 047193898X,
%    LC: QA299.3E93.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(2,N), the evaluation points.
%
%    Output, real FX(N,1), the integrand values.
%
  fx(1:n,1) = 1.0 ./ sqrt ( 1.0 - x(1,1:n).^2 .* x(2,1:n).^2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function result = ssc_mc_omp(X, D, K)
    solver = spx.cluster.ssc.SSC_TOMP(X, D, K);
    result = solver.solve();
end

author: indigits

repo: sparse-plex


text:
options.filter = 'linear';
options.filter = '9-7';

%% 1D %%
n = 1024;
x = linspace(0,1,n)';
x = cos(10*pi*x) + (x>.4);

Jmin = 6;

options.ti = 0;
y = perform_lifting_transform(x, Jmin, +1, options);
x1 = perform_lifting_transform(y, Jmin, -1, options);
disp(['Error(should be 0)=' num2str(norm(x-x1)/norm(x))]);

options.ti = 1;
y = perform_lifting_transform(x, Jmin, +1, options);
x1 = perform_lifting_transform(y, Jmin, -1, options);
disp(['Error(should be 0)=' num2str(norm(x-x1)/norm(x))]);

%% test for RWT
options.ti = 1;
options.use_mex = 1;
options.wavelet_vm = 3;
options.wavelet_type = 'daubechies';
y = perform_wavelet_transform(x, Jmin, +1, options);
x1 = perform_wavelet_transform(y, Jmin, -1, options);
disp(['Error(should be 0)=' num2str(norm(x-x1)/norm(x))]);

%% test for wavelab
options.ti = 1;
options.use_mex = 0;
options.wavelet_vm = 3;
options.wavelet_type = 'daubechies';
y = perform_wavelet_transform(x, Jmin, +1, options);
x1 = perform_wavelet_transform(y, Jmin, -1, options);
disp(['Error(should be 0)=' num2str(norm(x-x1)/norm(x))]);

return;

%% 2D %%
n = 256;
Jmin = 3;
x = load_image('lena', n);
x = rescale(x);

options.ti = 0;
y = perform_lifting_transform(x, Jmin, +1, options);
x1 = perform_lifting_transform(y, Jmin, -1, options);
disp(['Error(should be 0)=' num2str(norm(x-x1)/norm(x))]);

options.ti = 1;
y = perform_lifting_transform(x, Jmin, +1, options);
x1 = perform_lifting_transform(y, Jmin, -1, options);
disp(['Error(should be 0)=' num2str(norm(x-x1)/norm(x))]);

author: gpeyre

repo: matlab-toolboxes


text:
%  internal function

%  'xform_nii.m' is an internal function called by "load_nii.m", so
%  you do not need run this program by yourself. It does simplified
%  NIfTI sform/qform affine transform, and supports some of the
%  affine transforms, including translation, reflection, and
%  orthogonal rotation (N*90 degree).
%
%  For other affine transforms, e.g. any degree rotation, shearing
%  etc. you will have to use the included 'reslice_nii.m' program
%  to reslice the image volume. 'reslice_nii.m' is not called by
%  any other program, and you have to run 'reslice_nii.m' explicitly
%  for those NIfTI files that you want to reslice them.
%
%  Since 'xform_nii.m' does not involve any interpolation or any
%  slice change, the original image volume is supposed to be
%  untouched, although it is translated, reflected, or even
%  orthogonally rotated, based on the affine matrix in the
%  NIfTI header.
%
%  However, the affine matrix in the header of a lot NIfTI files
%  contain slightly non-orthogonal rotation. Therefore, optional
%  input parameter 'tolerance' is used to allow some distortion
%  in the loaded image for any non-orthogonal rotation or shearing
%  of NIfTI affine matrix. If you set 'tolerance' to 0, it means
%  that you do not allow any distortion. If you set 'tolerance' to
%  1, it means that you do not care any distortion. The image will
%  fail to be loaded if it can not be tolerated. The tolerance will
%  be set to 0.1 (10%), if it is default or empty.
%
%  Because 'reslice_nii.m' has to perform 3D interpolation, it can
%  be slow depending on image size and affine matrix in the header.
%
%  After you perform the affine transform, the 'nii' structure
%  generated from 'xform_nii.m' or new NIfTI file created from
%  'reslice_nii.m' will be in RAS orientation, i.e. X axis from
%  Left to Right, Y axis from Posterior to Anterior, and Z axis
%  from Inferior to Superior.
%
%  NOTE: This function should be called immediately after load_nii.
%
%  Usage: [ nii ] = xform_nii(nii, [tolerance], [preferredForm])
%
%  nii	- NIFTI structure (returned from load_nii)
%
%  tolerance (optional) - distortion allowed for non-orthogonal rotation
%	or shearing in NIfTI affine matrix. It will be set to 0.1 (10%),
%	if it is default or empty.
%
%  preferredForm (optional)  -  selects which transformation from voxels
%	to RAS coordinates; values are s,q,S,Q.  Lower case s,q indicate
%	"prefer sform or qform, but use others if preferred not present".
%	Upper case indicate the program is forced to use the specificied
%	tranform or fail loading.  'preferredForm' will be 's', if it is
%	default or empty.	- Jeff Gunter
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function nii = xform_nii(nii, tolerance, preferredForm)

   %  save a copy of the header as it was loaded.  This is the
   %  header before any sform, qform manipulation is done.
   %
   nii.original.hdr = nii.hdr;

   if ~exist('tolerance','var') | isempty(tolerance)
      tolerance = 0.1;
   elseif(tolerance<=0)
      tolerance = eps;
   end

   if ~exist('preferredForm','var') | isempty(preferredForm)
      preferredForm= 's';				% Jeff
   end

   %  if scl_slope field is nonzero, then each voxel value in the
   %  dataset should be scaled as: y = scl_slope * x + scl_inter
   %  I bring it here because hdr will be modified by change_hdr.
   %
   if nii.hdr.dime.scl_slope ~= 0 & ...
	ismember(nii.hdr.dime.datatype, [2,4,8,16,64,256,512,768]) & ...
	(nii.hdr.dime.scl_slope ~= 1 | nii.hdr.dime.scl_inter ~= 0)

      nii.img = ...
	nii.hdr.dime.scl_slope * double(nii.img) + nii.hdr.dime.scl_inter;

      if nii.hdr.dime.datatype == 64

         nii.hdr.dime.datatype = 64;
         nii.hdr.dime.bitpix = 64;
      else
         nii.img = single(nii.img);

         nii.hdr.dime.datatype = 16;
         nii.hdr.dime.bitpix = 32;
      end

      nii.hdr.dime.glmax = max(double(nii.img(:)));
      nii.hdr.dime.glmin = min(double(nii.img(:)));

      %  set scale to non-use, because it is applied in xform_nii
      %
      nii.hdr.dime.scl_slope = 0;

   end

   %  However, the scaling is to be ignored if datatype is DT_RGB24.

   %  If datatype is a complex type, then the scaling is to be applied
   %  to both the real and imaginary parts.
   %
   if nii.hdr.dime.scl_slope ~= 0 & ...
	ismember(nii.hdr.dime.datatype, [32,1792])

      nii.img = ...
	nii.hdr.dime.scl_slope * double(nii.img) + nii.hdr.dime.scl_inter;

      if nii.hdr.dime.datatype == 32
         nii.img = single(nii.img);
      end

      nii.hdr.dime.glmax = max(double(nii.img(:)));
      nii.hdr.dime.glmin = min(double(nii.img(:)));

      %  set scale to non-use, because it is applied in xform_nii
      %
      nii.hdr.dime.scl_slope = 0;

   end

   %  There is no need for this program to transform Analyze data
   %
   if nii.filetype == 0 & exist([nii.fileprefix '.mat'],'file')
      load([nii.fileprefix '.mat']);	% old SPM affine matrix
      R=M(1:3,1:3);
      T=M(1:3,4);
      T=R*ones(3,1)+T;
      M(1:3,4)=T;
      nii.hdr.hist.qform_code=0;
      nii.hdr.hist.sform_code=1;
      nii.hdr.hist.srow_x=M(1,:);
      nii.hdr.hist.srow_y=M(2,:);
      nii.hdr.hist.srow_z=M(3,:);
   elseif nii.filetype == 0
      nii.hdr.hist.rot_orient = [];
      nii.hdr.hist.flip_orient = [];
      return;				% no sform/qform for Analyze format
   end

   hdr = nii.hdr;

   [hdr,orient]=change_hdr(hdr,tolerance,preferredForm);

   %  flip and/or rotate image data
   %
   if ~isequal(orient, [1 2 3])

      old_dim = hdr.dime.dim([2:4]);

      %  More than 1 time frame
      %
      if ndims(nii.img) > 3
         pattern = 1:prod(old_dim);
      else
         pattern = [];
      end

      if ~isempty(pattern)
         pattern = reshape(pattern, old_dim);
      end

      %  calculate for rotation after flip
      %
      rot_orient = mod(orient + 2, 3) + 1;

      %  do flip:
      %
      flip_orient = orient - rot_orient;

      for i = 1:3
         if flip_orient(i)
            if ~isempty(pattern)
               pattern = flipdim(pattern, i);
            else
               nii.img = flipdim(nii.img, i);
            end
         end
      end

      %  get index of orient (rotate inversely)
      %
      [tmp rot_orient] = sort(rot_orient);

      new_dim = old_dim;
      new_dim = new_dim(rot_orient);
      hdr.dime.dim([2:4]) = new_dim;

      new_pixdim = hdr.dime.pixdim([2:4]);
      new_pixdim = new_pixdim(rot_orient);
      hdr.dime.pixdim([2:4]) = new_pixdim;

      %  re-calculate originator
      %
      tmp = hdr.hist.originator([1:3]);
      tmp = tmp(rot_orient);
      flip_orient = flip_orient(rot_orient);

      for i = 1:3
         if flip_orient(i) & ~isequal(tmp(i), 0)
            tmp(i) = new_dim(i) - tmp(i) + 1;
         end
      end

      hdr.hist.originator([1:3]) = tmp;
      hdr.hist.rot_orient = rot_orient;
      hdr.hist.flip_orient = flip_orient;

      %  do rotation:
      %
      if ~isempty(pattern)
         pattern = permute(pattern, rot_orient);
         pattern = pattern(:);

         if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 | ...
		hdr.dime.datatype == 128 | hdr.dime.datatype == 511

            tmp = reshape(nii.img(:,:,:,1), [prod(new_dim) hdr.dime.dim(5:8)]);
            tmp = tmp(pattern, :);
            nii.img(:,:,:,1) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);

            tmp = reshape(nii.img(:,:,:,2), [prod(new_dim) hdr.dime.dim(5:8)]);
            tmp = tmp(pattern, :);
            nii.img(:,:,:,2) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);

            if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
               tmp = reshape(nii.img(:,:,:,3), [prod(new_dim) hdr.dime.dim(5:8)]);
               tmp = tmp(pattern, :);
               nii.img(:,:,:,3) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);
            end

         else
            nii.img = reshape(nii.img, [prod(new_dim) hdr.dime.dim(5:8)]);
            nii.img = nii.img(pattern, :);
            nii.img = reshape(nii.img, [new_dim       hdr.dime.dim(5:8)]);
         end
      else
         if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 | ...
		hdr.dime.datatype == 128 | hdr.dime.datatype == 511

            nii.img(:,:,:,1) = permute(nii.img(:,:,:,1), rot_orient);
            nii.img(:,:,:,2) = permute(nii.img(:,:,:,2), rot_orient);

            if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
               nii.img(:,:,:,3) = permute(nii.img(:,:,:,3), rot_orient);
            end
         else
            nii.img = permute(nii.img, rot_orient);
         end
      end
   else
      hdr.hist.rot_orient = [];
      hdr.hist.flip_orient = [];
   end

   nii.hdr = hdr;

   return;					% xform_nii

%-----------------------------------------------------------------------
function [hdr, orient] = change_hdr(hdr, tolerance, preferredForm)

   orient = [1 2 3];
   affine_transform = 1;

   %  NIFTI can have both sform and qform transform. This program
   %  will check sform_code prior to qform_code by default.
   %
   %  If user specifys "preferredForm", user can then choose the
   %  priority.					- Jeff
   %
   useForm=[];					% Jeff

   if isequal(preferredForm,'S')
       if isequal(hdr.hist.sform_code,0)
           error('User requires sform, sform not set in header');
       else
           useForm='s';
       end
   end						% Jeff

   if isequal(preferredForm,'Q')
       if isequal(hdr.hist.qform_code,0)
           error('User requires qform, qform not set in header');
       else
           useForm='q';
       end
   end						% Jeff

   if isequal(preferredForm,'s')
       if hdr.hist.sform_code > 0
           useForm='s';
       elseif hdr.hist.qform_code > 0
           useForm='q';
       end
   end						% Jeff

   if isequal(preferredForm,'q')
       if hdr.hist.qform_code > 0
           useForm='q';
       elseif hdr.hist.sform_code > 0
           useForm='s';
       end
   end						% Jeff

   if isequal(useForm,'s')
      R = [hdr.hist.srow_x(1:3)
           hdr.hist.srow_y(1:3)
           hdr.hist.srow_z(1:3)];

      T = [hdr.hist.srow_x(4)
           hdr.hist.srow_y(4)
           hdr.hist.srow_z(4)];

      if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
         hdr.hist.old_affine = [ [R;[0 0 0]] [T;1] ];
         R_sort = sort(abs(R(:)));
         R( find( abs(R) < tolerance*min(R_sort(end-2:end)) ) ) = 0;
         hdr.hist.new_affine = [ [R;[0 0 0]] [T;1] ];

         if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
            msg = [char(10) char(10) '   Non-orthogonal rotation or shearing '];
            msg = [msg 'found inside the affine matrix' char(10)];
            msg = [msg '   in this NIfTI file. You have 3 options:' char(10) char(10)];
            msg = [msg '   1. Using included ''reslice_nii.m'' program to reslice the NIfTI' char(10)];
            msg = [msg '      file. I strongly recommand this, because it will not cause' char(10)];
            msg = [msg '      negative effect, as long as you remember not to do slice' char(10)];
            msg = [msg '      time correction after using ''reslice_nii.m''.' char(10) char(10)];
            msg = [msg '   2. Using included ''load_untouch_nii.m'' program to load image' char(10)];
            msg = [msg '      without applying any affine geometric transformation or' char(10)];
            msg = [msg '      voxel intensity scaling. This is only for people who want' char(10)];
            msg = [msg '      to do some image processing regardless of image orientation' char(10)];
            msg = [msg '      and to save data back with the same NIfTI header.' char(10) char(10)];
            msg = [msg '   3. Increasing the tolerance to allow more distortion in loaded' char(10)];
            msg = [msg '      image, but I don''t suggest this.' char(10) char(10)];
            msg = [msg '   To get help, please type:' char(10) char(10) '   help reslice_nii.m' char(10)];
            msg = [msg '   help load_untouch_nii.m' char(10) '   help load_nii.m'];
            error(msg);
         end
      end

   elseif isequal(useForm,'q')
      b = hdr.hist.quatern_b;
      c = hdr.hist.quatern_c;
      d = hdr.hist.quatern_d;

      if 1.0-(b*b+c*c+d*d) < 0
         if abs(1.0-(b*b+c*c+d*d)) < 1e-5
            a = 0;
         else
            error('Incorrect quaternion values in this NIFTI data.');
         end
      else
         a = sqrt(1.0-(b*b+c*c+d*d));
      end

      qfac = hdr.dime.pixdim(1);
      if qfac==0, qfac = 1; end
      i = hdr.dime.pixdim(2);
      j = hdr.dime.pixdim(3);
      k = qfac * hdr.dime.pixdim(4);

      R = [a*a+b*b-c*c-d*d     2*b*c-2*a*d        2*b*d+2*a*c
           2*b*c+2*a*d         a*a+c*c-b*b-d*d    2*c*d-2*a*b
           2*b*d-2*a*c         2*c*d+2*a*b        a*a+d*d-c*c-b*b];

      T = [hdr.hist.qoffset_x
           hdr.hist.qoffset_y
           hdr.hist.qoffset_z];

      %  qforms are expected to generate rotation matrices R which are
      %  det(R) = 1; we'll make sure that happens.
      %
      %  now we make the same checks as were done above for sform data
      %  BUT we do it on a transform that is in terms of voxels not mm;
      %  after we figure out the angles and squash them to closest
      %  rectilinear direction. After that, the voxel sizes are then
      %  added.
      %
      %  This part is modified by Jeff Gunter.
      %
      if det(R) == 0 | ~isequal(R(find(R)), sum(R)')

         %  det(R) == 0 is not a common trigger for this ---
         %  R(find(R)) is a list of non-zero elements in R; if that
         %  is straight (not oblique) then it should be the same as
         %  columnwise summation. Could just as well have checked the
         %  lengths of R(find(R)) and sum(R)' (which should be 3)
         %
         hdr.hist.old_affine = [ [R * diag([i j k]);[0 0 0]] [T;1] ];
         R_sort = sort(abs(R(:)));
         R( find( abs(R) < tolerance*min(R_sort(end-2:end)) ) ) = 0;
         R = R * diag([i j k]);
         hdr.hist.new_affine = [ [R;[0 0 0]] [T;1] ];

         if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
            msg = [char(10) char(10) '   Non-orthogonal rotation or shearing '];
            msg = [msg 'found inside the affine matrix' char(10)];
            msg = [msg '   in this NIfTI file. You have 3 options:' char(10) char(10)];
            msg = [msg '   1. Using included ''reslice_nii.m'' program to reslice the NIfTI' char(10)];
            msg = [msg '      file. I strongly recommand this, because it will not cause' char(10)];
            msg = [msg '      negative effect, as long as you remember not to do slice' char(10)];
            msg = [msg '      time correction after using ''reslice_nii.m''.' char(10) char(10)];
            msg = [msg '   2. Using included ''load_untouch_nii.m'' program to load image' char(10)];
            msg = [msg '      without applying any affine geometric transformation or' char(10)];
            msg = [msg '      voxel intensity scaling. This is only for people who want' char(10)];
            msg = [msg '      to do some image processing regardless of image orientation' char(10)];
            msg = [msg '      and to save data back with the same NIfTI header.' char(10) char(10)];
            msg = [msg '   3. Increasing the tolerance to allow more distortion in loaded' char(10)];
            msg = [msg '      image, but I don''t suggest this.' char(10) char(10)];
            msg = [msg '   To get help, please type:' char(10) char(10) '   help reslice_nii.m' char(10)];
            msg = [msg '   help load_untouch_nii.m' char(10) '   help load_nii.m'];
            error(msg);
         end

      else
         R = R * diag([i j k]);
      end					% 1st det(R)

   else
      affine_transform = 0;	% no sform or qform transform
   end

   if affine_transform == 1
      voxel_size = abs(sum(R,1));
      inv_R = inv(R);
      originator = inv_R*(-T)+1;
      orient = get_orient(inv_R);

      %  modify pixdim and originator
      %
      hdr.dime.pixdim(2:4) = voxel_size;
      hdr.hist.originator(1:3) = originator;

      %  set sform or qform to non-use, because they have been
      %  applied in xform_nii
      %
      hdr.hist.qform_code = 0;
      hdr.hist.sform_code = 0;
   end

   %  apply space_unit to pixdim if not 1 (mm)
   %
   space_unit = get_units(hdr);

   if space_unit ~= 1
      hdr.dime.pixdim(2:4) = hdr.dime.pixdim(2:4) * space_unit;

      %  set space_unit of xyzt_units to millimeter, because
      %  voxel_size has been re-scaled
      %
      hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,1,0));
      hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,2,1));
      hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,3,0));
   end

   hdr.dime.pixdim = abs(hdr.dime.pixdim);

   return;					% change_hdr

%-----------------------------------------------------------------------
function orient = get_orient(R)

   orient = [];

   for i = 1:3
      switch find(R(i,:)) * sign(sum(R(i,:)))
      case 1
         orient = [orient 1];		% Left to Right
      case 2
         orient = [orient 2];		% Posterior to Anterior
      case 3
         orient = [orient 3];		% Inferior to Superior
      case -1
         orient = [orient 4];		% Right to Left
      case -2
         orient = [orient 5];		% Anterior to Posterior
      case -3
         orient = [orient 6];		% Superior to Inferior
      end
   end

   return;					% get_orient

%-----------------------------------------------------------------------
function [space_unit, time_unit] = get_units(hdr)

   switch bitand(hdr.dime.xyzt_units, 7)	% mask with 0x07
   case 1
      space_unit = 1e+3;		% meter, m
   case 3
      space_unit = 1e-3;		% micrometer, um
   otherwise
      space_unit = 1;			% millimeter, mm
   end

   switch bitand(hdr.dime.xyzt_units, 56)	% mask with 0x38
   case 16
      time_unit = 1e-3;			% millisecond, ms
   case 24
      time_unit = 1e-6;			% microsecond, us
   otherwise
      time_unit = 1;			% second, s
   end

   return;					% get_units

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA