text:
function scalar2d = slice_scalar3d(scalar3d, normal_axis, intercept)

chkarg(istypesizeof(scalar3d, 'Scalar3d'), '"scalar3d" should be instance of Scalar3d.');
chkarg(istypesizeof(normal_axis, 'Axis'), '"normal_axis" should be instance of Axis.');
chkarg(istypesizeof(intercept, 'real'), '"intercept" should be real.');

grid3d = scalar3d.grid3d;
[h, v, n] = cycle(normal_axis);
chkarg(grid3d.comp(n).contains(intercept), ...
	'"intercept" should be contained in %s-axis range.', char(n));

lall = grid3d.lall(Axis.elems + Axis.count*subsindex(scalar3d.gt_array));
nall = lall{n};

ind = {':', ':', ':'};
indn = ismembc2(intercept, nall);
if indn ~= 0
	ind{n} = indn;
	array = scalar3d.array(ind{:});
else
% 	warning('Maxwell:interp', 'slice at %s = %s is not grid plane where %s is defined; fields are interpolated.', ...
% 		char(n), num2str(intercept), scalar3d.name);

	indn = find(nall < intercept, 1, 'last');
	ind{n} = [indn indn+1];

	V = scalar3d.array(ind{:});
	l = lall;
	l{n} = nall(ind{n});
	[Xh, Yv, Zn] = ndgrid(l{:});

	li = lall;
	li{n} = intercept;
	[XIh, YIv, ZIn] = ndgrid(li{:});

	array = interpn(Xh, Yv, Zn, V, XIh, YIv, ZIn);
end

array = permute(array, int([h v n]));
assert(ndims(array) == Dir.count);

grid2d = Grid2d(grid3d, n);
scalar2d = Scalar2d(array, grid2d, scalar3d.gt_array([h v]), scalar3d.osc, scalar3d.physQcell, scalar3d.name, intercept);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%
% Author: Marius Drulea
% http://www.cv.utcluj.ro/optical-flow.html

% References
% M. Drulea and S. Nedevschi, "Total variation regularization of
% local-global optical flow," in Intelligent Transportation Systems (ITSC),
% 2011 14th International IEEE Conference on, 2011, pp. 318-323.

% Copyright (C) 2011 Technical University of Cluj-Napoca

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

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

function [u, v, pu, pv] = solve_clg_tv_equation(u, v, u0, v0, pu, pv, ...
    I1, Ix, Iy, It, D, settings)
%%
lambda = settings.lambda;

theta = 0.5;
tau = 0.5;
% tau = 1.0/(4.0*theta + epsilon);

der_mask = [0 -1 1];
adjoint_der_mask = [-1 1 0];

if settings.use_bilateral == 1
     r0 = It - u0.*Ix - v0.*Iy;
    % apply the bilateral filter to the data terms
    [W_Ix_2, W_Ixy, W_Iy_2, W_Ix_r0, W_Iy_r0] = ...
    applyBilateralFilterToDataTerms(I1, Ix.*Ix, Ix.*Iy, Iy.*Iy, Ix.*r0, Iy.*r0, ...
        settings.wSize, settings.sigma_d, settings.sigma_r);
else
    gauss = fspecial('gaussian', [settings.wSize settings.wSize], settings.wSize/6);

    W_Ix_2 = imfilter(Ix.^2, gauss, 'replicate');
    W_Ixy = imfilter(Ix.*Iy, gauss, 'replicate');
    W_Iy_2 = imfilter(Iy.^2, gauss, 'replicate');

    r0 = It - u0.*Ix - v0.*Iy;
    W_Ix_r0 = imfilter(Ix.*r0, gauss, 'replicate');
    W_Iy_r0 = imfilter(Iy.*r0, gauss, 'replicate');
end

a11 = 1 + 2*lambda*theta*W_Ix_2;
a12 = 2*lambda*theta*W_Ixy;
a21 = a12;
a22 = 1 + 2*lambda*theta*W_Iy_2;
l_t_2_W_Ix_r0 = 2*lambda*theta*W_Ix_r0;
l_t_2_W_Iy_r0 = 2*lambda*theta*W_Iy_r0;

delta = a11.*a22 - a21.*a12;
%%
for k = 1:settings.its
    %%
    % 1. update the coupling variable

    b1 = u - l_t_2_W_Ix_r0;
    b2 = v - l_t_2_W_Iy_r0;

    % update u_ and v_
    deltaU_ = b1.*a22 - b2.*a12;
    deltaV_ = b2.*a11 - b1.*a21;

    u_ = deltaU_./delta;
    v_ = deltaV_./delta;

    %%
  % compute the divergence of the dual variable
  % the adjoint of the nabla (derivative) operator = (- divergence) operator

  div_u = imfilter(pu(:, :, 1), adjoint_der_mask, 'replicate') + ...
          imfilter(pu(:, :, 2), adjoint_der_mask', 'replicate');
  div_v = imfilter(pv(:, :, 1), adjoint_der_mask, 'replicate') + ...
          imfilter(pv(:, :, 2), adjoint_der_mask', 'replicate');

  % update primal variable u
  u = u_ + theta*div_u;

  % update primal variable v
  v = v_ + theta*div_v;

  % compute nabla(u); the derivative operator
  ux = imfilter(u, der_mask, 'replicate');
  uy = imfilter(u, der_mask', 'replicate');
  vx = imfilter(v, der_mask, 'replicate');
  vy = imfilter(v, der_mask', 'replicate');

  % update dual variable; gradient descent
  % p = p_k + tau*nabla(u);
  pu(:, :, 1) = pu(:, :, 1) + tau * ux;
  pu(:, :, 2) = pu(:, :, 2) + tau * uy;
  pv(:, :, 1) = pv(:, :, 1) + tau * vx;
  pv(:, :, 2) = pv(:, :, 2) + tau * vy;

  % project the dual variable to ensure the inequality |p| <= D
  % p = p./max(|p|, D) .* D;
  pu(:, :, 1) = pu(:, :, 1)./max(abs(pu(:, :, 1)), D) .* D;
  pu(:, :, 2) = pu(:, :, 2)./max(abs(pu(:, :, 2)), D) .* D;
  pv(:, :, 1) = pv(:, :, 1)./max(abs(pv(:, :, 1)), D) .* D;
  pv(:, :, 2) = pv(:, :, 2)./max(abs(pv(:, :, 2)), D) .* D;
end

author: thomas-koehler

repo: SupER


text:
function [DMIN,DMIN_aerr] = dicomrt_MINcal(inputdose,inputerror,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use)
% dicomrt_MINcal(inputdose,inputerror,dose_xmesh,dose_ymesh,dose_zmesh,VOI,voi2use)
%
% Calculate the MIN Dose for a given 3D dose distribution within a given VOI
%
% inputdose is the input 3D dose (e.g. RTPLAN or MC generated)
% inputerror is the relative error asociated with the calculated dose. If inputerror=0 error
% calculation is not performed (e.g. for TPS dose matrices), otherwise
% inputerror dimensions must match those of inputdose (e.g. for MC dose matrices).
% dose_xmesh,dose_ymesh,dose_zmesh are x-y-z coordinates of the center of the dose-pixel
%
% VOI is a cell array which contain the patients VOIs as read by dicomrt_loadvoi
% voi2use is a vector pointing to the number of VOIs to be used ot the analysis and for the display.
%
% Example:
%
% [A]=dicomrt_MINcal(B,0,dose_xmesh,dose_ymesh,dose_zmesh,demo_voi,9);
%
% returns in A the min dose for the dose matrix B and the voi # 9 within demo_voi.
%
% NOTE:
% due to the algorithm that MATLAB uses to filter matrices (roipoly), it can happen that
% some of the voxels which are partially outside the patient are included in the VOI.
% This can happen when the border of the VOI are too close to the patient outline with
% respect to the size of the calculation matrix used.
% This is not a problem for MC which calculate dose also outside the patient outline.
% It may represent a problem when dealing with TPS generated data, which do not calculate
% dose outside the patient outline, because voxels with zero dose can be counted in the target
% volume. This event is rare (1/10000) and should not influence dose analysis.
% A simple way to overcome this problem is to define VOIs as separate as possible.
% In this function the minimum dose ~=0 is reported.
%
% See also dicomrt_MDNcal, dicomrt_MDcal, dicomrt_MAXcal
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check case and set-up some parameters and variables
[dose_temp,derror_temp,type_dose,doselabel]=dicomrt_checkinputanderr(inputdose,inputerror,1);
dose=dicomrt_varfilter(dose_temp);
derror=dicomrt_varfilter(derror_temp);
[VOI_temp,type,label]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(VOI_temp);

% Initialise parameters
count=0;

[mask_dose,volume_VOI,mask4VOI]=dicomrt_mask(VOI_temp,dose_temp,dose_xmesh,dose_ymesh,dose_zmesh,voi2use,'nan','n');
mask=dicomrt_varfilter(mask_dose);

if size(derror)~=1
    [mask_error,volume_VOI,mask4VOI]=dicomrt_mask(VOI_temp,derror_temp,dose_xmesh,dose_ymesh,dose_zmesh,voi2use,'nan','n');
    msk_error=dicomrt_varfilter(mask_error);
end

for k=1:size(mask,3)
    [ilocation,jlocation]=find(isnan(mask(:,:,k))~=1);
    if isnumeric(ilocation)==1 | isnumeric(jlocation)==1
        for i=1:length(ilocation)
            count=count+1;
            dose_VOI(count)=mask(ilocation(i),jlocation(i),k);
            if size(derror)~=1
                % this is absolute error now
                error_VOI(count)=msk_error(ilocation(i),jlocation(i),k).*mask(ilocation(i),jlocation(i),k);
            else
                error_VOI=0;
            end
        end
    end
end

if size(derror)~=1
    % Mean is calculated as weighted average of data points
    [DMIN,I]=min(dose_VOI);
    % Dose min cannot be 0.
    if DMIN==0
        dose_VOI(I)=[];
        [DMIN,I]=min(dose_VOI);
    end
    % error is calulated as standard error of the mean
    DMIN_aerr=error_VOI(I);
else
    % Mean is calculated as average of data points
    [DMIN,I]=min(dose_VOI);
    % Dose min cannot be 0.
    if DMIN==0
        dose_VOI(I)=[];
        [DMIN,I]=min(dose_VOI);
    end
    % error is calulated as standard error of the mean
    DMIN_aerr=0;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [pro, res] = transferRT0toCR3(node,elem,fixedface)
% Reference: https://arxiv.org/abs/2102.03396
%
%
% See also transferRT0toCR

% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Data Structure and Basis for RT0
[elem2face,face,dofSign] = dof3face(elem);
Dlambda = gradbasis3(node,elem);
NF = size(face,1);
locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];

%% Assembling
T = sparse(3*NF,NF);
for i = 1:4
    for j = 1:4
        ii = double(elem2face(:,i));
		jj = double(elem2face(:,j));
        j1 = locFace(j,1); j2 = locFace(j,2);  j3 = locFace(j,3);
        % [j1,j2] is the edge opposite to vertex j.
%         sT = (1/2*(i~=j1)*Clambda(:,:,j2) - 1/2*(i~=j2)*Clambda(:,:,j1))...
%              .*repmat(double(dofSign(:,j)),1,2);      % phi_j evaluates at m_i.
        sT = -(Dlambda(:,:,j1) + Dlambda(:,:,j2) + Dlambda(:,:,j3))...
             .*repmat(double(dofSign(:,j)),1,3)/3;      % phi_j evaluates at c_t.
        sT = 1/2*sT;
        T = T + sparse([ii NF+ii 2*NF+ii], [jj jj jj], sT(:), 3*NF, NF);
    end
end
pro = T;
pro([fixedface;fixedface;fixedface],:) = 2*pro([fixedface;fixedface;fixedface],:);
res = pro';
end

repo: ifem


text:
classdef ChTargetFactory < handle

    methods (Access = public, Static)

        function Ch_star = create(cParams)

            switch cParams.type
                case 'negative_poisson'

                    E1  = cParams.E_plus;
                    E0  = cParams.E_minus;
                    nu1 = cParams.nu_plus;
                    nu0 = cParams.nu_minus;
                    kappa_f = @(E,nu) E/2*(1-nu);
                    mu_f = @(E,nu) E/2*(1-nu);

                    k_plus = kappa_f(E1,nu1);
                    mu_plus = mu_f(E1,nu1);

                    k_minus = kappa_f(E0,nu0);
                    mu_minus = mu_f(E0,nu0);

                    nu = @(k,mu) (k-mu)/(k+mu);
                    E = @(k,mu) (4*k*mu)/(k+mu);
                    C = @(E,nu) (E/(1-nu*nu)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2]);

                    kappa_nu_min = k_minus;
                    mu_nu_min = mu_plus;

                    nu_min = nu(kappa_nu_min,mu_nu_min);
                    E_nu_min = E(kappa_nu_min,mu_nu_min);
                    C_nu_min = C(E_nu_min,nu_min);
                    Ch_star = C_nu_min;

                case 'nu_0_6' %From Sigmund Thesis% rho = 0.38
                    C = @(E,nu) (E/(1-nu*nu)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2]);
                    nu = -0.6;
                    E = (1-nu*nu)*0.04;
                    Ch_star = C(E,nu);

                case 'Seba' % Es=0.08; nus=-0.25
                    Ch_star = [0.0853    -0.0213       0;
                               -0.0213    0.0853       0;
                                   0         0    0.0533];

                case 'Nu0_2'
                    %Ch_star = [0.0287   -0.0069    0.0068;
                    %           -0.0069    0.0287    0.0068;
                    %           0.0068    0.0068    0.0137];

                    %Ch_star = [0.0129   -0.0158    0.0000;
                    %           -0.0158    0.0812   -0.0000;
                    %           0.0000   -0.0000    0.0021];

                    Ch_star = [  0.0157   -0.0174   -0.0000;
                                -0.0174    0.0815   -0.0000;
                                -0.0000   -0.0000    0.0025];

                case 'nu_0_8' %From Sigmund Thesis% rho = 0.25
                    C = @(E,nu) (E/(1-nu*nu)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2]);
                    nu = -0.8;
                    E = (1-nu*nu)*0.02;
                    Ch_star = C(E,nu);

                case 'Vfrac07'
                    Ch_star =[
                        0.4256    0.2837         0
                        0.2837    0.4256         0
                        0         0    0.1419];

                case 'Vfrac06'
                    Ch_star =[
                        0.2909    0.1940         0
                        0.1940    0.2909         0
                        0         0    0.0970];

                case 'Vfrac05'
                    Ch_star =[
                        0.1892    0.1261         0
                        0.1261    0.1892         0
                        0         0    0.0631];

                case 'Vfrac04'
                    Ch_star =[
                        0.1141    0.0761         0
                        0.0761    0.1141         0
                        0         0    0.0380];

                case 'Vfrac04b' %Circle of 0.4 radius
                         Ch_star = [ 0.7519    0.2237   -0.0000;
                                  0.2237    0.7519   -0.0000;
                                 -0.0000   -0.0000    0.2245];

                case 'Vfrac03'
                    Ch_star =[
                        0.0611    0.0407         0
                        0.0407    0.0611         0
                        0         0    0.0204];

                case 'HorizontalRectangleInclusion'
                    Ch_star = [0.4637    0.0010    0.0000
                               0.0010    0.0031    0.0000
                               -0.0000    0.0000    0.0011];

                case 'Composite'
                    Ch_star =[1 0.15 0;
                        0.15 0.5 0;
                        0 0 0.2];
                case 'HoneyComb'
                    Ch_star =0.094*[1 0.75 0
                        0.75 1 0
                        0 0 0.125];
                case 'InvertedHoneyComb'
                    Ch_star =0.08*[1 -0.5 0
                        -0.5 1 0
                        0 0 0.06];
                case 'AcousticZeroShearA'
                    Ch_star =[1 1 0;
                        1 1 0;
                        0 0 0];
                case 'NegativePoiss06'
                    Ch_star =0.04*[1 -0.6 0;
                        -0.6 1 0;
                        0 0 0.8];
                case 'IsotropyHexagon'
                    Ch_star = [0.5931 0.1882 0.0000;
                               0.1882 0.5931 0.0000;
                               0.0000 0.0000 0.2025];
            end
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% SPEC_FREQ_AVERAGE Calulates the frequency-averaged spectrogram
%
% Usage
%    spec = SPEC_FREQ_AVERAGE(x, filters, options)
%
% Input
%    x (numeric): The signal to be transformed.
%    filters (cell): The set of filters used to calculate the spectrogram
%       and the frequency averaging.
%    options (struct): Various options to the function:
%       options.oversampling (int): The amount of oversampling with respect to
%          the critical bandwidth, as a power of 2 (default 1).
%
% Output
%    spec (numeric): The frequency-averaged spectrogram, arranged in as a
%       first-order scattering transform.
%
% Description
%    First, the spectrogram of x is computed with respect to the second-order
%    lowpass filter filters{2}.phi, if it is present, otherwise filters{1}.phi
%    is used. Once a spectrogram is obtained, the filters in filters{1}.psi are
%    averaged with the spectrogram to each yield a time-varying coefficient.
%    The resulting signals are arranged in spec in the format of a first-order
%    scattering transform. As a result, spec{1}.signal is empty, since no
%    zeroth-order coefficients are present. However, spec{2}.signal contains
%    the frequency-averaged spectrogram coefficients, arranged in order of
%    decreasing frequeny.
%
%    The accompanying meta structure is also analogous to the scattering
%    transform. Notably, spec{2}.meta.order is equal to 1 for each coefficient
%    and spec{2}.meta.j gives the index of the filter used to compute the
%    frequency averaging.
%
% See also
%   MORLET_FILTER_BANK_1D, WAVELET_FACTORY_1D

function [out,meta] = spec_freq_average(in,filters,options)
	options = fill_struct(options,'oversampling',1);

	filters1 = filters{1};
	filters2 = filters{min(2,numel(filters))};

	[temp1,temp2,phi2_bw] = filter_freq(filters{2}.meta);

	supp_mult = 4;

	N = size(in,1);
	sig_count = size(in,3);
	Nfilt = filters1.meta.size_filter;
	N1 = 2^round(log2(2*pi/phi2_bw));

	fs = zeros(N1*supp_mult,length(filters1.psi.filter));

	for k1 = 0:length(filters1.psi.filter)-1
		f_temp = realize_filter(filters1.psi.filter{k1+1});
		fs(:,k1+1) = abs(f_temp(1:Nfilt/(N1*supp_mult):Nfilt));
	end

	window = ifft(realize_filter(filters2.phi.filter));
	window = [window(Nfilt-N1*supp_mult/2+1:Nfilt); window(1:N1*supp_mult/2)];

	frames = zeros(length(window),round(N/N1*2^options.oversampling),sig_count);

	out = zeros(size(frames,2),size(fs,2),sig_count);

	for t = 0:size(out,1)-1
		ind = round(t*N1/2^options.oversampling+[-N1*supp_mult/2:N1*supp_mult/2-1]+1);

		% symmetric extension
		ind(ind<1) = 1-ind(ind<1);
		ind(ind>N) = 2*N+1-ind(ind>N);

		frames(:,t+1,:) = in(ind,:);
	end

	frames = bsxfun(@times,frames,window);

	frame_fm = abs(fft(frames,[],1));

	out = fs.'*reshape(frame_fm,[size(frame_fm,1) size(frame_fm,2)*size(frame_fm,3)]);
	out = reshape(out,[size(fs,2) size(frame_fm,2) size(frame_fm,3)]);

	resolution = round(log2(N1))-options.oversampling;

	meta.order = ones(1,size(fs,2));
	meta.scale = [0:size(fs,2)-1];
	meta.bandwidth = phi2_bw*ones(1,size(fs,2));
	meta.resolution = resolution*ones(1,size(fs,2));

	X = cell(1,2);
	X{1}.signal = {};
	X{1}.meta.bandwidth = zeros(1,0);
	X{1}.meta.scale = zeros(1,0);

	X{2}.signal = cell(1,size(out,1));
	X{2}.meta = meta;

	for k0 = 0:size(out,1)-1
		X{2}.signal{k0+1} = reshape(out(k0+1,:,:),[size(out,2) 1 size(out,3)]);
	end

	out = X;
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [R, singular_values, determinants] = round_solution(Yopt, problem_data)
%function [R, singular_values, determinants] = round_solution(Yopt, problem_data)
%
% Given an element Yopt in St(d, r)^n, this function rounds Yopt to an
% element of SO(d)^n

% Copyright (C) 2016 by David M. Rosen

r = size(Yopt, 1);

[U, Xi, V] = svd(Yopt, 'econ');
singular_values = diag(Xi)';

Xi_d = Xi(1:problem_data.d, 1:problem_data.d);  %Xi_d is the upper-left dxd submatrix of Xi
V_d = V(:, 1:problem_data.d);  %V_d contains the first d columns of V

R = Xi_d*V_d';

determinants = zeros(1, problem_data.n);

for k = 1:problem_data.n
    determinants(k) = det(R(:, problem_data.d*(k-1) + 1 : problem_data.d*(k-1) + problem_data.d));
end
ng0 = sum(determinants > 0);

reflector = diag([ones(1, problem_data.d - 1), -1]);  % Orthogonal matrix that we can use for reversing the orientations of the orthogonal matrix subblocks of R

if ng0 == 0
    % This solution converged to a reflection of the correct solution
    R = reflector*R;
    determinants = -determinants;
elseif ng0 < problem_data.n
    disp('WARNING: SOLUTION HAS INCONSISTENT ORIENTATIONS!');

    % If more than half of the determinants have negative sign, reverse
    % them
    if ng0 < problem_data.n / 2
        determinants = -determinants;
        R = reflector * R;
    end
end

% Finally, project each element of R to SO(d)
for i = 1:problem_data.n
    R(:, problem_data.d * (i-1) + 1 : problem_data.d *i) = project_to_SOd( R(:, problem_data.d * (i-1) + 1 : problem_data.d *i) );
end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function IR = readIR(data,frameID)

IR = imread(data.ir{frameID});

IR = IR(:,end:-1:1,:);

author: jianxiongxiao

repo: ProfXkit


text:
function A = four_corners(m,n)
%FOUR_CORNERS Construct a m x n sparse matrix with ones in each of its four
%corners.
%
% t = four_corners(m,n);
%
% Inputs:
%  m,n  the dimension of the matrix requested
% Outputs:
%  A  the sparse matrix with ones in their four corners
%

A = ...

end

author: odedstein

repo: sgi-introduction-course


text:
function mask = PTKGetLungROIFromLeftAndRightLungs(left_and_right_lung_mask, context, reporting)
    % PTKGetLungROIFromLeftAndRightLungs. Extracts a region of interest for the left or right lung given the original image and the left and right lung mask.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    mask = left_and_right_lung_mask.BlankCopy;

    if context == PTKContext.RightLung
        colormap_index = PTKColormapLabels.RightLung;
    elseif context == PTKContext.LeftLung
        colormap_index = PTKColormapLabels.LeftLung;
    else
        reporting.Error('PTKGetLungROIFromLeftAndRightLungs:UnknownContext', ['The context ' char(context) ' is unknown.']);
    end
    mask.ChangeRawImage(left_and_right_lung_mask.RawImage == colormap_index);
    mask.CropToFitWithBorder(5);
end

author: tomdoel

repo: pulmonarytoolkit


text:
clc;

%This code runs a simulation of a single spring pendulum, running all
%dynamics and integration in C++ code. Matlab is only used for plotting and
%compiling the code (to a mex function).

%First, compile the most recent version of the file
Recompile_Code = true;
if Recompile_Code
    mex Cpp_Integrator.cpp;
else
    disp('Warning - Did not recompile C++ files!')
    disp(' ')
end

% Initial state for the system
% [x,y,dx,dy]';
x0 = [  0;
        -0.4;
        0.6;
        0];   %Must be a column vector
Tspan = [0,20];
N = 1e5;

Params = zeros(1,6);  %Must be a row vector!
Params(1) = 0.8;    %(kg)   Mass
Params(2) = 89;    %(N/m)  Spring Constant
Params(3) = 0;    %(Ns/m)  Damping
Params(4) = 9.81;   %(m/s^2)  Gravity
Params(5) = 0.2;    %(m) Rest length of spring
Params(6) = 0.5;    %(Nm)    Forcing Amplitude
Params(7) = 20;      %(Hz)   Forcing frequency (Hz)

tic
[t, X] = Cpp_Integrator(Tspan,x0,N,Params);
disp(' ')
disp(['C++ Time: (Normal) ' num2str(toc)]);

x = X(1,:);
y = X(2,:);
dx = X(3,:);
dy = X(4,:);

figure(1)
plot(t,X)
title('All States')

figure(2)
plot(x,y)
title('Position Trace')

author: MatthewPeterKelly

repo: dscTutorials


text:
function volume = ellipsoid_volume ( m, a, v, r )

%*****************************************************************************80
%
%% ELLIPSOID_VOLUME returns the volume of an ellipsoid.
%
%  Discussion:
%
%    The points X in the ellipsoid are described by an M by M
%    positive definite symmetric matrix A, an M-dimensional point V,
%    and a "radius" R, such that
%      (X-V)' * A * (X-V) <= R * R
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, real A(M,M), the matrix that describes
%    the ellipsoid.  A must be symmetric and positive definite.
%
%    Input, real V(M), the "center" of the ellipse.
%    The value of V is not actually needed by this function.
%
%    Input, real R, the "radius" of the ellipse.
%
%    Output, real VOLUME, the volume of the ellipsoid.
%
  [ u, info ] = r8po_fa ( m, a );

  sqrt_det = 1.0;
  for i = 1 : m
    sqrt_det = sqrt_det * u(i,i);
  end

  volume = r ^ m * hypersphere_unit_volume ( m ) / sqrt_det;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Y = spm_ho_poly(P,M,U,varargin)
% General polynomial mapping with derivatives
% FORMAT Y = spm_ho_poly(P,M,U)
%
% P    - polynomial parameters (P{i} = i-th order coefficients)
% M    - model structure
% U    - (m,n) inputs
%
% Y(i) =  P{1} + P{2}*U(:,i) + P{3}*kron(U(:,i),U(:,i)) + ...
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_ho_poly.m 5709 2013-10-22 11:07:29Z guillaume $

% evaluate
%--------------------------------------------------------------------------
nu  = size(U,2);
if nargin > 3

    % evaluate
    %----------------------------------------------------------------------
    Y     = sparse(1,nu) + P{1};
    for i = 1:nu
        X     = 1;
        for j = 2:length(P)
            X      = kron(X,U(:,i));
            Y(1,i) = Y(1,i) + spm_vec(P{j})'*X;
        end
    end

else

    % evaluate with derivatives dFdu
    %----------------------------------------------------------------------
    for i = 1:nu
        [dFdu,F] = spm_diff('spm_ho_poly',P,M,U(:,i),'no diff',3);
        Y(:,i)   = spm_vec(F,dFdu);
    end

end

% place samples in rows
%----------------------------------------------------------------------
Y  = Y';

author: spm


text:
function [dP_dV, dP_dlam] = cpf_p_jac(parameterization, z, V, lam, Vprv, lamprv, pv, pq)
%CPF_P_JAC Computes partial derivatives of CPF parameterization function.
%   [DP_DV, DP_DLAM ] = CPF_P_JAC(PARAMETERIZATION, Z, V, LAM, ...
%                                                   VPRV, LAMPRV, PV, PQ)
%
%   Computes the partial derivatives of the continuation power flow
%   parameterization function w.r.t. bus voltages and the continuation
%   parameter lambda.
%
%   Inputs:
%       PARAMETERIZATION : Value of cpf.parameterization option.
%       Z : normalized tangent prediction vector from previous step
%       V : complex bus voltage vector at current solution
%       LAM : scalar lambda value at current solution
%       VPRV : complex bus voltage vector at previous solution
%       LAMPRV : scalar lambda value at previous solution
%       PV : vector of indices of PV buses
%       PQ : vector of indices of PQ buses
%
%   Outputs:
%       DP_DV : partial of parameterization function w.r.t. voltages
%       DP_DLAM : partial of parameterization function w.r.t. lambda
%
%   See also CPF_PREDICTOR, CPF_CORRECTOR.

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

if parameterization == 1        %% natural
    npv = length(pv);
    npq = length(pq);
    dP_dV = zeros(1, npv+2*npq);
    if lam >= lamprv
        dP_dlam = 1.0;
    else
        dP_dlam = -1.0;
    end
elseif parameterization == 2    %% arc length
    Va = angle(V);
    Vm = abs(V);
    Vaprv = angle(Vprv);
    Vmprv = abs(Vprv);
    dP_dV = 2*([Va([pv; pq]); Vm(pq)] - [Vaprv([pv; pq]); Vmprv(pq)])';
    if lam == lamprv    %% first step
        dP_dlam = 1.0;  %% avoid singular Jacobian that would result
                        %% from [dP_dV, dP_dlam] = 0
    else
        dP_dlam = 2*(lam-lamprv);
    end
elseif parameterization == 3    %% pseudo arc length
    nb = length(V);
    dP_dV = z([pv; pq; nb+pq])';
    dP_dlam = z(2*nb+1);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function submit()
  addpath('./lib');

  conf.assignmentSlug = 'linear-regression';
  conf.itemName = 'Linear Regression with Multiple Variables';
  conf.partArrays = { ...
    { ...
      '1', ...
      { 'warmUpExercise.m' }, ...
      'Warm-up Exercise', ...
    }, ...
    { ...
      '2', ...
      { 'computeCost.m' }, ...
      'Computing Cost (for One Variable)', ...
    }, ...
    { ...
      '3', ...
      { 'gradientDescent.m' }, ...
      'Gradient Descent (for Multiple Variables)', ...
    }, ...
    { ...
      '4', ...
      { 'featureNormalize.m' }, ...
      'Feature Normalization', ...
    }, ...
    { ...
      '5', ...
      { 'computeCostMulti.m' }, ...
      'Computing Cost (for Multiple Variables)', ...
    }, ...
    { ...
      '6', ...
      { 'gradientDescentMulti.m' }, ...
      'Gradient Descent (for Multiple Variables)', ...
    }, ...
    { ...
      '7', ...
      { 'normalEqn.m' }, ...
      'Normal Equations', ...
    }, ...
  };
  conf.output = @output;

  submitWithConfiguration(conf);
end

function out = output(partId)
  % Random Test Cases
  X1 = [ones(20,1) (exp(1) + exp(2) * (0.1:0.1:2))'];
  Y1 = X1(:,2) + sin(X1(:,1)) + cos(X1(:,2));
  X2 = [X1 X1(:,2).^0.5 X1(:,2).^0.25];
  Y2 = Y1.^0.5 + Y1;
  if partId == '1'
    out = sprintf('%0.5f ', warmUpExercise());
  elseif partId == '2'
    out = sprintf('%0.5f ', computeCost(X1, Y1, [0.5 -0.5]'));
  elseif partId == '3'
    out = sprintf('%0.5f ', gradientDescent(X1, Y1, [0.5 -0.5]', 0.01, 10));
  elseif partId == '4'
    out = sprintf('%0.5f ', featureNormalize(X2(:,2:4)));
  elseif partId == '5'
    out = sprintf('%0.5f ', computeCostMulti(X2, Y2, [0.1 0.2 0.3 0.4]'));
  elseif partId == '6'
    out = sprintf('%0.5f ', gradientDescentMulti(X2, Y2, [-0.1 -0.2 -0.3 -0.4]', 0.01, 10));
  elseif partId == '7'
    out = sprintf('%0.5f ', normalEqn(X2, Y2));
  end
end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function [X, sigma2, W] = ppcaEmbed(Y, dims)

% PPCAEMBED Embed data set with probabilistic PCA.
% FORMAT
% DESC returns latent positions for a given data set via probabilistic
% PCA.
% ARG Y : the data set which you want the latent positions for.
% ARG dims : the dimensionality of the latent space.
% RETURN X : the latent positions.
% RETURN sigma2 : the variance not explained by the latent positions.
% RETURN W : the matrix required to invert the transformation, Y=X*W'.
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% SEEALSO : lleEmbed, isomapEmbed

% MLTOOLS

if ~any(any(isnan(Y)))
  if size(Y, 1)<size(Y, 2)
    Ymean = mean(Y);
    Ycentre = zeros(size(Y));
    for i = 1:size(Y, 1)
      Ycentre(i, :) = Y(i, :) -Ymean;
    end
    if size(Ycentre, 2)>30000
      % Bug in MATLAB 7.0 means you have to do this.
      innerY = zeros(size(Ycentre, 1));
      for i = 1:size(Ycentre, 1)
        innerY(i, :) = Ycentre(i, :)*Ycentre';
      end
    else
      innerY = Ycentre*Ycentre';
    end
    [v, u] = eigdec(innerY, dims);
    v(find(v<0))=0;
    X = u(:, 1:dims)*sqrt(size(Y, 1));
    v = v/sqrt(size(Y, 1));
    sigma2 = (trace(innerY) - sum(v))/(size(Y, 2)-dims);
    W = X'*Ycentre;

  else
    [v, u] = pca(Y);
    v(find(v<0))=0;
    Ymean = mean(Y);
    Ycentre = zeros(size(Y));
    for i = 1:size(Y, 2);
      Ycentre(:, i) = Y(:, i) - Ymean(i);
    end
    X = Ycentre*u(:, 1:dims)*diag(1./sqrt(v(1:dims)));
    sigma2 = mean(v(dims+1:end));
    W = X'*Ycentre;
  end
else
  % Hacky implementation of Probabilistic PCA for when there is missing data.
  iters = 100;
  % Initialise W randomly
  d = size(Y, 2);
  q = dims;
  N = size(Y, 1);
  W = randn(d, q)*1e-3;
  sigma2 = 1;
  mu = zeros(d, 1);
  for i = 1:d
    obs = find(~isnan(Y(:, i)));
    if length(obs)>0
      mu(i) = mean(Y(obs, i));
    else
      mu(i) = 0;
    end
  end
  numObs = sum(sum(~isnan(Y)));
  for i = 1:iters
    M = W'*W + sigma2*eye(q);
    invM = inv(M);
    exp_xxT = zeros(q);
    exp_x = zeros(N, q);
    for n = 1:N
      obs = find(~isnan(Y(n, :)));
      exp_x(n, :) = (invM*W(obs, :)'*(Y(n, obs)' - mu(obs)))';
    end
    exp_xxT = N*sigma2*invM + exp_x'*exp_x;
    s = zeros(d, q);
    s2 = 0;
    for n = 1:N
      obs = find(~isnan(Y(n, :)));
      subY = zeros(size(Y(n, :)))';
      subY(obs) = Y(n, obs)' - mu(obs);
      s = s + (subY)*exp_x(n, :);
      s2 = s2 + sum((Y(n, obs)' - mu(obs)).^2) - 2*exp_x(n, :)*W(obs, :)'*(Y(n, obs)' - mu(obs));
    end
    W = s*inv(exp_xxT);
    sigma2 = 1/(numObs)*(s2 + trace(exp_xxT*W'*W));
  end

  X = exp_x;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Earth Observing System Data Visualization
%
%   Part 1 (Main Example)
%
%   This M-file contains the sequence of MATLAB commands listed in the
%   November 2002 MATLAB Digest article on "Reading and Visualizing Data
%   from the Earth Observing System (EOS)," up to the section that
%   requires the Mapping Toolbox.
%
%   Rob Comer and Chris Lawton
%   Copyright 2002 The MathWorks, Inc.

%== DATA SELECTION AND HDFTOOL============================================

% Load the sea surface temperature (SST) swath data set, the corresponding
% geolocation data set, and the NDVI grid into HDFTOOL.
hdftool('MOD28L2.A2001185.0830.003.2001308112641.hdf')
hdftool('MOD03.A2001185.0830.003.2001305034045.hdf')
hdftool('PAL_CLIMATE_JUL_01-10_2001.HDF')

%== COLOR SCALES AND SHARING A COLORMAP ==================================

% We'll want separate colorbars for SST and NDVI, and we'll want tick marks
% labelled in physical units. So instead of using MTALAB's COLORBAR
% command, we'll use our COLORSCALE function. Because each figure can only
% have one colormap, we concatenate a 128-color SST colormap with a
% 128-color NDVI colormap.

% Load the colormap and display a standard colorbar.
load colormaps
cmap = [sstCmap; ndviCmap];   % We saved our colormap in two equal parts
figure
colormap(cmap)
set(axes,'Visible','off')
h = colorbar('horiz');
set(get(h,'Title'),'String','Standard Colorbar')

% Examine the colormap with COLORMAPEDITOR.
colormapeditor

% Assign colormap ranges to SST and NDVI, and relate them to physical
% units.
sstCmapLim  = [3 128];
sstDataLim  = [0  35];
ndviCmapLim = [  131   256];
ndviDataLim = [-0.25  1.00];

% Add color scales for SST and NDVI to the current figure.
pos = get(h,'Position'); % Position color scales relative to the color bar
colorscale(sstCmapLim, sstDataLim,5,'horiz','Position',pos + [0 0.70 0 0])
title('SST Colorscale')
xlabel('degrees Celsius')
colorscale(ndviCmapLim,ndviDataLim,0.2,'horiz','Position',pos + [0 0.35 0 0])
title('NDVI Colorscale')

%== READING MODIS SEA SURFACE TEMPERATURE SWATH ==========================

% Read MODIS Sea-Surface Temperature (SST).
mod28L2File = 'MOD28L2.A2001185.0830.003.2001308112641.hdf';
sstInfo = hdfinfo(mod28L2File,'eos');
sstSwath = sstInfo.Swath;
idx = {[1 1],[4 4],[]}; % Subset by 4 in each dimension
sst = hdfread(sstSwath,'Fields','sst','Index',idx);

% Read flags and construct land mask.
commonFlags = hdfread(sstSwath,'Fields','common_flags','Index',idx);
land = (bitget(commonFlags,8) == 1);

% Read SST attributes: scale and units.
sstInfoHDF = hdfinfo(mod28L2File);
datasets = sstInfoHDF.Vgroup.Vgroup(2).SDS;
sstDataset = datasets(strmatch('sst',{datasets.Name},'exact'));
sstAttr    = sstDataset.Attributes;
sstSlope     = double(sstAttr(strmatch('Slope',    {sstAttr.Name})).Value)
sstIntercept = double(sstAttr(strmatch('Intercept',{sstAttr.Name})).Value)
sstUnits     = sstAttr(strmatch('Units',{sstAttr.Name})).Value

% Read MODIS geolocation data fields.
mod03File = 'MOD03.A2001185.0830.003.2001305034045.hdf';
geoInfo = hdfinfo(mod03File,'eos');
sstLat = double(hdfread(geoInfo.Swath,'Fields','Latitude', 'Index',idx));
sstLon = double(hdfread(geoInfo.Swath,'Fields','Longitude','Index',idx));
sstLat(sstLat == -999) = NaN;
sstLon(sstLon == -999) = NaN;

%== MAPPING SST VALUES TO THE COLORMAP ===================================

% The SST values are encoded in the HDF-EOS data file.  To convert an encoded
% data value to absolute temperature units (degrees Celsius), one must
% multiply by sst_slope, then add sst_intercept. However, instead of
% converting to an intermediate double-valued array in degrees Celsius,
% we directly scale and shift the encoded SST values to correspond to the
% SST section of the colormap, in preparation for visualizing them via
% direct color mapping (i.e., with 'CDataMapping' set to 'direct'). First
% determine coefficients a(1) and a(2) such that a(1) * sst + a(2) will
% linearly map the encoded SST values directly onto the portion of the
% SST section of our colormap.

% Find the encoded SST limits that correspond to the temperature limits
% in sstDataLim.
sstCodelim = (sstDataLim - sstIntercept)/sstSlope;

% Solve a 2-by-2 linear system:
%   a(1) * sstCodelim + a(2) = sstCmapLim - 1
a = (sstCmapLim - 1)/[sstCodelim; [1 1]];

% Finally, scale and shift, clip to stay within the appropriate part
% of the colormap, and convert to uint8 (for economy of storage).
% Also, be sure to remap the special value of zero and apply the land mask.
sstMapped = a(1) * double(sst) + a(2);  % Scale and shift
sstMapped(sstMapped < sstCmapLim(1)-1) = sstCmapLim(1)-1;  % Clip bottom
sstMapped(sstMapped > sstCmapLim(2)-1) = sstCmapLim(2)-1;  % Clip top
sstMapped(sst == 0) = 0;   % Remap special value (displays as white)
sstMapped(land) = 1;       % Apply land mask (displays as dark gray)
sstMapped = uint8(round(sstMapped)); % Save storage with uint8

%== VISUALIZING THE MODIS SST SWATH ======================================

% Display the re-mapped SST as an indexed-color image using the colormap.
figure('Position',[100 100 size(sstMapped,2) size(sstMapped,1)])
colormap(cmap)
image(sstMapped)
set(gca,'Position',[0 0 1 1],'DataAspectRatio',[1 1 1]);

% Apply CONTOUR to the latitude and longitude geolocation arrays to
% create a lat-lon grid that illustrates the location of the swath and its
% spatial distortions. Use a contour interval of 5 degrees.
interval = 5;
latContourValues = interval ...
    * (ceil(min(sstLat(:))/interval) : floor(max(sstLat(:))/interval));
lonContourValues = interval ...
    * (ceil(min(sstLon(:))/interval) : floor(max(sstLon(:))/interval));
hold on
[c,h] = contour(sstLat,latContourValues,'k'); clabel(c,h)
[c,h] = contour(sstLon,lonContourValues,'k'); clabel(c,h)
hold off

% Use the geolocation fields as inputs to SURFACE to re-display the swath in
% latitude-longitude coordinates, the use COLORSCALE to add a color scale bar
% with units indicated in degrees Celsius.
figure
colormap(cmap)
axes('DataAspectRatio',[1 1 1],...
    'XTick',-180:5:180,'YTick',-90:5:90,...
    'XLim',[20 50],'YLim',[12 37],...
    'XGrid','on','YGrid','on')
xlabel('Longitude, degrees')
ylabel('Latitude, degrees')
title('Geolocated SST swath, July 4, 2001')
surface(sstLon, sstLat, double(sstMapped) + 1,...
    'Linestyle','none','CDataMapping','Direct')
% Add color scale
pos = get(gca,'Position');
set(gca,'Position',pos + [-0.06 0 0 0])
colorscale(sstCmapLim, sstDataLim, 5, 'vert',...
           'Position',[0.92 pos(2) 0.025 pos(4)])
xlabel('SST')
ylabel('degrees Celsius')

%== READING AVHRR GLOBAL NDVI GRID =======================================

% Read NDVI data from the AVHRR Climate file.
ndviFile = 'PAL_CLIMATE_JUL_01-10_2001.HDF';
ndvi = hdfread(ndviFile,'Data-Set-2');

% Read NDVI data description and attributes.
ndviInfo = hdfinfo(ndviFile);
ndviAttr = ndviInfo.SDS.Attributes;
ndviAddOffset   = double(ndviAttr(...
    strmatch('add_offset',  {ndviAttr.Name},'exact')).Value)
ndviScaleFactor = double(ndviAttr(...
    strmatch('scale_factor',{ndviAttr.Name},'exact')).Value)

%== MAPPING NDVI VALUES TO THE COLORMAP ==================================

% The approach here is the same as for the sea surface temperature data
% mappings, except that the linear relationship between the encoded NDVI
% values and actual NDVI levels is expressed in terms of an offset and
% scale factor rather than an intercept and slope. To convert an encoded
% data value to absolute NDVI units, one must subtract the offset, then
% multiply by the scale factor.

% Find the encoded NDVI limits that correspond to the absolute NDVI limits
% in ndviDataLim.
ndviCodelim = ndviDataLim/ndviScaleFactor + ndviAddOffset;

% Solve a 2-by-2 linear system:
%   b(1) * ndviCodelim + b(2) = ndviCmapLim - 1
b = (ndviCmapLim - 1)/[ndviCodelim; [1 1]];

% Finally, scale and shift, clip to stay within the appropriate part
% of the colormap, and convert to uint8. Also, be sure to individually
% remap the special values (1 denoting water and 0 missing data).
ndviMapped = b(1) * double(ndvi) + b(2);  % Scale and shift
ndviMapped(ndviMapped < ndviCmapLim(1)-1) = ndviCmapLim(1)-1;  % Clip bottom
ndviMapped(ndviMapped > ndviCmapLim(2)-1) = ndviCmapLim(2)-1;  % Clip top
ndviMapped(ndvi==0) = 130-1;   % Missing data (displays as gray)
ndviMapped(ndvi==1) = 129-1;   % Water (displays as blue)
ndviMapped = uint8(round(ndviMapped));  % Save storage with uint8

%== VISUALIZING THE NDVI GRID ============================================

% Show the NDVI grid over global latitude-longitude axes, and add a color
% scale.
figure
colormap(cmap)
image(ndviMapped,'XData',[-179.5 179.5],'YData',[89.5 -89.5])
set(gca,'YDir','normal','DataAspectRatio',[1 1 1],...
    'XTick',[-180 : 30 : 180],'YTick',[-90 : 30 : 90],...
    'XLim',[-180 180],'YLim',[-90 90])
xlabel('Longitude, degrees')
ylabel('Latitude, degrees')
title('NDVI from AVHRR, July 1-10, 2001')
% Add color scale
pos = get(gca,'Position');
set(gca,'Position',pos + [0 0.075 0 0])
colorscale( ndviCmapLim, ndviDataLim, 0.2, 'horiz','Position',[pos(1:3) 0.05])
xlabel('NDVI')

%== COMBINED SST AND NDVI DISPLAY ========================================

% Now show the NDVI and SST together on a latitude-longitude grid
% after zooming in on the region containing the SST swath. We add a
% color scale for each data set.
figure('Renderer','zbuffer')
colormap(cmap)
image(ndviMapped,'XData',[-179.5 179.5],'YData',[89.5 -89.5])
set(gca,'YDir','normal','DataAspectRatio',[1 1 1],...
    'XTick',-180:5:180,'YTick',-90:5:90,...  % Ticks every 5 degrees
    'XLim',[15 60],'YLim',[0 45]) % Axes limits in degrees
xlabel('Longitude, degrees')
ylabel('Latitude, degrees')
title('SST swath plotted over one-degree NDVI grid')
surface(sstLon, sstLat, ones(size(sstLon)), double(sstMapped) + 1,...
    'Linestyle','none','CDataMapping','Direct');
% Add color scales
pos = get(gca,'Position');
set(gca,'Position',pos + [-0.12 0 0 0])
colorscale(ndviCmapLim, ndviDataLim, 0.2, 'vert',...
           'Position',[0.9 pos(2) 0.025 pos(4)])
xlabel('NDVI')
colorscale(sstCmapLim, sstDataLim, 5, 'vert',...
           'Position',[0.8 pos(2) 0.025 pos(4)])
xlabel('SST')
ylabel('degrees Celsius')

author: Sable

repo: mcbench-benchmarks


text:
function r82poly2_print ( a, b, c, d, e, f )

%*****************************************************************************80
%
%% R82POLY2_PRINT prints a second order polynomial in two variables.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, D, E, F, the coefficients.
%
  fprintf ( 1, '  p(x,y) = %g * x^2 + %g * y^2 + %g * xy\n', a, b, c );
  fprintf ( 1, '         + %g * x + %g * y + %g\n', d, e, f );

  return
end


text:
% Mathematics Q2886713
% https://math.stackexchange.com/questions/2886713
% Prox Operator of L1 Norm with Linear Equality Constraint (Sum of Elements)
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     18/08/2018
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

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

generateFigures = OFF;

%% Simulation Parameters

numElements = 10;

%% Generate Data

vY = 10 * randn(numElements, 1);

paramGamma  = 0;
paramB      = 1;

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable vX(numElements);
    minimize( 0.5 * pow_pos(norm(vX - vY, 2), 2) + (paramGamma * norm(vX ,1)) );
    subject to
        sum(vX) == paramB;
cvx_end

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

%% Solution by Prox Function (Closed Form Solution)
%{
Solving
%}

hObjFun = @(vX) (0.5 * sum((vX - vY) .^ 2)) + (paramGamma * norm(vX ,1));

vXRef = vX;

vX = ProxL1NormSum(vY, paramGamma, paramB);

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

%% Display Results

figureIdx = figureIdx + 1;

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

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

author: RoyiAvital

repo: StackExchangeCodes


text:
% evaluate simulated data - ICA and NMF
function separ(sep0, offset0, path_data, path_res, prename, niter, savethis, hint, sep_how)
% sep_how: string i - ica, n - nmf -> method for separatdion
if ~exist('hint', 'var')
    hint = 0;
end

if ~exist('sep_how', 'var')
    sep_how = 'in';
end

fprintf('Separation of components... \n')
for rr = 1: length(offset0)
    fprintf('\n%g: ',rr)
    for ll=1 : length(sep0)
        fprintf('.')
        p.namedir = [prename num2str(100*sep0(ll)) 'offset_' num2str(offset0(rr))];
        cd ([path_data p.namedir])
        for mm=1:niter

            namefile = [p.namedir '-iter_' num2str(mm)];
            load ([namefile '.mat'])
            %         ims(psf);
            %         SaveImageFULL('psf', 'pf');

            if sum(sep_how == 'i')>0 %ICA
                [icasig{mm}, A{mm}, W{mm}] = fastica (dveccr, 'numOfIC', 2, 'g', 'tanh');
                icapixICA{mm} = reshape(A{mm},32, 32, 2);
            end
            if sum(sep_how == 'n')>0 %NMF
                ncomp = 2; %number of components to be separated
                if hint
                    dvec_ind = squeeze(reshape(double(array2im(dpixc_ind)), p.nx*p.ny, 1, 2)); % vectors of resized images
%                     [out, bg(mm), bg_im]=backgroundoffset(dpixc);
                    [out, bg(mm), bg_im]=backgroundoffset(dpixc, 'no', 5, 20, 8); %empirical values...
                    dvec_bg = bg(mm)*ones(1, p.nx*p.ny);
%                     dvec_bg = p.offset*ones(1, p.nx*p.ny); %changed for offset 10...

                    blinkmatrand = rand(p.Nt, ncomp);
                    winit = [blinkmatrand,ones(p.Nt,1)];             %random weights will be assigned to firts two and bg fixed

                    hinit = [dvec_ind'; dvec_bg];       %original 'true' points + background
%                     hinit = [rand(ncomp, p.nx*p.ny); dvec_bg];
                    ncomp = ncomp+1; %background added
%%%                    [w{mm},h{mm}, wtrace{mm},htrace{mm}]=nmf_test(double(dveccr'),ncomp+1,1,winit,hinit, [3], [3]);
[w{mm},h{mm}, wtrace,htrace]=nmf_test(double(dveccr'),ncomp+1,1,winit,hinit, [3], [3]);

                else
                    [w{mm},h{mm}]=nmf(double(dveccr'),ncomp,1);
                end
                icapixNMF{mm} = reshape(h{mm}',32,32,ncomp);
            end
            %             imstiled(icapixICA{mm});
            %             SaveImageFULL([p.namedir 'ICA_' num2str(mm)], 'p');
            %             imstiled(icapixNMF{mm});
            %             SaveImageFULL([p.namedir 'NMF_' num2str(mm)], 'p');
            close all

        end

        p.path_data = path_data;
        p.path = path_res;

        if savethis == 1
            fprintf('saving data \n');
            if ~(strcmp(p.path, p.path_data)) %not identical
                mkdir ([p.path p.namedir]);
                cd ([p.path p.namedir]);
            end
            save ([p.namedir '_separ'])
            writedata([],[],p,[p.namedir '_param'])
        end
    end
end

fprintf('\n')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [Next,Range] = EnvironmentalSelection(Obj,RefPoint,Range,N)
% The environmental selection of AR-MOEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

        %% Selection among feasible solutions
        % Non-dominated sorting
        [FrontNo,MaxFNo] = NDSort(Obj,N);%population 2*N
        Next = FrontNo < MaxFNo;
        % Select the solutions in the last front
        Last   = find(FrontNo==MaxFNo);
        Choose = LastSelection(Obj(Last,:),RefPoint,Range,N-sum(Next));
        Next(Last(Choose)) = true;
        %Population = Population(Next,:);
        % Update the range for normalization
        Range(2,:) = max(Obj,[],1);
        Range(2,Range(2,:)-Range(1,:)<1e-6) = 1;
end

function Remain = LastSelection(PopObj,RefPoint,Range,K)
% Select part of the solutions in the last front

    N  = size(PopObj,1);
    NR = size(RefPoint,1);

    %% Calculate the distance between each solution and point
    Distance    = CalDistance(PopObj-repmat(Range(1,:),N,1),RefPoint);
    Convergence = min(Distance,[],2);

    %% Delete the solution which has the smallest metric contribution one by one
    [dis,rank] = sort(Distance,1);
    Remain     = true(1,N);
    while sum(Remain) > K
        % Calculate the fitness of noncontributing solutions
        Noncontributing = Remain;
        Noncontributing(rank(1,:)) = false;
        METRIC = sum(dis(1,:)) + sum(Convergence(Noncontributing));
        Metric = inf(1,N);
        Metric(Noncontributing) = METRIC - Convergence(Noncontributing);
        % Calculate the fitness of contributing solutions
        for p = find(Remain & ~Noncontributing)
            temp = rank(1,:) == p;
            noncontributing = false(1,N);
            noncontributing(rank(2,temp)) = true;
            noncontributing = noncontributing & Noncontributing;
            Metric(p) = METRIC - sum(dis(1,temp)) + sum(dis(2,temp)) - sum(Convergence(noncontributing));
        end
        % Delete the worst solution and update the variables
        [~,del] = min(Metric);
        temp    = rank ~= del;
        dis     = reshape(dis(temp),sum(Remain)-1,NR);
        rank    = reshape(rank(temp),sum(Remain)-1,NR);
        Remain(del) = false;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function analyse_dev_lin_and_qual_fusion(key,scores_obj_array,scr_lin,scr_qual,prior)
% Analyses the dev scores for all input systems (including fusions)
% to so that there relative strength can be compared.  It
% calculates the percentage of target trials calculated by each
% system and the actual dcf, min dcf and prbep values for each
% system.
% Inputs:
%   key: The Key for the dev trials.
%   scores_obj_array: An array of score objects (one object for
%     each system) that was fused.
%   scr_lin: An object of type Scores containing the result of the
%     linear fusion.
%   scr_qual: An object of type Scores containing the result of the
%     quality fusion.
%   prior: The effective target prior.

assert(nargin==5)
assert(isa(key,'Key'))
assert(isa(scr_lin,'Scores'))
assert(isa(scr_qual,'Scores'))
assert(key.validate())
assert(scr_lin.validate())
assert(scr_qual.validate())

trialmask = key.to_ndx().trialmask;
dev_stacked_scores = stackScores(key,scores_obj_array);
dev_lin_fusion = scr_lin.scoremat(trialmask(:))';
dev_qfusion = scr_qual.scoremat(trialmask(:))';
dev_scores = [dev_stacked_scores;dev_lin_fusion;dev_qfusion];

true_prop = sum(key.tar(:))/(sum(key.tar(:))+sum(key.non(:)));
logprint(Logger.Info,'True Dev target proportion is %g%%. System estimates are:\n',100*true_prop);
logprint(Logger.Info,'%s',sprintfmatrix(100*estimate_target_proportions(dev_scores,true_prop)'));

numsystems = size(dev_stacked_scores,1);
mf = size(dev_scores,1);
logprint(Logger.Info,'\n\nAnalysing Dev: %i subsystems and %i fusions:\n',numsystems,mf-numsystems);
E = evalAll(dev_scores,key,prior);
logprint(Logger.Info,'\nDev all: dcf*100, mindcf*100, prbep:\n%s',sprintfmatrix(E,1:mf));
end

author: nesl


text:
function [status,varargout] = spellcheck(text)

%SPELLCHECK checks the spelling of word(s) and returns suggestions if misspelled.
%
% [status,suggestions] = spellcheck(text);
% status = spellcheck(text);
%
%     text:         a word or string of words separated by space(s).
%     status:       returns '1' if words exist in dictionary or '0' otherwise.
%     suggestions:  an array of suggested corrections of the misspelled word
%                   or otherwise returns a message 'Correct Spelling'.
%
% Examples:
%      [status,suggestions] = spellcheck('Hellow');
%      [status,suggestions] = spellcheck('My Name is Fahad!');
%      status = spellcheck('Emirates');
%
%   Copyright 2004 Fahad Al Mahmood
%   Version: 1.0 $  $Date: 30-Jun-2004
%   Version: 1.5 $  $Date: 01-Jul-2004 (Multiple words check & status,
%                                       thanks Michael Kleder!)
%   Version: 1.6 $  $Date: 11-Jul-2004 (Fixed the 'no suggestions' result,
%                                       thanks Sander Stepanov!)

% Separating string of words into arrays of words.
k=1;
temp='';
for n=1:length(text);
    if ~isspace(text(n))
        temp = [temp text(n)];
    else
        if ~isspace(text(n-1))
            words{k} = temp;
        end
        temp='';
        k=k+1;
    end
end
words{k} = temp;

% Opening MS Word and Starting the Spelling Check

Doc = actxserver('Word.Application');
m=1;
for n=1:length(words)
    if ~isempty(words{n})
        status(m) = invoke(Doc,'CheckSpelling',words{n},[],1,1);
        if nargout==2
            invoke(Doc.Documents,'Add');
            X = invoke(Doc,'GetSpellingSuggestions',words{n});
            count = get(X,'Count');
            suggestions{m,1} = words{n};
            if count~=0
                for k=2:count+1
                    suggestion = invoke(X,'Item',k-1);
                    suggestions{m,k} = get(suggestion,'Name');
                end
            elseif count==0 & status(m)==1
                suggestions{m,2} = 'Correct Spelling';
            elseif count==0 & status(m)==0
                suggestions{m,2} = 'No Suggestions!';
            end
            varargout = {suggestions};
        end
        m=m+1;
    end
end

status = all(status);

invoke(Doc,'Quit');
delete(Doc);

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
classdef MimUserSavedItem < GemListItem
    % MimUserSavedItem. Part of the gui for the TD MIM Toolkit.
    %
    %     This class is used internally within the TD MIM Toolkit to help
    %     build the user interface.
    %
    %     MimUserSavedItem is used in a list box to show a list of cached
    %     items that are available for the user to select, such as markers or
    %     manual segmentations.
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        CacheFileName
        VisibleName
    end

    properties (Constant)
        UserSavedItemTextHeight = 21
    end

    properties (Access = private)
        TextControl

        LoadCallback
        DeleteCallback
        RenameCallback
        DuplicateCallback
        AddCallback
    end

    properties (Constant, Access = private)
        UserSavedItemFontSize = 14
    end

    methods
        function obj = MimUserSavedItem(parent, name, tooltip, load_callback, delete_callback, rename_callback, duplicate_callback, add_callback)
            obj = obj@GemListItem(parent, name);
            obj.CacheFileName = name;
            obj.VisibleName = name;
            obj.TextHeight = obj.UserSavedItemTextHeight;

            if nargin > 0
                obj.LoadCallback = load_callback;
                obj.DeleteCallback = delete_callback;
                obj.RenameCallback = rename_callback;
                obj.DuplicateCallback = duplicate_callback;
                obj.AddCallback = add_callback;

                obj.TextControl = GemText(obj, obj.VisibleName, tooltip, 'ItemSelect');
                obj.TextControl.FontSize = obj.UserSavedItemFontSize;
                obj.TextControl.HorizontalAlignment = 'left';
                obj.AddTextItem(obj.TextControl);
            end
        end

        function Resize(obj, location)
            size_changed = ~isequal(location, obj.Position);

            % Don't call the parent class
            Resize@GemVirtualPanel(obj, location);

            obj.TextControl.Resize(location);

            % A resize may change the location of the highlighted item
            if size_changed
                obj.Highlight(false);
            end
        end
    end

    methods (Access = protected)

        function ItemLeftClicked(obj, src, eventdata)
            ItemLeftClicked@GemListItem(obj, src, eventdata);
            obj.LoadCallback(obj.CacheFileName);
        end

        function ItemRightClicked(obj, src, eventdata)
            ItemRightClicked@GemListItem(obj, src, eventdata);

            if isempty(get(obj.TextControl.GraphicalComponentHandle, 'uicontextmenu'))
                context_menu = uicontextmenu;
                context_menu_delete = uimenu(context_menu, 'Label', 'Delete', 'Callback', @obj.DeleteCachedItem);
                context_menu_rename = uimenu(context_menu, 'Label', 'Rename', 'Callback', @obj.RenameCachedItem);
                context_menu_duplicate = uimenu(context_menu, 'Label', 'Duplicate', 'Callback', @obj.DuplicateCachedItem);
                context_menu_add = uimenu(context_menu, 'Label', 'Add', 'Callback', @obj.AddCachedItem);
                obj.SetContextMenu(context_menu);
            end
        end
    end

    methods (Access = private)
        function DeleteCachedItem(obj, ~, ~)
            parent_figure = obj.GetParentFigure;
            parent_figure.ShowWaitCursor;
            obj.DeleteCallback(obj.CacheFileName);

            % Note that at this point obj may have been deleted, so we can no longer use it
            parent_figure.HideWaitCursor;
        end

        function RenameCachedItem(obj, ~, ~)
            parent_figure = obj.GetParentFigure;
            parent_figure.ShowWaitCursor;
            obj.RenameCallback(obj.CacheFileName);
            parent_figure.HideWaitCursor;
        end

        function DuplicateCachedItem(obj, ~, ~)
            parent_figure = obj.GetParentFigure;
            parent_figure.ShowWaitCursor;
            obj.DuplicateCallback(obj.CacheFileName);
            parent_figure.HideWaitCursor;
        end

        function AddCachedItem(obj, ~, ~)
            parent_figure = obj.GetParentFigure;
            parent_figure.ShowWaitCursor;
            obj.AddCallback();
            parent_figure.HideWaitCursor;
        end
    end
end


text:
function [D,DJacob,DHess,pDpt]=DCoordTurn2D(x,q0,qTurn,qLin)
%DCOORDTURN2D The continuous-time diffusion matrix function for a 2D
%             coordinated turn model with a Cartesian state. The turn rate
%             can be specified in terms of a turn rate in radians per
%             second, or in terms of a transversal acceleration.
%             Additionally, a linear acceleration can be given. This
%             diffusion matrix goes with the drift functions
%             aCoordTurn2DOmega and aCoordTurn2DTrans.
%
%INPUTS: x The target state for 2D motion. If there is no linear
%          acceleration (acceleration along the direction of motion), then
%          x can either be x=[x;y;xdot;ydot;omega], where omega is the turn
%          rate estimate in radians per second counterclockwise from the
%          x-axis or x=[x;y;xdot;ydot;at] where at is the transversal
%          acceleration, which is orthogonal to the velocity and is defined
%          such that positive values of at map to positive values of omega.
%          If there is a linear acceleration, then the target state is
%          either x=[x;y;xdot;ydot;omega;al] where omega is the turn rate
%          and al is the linear acceleration or the target state is
%          x=[x;y;xdot;ydot;at;al] if the turn is expressed in terms of a
%          transversal acceleration. The dimensionality of the state is
%          used to determine whether a linear acceleration component is
%          present. The linear acceleration component changes the speed.
%          That means that it acts in the direction of the velocity vector.
%          The number of columns in x determine how many copies of D are
%          returned.
%       q0 The power spectral density of the process noise of the velocity
%          components. It is assumed to be the same in both dimensions. It
%          covers perturbations from an ideal coordinated turn trajectory
%          and has units of m^2/s^3. If no process noise is desired for the
%          velocity (i.e. all perturbations should be covered by noise on
%          the turn component and on the linear acceleration), then a value
%          of zero should be passed.
%    qTurn If the turn is specified in terms of a turn rate in radians per
%          second, then this is the power spectral density of the turn rate
%          noise having units of radians squared per seconds cubed. If the
%          turn is expressed in terms of a transverse acceleration, then
%          this is the power spectral density of the transverse
%          acceleration noise, having units of m^2/s^5.
%     qLin This parameter is only needed if a linear acceleration is
%          present. It is the power spectral density of the linear
%          acceleration noise having units of m^2/s^5.
%
%OUTPUTS: D The diffusion matrix of a 2D continuous-time turning model
%           where the velocity is given as a Cartesian vector. If x has
%           N columns, then N copies of D are returned with D(:,:,i) being
%           the ith one.
%   DJacob, DHess The xDimXmXxDim, xDimXmXxDimXxDim matrices of first and
%           second partial derivatives of the elements of D with respect to
%           x. These are all zero, since D is a constant. it is the same
%           for all D and is not repeated N times.
%      pDpt The xDimX2 partial derivative of D with respect to time. This
%           is all zeros, because D is a constant.
%
%The basic 2D coordinated turn model in Cartesian coordinates is described
%in Section VA of [1]. When the turn rate is something that must be
%estimated, it is assumed that the continuous-time turn rate model is
%omegaDot=-(1/tau)*Omega+noise
%Note that the ordering of the state elements assumed by this function
%differs from the ordering of the state elements assumed in the paper.
%
%The 2D coordinates turn model in Cartesian coordinates is also described
%in  Chapter 4.2.3 of [2].
%
%The concept of using the transversal acceleration instead of the turn rate
%is not discussed in either of those references. It is, however, mentioned
%in [3], though no differential equations are given and a more detailed
%reference cited therein is a hard-to-get dissertation in French. The use
%of transversal acceleration is discussed in more detail in [4], though
%expressions are given when considering the 2D velocity are broken into
%components of heading and speed rather than in Cartesian space. The
%generalization to Cartesian space is not difficult and is done here.
%
%A starting point for setting q0 is to use processNoiseSuggest with
%'PolyKal-ROT' and order=1. The noise must cover small velocity deviations
%from the ideal turn model. A starting point for setting qTurn when a turn
%rate is given is to use processNoiseSuggest with 'PolyKal-ROT' and
%order=1. A starting point for setting qTurn when using a linear
%acceleration and for qLin is to use processNoiseSuggest with 'PolyKal-ROT'
%and order=2. The order chosen in the suggestion function just depends on
%the number of derivatives of time present.
%
%The corresponding drift functions are given by the functions
%aCoordTurn2DOmega and aCoordTurn2DTrans. The corresponding discrete-time
%functions are FCoordTurn2D and QCoordTurn. However, note that the
%discrete-time functions with unknown noise is a direct-discrete model and
%not a discretization of the continuous-time model.
%
%REFERENCES:
%[1] X. R. Li and V. P. Jilkov, "Survey of maneuvering target tracking.
%    Part I: Dynamic models," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 39, no. 4, pp. 1333-1364, Oct. 2003.
%[2] S. Blackman and R. Popoli, Design and Analysis of Modern Tracking
%    Systems. Norwood, MA: Artech House, 1999.
%[3] P. Vacher, I. Barret, and M. Gauvrit, "Design of a tracking algorithm
%    for an advanced ATC system," in Multitarget-Multisensor Tracking:
%    Applications and Advances, Y. Bar-Shalom, Ed. Norwood, MA: Artech
%    House, 1992, vol. II, ch. 1.
%[4] H. A. P. Blom, R. A. Hogendoorn, and B. A. van Doorn, "Design
%    of a multisensor tracking system for advanced air traffic control," in
%    Multitarget-Multisensor Tracking: Applications and Advances, Y. Bar-
%    Shalom, Ed. Norwood, MA: Artech House, 1992, vol. II, ch. 2.
%
%July 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(x,2);

rootQ0=sqrt(q0);
rootQTurn=sqrt(qTurn);

numDim=length(x);
switch(numDim)
    case 5%There is no linear acceleration
        %Equation 61 in Li's paper combined with Equation 67.
        D=[0,       0,      0;%Row for x-component noise.
           0,       0,      0;%Row for y-component noise.
           rootQ0,  0,      0;%Row for velocity-x noise.
           0,       rootQ0, 0;%Row for velocity-y noise.
           0,       0,      rootQTurn];%Row for turn noise.
       m=3;
    case 6%There is a linear acceleration component.
        rootQLin=sqrt(qLin);
        %Similar to Equation 61 in Li's paper combined with Equation 67,
        %but with an added row for noise in the linear acceleration
        %component.
        D=[0,       0,      0,  0;%Row for x-component noise.
           0,       0,      0,  0;%Row for y-component noise.
           rootQ0,  0,      0,  0;%Row for velocity-x noise.
           0,       rootQ0, 0,  0;%Row for velocity-y noise.
           0,       0, rootQTurn,0;%Row for turn rate noise.
           0,       0,      0,  rootQLin];%Row for linear accel. noise.
       m=4;
    otherwise
        error('The length of x is neither 5 nor 6.');
end

if(nargout>1)
    DJacob=zeros(numDim,m,numDim);
    if(nargout>2)
        DHess=zeros(numDim,m,numDim,numDim);
        if(nargout>3)
            pDpt=zeros(numDim,m);
        end
    end
end

if(N>1)
   D=repmat(D,[1,1,N]);
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [gx] = g_goNogo(x,P,u,in)
% generates the probability of 'go' choice: P(go) = sig(V/temperature+bias)
% function [gx] = g_goNogo(x,P,u,in)
% IN:
%   - x: value of the 'go' option
%   - P: log-temperature and bias
%   - u: [useless]
%   - in: [useless]
% OUT:
%   - gx: P(go) = sig(V/temperature+bias)

gx = 1./(1+exp(P-x));

author: MBB-team

repo: VBA-toolbox


text:
function writeNewtExperiment(model,metData,metSampleNames,fileName,param)
%writes out an experimental data file for an model exported to SBML and
%imported into newteditor

% version\t1.0
% name\tsample experiment data
% description\tAdenoid Cystic Carcinoma 2014 vs 2019
% color\t0\t#FFFFFF\t100\t#FF0000
% gene\t2014\t2019
% RB1\t36\t12
% TP53\t36\t72
% CDKN2A\t0\t14
% MDM2\t0\t5
% CCNE\t0\t7

[nNodes,nSamples]=size(metData);
[nMets,nRxns]=size(model.S);
if nNodes ~= nMets
    error('size(metData,1) must equal size(model.S,1)')
end

if ~exist('fileName','var')
    fileName = [pwd filesep 'newtExpData.txt'];
end

fid = fopen(fileName,'w');

if ~exist('param','var')
    param=struct;
end

if ~isfield(param,'version')
    param.version = datestr(now,30);
end

if ~isfield(param,'name')
    if isfield(model,'modelID')
        param.name = model.modelID;
    else
        param.name = 'aModel';
    end
end

if ~isfield(param,'description')
    if isfield(model,'description')
        param.description = model.description;
    else
        param.description = 'aDescription';
    end
end

% color\t0\t#FFFFFF\t100\t#FF0000
%color -100  0 #00FF00 100 #0000FF
if ~isfield(param,'color')
   param.color.minValue = -100;
   param.color.minColor = '#FF0000';
   param.color.zeroValue = 0;
   param.color.zeroColor = '#FF0000';
   param.color.maxValue =  100;
   param.color.maxColor = '#0000FF';
%    if ~isfield(param.color,'string')
%         param.color.string = [param.color.minValue '\t' param.color.minColor '\t' param.color.zeroValue '\t' param.color.zeroColor '\t' param.color.maxValue '\t' param.color.maxColor];
%    end
end

%fprintf('name\tsample experiment data\r\ndescription\tAdenoid Cystic Carcinoma 2014 vs 2019\r\nel\t2014\t2019\r\nRB1\t36\t12\r\nTP53\t36\t72\r\nCDKN2A\t0\t14\r\nMDM2\t0\t5\r\nCCNE\t0\t7\r')

%header
fprintf(fid,'%s\t%s\n','version', param.version);
fprintf(fid,'%s\t%s\n','name', param.name);
fprintf(fid,'%s\t%s\n','description', param.description);
fprintf(fid,'%s\t%i\t%s\t%i\t%s\t%i\t%s\n','color', param.color.minValue,param.color.minColor,param.color.zeroValue,param.color.zeroColor,param.color.maxValue ,param.color.maxColor);

%experiment names
fprintf(fid,'%s\t','node');
if ~exist('metSampleNames','var')
    for n=1:nSamples
        if n~=nSamples
            fprintf(fid,'%s\t',['exp' int2str(n)]);
        else
            fprintf(fid,'%s',['exp' int2str(n)]);
        end
    end
    fprintf(fid,'\n');
else
    if ischar(metSampleNames)
        metSampleNames={metSampleNames};
    end
    for n=1:nSamples
        if n~=nSamples
            fprintf(fid,'%s\t',metSampleNames{n});
        else
            fprintf(fid,'%s',metSampleNames{n});
        end
    end
    fprintf(fid,'\n');
end

%assume that the map in newt was exported via an sbml file, in which case
%the node ids are in sbml format
convertedMets = convertSBMLID(model.mets);

%experimental data
useM=0;
for i=1:nMets
    if useM
        fprintf(fid,'%s',['M_' convertedMets{i}]);
    else
        fprintf(fid,'%s\t',model.metNames{i});
    end

    for j=1:nSamples
        if j~=nSamples
            fprintf(fid,'%i\t', full(metData(i,j)));
        else
            fprintf(fid,'%i', full(metData(i,j)));
        end
    end
    fprintf(fid,'\n') ;
end

fclose(fid);

author: opencobra

repo: cobratoolbox


text:
% load sample EEG data
load sampleEEGdata

% wavelet parameters
min_freq = 2;
max_freq = 48;
num_frex = 30;

% other wavelet parameters
frequencies = logspace(log10(min_freq),log10(max_freq),num_frex);
time = -1:1/EEG.srate:1;
half_of_wavelet_size = (length(time)-1)/2;

% FFT parameters (use next-power-of-2)
n_wavelet     = length(time);
n_data        = EEG.pnts*EEG.trials;
n_convolution = n_wavelet+n_data-1;
n_conv_pow2   = pow2(nextpow2(n_convolution));
wavelet_cycles= 4;

chan2plot = 'fz';

% define baseline period
baselinetime = [ -500 -200 ]; % in ms

% convert baseline window time to indices
[junk,baselineidx(1)]=min(abs(EEG.times-baselinetime(1)));
[junk,baselineidx(2)]=min(abs(EEG.times-baselinetime(2)));

tf_data = zeros(2,length(frequencies),EEG.pnts);

fft_data = fft(reshape(EEG.data(strcmpi(chan2plot,{EEG.chanlocs.labels}),:,:),1,[]),n_conv_pow2);

for fi=1:length(frequencies)

    % create wavelet and get its FFT
    wavelet = exp(2*1i*pi*frequencies(fi).*time) .* exp(-time.^2./(2*( wavelet_cycles /(2*pi*frequencies(fi)))^2));
    fft_wavelet = fft(wavelet,n_conv_pow2);
    fft_wavelet = fft_wavelet./max(fft_wavelet);

    % run convolution
    convolution_result_fft = ifft(fft_wavelet.*fft_data,n_conv_pow2);
    convolution_result_fft = convolution_result_fft(1:n_convolution);
    convolution_result_fft = convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size);
    convolution_result_fft = reshape(convolution_result_fft,EEG.pnts,EEG.trials);

    % put power data into time-frequency matrix
    tf_data(1,fi,:) = mean(abs(convolution_result_fft).^2,2);
    tf_data(2,fi,:) = median(abs(convolution_result_fft).^2,2);
end

baseline_power = squeeze(mean(tf_data(1,:,baselineidx(1):baselineidx(2)),3));
dbconverted = 10*log10( squeeze(bsxfun(@rdivide,tf_data(1,:,:),baseline_power) ));

%%

surf(dbconverted(4:end,155:539))
h=findobj('Type','surface');
%set(h,'CData',double(asdf(4:end,155:539)))
set(gca,'clim',[-3 3],'yscale','log'), shading interp, axis off
set(gcf,'color','k')
cmap=(1+[cos(linspace(0,pi*2,100)); sin(linspace(0,pi*2,100)); cos(linspace(0,pi*2,100))])/2;
colormap(cmap')

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
% Faust Synthetic
test_idx=[80:99];
pairs_list = [];
for i=1:20
    for j = 1:20
        pairs_list = [pairs_list; test_idx(i),test_idx(j)];
    end
end

% delete_idx = floor(pairs_list(:,1)/10) == floor(pairs_list(:,2)/10);
% pairs_list(delete_idx,:)=[];

fileID = fopen('test_pairs.txt','w');
for i=1:size(pairs_list,1)
    id1 = sprintf('%03d', pairs_list(i,1));
    id2 = sprintf('%03d', pairs_list(i,2));
    fprintf(fileID,['tr_reg_', id1, '.mat tr_reg_' ,id2, '.mat\n']);
end
fclose(fileID)

author: OshriHalimi


text:
function Speciation
% Speciation
%    using MATLAB
%
%   $Ekkehard Holzbecher  $Date: 2006/05/03 $
%--------------------------------------------------------------------------
toll = 1.e-15;                    % tolerance
nmax = 50;                        % max. no. of iterations
Se = [-1 -1 1];                   % reaction matrix
logc = [1.e-10; 1.e-10; 0];       % initial guess (log)
logK = [-0.93];                   % equilibrium constants (log)
logu = [-0.301; 0];               % total concentrations (log)

ln10 = 2.3026;
n=size(Se,1); m=size(Se,2);
S1 = Se(:,1:m-n);
S2 = Se(:,m-n+1:m);
S1star = -S2\S1;
U = [eye(m-n),S1star'];

c=exp(ln10*logc);
u(1:m-n,1) = exp(ln10*logu);
err = toll+1; nit = 0;
while (nit < nmax && err > toll),
    nit = nit+1;
    F = [U*c-u; Se*logc-logK];
    DF = [U; Se*diag((1/ln10)./c)];
    dc = -DF\F;
    cn = max(c+dc,0.005*abs(c));
    err = max(abs(cn-c));
    c = cn;
end
display (['Species concentrations obtained after ' num2str(nit) ' iterations:']);
c
exp(ln10*logK)-c(3)/c(1)/c(2)


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the column row of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);

% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

epsilon = 0.12;
W = rand(L_out, 1 + L_in) * 2 * epsilon - epsilon;

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

end

author: AvaisP


text:
function stroud_test163 ( )

%*****************************************************************************80
%
%% STROUD_TEST163 tests the rules for CN with Gegenbauer weight on monomials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2010
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;

  alpha_test = [ -0.5, 0.0, 0.5, 1.0, 1.5 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'STROUD_TEST163\n' );
  fprintf ( 1, '  Demonstrate the use of quadrature rules for the region\n' );
  fprintf ( 1, '  CN_GEG, that is, the hypercube [-1,+1]^N, with the\n' );
  fprintf ( 1, '  weight W(ALPHA;X) = product ( 1 <= I <= N )\n' );
  fprintf ( 1, '    (1-X(I)^2)^ALPHA\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use the formulas to integrate various monomials of\n' );
  fprintf ( 1, '  the form X(1)^E(1) * X(2)^E(2) * ... X(N)^E(N)\n' );
  fprintf ( 1, '  and compare to the exact integral.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The precision of each formula is known, and we only use\n' );
  fprintf ( 1, '  a formula if its precision indicates it should be able to\n' );
  fprintf ( 1, '  produce an exact result.\n' );

  for n = 1 : 6

    for test = 1 : test_num

      alpha = alpha_test(test);

      expon(1:n) = 0;
      cn_geg_test ( n, alpha, expon );

    end

    for test = 1 : test_num

      alpha = alpha_test(test);

      expon(1:n) = 0;
      expon(n) = 1;
      cn_geg_test ( n, alpha, expon );

    end

    if ( 2 <= n )

      for test = 1 : test_num

        alpha = alpha_test(test);

        expon(1:n) = 0;
        expon(1) = 1;
        expon(2) = 1;
        cn_geg_test ( n, alpha, expon );

      end

    end

    for test = 1 : test_num

      alpha = alpha_test(test);

      expon(1:n) = 0;
      expon(1) = 2;
      cn_geg_test ( n, alpha, expon );

    end

  end

  return
end

author: johannesgerer


text:
function [it,Rt,Pc] = splitInvFrame(F)

% SPLITINVFRAME  Split inverse frame information.
%   [iT,Rt] = SPLITINVFRAME(F), for a frame F, returns the inverse
%   translation vector iT and the transposed rotation matrix Rt. The frame
%   F can be either a 7-vector F=[T;Q] or a structure containing, at least,
%   the fields F.it and F.Rt.
%
%   [iT,Rt,Pc] = SPLITINVFRAME(F) returns in addition the conjugated Pi
%   matrix Pc. See toFrame for explanations on matrices Pi and Pc. If F is
%   a frame structure, it needs to contain the field F.Pc.
%
%   See also FRAME, SPLITFRAME, UPDATEFRAME, TOFRAME.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if isa(F,'struct')
    it = F.it;
    Rt = F.Rt;

    if nargout > 2
        Pc = F.Pc;
    end

else % F is a 7-vector
    Rt = q2R(F(4:7))';
    it = -Rt*F(1:3);

    if nargout > 2
        Pc = q2Pi(q2qc(F(4:7)));
    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function write_sicd_security_tags(obj)
%WRITE_SICD_SECURITY_TAGS Writes the NITF security tags
%
% It's not expected that the typical user would call this function, rather
% it is called as part of a bigger NITF/SICD writer.
%
% Inputs:
%       obj:       This is the 'hidden' object-specific handle.  It is
%                  similar to 'this' in Java.
%
% References:
%       MIL-STD-2500C, Department of Defense Interface Standard
%                      National Imagery Transmission Format version 2.1
%                      01 May 2006
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Security is limited in the SICD XML structure.  Without explicit tags
% passed, we make a guess from XML and fill in what we can.
if isfield(obj.sicdmeta,'CollectionInfo') && ...
        isfield(obj.sicdmeta.CollectionInfo,'Classification')&& ...
        ~isempty(obj.sicdmeta.CollectionInfo.Classification)
    classification = obj.sicdmeta.CollectionInfo.Classification;
else
    classification = ' ';
end
codes = strsplit(classification,'/');  % Cell array
codes = strjoin(cellfun(@(x) x(1:min(2,end)), codes(2:end), 'UniformOutput', false));

% Tag name, field length, default value
tag_info = { ...
    'CLAS', 1, classification(1); ...
    'CLSY', 2, 'US'; ...
    'CODE', 11, codes; ...
    'CTLH', 2, ''; ...
    'REL', 20, ''; ...
    'DCTP', 2, ''; ...
    'DCDT', 8, ''; ...
    'DCXM', 4, ''; ...
    'DG', 1, ''; ...
    'DGDT', 8, ''; ...
    'CLTX', 43, ''; ...
    'CATP', 1, ''; ...
    'CAUT', 40, ''; ...
    'CRSN', 1, ''; ...
    'SRDT', 8, ''; ...
    'CTLN', 15, ''};

if isfield(obj.sicdmeta,'NITF') && isfield(obj.sicdmeta.NITF,'security')
    sec_tags = obj.sicdmeta.NITF.security;
else
    sec_tags = struct();
end

for i=1:size(tag_info,1)
    if isfield(sec_tags, tag_info{i,1})
        fwriten(obj.FID, sec_tags.(tag_info{i,1}), tag_info{i,2});
    else
        fwriten(obj.FID, tag_info{i,3}, tag_info{i,2});
    end
end

end

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

author: ngageoint


text:
function [ n_data, x, fx ] = erf_values ( n_data )

%*****************************************************************************80
%
%% ERF_VALUES returns some values of the ERF or "error" function.
%
%  Discussion:
%
%    The error function is defined by:
%
%      ERF(X) = ( 2 / sqrt ( PI ) * integral ( 0 <= T <= X ) exp ( - T^2 ) dT
%
%    In Mathematica, the function can be evaluated by:
%
%      Erf[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 21;

  fx_vec = [ ...
     0.0000000000000000E+00, ...
     0.1124629160182849E+00, ...
     0.2227025892104785E+00, ...
     0.3286267594591274E+00, ...
     0.4283923550466685E+00, ...
     0.5204998778130465E+00, ...
     0.6038560908479259E+00, ...
     0.6778011938374185E+00, ...
     0.7421009647076605E+00, ...
     0.7969082124228321E+00, ...
     0.8427007929497149E+00, ...
     0.8802050695740817E+00, ...
     0.9103139782296354E+00, ...
     0.9340079449406524E+00, ...
     0.9522851197626488E+00, ...
     0.9661051464753107E+00, ...
     0.9763483833446440E+00, ...
     0.9837904585907746E+00, ...
     0.9890905016357307E+00, ...
     0.9927904292352575E+00, ...
     0.9953222650189527E+00 ];

  x_vec = [ ...
     0.0E+00, ...
     0.1E+00, ...
     0.2E+00, ...
     0.3E+00, ...
     0.4E+00, ...
     0.5E+00, ...
     0.6E+00, ...
     0.7E+00, ...
     0.8E+00, ...
     0.9E+00, ...
     1.0E+00, ...
     1.1E+00, ...
     1.2E+00, ...
     1.3E+00, ...
     1.4E+00, ...
     1.5E+00, ...
     1.6E+00, ...
     1.7E+00, ...
     1.8E+00, ...
     1.9E+00, ...
     2.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
%% This is for testing the SerialLink models in the robotics Toolboxg
function tests = SerialLinkModelsTest
  tests = functiontests(localfunctions);
end

function models_test(tc)
    models

    names = models();
    verifyTrue(tc, iscell(names) );

    models('6dof')

    names = models('6dof');
    verifyTrue(tc, iscell(names) );
end

function all_models_test(tc)

    all = models();

    for model = all'
        try
            eval( [model{1} ';'] );
        catch
            verifyFail(tc, sprintf('model %s failed to load', model{1}) );
        end

        % find all the SerialLink objects
        vars = whos;
        k = find( cellfun(@(s) strcmp(s,'SerialLink'), {vars.class}) );
        if isempty(k)
            verifyFail(tc, sprintf('no SerialLink models loaded by %s', model{1}) );
        end

        for kk = k
            clear(vars(kk).name)
        end

%         if exist('qz', 'var') ~= 1
%             verifyFail(tc, sprintf('model %s doesn''t set qz', model{1}));
%         end
        clear qz
    end
end

function query_test(tc)
    names = models('6dof');
    for model = names'
        % attempt to load the model
        try
            eval( [model{1} ';'] );
        catch
            verifyFail(tc, 'model %s failed', model{1});
        end

        % find the number of DOF in the model
        vars = whos;

        k = find( cellfun(@(s) strcmp(s,'SerialLink'), {vars.class}) );

        if ~isempty(k)
            if eval( [vars(k).name '.n'] ) ~= 6
                verifyFail(tc, 'query failed');
            end
        end

        clear(vars(k).name)

    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function createReadmeSession(info)
% function createReadmeSession
%
% Dialog box to get basic descriptors of session.
% Writes these to Readme file.
% Called by mrCreateReadme.m
%
% djh, 9.4.01
% ras, 01/12/04: added 'RessCoil' and 'fMRICoil' to list of available
% coils.
% ras, 04/18/04: also added 'Posthead-ReceiveOnly' to coil list.
% ras, 03/30/05: uses info from existing readme if supplied.
% ras, 08/11/06: tries to get protocl from headers.
global mrSESSION dataTYPES

if notDefined('info')
    % set blank default vals, or vals from mrSESSION
    defaults = {...
            mrSESSION.sessionCode, mrSESSION.description,...
            mrSESSION.subject,'','','','',...
            num2str(mrSESSION.examNum),'',''...
        };

    % try to get protocol from header
    try
        w = dir('Raw/Anatomy/Inplane/I*');
        firstIfile = fullfile('Raw', 'Anatomy', 'Inplane', w(1).name);
        [img hdr] = mrReadIfile(firstIfile);
        defaults{10} = hdr.series.prtcl;
    end

else
    % get defaults from existing Readme.txt file
    defaults = {...
            info.session, info.description,...
            info.subject, info.operator, ...
            info.operator, info.magnet, info.coil...
            info.examNum, info.sliceOrientation,...
            info.protocol...
    };
end

c=0;

title = 'SESSION IDENTIFIERS';

c=c+1;
dlg(c).string = 'Session code:';
dlg(c).fieldName = 'sessionCode';
dlg(c).style = 'edit';
dlg(c).value = defaults{1};

c=c+1;
dlg(c).string = 'Description:';
dlg(c).fieldName = 'description';
dlg(c).style = 'edit';
dlg(c).value = defaults{2};
c=c+1;
dlg(c).string = 'Subject:';
dlg(c).fieldName = 'subject';
dlg(c).style = 'edit';
dlg(c).value = defaults{3};

c=c+1;
dlg(c).string = 'Operator:';
dlg(c).fieldName = 'operator';
dlg(c).style = 'edit';
dlg(c).value = defaults{4};

c=c+1;
dlg(c).string = 'Recon/Readme by:';
dlg(c).fieldName = 'scribe';
dlg(c).style = 'edit';
dlg(c).value = defaults{5};

c=c+1;
dlg(c).string = 'Magnet:';
dlg(c).fieldName = 'magnet';
dlg(c).style = 'popupmenu';
dlg(c).list = {'1.5 T', '3 T'};
dlg(c).choice = 2;
dlg(c).value = defaults{6};

c=c+1;
dlg(c).string = 'Coil:';
dlg(c).fieldName = 'coil';
dlg(c).style = 'popupmenu';
dlg(c).list = {'Posthead',  'RessCoil', 'fMRICoil', 'Posthead-ReceiveOnly','Helmet-head','Small-head', 'GE-head', 'Flex-quad'};
dlg(c).choice = 1;
dlg(c).value = defaults{7};

c=c+1;
dlg(c).fieldName = 'examNumber';
dlg(c).string = 'Exam number:';
dlg(c).style = 'edit';
dlg(c).value = defaults{8};

c=c+1;
dlg(c).string = 'Slice orientation:';
dlg(c).fieldName = 'sliceOrientation';
dlg(c).style = 'popupmenu';
dlg(c).list = {'Oblique', 'Coronal', 'Axial', 'Saggital'};
dlg(c).choice = 1;
dlg(c).value = defaults{9};

c=c+1;
dlg(c).string = 'Protocol name:';
dlg(c).fieldName = 'protocol';
dlg(c).style = 'edit';
dlg(c).value = defaults{10};

% Prepare to make dialog box, make it
height = 1;
vSkip = 0.12;
pos = [35,10,80,length(dlg)*(height+vSkip)+3];
x = 1;
y = length(dlg)*(height+vSkip)+1;
editWidth = 40;
stringWidth = 30;
for uiNum = 1:length(dlg)
    dlg(uiNum).stringPos = [x,y,stringWidth,height];
    dlg(uiNum).editPos = [x+stringWidth-1,y,editWidth,height];
    y = y-(height+vSkip);
end
resp = generaldlg(dlg, pos, title);

% Write to Readme.txt
[fid, message] = fopen('Readme.txt', 'w');
if (~exist('message','var'))
    message='default message';
end

if fid == -1
    warndlg(message);
    return
end

for i = 1:length(dlg)
    fprintf(fid, '%s %s\n', dlg(i).string, resp.(dlg(i).fieldName));
end

status = fclose(fid);
if status == -1
    warndlg(messsage);
    return
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
global planC
[SUVmax,SUVmean,BKG,BSL,Thresh,equiVol] = calcBSL(1,planC,1,1);
[SUVmax,SUVmean,BKG,BSL,Thresh,equiVol] = calcBSL(2,planC,1,1);
[SUVmax,SUVmean,BKG,BSL,Thresh,equiVol] = calcBSL(3,planC,1,1);

author: cerr

repo: CERR


text:
function n = p01_n ( )

%*****************************************************************************80
%
%% P01_N returns the number of variables for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 3;

  return
end

repo: jburkardt-m