text:
% The original source code is from https://github.com/zhengliu6699/imageFusionMetrics/blob/master/metricChenBlum.m
% The interface is modified by the authors of VIFB to integrate it into VIFB.

function res=metricsQcb(img1,img2,fused)

    fused = double(fused);
    img1 = double(img1);
    img2 = double(img2);
    % Get the size of img
    [m,n,b] = size(fused);
    [m1,n1,b1] = size(img2);

    if b == 1
        g = Qcb(img1,img2,fused);
        res = g;
    elseif b1 == 1
        for k = 1 : b
            g(k) = Qcb(img1(:,:,k),img2,fused(:,:,k));
        end
        res = mean(g);
    else
        for k = 1 : b
            g(k) = Qcb(img1(:,:,k),img2(:,:,k),fused(:,:,k));
        end
        res = mean(g);
    end

end

function output = Qcb(im1, im2, fused)

    % function res=metricChenBlum(im1,im2,fused)
    %
    % This function implements Yin Chen's algorithm for fusion metric.
    % im1, im2 -- input images;
    % fused      -- fused image;
    %
    % IMPORTANT: The size of the images need to be 2X.
    % See also: evalu_fusion.m
    %
    % Z. Liu [July 2009]    %

    % Ref: A new automated quality assessment algorithm for image fusion, Image and Vision Computing, 27 (2009) 1421-1432
    % By Yin Chen et al.
    %

    im1 = im2double(im1);
    im2 = im2double(im2);
    fused = im2double(fused);

    im1=normalize1(im1);
    im2=normalize1(im2);
    fused=normalize1(fused);

    %% set up some constant values for experiment

    f0=15.3870;
    f1=1.3456;
    a=0.7622;

    % parameters for local constrast computation
    k=1;
    h=1;
    p=3; %2.4;
    %p=2.4;
    q=2;
    Z=0.0001;
    sigma=2;
    %% caculate the quality Q

    [hang,lie]=size(im1);

    %DoG filter
    %DoG1
    %HH=hang/2; LL=lie/2;
    HH=hang/30; LL=lie/30;

    %DoG2
    %HH=hang/4; LL=lie/4;

    %DoG3
    %HH=hang/8; LL=lie/8;

    [u,v]=freqspace([hang,lie],'meshgrid');
    u=LL*u; v=HH*v;
    r=sqrt(u.^2+v.^2);

    Sd=exp(-(r/f0).^2)-a*exp(-(r/f1).^2);

    % constrast sensitivity filtering
    fused1=ifft2(ifftshift(fftshift(fft2(im1)).*Sd));
    fused2=ifft2(ifftshift(fftshift(fft2(im2)).*Sd));
    ffused=ifft2(ifftshift(fftshift(fft2(fused)).*Sd));

    %--------------------
    %fused1=normalize1(fused1);
    %fused2=normalize1(fused2);
    %ffused=normalize1(ffused);

    % local contrast computation
    % one level of contrast
    G1=gaussian2d(hang,lie,2);
    G2=gaussian2d(hang,lie,4);

    % filtering in frequency domain
    C1=contrast(G1,G2,fused1);
    C1=abs(C1); % I add this. (see your notes)
    C1P=(k*(C1.^p))./(h*(C1.^q)+Z);

    C2=contrast(G1,G2,fused2);
    C2=abs(C2); % I add this.
    C2P=(k*(C2.^p))./(h*(C2.^q)+Z);

    Cf=contrast(G1,G2,ffused);
    Cf=abs(Cf); % I add this.
    CfP=(k*(Cf.^p))./(h*(Cf.^q)+Z);

    % contrast preservation calculation
    mask=(C1P<CfP);
    mask=double(mask);
    Q1F=(C1P./CfP).*mask+(CfP./C1P).*(1-mask);

    mask=(C2P<CfP);
    mask=double(mask);
    Q2F=(C2P./CfP).*mask+(CfP./C2P).*(1-mask);

    % Saliency map generation
    ramda1=(C1P.*C1P)./(C1P.*C1P+C2P.*C2P);
    ramda2=(C2P.*C2P)./(C1P.*C1P+C2P.*C2P);

    % global quality map

    Q=ramda1.*Q1F+ramda2.*Q2F;

    output=mean2(Q);
end

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% sub-functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res=gaussian2d(n1,n2,sigma)

% creat a 2D Gaussian filter in spatial domain
%

% hang (H)-> y; lie (L) -> x

H=floor((n1-1)/2);
L=floor((n2-1)/2);

[x,y]=meshgrid(-15:15,-15:15);
G=exp(-(x.*x+y.*y)/(2*sigma*sigma))/(2*pi*sigma*sigma);

%This is to normalize
%G=G/sum(G(:));
res=G;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res=contrast(G1,G2,im)

%[hang,lie]=size(im);

%FG1=fft2(G1,hang,lie);
%FG2=fft2(G2,hang,lie);
%fused=fft2(im);

%buff=real(ifft2(FG1.*fused));
%buff1=real(ifft2(FG2.*fused));

buff=filter2(G1,im,'same');
buff1=filter2(G2,im,'same');

res=buff./buff1-1;
end

function RES=normalize1(data)

    % function RES=normalize1(data)
    %
    % This function is to NORMALIZE the data.
    % The data will be in the interval 0-255 (gray level) and pixel value has
    % been rounded to an integer.
    %
    % See also: normalize.m
    %
    % Z. Liu @NRCC (Aug 24, 2009)

    data=double(data);
    da=max(data(:));
    xiao=min(data(:));
    if (da==0 & xiao==0)
        RES=data;
    else
        newdata=(data-xiao)/(da-xiao);
        RES=round(newdata*255);
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
clc;
clearvars;
close all;
rng default;

% subspace dimension
D = 4;
% ambient dimension
M = 10;
% Number of subspaces
K = 2;
import spx.data.synthetic.subspaces.random_subspaces;
bases = random_subspaces(M, K, D);
A = bases{1};
B = bases{2};
G = A' * B
sigmas = svd(G)'
largest_product = sigmas(1)
smallest_angle_rad  = acos(largest_product)
smallest_angle_deg = rad2deg(smallest_angle_rad)

K = 6;
bases = random_subspaces(M, K, D);
angles = spx.la.spaces.smallest_angles_deg(bases)
angles = spx.matrix.off_diag_upper_tri_elements(angles)'

angles = spx.la.spaces.smallest_angles_rad(bases)
angles = spx.la.spaces.smallest_angles_cos(bases)

author: indigits

repo: sparse-plex


text:
% Sigmoid Belief Hidden Markov Decision Tree    (Jordan/Gharhamani 1996)
%
clear all;
%clc;
rand('state',0); randn('state',0);
X = 1; Q1 = 2; Q2 = 3; Y = 4;
% intra time-slice graph
intra=zeros(4);
intra(X,[Q1 Q2 Y])=1;
intra(Q1,[Q2 Y])=1;
intra(Q2, Y)=1;
% inter time-slice graph
inter=zeros(4);
inter(Q1,Q1)=1;
inter(Q2,Q2)=1;

ns = [1 2 3 1];
dnodes = [2 3];
eclass1 = [1 2 3 4];
eclass2 = [1 5 6 4];
bnet = mk_dbn(intra, inter, ns, dnodes, eclass1, eclass2);

bnet.CPD{1} = root_CPD(bnet, 1);
% =========================================
bnet.CPD{2} = softmax_CPD(bnet, 2);
bnet.CPD{3} = softmax_CPD(bnet, 3, 'discrete', [2]);
bnet.CPD{5} = softmax_CPD(bnet, 6);
bnet.CPD{6} = softmax_CPD(bnet, 7, 'discrete', [3 6]);
% =========================================
bnet.CPD{4} = gaussian_CPD(bnet, 4);

% make some data
T=20;
cases = cell(4, T);
cases(1,:)=num2cell(round(rand(1,T)*2)+1);
%cases(2,:)=num2cell(round(rand(1,T))+1);
%cases(3,:)=num2cell(round(rand(1,T)*2)+1);
cases(4,:)=num2cell(rand(1,T));

engine = bk_inf_engine(bnet, 'exact', [1 2 3 4]);

% log lik before learning
[engine, loglik] = enter_evidence(engine, cases);

% do learning
ev=cell(1,1);
ev{1}=cases;
[bnet2, LL2] = learn_params_dbn_em(engine, ev, 10);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [spec, headerinfo] = read_caret_spec(specfile)

% READ_CARET_SPEC reads in a caret .spec file.
%
% Use as
%   [spec, headerinfo] = read_caret_spec(specfile)
%
% Output arguments:
%   spec       = structure containing per file type the files listed
%   headerinfo = structure containing the specfile header
%
% The file can be an xml-file or an ascii formatted file

% Copyright (C) 2013, Jan-Mathijs Schoffelen

try,
  % read xml-file that contains a description to a bunch of files
  % belonging together
  ft_hastoolbox('gifti', 1);
  g = xmltree(specfile);

  % convert into a structure
  s = convert(g);

  if isfield(s, 'FileHeader')
    headerinfo = s.FileHeader;
    spec       = rmfield(s, 'FileHeader');
  else
    headerinfo = [];
    spec       = s;
  end

  % process the headerinfo
  if ~isempty(headerinfo)
    if isfield(headerinfo, 'Element')
      tmp = headerinfo.Element;
      tmp2 = struct([]);
      for k = 1:numel(headerinfo.Element)
        tmp2(1).(strrep(headerinfo.Element{k}.Name, '-', '_')) = headerinfo.Element{k}.Value;
      end
      headerinfo = tmp2;
    end
  end

  % further process the fields in spec
  f = fieldnames(spec);

  for k = 1:numel(f)
    if isempty(strfind(f{k}, 'study_metadata'))
      if iscell(spec.(f{k}))
        tmp = spec.(f{k});
        tmp2 = {};
        for m = 1:numel(tmp)
          tmpx = tmp{m};
          if isstruct(tmpx)
            fn = fieldnames(tmpx)
            for i = 1:numel(fn)
              tmp2{end+1,1} = tmpx.(fn{i});
            end
          end
        end
        spec.(f{k}) = tmp2;
      elseif isstruct(spec.(f{k}))
        tmp = spec.(f{k});
        fn  = fieldnames(tmp);
        tmp2 = {};
        for m = 1:numel(fn)
          tmp2{end+1,1} = tmp.(fn{m});
        end
        spec.(f{k}) = tmp2;
      else
        % don't know what to do with it
        spec = rmfield(spec, f{k});
      end
    else
      % don't know what to do with it
      spec = rmfield(spec, f{k});
    end
  end

catch

  % process as ASCII-file
  fid  = fopen_or_error(specfile);
  line = 'some text';
  while isempty(strfind(line, 'EndHeader'))
    line = fgetl(fid);
    if isempty(strfind(line, 'BeginHeader')) && isempty(strfind(line, 'EndHeader'))
      tok = tokenize(line, ' ');
      headerinfo.(strrep(tok{1},'-','_')) = tok{2};
    end
  end
  line = fgetl(fid); % empty line

  spec = struct([]);
  while 1
    line = fgetl(fid);
    if ~ischar(line), break, end
    tok = tokenize(line, ' ');
    if ~isempty(tok{1})
    if isfield(spec, tok{1})
      spec(1).(tok{1}){end+1,1} = tok{2};
    else
      spec(1).(tok{1}){1} = tok{2};
    end
    end
  end
  fclose(fid);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function restorefig(h,old)
%RESTOREFIG  Restore a figure's properties
%   RESTOREFIG(H,OLD) restores the properties of H specified in
%   OLD. The state-difference structure OLD is the output of the
%   function APPLYTOFIG. Deleting object handles after the call to
%   APPLYFIG might cause undesired behavior.
%
%   See also EXPORTFIG, PREVIEWFIG, APPLYTOFIG.

%  Copyright 2000-2009 The MathWorks, Inc.

for n=1:length(old.objs)
  try
    if ~iscell(old.values{n}) & iscell(old.prop{n})
      old.values{n} = {old.values{n}};
    end
    set(old.objs{n}, old.prop{n}, old.values{n});
  end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [Va, Mask] = reslice_vol(Vo, M, interp)
% function [Va, Mask] = reslice_vol(Vo, M, interp)

% Ripped out of SPM 8 and modified (2010, S.Klanke)

dim = size(Vo);
wrap = [1;1;0];
d = [interp*[1;1;1] wrap];

[x1,x2] = ndgrid(1:dim(1), 1:dim(2));
C = spm_bsplinc(Vo, d);
Va = zeros(dim);
Mask = zeros(dim);

for x3 = 1:dim(3)
	[msk_x3,y1,y2,y3] = getmask(M,x1,x2,x3,dim,wrap);
	Mask(:,:,x3) = msk_x3;
    Va(:,:,x3) = spm_bsplins(C, y1,y2,y3, d);
end
return;

function [Mask,y1,y2,y3] = getmask(M,x1,x2,x3,dim,wrp)
tiny = 5e-2; % From spm_vol_utils.c
y1   = M(1,1)*x1+M(1,2)*x2+(M(1,3)*x3+M(1,4));
y2   = M(2,1)*x1+M(2,2)*x2+(M(2,3)*x3+M(2,4));
y3   = M(3,1)*x1+M(3,2)*x2+(M(3,3)*x3+M(3,4));
Mask = true(size(y1));
if ~wrp(1), Mask = Mask & (y1 >= (1-tiny) & y1 <= (dim(1)+tiny)); end
if ~wrp(2), Mask = Mask & (y2 >= (1-tiny) & y2 <= (dim(2)+tiny)); end
if ~wrp(3), Mask = Mask & (y3 >= (1-tiny) & y3 <= (dim(3)+tiny)); end
return;


text:
function volume = ball01_volume ( )

%*****************************************************************************80
%
%% BALL01_VOLUME returns the volume of the unit ball.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VOLUME, the volume of the unit ball.
%
  r = 1.0;
  volume = 4.0 * pi * r ^ 3 / 3.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function obj = changeFont(obj,varargin)
	% Change figure font size
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% DESCRIPTION
	options.fontSize = [];
	% get options
	options = getOptions(options,varargin);
	% disp(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	if isempty(options.fontSize)
		userInput = inputdlg('New font');
		userInput = str2num(userInput{1});
		set(findall(gcf,'-property','FontSize'),'FontSize',userInput);
	else
		set(findall(gcf,'-property','FontSize'),'FontSize',options.fontSize);
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% In the field of robotics there are many possible ways of representing
% orientations of which the most common are:
% - orthonormal rotation matrices (3x3),
% - three angles (1x3), and
% - quaternions.

% A rotation of pi/2 about the x-axis can be represented as an orthonormal rotation
% matrix

R = rotx(pi/2)
% which we can see is a 3x3 matrix.

% Such a matrix has the property that it's columns (and rows) are sets of orthogonal
% unit vectors.  The determinant of such a matrix is always 1

det(R)

% Let's create a more complex rotation

R = rotx(30, 'deg') * roty(50, 'deg') * rotz(10, 'deg')
% where this time we have specified the rotation angle in degrees.

% Any rotation can be expressed in terms of a single rotation about some axis
% in space

[theta,vec] = tr2angvec(R)
% where theta is the angle (in radians) and vec is unit vector representing the
% direction of the rotation axis.

% Commonly rotations are represented by Euler angles

eul = tr2eul(R)
% which are three angles such that R = rotz(a)*roty(b)*rotz(c), ie. the rotations
% required about the Z, then then the Y, then the Z axis.

% Rotations are also commonly represented by roll-pitch-yaw  angles

rpy = tr2rpy(R)
% which are three angles such that R = rotx(r)*roty(p)*rotz(y), ie. the rotations
% required about the X, then then the Y, then the Z axis.

% We can investigate the effects of rotations about different axes
% using this GUI based demonstration.  The menu buttons allow the rotation
% axes to be varied
%  *** close the window when you are done.
tripleangle('rpy', 'wait')

% The final useful form is the quaternion which comprises 4 numbers.  We can create
% a quaternion from an orthonormal matrix

q = Quaternion(R)
% where we can see that it comprises a scalar part and a vector part.  To convert back

q.R
% which is the same of the value of R we determined above.

% Quaternions are a class and the orientations they represent can be compounded, just
% as we do with rotation matrices by multiplication.

% First we create two quaternions

q1 = Quaternion( rotx(pi/2) )
q2 = Quaternion( roty(pi/2) )

% then the rotation of q1 followed by q2 is simply

q1 * q2

% We can also take the inverse of a Quaternion

q1 * inv(q1)
% which results in a null rotation (zero vector part)

author: Allopart

repo: rbpf-gmapping


text:
function [Yest, results] = lle(Y, ssub, rr, ACTIVE_PX)
%% approximate the background with locally-linear embedding
%% inputs:
%   Y: d1*d2*T 3D matrix, video data
%   rr: scalar, average neuron size
%   ssub: spatial downsampling factor
%   ACTIVE_PX:  indicators of pixels to be approximated

%% outputs:
%   Yest: d1*d2*T 3D matrix, reconstructed video data
%   results: struct variable {weights, ssub}
%       weights: d1*d2 cell, each element is a 2*J matrix. Row 1 has the indice of the
%       ring neighbors and row 2 has the corresponding weights.
%       ssub:    scalar, spatial downsampling factor

%% Author: Pengcheng Zhou, Carnegie Mellon University,2016

%% input arguments
[d1, d2, T] = size(Y);

% center the fluorescence intensity by its mean
Ymean = mean(Y, 3);
Y = Y - bsxfun(@minus, Ymean, ones(1, 1, T));

% average neuron size
if ~exist('rr', 'var')|| isempty(rr)
    rr = 15;
end
% spatial downsampling
if ~exist('ssub', 'var') || isempty(ssub)
    ssub = 1;
end

%downsample the data
if ssub>1
    Y = imresize(Y, 1./ssub);
    [d1s, d2s, ~] = size(Y);
    rr = round(rr/ssub)+1;
else
    d1s = d1;
    d2s = d2;
end

% pixels to be approximated
if exist('ACTIVE_PX', 'var') && ~isempty(ACTIVE_PX)
    ACTIVE_PX = reshape(double(ACTIVE_PX), d1, d2);
    ACTIVE_PX = (imresize(ACTIVE_PX, 1/ssub)>0);
    ind_px = find(ACTIVE_PX(:));  % pixels to be approxiamted
else
    ind_px = (1:(d1s*d2s));
end

%% determine neibours of each pixel
c_shift = (-rr:rr);
r_shift = round(sqrt(rr^2-c_shift.^2));
c_shift = [c_shift, c_shift(2:end)];
r_shift = [-r_shift, r_shift(2:end)];

[csub, rsub] = meshgrid(1:d2s, 1:d1s);
csub = reshape(csub, [], 1);
rsub = reshape(rsub, [], 1);
csub = bsxfun(@plus, csub, c_shift);
rsub = bsxfun(@plus, rsub, r_shift);
% remove neighbors that are out of boundary
ind = or(or(csub<1, csub>d2s), or(rsub<1, rsub>d1s));
csub(ind) = nan;
rsub(ind) = nan;
% options = optimoptions('lsqlin','Algorithm','active-set', 'Display', 'none');

%% run approximation
 gamma = 0.001; % add regularization
Y = reshape(Y, d1s*d2s, []);
Yest = zeros(size(Y));
weights = cell(d1s, d2s);
for m=1:length(ind_px)
    px = ind_px(m);
    ind_nhood = sub2ind([d1s,d2s], rsub(px, :), csub(px, :));
    ind_nhood(isnan(ind_nhood)) = [];
    J = length(ind_nhood);

    G = bsxfun(@minus, Y(ind_nhood, 1:3:end), Y(px, 1:3:end));
    w = (G*G'+gamma*sum(G(:).^2)*eye(J))\ones(length(ind_nhood), 1);
    w = w/sum(w);
%     w = lsqlin(Y(ind_nhood,:)', Y(px,:)', [], [], ones(1, J), ...
%         1, ones(J, 1)*0, ones(J, 1)*3/J, [], options); % weights have box constraints [min_w, max_w]
    Yest(px, :) = w'*Y(ind_nhood, :);
    weights{px} = [ind_nhood; w'];
end
results.weights = weights;
results.ssub = ssub;

ind = 1:(d1s*d2s);
ind(ind_px) = [];
if ~isempty(ind)
    temp = imfilter(Y, ones(3, 3)/9, 'replicate');
    Yest(ind, :) = temp(ind, :); % without approximation
end
Yest = reshape(Yest, d1s, d2s, []);

    %% return the result
if ssub>1 %up sampling
    Yest = imresize(Yest, [d1, d2]);
end

clear Y;
Ybaseline = Ymean - median(Yest, 3);
Yest = bsxfun(@plus, Yest, Ybaseline);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function r = eq(p,q)
% this function compares the equality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

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

if (isa(p, 'CellVariable')&&isa(q, 'CellVariable'))
    r=p;
    r.value = p.value==q.value;
elseif isa(p, 'CellVariable')
    r=p;
    r.value = p.value==q;
else
    r=q;
    r.value = p==q.value;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function DataSet = prtDataGenNoisyLine(varargin)
% prtDataGenNoisyLine Generates noisy line example data
%
%   DATASET = prtDataGenNoisyLine returns a prtDataSetRegress with 100
%   samples of a line on x = [-1, 1] with zero-mean additive Gaussian
%   noise. By default, the slope is 1 and the y-intercept is 0. The noise
%   variance is .1.
%
%   DATASET = prtDataGenNoisyLine(param,val) enables setting the following
%   parameters:
%       slope: 1
%       yIntercept: 0
%       xRange: [-1, 1]
%       xRange: [-1, 1]
%       nSamples: 1000
%       stdev: 0.1
%
%   Example:
%
%   ds = prtDataGenNoisyLine;
%   plot(ds)
%
%   See also: prtDataSetRegress, prtDataGenNoisySinc

p = inputParser;
p.addParameter('slope',1);
p.addParameter('yIntercept',0);
p.addParameter('xRange',[-1 1]);
p.addParameter('nSamples',100);
p.addParameter('stdev',0.1);

p.parse(varargin{:});

x = rand(p.Results.nSamples,1)*range(p.Results.xRange)+p.Results.xRange(1);
t = x*p.Results.slope + p.Results.yIntercept;
y = t + p.Results.stdev*randn(size(x));

DataSet = prtDataSetRegress(x,y,'name','Noisy Line');

% function x = prtUtilRandSample(vector,nSamples)
%
% ind = randperm(length(vector));
% ind = ind(1:nSamples);
% x = vector(sort(ind));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function mv = mv_reorderCondNums(mv,newOrder);
%
% mv = mv_reorderCondNums(mv,[newOrder]);
%
% Specify a new order for non-null conditions.
% If not specified, will pop up a UI.
%
%
%
% ras 03/05.
if notDefined('mv'),    mv = get(gcf,'UserData');		end
if ishandle(mv),		mv = get(mv, 'UserData');		end

if notDefined('newOrder')
    % put up a dialog
    nums = mv.trials.condNums;
    names = mv.trials.condNames;
    def = {num2str(nums)};
    prompt{1} = 'New Condition # Order [Current Order: ';
    for i = 1:length(names)
        prompt{1} = [prompt{1} num2str(nums(i)) ') ' names{i} ' '];
    end
    prompt{1} = [prompt{1} ']'];
	dlgTitle='Reorder Conditions';
	answer=inputdlg(prompt,dlgTitle,1,def);
    newOrder = str2num(answer{1});
end

% if the null isn't included, put it first:
if ~ismember(newOrder, 0)
	newOrder = [0 newOrder];
end

%% need to actually change the condition numbers in the event
%% specification...tricky...
% first, make all the conditions negative. We don't allow negative
% condition specification, so these won't be confused in the shuffling
mv.trials.cond = -(mv.trials.cond + 1); % null is -1, down from there

% now, map each original value to the new value
for c = 1:length(newOrder)
	I = find(mv.trials.cond == -c);
	mv.trials.cond(I) = newOrder(c);
end

% get the indexing needed to produce this
% new order:
% [ignore ind] = sort(newOrder);
newOrder = newOrder + 1; % just account for null

% now change the labels, colors emv.
mv.trials.condNums = mv.trials.condNums(newOrder);
mv.trials.condNames = mv.trials.condNames(newOrder);
mv.trials.condColors = mv.trials.condColors(newOrder);

set(gcf,'UserData',mv);

multiVoxelUI;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function p = circle_imp_points_2d ( r, center, n )

%*****************************************************************************80
%
%% CIRCLE_IMP_POINTS_2D returns N equally spaced points on an implicit circle in 2D.
%
%  Discussion:
%
%    The first point is always ( CENTER(1) + R, CENTER(2) ), and subsequent
%    points proceed counterclockwise around the circle.
%
%    An implicit circle in 2D satisfies the equation:
%
%      ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 = R^2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, integer N, the number of points desired.  N must be at least 1.
%
%    Output, real P(2,N), the coordinates of points
%    on the circle.
%
  for i = 1 : n
    theta = ( 2.0 * pi * ( i - 1 ) ) / n;
    p(1,i) = center(1,1) + r * cos ( theta );
    p(2,i) = center(2,1) + r * sin ( theta );
  end

  return
end


text:
function [KMZFullFile] = kml2kmz(varargin)
%kml2kmz converts .kml file types into self-contained .kmz filetypes.
%   kml2kmz(KMLFullFile) copies all the support documents used by the kml
%   into a single 'files' folder, and updates the file path in the kml to
%   the new relative path accordingly.  The kml is also copied, then zipped
%   up with the new copies of supporting files.  The zip file is then
%   renamed with a .kmz extension, and put in the same folder as the
%   original kml.  The new copies of supporting documents are then deleted
%   by default.  Supporting documents referenced by KML attributes are not yet
%   supported.
%
%   kml2kmz(KMLFullFile,'keepsupportfiles') keeps the extra copies of the
%   supporting documents.
%
%   kml2kmz(KMLFullFile, 'deleteoriginalkml') deletes the original kml.
%
%   Both 'keepsupportfiles' and 'deleteoriginalkml' can be used at the same time.
%
%%%%%%%%%%%%%%%%%%%%
%
%   Still to be Considered: Please contact me if you can help me solve these
%   questions/points.
%       1:  Currently checks for addresses in kml elements, but we need to deal with
%       the kml attribute case too.  I don't have an example to work with.
%       2:  Haven't checked fucntionality on Linux or Mac, but should work
%       as far as I know
%       3:  Haven't checked if ".." notation is handled properly
%       4:  I haven't worked with Alias' before.  Are Targethref and
%       Sourcehref handled properly?
%       5:  Are there other element types I've missed?
%
%
%   Author(s): Ryan Bell, ryan.bell@viu.ca
%   Vancouver Island University
%   Version 2.1, December 05, 2012
%
%
%   version 1.0  - intial coding
%   version 2.0  - added internet file download
%                - added check to prevent downloading same file twice
%                - added check to ensure file names don't conflict
%                - removed platform check
%                - added <a href> to catch description text
%                - improved description text and added a todo list
%   version 2.1  - added conditional check for existing, but empty, elements

%% Setup Inputs

%     % Defaults;
    keep_support = 0;
    delete_original = 0;
    KMLFullFile = varargin{1};
    % holds a list of files accessed, so we don't copy the same one more than once
    FileList={'source files','destination files'};
    nin = length(varargin);

%% check for incorrect arguments
    % if the input args is more than 1 - it should be either
    % 'keeplimits' or 'keepticks' or both.
    if nin > 1
        for i = nin:-1:2
            if ~(strcmpi(varargin{i},'keepsupportfiles') || strcmpi(varargin{i},'deleteoriginalkml'))
                error('Incorrect input arguments');
            end
        end
    end

    % Look for 'keeplimits'
    for i=2:nin
        if strcmpi(varargin{i},'keepsupportfiles')
            keep_support = 1;
        elseif strcmpi(varargin{i},'deleteoriginalkml')
            delete_original = 1;
        end
    end

%% Setup support file/folder names
   % Decide where to keep the temporary support files
    if keep_support
        %put the support files where they can be found
        parentfolder = fileparts(KMLFullFile);
    else
        %put the temporary support files where they won't be noticed
        parentfolder = tempname;
    end

    %location of the temporary support documents
    filesPath = fullfile(parentfolder,'files');
    %location of the temporary kml
    tempKMLFullFile = fullfile(parentfolder,'doc.kml');
    %location of the new kmz
    KMZFullFile = [KMLFullFile(1:end-3) 'kmz'];
    %make sure folder exists
    if ~isdir(filesPath)
        mkdir(filesPath);
    end

%% Open file and look for references to local and internet addresses

    %open the original kml as an xDocument
    xDoc = xmlread(KMLFullFile);

    % Go through the kml and find all the referenced files, copy them in to
    % support folder, and update the file names in the kml to relative paths

    % checks for address in elements, but
    for tag = {'href','targetHref','sourceHref','description'}  %others?  https://developers.google.com/kml/documentation/kmzarchives
        allListitems = xDoc.getElementsByTagName(tag);  %find all the references to external files
        for k = 0:allListitems.getLength-1 %go through each refence (Java index starts at 0)
            thisListitem = allListitems.item(k);
            %make sure xml element isn't empty
            if ~isempty(thisListitem.getFirstChild)
                Contents = char(thisListitem.getFirstChild.getData);  %this is the location of support file
                if strcmp(tag,'description')
                    Description = Contents; %Contents are the description html text
                    [Description, FileList] = handle_description(Description, filesPath,FileList);
                    thisListitem.getFirstChild.setData(Description);
                else
                    SourceFile = Contents; %Contents are the source file
                    [DestinationFile_GE, FileList] = handle_link(SourceFile, filesPath,FileList);
                    thisListitem.getFirstChild.setData(DestinationFile_GE);
                end
            end
        end
    end

%% Save the work and Tidy up

    % Rewrite the kml with its new file paths
    write_XML_no_extra_lines(tempKMLFullFile,xDoc) %bug fix see:http://www.mathworks.com/matlabcentral/newsreader/view_thread/245555
    % xmlwrite(tempKMLFullFile,xDoc);  %write the updated xml file into out temp directory

    % This creates a zipped file named 'myKmz.kmz.zip'
    zip(KMZFullFile,{filesPath tempKMLFullFile});
    %renames the zipped file to 'mykmz.kmz'.
    movefile([KMZFullFile '.zip'], KMZFullFile);

    if keep_support == 0
        % Delete supporting files
        delete(tempKMLFullFile) %Delete the temp kml
        rmdir(filesPath,'s');  %delete the files folder
    end
    if delete_original == 1
        delete(KMLFullFile);  %delete original kml
    end

end

%% Functions...

%% Check description text for instances of images of
function [Description, FileList] = handle_description(Description, filesPath,FileList)
    addresses = [strfind(Description,'img src'),...
                 strfind(Description,'a href')];
    for i = addresses
        j = strfind(Description(i:end),'"')+i;
        SourceFile = Description(j(1):j(2)-2);
        DestinationFile_GE = handle_link(SourceFile, filesPath,FileList);
        Description = [Description(1:j(1)-1) DestinationFile_GE Description(j(2)-1:end)];
    end
end

% Move support file into 'files' folder and convert address relative
function [DestinationFile_GE, FileList] = handle_link(SourceFile, filesPath,FileList)

    if ~any(strcmp(SourceFile,FileList(:,1)))  %Check if we've copied or download this file already
        [~,FileName,FileExt] = fileparts(SourceFile);
        DestinationFile_GE = ['files/' FileName FileExt]; %relative GE filepath, using GE file separator (/)
        n=1;
        while any(strcmp(DestinationFile_GE,FileList(:,2)))  %Check if this path is already used (eg the file shares a file name, but was in a different folder from a previous support file)
            FileNameOld = FileName;
            FileName = [FileNameOld ' ' num2str(n,'%.0f')];
            DestinationFile_GE = ['files/' FileName FileExt]; %add a numeric extension to this file that shared the same name as a previous support file
            n=n+1;
        end
        DestinationFile_Local = fullfile(filesPath,[FileName FileExt]);  %destination file path using local file separator

        if strcmp(SourceFile(1:4),'http') || strcmp(SourceFile(1:3),'ftp')  %file is on the internet
            [~,status] = urlwrite(SourceFile,DestinationFile_Local); %download the file this can be painfully slow for some reason...
            if status == 0
                error(['Problem downloading ' SourceFile]);
            end

        else  %file is local
            if ~strcmp(SourceFile(2),':')  %file path relative, add the current folder to the filepath
                SourceFile = fullfile(currentFolder, SourceFile); %change file path from relative to absolute
            end

            if exist(SourceFile,'file') %check if file to be copied exists
                status = copyfile(SourceFile,DestinationFile_Local);  % copyfile into the 'files' folders
                if status == 0
                    error(['Problem copying ' SourceFile]);
                end
            else
                error(['Can not find file:' SourceFile]);
            end
        end
        FileList{size(FileList(:,1),2)+1,1} = SourceFile;                %Store the source file for later checking
        FileList{size(FileList(:,1),2)+1,2} = DestinationFile_GE;        %Store the GE relative path for possible, this path may be reused later
    else
        DestinationFile_GE = FileList{strcmp(SourceFile,FileList(:,1)),2}; % This file is being used again, get the GE relative path
    end
%     save('kmldata')  %for debug
end

%% bug fix.  See: http://www.mathworks.com/matlabcentral/newsreader/view_thread/245555
function write_XML_no_extra_lines(save_name,XDOC)

    XML_string = xmlwrite(XDOC); %open KML as string
    XML_string = regexprep(XML_string,'\n[ \t\n]*\n','\n'); %removes extra tabs, spaces and extra lines

    %Write to file
    fid = fopen(save_name,'w');
    fprintf(fid,'%s\n',XML_string);
    fclose(fid);
end

author: Sable

repo: mcbench-benchmarks


text:
function tagS = code_sequence_macro_tags
%"code_sequence_macro_tags"
%   Returns the tags associated with a code sequence, specified by section
%   8.8 in PS3.3 of 2006 DICOM.
%
%WARNING: The types of the tags specified in this function DO NOT MATCH the
%DICOM specification, since Type 1C and 2C elements whose only condition
%for existance is that the sequence is present have already had that
%requirement implicitly met; hence they are changed to Type 1 and 2
%respectively.
%
%JRA 06/06/06
%
%Usage:
%   tagS = code_sequence_macro_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% Define tags
tagsC = {'00080100','00080102','00080103','00080104','0008010F','00080105',...
    '00080106','0008010B','00080107','0008010D'};
tagDecC = {524544, 524546, 524547, 524548, 524559, 524549, 524550,...
    524555, 524551, 524557};
% Define tag type
tagTypeC = {'1','1','1C','1','3','1C','1C','3','1C','1C'};

% Define children
childC = {[],[],[],[],[],[],[],[],[],[]};

tagS = struct('tag', tagsC, 'tagdec', tagDecC, 'type', tagTypeC, 'children', childC);

% %Initialize the tagS structure.
% tagS = struct('tag', {}, 'type', {}, 'children', {]);
%
% %Add tags based on PS3.3 attribute lists.
%
% %Code Value
% tagS(end+1) = struct('tag', ['00080100'], 'type', ['1'], 'children', []);
%
% %Coding Scheme Designator
% tagS(end+1) = struct('tag', ['00080102'], 'type', ['1'], 'children', []);
%
% %Coding Scheme Version
% tagS(end+1) = struct('tag', ['00080103'], 'type', ['1C'], 'children', []);
%
% %Code Meaning
% tagS(end+1) = struct('tag', ['00080104'], 'type', ['1'], 'children', []);
%
% %Context Identifer
% tagS(end+1) = struct('tag', ['0008010F'], 'type', ['3'], 'children', []);
%
% %Mapping Resource
% tagS(end+1) = struct('tag', ['00080105'], 'type', ['1C'], 'children', []);
%
% %Context Group Version
% tagS(end+1) = struct('tag', ['00080106'], 'type', ['1C'], 'children', []);
%
% %Context Group Extension Flag
% tagS(end+1) = struct('tag', ['0008010B'], 'type', ['3'], 'children', []);
%
% %Context Group Local Version
% tagS(end+1) = struct('tag', ['00080107'], 'type', ['1C'], 'children', []);
%
% %Context Group Extension Creator UID
% tagS(end+1) = struct('tag', ['0008010D'], 'type', ['1C'], 'children', []);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [H] = scale(f)

% SCALE returns the homogenous coordinate transformation matrix
% corresponding to a scaling along the x, y and z-axis
%
% Use as
%   [H] = translate(S)
% where
%   S       [sx, sy, sz] scaling along each of the axes
%   H   corresponding homogenous transformation matrix

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

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

if numel(f)~=3
  ft_error('incorrect input vector');
end

H = [
  f(1) 0    0    0
  0    f(2) 0    0
  0    0    f(3) 0
  0    0    0    1
  ];

author: spm


text:
function hash = string2hash(str, len, type)
% This function generates a hash value from a text string
%
% hash = string2hash(str,type);
%
% inputs,
%   str : The text string, or array with text strings.
% outputs,
%   hash : The hash value, integer value between 0 and 2^32-1
%   type : Type of has 'djb2' (default) or 'sdbm'
%
% From c-code on : http://www.cse.yorku.ca/~oz/hash.html
%
% djb2
%  this algorithm was first reported by dan bernstein many years ago
%  in comp.lang.c
%
% sdbm
%  this algorithm was created for sdbm (a public-domain reimplementation of
%  ndbm) database library. it was found to do well in scrambling bits,
%  causing better distribution of the keys and fewer splits. it also happens
%  to be a good general hashing function with good distribution.
%
% example,
%
%  hash=string2hash('hello world');
%  disp(hash);
%
% Function is written by D.Kroon University of Twente (June 2010)
% From string to double array
str = double(str);
if ~exist('len','var')
    len = [];
end
if ~exist('type','var')
    type = 'djb2';
end
N = size(str, 2);
k = 0;
switch(type)
    case 'djb2'
        hash = 5381*ones(size(str,1),1);
        for i = 1:N
            hash = mod(hash * 33 + str(:,i), 2^32-1);
            k = k + str(:,i)*2^i;
        end
    case 'sdbm'
        hash = zeros(size(str,1),1);
        for i = 1:N
            hash = mod(hash * 65599 + str(:,i), 2^32-1);
            k = k + str(:,i)*2^i;
        end
    otherwise
        error('string_hash:inputs','unknown type');
end

hash = hash+k;
%fprintf('Initial hash:  %d\n', hash)

if ~isempty(len) && len<6
    N = round(log10(hash));
    if N <= len
        return;
    end

    d = isolateDigits(hash);
    hash = 0;
    for ii = 1:len
        hash = hash + d(ii)*10^(ii-1);
    end
end

% --------------------------------------------------------------
function digits = isolateDigits(x, base)
digits = [];

if ~exist('base','var') || (base < 2 || base > 16)
    base = 10;
end
if base ~= uint32(base)
    return;
end
if x == 0
    digits = 0;
    return;
end

kk = 1;
while 1

    d = mod(x,base);
    x = floor(x/base);

    if (x == 0) && (d == 0)
        break;
    end

    digits(kk) = d;
    kk = kk+1;

end

author: BUNPC


text:
function [vs, output, v0] = computeRatioConfidenceInterval(v0, expdata, model, max_score, ratio)

% v0 - set of flux vectors to be used as initial guesses.  They may be
% valid or not.
% expdata - experimental data.
% model - The standard model.  Additional field .N (= null(S)) should also
% be provided.  This is a basis of the flux space.
% max score - maximum allowable data fit error.
% ratio - which reactions to take the ratio of.  positive values indicate
% reactions in the numerator and negative values, reactions in the
% denominator.

majorIterationLimit = 2000;  %max number of iterations
minorIterationLimit = 1e7;   % essentially infinity
diffInterval = 1e-5;         %gradient step size.
feasibilityTolerance = max_score/20; % how close you need to be to the max score.

x0 = model.N\v0; % back substitute
numpoints = size(v0,2);
numratios = size(ratio,2);

scores = zeros(numpoints,1);

tProb.user.expdata = expdata;
tProb.user.model = model;
for i = 1:numpoints
    scores(i) = errorComputation2(x0(:,i),tProb);
end

v0(:,scores > max_score) = fitC13Data(v0(:,scores > max_score),expdata,model);

if ~isfield(model, 'N')
   model.N = null(model.S);
end

x0 = model.N\v0; % back substitute

% safety check:
if (max(abs(model.S*v0))> 1e-6)
    display('v0 not quite in null space');
    pause;
end
if(max(abs(model.N*x0 - v0)) > 1e-6)
    display('null basis is weird');
    pause;
end

nalpha = size(model.N, 2);
x_L = -1000*ones(nalpha,1);
x_U = 1000*ones(nalpha,1);
Name = 't2';
[A, b_L, b_U] = defineLinearConstraints(model);

scores = zeros(numpoints,1);
% compute scores for all points.
for i = 1:numpoints
    scores(i) = errorComputation2(x0(:,i),tProb);
end
valid_index = scores < max_score + feasibilityTolerance;
fprintf('found %d valid points\n', sum(valid_index));

c = 'errorComputation2';
c_L = 0; c_U = max_score;

dc = 'errorComputation2_grad';
d2c = []; ConsPattern = [];
%pSepFunc = [];
x_min = []; x_max = []; f_opt = []; x_opt = [];
H = [];  HessPattern = []; pSepFunc = [];  fLowBnd = [];
Solver = 'snopt';

outputminv = 222*ones(numratios, numpoints);
outputminexitflag = -222*ones(numratios, numpoints);
outputminfinalscore = -222*ones(numratios, numpoints);
outputmaxv = -222*ones(numratios, numpoints);
outputmaxexitflag = -222*ones(numratios, numpoints);
outputmaxfinalscore = -222*ones(numratios, numpoints);
outputminstruct = cell(numratios, numpoints);
outputmaxstruct = cell(numratios, numpoints);

%iterate through directions
for i = 1:numratios
%for i = 10:200
    for j = -1:2:1 % max and min
        ration = zeros(size(objective_coefficient(1,model)));
        ratiod = zeros(size(objective_coefficient(1,model)));
        for m = 1:size(ratio,1);
            if(ratio(m,i) > 0)
                ration = ration + objective_coefficient(m,model)*j;
            elseif(ratio(m,i) < 0)
                ratiod = ratiod + objective_coefficient(m,model);
            end
        end

        % initialize temp variables to make parfor work.
        v = j*222*ones(1,numpoints);
        exitflag = -222*ones(1,numpoints);
        finalscore = 222*ones(1,numpoints);
        ostruct = cell(1,numpoints);

        for k = 1:numpoints
            if exist('ttt.txt', 'file')
                fprintf('quitting due to file found\n');
                continue;
            end
            xinitial = x0(:,k);
%                 Prob  = lpconAssign(d, x_L, x_U, Name, xinitial,...
%                                       A, b_L, b_U,...
%                                       c, dc, d2c, ConsPattern, c_L, c_U,...
%                                       fLowBnd, x_min, x_max, f_opt,
%                                       x_opt);

            Prob = conAssign('ratioScore', 'ratioScore_grad', H, HessPattern, x_L, x_U, Name, xinitial, ...
                            pSepFunc, fLowBnd, ...
                            A, b_L, b_U, c, dc, d2c, ConsPattern, c_L, c_U, ...
                            x_min, x_max, f_opt, x_opt);
            %Prob.NumDiff = 2; % central diff
            %Prob.optParam.CentralDiff = 1e-5;
            %pause;
            Prob.user.expdata = expdata;
            Prob.user.model = model;
            Prob.user.ration = ration;
            Prob.user.ratiod = ratiod;

            Prob.user.max_error = max_score;
            Prob.user.diff_interval = diffInterval;
            Prob.user.useparfor = true;

            Prob.optParam.IterPrint = 0;
            %Prob.optParam.cTol = .1*feasibilityTolerance;

            Prob.PriLevOpt = 1;
            Prob.SOL.PrintFile = strcat('temp/snoptp', num2str(i), 'x', num2str(j), 'x', num2str(k),'.txt');
            Prob.SOL.SummFile = strcat('temp/snopts', num2str(i), 'x', num2str(j), 'x', num2str(k),'.txt');
            Prob.SOL.optPar(35) = majorIterationLimit; %This is major iteration count.
            Prob.SOL.optPar(30) = minorIterationLimit; %total iteration limit;
            %Prob.SOL.optPar(11) = feasibilityTolerance; % feasibility tolerance

            Result = tomRun(Solver, Prob, 5);
            tscore = errorComputation2(Result.x_k, tProb);
            tbest = Result.f_k;

            fprintf('reaction %d (%d), x %d; x=%f (%f); score=%f (%f)\n', i,length(model.lb),j, tbest,fLowBnd, tscore, max_score)

            v(k) = j*tbest;
            exitflag(k) = Result.Inform;
            finalscore(k) = tscore;
            ostruct{k} = Result;
        end
        if (j > 0) %minimizing
            outputminv(i,:) = v; % multiply by j to correct sign.
            outputminexitflag(i,:) = exitflag;
            outputminfinalscore(i,:) = finalscore;
            outputminstruct(i,:) = ostruct;
        else
            outputmaxv(i,:) = v; % multiply by j to correct sign.
            outputmaxexitflag(i,:) = exitflag;
            outputmaxfinalscore(i,:) = finalscore;
            outputmaxstruct(i,:) = ostruct;
        end
    end
end

output.minv = outputminv;
output.maxv = outputmaxv;
output.minexitflag = outputminexitflag;
output.maxexitflag = outputmaxexitflag;
output.minfinalscore = outputminfinalscore;
output.maxfinalscore = outputmaxfinalscore;
output.minstruct = outputminstruct;
output.maxstruct = outputmaxstruct;

vs = zeros(numratios, 2);
for i = 1:numratios
    validindex = output.minfinalscore(i,:) < max_score + feasibilityTolerance;
    if any(validindex)
        vs(i,1) = min(output.minv(i,validindex));
    else
        vs(i,1) = 222;
    end
    validindex = output.maxfinalscore(i,:) < max_score + feasibilityTolerance;
    if any(validindex)
        vs(i,2) = max(output.maxv(i,validindex));
    else
        vs(i,2) = -222;
    end
end

return;

% function that returns the proper objective coefficient for each reaction
% takes into account the reversibility of reactinos etc.
function [d] =  objective_coefficient(i, model)
d = zeros(length(model.lb),1);
d(i) = 1;
if (model.match(i))
    d(model.match(i)) = -1;
end
d = (d'*model.N)'; % transform to null space;
return

author: opencobra

repo: cobratoolbox


text:
%% pdf_Uniform: Uniform Probability Density Function
%
% Usage:
%   x = pdf_Uniform(N, range, seed)
%
% Inputs:
%     N                scalar, number of samples
%     range            vector [min max] range of the random variable
%     seed             optional, seed for the random number generator
%
% Output:
%     x                vector with the sampled data
%                      if N==0 -> x = range
%
% ------------------------------------------------------------------------
% See also
%
% Author : Flavio Cannavo'
% e-mail: flavio(dot)cannavo(at)gmail(dot)com
% Release: 1.0
% Date   : 28-01-2011
%
% History:
% 1.0  28-01-2011  First release.
%%

function x = pdf_Uniform(N, range, seed)

if nargin>2
    rand('twister',seed);
else
    rand('twister', round(cputime*1000 + 100000*rand));
end

if N==0
    x = range(:)';
else
    x = rand(N,1)*diff(range) + range(1);
end


text:
%compute the xcoordinate of the head in the tails coordinate system

function [data,units]=compute_xtailhead_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
xtailehad_mm=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    xtailehad_mm{1,i}=trx(larva).xhead_mm-trx(larva).xtail_mm;
end

units=parseunits('mm');
data=xtailehad_mm;

author: kristinbranson

repo: JAABA


text:
function plotFeatureImportanceRFvolClin_HN(pathRF)

startpath = pwd;

cd(pathRF), load('training'), load('testingVariableImportance')
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
% fSetNames = fieldnames(training.textures.(nameOutcomes{1})); nFset = numel(fSetNames);
fSetNames = {'volClinical'}; nFset = numel(fSetNames);

for o = 1:nOutcomes
    for f = 1:nFset
        percentAUCdecrease = variableImportance.(nameOutcomes{o}).(fSetNames{f}).percentAUCdecrease;
        varNames = variableImportance.(nameOutcomes{o}).(fSetNames{f}).varNames; nVar = numel(varNames);
        figure
        barh(1:nVar,percentAUCdecrease*100)
        for i = 1:nVar
            ind = strfind(varNames{i},'_'); nInd = numel(ind);
            if ~isempty(ind)
                for n = 1:nInd
                    varNames{i} = [varNames{i}(1:(ind(n)-1+(n-1))),'\',varNames{i}((ind(n)+(n-1)):end)];
                end
            end
        end
        set(gca,'yticklabel',varNames)
        nameOutcome = nameOutcomes{o};
        ind = strfind(nameOutcome,'Death');
        if ~isempty(ind)
            nameOutcome(ind:ind+4) = [];
            nameOutcome = ['Survival',nameOutcome];
        end
        title(['RANDOM PERMUTATIONS (RF):',nameOutcome,' -- ',fSetNames{f}])
        xlabel('Percent AUC decrease')
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function test_failed = test_libltfat_maxtree(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr =    [100,101,500,501,1025];
darr =    0;

for repeat = 1:1
  for depth = darr

for Lidx = 1:numel(Larr)

    L = Larr(Lidx);

    f = cast(randn(L,1)',flags.complexity);

    fPtr = libpointer(dataPtr,f);

    funname = makelibraryname('maxtree_initwitharray',flags.complexity,0);

    p = libpointer();

    calllib('libltfat',funname,L,depth,fPtr,p);

    maxPtr = libpointer(dataPtr,5);
    maxposPtr = libpointer('int64Ptr',cast(5,'int64'));

    funname = makelibraryname('maxtree_findmax',flags.complexity,0);
    status=calllib('libltfat',funname,p,maxPtr,maxposPtr);

    fprintf('max=%.3f, maxPos=%d\n',maxPtr.value, maxposPtr.value);

    fPtr.value(1) = 1000;

    funname = makelibraryname('maxtree_setdirty',flags.complexity,0);
    status=calllib('libltfat',funname,p,0,1);

    funname = makelibraryname('maxtree_findmax',flags.complexity,0);
    status=calllib('libltfat',funname,p,maxPtr,maxposPtr);

    fprintf('max=%.3f, maxPos=%d\n',maxPtr.value, maxposPtr.value);

    [fmax,fIdx] = max(fPtr.value);

    fprintf('max=%.3f, maxPos=%d\n', fmax, fIdx -1);
    fprintf('max=%.3f, maxPos=%d\n',maxPtr.value, maxposPtr.value);

    [test_failed,fail]=ltfatdiditfail(maxPtr.value-fmax + maxposPtr.value - (fIdx -1) ,test_failed,0);
    fprintf(['MAXTREE L:%3i, %s %s %s\n'],L,flags.complexity,ltfatstatusstring(status),fail);

    funname = makelibraryname('maxtree_done',flags.complexity,0);
    calllib('libltfat',funname,p);

end
  end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function a = r8mat_diag_add_scalar ( n, a, s )

%*****************************************************************************80
%
%% R8MAT_DIAG_ADD_SCALAR adds a scalar to the diagonal of a matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows and columns of the matrix.
%
%    Input, real A(N,N), the matrix to be modified.
%
%    Input, real S, the value to be added to the diagonal
%    of the matrix.
%
%    Output, real A(N,N), the modified matrix.
%
  for i = 1 : n
    a(i,i) = a(i,i) + s;
  end

  return
end


text:
function varargout=cent2frq(varargin)
%
% For calling details please see v_cent2frq.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_cent2frq.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: cent2frq.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_cent2frq(varargin{:});
else
    v_cent2frq(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function x = Adj_USFFT(n,y,omega,D,L,center)
% USFFT_simple -- 1d adjoint unequispaced Fourier transform
% Usage:
%   x = USFFT(n,y,omega,D,L,center)
% Inputs:
%   n      length of output x
%   y	   vector of length m,
%   omega  vector of sampled frequencies (length m)
%   center 0 for unbiased FT, 1 otherwise
% Outputs:
%   x      vector of length n
% Description:
%  Evaluates the adjoint of the FT
%
%  If center = 0,
%
%  x(t) = sum_{k} exp(i omega_k t) y(k), -n/2 <= t < n/2
%
%  If center = 1,
%
%  x(t) = sum_{k} exp(i omega_k t) y(k), 0 <= t < n
%
%  See Also
%   USFFT, USFT_simple, Evaluate_FT
%  Notes:
%   To work properly, D*n must be even
%
% By Emmanuel candes, 2003-2004

  if nargin < 6,
    center = 0;
  end

  if nargin < 5,
    L = 4;
  end

  if nargin < 4,
    D = 16;
  end

  n2 = n/2;
  N = D*n; N2 = N/2;
  m = length(omega);

  if rem(L,2) == 1,
    L = L + 1;
  end

  % Nearest point calculations

  near_index = round(N *omega/(2*pi));
  row = near_index + N/2 + 1;

  delta = omega - 2*pi*near_index/N;

  % Make matrix y(k) * delta(k)^l, l = 0, 1, ..., L-1

  FA = repmat(y(:),1,L);
  for l = 2:L;
    FA(:,l) = FA(:,l-1) .* delta/(l-1);
  end

  % Sum entries corresponding to the same index

  F = zeros(N,L);
  for k = 1:m,
    F(row(k),:) = F(row(k),:) + FA(k,:);
  end

  % Take IFFT along columns
  X = ifft_mid0(F) * N;
  X = X((N2-n2+1):(N2+n2),:);

  t = -n2:(n2-1);
  tp = ones(n,L);

  for k = 2:L;
    tp(:,k) = tp(:,k-1) .* (i*t.');
  end

  x = sum(X.*tp,2);

author: thomaskuestner


text:
function [values,modes] = max(SO3F,varargin)
% global, local and pointwise maxima of functions on SO(3)
%
% Syntax
%   [v,pos] = max(SO3F) % the position where the maximum is atained
%
%   [v,pos] = max(SO3F,'numLocal',5) % the 5 largest local maxima
%
%   SO3F = max(SO3F, c) % maximum of a rotational functions and a constant
%   SO3F = max(SO3F1, SO3F2) % maximum of two rotational functions
%   SO3F = max(SO3F1, SO3F2, 'bandwidth', bw) % specify the new bandwidth
%
%   % compute the maximum of a multivariate function along dim
%   SO3F = max(SO3Fmulti,[],dim)
%
% Input
%  SO3F, SO3F1, SO3F2 - @SO3Fun
%  SO3Fmulti          - a multivariate @SO3Fun
%  c                  - double
%
% Output
%  v - double
%  pos - @rotation / @orientation
%
% Options
%  kmax          - number of iterations
%  numLocal      - number of peaks to return
%  startingNodes - @rotation / @orientation
%  tolerance     - minimum distance between two peaks
%  resolution    - minimum step size
%  maxStepSize   - maximm step size
%
% Example
%
%   %find the local maxima of the <SantaFe.html SantaFe> ODF
%   mode = calcModes(SantaFe)
%   plotPDF(SantaFe,Miller(0,0,1,mode.CS))
%   annotate(mode)
%
% See also
% SO3Fun/min SO3Fun/max

if isa(SO3F,'SO3FunHarmonic') && ~SO3F.isReal
  SO3F = SO3F.isReal;
  warning('By taking the maxima of SO3Funs, the functions should be real valued.')
end
if nargin>1 && isa(varargin{1},'SO3FunHarmonic') && ~varargin{1}.isReal
  varargin{1}.isReal = 1;
  warning('By taking the maxima of SO3Funs, the functions should be real valued.')
end

if numel(SO3F)==1
  [values,modes] = max@SO3Fun(SO3F,varargin{:});
  return
end

% multivariate functions
s = size(SO3F);

if nargin>1 && (isa(varargin{1},'SO3FunHarmonic') || isnumeric(varargin{1}))
  t = size(varargin{1});
  SO3F1 = SO3F.*ones(t);
  SO3F2 = varargin{1}.*ones(s);
  values = [];
  for k=1:numel(SO3F1)
    if isa(SO3F2,'SO3FunHarmonic')
      A = max@SO3Fun(SO3F1.subSet(k),SO3F2.subSet(k),varargin{:});
    else
      A = max@SO3Fun(SO3F1.subSet(k),SO3F2(k),varargin{:});
    end
    values = [values,A];
  end
  values = reshape(values,size(SO3F1));
  return
end

len = get_option(varargin,'numLocal',1);
values = zeros(len,prod(s));
modes = rotation.id(len,prod(s));
for k=1:numel(SO3F)
  [v,m] = max@SO3Fun(SO3F.subSet(k),varargin{:});
  values(:,k)=v; modes(:,k)=m;
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [balanced_ds,idxs,classes]=cosmo_balance_dataset(ds,varargin)
% sub-sample a dataset to have an equal number of samples for each target
%
% [balanced_ds,idxs,classes]=cosmo_balance_dataset(ds)
%
% Inputs:
%   ds                      dataset struct with fields .samples,
%                           .sa.targets and .sa.chunks. All values in
%                           .sa.chunks must be different from each other.
%   'sample_balancer', f    (optional)
%                           function handle with signature
%                               [idxs,classes]=f(targets,seed)
%                           where idxs is a SxC vector with indices for C
%                           classes and S targets per class. If omitted a
%                           builtin function is used.
%   'seed', s               (optional, default=1)
%                           Seed to use for pseudo-random number generation
%
% Output:
%   balanced_ds             dataset with a subset of the samples from ds
%                           so that each target occurs equally often.
%                           Selection is (by default) done in a
%                           pseudo-determistic manner.
%   idxs                    SxC vector indicating which
%   classes                 Cx1 vector containing unique class labels
%
% Notes:
%   - this function is to be used with MEEG datasets. it is not intended
%     for fMRI data.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    defaults=struct();
    defaults.seed=1;
    defaults.sample_balancer=[];

    opt=cosmo_structjoin(defaults,varargin{:});

    check_inputs(ds,opt)

    [idxs,classes]=select_subset(ds.sa.targets,opt);
    balanced_ds=cosmo_slice(ds,idxs(:),1);

function [idxs,classes]=select_subset(targets,opt)
    sample_balancer=opt.sample_balancer;
    if isempty(sample_balancer)
        sample_balancer=@default_sample_balancer;
    end

    [idxs,classes]=sample_balancer(targets,opt.seed);

function [idxs,classes]=default_sample_balancer(targets,seed)
    [all_idxs,classes]=cosmo_index_unique(targets);
    class_counts=cellfun(@numel,all_idxs);

    min_count=min(class_counts);
    max_count=max(class_counts);

    nclasses=numel(class_counts);
    % invoke PRNG only once
    rand_vals=cosmo_rand(max_count,nclasses,'seed',seed);

    idxs=zeros(min_count,nclasses);
    for k=1:nclasses
        % set class_idxs to have values in the range 1:class_counts in
        % random order
        [unused,class_idxs]=sort(rand_vals(1:class_counts(k),k));

        % select min_count values from the indices in class_idxs
        idxs(:,k)=all_idxs{k}(class_idxs(1:min_count));
    end

function check_inputs(ds,opt)
    cosmo_check_dataset(ds);

    % chunks and targets must be present
    raise_exception=true;
    cosmo_isfield(ds,{'sa.targets','sa.chunks'},raise_exception);

    chunks=ds.sa.chunks;
    if numel(sort(chunks)) ~= numel(unique(chunks))
        error(['All values in .sa.chunks must be unique. If '...
                '*and only if* all '...
                'observations in .samples can be assumed to be '...
                'independent, for a dataset ds you can set '...
                '  ds.sa.chunks(:)=1:numel(ds.sa.chunks,1)'...
                'to indicate independence. This assumption typically '...
                'only applies to M/EEG datasets; this function should '...
                'not be used for typical fMRI datasets']);
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
clc;
clear all;
close all;
A=imread('cameraman.tif');
figure,imshow(uint8(A))
title('Original Image');
A=double(A);
[s1 s2]=size(A);
% bs=input('Enter the block sizes for division of the image: '); % Block Size
bs=256;

% Slant
temp=double(zeros(size(A)));
for y=1:bs:s1-bs+1
    for x=1:bs:s2-bs+1
        croppedImage = A((y:y+bs-1),(x:x+bs-1));
        t=getSlantTransform(croppedImage,bs);
        temp((y:y+bs-1),(x:x+bs-1))=t;
    end
end
figure,imshow(uint8(temp))

% Inverse Slant
temp1=double(zeros(size(A)));
for y=1:bs:s1-bs+1
    for x=1:bs:s2-bs+1
        croppedImage = temp((y:y+bs-1),(x:x+bs-1));
        t=getInvSlantTransform(croppedImage,bs);
        temp1((y:y+bs-1),(x:x+bs-1))=t;
    end
end
figure,imshow(uint8(temp1))


text:
HeadModelPath = [hlp_getSiftRoot filesep 'resources' filesep 'headmodels' filesep 'standard-Colin27-385ch.mat'];
hmObj = hlp_validateHeadModelObject(HeadModelPath);

%% Simulate VAR model
[EEGsim] = pop_sim_varmodel([]);
if length(EEGsim)>1
    EEGtrue = EEGsim(2);
    EEGsim  = EEGsim(1);
end

%% Simulate scalp activity using forward model

% these are the source locations
sourceRois = {'Cingulum_Mid_L','Occipital_Mid_L','Parietal_Sup_L','Frontal_Sup_Medial_R','Precentral_R'};

cfg = arg_guipanel('Function',@sim_simulateSources, ...
                   'Parameters',{ ...
                        'hmObj',HeadModelPath, ...
                        'sourceAtlasLabels', setdiff_bc(hmObj.atlas.label,{'Thalamus_L','Thalamus_R'}), ...
                        'Channels',hmObj.label, ...
                        'sourceShape' {'gausspatch' 'roiOrdered' sourceRois, 'sigma', 10}, ...
                        'addNoise', {'SignalToNoise' 20} ...
                        },'PanelOnly',false);
[scalpData srcData LFM] = sim_simulateSources('sourceAmps',EEGsim.data,cfg);

%%

chanlabels = hmObj.getChannelLabels();
elocs = hmObj.channelSpace;
EEG = EEGsim;
EEG.data       = scalpData;
EEG.nbchan     = size(EEG.data,1);
EEG.setname    = 'VAR Simulation';
EEG.condition  = 'VAR Simulation';
EEG.srcpot_all = srcData;
EEG = rmfield(EEG,'chanlocs');
for k=1:EEG.nbchan
    EEG.chanlocs(k) = struct(...
        'X',elocs(k,1), ...
        'Y',elocs(k,2), ...
        'Z',elocs(k,3), ...
        'labels',chanlabels{k}, ...
        'type','EEG', ...
        'theta',[], ...
        'radius',0.5, ...
        'sph_theta',[],...
        'sph_phi',[],...
        'sph_radius',[],...
        'urchan',k, ...
        'ref','');
end
EEG.chanlocs = convertlocs(EEG.chanlocs,'cart2all');

%% visualize
latency = [];
gobj=vis_csd('hmObj',HeadModelPath, ...
         'signal',EEG,  ...
         'cortexMesh',[], ...
         'times',latency,        ...
         'avgtimes',false,       ...
         'frameskip',10,          ...
         'cortexlims',99,       ...
         'scalplims',[99], ...
         'showpower',false, ...
         'title','True');

repo: REST


text:
function Omega = logm(ori,varargin)
% the logarithmic map that translates a rotation into a spin tensor
%
% Syntax
%   Omega = logm(mori) % spin tensor of a misorientation
%
%   Omega = logm(ori,ori_ref) % spin tensor in crystal coordinates
%   Omega = logm(ori,ori_ref,'left') % spin tensor in specimen coordinates
%
% Input
%  mori - misorientation
%  ori - @orientation
%  ori_ref - @orientation
%
% Output
%  Omega - @spinTensor
%
% See also
% orientation/log spinTensor/exp

m = log(ori,varargin{:});
Omega = spinTensor(m);

%T = logm@quaternion(ori,varargin{:});

%if check_option(varargin,'left')
%  T.CS = ori.SS;
%else
%  T.CS = ori.CS;
%end

end


text:
function A = construct_incidence_matrix(measurements)
%function A = construct_incidence_matrix(measurements)
%
% This function computes and returns the oriented incidence matrix of the
% underlying directed graph of measurements:
%
% A_{ie} =  -1, if edge e *leaves* node i,
%           +1, if edge e *enters* node i,
%           0, otherwise.
%
% (see eq. (7) in the paper).

% Copyright (C) 2016 by David M. Rosen

N = max(max(measurements.edges));  % Number of nodes in the pose graph
M = size(measurements.edges, 1);  % Number of edges in the pose graph

out_nodes = measurements.edges(:, 1)';  %out_nodes(e) = i if edge e leaves node i
in_nodes = measurements.edges(:, 2)';  %out_nodes(e) = j if edge e enters node i

node_indices = [out_nodes, in_nodes];
edge_indices = [ [1:M], [1:M] ];
vals = [-ones(1, M), ones(1, M)];

A = sparse(node_indices, edge_indices, vals, N, M);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
classdef SettingsElementalDensity < AbstractSettings
    properties
        filterParams = SettingsFilterP0()
        levelSetCreatorParams = SettingsLevelSetCreator()
    end
    methods
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef test_ChanDetails < matlab.unittest.TestCase
   %UNTITLED Summary of this class goes here
   %   Detailed explanation goes here

   properties
   end

   methods(Test)
      function TestRetrieve(testCase)
         ANMO = 'IU.ANMO.00.BHZ';
         ANTO = 'IU.ANTO.00.BHZ';
         % * test retrieve with default source *
         cdKey = ChannelDetails.retrieve([], 'station','ANMO',...
            'channel','BHZ','location','00','network','IU',...
            'starttime','2015-10-21');
         % make sure the first one works, so we can trust use it to
         % evaluate the other ones
         testCase.assertLength(cdKey, 1);
         testCase.assertEqual(cdKey.channelinfo,ChannelTag(ANMO));
         testCase.assertEqual(cdKey.samplerate, 20);
         testCase.assertEqual(cdKey.elevation, 1671);
         testCase.assertEqual(cdKey.latitude, 34.9459, 'absTol', 0.0001);
         testCase.assertEqual(cdKey.depth, 145);
         testCase.assertEqual(cdKey.azimuth, 0);
         testCase.assertEqual(cdKey.dip, -90);
         testCase.assertEqual(cdKey.sensordescription,'Geotech KS-54000 Borehole Seismometer');
         testCase.assertEqual(cdKey.scalefreq,0.0200);
         testCase.assertEqual(cdKey.scaleunits,'M/S');
         testCase.assertEqual(datestr(cdKey.starttime),'17-Dec-2014 18:40:00');
         testCase.assertGreaterThan(cdKey.endtime, cdKey.starttime);

         % test retrieve with a single N.S.L.C string
         chd = ChannelDetails.retrieve([], ANMO);
         testCase.verifyEqual(chd(end),cdKey);
         % test retrieve with ChannelTag
         chaninfo = ChannelTag(ANMO);
         chd = ChannelDetails.retrieve([], chaninfo);
         testCase.verifyEqual(chd(end),cdKey);
         % test retrieve with multiple channeltags
         chansinfo = ChannelTag({ANMO, ANTO});
         chd = ChannelDetails.retrieve([], [chansinfo; chansinfo]);
         testCase.verifySize(chd, [2 2], 'Wrong size returned for multiple Traces');
         % test retrieve with SeismicTrace
         T = SeismicTrace; T.name = ANMO;T.start = datenum('2015-10-21');
         chd = ChannelDetails.retrieve([],T);
         testCase.verifyEqual(chd,cdKey);
         % test retrieve with multiple Traces
         chd = ChannelDetails.retrieve([], [T T; T T]);
         testCase.verifySize(chd, [2, 2], 'Wrong size returned for multiple Traces');
      end
   end

end

author: geoscience-community-codes

repo: GISMO


text:
function score = IGD(Population,optimum)
% <min> <multi/many> <real/integer/label/binary/permutation> <large/none> <constrained/none> <expensive/none> <multimodal/none> <sparse/none> <dynamic/none> <robust/none>
% Inverted generational distance

%------------------------------- Reference --------------------------------
% C. A. Coello Coello and N. C. Cortes, Solving multiobjective optimization
% problems using an artificial immune system, Genetic Programming and
% Evolvable Machines, 2005, 6(2): 163-190.
%------------------------------- 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".
%--------------------------------------------------------------------------

    PopObj = Population.best.objs;
    if size(PopObj,2) ~= size(optimum,2)
        score = nan;
    else
        score = mean(min(pdist2(optimum,PopObj),[],2));
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function store = getStore(problem, x, storedb)
% Extracts a store struct. pertaining to a point from the storedb database.
%
% function store = getStore(problem, x, storedb)
%
% Queries the storedb database of structures (itself a structure) and
% returns the store structure corresponding to the point x. If there is no
% record for the point x, returns an empty structure.
%
% See also: setStore purgeStoredb

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    % Construct the fieldname (key) associated to the queried point x.
    key = problem.M.hash(x);

    % If there is a value stored for this key, return it.
    % Otherwise, return an empty structure.
    if isfield(storedb, key)
        store = storedb.(key);
    else
        store = struct();
    end

end

author: andrewssobral

repo: lrslibrary


text:
function sample()
% Create an output file with sys command
!touch test_file.txt
fid = fopen('test_file.txt', 'w');
for i=1:20
    fprintf(fid, '%d unterminated\n, i);
end
fclose(fid);
end

function unusedFunc(passedInput)
persistent global_value;
global_value = global_value + 1;
%% Title of cell
printed_var = subFunc(global_value)
unused_var = printed_var;
    function out = subFunc(scaleFactor)
        fprintf('Scale %.4f\n', scaleFactor);
        out = passedInput * scaleFactor;
    end
end

author: scottclowe

repo: matlab-schemer


text:
function [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = plot_McCdfnormal(mCatalog, fBinning)
% function [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = plot_McCdfnormal(mCatalog, fBinning);
% -----------------------------------------------------------------------------------------------------
% Same as calc_McCdfnormal with plotting the fitting steps and the final result
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Binning interval, usually 0.1
%
% Outgoing variables:
% mResult     : Solution matrix including
%               vProbability: maximum likelihood score
%               vMc         : Mc values
%               vX_res      : mu (of normal CDF), sigma (of normal CDF), residuum, exitflag
%               vNmaxBest   : Number of events in lowest magnitude bin considered complete
%               vABValue    : a and b-value
% fMls       : minimum maximum likelihood score --> best Mc
% fMc        : Best estimated magnitude of completeness
% mDatPredBest   : Matrix of non-cumulative FMD [Prediction, magnitudes, original distribution]
% vPredBest      : Matrix of non-cumulative FMD below Mc [magnitude, prediction, uncertainty of prediction]
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 10.02.03

% Initialize
vProbability = [];
vMc = [];
vABValue =[];
mFitRes = [];
vX_res = [];
vNCumTmp = [];
mDataPred = [];
vPredBest = [];
vDeltaBest = [];
vX_res = [];
vNmaxBest = [];
mResult=[];
mDatPredBest = [];

% Determine exact time period
fPeriod1 = max(mCatalog(:,3)) - min(mCatalog(:,3));

% Determine max. and min. magnitude
fMaxMag = ceil(10 * max(mCatalog(:,6))) / 10;

% Set starting value for Mc loop and LSQ fitting procedure
fMcTry= calc_Mc(mCatalog,1);
fSmu = fMcTry/2;
fSSigma = fMcTry/4;
if (fSmu > 1)
    fSmu = fMcTry/10;
    fSSigma = fMcTry/20;
end
fMcBound = fMcTry;

% Calculate FMD for original catalog
[vFMDorg, vNonCFMDorg] = calc_FMD(mCatalog);
fMinMag = min(vNonCFMDorg(1,:));

%% Shift to positive values
if fMinMag ~= 0
    fMcBound = fMcTry-fMinMag;
end

% Loop over Mc-values
for fMc = fMcBound-0.4:0.1:fMcBound+0.4
    fMc = round(fMc*10)/10;
    vFMD = vFMDorg;
    vNonCFMD = vNonCFMDorg;
    vNonCFMD = fliplr(vNonCFMD);
    % Calculate a and b-value for GR-law and distribution vNCum
    [nIndexLo, fMagHi, vSel, vMagnitudes] = fMagToFitBValue(mCatalog, vFMD, fMc);
    if (length(mCatalog(vSel,1)) >= 20)
        [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog(vSel,:), fBinning);
        % Normalize to time period
        vFMD(2,:) = vFMD(2,:)./fPeriod1;
        vNonCFMD(2,:) = vNonCFMD(2,:)./fPeriod1;
        % Compute quantity of earthquakes by power law
        fMaxMagFMD = max(vNonCFMD(1,:));
        fMinMagFMD = min(vNonCFMD(1,:));
        vMstep = [fMinMagFMD:0.1:fMaxMagFMD];
        vNCum = 10.^(fAValue-fBValue.*vMstep); % Cumulative number

        % Compute non-cumulative numbers vN from GR-law
        fNCumTmp = 10^(fAValue-fBValue*(fMaxMag+0.1));
        vNCumTmp  = [vNCum fNCumTmp ];
        vN = abs(diff(vNCumTmp));

        % Normalize vN
        vN = vN./fPeriod1;
        % Data selection
        % mData = Non-cumulative FMD values from GR-law and original data
        mData = [vN' vNonCFMD'];
        vSel = (mData(:,2) >= fMc);
        mDataTest = mData(~vSel,:);
        mDataTmp = mData.subset(vSel);
        % Check for zeros in observed data
%         vSelCheck = (mDataTest(:,3) == 0);
%         mDataTest = mDataTest(~vSelCheck,:);
        % Choices of normalization
        fNmax = mDataTmp(1,3); % Frequency of events in Mc bin
        %fNmax = max(mDataTest(:,3));  % Use maximum frequency of events in bins below Mc
        %fNmax = mDataTest(length(mDataTest(:,1)),3); % Use frequency of events at bin Mc-0.1 -> best fit
        %fNmax = (mDataTest(length(mDataTest(:,1)),3)+mDataTmp(1,3))/2;
        if (~isempty(fNmax) & ~isnan(fNmax) & fNmax ~= 0 & length(mDataTest(:,1)) > 4)
            mDataTest(:,3) = mDataTest(:,3)/fNmax; % Normalize datavalues for fitting with CDF
            % Move to M=0 to fit with lsq-algorithm
            fMinMagTmp = min(mDataTest(:,2));
            mDataTest(:,2) = mDataTest(:,2)-fMinMagTmp;
            % Curve fitting: Non cumulative part below Mc
            options = optimset;
            %options = optimset('Display','off','Tolfun',1e-5,'TolX',0.0001,'MaxFunEvals', 100000,'MaxIter',10000);
            options = optimset('Display','off','Tolfun',1e-5,'TolX',0.001,'MaxFunEvals', 1000,'MaxIter',1000);
            [vX, resnorm, resid, exitflag, output, lambda, jacobian]=lsqcurvefit(@calc_normalCDF,[fSmu  fSSigma], mDataTest(:,2), mDataTest(:,3),[],[],options);
            mDataTest(:,1) = normcdf(mDataTest(:,2), vX(1), vX(2))*fNmax;
            if (length(mDataTest(:,2)) > length(vX(1,:)))
                %% Confidence interval determination
                % vPred : Predicted values of lognormal function
                % vPred+-delta : 95% confidence level of true values
                [vPred,delta] = nlpredci(@calc_normalCDF,mDataTest(:,2),vX, resid, jacobian);
            else
                vPred = NaN;
                delta = NaN;
            end % END: This section is due for errors produced with datasets less long than amount of parameters in vX
            % Results of fitting procedure
            mFitRes = [mFitRes; vX resnorm exitflag];
            % Move back to original magnitudes
            mDataTest(:,2) = mDataTest(:,2)+fMinMagTmp;
            %% Set data together
            mDataTest(:,3) = mDataTest(:,3)*fNmax;
            mDataPred = [mDataTest; mDataTmp];
            % Denormalize to calculate probabilities
            mDataPred(:,1) = round(mDataPred(:,1).*fPeriod1);
            mDataPred(:,3) = mDataPred(:,3).*fPeriod1;
            vProb_ = calc_log10poisspdf2(mDataPred(:,3), mDataPred(:,1)); % Non-cumulative
            %vPro_ = calc_log10poisspdf(mDataPred(:,3), mDataPred(:,1)); % Non-cumulative
            % Sum the probabilities
            fProbability = (-1) * sum(vProb_);
            vProbability = [vProbability; fProbability];
            % Move magnitude back
            mDataPred(:,2) = mDataPred(:,2)+fMinMag;
            vMc = [vMc; fMc+fMinMag];
            vABValue = [vABValue; fAValue fBValue];

             % Keep values
            vDeltaBest = [vDeltaBest; delta];
            vX_res = [vX_res; vX resnorm exitflag];
            vNmaxBest = [vNmaxBest; fNmax];

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Plot routines
            if exist('cdfit_fig','var') &  ishandle(cdfit_fig)
                set(0,'Currentfigure',cdfit_fig);
            else
                cdfit_fig=figure_w_normalized_uicontrolunits('tag','cdfit','Name','Fitting CDF','Units','normalized','Nextplot','add',...
                    'Numbertitle','off','visible','on');
                cdfit_axs=axes('tag','ax_cdfit','Nextplot','add','box','off');
            end

            set(gcf,'tag','cdfit');
            subplot(3,1,1); % Fitting curve and original data => Mc
            plot(mDataTest(:,2)+fMinMag, mDataTest(:,1),'-r', mDataTest(:,2)+fMinMag, mDataTest(:,3), '*')
            hold on;
            plot(mDataTest(:,2)+fMinMag, mDataTest(:,1)+delta*fNmax,'--g',mDataTest(:,2)+fMinMag, mDataTest(:,1)-delta*fNmax,'--g');
            hold off;
            ylabel('CDF fit');
            sTitle = ['Test: Mc = ', num2str(fMc+fMinMag)];
            title(sTitle);
            subplot(3,1,2); % Non-cumulative data, Non-cumulative data from GR-law, Non-cumulative data predicted with normal CDF
            plot(vNonCFMD(1,:)'+fMinMag, vNonCFMD(2,:)'.*fPeriod1, '^b', vNonCFMD(1,:)'+fMinMag, vN.*fPeriod1, '*g', mDataPred(:,2),mDataPred(:,1),'or')
            ylabel('NonCumFMD fit')
            subplot(3,1,3); % Cumlative data, cumulative data from GR-law, cumulative data predicted with normal CDF
            mTmp = flipud(mDataPred);
            mTmp(:,1) = cumsum(mTmp(:,1));
            mTmp(:,3) = cumsum(mTmp(:,3));
            semilogy(mTmp(:,2), mTmp(:,1),'ro',mTmp(:,2), mTmp(:,3),'g*',vMstep+fMinMag, vNCum./fPeriod1)
            xlabel('Magnitude')
            ylabel('FMD fit')
            drawnow;
            %pause
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

             % Keep best fitting model for plotting
            if (fProbability == nan(vProbability))
                vDeltaBest = delta;
                vPredBest = [mDataTest(:,2)+fMinMag vPred*fNmax*fPeriod1 delta*fNmax*fPeriod1]; % Gives back uncertainty
                vNBest = vN;
                %fMc+fMinMag : Test procedure
                mDatPredBest = [mDataPred];
           end
        else
            %disp('Not enough data');
            % Setting values
            fProbability = NaN;
            fMc = NaN;
            vX(1) = NaN;
            vX(2) = NaN;
            resnorm = NaN;
            exitflag = NaN;
            delta = NaN;
            vPred = [NaN NaN NaN];
            fNmax = NaN;
            fAValue = NaN;
            fBValue = NaN;
            vProbability = [vProbability; fProbability];
            vMc = [vMc; fMc+fMinMag];
            vX_res = [vX_res; vX resnorm exitflag];
            vDeltaBest = [vDeltaBest; delta];
            vPredBest = [vPredBest; vPred*fNmax];
            vNmaxBest = [vNmaxBest; fNmax];
            vABValue = [vABValue; fAValue fBValue];
        end % END of IF fNmax
    end % END of IF length(mCatalog(vSel,1))

    % Clear variables
    vNCumTmp = [];
    mModelDat = [];
    vNCum = [];
    vSel = [];
    mDataTest = [];
    mDataPred = [];
end % END of FOR fMc

if isempty(vProbability)
    disp('Not enough data');
    return;
end

mResult = [mResult; vProbability vMc vX_res vNmaxBest vABValue];

% Find best estimate, excluding the case of mResult all NAN
if  ~isempty(nan(mResult))
    if ~isnan(nan(mResult(:,1)))
        vSel = find(nan(mResult(:,1)) == mResult(:,1));
        fMc = min(mResult(vSel,2));
        fMls = min(mResult(vSel,1));
        fMu = min(mResult(vSel,3));
        fSigma = min(mResult(vSel,4));
    else
        fMc = NaN;
        fMls = NaN;
        fMu = NaN;
        fSigma = NaN;

    end
else
    fMc = NaN;
    fMls = NaN;
    fMu = NaN;
    fSigma = NaN;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plot result
% Plot MLS vs. Mc
figure_w_normalized_uicontrolunits('tag','mlsmc','Name','MLS vs. Mc','Units','normalized','Nextplot','add',...
                    'Numbertitle','off','visible','on');
plot(vMc, vProbability,'Marker','s','MarkerFaceColor',[0.5 0.5 0.5],'Markersize',8');
xlabel('Mc');
ylabel('MLS');
sTitlestr = ['Minimum MLS at Mc = ' num2str(fMc)];
title(sTitlestr);

% Plot Best fitting CDF vs. magnitude
figure_w_normalized_uicontrolunits('tag','best_cdfit','Name','CDF fit','Units','normalized','Nextplot','add',...
                    'Numbertitle','off','visible','on');
plot(vPredBest(:,1),vPredBest(:,2))
hold on;
%plot(mDat(:,2)+fMinMag, mDat(:,3)*fNmax*fPeriod1,'+r')
plot(mDatPredBest(:,2),mDatPredBest(:,1),'o',mDatPredBest(:,2), mDatPredBest(:,3),'<');
sTitlestr = ['mu = ' num2str(fMu) ', sigma = ' num2str(fSigma)];
title(sTitlestr)
plot(vPredBest(:,1),vPredBest(:,2)+vPredBest(:,3),'--g')
plot(vPredBest(:,1),vPredBest(:,2)-vPredBest(:,3),'--g')
xlabel('Magnitude');
ylabel('Non-cumulative FMD');
hold off;

% Plot Non-cumulative distribution, original and predicted
figure_w_normalized_uicontrolunits('tag','ncumdist','Name','Best model','Units','normalized','Nextplot','add',...
                    'Numbertitle','off','visible','on');
semilogy(vNonCFMDorg(1,:)', vNonCFMDorg(2,:)', '^', vNonCFMD(1,:)', vNBest.*fPeriod1, '*', mDatPredBest(:,2),mDatPredBest(:,1),'o')
xlim = ([min(mDatPredBest(:,2)) max(mDatPredBest(:,2))]);
sTitlestr = ['Mc = ' num2str(fMc) ' using Normal CDF fitting'];
title(sTitlestr)
xlabel('Magnitude');
ylabel('Non-cumulative FMD');
drawnow;

author: CelsoReyes


text:
function Y = linsolve_cov(C,Y,type)
% X = linsolve_cov(C,Y)
%
% Solves C*X = Y, where C is symmetric positive definite matrix.

Y = linsolve_chol(chol(C),Y);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function tf = istype(obj,  typeName)
    % ISTYPE Determine if input object is of a particular type
    %   ISTYPE(obj, 'TypeName') returns true if obj contains a property named 'Type', and
    %   the contents match.  This is designed for use with graphical objects.
    %
    %   also, will work with structs that have fields named 'Type'
    %
    % SEE ALSO isa

    tf = (isobject(obj) & isprop(obj,'Type')) | ...% object with a Type property OR
        (isstruct(obj) & isfield(obj,'Type') ); % struct with a Type field
    tf(tf)= arrayfun(@(x)ischar(x.Type)||isstring(x.Type)&&strcmp(x.Type,typeName),obj(tf));
end