text:
function [V,F3,F4] = load_obj_file(modelname, opts)
% loadawobj load a Wavefront/Alias obj style model.
%    It will only consider polygons with 3 or 4 vertices.
%    [v,f3,f4]=loadawobj(modelname)
%    s=loadawobj(modelname)
%    s=loadawobj(modelname,1) % will load vertex normals and textures
% Programme will normally ignore normal and texture data unless the
% option is given.
% It will also ignore any part of the obj specification that is not a
% polygon mesh, ie nurbs and the materials file (.mtl). Please let me
% know if these are needed.
%
% Examples:
%   loadawobj('name.obj') will load and display name.obj
%   S=loadawobj('name.obj') will load obj into S
%   patch('Vertices',S.v','Faces',S.f3','FaceColor','g');
%
%
% s=loadawobj(modelname) gives a structure
%  version
%  vertices
%  f3 (triangles)
%  f4 (quadrilaterals)
%  groups
%  g3 faces per group
%  g4 faces per group
%  Vn3 Vertex normals indices for f3
%  Vn4 Vertex normals indices for f4
%  Vertex normals (if option is given)
%  Vertex textures (if option is given)
%
% where the file is grouped the following code fragment can be adapted
% to split out groups
%
% clf; axis('equal'); col='mybgrcmy';
% for jj=3:length(S.g)+2
%     patch('Vertices',S.v','Faces',S.f3(:,S.g3(jj-1):S.g3(jj)-1)','FaceColor',col(jj));
%     pause;
% end;
%
% See also Anders Sandberg's vertface2obj.m and saveobjmesh.m
%
% there is still a good chance that obj files will not load. I
% would be grateful for any reports and examples of those that fail.
%
% W.S. Harwin, University Reading, 2006,2010,2015.
% Matlab BSD license
% thanks also to Doug Hackett

version=0.3;
if nargin <1
  disp('specify model name')
end

fid = fopen(modelname,'r');
if (fid<0)
  error(['can not open file: ' modelname]);
  return ;
end

vnum=1;
f3num=1; % triangle (f3) faces so far
f4num=1; % quad (f4) faces so far
vtnum=1; % vertex textures so far
vnnum=1; % vertex normals so far
gnum=0;  % groups so far
g3num=1; %
g4num=1;
Vtmp=[]; % set to null so that it can be assigned in the structure

% Line by line passing of the obj file

while ~feof(fid)
  Ln=fgets(fid);
  Ln=removespace(Ln);
  objtype=sscanf(Ln,'%s',1);
  l=length(Ln);
  if l==0  % isempty(s) ;
%    disp(['empty' Ln]);
    continue
  end
%  disp(Ln)
  Lyn=Ln; %temp hack
  switch objtype
    case '#' % comment
%       disp(Ln);
    case 'v' % vertex
      v=sscanf(Ln(2:end),'%f');
      Vtmp(:,vnum)=v;
      vnum=vnum+1;
    case 'vt'	% textures
      if nargin>=2; % option to collect
          vtxtr=sscanf(Ln(3:end),'%f');
          Vtexture(:,vtnum)=vtxtr;
          vtnum=vtnum+1;
      end
    case 'g' % sub mesh
%       disp(Ln);
%      if f3num > g3num(end);
      g3num=[g3num f3num];
      g4num=[g4num f4num];
      gnum=gnum+1;
      g{gnum}=Ln(3:end);
    case 'mtllib' % material library
%         disp(Ln);
    case 'usemtl' % use this material name
%         disp(Ln);
    case 'l' % Line
%         disp(Ln);
    case 's' %smooth shading across polygons
%         disp(Ln);
    case 'vn' % normals
      if nargin>=2; % option to collect
          vnorms=sscanf(Ln(3:end),'%f');
          Vnorms(:,vnnum)=vnorms;
          vnnum=vnnum+1;
      end
    case 'f' % faces
      nvrts=length(findstr(Ln,' ')); % spaces as a predictor of n vertices
      slashpat=findstr(Ln,'/');
      nslash=length(slashpat);
      if nslash >1 % dblslash can be 0, 1 or >1
        dblslash=slashpat(2)-slashpat(1); else dblslash=0;
      end
      Ln=Ln(3:end); % get rid of the f
      if nslash == 0 % Face = vertex
%        disp('-v');
        f1=sscanf(Ln,'%f');
      elseif nslash == nvrts && dblslash>1 % Face = v/tc
%        disp('-v/tc');
        data1=sscanf(Ln,'%f/%f');
        if nvrts == 3
          f1=data1([1 3 5]);
          tc1=data1([2 4 6]);
        end
        if nvrts == 4;
          f1=data1([1 3 5 7]);
          tc1=data1([2 4 6 8]);
        end
      elseif nslash == 2*nvrts && dblslash==1 % v//n
%        disp('-v//n');
        data1=sscanf(Ln,'%f//%f');
        if nvrts == 3
          f1=data1([1 3 5]);
          vn1=data1([2 4 6]);
          Vn3(:,f3num)=f1;
        end
        if nvrts == 4;
          f1=data1([1 3 5 7]);
          vn11=data1([2 4 6 8]);
          Vn4(:,f4num)=f1;
        end
      elseif nslash == 2*nvrts && dblslash>1 % v/tc/n
%        disp('-v/tc/n');
        data1=sscanf(Ln,'%f/%f/%f');
        if nvrts == 3
          f1=data1([1 4 7]);
          tc1=data1([2 5 8]);
          vn1=data1([3 6 9]);
          Vn3(:,f3num)=f1;
        end
        if nvrts == 4;
          f1=data1([1 4 7 10]);
          tc1=data1([2 5 8 11]);
          vn1=data1([3 6 9 12]);
          Vn4(:,f4num)=f1;
        end
      end
% Now put the data into the array(s)
      if nvrts == 3
          F3(:,f3num)=f1;
          f3num=f3num+1;
      elseif nvrts ==4
          F4(:,f4num)=f1;
          f4num=f4num+1;
      else
          warning(sprintf('v nvrts=%d %s',nvrts, Ln));
      end
    otherwise
%      if ~strcmp(Lyn,char([13 10])) % carriage return
%        if ~ all(Lyn == ' ') % ignore lines with only space
%           disp(['unprocessed-' Ln '-']); % see what has not been processed
%          double(Ln)
          %        end
%      end
    end

end

fclose(fid);

% plot if no output arguments are given
if nargout ==0
  if exist('F3','var')
    patch('Vertices',Vtmp','Faces',F3','FaceColor','g');
  end
  if exist('F4','var')
    patch('Vertices',Vtmp','Faces',F4','FaceColor','b');
  end
  axis('equal')
  clear Vtmp F3 F4
end

if nargout >=2
  V=Vtmp;
  if ~ exist('F3','var')
    warning('No 3 element faces')
    F3=[];
  end
  if nargout ==3
    if ~ exist('F4','var')
      warning('No 4 element faces')
      F4=[];
    end
  end
end

if nargout ==1
  V.version=version;
  V.v=Vtmp;
  if exist('F3','var')
      V.f3=F3;  end
  if exist('F4','var')
    V.f4=F4;  end
  if gnum>0
      V.g=g;  end
  V.g3=[g3num f3num];
  V.g4=[g4num f4num];
  if exist('Vn3','var')
    V.Vn3=Vn3;  end
  if exist('Vn4','var')
    V.Vn4=Vn4;  end
  if exist('Vnorms','var')
    V.Vnorms=Vnorms;  end
  if exist('Vtexture','var')
    V.Vtexture=Vtexture;  end
end

function Lyn=removespace(Lyn)
% A not an elegant way to remove
% surplus space
Lyn=strtrim(Lyn);
Lyn=strrep(Lyn,'       ',' '); % 8-2 .. 12-6
Lyn=strrep(Lyn,'    ',' '); % 5-2 6-3 4-1
Lyn=strrep(Lyn,'  ',' '); % 3-2 2-1
Lyn=strrep(Lyn,'  ',' ');
Lyn=strrep(Lyn,char([13 10]),''); % remove cr/lf
Lyn=strrep(Lyn,char([10]),''); % remove lf

author: yuxng


text:
function imgOut = SchlickTMO(img, s_mode, p, nBit, L0, k)
%
%       imgOut = SchlickTMO(img, s_mode, p, nBit, L0, k)
%
%
%       Input:
%           -img: input HDR image.
%           -s_mode: parameters selection: 'manual', 'automatic', and 'nonuniform'.
%           -p: a model parameter which takes values in [1,+inf].
%           -nBit: number of bit for the quantization step.
%           -L0: lowest value of the LDR monitor that can be perceived
%           by the HVS.
%           -k: a value in [0,1].
%
%       Output
%           -imgOut: tone mapped image
%
%     Copyright (C) 2010-15  Francesco Banterle
%
%     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/>.
%
%     The paper describing this technique is:
%     "Quantization Techniques for Visualization of High Dynamic Range Pictures"
% 	  by Christophe Schlick
%     in "Photorealistic Rendering Techniques" 1995
%

checkNegative(img);

check13Color(img);

check3Color(img);

if(~exist('s_mode', 'var'))
    s_mode = 'automatic';
end

if(~exist('nBit', 'var'))
    nBit = 8;
end

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

if(~exist('k','var'))
    k = 0.5;
end

if(~exist('p', 'var'))
    p = 1 / 0.005;
end

%compute the luminance channel
L = lum(img);

%compute the min luminance
LMin = min(L(L > 0.0));

%compute the max luminance
LMax = max(L(L>0.0));

%mode selection
switch s_mode
    case 'manual'
        p = max([p, 1]);

    case 'automatic'
        p = L0 * LMax / (2^nBit * LMin);

    case 'nonuniform'
        p = L0 * LMax / (2^nBit * LMin);
        p = p * (1 - k + k * L / sqrt(LMax * LMin));
end

%dynamic range reduction
Ld = p .* L ./ ((p - 1) .* L + LMax);

%change luminance
imgOut = ChangeLuminance(img, L, Ld);

end

author: banterle


text:
function c=ref_dgt_4(f,g,a,M)
%REF_DGT_4  DGT algorithm 4
%
%  This algorithm makes the r-loop be the outermost loop, in order to
%  reduce the size of the intermediate buffers, and hopefully better
%

  L=size(g,1);
N=L/a;
b=L/M;

[c,h_a,h_m]=gcd(-a,M);
p=a/c;
q=M/c;
d=N/q;

w=zeros(M,N);

% This version uses matrix-vector products and ffts

F=zeros(d,p,q);

C=zeros(d,q,q);

% Setup G before we start
G=zeros(c,d,p,q);
for r=0:c-1
  for s=0:d-1
    for k=0:p-1
      for l=0:q-1
        G(r+1,s+1,k+1,l+1)=sqrt(M*d)*g(mod(r+k*M-l*a+s*p*M,L)+1);
      end;
    end;
  end;
end;
G=dft(G,[],2);

% Set up the matrices
for r=0:c-1

  for s=0:d-1
    for k=0:p-1
      for l=0:q-1
        F(s+1,k+1,l+1)=f(mod(r+k*M+s*p*M-l*h_a*a,L)+1);
      end;
    end;
  end;

  % fft them
  F=dft(F,[],1);

  % Multiply them
  for s=0:d-1
    GM=reshape(G(r+1,s+1,:,:),p,q);
    FM=reshape(F(s+1,:,:),p,q);
    C(s+1,:,:)=GM'*FM;
  end;

  % Inverse fft
  C=idft(C,[],1);

  % Place the result
  for l=0:q-1
    for u=0:q-1
      for s=0:d-1
        w(r+l*c+1,mod(u+s*q-l*h_a,N)+1)=C(s+1,u+1,l+1);
      end;
    end;
  end;
end;

c=dft(w);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% barvar(X,Y,W)
%
% using patches to re-implement the bar-plot function, allowing different
% width for each bar.
%
% columns in Y are regarded as a data series (same color). rows are grouped
% into adjacent bars having different colors.
%
% to set the width for each data series input a row vector:
% W = [w1,w2,w3...]. to a set a different width for each bars group input a
% column vector: W = [w1;w2;w3...].

function hBar = barvar( X, Y, W)

if nargin < 3
    fprintf( 'this is a demo of barvar(X,Y,Z). \nnext time give me all the required inputs.\n');
    W = randi(3,1,5);
    X = (1:5)' * sum(W) * 1.5;
    Y = 10*(rand(5,5) - 0.5);
end

% rectify input, if possible
for dim = 1:2
    Repeat = ones(1,2);
    if size(X,dim) ~= size(Y,dim)
        if size(X,dim) == 1
            Repeat(dim) = size( Y,dim);
            X = repmat( X, Repeat);
        elseif size(Y,dim) == 1
            Repeat(dim) = size( X,dim);
            Y = repmat( Y, Repeat);
        else
            error( 'each input dimension should be either 1 or equal to other inputs.');
        end
    end
    if size(W,dim) ~= size(Y,dim)
        if size(W,dim) == 1
            Repeat(dim) = size( Y,dim);
            W = repmat( W, Repeat);
        else
            error( 'each input dimension should be either 1 or equal to other inputs.');
        end
    end
end

%% draw bars

nSets = size( Y,2);
nBars = size( Y,1);
hBar = zeros( 1, nSets);
SetWidth = sum( W,2); % width of combined sets per bar
cdata = (1:nSets)';

for s = 1:nSets

    % 4 vertices for each bar
    Xdata = zeros(4,nBars);
    Ydata = zeros(4,nBars);
    for b = 1:nBars
        % X-left
        Xdata(1,b) = X(b,s); % ref data point
        Xdata(1,b) = Xdata(1,b) - SetWidth(b) / 2; % begining of set (centered on Xdata)
        Xdata(1,b) = Xdata(1,b) + sum( W(b,1:s-1)); % current set position
        Xdata(2,b) = Xdata(1,b);
        % X-right = Xl + W
        Xdata(3,b) = Xdata(1,b) + W(b,s);
        Xdata(4,b) = Xdata(3,b);
        % Y-down
        Ydata(1,b) = 0;
        Ydata(4,b) = 0;
        % Y-up
        Ydata(2,b) = Y(b,s);
        Ydata(3,b) = Y(b,s);
    end

    hBar(s) = patch( Xdata, Ydata, s);
%     set(gca,'CLim',[0 40])
%     set( hBar(s), 'FaceColor', 'flat', 'FaceVertexCData', cdata(b), 'CDataMapping', 'direct');
    hold on
end

% x-axis
line( [ min(X(:)) - 100, max(X(:)) + 100], [0,0], 'Color', 'k');
hold off
AxisMargin = max( sum(W,2));
set( gca, 'Box', 'on', 'Xlim', [ min(X(:)) - AxisMargin, max(X(:)) + AxisMargin]);
figure(gcf)

author: Sable

repo: mcbench-benchmarks


text:
% eeg_anova_script - example task ERP analysis

%  Analysis script for ERP components in
%  fixed and variable task waveforms for PTSD-PET project

% Define sample rate & epoch parameters
fprintf('Defining ERP parameters\n');
sample_rate = 2.5;
epoch_start = -200;
epoch_end = 1500;
points = 680;
times = linspace(epoch_start,(epoch_end - sample_rate),points);

% Define component timing & rows
comp = [145 195 250 385 550 750 950 1150 1350];
comp_rows = (comp - epoch_start) / sample_rate;

% Load data (680 rows x 124 cols)
path = 'd:\data_emse\ptsdpet\link14hz\';
for i=1:10,
    file = sprintf('%sc%02do_link14hz.txt',path,i);
    var  = sprintf('c%02do_link14hz',i);
    if ~exist(var,'var')
        fprintf('loading %s\n',file);
        load(file);
    end
    file = sprintf('%sc%02dt_link14hz.txt',path,i);
    var  = sprintf('c%02dt_link14hz',i);
    if ~exist(var,'var')
        fprintf('loading %s\n',file);
        load(file);
    end
    file = sprintf('%sp%02do_link14hz.txt',path,i);
    var  = sprintf('p%02do_link14hz',i);
    if ~exist(var,'var')
        fprintf('loading %s\n',file);
        load(file);
    end
    file = sprintf('%sp%02dt_link14hz.txt',path,i);
    var  = sprintf('p%02dt_link14hz',i);
    if ~exist(var,'var')
        fprintf('loading %s\n',file);
        load(file);
    end
end

% Select data rows for each component

cont = zeros(1,2,124);
ptsd = zeros(1,2,124);

for i = 1:max(size(comp))
    c = comp(i);
    t = comp_rows(i);
    for j=1:10,
        sub = sprintf('c%02do_link14hz',j);
        fprintf('Component %d, selecting subject %s\n',c,sub);
        data = eval(sub);
        cont(j,1,:) = data(t,:);

        sub = sprintf('c%02dt_link14hz',j);
        fprintf('Component %d, selecting subject %s\n',c,sub);
        data = eval(sub);
        cont(j,2,:) = data(t,:);

        sub = sprintf('p%02do_link14hz',j);
        fprintf('Component %d, selecting subject %s\n',c,sub);
        data = eval(sub);
        ptsd(j,1,:) = data(t,:);

        sub = sprintf('p%02dt_link14hz',j);
        fprintf('Component %d, selecting subject %s\n',c,sub);
        data = eval(sub);
        ptsd(j,2,:) = data(t,:);
    end

    % Combine component data into single matrix
    data = [cont;ptsd]; % 20x2x124

    % Output ANOVA analysis.
    file = sprintf('%scomp%d_anova.txt',path,c);
    ANOVA = fopen(file,'w');
    file = sprintf('%scomp%d_p_task.txt',path,c);
    TASK = fopen(file,'w');
    file = sprintf('%scomp%d_p_group.txt',path,c);
    GROUP = fopen(file,'w');
    file = sprintf('%scomp%d_p_taskxgroup.txt',path,c);
    TASKxGROUP = fopen(file,'w');

    formatstr = '';

    fprintf('Doing ANOVA on electrode:\n');

    x = 0;
    for k=1:124,

        % Extract data for analysis
        anova_data = data(:,:,k);

        if (x==24),
            x=0;
            fprintf('%3d\n',k);
        else,
            x = x + 1;
            fprintf('%3d ',k);
        end
        fprintf(ANOVA,'\n\nELECTRODE %3d\n',k);

        % Calculate & output summary stats
        avg.cont = mean(anova_data( 1:10,:));
        avg.ptsd = mean(anova_data(10:20,:));
        sd.cont = std(anova_data( 1:10,:));
        sd.ptsd = std(anova_data(10:20,:));
        se.cont = sd.cont / sqrt(10);
        se.ptsd = sd.ptsd / sqrt(10);
        fprintf(ANOVA,'%-12s\t%-12s\t%12s\t%12s\t%12s\n','Group','Task','Mean','StDev','SE');
        fprintf(ANOVA,'%-12s\t%-12s\t%12s\t%12s\t%12s\n','-----','----','----','-----','--');
        fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n','Cont','Fixed',   avg.cont(1),sd.cont(1),se.cont(1));
        fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n','Cont','Variable',avg.cont(2),sd.cont(2),se.cont(2));
        fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n','PTSD','Fixed',   avg.ptsd(1),sd.ptsd(1),se.ptsd(1));
        fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n\n','PTSD','Variable',avg.ptsd(2),sd.ptsd(2),se.ptsd(2));

        % Run 2way ANOVA
        [p, tab] = anova2(anova_data,10,'off');

        % Collate the p values
        task(k) = p(1);
        group(k) = p(2);
        taskxgroup(k) = p(3);

        % Output ANOVA table
        tab(2,1) = {'Task'};
        tab(3,1) = {'Group'};
        tab(4,1) = {'Task x Group'};
        for m=1:6,
            for n=1:6,
                if n==1, fprintf(ANOVA,'%-12s\t',tab{m,n});
                else
                    if m==1, fprintf(ANOVA,'%12s\t',tab{m,n});
                    else     fprintf(ANOVA,'%12.6f\t',tab{m,n});
                    end
                end
            end
            fprintf(ANOVA,'\n');
        end

        % This string is used to output all the p values below
        formatstr = strcat(formatstr,'%12.6f\t');

    end
    fclose(ANOVA);

    % Output the p values
    fprintf(TASK,formatstr,task); fclose(TASK);
    fprintf(GROUP,formatstr,group); fclose(GROUP);
    fprintf(TASKxGROUP,formatstr,taskxgroup); fclose(TASKxGROUP);

end

author: PatternRecognition

repo: OpenBMI


text:
%% calculation of RMS error
%After putting breakpoint in DPMpcOneBeam_aavc line 284, run the following code:
[xGrid,yGrid,weightMap] = plotIntensityMap(whichBeam,planC);
wtAtPB = [];
for i=1:length(xPosV)
    wtAtPB(i) = interp2(xGrid,yGrid,weightMap,xPosV(i)/10,yPosV(i)/10,'nearest',0);
end
MUweight = sum(planC{indexS.IM}(end).IMDosimetry.LeafSeq.MU{whichBeam});
disp('RMS ERROR:')
disp(sqrt(sum((wtAtPB-w_field*MUweight).^2)/length(w_field))/max(w_field*MUweight))

disp(sqrt(sum((wtAtPB-w_field*MUweight).^2)/length(w_field))/max(w_field*MUweight))

indKeep = find(wtAtPB>max(wtAtPB)*0.1);
x = wtAtPB(indKeep);
y = w_field(indKeep)*MUweight;
sqrt(1/length(x)*sum(((x-y)./x).^2))

%% Calculate RMS error for dose
baseIndex = 5;
newIndex = 6;
indKeep = find(planC{indexS.dose}(baseIndex).doseArray>max(planC{indexS.dose}(baseIndex).doseArray(:))*0.5);
x = planC{indexS.dose}(baseIndex).doseArray(indKeep);
y = planC{indexS.dose}(newIndex).doseArray(indKeep);
sqrt(1/length(x)*sum(((x-y)./x).^2))

%% Scale beamletWeights
% for i=1:length(planC{indexS.IM}(end).IMDosimetry.beams)
%     for j=1:size(planC{indexS.IM}(end).IMDosimetry.beams(i).beamlets,1)
%         for k=1:size(planC{indexS.IM}(end).IMDosimetry.beams(i).beamlets,2)
%             planC{indexS.IM}(end).IMDosimetry.beams(i).beamlets(j,k).maxInfluenceVal = planC{indexS.IM}(end).IMDosimetry.beams(i).beamlets(j,k).maxInfluenceVal * 11.2/23.668;
%         end
%     end
% end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef GemLabelEditBox < GemVirtualPanel
    % GemLabelEditBox GEM class for an edit box with label text
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = protected)
        Text
        EditBox
    end

    properties
        LabelWidth = 50
        LabelHeight = 20
        EditBoxWidth = 50
        EditBoxHeight = 20

        TextHeight = 15
        VerticalSpacing = 2
        LabelFontSize = 9
        EditBoxFontSize = 9
        EditBoxSpacing = 6
        EditBoxPosition

        StackVertically = false % When multiple controls are grouped, should they be stacked vertically or horizontally?
    end

    methods
        function obj = GemLabelEditBox(parent, text, tooltip, tag)
            obj = obj@GemVirtualPanel(parent);

            label_horizontal_alignment = 'center';

            obj.EditBox = GemEditBox(parent, tooltip);
            obj.AddChild(obj.EditBox);
            obj.AddEventListener(obj.EditBox, 'TextChanged', @obj.EditBoxCallback);
            obj.EditBox.HorizontalAlignment = 'left';
            obj.EditBox.FontSize = obj.EditBoxFontSize;
            text = [text ':'];

            obj.Text = GemText(parent, text, tooltip, tag);
            obj.Text.HorizontalAlignment = label_horizontal_alignment;
            obj.Text.FontSize = obj.LabelFontSize;
            obj.AddChild(obj.Text);
        end

        function Resize(obj, new_position)
            Resize@GemVirtualPanel(obj, new_position);

            control_width = round(new_position(3)/2);
            label_y_position = new_position(2) + obj.VerticalSpacing;
            label_x_position = new_position(1) + control_width - obj.LabelWidth;
            obj.Text.Resize([label_x_position, label_y_position, obj.LabelWidth, obj.LabelHeight]);

            edit_y_position = new_position(2) + obj.VerticalSpacing;
            edit_x_position = new_position(1) + control_width;
            obj.EditBox.Resize([edit_x_position, edit_y_position, obj.EditBoxWidth, obj.EditBoxHeight]);
        end

        function width = GetWidth(obj)
            width = obj.LabelWidth + obj.LabelWidth;
        end

        function height = GetRequestedHeight(obj, width)
            height = max(obj.EditBoxHeight, obj.LabelHeight);
        end

        function Enable(obj)
            Enable@GemVirtualPanel(obj);
            obj.Text.Enable;
            obj.EditBox.Enable;
        end

        function Disable(obj)
            Disable@GemVirtualPanel(obj);
            obj.Text.Disable;
            obj.EditBox.Disable;
        end

        function Select(obj, is_selected)
        end

    end

    methods (Access = protected)
        function EditBoxCallback(obj, hObject, ~)
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% demo of v_windows
close all;
wty=[1:5 7:15]; % list of window types to plot
for i=wty
    figure(i);
    v_windows(i,2520,'usw');
end
wty=[2,3,8,11]; % list of square root window types to plot
for i=wty
    figure(i+100);
    v_windows(i,2520,'usqw');
    title(['sqrt ' get(get(gca,'Title'),'String')]);
end
tilefigs;

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [ out ] = laplaceSP(a,rho)

if length(a)==1
    b= 1;
else
    b = a(2);
    a = a(1);
end
F = @(r)(Cp(rho).*besselj(0,rho*r));
out = 2*pi*(F(b) - F(a));

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%PRINT2ARRAY  Exports a figure to an image array
%
% Examples:
%   A = print2array
%   A = print2array(figure_handle)
%   A = print2array(figure_handle, resolution)
%   A = print2array(figure_handle, resolution, renderer)
%   [A bcol] = print2array(...)
%
% This function outputs a bitmap image of the given figure, at the desired
% resolution.
%
% If renderer is '-painters' then ghostcript needs to be installed. This
% can be downloaded from: http://www.ghostscript.com
%
% IN:
%   figure_handle - The handle of the figure to be exported. Default: gcf.
%   resolution - Resolution of the output, as a factor of screen
%                resolution. Default: 1.
%   renderer - string containing the renderer paramater to be passed to
%              print. Default: '-opengl'.
%
% OUT:
%   A - MxNx3 uint8 image of the figure.
%   bcol - 1x3 uint8 vector of the background color

% Copyright (C) Oliver Woodford 2008-2011

% 5/9/2011 Set EraseModes to normal when using opengl or zbuffer renderers.
% Thanks to Pawel Kocieniewski for reporting the issue.

% 21/9/2011 Bug fix: unit8 -> uint8!
% Thanks to Tobias Lamour for reporting the issue.

% 14/11/2011 Bug fix: stop using hardcopy(), as it interfered with figure
% size and erasemode settings. Makes it a bit slower, but more reliable.
% Thanks to Phil Trinh and Meelis Lootus for reporting the issues.

% 9/12/2011 Pass font path to ghostscript.

% 27/1/2012 Bug fix affecting painters rendering tall figures. Thanks to
% Ken Campbell for reporting it.

function [A bcol] = print2array(fig, res, renderer)
% Generate default input arguments, if needed
if nargin < 2
    res = 1;
    if nargin < 1
        fig = gcf;
    end
end
% Warn if output is large
old_mode = get(fig, 'Units');
set(fig, 'Units', 'pixels');
px = get(fig, 'Position');
set(fig, 'Units', old_mode);
npx = prod(px(3:4)*res)/1e6;
if npx > 30
    % 30M pixels or larger!
    warning('MATLAB:LargeImage', 'print2array generating a %.1fM pixel image. This could be slow and might also cause memory problems.', npx);
end
% Retrieve the background colour
bcol = get(fig, 'Color');
% Set the resolution parameter
res_str = ['-r' num2str(ceil(get(0, 'ScreenPixelsPerInch')*res))];
% Generate temporary file name
tmp_nam = [tempname '.tif'];
if nargin > 2 && strcmp(renderer, '-painters')
    % Print to eps file
    tmp_eps = [tempname '.eps'];
    print2eps(tmp_eps, fig, renderer, '-loose');
    try
        % Initialize the command to export to tiff using ghostscript
        cmd_str = ['-dEPSCrop -q -dNOPAUSE -dBATCH ' res_str ' -sDEVICE=tiff24nc'];
        % Set the font path
        fp = font_path();
        if ~isempty(fp)
            cmd_str = [cmd_str ' -sFONTPATH="' fp '"'];
        end
        % Add the filenames
        cmd_str = [cmd_str ' -sOutputFile="' tmp_nam '" "' tmp_eps '"'];
        % Execute the ghostscript command
        ghostscript(cmd_str);
    catch
        % Delete the intermediate file
        delete(tmp_eps);
        rethrow(lasterror);
    end
    % Delete the intermediate file
    delete(tmp_eps);
    % Read in the generated bitmap
    A = imread(tmp_nam);
    % Delete the temporary bitmap file
    delete(tmp_nam);
    % Set border pixels to the correct colour
    if isequal(bcol, 'none')
        bcol = [];
    elseif isequal(bcol, [1 1 1])
        bcol = uint8([255 255 255]);
    else
        for l = 1:size(A, 2)
            if ~all(reshape(A(:,l,:) == 255, [], 1))
                break;
            end
        end
        for r = size(A, 2):-1:l
            if ~all(reshape(A(:,r,:) == 255, [], 1))
                break;
            end
        end
        for t = 1:size(A, 1)
            if ~all(reshape(A(t,:,:) == 255, [], 1))
                break;
            end
        end
        for b = size(A, 1):-1:t
            if ~all(reshape(A(b,:,:) == 255, [], 1))
                break;
            end
        end
        bcol = median([reshape(A(:,[l r],:), [], size(A, 3)); reshape(A([t b],:,:), [], size(A, 3))], 1);
        for c = 1:size(A, 3)
            A(:,[1:l-1, r+1:end],c) = bcol(c);
            A([1:t-1, b+1:end],:,c) = bcol(c);
        end
    end
else
    if nargin < 3
        renderer = '-opengl';
    end
    err = false;
    % Set paper size
    old_mode = get(fig, 'PaperPositionMode');
    set(fig, 'PaperPositionMode', 'auto');
    try
        % Print to tiff file
        print(fig, renderer, res_str, '-dtiff', tmp_nam);
        % Read in the printed file
        A = imread(tmp_nam);
        % Delete the temporary file
        delete(tmp_nam);
    catch ex
        err = true;
    end
    % Reset paper size
    set(fig, 'PaperPositionMode', old_mode);
    % Throw any error that occurred
    if err
        rethrow(ex);
    end
    % Set the background color
    if isequal(bcol, 'none')
        bcol = [];
    else
        bcol = bcol * 255;
        if isequal(bcol, round(bcol))
            bcol = uint8(bcol);
        else
            bcol = squeeze(A(1,1,:));
        end
    end
end
% Check the output size is correct
if isequal(res, round(res))
    px = [px([4 3])*res 3];
    if ~isequal(size(A), px)
        % Correct the output size
        A = A(1:min(end,px(1)),1:min(end,px(2)),:);
    end
end
return

% Function to return (and create, where necessary) the font path
function fp = font_path()
fp = user_string('gs_font_path');
if ~isempty(fp)
    return
end
% Create the path
% Start with the default path
fp = getenv('GS_FONTPATH');
% Add on the typical directories for a given OS
if ispc
    if ~isempty(fp)
        fp = [fp ';'];
    end
    fp = [fp getenv('WINDIR') filesep 'Fonts'];
else
    if ~isempty(fp)
        fp = [fp ':'];
    end
    fp = [fp '/usr/share/fonts:/usr/local/share/fonts:/usr/share/fonts/X11:/usr/local/share/fonts/X11:/usr/share/fonts/truetype:/usr/local/share/fonts/truetype'];
end
user_string('gs_font_path', fp);
return

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function T_hat = run_tc_full(params)
  T = params.T;
  Idx = params.Idx;

  DIM = size(T);
  R = 1;
  T(T == 0) = 1e-3;
  T = tensor(T);
  T = T.*Idx;
  %T = sptensor(T);

  % Bayes CP algorithm for incomplete tensor and tensor completion
  [model] = BCPF_TC(T, 'obs', Idx, 'init', 'ml', ...
    'maxRank', max([DIM 2*R]), 'dimRed', 1, 'tol', 1e-6, 'maxiters', 100, 'verbose', 0);

  T_hat = double(model.X);
end

author: andrewssobral


text:
function out = cplexoptimget (options, name, default)
%%
% Purpose
% Retrieve optimization options values.
%
% Syntax
%    out = cplexoptimget(options,'param')
%    out = cplexoptimget(options,'param',default)
%
% Description
% val = cplexoptimget(options,'param') returns the value of the specified
% option in the optimization options structure options.
%
% val = cplexoptimget(options,'param',default) returns default if the
% specified option is not defined in the optimization options structure
% options.
%
%  See also cplexoptimset
%

% ---------------------------------------------------------------------------
% File: cplexoptimget.m
% ---------------------------------------------------------------------------
% Licensed Materials - Property of IBM
% 5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55
% Copyright IBM Corporation 2008, 2010. All Rights Reserved.
%
% US Government Users Restricted Rights - Use, duplication or
% disclosure restricted by GSA ADP Schedule Contract with
% IBM Corp.
% ---------------------------------------------------------------------------

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


function [price, opt, L] = PROJ_GMDB_DCA(proj_params, S_0, gmdb_params, r, q, modelInput)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for DCA-Style Garuanteed Minimum Withdraw Benefit (GMWB) using PROJ method
%
% Terminal Payoff:  Payoff(tau) = L*exp(g*tau) + (Gam(tau) - L*exp(g*tau))^+
%                      Gam(tau) = S_M * sum_{m=0}^M(alpha*gamma / S_m)
%                          tau  = time of death (discrete periods)
%
% Models Supported: Levy Processes, including jump diffusions and Black-Scholes model
% Returns: price of contract
%
% NOTE: this is the SLOW "Direct" version for testing purpose. In general, use PROJ_GMDB_DCA_Fast
%
% Author: Justin Lars Kirkby
% References: 1) Equity-Linked  Guaranteed Minimum Death Benefits with Dollar Cost Averaging, J.L.Kirkby & D.Nguyen, 2021
%
% ----------------------
% Contract/Model Params
% ----------------------
% S_0 = initial stock price (e.g. 100)
% r   = interest rate (e.g. 0.05)
% q   = dividend yield (e.g. 0.05)
% M   = number of subintervals of [0,T] (total of M+1 monitoring points in time grid, including S_0)
% gmdb_params = container of GMDB contract params, see below
% modelInput =  model inputs, see below
%
% ----------------------
% Numerical (PROJ) Params
% ----------------------
% proj_params = numerical params
%   proj_params.N = number of basis elements, e.g. N = 2^10
%   proj_params.L1 = gridwidth param, e.g. L1 = 8
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ------------------
% GMDB Contract Params
% ------------------
L = gmdb_params.L;   % Guarantee Level: Set L = -1 to use ATMF value for L
alpha = gmdb_params.alpha;  % Period premium payment, paid every dt time units
gamma = gmdb_params.gamma;  % Proportion of investment retained by policyholder (fee is 1-gamma)
contract_type = gmdb_params.contract_type;  % Contract type: 1 = GMDB, 2 = GMDB-RS (Ratchet strike)
p = gmdb_params.death_prob;  % death probability distribution, must be consistent with dt
g = gmdb_params.g;

% ------------------
% Model Inputs
% ------------------
dt = modelInput.dt;  % Time increment, premiums paid / underlying S is monitored every dt
phiR = modelInput.rnCHF;  % Risk neutral CHF for time period dt

Z = gen_func(-r, dt, p);
if g == 0
    Zrg = Z;
else
    Zrg = gen_func(-(r-g), dt, p);
end

if L == -1
    MF = gen_func(r - q - g, dt, p);
    Zg = gen_func(-g, dt, p);
    L = alpha * gamma * (exp((r-q)*dt)*MF - Zg) / (exp((r-q)*dt) - 1);
end

call = 1;
N = proj_params.N;

s = 0;

for n = 1 : length(p)

    M = n;
    T = n*dt;
    if contract_type == 2  % GMDB-RS
       W = S_0;
    else
       W = S_0*L*exp(g*T) / (alpha*gamma*(M+1));
    end

    pr_alpha = getTruncationAlpha(T, proj_params.L1, modelInput, proj_params.model);

    if n == 1
        W = 2*W - S_0;
        opt_v = 0.5*PROJ_European(3, N, 2*pr_alpha, r, q, T, S_0, W, call, phiR, modelInput.c1);
    else
        if contract_type == 1 || contract_type == 2  % GMDB / GMDB-RS
            ER = 0;
            opt_v = PROJ_Asian(N, pr_alpha, S_0, M, W, call, T, r, q, phiR, ER);
        elseif contract_type == 3  % European (for upper bound)
            opt_v = PROJ_European(3, N, 2*pr_alpha, r, q, T, S_0, W, call, modelInput.rnCHF_T, modelInput.c1);
            if r ~= 0  % else there is no multiplier)
                opt_v = opt_v * ((1 - exp(-r*(n+1)*dt)) / (1 - exp(-r*dt)))/(n+1);
            end
        elseif contract_type == 4  % geometric Asian (for lower bound)
            opt_v = PROJ_Geometric_Asian(N, pr_alpha, S_0, M, W, call, T, r, q, modelInput.rnSYMB);
        end
    end

    % fprintf('%.12f \n', opt_v);
    s = s + p(n) * (n + 1) * opt_v;

end

opt = s * alpha * gamma / S_0;
price = L*Zrg - alpha * (exp(r*dt) - Z) / (exp(r*dt) - 1) + opt;

end


text:
function [chi,omega_b,a_b,S] = lukfUpdate(chi,omega_b,a_b,...
    S,y,param,R,ParamFilter)
param.Pi = ParamFilter.Pi;
param.chiC = ParamFilter.chiC;

k = length(y);
q = length(S);
N_aug = q+k;
Rc = chol(kron(eye(k/2),R));
S_aug = blkdiag(S,Rc);

% scaled unsented transform
W0 = 1-N_aug/3;
Wj = (1-W0)/(2*N_aug);
gamma = sqrt(N_aug/(1-W0));
alpha = 1;
beta = 2;

% Compute transformed measurement
X = gamma*[zeros(N_aug,1) S_aug' -S_aug'];% sigma-points
Y = zeros(k,2*N_aug+1);
Y(:,1) = h(chi,zeros(q-6,1),param,zeros(N_aug-q,1));
for j = 2:2*N_aug+1
    xi_j = X([1:9 16:q],j);
    v_j = X(q+1:N_aug,j);
    Y(:,j) = h(chi,xi_j,param,v_j);
end
ybar = W0*Y(:,1) + Wj*sum(Y(:,2:end),2);% Measurement mean
Y(:,1) = sqrt(abs(W0+(1-alpha^2+beta)))*(Y(:,1)-ybar);
YY = sqrt(Wj)*(Y(:,2:2*N_aug+1)-ybar*ones(1,2*N_aug));
[~,Rs] = qr(YY');
Ss = Rs(1:k,1:k);
[Sy,~] = cholupdate(Ss,Y(:,1),'-'); % Sy'*Sy = Pyy
Pxy = zeros(q,k);
for j = 2:2*N_aug+1
    Pxy = Pxy + Wj*X(1:q,j)*(Y(:,j)-ybar)';
end

K = Pxy*Sy^-1*Sy'^-1; % Gain

xibar = K*(y-ybar);
omega_b = omega_b + xibar(10:12);
a_b = a_b + xibar(13:15);
xibar = xibar([1:9 16:q]);

% Covariance update
A = K*Sy';
for n = 1:k
    S = cholupdate(S,A(:,n),'-');
end

% Update mean state
chi = chi*exp_multiSE3(xibar);
J = xi2calJr(xibar);
S = S*J;
end

%--------------------------------------------------------------------------
function y = h(chi,xi,param,v)
Pi = param.Pi;
chiC = param.chiC;
RotC = chiC(1:3,1:3);
xC = chiC(1:3,4);

yAmers = param.yAmers;
NbAmers = length(yAmers);

chi_j = chi*exp_multiSE3(xi);
Rot = chi_j(1:3,1:3);
x = chi_j(1:3,5);
PosAmers = chi_j(1:3,6:end);
posAmers = PosAmers(:,yAmers);
z = Pi*( (Rot*RotC)'*(posAmers-kron(x,ones(1,NbAmers))) ...
    - kron(xC,ones(1,NbAmers)));
y = z(1:2,:)./z(3,:);
y = y(:) + v;
end

author: mbrossar


text:
classdef MimMarkerGreen < MimGuiPlugin
    % MimMarkerGreen. Gui Plugin for switching transparency of zero
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimMarkerGreen is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     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
        ButtonText = 'Green'
        SelectedText = 'Green'
        ToolTip = 'Place green markers'
        Category = 'Marker colour'
        Visibility = 'Dataset'
        Mode = 'Markers'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'single_marker.png'
        IconColour = GemMarkerPoint.DefaultColours{2}
        Location = 2
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.NewMarkerColour = 2;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && isequal(gui_app.GetCurrentModeName, MimModes.MarkerMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.NewMarkerColour == 2;
        end
    end
end


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


text:
% SPECTRAL EFFICIENCY CURVES - UPPER BOUNDS ON RATE

if(0) % set to (1) if you wish to run the converse_mc function

    n = 1e3;
    Pe = 1e-5;
    SNRdB = -15:7;
    uno = ones(size(SNRdB));tic

    % normal approximation - expected execution time 0.5 sec
    tic; rhoNA = 2*converse_mc(n*uno,Pe*uno,SNRdB,'normal'); % spectral efficiency
    EbN0NA = SNRdB-10*log10(rhoNA); toc
    % O(n^-2) approximation - expected execution time 10 sec
    tic; rhoPPV2 = 2*converse_mc(n*uno,Pe*uno,SNRdB,'On2'); % spectral efficiency
    EbN0PPV2 = SNRdB-10*log10(rhoPPV2); toc
    % O(n^-3) approximation - expected execution time 20 sec
    tic; rhoPPV3 = 2*converse_mc(n*uno,Pe*uno,SNRdB,'On3'); % spectral efficiency
    EbN0PPV3 = SNRdB-10*log10(rhoPPV3); toc
    % full approximation - expected execution time 10 min
    tic; rhoPPV = 2*converse_mc(n*uno,Pe*uno,SNRdB,'full'); % spectral efficiency
    EbN0PPV = SNRdB-10*log10(rhoPPV); toc

    % save('example_speff.mat')
else
    load('example_speff.mat')
end

% display results
close all
figure(1)
set(0,'defaulttextinterpreter','latex')
semilogy(EbN0PPV3,rhoPPV3,'--', EbN0PPV,rhoPPV,'-',EbN0PPV2,rhoPPV2,'--',EbN0NA,rhoNA,'-.x')
xlabel('SNR $E_b/N_0$')
ylabel('spectral efficiency $\rho$ [bit/s/Hz]')
title(['n = ' num2str(n) ', Pe = ', num2str(Pe)])
legend({},'interpreter','latex')
legend('$O(n^{-3}\,)$ approx $\quad$','full approx','$O(n^{-2}\,)$ approx',...
    'normal approx','Location','Best')
grid

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function z = r8vec_convolution ( m, x, n, y )

%*****************************************************************************80
%
%% R8VEC_CONVOLUTION returns the convolution of two R8VEC's.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    The I-th entry of the convolution can be formed by summing the products
%    that lie along the I-th diagonal of the following table:
%
%    Y3 : 3   4   5   6   7
%    Y2 : 2   3   4   5   6
%    Y1 : 1   2   3   4   5
%       +------------------
%        X1  X2  X3  X4  X5
%
%    which will result in:
%
%    Z = ( X1 * Y1,
%          X1 * Y2 + X2 * Y1,
%          X1 * Y3 + X2 * Y2 + X3 * Y1,
%                    X2 * Y3 + X3 * Y2 + X4 * Y1,
%                              X3 * Y3 + X4 * Y2 + X5 * Y1,
%                                        X4 * Y3 + X5 * Y2,
%                                                  X5 * Y3 )
%
%  Example:
%
%    Input:
%
%      X = (/ 1, 2, 3, 4 /)
%      Y = (/ -1, 5, 3 /)
%
%    Output:
%
%      Z = (/ -1, 3, 10, 17, 29, 12 /)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the dimension of X.
%
%    Input, real X(M), the first vector to be convolved.
%
%    Input, integer N, the dimension of Y.
%
%    Input, real Y(N), the second vector to be convolved.
%
%    Output, real Z(M+N-1), the convolution of X and Y.
%
  z(1:m+n-1,1) = 0.0;

  for j = 1 : n
    z(j:j+m-1,1) = z(j:j+m-1,1) + x(1:m,1) * y(j,1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%FEATRANK Feature ranking on individual performance for classification
%
% 	[I,F] = FEATRANK(A,CRIT,T)
% 	  I   = A*FEATRANK([],CRIT,T)
% 	  I   = A*FEATRANK(CRIT,T)
%
% INPUT
%   A      input dataset
%   CRIT   string name of a method or untrained mapping, default 'NN'
%   T      validation dataset (optional)
%
% OUTPUT
%   I      vector with sorted feature indices
%   F      vector with criteria values
%
% DESCRIPTION
% Feature ranking based on the training dataset A. CRIT determines
% the criterion used by the feature evaluation routine feateval. If
% the dataset T is given, it is used as test set for feateval. In I
% the features are returned in decreasing performance. In F the
% corresponding values of feateval are given. Default: crit='NN'.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FEATEVAL, FEATSELO, FEATSELB, FEATSELF,
% FEATSELP, FEATSELM

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function [I,F] = featrank(varargin)

argin = shiftargin(varargin,'char');
argin = setdefaults(argin,[],'NN',[]);
if mapping_task(argin,'definition')
  I = define_mapping(argin,'fixed');
else
  [a,crit,t] = deal(argin{:});
	[m,k,c] = getsize(a);
	F = zeros(1,k);
	isvaldfile(a,1,2); % at least 1 object per class, 2 classes
	a = testdatasize(a);
	iscomdset(a,t);

	if isempty(t)
		for j = 1:k
			F(j) = feateval(a(:,j),crit);
		end
	else
		% run the criterion on the validation set
		for j = 1:k
			F(j) = feateval(a(:,j),crit,t(:,j));
		end
	end

	[F,I] = sort(-F);
	F = -F;
end
return

author: marianux

repo: ecg-kit

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


text:
function scMove(Pos1,length,LinCol)
%scMove : Move along the transmission line on smith chart to transform impedance
%
%
%  SYNOPSIS:
%     Transform a point to another one as one moves along a lossless transmission line on
%     smith chart.
%
%     See also scDraw, scInv, scConCirc, scMatchCirc
%
%  SYNTAX:
%     scMove(Pos1, length, Direction, LinCol)
%
%  INPUT ARGUMENTS:
%     Pos1      : Starting point coordinates [r x] normalized
%     length    : Length of the transmission line normalized to wavelength
%     LinCol     : Color of the end-ray emanating from origin
%     Direction : Towards(+1) or away from the generator(-1) %%%%TO BE INCORPORATED YET.
%
%  OUTPUT ARGUMENT:
%     none
%
%  EXAMPLE:
%         The Command sequence
%         scDraw;
%         scMove([2 3],.3)
%         will draw a blank smith chart and will plot the point
%         corresponding to impedance (2+j*3)*Z_L and its transformed value
%         after moving 0.3*lambda towards the generator along the
%         transmission line
%
%     Mohammad Ashfaq - (31-05-2000)
%     Mohammad Ashfaq - (13-04-2006) Modified (example included)
%

if nargin == 2
    LinCol = 'm';
end

 r1       = Pos1(1);
 x1       = Pos1(2);
 [u1, v1] = scPOI(r1,x1);

 % MARK STARTING POINT
 scRay(Pos1);

 r        = sqrt(u1^2+v1^2);
 Theta1   = atan2(v1,u1);
 length1  = (1-Theta1/pi)/4;
 length2  = length1 + length;
 Theta2   = (1-4*length2)*pi;
 u2       = r * cos(Theta2);
 v2       = r * sin(Theta2);

 % MARK END POINT
 plot(u2, v2,'r*')

 Theta2d  = Theta2*180/pi;
 if (Theta2d<-180)
    Theta2d  = 360 + Theta2d;
 elseif (Theta2d>180)
    Theta2d  = 360 - Theta2d;
 end
 length2 = length2-floor(length2/.5)*.5;

 plot([0 u2],[0 v2],LinCol);
 plot([1.25*cos(Theta2) u2],[1.25*sin(Theta2) v2],'k');
 string = ['\theta=', num2str(Theta2d,'%3.2f'), '  l/\lambda=', num2str(length2, '%0.3f')];

 if abs(Theta2d)>90
    Theta2d  = Theta2d+180;
    h = text(1.25*cos(Theta2),1.25*sin(Theta2), string);
    set(h,'rotation', Theta2d, 'Fontsize', 8, 'HorizontalAlignment','right');
 else
    h = text(1.25*cos(Theta2),1.25*sin(Theta2), string);
    set(h,'rotation', Theta2d, 'Fontsize', 8);
 end


text:
function [efeaturesdj, adjlistdj] = mcmcGetAllEdgeData(spfeatures, imsegs, adjlist)

for f = 1:numel(imsegs)
    disp([num2str(f) ': getting edge data'])
    npertype = imsegs(f).nseg*4;
    [efeaturesdj{f}, adjlistdj{f}] = mcmcGetDisjointEdgeData(imsegs(f), spfeatures{f},  npertype, adjlist{f});
    disp(num2str([imsegs(f).nseg size(adjlist{f}, 1) size(adjlistdj{f},1)]))
end

author: zouchuhang

repo: LayoutNet


text:
classdef TestAllDemos < matlab.unittest.TestCase

    properties (TestParameter)
      demo = findDemos ();
    end

    methods (Test)
        function testDemo (testCase, demo)
            setup = 'pause off; warning off; ' ;
            try
                evalin ('base', [setup demo ' ; close all'])
                %evalc ([setup demo ' ; close all']);
            catch err
                testCase.verifyFail (err.message);
            end
        end
    end
end

function demos = findDemos ()

    vba_info = VBA_version();
    d = dir(fullfile(vba_info.path, '**', 'demo_*.m'));

    demos = {};
    for i = 1 : numel (d)
        demos{end+1} = d(i).name(1:end-2);
    end

end

author: MBB-team

repo: VBA-toolbox


text:
function y = size( x, varargin )
y = size( cvxaff( x ), varargin{:} );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function varargout = milpsubsref(varargin)
%MILPSUBSREF

switch class(varargin{1})
    case 'double'
        varargin{1} = double(varargin{1});
        for i = 1:length(varargin{2}.subs)
            varargin{2}.subs{i} = round(double(varargin{2}.subs{i}));
        end
        varargout{1} = subsref(varargin{:});

    case 'sdpvar'
        switch length(varargin{2}.subs)
            case 1
                index = varargin{2}.subs{1};
                if length(index) > 1
                    y = [];
                    for i = 1:length(index)
                        varargin{2}.subs{1} = index(i);
                        y = [y;yalmip('define',mfilename,varargin{:})];
                    end
                    % Figure out dims of variable
                    X = randn(size(varargin{1}));
                    varargin{2}.subs{1} = ones(size(index));
                    X = subsref(X,varargin{2});
                    y = reshape(y,size(X,1),size(X,2));
                    varargout{1} = y;
                else
                    varargout{1} = yalmip('define',mfilename,varargin{:});
                end
            case 2
                index1 = varargin{2}.subs{1};
                if isa(index1,'char')
                    index1 = 1:size(varargin{1},1);
                end
                index2 = varargin{2}.subs{2};
                if isa(index2,'char')
                    index2 = 1:size(varargin{1},2);
                end
                y = [];
                if length(index1)*length(index2) == 1
                    varargin{2}.subs{1} = index1;
                    varargin{2}.subs{2} = index2;
                    y = yalmip('define',mfilename,varargin{:});
                else
                    for i = 1:length(index1)
                        temp = [];
                        for j = 1:length(index2)
                            varargin{2}.subs{1} = index1(i);
                            varargin{2}.subs{2} = index2(j);
                            temp = [temp yalmip('define',mfilename,varargin{:})];
                        end
                        y = [y;temp];
                    end
                end
                % Figure out dims of variable
                X = randn(size(varargin{1}));
                varargin{2}.subs{1} = ones(size(index1));
                varargin{2}.subs{2} = ones(size(index2));
                X = subsref(X,varargin{2});
                y = reshape(y,size(X,1),size(X,2));
                varargout{1} = y;

            otherwise
                error('Only 1D and 2D variable subsref implemented');
        end

    case 'char'

        X = varargin{3};
        Y = varargin{2};
        if length(varargin{4}.subs) == 1
            X = X(:);
            i = varargin{4}.subs{1};
            M = length(X);
            m = 1;
            F = (integer(i)); % just to be sure
            d = binvar(length(X),1);
            [Mx,mx]=derivebounds(X);
            for j = m:M
                di = d(j);
                %                        F = F + (mx*(1-di) <= Y-X(j) <= Mx*(1-di));
                F = F + (-(max(Mx)-min(mx))*(1-di) <= Y-X(j) <= (max(Mx)-min(mx))*(1-di));
                F = F + (-(1+M-m)*(1-di) <= i-j <= (1+M-m)*(1-di));
            end
            F = F + (sum(d)==1);
        else
            i1 = varargin{4}.subs{1};
            i2 = varargin{4}.subs{2};
            M1 = size(X,1);
            M2 = size(X,2);
            m1 = 1;
            m2 = 1;
            if isa(i1,'sdpvar')
                F = (integer(i1)); % just to be sure
            end
            if isa(i2,'sdpvar')
                F = (integer(i2)); % just to be sure
            end
            d = binvar(size(X,1),size(X,2),'full');
            [Mx,mx]=derivebounds(X);
            for i = m1:M1
                for j = m2:M2
                    di = d(i,j);
                    F = F + (-(max(Mx)-min(mx))*(1-di) <= Y-X(i,j) <= (max(Mx)-min(mx))*(1-di));
                    if isa(i1,'sdpvar')
                        F = F + (-(1+M1-m1)*(1-di) <= i1-i <= (1+M1-m1)*(1-di));
                    elseif i1~=i
                        F = [F, sum(d(i,:))==0];
                    end
                    if isa(i2,'sdpvar')
                        F = F + (-(1+M2-m2)*(1-di) <= i2-j <= (1+M2-m2)*(1-di));
                    elseif i2~=j
                        F = [F, sum(d(:,j))==0];
                    end
                end
            end
            F = F + (sum(sum(d))==1);

        end
        varargout{1} = F;
        varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none','model','integer');
        varargout{3} = [X(:);i(:)];

    otherwise
        error('Strange type on first argument in SDPVAR/ABS');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function compile_mex( cudaRoot )
%compile_mex compiles all the MEX-functions included in this package

rootDir = fileparts( mfilename( 'fullpath' ) );

% image cropping on a GPU
cd( fullfile(rootDir, 'utils', 'cropRectanglesMex') );
if exist('cudaRoot', 'var')
	build_cropRectanglesMex( cudaRoot );
else
	build_cropRectanglesMex;
end
cd(rootDir);

% computation of the min-marginals
cd( fullfile(rootDir, 'pairwiseModel', 'computeMinMarginalsBinaryMex') );
build_computeMinMarginalsBinaryMex;
cd(rootDir);

% brute force energy minimization
cd( fullfile(rootDir, 'pairwiseModel', 'energyMinimization', 'bruteForceBinaryPairwiseMex') );
build_bruteForceBinaryPairwiseMex;
cd( rootDir );

% QPBO
cd( fullfile(rootDir, 'pairwiseModel', 'energyMinimization', 'qpboMex') );
build_qpboMex;
cd( rootDir );

% TRW-S
cd( fullfile(rootDir, 'pairwiseModel', 'energyMinimization', 'trwsMex') );
build_trwsMex;
cd( rootDir );

end

author: aosokin


text:
function asa047_test02 ( )

%*****************************************************************************80
%
%% ASA047_TEST02 demonstrates the use of NELMIN on POWELL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2008
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA047_TEST02\n' );
  fprintf ( 1, '  Apply NELMIN to POWELL quartic function.\n' );

  start(1:n) =  [ 3.0; - 1.0; 0.0; 1.0 ];
  reqmin = 1.0E-08;
  step(1:4) = 1.0;
  konvge = 10;
  kcount = 500;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Starting point X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %f\n', start(i) );
  end

  ynewlo = powell ( start );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X) = %f\n', ynewlo );

  [ xmin, ynewlo, icount, numres, ifault ] = nelmin ( @powell, n, start, ...
  reqmin, step, konvge, kcount );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Return code IFAULT = %d\n', ifault );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate of minimizing value X*:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %f\n', xmin(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %f\n', ynewlo );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', icount );
  fprintf ( 1, '  Number of restarts =   %d\n', numres );

  return
end


text:
function [xc,yc,radius] = fit_circle_to_points(x,y)

x = x(:);
y = y(:);

% solve for parameters a, b, and c in the least-squares sense by
% using the backslash operator
abc = [x y ones(numel(x),1)] \ -(x.^2+y.^2);
a = abc(1); b = abc(2); c = abc(3);

% calculate the location of the center and the radius
xc = -a/2;
yc = -b/2;
radius  =  sqrt((xc^2+yc^2)-c);

author: kristinbranson

repo: JAABA


text:
%IMPEXP_GID Import/export grid in GID msh format.
%
%   [ GRID ] = IMPEXP_GID( FILE_NAME, MODE, GRID, FID_LOG ) Import or export of
%   a GiD mesh ASCII format .msh file. FILE_NAME is a string specifying the (root)
%   file name to process. MODE can either be a string indicating import (no boundary
%   reconstruction), import_bdr (boundary reconstruction with gridbdr), or export.
%   GRID is a grid struct input for export and output for import. FID_LOG is an
%   optional log file handle for message output (negative for gui output or empty
%   for no output).

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function testFinalModelsAerts_HN(pathExperiment,fSetNames)
% - paramCells: One cell of parameters for each fSet
% - nameCells: One cell of names for each fSet

startpath = pwd;

cd(pathExperiment), load('training'), load('testing')
cd('FINAL_MODELS'), pathFinalModels = pwd;
nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
nFset = numel(fSetNames);

for o = 1:nOutcomes
    for f = 1:nFset
        results = []; results = struct;
        cd(fullfile(pathFinalModels,nameOutcomes{o},[fSetNames{f},'sign']))
        load('finalModel'), load('coeff'), load('response'), order = 4;
        results.model.Name = finalModel.Name; results.model.coeff = coeff; results.model.order = order;
        results.trainData.data = finalModel.Data; results.trainData.response = response; results.trainData.outcome = training.(nameOutcomes{o}).outcome;

        % TESTING MODEL
        outcome = testing.(nameOutcomes{o}).outcome;
        resp = zeros(numel(outcome),1);
        data = zeros(numel(outcome),order);
        for n = 1:order
            data(:,n) = testing.(nameOutcomes{o}).sign.(fSetNames{f})(:,n);
            resp = resp + coeff(n)*data(:,n);
        end
        resp = resp + coeff(end);
        results.testData.data = data; results.testData.response = resp; results.testData.outcome = outcome;
        [AUC,sensitivity,specificity,accuracy] = calcPerformMetrics(resp,outcome,0);
        results.AUC = AUC;
        results.Sensitivity = sensitivity;
        results.Specificity = specificity;
        results.Accuracy = accuracy;
        cd(pathExperiment), save(['testResultsLR_',[fSetNames{f},'sign'],'_',nameOutcomes{o}],'results')
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function u=transsplit(u0,a,b,x,y,T,nstep)
%%
% Computes an approximation to the solution of
%
% $$ u_t+  a u_x + b u_y=0 $$
%
% with initial data u0 on the rectangle spanned by the vectors x and y. u0 must be
% of the size |[length(x),length(y)]|.
%
% Output is a 3d matrix of size |[lentgth(x),length(y),((T/dt)+1)]|.
%
%% Initial setup
dt=T/nstep; Nt=ceil(T/dt); dt=T/Nt;
dx=x(2)-x(1); dy=y(2)-y(1);
u=zeros(length(x),length(y),Nt+1);

%% Solving by Strang splitting
u(:,:,1) = u0;
u1 = transport(b,u(:,:,1),0.5*dt,dx,1);   % Transport in the x-direction
for i=2:Nt,
   u(:,:,i) = transport(a,u1,      dt,dy,2);   % Transport in the y-direction
   u1       = transport(b,u(:,:,i),dt,dx,1);   % Transport in the x-direction
end
u(:,:,Nt+1) = transport(a,u1,            dt,dy,2);   % Transport in the y-direction
u(:,:,Nt+1) = transport(b,u(:,:,Nt+1),0.5*dt,dx,1);   % Transport in the x-direction

repo: Aero-matlab


text:
function result = is_symmetric_matrix(A)
%
% Symmetric Matrices
%
% is_symmetric_matrix(A) determines if the matrix A is a symmetric
% matrix. An error is returned if a matrix that is not square is attempted
% to be determined for symmetry.

matrix_size = size(A);

m = matrix_size(1,1);
n = matrix_size(1,2);

if m ~= n
    error('Only square matrices can be symmetric.');
else
    if A == A'
        result = 1;
    else
        result = 0;
    end
end

author: xiuyechen

repo: FishExplorer


text:
function varargout = LabBook(varargin)
%The labbook takes screenshots of the current matlab figure and generates
%html to display a collection of shots + annotations.
%
% LM:   4 Feb 2008   KU    Fixed distortion in screenshot thumbnails (function 'smartresize')
%
% 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/>.

persistent lastCapturedImageName;
persistent currentLabBook;

thumbnailWidth = 300;

if(nargin > 0)
	switch(upper(varargin{1}))

        case 'CAPTURE' %Capture passed figure, thumbnail and return both, [image, thumb]
            if(nargin > 1)
                figToCapture = varargin{2};
                F = getframe(figToCapture);
                [y,x,color] = size(F.cdata);
                resizeFactor = thumbnailWidth/x;
                thumb = F;
                thumb.cdata = imresize(F.cdata,resizeFactor,'nearest');
                varargout{1} = F;
                varargout{2} = thumb;
                return;
            else
                error('LabBook: No figure specified.');
            end

        case 'SAVE' %Save the passed image, thumbnail, and annotation to the passed filename, with the proper extentions.
            if(nargin == 5)
                if ~exist (currentLabBook,'dir')
                    error('LabBook: Specified LabBook Directory does not exist');
                    return;
                end
                large = varargin{2};
                thumb = varargin{3};
                imageName = varargin{4};
                annotation = varargin{5};
                filename = fullfile(currentLabBook, imageName);
                imwrite(large.cdata, [filename '.png'], 'png');
        %        imwrite(thumb.cdata, [filename 'Thumb' '.png'], 'png');
                fid = fopen([filename '.txt'],'w');
                fwrite(fid, annotation, 'char');
                fclose(fid);
                varargout{1} = [filename '.png'];
            else
                error('LabBook: Improper Number of Arguments');
            end

        case 'ANNOTATE' %Annotate figure stored in filename with text in annotation.  LabBook('Annotate', filename, annotation);
            if(nargin > 1)
                filename = varargin{2};
            else
                filename = lastCapturedImageName
            end

            if(nargin > 2)
                annotation = varargin{3};
            else
                error('No annotation specified.');
            end
            fid = fopen([filename '.txt'],'w');
            fwrite(fid, annotation, 'char');
            fclose(fid);

        case 'GENERATEREPORT' %generate an html report
            %Ideas: an Index.html page, using thumbs, 3x2 per page
            %w/scrolling, to display all possible images.  Thumbs are
            %linked to the corresponding large images.  When large images
            %is called up, txt from annotation report is displayed in a
            %frame at bottom of doc.
            if isempty(currentLabBook)
                currentLabBook = fullfile(getCERRPath, 'LabBook', 'Default');
            end
            fid = fopen(fullfile(currentLabBook, 'index.html'), 'w');
            %html = ['

 \r'...
					'\r'...
					'\r'...
					'
\r'...
					'
\r'...
					'\r'...
					'
\r'...
					'\r'...
				'
\r'...
				'
\r'...
				'
\r'...
				'
\r'...
				'	\r'...
				'

Back
Filename
	Prev Next

\r'...
				'
\r'...
				'	\r'...
				'

Comment:  \r'...
				'

\r'...
				'
\r'...
				'\r\'...
				''];

       %%%%%%

            fprintf(fid, html, 'char');
%            fprintf(fid, ['   document.write(''
'');for (var i = 0; i < (aImage.length); i = i + 1){	if (i %% 3 == 0)	{ 			document.write(''

'');   	}	document.write('''');	document.write("
" + aImage[i] + ".png");	document.write(''	'');}document.write(''
'');</script>

Caption

  ']);
            fclose(fid);

        case 'DISPLAYREPORT'
            web(fullfile(currentLabBook, 'index.html'), '-browser');

        case 'GETCURRENTBOOK'
            if isempty(currentLabBook) | ~exist(currentLabBook,'dir')
                currentLabBook = fullfile(getCERRPath, 'LabBook', 'Default', '');
                [success,msg,msgid] = mkdir(fullfile(getCERRPath, 'LabBook', ''), 'Default');
                if ~success
                    error('Unable to create default labbook directory.');
                end
            end
            varargout{1} = currentLabBook;
            return;
        case 'SETCURRENTBOOK'
            if(nargin > 1)
                currentLabBook = varargin{2};
            else
                error('No Book specified');
            end
            return;
        case 'CHECKDIRECTORY'
            if(nargin == 2)
                if(~exist(varargin{2}, 'dir'))
                    error('Specified LabBook Directory does not exist');
                end
            else
                error('Incorrect number of checkdirectory Parameters');
            end
	end
else
    warning('not enough parameters');
end

author: cerr

repo: CERR


text:
classdef SVMData < handle
    properties
        data
        conditions
        labels
        trials
        nTrials
        groups
        groupLabels
    end

    methods
        %%
        function obj = SVMData(raw)
        % SVMData(raw)
        %   Constructor for SVMData.
        %
        % @param rawDataClass raw
        % @return SVMData obj
        %
            if (exist('raw','var'))
                type = class(raw);
                switch (type)
                    case {'TimeSeriesData'}
                        obj.(['convert' type])(raw);
                    otherwise
                        fprintf(1,'\nUnrecognized data type.\n');
                end
            else
                fprintf(1,'\nObject created.  No data loaded.\n');
            end
        end

        %%
        function plot(obj, selectedData, color)
        % plot(selectedData, colorMap)
        %   Plots N rows of images depicting the data values with a
        %   colormap for ease of visual inspection.
        %
        %   Usage (Function):
        %       obj.plot({2, 3, 'first'; 3:4, 7:9, 'second'}, 'autumn');
        %       The first row will be condition 2, trial 3, label 'first',
        %       while the second row will be the mean of conditions 3 and
        %       4, trials 7 through 9, label 'second'.
        %
        %   Usage (Plot):
        %       Click and drag the mouse upon the plot to zoom in.  Click
        %       without dragging to zoom out one step from where you were.
        %
        % @param Nx2_cellarray<vector<int>> selectedData
        % @param string colorMap
        % @return void
        %
            figureHandle = figure;
            isDataSelected = (exist('selectedData','var'));
            if (isDataSelected)
                nRows = size(selectedData, 1);
                dataToPlot = zeros(nRows,size(obj.data,2));
                for i = 1:nRows
                    dataToPlot(i,:) = mean(obj.data(...
                        ismember(obj.conditions, selectedData{i,1}) ...
                        & ismember(obj.trials, selectedData{i,2}),:),1);
                end
            else
                nRows = length(obj.labels);
                dataToPlot = zeros(nRows,size(obj.data,2));
                for i = 1:nRows
                    dataToPlot(i,:) = mean(obj.data(...
                        ismember(obj.conditions, i) ...
                        & ismember(obj.trials, 1:obj.nTrials(i)),:),1);
                end
            end

            if (~exist('color','var')), color = 'bone'; end
            colormap(color);

            imagesc(dataToPlot);
            set(gca, 'YTick', 1:nRows);
            if (isDataSelected)
                set(gca, 'YTickLabel', selectedData(:,3));
            else
                set(gca,'YTickLabel', obj.labels);
            end

            startLine = line(0,0);
            endLine = line(0,0);
            xClick = [0, size(dataToPlot,2)];
            xLimits = xClick;
            clickHistory = [0, size(dataToPlot,2)];
            set(figureHandle,...
                'WindowButtonUpFcn', @stopDragCallback, ...
                'WindowButtonDownFcn', @startDragCallback);

            function startDragCallback(varargin)
                point = get(gca,'CurrentPoint');
                xValue = point(1);
                if (xValue > 0 || xValue < size(dataToPlot,2))
                    set(startLine, ...
                        'XData', [xValue xValue], ...
                        'YData', [0 (nRows + .5)]);
                    set(endLine, ...
                        'XData', [xValue xValue], ...
                        'YData', [0 (nRows + .5)]);
                    xClick = ones(1,2)*round(xValue);
                    set(figureHandle, 'WindowButtonMotionFcn', @draggingCallback);
                end
            end

            function stopDragCallback(varargin)
                set(figureHandle, 'WindowButtonMotionFcn', '');
                set(startLine, 'XData', 0, 'YData', 0);
                set(endLine, 'XData', 0, 'YData', 0);
                if (xClick(1) < xClick(2))
                    clickHistory = [clickHistory; xClick];
                elseif (xClick(1) > xClick(2))
                    xClick = fliplr(xClick);
                    clickHistory = [clickHistory; xClick];
                elseif (xClick(1) == xClick(2))
                    lastClick = size(clickHistory,1);
                    if (lastClick > 1)
                        clickHistory(lastClick,:) = [];
                        xClick = clickHistory(lastClick - 1,:);
                    elseif (lastClick == 1)
                        xClick = clickHistory(1,:);
                    end
                end
                xLimits = xClick;
                set(gca, 'XLim', xClick);

            end

            function draggingCallback(varargin)
                point = get(gca,'CurrentPoint');
                xValue = point(1);
                if (xValue > xLimits(1) && xValue <= xLimits(2))
                    set(endLine, ...
                        'XData', [xValue xValue]);
                    xClick(2) = round(xValue);
                end

            end
        end

        %%
        function convertTimeSeriesData(obj, raw)
        % convertTimeSeriesData(raw)
        %   Read a TimeSeriesData object and get the data into SVM format.
        %
        % @param TimeSeriesData raw
        % @return void
        %
                if (isempty(raw.analysis))
                    raw.loadTimeCourses();
                end
                nVoxels                     = length(raw.analysis);
                [nTimePoints nConds]        = size(raw.analysis(1).meanTcs);
                obj.labels                  = strtrim(raw.analysis(1).labels)';
                obj.nTrials                 = zeros(1,nConds);

                for i = 1:nConds
                    obj.nTrials(i) = length(raw.analysis(1).allTcs(1,~isnan(raw.analysis(1).allTcs(1,:,i)),i));
                end

                temp = zeros([size(raw.analysis(1).allTcs) nVoxels]);
                for i = 1:nVoxels
                    temp(:,:,:,i) = raw.analysis(i).allTcs;
                end
                temp = permute(temp, [2 1 3 4]);

                totalTrials              = sum(obj.nTrials);
                obj.data                 = zeros(totalTrials,nVoxels*nTimePoints);
                obj.conditions           = zeros(totalTrials, 1);
                obj.trials               = zeros(totalTrials, 1);
                obj.groups               = zeros(totalTrials, 1);

                for i = 1:nConds
                    startIndex    = sum(obj.nTrials(1 : (i - 1))) + 1;
                    endIndex      = startIndex + obj.nTrials(i) - 1;
                    obj.data(startIndex : endIndex, :) = ...
                        reshape(temp(1 : obj.nTrials(i), :, i, :), obj.nTrials(i), nVoxels * nTimePoints);
                    obj.conditions(startIndex : endIndex) = i;
                    obj.trials(startIndex : endIndex) = 1:obj.nTrials(i);
                end
        end

        %% I'm unhappy with the way labels and groups are handled...
        function group(obj, conditions, groupName)
        % group(conditions, groupName)
        %   Groups a set of conditions for SVM training.
        %
        % @param vector<int> conditions
        % @param string groupName
        % @return void
        %
            if (~exist('conditions','var'))
                conditions = obj.promptUser('Enter condition(s): ', {1:length(obj.labels)});
            end

            if (length(conditions) == 1)
                tellUser('Please select more than one condition when grouping.');
                return;
            end

            % Checks if we've already grouped these conditions
            if (sum(ismember(obj.conditions(find(obj.groups)),conditions)) ~= 0)
                tellUser(['One or more of the specified conditions have already '...
                    'been grouped.\nUngroup and rerun, or choose different conditions.']);
                return;
            end

            targetIndices = ismember(obj.conditions, conditions);

            obj.groups(targetIndices) = conditions(1);

            if (exist('groupName','var'))
                obj.groupLabels = [obj.groupLabels; {conditions(1), groupName}];
            end

        end

        %% Ditto from the groupSVM function
        function ungroup(obj)
        % ungroup()
        %   Run to remove groupings of training groups.
        %
        % @return void
        %
            obj.groups(:,1) = 0;
        end

        %%
        function clonedObject = clone(obj)
        % clonedObject = clone()
        %   Generate a deep copy of the object.
        %
        % @return SVMData clonedObject
        %
            clonedObject = SVMData();
            classInfo = metaclass(obj);
            nProperties = length(classInfo.Properties);
            for i = 1:nProperties
                if (~classInfo.Properties{i}.Constant)
                    name = classInfo.Properties{i}.Name;
                    clonedObject.(name) = obj.(name);
                end
            end
        end
    end

    methods (Static)
        %%
        function userResponse = promptUser(userMessage, acceptableInput)
        % userResponse = promptUser(userMessage, acceptableInput)
        %   Prompt user with message and limit acceptable responses.
        %
        %   Usage:
        %       resp = promptUser('Give me 1, 2, or 3: ', [1 2 3]);
        %       resp = promptUser('What is your name? ', 'string');
        %       resp = promptUser('How about a vector?', {1:8});
        %
        % @param string userMessage
        % @param string|vector<int>|cellarray acceptableInput
        % @return string|int|vector<int> userResponse
        %
            while (true)
                try
                    userResponse = input(userMessage);
                catch %#ok<*CTCH>
                    tellUser('Erroneous input: Could not process input.');
                    continue;
                end

                if (isa(acceptableInput, 'char'))
                    if (isa(userResponse, 'char'))
                        break;
                    else
                        tellUser('Erroneous input: Please enter a string.');
                    end
                elseif (isa(acceptableInput,'cell'))
                    if (~isa(userResponse, 'char')) && ...
                            (length(intersect(acceptableInput{1}, userResponse)) ...
                            == length(userResponse))
                        break;
                    else
                        tellUser(['Erroneous input: One or more of the ' ...
                            'integers entered was not valid.']);
                    end
                else
                    if (~isa(userResponse, 'char') && ...
                            (~sum(acceptableInput==userResponse)==0))
                        break;
                    else
                        tellUser('Erroneous input: Please enter a valid integer.');
                    end
                end
            end
        end
    end
end

%%
function tellUser(message)
% tellUser(message)
%   Shorthand to print a message to the user with some convenient
%   line breaks.
%
% @param string message
% @return void
%
    fprintf(1,['\n' message '\n']);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function se = cross3d
%CROSS3D Return a 3D structuring element with cross shape
%
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 30/09/2004.
%

%   HISTORY

se = zeros([3 3 3]);
se(1,2,2) = 1;
se(3,2,2) = 1;
se(2,:,:) = [0 1 0;1 1 1;0 1 0];

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function h = huber_pot(t, d)
%| function h = huber_pot(t, d)
%| huber potential function

if nargin < 1, ir_usage, end
if nargin < 2, huber_pot_test, return, end

h = t.^2 / 2;
ii = abs(t) > d;
h(ii) = d * abs(t(ii)) - d.^2/2;

function huber_pot_test
t = linspace(-9,9,101);
delta = 4;
plot(t, huber_pot(t, delta), '-', delta, huber_pot(delta, delta), 'o')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [info,meas] = fiff_read_meas_info(source,tree)
%
% [info,meas] = fiff_read_meas_info(source,tree)
%
% Read the measurement info
%
% If tree is specified, source is assumed to be an open file id,
% otherwise a the name of the file to read. If tree is missing, the
% meas output argument should not be specified.
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.14  2009/03/31 01:12:30  msh
%   Improved ID handling
%
%   Revision 1.13  2008/04/16 22:24:57  msh
%   Added megacq parameters to the measurement info
%
%   Revision 1.12  2008/03/13 19:18:06  msh
%   Read and write FIFF_MEAS_DATE from/to FIFFB_MEAS_INFO as appropriate
%
%   Revision 1.11  2007/11/13 10:55:32  msh
%   Specify the second argument to all calls to the exist function
%
%   Revision 1.10  2006/09/24 18:52:43  msh
%   Added FIFFV_REF_MEG_CH to fiff_define_constants.
%   Added coord_frame field to dig point structure.
%
%   Revision 1.9  2006/09/08 19:27:13  msh
%   Added KIT coil type to mne_load_coil_def
%   Allow reading of measurement info by specifying just a file name.
%
%   Revision 1.8  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.7  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.6  2006/04/14 15:49:49  msh
%   Improved the channel selection code and added ch_names to measurement info.
%
%   Revision 1.5  2006/04/13 23:09:46  msh
%   Further streamlining of the coordinate transformations.
%
%   Revision 1.4  2006/04/13 22:37:03  msh
%   Added head_head_trans field to info.
%
%   Revision 1.3  2006/04/13 17:05:45  msh
%   Added reading of bad channels to fiff_read_meas_info.m
%   Added mne_pick_channels_cov.m
%
%   Revision 1.2  2006/04/12 10:29:02  msh
%   Made evoked data writing compatible with the structures returned in reading.
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_read_meas_info';

if nargin ~= 2 & nargin ~= 1
    error(me,'Incorrect number of arguments');
end

if nargin == 1 & nargout == 2
    error(me,'meas output argument is not allowed with file name specified');
end

if nargin == 1
    [ fid, tree ] = fiff_open(source);
    open_here = true;
else
    fid = source;
    open_here = false;
end

%
%   Find the desired blocks
%
meas = fiff_dir_tree_find(tree,FIFF.FIFFB_MEAS);
if length(meas) == 0
    if open_here
        fclose(fid);
    end
    error(me,'Could not find measurement data');
end
%
meas_info = fiff_dir_tree_find(meas,FIFF.FIFFB_MEAS_INFO);
if length(meas_info) == 0
    if open_here
        fclose(fid);
    end
    error(me,'Could not find measurement info');
end
%
%   Read measurement info
%
dev_head_t=[];
ctf_head_t=[];
meas_date=[];
p = 0;
for k = 1:meas_info.nent
    kind = meas_info.dir(k).kind;
    pos  = meas_info.dir(k).pos;
    switch kind
        case FIFF.FIFF_NCHAN
            tag = fiff_read_tag(fid,pos);
            nchan = tag.data;
        case FIFF.FIFF_SFREQ
            tag = fiff_read_tag(fid,pos);
            sfreq = tag.data;
        case FIFF.FIFF_CH_INFO
            p = p+1;
            tag = fiff_read_tag(fid,pos);
            chs(p) = tag.data;
        case FIFF.FIFF_LOWPASS
            tag = fiff_read_tag(fid,pos);
            lowpass = tag.data;
        case FIFF.FIFF_HIGHPASS
            tag = fiff_read_tag(fid,pos);
            highpass = tag.data;
        case FIFF.FIFF_MEAS_DATE
            tag = fiff_read_tag(fid,pos);
            meas_date = tag.data;
        case FIFF.FIFF_COORD_TRANS
            tag = fiff_read_tag(fid,pos);
            cand = tag.data;
            if cand.from == FIFF.FIFFV_COORD_DEVICE && ...
                    cand.to == FIFF.FIFFV_COORD_HEAD
                dev_head_t = cand;
            elseif cand.from == FIFF.FIFFV_MNE_COORD_CTF_HEAD && ...
                    cand.to == FIFF.FIFFV_COORD_HEAD
                ctf_head_t = cand;
            elseif cand.from == FIFF.FIFFV_COORD_HEAD && ...
                    cand.to == FIFF.FIFFV_COORD_DEVICE
                dev_head_t = fiff_invert_transform(cand);
            end
    end
end
%
%   Check that we have everything we need
%
if ~exist('nchan','var')
    if open_here
        fclose(fid);
    end
    error(me,'Number of channels in not defined');
end
if ~exist('sfreq','var')
    if open_here
        fclose(fid);
    end
    error(me,'Sampling frequency is not defined');
end
if ~exist('chs','var')
    if open_here
        fclose(fid);
    end
    error(me,'Channel information not defined');
end
if length(chs) ~= nchan
    if open_here
        fclose(fid);
    end
    error(me,'Incorrect number of channel definitions found');
end

if isempty(dev_head_t) || isempty(ctf_head_t)
    hpi_result = fiff_dir_tree_find(meas_info,FIFF.FIFFB_HPI_RESULT);
    if length(hpi_result) == 1
        for k = 1:hpi_result.nent
            kind = hpi_result.dir(k).kind;
            pos  = hpi_result.dir(k).pos;
            if kind == FIFF.FIFF_COORD_TRANS
                tag = fiff_read_tag(fid,pos);
                cand = tag.data;
                if cand.from == FIFF.FIFFV_COORD_DEVICE && ...
                        cand.to == FIFF.FIFFV_COORD_HEAD
                    dev_head_t = cand;
                elseif cand.from == FIFF.FIFFV_MNE_COORD_CTF_HEAD && ...
                        cand.to == FIFF.FIFFV_COORD_HEAD
                    ctf_head_t = cand;
                end
            end
        end
    end
end
%
%   Locate the Polhemus data
%
isotrak = fiff_dir_tree_find(meas_info,FIFF.FIFFB_ISOTRAK);

dig=struct('kind',{},'ident',{},'r',{},'coord_frame',{});
coord_frame = FIFF.FIFFV_COORD_HEAD;
if length(isotrak) == 1
    p = 0;
    for k = 1:isotrak.nent
        kind = isotrak.dir(k).kind;
        pos  = isotrak.dir(k).pos;
        if kind == FIFF.FIFF_DIG_POINT
            p = p + 1;
            tag = fiff_read_tag(fid,pos);
            dig(p) = tag.data;
        else
            if kind == FIFF.FIFF_MNE_COORD_FRAME
                tag = fiff_read_tag(fid,pos);
                coord_frame = tag.data;
            elseif kind == FIFF.FIFF_COORD_TRANS
                tag = fiff_read_tag(fid,pos);
                dig_trans = tag.data;
            end
        end
    end
end
for k = 1:length(dig)
    dig(k).coord_frame = coord_frame;
end

if exist('dig_trans','var')
    if (dig_trans.from ~= coord_frame && dig_trans.to ~= coord_frame)
        clear('dig_trans');
    end
end
%
%   Locate the acquisition information
%
acqpars = fiff_dir_tree_find(meas_info,FIFF.FIFFB_DACQ_PARS);
acq_pars = [];
acq_stim = [];
if length(acqpars) == 1
    for k = 1:acqpars.nent
        kind = acqpars.dir(k).kind;
        pos  = acqpars.dir(k).pos;
        if kind == FIFF.FIFF_DACQ_PARS
            tag = fiff_read_tag(fid,pos);
            acq_pars = tag.data;
        else if kind == FIFF.FIFF_DACQ_STIM
                tag = fiff_read_tag(fid,pos);
                acq_stim = tag.data;
            end
        end
    end
end
%
%   Load the SSP data
%
projs = fiff_read_proj(fid,meas_info);
%
%   Load the CTF compensation data
%
comps = fiff_read_ctf_comp(fid,meas_info,chs);
%
%   Load the bad channel list
%
bads = fiff_read_bad_channels(fid,meas_info);
%
%   Put the data together
%
if ~isempty(tree.id)
    info.file_id = tree.id;
else
    info.file_id = [];
end
%
%  Make the most appropriate selection for the measurement id
%
if isempty(meas_info.parent_id)
    if isempty(meas_info.id)
        if isempty(meas.id)
            if isempty(meas.parent_id)
                info.meas_id = info.file_id;
            else
                info.meas_id = meas.parent_id;
            end
        else
            info.meas_id = meas.id;
        end
    else
        info.meas_id = meas_info.id;
    end
else
    info.meas_id = meas_info.parent_id;
end
if isempty(meas_date)
    info.meas_date = [ info.meas_id.secs info.meas_id.usecs ];
else
    info.meas_date = meas_date;
end
info.nchan     = nchan;
info.sfreq     = sfreq;
if exist('highpass','var')
    info.highpass = highpass;
else
    info.highpass = 0;
end
if exist('lowpass','var')
    info.lowpass = lowpass;
else
    info.lowpass = info.sfreq/2.0;
end
%
%   Add the channel information and make a list of channel names
%   for convenience
%
info.chs = chs;
for c = 1:info.nchan
    info.ch_names{c} = info.chs(c).ch_name;
end
%
%  Add the coordinate transformations
%
info.dev_head_t = dev_head_t;
info.ctf_head_t = ctf_head_t;
if ~isempty(info.dev_head_t) && ~isempty(info.ctf_head_t)
    info.dev_ctf_t    = info.dev_head_t;
    info.dev_ctf_t.to = info.ctf_head_t.from;
    info.dev_ctf_t.trans = inv(ctf_head_t.trans)*info.dev_ctf_t.trans;
else
    info.dev_ctf_t = [];
end
%
%   All kinds of auxliary stuff
%
info.dig   = dig;
if exist('dig_trans','var')
    info.dig_trans = dig_trans;
end
info.bads  = bads;
info.projs = projs;
info.comps = comps;
info.acq_pars = acq_pars;
info.acq_stim = acq_stim;

if open_here
    fclose(fid);
end

return;

    function [tag] = find_tag(node,findkind)

        for p = 1:node.nent
            if node.dir(p).kind == findkind
                tag = fiff_read_tag(fid,node.dir(p).pos);
                return;
            end
        end
        tag = [];
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Cb2n = ch_rotz(theta)
% 3D初等旋转， theta为旋转角度，rad
Cb2n = [cos(theta) -sin(theta) 0; sin(theta) cos(theta) 0; 0 0 1];

end

author: yandld


text:
%
% Example #1  Best-Fit Circle
%
%{

% pseudo-data setting
t=-pi:.01:pi;
x=sin(t);
y=cos(t);
xn=.05*randn(1,629);
yn=.05*randn(1,629);
x=x+xn;
y=y+yn;

% additional data for objective function
mydata.x=x;
mydata.y=y;

ds(1,'circlefit',mydata,10,3,-10,10,2000)

plot(x,y,'o','markersize',2);
hold on
plotcircle(globalminimizer(1),globalminimizer(2),globalminimizer(3),'r')
daspect([1 1 1])
shg

%}
function out=circlefit(X,mydata)

xdata=mydata.x;
ydata=mydata.y;

[N,~]=size(X);
out=ones(N,1); % pre-memory
for i=1:N
    x=X(i,:);
    a=x(1);
    b=x(2);
    r=x(3);

    out(i)=sum(abs((xdata-a).^2+(ydata-b).^2-repmat(r.^2,[size(xdata,1),1])));
end

% out---> Nx1 sized,  where N:size of superorganism (i.e., population size; pattern-matrix size)

author: Sable

repo: mcbench-benchmarks


text:
function CreateMovie(T, fileName)
[height, width, layer] = size(T{1});
aviobj = avifile(fileName, 'fps', 30, 'compression', 'None');
for i = 1:layer
    img(:,:,1) = T{1}(:,:,i);
    img(:,:,2) = T{2}(:,:,i);
    img(:,:,3) = T{3}(:,:,i);
    aviobj = addframe(aviobj, img);
end
aviobj = close(aviobj);

author: andrewssobral


text:
function out = nrlremap(x, a, c)
% NRLREMAP Lin-log style remap
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if nargin<3
    c = 220; % Output "knee" in lin-log curve
end
if nargin<2
    a = 1; % Scale factor of 99th percentile for input "knee"
end

A = abs(single(x));
Amin = min(A(:));
P99 = pctl(A(:),99);
b = (255-c)/log10((max(A(:))-Amin)/((a*P99)-Amin));

out = uint8(zeros(size(x)));
linear_region = (A<=a*P99);
out(linear_region) = (A(linear_region) - Amin)*c/((a*P99)-Amin);
out(~linear_region) = b*log10((A(~linear_region)-Amin)/((a*P99)-Amin))+c;

end

function ptile = pctl(v, p)
ptile = interp1(linspace(0.5/numel(v), 1-0.5/numel(v), numel(v))', sort(v(:)), p*0.01, 'spline');
end
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [y,t]=schmitt(x,thresh,minwid)
% Pass input signal X through a schmitt trigger
% SCHMITT(X,[LOW HIGH]) gives low and high thresholds. LOW and HIGH can be
% scalars or can be vectors specifiying different thresholds for each X element.
%
% SCHMITT(X,HYSTERESIS) specifies the thresholds as MAX-DELTA and MIN+DELTA where
% DELTA=(MAX-MIN)*(1-HYSTERESIS)/2 and MAX and MIN are the max and min of X.
% HYSTERESIS must be in the range 0 to 1 and represents the fraction of MAX-MIN between
% the two threshold values; default is 0.5.
%
% MINWID specifies the minimum width of a pulse (in samples). Pulses thinner than this will
% be ignored.
%
% Output Y takes values -1, +1 according to whether X<LOW or X>HIGH most recently.
% Y may be 0 for an initial segment if neither condition is initially true.
% For [Y,T]=SCHMITT(...) Y contains alternate +1 and -1 values and T contains
% the sample numbers at which X crossed the thresholds.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: schmitt.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
    minwid=0;
    if nargin<2
        thresh=0.5;
    end
end
if length(thresh)<2
    xmax=max(x);
    xmin=min(x);
    low=(xmax-xmin)*(1-thresh)/2;
    high=xmax-low;
    low=xmin+low;
else
    low=thresh(1);
    high=thresh(2);
end
c=(x>high)-(x<low);
c(2:end)=c(2:end).*(c(2:end)~=c(1:end-1));
t=find(c);
t(1+find(c(t(2:end))==c(t(1:end-1))))=[]; % remove duplicates
if minwid>=1
    t(t(2:end)-t(1:end-1)<minwid)=[];
    t(1+find(c(t(2:end))==c(t(1:end-1))))=[]; % remove duplicates
end
if nargout>1 y=c(t);
else
    y=zeros(size(c));
    if ~isempty(t)
        y(t)=2*c(t);
        y(t(1))=c(t(1));
        y=cumsum(y);
    end
end
if ~nargout
        xmax=max(x);
    xmin=min(x);
    if high-low<0.1*(xmax-xmin)
        high=xmax;
        low=xmin
    end
    plot([x(:) low+(y(:)+1)*(high-low)/2]);
end

author: decouples


text:
function r8_power_fast_test ( )

%*****************************************************************************80
%
%% R8_POWER_FAST_TEST tests R8_POWER_FAST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_POWER_FAST_TEST\n' );
  fprintf ( 1, '  R8_POWER_FAST computes R^P, economizing on\n' );
  fprintf ( 1, '  multiplications.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      R          P       R^P        Mults\n' );
  fprintf ( 1, '\n' );

  for i = -10 : 40

    r = 2.0;
    p = i;
    [ rp, mults ] = r8_power_fast ( r, p );
    fprintf ( 1, '  %12f  %5d  %12f  %5d\n', r, p, rp, mults );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Gausspoint,Gaussweight] = GaussQuadrature(ngl)
%-------------------------------------------------------------------
%  Purpose:
%     determine the integration points and weighting coefficients
%     of Gauss-Legendre quadrature for two-dimensional integration
%
%  Synopsis:
%     [point,weight]=GaussQuadrature(nglx,ngly)
%
%  Variable Description:
%     ngl - number of integration points
%     point - vector containing integration points
%     weight - vector containing weighting coefficients
%-------------------------------------------------------------------
%  initialization

   Gausspoint=zeros(ngl,1);
   Gaussweight=zeros(ngl,1);

%  corresponding integration points and weights
    % 2-point quadrature rule
    Gausspoint(1)=-0.577350269189626;
    Gausspoint(2)=-Gausspoint(1);
    Gaussweight(1)=1.0;
    Gaussweight(2)=Gaussweight(1);


text:
function stroud_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests BALL_MONOMIAL_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  r = 2.0;
  xc(1:n) = [ 0.0, 0.0, 0.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  For the integral of a monomial in a ball in ND:\n' );
  fprintf ( 1, '  BALL_MONOMIAL_ND approximates the integral.\n' );
  fprintf ( 1, '  BALL_F1_ND, which can handle general integrands,\n' );
  fprintf ( 1, '    will be used for comparison.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension N = %d\n', n );
  fprintf ( 1, '  Ball radius = %f\n', r )
  fprintf ( 1, '  Ball volume = %f\n', ball_volume_nd ( n, r ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        Rule:	   MONOMIAL	      BALL_F1_ND\n' );
  fprintf ( 1, '        F(X)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 4

    if ( i == 1 )
      string = '1';
      p = [ 0, 0, 0 ];
      result2 = ball_f1_nd ( 'mono_000_3d', n, xc, r );
    elseif ( i == 2 )
      string = 'xyz';
      p = [ 1, 1, 1 ];
      result2 = ball_f1_nd ( 'mono_111_3d', n, xc, r );
    elseif ( i == 3 )
      string = 'x^2z^2';
      p = [ 2, 0, 2 ];
      result2 = ball_f1_nd ( 'mono_202_3d', n, xc, r );
    elseif ( i == 4 )
      string = 'x^4y^2z^2';
      p = [ 4, 2, 2 ];
      result2 = ball_f1_nd ( 'mono_422_3d', n, xc, r );
    end

    result1 = ball_monomial_nd ( n, p, r );

    fprintf ( 1, '  %10s  %14f  %14f\n', string, result1, result2 );

  end

  return
end


text:
function sgmga_importance_to_aniso_test ( dim_num, importance )

%*****************************************************************************80
%
%% SGMGA_IMPORTANCE_TO_ANISO_TEST calls SGMGA_IMPORTANCE_TO_ANISO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_IMPORTANCE_TO_ANISO_TEST\n' );
  fprintf ( 1, '  Importances:\n' );
  for dim = 1 : dim_num
    fprintf ( 1, '  %12f', importance(dim) );
  end
  fprintf ( 1, '\n' );

  level_weight = sgmga_importance_to_aniso ( dim_num, importance );

  fprintf ( 1, '  Anistropic coefficients:\n' );
  for dim = 1 : dim_num
    fprintf ( 1, '  %12f', level_weight(dim) );
  end
  fprintf ( 1, '\n' );

  return
end