text:
function [trl] = ft_trialfun_twoclass_classification(cfg)

% FT_TRIALFUN_TWOCLASS_CLASSIFICATION can be used to train and test a real-time
% classifier in offline and online mode. It selects pieces of data in the two classes
% based on two trigger values. The first N occurences in each class are marked as
% training items. All subsequent occurrences are marked as test items.
%
% This function can be used in conjunction with FT_REALTIME_CLASSIFICATION. The
% configuration structure should contain
%   cfg.dataset              = string with the filename
%   cfg.trialfun             = 'ft_trialfun_twoclass_classification'
%   cfg.trialdef.numtrain    = number of training items, e.g. 20
%   cfg.trialdef.eventvalue1 = trigger value for the 1st class
%   cfg.trialdef.eventvalue2 = trigger value for the 2nd class
%   cfg.trialdef.eventtype   = string, e.g. 'trigger'
%   cfg.trialdef.prestim     = latency in seconds, e.g. 0.3
%   cfg.trialdef.poststim    = latency in seconds, e.g. 0.7
%
% See also FT_DEFINETRIAL, FT_TRIALFUN_GENERAL

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

% these are used to count the number of training items in each class
persistent numtrain1
persistent numtrain2

if isempty(numtrain1)
  numtrain1 = 0;
end
if isempty(numtrain2)
  numtrain2 = 0;
end

if isfield(cfg, 'hdr')
  hdr = cfg.hdr;
else
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat);
end

if isfield(cfg, 'event')
  event = cfg.event;
else
  event = ft_read_event(cfg.headerfile, 'headerformat', cfg.headerformat);
end

baseline = round(cfg.trialdef.prestim*hdr.Fs);
duration = round(cfg.trialdef.prestim*hdr.Fs + cfg.trialdef.poststim*hdr.Fs);

% make a subset of the interesting events
sel   = strcmp(cfg.trialdef.eventtype, {event.type});
event = event(sel);
num   = length(event);
trl   = zeros(num,5);

for i=1:num
  % determine the location of this trial in the data stream
  begsample = event(i).sample - baseline;
  endsample = begsample + duration - 1;
  offset    = baseline;
  % determine the class and wether this trial is eligeable for training
  if event(i).value==cfg.trialdef.eventvalue1
    class = 1;
    train = (numtrain1 < cfg.trialdef.numtrain); % boolean
    numtrain1 = numtrain1 + train;  % increment the counter
  elseif event(i).value==cfg.trialdef.eventvalue2
    class = 2;
    train = (numtrain2 < cfg.trialdef.numtrain); % boolean
    numtrain2 = numtrain2 + train;  % increment the counter
  else
    % the class is unknown and therefore irrelevant
    class = nan;
    train = false;
  end
  % remember this trial, the class and whether it should be used for training
  trl(i,:) = [begsample endsample offset class train];
end


text:
function [output] = getCompareBinsOfFluxes(xglc, model, samplo, samphi, metabolites)
% Compares the bins of fluxes between `samplo` and `samphi`,
% calls `compareBinsOfFluxes(xglc, model, sammin, sammax, metabolites).
% `samplo` and `samphi` each contain samples in `x.points`
%
% USAGE:
%
%    [output] = getCompareBinsOfFluxes(xglc, model, samplo, samphi, metabolites)
%
% INPUTS:
%    xglc:           sugar distribution
%    model:          model structure
%    samplo:         samples containing bins of fluxes
%    samphi:         samples containing bins of fluxes
%
% OPTIONAL INPUT:
%    metabolites:    list of metabolites
%
% OUTPUT:
%    output:         result of comparison
%
% .. Author: - Wing Choi 3/7/08

output = 0;

if (nargin < 4)
    disp '[output] = getCompareBinsOfFluxes(xglc,model,samplo,samphi,metabolites)'
    return;
end

if (nargin < 5)
    metabolites = [];
end

if (isempty(xglc))
    % random glucose
    %xglc = rand(64,1);
    %xglc = xglc/sum(xglc);
    %xglc = idv2cdv(6)*xglc;

    glucose = rand(8,1);
    glucose = glucose/sum(glucose);
    %glc = idv2cdv(6)*glc;

    % glc 1-6 = carbon 1-6
    % glc 7 = carbon 1+2 (really 5 and 6)
    % glc 8 = unlabeled
    % glc 9 = fully labeled
    glc = zeros(64,9);
    glc(1+1,1) = 1;
    glc(2+1,2) = 1;
    glc(4+1,3) = 1;
    glc(8+1,4) = 1;
    glc(16+1,5) = 1;
    glc(32+1,6) = 1;
    glc(32+16+1,7) = 1;
    glc(0+1,8) = 1;
    glc(63+1,9) = 1;

    xGlc = zeros(64,1);
    for i = 1:8
        xGlc = xGlc + glucose(i)*glc(:,i);
    end

    xglc = idv2cdv(6)*xGlc;

end

% specify the #fluxes and #bins to work with.
bins = [ 100 40 ; 125 32 ; 200 20 ; 250 16 ; 300 13 ; 400 10 ; 1000 4 ];

for i = 1:7
    nflux = bins(i,1);
    nbins = bins(i,2);
    ol = getBinsOfFluxes(samplo,nflux,nbins);
    ou = getBinsOfFluxes(samphi,nflux,nbins);
    ome = compareBinsOfFluxes(xglc,model,ol,ou,metabolites);
    name = sprintf('n%d',nbins);
    output.data(i,1:length(ome.totalz)) = ome.totalz';
    output.atotalz(i,1) = mean(ome.totalz);
    output.std(i,1) = std(ome.totalz);
end

% get the single bin of 4000 run for compare.
sl.samps(1,1).points = samplo.points(:,1:4000);
su.samps(1,1).points = samphi.points(:,1:4000);
ome = compareBinsOfFluxes(xglc,model,sl,su,metabolites);
output.data(8,1) = ome.totalz(1,1);
output.atotalz(8,1) = ome.totalz(1,1);
output.std(8,1) = 0;
%[totalz,zscore,mdvs] = compareMultSamp(xGlc,x,samps);

%[totalz] = compareTwoSamp(xGlc,x,samp1,samp2,metabolite);

return
end

author: opencobra

repo: cobratoolbox


text:
function singleMetaboliteThermoFigure(model, miliMolarStandard, metAbbr)
% Creates a vertical errorbar figure of the reactions involving the
% metabolite given by the abbreviation
%
% .. todo = fix this up
%
% USAGE:
%
%    singleMetaboliteThermoFigure(model, miliMolarStandard, metAbbr)
%
% INPUTS:
%    model:                structure
%    miliMolarStandard:
%    metAbbr:
%
% .. Author: - Ronan M.T. Fleming

[nMet,nRxn]=size(model.S); %%%%%%%%%vertical%%%%%%%%%%%%%
if miliMolarStandard
    dGt0Min=model.dGt0Min;
    dGt0Max=model.dGt0Max;
    for n=1:nRxn
        model.dGt0Min(n)=model.rxn(n).dGtmMMin;
        model.dGt0Max(n)=model.rxn(n).dGtmMMax;
    end
end
%forward possibly reversible
dGfGCforwardReversibleBool=directions.ChangeForwardReversibleBool_dGfGC;
%forward, probably reverse
fwdProbReverse=directions.ChangeForwardReversibleBool_dGfGC_byConcRHS | ...
    directions.ChangeForwardReversibleBool_dGfGC_bydGt0RHS | ...
    directions.ChangeForwardReversibleBool_dGfGC_byConc_No_dGt0ErrorRHS;

%plot regions 5,6,7, with region 6 shaded
X1=1:nRxn;
%dGrt0
Y0=(model.dGt0Min+model.dGt0Max)/2;
L0=Y0-model.dGt0Min;
U0=model.dGt0Max-Y0;
%dGrt
Y=(model.dGtMin+model.dGtMax)/2;
L=Y-model.dGtMin;
U=model.dGtMax-Y;
%find the amount of reactions with normal cumulative distribution over
%range of dGt0
P = normcdf(0,Y0,L0);
%sort by probability that a reaction is forward (puts any NaN first)
[tmp,xip]=sort(P,'descend');
%     only take the indices of the problematic reactions, but be sure to
%     take them in order of descending P
xip2=zeros(nnz(fwdProbReverse),1);
p=1;
for n=1:nRxn
    if fwdProbReverse(xip(n))
        xip2(p)=xip(n);
        p=p+1;
    end
end
xip=xip2;
X1=1:length(xip);

%replace the NaN due to zero st dev
nNaNpRHS=nnz(directions.ChangeForwardReversibleBool_dGfGC_byConc_No_dGt0ErrorRHS);
if nNaNpRHS~=nnz(isnan((P(fwdProbReverse))))
    warning('ExtraCategory');
end
%nans are first in the ordering of indexes
NaNPInd=xip(1:nNaNpRHS);
%sorts indices of the zero std dev met by their mean dG0t
[tmp,xipNaNPInd]=sort(Y0(NaNPInd));
%new ordering with NaNs sorted to right
xip=[xip(nNaNpRHS+1:end); NaNPInd(xipNaNPInd(1:nNaNpRHS))];

% close all
figure1 = figure;
% Create axes
axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
hold on;

%upper and lower X
minX=min(model.dGtMin(fwdProbReverse));
maxX=max(model.dGtMax(fwdProbReverse));
%baselines
PreversibleBar_byConcRHS=ones(1,nRxn)*minX;
%bar for 6
PreversibleBar_byConcRHS(directions.ChangeForwardReversibleBool_dGfGC_byConcRHS)=maxX;
bar_handle6=barh(X1,PreversibleBar_byConcRHS(xip),1,'BaseValue',minX,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

%dGrt errorbar
%   PLOTERR(X,Y,{LX,UX},{LY,UY}) plots the graph of vector X vs. vector Y
%   with error bars specified by LX and UX in horizontal direction and
%   with error bars specified by LY and UY in vertical direction.
%   H = PLOTERR(...) returns a vector of line handles in this order:
%      H(1) = handle to datapoints
%      H(2) = handle to errorbar y OR errorbar x if error y not specified
%      H(3) = handle to errorbar x if error y specified
% LineSpec=[''LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r'']
LineSpec='r.';
hE = ploterr(Y(xip),X1,{model.dGtMin(xip) model.dGtMax(xip)},[],LineSpec,'hhxy',0);
set(hE(2),'LineWidth',10);
% hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%dGrt0 errorbar on top and inside dGrt
LineSpec='b.';
hE = ploterr(Y0(xip),X1,{model.dGt0Min(xip) model.dGt0Max(xip)},[],LineSpec,'hhxy',0);
set(hE(2),'LineWidth',10);
% hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');

%mean dGrt0
plot(Y0(xip),X1,'.','LineStyle','none','Color',[0.3412 0.7961 0.1922]);
%zero line
%cumulative probability that reaction is really forward, assuming a
%normal distribution about the mean dGt0
[AX,H1,H2]=plotyy(zeros(1,length(X1)),X1,P(xip),X1);
set(gca,'XTick',[]);
set(AX(1),'XTickMode','manual','YTick',floor(minX/100)*200:50:maxX);%,'TickDirMode','manual','TickDir','out');
set(AX(2),'YTick',1:length(xip));
set(H1,'LineStyle','none')
set(H2,'LineStyle','-','LineWidth',2,'Color','k')
plot(zeros(1,length(X1)),X1,'w','LineWidth',2,'LineStyle','--');
%axis limits
set(AX(1),'FontSize',16,'YColor','k')
set(AX(1),'YTick',[],'YTickLabelMode','manual');
set(AX(1),'XTickMode','manual','XTick',floor(minX/100)*200:50:maxX);
set(get(AX(1),'Xlabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
set(get(AX(1),'Xlabel'),'FontSize',16)
set(AX(1),'Position',[0.05 0.1 0.4 0.8]);
axis(AX(1),[minX maxX 0.5 length(X1)+0.5 ]);
% axis(AX(1),[-250 250 0 length(X1)+0.5 ]);

set(AX(2),'FontSize',16,'YColor','k')
set(AX(2),'YTick',1:length(X1),'YTickLabelMode','manual');
set(AX(2),'XAxisLocation','top');
set(get(AX(2),'Xlabel'),'String','P(\Delta_{r}G^{\primem}<0)');
set(get(AX(2),'Xlabel'),'FontSize',16);
axis(AX(2),[0 1 0.5 length(X1)+0.5 ]);
set(AX(2),'Position',[0.05 0.1 0.4 0.8]);

% ylabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0) (blue)');
title('Qualitatively forward, but probably quantitatively reverse, using group contribution estimates of \Delta_{r}G^{\primem}.','FontSize',16)

for n=1:length(xip)
%     textLen(n)=length(model.rxn(xip(n)).officialName)+length(model.rxn(xip(n)).equation);
%     textLen(n)=length(model.rxn(xip(n)).equation);
   textLen(n)=length(model.rxn(xip(n)).officialName);
end
textLenMax=max(textLen);
for n=1:length(xip)
%     YTickLabel{n}=[model.rxn(xip(n)).officialName blanks(textLenMax-textLen(n)+5) model.rxn(xip(n)).equation];
    YTickLabel{n}=[model.rxn(xip(n)).officialName ' ' model.rxn(xip(n)).equation];
    %get rid of the cytoplasmic compartment abbreviations
    YTickLabel{n} = strrep(YTickLabel{n}, '[c]', '');

   % no idea why this does not work
   %     fprintf('%s\n',[model.rxn(xip(n)).officialName blanks(textLenMax-textLen(n)+5) model.rxn(xip(n)).equation]);
end
set(AX(2),'YTickLabel',YTickLabel,'FontSize',16);
saveas(figure1 ,'fwdProbReverseGC','fig');

%change back
if miliMolarStandard
    model.dGt0Min=dGt0Min;
    model.dGt0Max=dGt0Max;
end

%horizontal bar
horizontal=0;
if horizontal
    [nMet,nRxn]=size(model.S);
    if miliMolarStandard
        dGt0Min=model.dGt0Min;
        dGt0Max=model.dGt0Max;
        for n=1:nRxn
            model.dGt0Min(n)=model.rxn(n).dGtmMMin;
            model.dGt0Max(n)=model.rxn(n).dGtmMMax;
        end
    end
    %forward possibly reversible
    dGfGCforwardReversibleBool=directions.ChangeForwardReversibleBool_dGfGC;
    %forward, probably reverse
    fwdProbReverse=directions.ChangeForwardReversibleBool_dGfGC_byConcRHS | ...
        directions.ChangeForwardReversibleBool_dGfGC_bydGt0RHS | ...
        directions.ChangeForwardReversibleBool_dGfGC_byConc_No_dGt0ErrorRHS;

    %plot regions 5,6,7, with region 6 shaded
    X1=1:nRxn;
    %dGrt0
    Y0=(model.dGt0Min+model.dGt0Max)/2;
    L0=Y0-model.dGt0Min;
    U0=model.dGt0Max-Y0;
    %dGrt
    Y=(model.dGtMin+model.dGtMax)/2;
    L=Y-model.dGtMin;
    U=model.dGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');
    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(fwdProbReverse),1);
    p=1;
    for n=1:nRxn
        if fwdProbReverse(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    nNaNpRHS=nnz(directions.ChangeForwardReversibleBool_dGfGC_byConc_No_dGt0ErrorRHS);
    if nNaNpRHS~=nnz(isnan((P(fwdProbReverse))))
        warning('ExtraCategory');
    end
    %nans are first in the ordering of indexes
    NaNPInd=xip(1:nNaNpRHS);
    %sorts indices of the zero std dev met by their mean dG0t
    [tmp,xipNaNPInd]=sort(Y0(NaNPInd));
    %new ordering with NaNs sorted to right
    xip=[xip(nNaNpRHS+1:end); NaNPInd(xipNaNPInd(1:nNaNpRHS))];

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.dGtMin(fwdProbReverse));
    maxY=max(model.dGtMax(fwdProbReverse));
    %baselines
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    %bar for 6
    PreversibleBar_byConcRHS(directions.ChangeForwardReversibleBool_dGfGC_byConcRHS & directions.ChangeForwardReversibleBool_dGfGC_bydGt0RHS)=maxY;
    bar_handle6=bar(X1,PreversibleBar_byConcRHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
    % adjust error bar width
    hE_c=get(hE, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
    % adjust error bar width
    hE_c=get(hE2, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);

    %mean dGrt0
    plot(X1,Y0(xip),'.','LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:50:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none')
    set(H2,'LineStyle','-','LineWidth',2,'Color','k')
    plot(X1,zeros(1,length(X1)),'w','LineWidth',2,'LineStyle','--');
    %axis limits
    axis(AX(1),[0 length(X1) minY maxY])
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    axis(AX(2),[0 length(X1) 0 1])
    title({'Qualitatively forward, but quantitatively reversible' ; 'using group contribution estimates of \Delta_{r}G^{\primem}.'},'FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0) (blue)');
    saveas(figure1 ,'GCfwdProbFwd','fig');

    %change back
    if miliMolarStandard
        model.dGt0Min=dGt0Min;
        model.dGt0Max=dGt0Max;
    end
end


text:
function [ unique_num, seed ] = point_radial_tol_unique_count ( m, n, a, ...
  tol, seed )

%*****************************************************************************80
%
%% POINT_RADIAL_TOL_UNIQUE_COUNT counts the tolerably unique points.
%
%  Discussion:
%
%    The input data is an M x N array A, representing the M-dimensional
%    coordinates of N points.
%
%    The output is the number of tolerably unique points in the list.
%
%    This program performs the same task as POINT_TOL_UNIQUE_COUNT.
%    But that program is guaranteed to use N^2 comparisons.
%
%    It is hoped that this function, on the other hand, will tend
%    to use O(N) comparisons after an O(NLog(N)) sort.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 July 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows.
%
%    Input, integer N, the number of columns.
%
%    Input, real A(M,N), the array of N columns of data.
%
%    Input, real TOL, a tolerance for equality.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, integer UNIQUE_NUM, the number of tolerably
%    unique points.
%
  if ( n <= 0 )
    unique_num = 0;
    return
  end
%
%  Assign a base point Z randomly in the convex hull.
%
  [ w, seed ] = r8vec_uniform_01 ( n, seed );

  z = zeros ( m, 1 );
  for i = 1 : m
    z(i) = a(i,1:n) * w(1:n,1) / sum ( w(1:n) );
  end
%
%  Compute the radial distance R of each point to Z.
%
  r = zeros ( n, 1 );
  for j = 1 : n
    r(j) = sqrt ( sum ( ( a(1:m,j) - z(1:m) ).^2 ) );
  end
%
%  Implicitly sort the R array.
%
  indx = r8vec_sort_heap_index_a ( n, r );
%
%  To determine if a point I is tolerably unique, we only have to check
%  whether it is distinct from all points J such that R(I) <= R(J) <= R(J)+TOL.
%
  unique_num = 0;
  unique(1:n) = 1;

  for i = 1 : n

    if ( unique(indx(i)) )
%
%  Point INDX(I) is unique, in that no earlier point is near it.
%
      unique_num = unique_num + 1;
%
%  Look for later points which are close to point INDX(I)
%  in terms of R.
%
      hi = i;

      while ( hi < n )
        if ( r(indx(i)) + tol < r(indx(hi+1)) )
          break
        end
        hi = hi + 1;
      end
%
%  Points INDX(I+1) through INDX(HI) have an R value close to
%  point INDX(I).  Are they truly close to point INDEX(I)?
%
      for j = i + 1 : hi
        if ( unique(indx(j)) )
          dist = sqrt ( sum ( ( a(1:m,indx(i)) - a(1:m,indx(j)) ).^2 ) );
          if ( dist <= tol )
            unique(indx(j)) = 0;
          end
        end
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%
%  Copyright (c) 2014, Facebook, Inc.
%  All rights reserved.
%
%  This source code is licensed under the BSD-style license found in the
%  LICENSE file in the root directory of this source tree. An additional grant
%  of patent rights can be found in the PATENTS file in the same directory.
%
%%
function  [fail_idx, patches] = extract_person_patches(config, hits, dims, ...
    width_zoom, height_zoom,y_offset)
%Cropping people bounds

assert(config.DATASET == config.DATASET_ICCV);
bounds = int32(hits.bounds);
bounds(3:4,:) = bounds(3:4,:) + bounds(1:2,:);
bounds(1,:) = max(1,bounds(1,:));
bounds(2,:) = max(1,bounds(2,:));
patches = zeros([dims 3 hits.size],'uint8');
tic;
parfor i=1:hits.size
    img = load_image(hits.image_id(i),config);
    img = img(bounds(2,i):min(size(img,1),bounds(4,i)), ...
        bounds(1,i):min(size(img,2),bounds(3,i)),:);
    patches(:,:,:,i) = imresize(img, dims);
end
toc;
fail_idx = [];
end

author: facebookarchive

repo: pose-aligned-deep-networks


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

t = 0;              % current time
T = 1;              % maturity

F = 80:.25:120;       % possible forwards in T
A = 0.01:0.01:0.4;  % possible volatilities in T

f = 100;            % current spot asset
alpha = 0.2;        % current spot volatility

beta = 0.5;         % SABR CEV coeff
nu = 0.2;           % Vol of vol
rho = 0;         % Correlation

d_base = @(x,y) psabr(t,T,f,x,alpha,y,beta,nu,rho);   % approx density

xret = ((F-f)/f);   % possible relative changes

low = 0;            % limit integration low
up = 1;             % limit integration high

y_base = density1d(F,d_base,low,up);    % compute density
legend_base = 'Base';
title_plot = 'SABR Density';

y2_base = density2d(F,A,d_base);

%% Change alpha
d_low = @(x,y) psabr(t,T,f,x,alpha-0.1,y,beta,nu,rho); % low density
d_high = @(x,y) psabr(t,T,f,x,alpha+0.1,y,beta,nu,rho);% high density

y_low = density1d(F,d_low,low,up);
y_high = density1d(F,d_high,low,up);
legend_low = 'Changing \alpha low';
legend_high = 'Changing \alpha high';

createfigure_density(xret,y_base',y_low',y_high',title_plot,legend_base,legend_low,legend_high);

%% change beta
d_low = @(x,y) psabr(t,T,f,x,alpha,y,beta-.2,nu,rho); % low density
d_high = @(x,y) psabr(t,T,f,x,alpha,y,beta+.2,nu,rho);% high density

y_low = density1d(F,d_low,low,up);
y_high = density1d(F,d_high,low,up);
legend_low = 'Changing \beta low';
legend_high = 'Changing \beta high';

createfigure_density(xret,y_base',y_low',y_high',title_plot,legend_base,legend_low,legend_high);

y2_low = density2d(F,A,d_low);
y2_high = density2d(F,A,d_high);

create_density_2d(F',A,y2_low', y2_base', y2_high', legend_low, legend_base, legend_high, title_plot);

%% change nu
d_low = @(x,y) psabr(t,T,f,x,alpha,y,beta,nu-.1,rho); % low density
d_high = @(x,y) psabr(t,T,f,x,alpha,y,beta,nu+.1,rho);% high density

y_low = density1d(F,d_low,low,up);
y_high = density1d(F,d_high,low,up);
legend_low = 'Changing \nu low';
legend_high = 'Changing \nu high';

createfigure_density(xret,y_base',y_low',y_high',title_plot,legend_base,legend_low,legend_high);

%% change rho
d_low = @(x,y) psabr(t,T,f,x,alpha,y,beta,nu,rho-.6); % low density
d_high = @(x,y) psabr(t,T,f,x,alpha,y,beta,nu,rho+.6);% high density

y_low = density1d(F,d_low,low,up);
y_high = density1d(F,d_high,low,up);
legend_low = 'Changing \rho low';
legend_high = 'Changing \rho high';

createfigure_density(xret,y_base',y_low',y_high',title_plot,legend_base,legend_low,legend_high);

author: Sable

repo: mcbench-benchmarks


text:
function S = similaritymat(idx, samples, n)
%SIMILARITYMAT Similarity matrix.
% S = SIMILARITYMAT(IDX, SAMPLES, N) returns the N-by-N similarity matrix
% associated with the similarity of elements in the Ns-by-1 vector, IDX.
% The similarity value (0 or 1) is stored in the matrix S at the indices
% associated with the Ns-by-1 vector, SAMPLES
%
% Copyright (2009) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
% certain rights in this software.

S = zeros(n,n);

for i = 1:size(samples,1)
    S(samples(i),samples(i)) = 1;
    for j = (i+1):size(samples,1)
        if idx(i) == idx(j)
            S(samples(i),samples(j)) = 1;
            S(samples(j),samples(i)) = 1;
        end
    end
end


text:
function [X, optinf] = cbpdnms(D, S, lambda, opt)

% cbpdnms -- Convolutional Basis Pursuit DeNoising (Mask Simulation)
%
%         argmin_{x_k} (1/2)||W (\sum_k d_k * x_k - s)||_2^2 +
%                           lambda \sum_k ||x_k||_1
%
%         The solution is computed using an ADMM approach (see
%         boyd-2010-distributed) with efficient solution of the main
%         linear systems (see wohlberg-2016-efficient and
%         wohlberg-2016-boundary).
%
% Usage:
%       [X, optinf] = cbpdnms(D, S, lambda, opt)
%
% Input:
%       D           Dictionary filter set (3D array)
%       S           Input image
%       lambda      Regularization parameter
%       opt         Algorithm parameters structure
%
% Output:
%       X           Dictionary coefficient map set (3D array)
%       optinf      Details of optimisation
%
%
% Options structure fields:
%   Verbose          Flag determining whether iteration status is displayed.
%                    Fields are iteration number, functional value,
%                    data fidelity term, l1 regularisation term, and
%                    primal and dual residuals (see Sec. 3.3 of
%                    boyd-2010-distributed). The value of rho is also
%                    displayed if options request that it is automatically
%                    adjusted.
%   MaxMainIter      Maximum main iterations
%   AbsStopTol       Absolute convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   RelStopTol       Relative convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   L1Weight         Weighting array for coefficients in l1 norm of X
%   Y0               Initial value for Y
%   U0               Initial value for U
%   rho              ADMM penalty parameter
%   AutoRho          Flag determining whether rho is automatically updated
%                    (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoRhoPeriod    Iteration period on which rho is updated
%   RhoRsdlRatio     Primal/dual residual ratio in rho update test
%   RhoScaling       Multiplier applied to rho when updated
%   AutoRhoScaling   Flag determining whether RhoScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, RhoScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier.
%   RhoRsdlTarget    Residual ratio targeted by auto rho update policy.
%   StdResiduals     Flag determining whether standard residual definitions
%                    (see Sec 3.3 of boyd-2010-distributed) are used instead
%                    of normalised residuals (see wohlberg-2015-adaptive)
%   RelaxParam       Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed)
%   NoBndryCross     Flag indicating whether all solution coefficients
%                    corresponding to filters crossing the image boundary
%                    should be forced to zero.
%   AuxVarObj        Flag determining whether objective function is computed
%                    using the auxiliary (split) variable
%   HighMemSolve     Use more memory for a slightly faster solution
%   W                Synthesis spatial weighting matrix
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2016-05-10
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 4,
  opt = [];
end
checkopt(opt, defaultopts([]));
opt = defaultopts(opt);

% Set up status display for verbose operation
hstr = 'Itn   Fnc       DFid      l1        r         s      ';
sfms = '%4d %9.2e %9.2e %9.2e %9.2e %9.2e';
nsep = 54;
if opt.AutoRho,
  hstr = [hstr '   rho  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.Verbose && opt.MaxMainIter > 0,
  disp(hstr);
  disp(char('-' * ones(1,nsep)));
end

% Start timer
tstart = tic;

% Insert impulse filter into dictionary
imp = zeros(size(D,1), size(D,2), 1);
imp(1,1,1) = 1.0;
Di = cat(3, D, imp);

% Collapsing of trailing singleton dimensions greatly complicates
% handling of both SMV and MMV cases. The simplest approach would be
% if s could always be reshaped to 4d, with dimensions consisting of
% image rows, image cols, a single dimensional placeholder for number
% of filters, and number of measurements, but in the single
% measurement case the third dimension is collapsed so that the array
% is only 3d.
if size(S,3) > 1 && size(S,4) == 1,
  xsz = [size(S,1) size(S,2) size(Di,3) size(S,3)];
  % Insert singleton 3rd dimension (for number of filters) so that
  % 4th dimension is number of images in input s volume
  S = reshape(S, [size(S,1) size(S,2) 1 size(S,3)]);
  if ~isscalar(opt.W) & ndims(opt.W) > 2,
    opt.W = reshape(opt.W, [size(opt.W,1) size(opt.W,2) 1 size(opt.W,3)]);
  end
else
  xsz = [size(S,1) size(S,2) size(Di,3) size(S,4)];
end
IYW = 1.0 - opt.W;

% Compute filters in DFT domain
Df = fft2(Di, size(S,1), size(S,2));
% Convolve-sum and its Hermitian transpose
Dop = @(x) sum(bsxfun(@times, Df, x), 3);
DHop = @(x) bsxfun(@times, conj(Df), x);
% Compute signal in DFT domain
Sf = fft2(S);
% S convolved with all filters in DFT domain
DSf = DHop(Sf);

% Default lambda is 1/10 times the lambda value beyond which the
% solution is a zero vector
if nargin < 3 | isempty(lambda),
  b = ifft2(DHop(Sf), 'symmetric');
  lambda = 0.1*max(vec(abs(b)));
end

% Set up algorithm parameters and initialise variables
rho = opt.rho;
if isempty(rho), rho = 50*lambda+1; end;
if isempty(opt.RhoRsdlTarget),
  if opt.StdResiduals,
    opt.RhoRsdlTarget = 1;
  else
    opt.RhoRsdlTarget = 1 + (18.3).^(log10(lambda) + 1);
  end
end
if opt.HighMemSolve,
  C = bsxfun(@rdivide, Df, sum(Df.*conj(Df), 3) + rho);
else
  C = [];
end
Nx = prod(xsz);
optinf = struct('itstat', [], 'opt', opt);
r = Inf;
s = Inf;
epri = 0;
edua = 0;

% Initialise main working variables
X = [];
if isempty(opt.Y0),
  Y = zeros(xsz);
else
  Y = opt.Y0;
end
Yprv = Y;
if isempty(opt.U0),
  if isempty(opt.Y0),
    U = zeros(xsz);
  else
    U = (lambda/rho)*sign(Y);
  end
else
  U = opt.U0;
end

% Main loop
k = 1;
while k <= opt.MaxMainIter && (r > epri | s > edua),

  % Solve X subproblem
  Xf = solvedbi_sm(Df, rho, DSf + rho*fft2(Y - U), C);
  X = ifft2(Xf, 'symmetric');

  % See pg. 21 of boyd-2010-distributed
  if opt.RelaxParam == 1,
    Xr = X;
  else
    Xr = opt.RelaxParam*X + (1-opt.RelaxParam)*Y;
  end

  % Solve Y subproblem
  Y(:,:,1:(end-1),:) = shrink(Xr(:,:,1:(end-1),:) + U(:,:,1:(end-1),:), ...
                              (lambda/rho)*opt.L1Weight);
  Y(:,:,end,:) = bsxfun(@times, IYW, Xr(:,:,end,:) + U(:,:,end,:));
  if opt.NoBndryCross,
    Y((end-size(D,1)+2):end,:,1:(end-1),:) = 0;
    Y(:,(end-size(D,2)+2):end,1:(end-1),:) = 0;
  end

  % Update dual variable
  U = U + Xr - Y;

  % Compute data fidelity term in Fourier domain (note normalisation)
  if opt.AuxVarObj,
    Yf = fft2(Y); % This represents unnecessary computational cost
    Jdf = sum(vec(abs(sum(bsxfun(@times,Df,Yf),3)-Sf).^2))/(2*xsz(1)*xsz(2));
    Jl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, Y(:,:,1:(end-1),:)))));
  else
    Jdf = sum(vec(abs(sum(bsxfun(@times,Df,Xf),3)-Sf).^2))/(2*xsz(1)*xsz(2));
    Jl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, X(:,:,1:(end-1),:)))));
  end
  Jfn = Jdf + lambda*Jl1;

  nX = norm(X(:)); nY = norm(Y(:)); nU = norm(U(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    r = norm(vec(X - Y));
    s = norm(vec(rho*(Yprv - Y)));
    epri = sqrt(Nx)*opt.AbsStopTol+max(nX,nY)*opt.RelStopTol;
    edua = sqrt(Nx)*opt.AbsStopTol+rho*nU*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    r = norm(vec(X - Y))/max(nX,nY);
    s = norm(vec(Yprv - Y))/nU;
    epri = sqrt(Nx)*opt.AbsStopTol/max(nX,nY)+opt.RelStopTol;
    edua = sqrt(Nx)*opt.AbsStopTol/(rho*nU)+opt.RelStopTol;
  end

  % Record and display iteration details
  tk = toc(tstart);
  optinf.itstat = [optinf.itstat; [k Jfn Jdf Jl1 r s epri edua rho tk]];
  if opt.Verbose,
    if opt.AutoRho,
      disp(sprintf(sfms, k, Jfn, Jdf, Jl1, r, s, rho));
    else
      disp(sprintf(sfms, k, Jfn, Jdf, Jl1, r, s));
    end
  end

  % See wohlberg-2015-adaptive and pp. 20-21 of boyd-2010-distributed
  if opt.AutoRho,
    if k ~= 1 && mod(k, opt.AutoRhoPeriod) == 0,
      if opt.AutoRhoScaling,
        rhomlt = sqrt(r/(s*opt.RhoRsdlTarget));
        if rhomlt < 1, rhomlt = 1/rhomlt; end
        if rhomlt > opt.RhoScaling, rhomlt = opt.RhoScaling; end
      else
        rhomlt = opt.RhoScaling;
      end
      rsf = 1;
      if r > opt.RhoRsdlTarget*opt.RhoRsdlRatio*s, rsf = rhomlt; end
      if s > (opt.RhoRsdlRatio/opt.RhoRsdlTarget)*r, rsf = 1/rhomlt; end
      rho = rsf*rho;
      U = U/rsf;
      if opt.HighMemSolve && rsf ~= 1,
        C = bsxfun(@rdivide, Df, sum(Df.*conj(Df), 3) + rho);
      end
    end
  end

  Yprv = Y;
  k = k + 1;

end

% Record run time and working variables
optinf.runtime = toc(tstart);
optinf.X = X;
optinf.Xf = Xf;
optinf.Y = Y;
optinf.U = U;
optinf.lambda = lambda;
optinf.rho = rho;

if opt.Verbose && opt.MaxMainIter > 0,
  disp(char('-' * ones(1,nsep)));
end

% Remove coefficient map for impulse filter
X = X(:,:,1:(end-1), :);

return

function u = vec(v)

  u = v(:);

return

function u = shrink(v, lambda)

  if isscalar(lambda),
    u = sign(v).*max(0, abs(v) - lambda);
  else
    u = sign(v).*max(0, bsxfun(@minus, abs(v), lambda));
  end

return

function opt = defaultopts(opt)

  if ~isfield(opt,'Verbose'),
    opt.Verbose = 0;
  end
  if ~isfield(opt,'MaxMainIter'),
    opt.MaxMainIter = 1000;
  end
  if ~isfield(opt,'AbsStopTol'),
    opt.AbsStopTol = 0;
  end
  if ~isfield(opt,'RelStopTol'),
    opt.RelStopTol = 1e-4;
  end
  if ~isfield(opt,'L1Weight'),
   opt.L1Weight = 1;
  end
  if ~isfield(opt,'Y0'),
    opt.Y0 = [];
  end
  if ~isfield(opt,'U0'),
    opt.U0 = [];
  end
  if ~isfield(opt,'rho'),
    opt.rho = [];
  end
  if ~isfield(opt,'AutoRho'),
    opt.AutoRho = 1;
  end
  if ~isfield(opt,'AutoRhoPeriod'),
    opt.AutoRhoPeriod = 1;
  end
  if ~isfield(opt,'RhoRsdlRatio'),
    opt.RhoRsdlRatio = 1.2;
  end
  if ~isfield(opt,'RhoScaling'),
    opt.RhoScaling = 100;
  end
  if ~isfield(opt,'AutoRhoScaling'),
    opt.AutoRhoScaling = 1;
  end
  if ~isfield(opt,'RhoRsdlTarget'),
    opt.RhoRsdlTarget = [];
  end
  if ~isfield(opt,'StdResiduals'),
    opt.StdResiduals = 1;
  end
  if ~isfield(opt,'RelaxParam'),
    opt.RelaxParam = 1.8;
  end
  if ~isfield(opt,'NoBndryCross'),
    opt.NoBndryCross = 0;
  end
  if ~isfield(opt,'AuxVarObj'),
    opt.AuxVarObj = 0;
  end
  if ~isfield(opt,'HighMemSolve'),
    opt.HighMemSolve = 0;
  end
  if ~isfield(opt,'W'),
    opt.W = 1.0;
  end

return

author: thfylsty


text:
function[blobs] = rpUijlings2013(image, varargin)
% [blobs] = rpUijlings2013(image, varargin)
%
% Create a list of region proposals (blobs) of an image using the method of Uijlings2013.
%
% For information on the returned blobs see blob().
%
% Copyright by Holger Caesar, 2014

% Parse input
p = inputParser;
% Default settings suggested by Jasper
addParameter(p, 'ks', [50, 100]);
addParameter(p, 'sigma', 0.8);
addParameter(p, 'colorTypes', {'Hsv', 'Lab'});
addParameter(p, 'simFunctionHandles', {@SSSimColourTextureSizeFillOrig, @SSSimTextureSizeFill});
parse(p, varargin{:});

ks = p.Results.ks;
sigma = p.Results.sigma;
colorTypes = p.Results.colorTypes;
simFunctionHandles = p.Results.simFunctionHandles;

% Init SS
hierarchyCount = numel(ks) * numel(colorTypes) * numel(simFunctionHandles);
blobs = cell(hierarchyCount, 1);
priorities = cell(hierarchyCount, 1);
hBlobsInd = 1;

% Perform Selective Search on all parameter combinations
for ksInd = 1 : numel(ks),
    for colorTypeInd = 1 : numel(colorTypes),
        for simHndsI = 1 : numel(simFunctionHandles),
            % Compute hierarchical grouping with current settings
            [~, blobIndIm, blobBoxesInIm, hierarchy, priority] = Image2HierarchicalGrouping(image, sigma, ks(ksInd), ks(ksInd), colorTypes{colorTypeInd}, simFunctionHandles(simHndsI));

            % Convert boxes to blobs (when there is only one blob, we can't
            % create a hierarchy, i.e. for k=400)
            if isempty(hierarchy),
                blobs{hBlobsInd} = RecreateBlobHierarchyIndIm(blobIndIm, blobBoxesInIm, []);
            else
                blobs{hBlobsInd} = RecreateBlobHierarchyIndIm(blobIndIm, blobBoxesInIm, hierarchy{1});
            end;

            % Store priority
            priorities{hBlobsInd} = priority;
            hBlobsInd = hBlobsInd + 1;
        end;
    end;
end;

% Put all parameter results in a single cell array
blobs = flattenCellArray(blobs);
priorities = flattenCellArray(priorities);

% Do pseudo random sorting as in paper
priorities = priorities .* rand(size(priorities));
[~, sortIds] = sort(priorities, 'ascend');
blobs = blobs(sortIds, :);

% Convert to col struct
blobs = [blobs{:}];
blobs = blobs(:);

author: nightrome

repo: matconvnet-calvin


text:
function y = norm_cdf(x,mu,sigma)
%NORM_CDF Normal cumulative probability density function (cdf)
%
%   Y = NORMCDF(X,MU,SIGMA) Returns the normal cdf with
%   mean, MU, and standard deviation, SIGMA, at the values in X.
%
%   The size of X is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.
%
% Copyright (c) 1998-2004,2011 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3,
  sigma = 1;
end

if nargin < 2;
  mu = 0;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y = 0.5 * erfc(-((x-mu) ./ sigma) ./ sqrt(2));

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
%% Lucas-Kanade Sparse Optical Flow
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/tapi/pyrlk_optical_flow.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/gpu/pyrlk_optical_flow.cpp>
%

function pyrlk_optical_flow_demo()
    % Prepare video source
    if true
        vid = 0;
    elseif true
        vid = fullfile(mexopencv.root(), 'test', '768x576.avi');
    elseif mexopencv.require('vision')
        vid = fullfile(toolboxdir('vision'), 'visiondata', 'visiontraffic.avi');
    end
    cap = cv.VideoCapture(vid);
    assert(cap.isOpened(), 'Failed to initialize capturing');

    % Grab first frame
    frame = cap.read();
    assert(~isempty(frame), 'Failed to read frame');
    gray0 = cv.cvtColor(frame, 'RGB2GRAY');

    % Plot
    hImg = imshow(frame);
    title('PyrLK [Sparse]')

    % Main loop
    while ishghandle(hImg)
        % Grab next frame
        frame = cap.read();
        if isempty(frame), break; end
        gray1 = cv.cvtColor(frame, 'RGB2GRAY');

        % Detect corners in previous frame
        pts0 = cv.goodFeaturesToTrack(gray0, ...
            'MaxCorners',100, 'QualityLevel',0.01, 'MinDistance',0);
        pts0 = cat(1, pts0{:});
        if isempty(pts0), continue; end

        % Compute sparse optical flow (track points from previous to current frame)
        [pts1, status] = cv.calcOpticalFlowPyrLK(gray0, gray1, pts0);
        pts1 = cat(1, pts1{:});
        status = logical(status);

        % Draw sparse flow (only good points for which flow was found)
        if any(status)
            frame = drawArrows(frame, pts0(status,:), pts1(status,:));
        end

        % Display result
        set(hImg, 'CData',frame);
        drawnow;

        % Next iteration
        gray0 = gray1;
        pts0 = pts1;
    end
    cap.release();
end

%% Helper function

function img = drawArrows(img, pts0, pts1)
    %DRAWARROWS  Draw sparse optical flow
    %
    % See also: quiver
    %

    % compute angle and hypotenuse
    uv = pts0 - pts1;
    ang = atan2(uv(:,2), uv(:,1));
    mag = hypot(uv(:,2), uv(:,1));

    % skip short ones
    idx = (mag < 1.0);
    if all(idx), return; end
    pts0(idx,:) = [];
    pts1(idx,:) = [];
    ang(idx,:) = [];
    mag(idx,:) = [];

    % Options for line drawing
    props = {'Color',[255 0 0], 'Thickness',1};

    % Here we lengthen the arrow by a factor of three,
    % and draw the main line of the arrow
    pts1 = pts0 - 3 * bsxfun(@times, mag, [cos(ang) sin(ang)]);
    img = cv.line(img, pts0, pts1, props{:});

    % Now we draw the tips of the arrow. We do some scaling so that the
    % tips look proportional to the main line of the arrow
    pts0 = pts1 + 9 * [cos(ang + pi/4) sin(ang + pi/4)];
    img = cv.line(img, pts0, pts1, props{:});
    pts0 = pts1 + 9 * [cos(ang - pi/4) sin(ang - pi/4)];
    img = cv.line(img, pts0, pts1, props{:});
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% std_selsubject() - Helper function for std_erpplot(), std_specplot()
%                    and std_erspplot() to select specific subject when
%                    plotting channel data.
% Usage:
%  >> data = std_selsubject( data, subject, setinds, allsubjects);
%
% Inputs:
%  data  -  [cell array] mean data for each subject group and/or data
%           condition. For example, to compute mean ERPs statistics from a
%           STUDY for epochs of 800 frames in two conditions from three
%           groups of 12 subjects,
%           >> data = { [800x12] [800x12] [800x12];... % 3 groups, cond 1
%                       [800x12] [800x12] [800x12] };  % 3 groups, cond 2
% subject - [string] subject name
% setinds - [cell array] set indices for each of the last dimension of the
%           data cell array.
%           >> setinds = { [12] [12] [12];... % 3 groups, cond 1
%                          [12] [12] [12] };  % 3 groups, cond 2
% allsubject - [cell array] all subjects (same order as in
%              STUDY.datasetinfo)
%
% Output:
%  data       - [cell array] data array with the subject or component selected
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: std_erpplot(), std_specplot() and std_erspplot()

function [data] = std_selsubject(data, subject, setinds, allsubjects, optndims);

if nargin < 2
    help std_selsubject;
    return;
end;

optndims = max(optndims, ndims(data{1}));
if isempty(strmatch(lower(subject), lower(allsubjects)))
    error(sprintf('Cannot select subject %s in list %s', subject, vararg2str({ allsubjects })));
end;

% plot specific subject
% ---------------------
if size(setinds{1},1) > 1 && size(setinds{1},2) > 1 % single trials
    % possible subject indices
    selectInds = strmatch(lower(subject), lower(allsubjects));
    for c = 1:size(data,1)
        for g = 1:size(data,2)
            selectCol = [];
            for ind = 1:length(selectInds)
                selectCol = [ selectCol find(setinds{c,g}  == selectInds') ];
            end;
            if optndims == 2
                data{c,g} = data{c,g}(:,selectCol); %2-D
            elseif optndims == 3
                data{c,g} = data{c,g}(:,:,selectCol); %3-D
            else
                data{c,g} = data{c,g}(:,:,:,selectCol); %4-D
            end;
        end;
    end;
else
    for c = 1:size(data,1)
        for g = 1:size(data,2)
            subjectind = strmatch(lower(subject), lower(allsubjects));
            l = zeros(size(setinds{c,g}));
            for iSubj = 1:length(subjectind), l = l | setinds{c,g} == subjectind(iSubj); end;
            if optndims == 2
                data{c,g}(:,~l) = []; %2-D
            elseif optndims == 3
                data{c,g}(:,:,~l) = []; %3-D
            else
                data{c,g}(:,:,:,~l) = []; %4-D
            end;
        end;
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function visualizeData( X, k, IDX, types, C )
% Project high dim. data unto principal components (PCA) for visualization.
%
% Optionally IDX can be specified to indicate different classes for the
% points; in this case points in different classes are displayed using
% different colors. Up to 12 types are handled (for technical reasons
% involving plot), any cluster with a label>12 is assigned the label 12.
%
% USAGE
%  visualizeData( X, k, [IDX], [types], [C] )
%
% INPUTS
%  X       - column vector of data - N vectors of dimension p (X is Nxp)
%  k       - dimension to which to reduce data (2 or 3)
%  IDX     - [] cluster membership [see kmeans2.m]
%  types   - [] cell array of length ntypes of text labels for each type
%  C       - [] cluster centers (Kxp)
%
% OUTPUTS
%
% EXAMPLE
%  X = [randn(100,5); randn(100,5)+4];
%  C = [mean(X(1:100,:)); mean(X(101:200,:))];
%  IDX = [ones(100,1); 2*ones(100,1)];
%  visualizeData( X, 2, IDX, {'type1','type2' }, C);
%
% See also KMEANS2, DEMOCLUSTER
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(IDX) ); IDX=[]; end
if( nargin<4 || isempty(types) ); types=[]; end
if( nargin<5 || isempty(C) ); C=[]; end

% apply PCA if necessary
if( size(X,2)~= k )
  [ U, mu ] = pca( X' );
  X = pcaApply( X', U, mu, k )';
  if(~isempty(C)); C = pcaApply( C', U, mu, k )'; end
end

%%% get k
k = size(X,2);
if( k==1 ); X = [X zeros(size(X))]; k = 2; end
if( k>3 ); error( 'k must be <= 3'); end

%%% show points
if( isempty(IDX) )
  if( k==2 )
    plot( X(:,1), X(:,2), '.' );
  elseif( k==3 )
    plot3( X(:,1), X(:,2), X(:,3), '.' );
  end;

else
  IDX(IDX>12)=12;  m=max(IDX);

  if( k==2)
    % plot points
    R = cell(1,3*m+3);
    for i=1:m;
      R((3*i-2):(3*i)) = {X(IDX==i,1), X(IDX==i,2), '.'};
    end;
    R((3*m+1):(3*m+3)) = {X(IDX==-1,1), X(IDX==-1,2), 'k.'};
    plot( R{:} );

    % plot centers
    if( ~isempty(C) )
      R=cell(1,3*m);
      for i=1:m;  R((3*i-2):(3*i)) = {C(i,1), C(i,2), 'x'}; end
      hold('on');  plot( R{:}, 'MarkerSize', 30 );  hold('off');
    end;

  elseif( k==3 )
    % plot points
    R = cell(1,4*m+4);
    for i=1:m;
      R((4*i-3):(4*i)) = {X(IDX==i,1), X(IDX==i,2), X(IDX==i,3), '.'};
    end;
    R((4*m+1):(4*m+4)) = {X(IDX==-1,1), X(IDX==-1,2), X(IDX==-1,3), 'k.'};
    plot3( R{:} );

    % plot centers
    if( ~isempty(C) )
      R=cell(1,4*m);
      for i=1:m;  R((4*i-3):(4*i)) = {C(i,1), C(i,2), C(i,3), 'x'}; end
      hold('on'); plot3( R{:}, 'MarkerSize', 30 );  hold('off');
    end
  end
end
axis('equal');

%%% show legend if types is provided
if(~isempty(types));  legend(types); end

author: zouchuhang

repo: LayoutNet


text:
function v = subsref( x, S )
try
    v = subsref( x.value_, S );
catch
    cvx_throw( 'Invalid tuple reference: %s%s', inputname(1), cvx_subs2str( S ) );
end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/mmv_phase_transition_snr_40db_s_2.mat';
N = 256;
S = 2;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 40;
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function edge = lineToEdge(line)
%LINETOEDGE Convert a straight line to a finite edge.
%
%   EDGE = lineToEdge(LINE)
%   Returns the edge with same origin as the line LINE, and with second
%   extremity corresponding to the addition of line origin and direction.
%   LINE is represented as [X0 Y0  DX DY]
%   EDGE is represented as [X1 Y1  X2 Y2]
%
%   Example
%     line = [3 4  1 2];
%     edge = lineToEdge(line)
%     edge =
%          3   4   4   6
%
%   See also
%     lines2d, edges2d, edgeToLine

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-05-07, using Matlab 9.6.0.1072779 (R2019a)
% Copyright 2019-2022 INRA - Cepia Software Platform

edge = [line(:, 1:2) line(:,1:2)+line(:,3:4)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function a = daub8 ( n )

%*****************************************************************************80
%
%% DAUB8 returns the DAUB8 matrix.
%
%  Discussion:
%
%    The DAUB8 matrix is the Daubechies wavelet transformation matrix
%    with 8 coefficients.
%
%  Properties:
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gilbert Strang, Truong Nguyen,
%    Wavelets and Filter Banks,
%    Wellesley-Cambridge Press, 1997,
%    ISBN: 0-9614088-7-1,
%    LC: TK7872.F5S79 / QA403.3.S87
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 8, and a multiple of 2.
%
%    Output, real A(N,N), the matrix.
%
  if ( n < 8 || mod ( n, 2 ) ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DAUB8 - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 6 and a multiple of 2.\n' );
    error ( 'DAUB8 - Fatal error!' );
  end

  a = zeros ( n, n );

  c = [ ...
    0.2303778133088964, ...
    0.7148465705529154, ...
    0.6308807679298587, ...
   -0.0279837694168599, ...
   -0.1870348117190931, ...
    0.0308413818355607, ...
    0.0328830116668852, ...
   -0.0105974017850690 ]';

  for i = 1 : 2 : n - 1

    a(i,i)                  =   c(1);
    a(i,i+1)                =   c(2);
    a(i,i4_wrap(i+2,1,n))   =   c(3);
    a(i,i4_wrap(i+3,1,n))   =   c(4);
    a(i,i4_wrap(i+4,1,n))   =   c(5);
    a(i,i4_wrap(i+5,1,n))   =   c(6);
    a(i,i4_wrap(i+6,1,n))   =   c(7);
    a(i,i4_wrap(i+7,1,n))   =   c(8);

    a(i+1,i)                =   c(8);
    a(i+1,i+1)              = - c(7);
    a(i+1,i4_wrap(i+2,1,n)) =   c(6);
    a(i+1,i4_wrap(i+3,1,n)) = - c(5);
    a(i+1,i4_wrap(i+4,1,n)) =   c(4);
    a(i+1,i4_wrap(i+5,1,n)) = - c(3);
    a(i+1,i4_wrap(i+6,1,n)) =   c(2);
    a(i+1,i4_wrap(i+7,1,n)) = - c(1);

  end

  return
end


text:
subjectDir = '/biac2/wandell2/data/reading_longitude/templates/child_new/SIRL54warp3';
fileFragment = '*sn*';
[s,subCode] = findSubjects(subjectDir,fileFragment,{'tk040817'});
N = length(s);
tdir = '/silver/scr1/data/templates/child_new';
tname = 'SIRL54';
avgdir = fullfile(tdir,[tname 'warp3']);
template = load(fullfile(avgdir,'average_dt6'));
xformDtToAcpc = template.xformToAcPc;
dtMmPerVox = template.mmPerVox;
avgBrain = template.anat.img;
avgBrain(template.anat.brainMask<0.25) = 0;
avgBrain(template.anat.brainMask<0.5) = avgBrain(template.anat.brainMask<0.5)*.5;
outDir = '/silver/scr1/readingGroupStats/';
if(~exist(outDir,'dir')) mkdir(outDir); end

[behData,colNames] = dtiGetBehavioralData(subCode);
readerType = strmatch('Type of Reader',colNames);
pr = find(behData(:,readerType)==-1);
gr = find(behData(:,readerType)==1);
npr = length(pr);
ngr = length(gr);

dt = load(s{pr(1)});
prDt6 = zeros([size(dt.dt6) npr]);
prDt6(:,:,:,:,1) = dt.dt6;
mnB0 =  double(dt.b0);
for(ii=2:npr)
  disp(['Loading ' s{pr(ii)} '...']);
  dt = load(s{pr(ii)});
  dt.dt6(isnan(dt.dt6)) = 0;
  prDt6(:,:,:,:,ii) = dt.dt6;
  mnB0 = mnB0 + double(dt.b0);
end
grDt6 = zeros([size(dt.dt6) ngr]);
for(ii=1:ngr)
  disp(['Loading ' s{gr(ii)} '...']);
  dt = load(s{gr(ii)});
  dt.dt6(isnan(dt.dt6)) = 0;
  grDt6(:,:,:,:,ii) = dt.dt6;
  mnB0 = mnB0 + double(dt.b0);
end
mnB0 = mnB0./(npr+ngr);
mask = mnB0>250 & all(squeeze(prDt6(:,:,:,1,:)),4)>0 & all(squeeze(grDt6(:,:,:,1,:)),4)>0;
prDt6_ind = dtiImgToInd(prDt6, mask);
clear prDt6;
grDt6_ind = dtiImgToInd(grDt6, mask);
clear grDt6;
% convert to log tensors
[prVec, prVal] = dtiEig(prDt6_ind);
prVal(prVal<0) = 0;
prLogDt6_ind = dtiEigComp(prVec, log(prVal));
[grVec, grVal] = dtiEig(grDt6_ind);
grVal(grVal<0) = 0;
grLogDt6_ind = dtiEigComp(grVec, log(grVal));

showSlices = [20:60];

% Watson test for principal eigenvector difference
[prDir.mean, prDir.stdev, prDir.n, prDir.sbar] = dtiDirMean(squeeze(prVec(:,:,1,:)));
[grDir.mean, grDir.stdev, grDir.n, grDir.sbar] = dtiDirMean(squeeze(grVec(:,:,1,:)));
[T, DISTR, df] = dtiDirTest(prDir.sbar, prDir.n, grDir.sbar, grDir.n);
Timg = dtiIndToImg(T, mask);
fThresh = finv(1-10^-4, df(1), df(2));
fMax = finv(1-10^-12, df(1), df(2));
Timg(Timg>fMax) = fMax;
fMax = max(Timg(:));
figure; imagesc(makeMontage(Timg,showSlices)); axis image; colormap hot; colorbar;
set(gcf,'Name','Dir test'); title(sprintf('fthresh (p<10^-^4) = %0.1f',fThresh));
dtiWriteNiftiWrapper(Timg, xformDtToAcpc, fullfile(outDir,['dir_' DISTR '-test_' num2str(df(1)) ',' num2str(df(2)) 'df.nii.gz']));

% Simple FDR analysis for Watson test
%
fdrVal = 0.10; fdrType = 'original';
T(isnan(T)) = 0;
pvals = 1-fcdf(T, df(1), df(2));
[n_signif,index_signif] = fdr(pvals,fdrVal,fdrType,'mean');
pThreshFDR = max(pvals(index_signif));
% Convert back to an fThreshold
fThreshFDR = finv(1-pThreshFDR, df(1), df(2));
disp(sprintf('f-threshold for FDR (%s method) of %0.3f: %0.2f (%0.3f).\n',...
             fdrType,fdrVal,fThreshFDR,fThreshFDR/fMax));
logPimg = dtiIndToImg(-log10(pvals), mask);
cmap = autumn(256);
maxLogP = 10;
%minLogP = -log10(pThreshFDR);
minLogP = -log10(0.01);
sl = [-20:2:50];

mrAnatOverlayMontage(logPimg, xformDtToAcpc, avgBrain, template.anat.mmPerVox, cmap, [minLogP maxLogP], sl);

% Log-normal tensor tests
%
% Test for eigenvector differences
[prLogDt.mean, prLogDt.stdev, prLogDt.n] = dtiLogTensorMean(prLogDt6_ind);
[grLogDt.mean, grLogDt.stdev, grLogDt.n] = dtiLogTensorMean(grLogDt6_ind);

[T, DISTR, df] = dtiLogTensorTest('vec', prLogDt.mean, prLogDt.stdev, prLogDt.n, grLogDt.mean, grLogDt.stdev, grLogDt.n);

Timg = dtiIndToImg(T, mask);
fThresh = finv(1-10^-4, df(1), df(2));
fMax = finv(1-10^-12, df(1), df(2));
Timg(Timg>fMax) = fMax;
dtiWriteNiftiWrapper(Timg, xformDtToAcpc, fullfile(outDir,['vec_' DISTR '-test_' num2str(df(1)) ',' num2str(df(2)) 'df.nii.gz']));

% Sqrt(F) is the standardized distance between the groups.
%figure; imagesc(makeMontage(sqrt(Timg),[20:55])); axis image; colormap hot; colorbar;
%set(gcf,'Name','Vec Standardized Distance');

% Simple FDR analysis for eigenvector differences
%
fdrVal = 0.05; fdrType = 'original';
T(isnan(T)) = 0;
pvals = 1-fcdf(T, df(1), df(2));
[n_signif,index_signif] = fdr(pvals,fdrVal,fdrType,'mean');
%pThreshFDR = max(pvals(index_signif));
% Convert back to an fThreshold
%fThreshFDR = finv(1-pThreshFDR, df(1), df(2));
%disp(sprintf('f-threshold for FDR (%s method) of %0.3f: %0.2f (%0.3f).\n',...
%             fdrType,fdrVal,fThreshFDR,fThreshFDR/fMax));

logPimg = dtiIndToImg(-log10(pvals), mask);
cmap = autumn(256);
maxLogP = 10;
%minLogP = -log10(pThreshFDR);
minLogP = -log10(0.01);
mrAnatOverlayMontage(logPimg, xformDtToAcpc, avgBrain, template.anat.mmPerVox, cmap, [minLogP maxLogP], sl);

% Test for eigenvalue differences
%
[T, DISTR, df] = dtiLogTensorTest('val', prLogDt.mean, prLogDt.stdev, prLogDt.n, grLogDt.mean, grLogDt.stdev, grLogDt.n);

T(isnan(T)) = 0;
pvals = 1-fcdf(T, df(1), df(2));
[n_signif,index_signif] = fdr(pvals,fdrVal,fdrType,'mean');
pThreshFDR = max(pvals(index_signif));
% Convert back to an fThreshold
fThreshFDR = finv(1-pThreshFDR, df(1), df(2));
disp(sprintf('f-threshold for FDR (%s method) of %0.3f: %0.2f (%0.3f).\n',...
             fdrType,fdrVal,fThreshFDR,fThreshFDR/fMax));
Timg = dtiIndToImg(T, mask);
fMax = finv(1-10^-12, df(1), df(2));
Timg(Timg>fMax) = fMax;
dtiWriteNiftiWrapper(Timg, xformDtToAcpc, fullfile(outDir,['val_' DISTR '-test_' num2str(df(1)) ',' num2str(df(2)) 'df.nii.gz']));
logPimg = dtiIndToImg(-log10(pvals), mask);
cmap = autumn(256);
maxLogP = 10;
minLogP = -log10(pThreshFDR);
%minLogP = -log10(0.001);
mrAnatOverlayMontage(logPimg, xformDtToAcpc, avgBrain, template.anat.mmPerVox, cmap, [minLogP maxLogP], sl);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


%==============================================================================
% This code is part of the Finite Element Method app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/FAIRFEM
%==============================================================================
%
% function [Dx1,Dx2,Dx3] = getGradientMatrixFEM(Mesh)
%
% builds gradient operator for tetrahedral finite element discretization
% with piecewise linear basis functions
%
% Input:
%
% Mesh          - description of mesh, struct
%
% Output:
%
% [Dx1,Dx2,Dx3] - discrete partial derivative operators
%
%==============================================================================

function [Dx1,Dx2,Dx3] = getGradientMatrixFEM(Mesh,matrixFree)
if nargin==0, help(mfilename); runMinimalExample; return; end;

if not(exist('matrixFree','var'))||isempty(matrixFree),
    matrixFree=0;
end

Dx3 = [];

xn  = Mesh.xn;
vol = Mesh.vol;
dim = Mesh.dim;

flag = [num2str(dim) 'D'];
if matrixFree
    flag = [flag '-mf'];
else
    flag = [flag '-mb'];
end

switch flag
    case '2D-mb'
        % compute edges
        e1 =  Mesh.P1*xn - Mesh.P3*xn;
        e2 =  Mesh.P2*xn - Mesh.P3*xn;
        % compute gradients of basis functions
        dphi1 =   [ e2(:,2) -e2(:,1)]./[2*vol,2*vol];
        dphi2 =   [-e1(:,2)  e1(:,1)]./[2*vol,2*vol];
        dphi3 = -dphi1-dphi2;

        Dx1 = sdiag(dphi1(:,1))*Mesh.P1 +  sdiag(dphi2(:,1))*Mesh.P2  + sdiag(dphi3(:,1))*Mesh.P3;
        Dx2 = sdiag(dphi1(:,2))*Mesh.P1 +  sdiag(dphi2(:,2))*Mesh.P2  + sdiag(dphi3(:,2))*Mesh.P3;

        if nargout==1, % return gradient operator for vector field
            B = [Dx1;Dx2];
            Dx1 = blkdiag(B,B);
        end

    case '2D-mf'
        % compute edges
        x1 =  Mesh.mfPi(xn,1);
        x2 =  Mesh.mfPi(xn,2);
        x3 =  Mesh.mfPi(xn,3);
        e1 =  x1 - x3;
        e2 =  x2 - x3;

        % compute gradients of basis functions
        dphi1 =   [ e2(:,2) -e2(:,1)]./[2*vol,2*vol];
        dphi2 =   [-e1(:,2)  e1(:,1)]./[2*vol,2*vol];
        dphi3 = -dphi1 - dphi2;

        Dxi    = @(x,i)   dphi1(:,i).*Mesh.mfPi(x,1) ...
                        + dphi2(:,i).*Mesh.mfPi(x,2) ...
                        + dphi3(:,i).*Mesh.mfPi(x,3);

        Dxiadj = @(x,i)   Mesh.mfPi(dphi1(:,i).*x,1) ....
                        + Mesh.mfPi(dphi2(:,i).*x,2) ...
                        + Mesh.mfPi(dphi3(:,i).*x,3);
        if nargout==1
            % matrix free mode, return D*xc and operators
            Dx1.D    = @(x) getGradientMF(x,Dxi,dim,'Du');
            Dx1.Dadj = @(x) getGradientMF(x,Dxiadj,dim,'DTu');
        else
            Dx1.D    = @(x) Dxi(x,1);
            Dx1.Dadj = @(x) Dxiadj(x,1);
            Dx2.D    = @(x) Dxi(x,2);
            Dx2.Dadj = @(x) Dxiadj(x,2);
        end
    case '3D-mf'
        % compute edges
        v1   = Mesh.mfPi(xn,1);
        v2   = Mesh.mfPi(xn,2);
        v3   = Mesh.mfPi(xn,3);
        v4   = Mesh.mfPi(xn,4);
        e1   =  v1 - v4;
        e2   =  v2 - v4;
        e3   =  v3 - v4;
        % compute inverse transformation to reference element
        cofA =  [
              e2(:,2).*e3(:,3)-e2(:,3).*e3(:,2), ...
            -(e1(:,2).*e3(:,3)-e1(:,3).*e3(:,2)), ...
              e1(:,2).*e2(:,3)-e1(:,3).*e2(:,2), ...
            -(e2(:,1).*e3(:,3)-e2(:,3).*e3(:,1)), ...
              e1(:,1).*e3(:,3)-e1(:,3).*e3(:,1), ...
            -(e1(:,1).*e2(:,3)-e1(:,3).*e2(:,1)), ...
              e2(:,1).*e3(:,2)-e2(:,2).*e3(:,1), ...
            -(e1(:,1).*e3(:,2)-e1(:,2).*e3(:,1)), ...
              e1(:,1).*e2(:,2)-e1(:,2).*e2(:,1), ...
            ];
        detA = e1(:,1).*cofA(:,1) + e2(:,1).*cofA(:,2)+ e3(:,1).*cofA(:,3);

        % compute gradients of basis functions
         dphi1 =   cofA(:,[1 4 7])./repmat(detA,[1 3]);
        dphi2 =   cofA(:,[2 5 8])./repmat(detA,[1 3]);
        dphi3 =   cofA(:,[3 6 9])./repmat(detA,[1 3]);
       dphi4 =   -dphi1 - dphi2 - dphi3;

        Dxi    = @(x,i) dphi1(:,i).*Mesh.mfPi(x,1) + dphi2(:,i).*Mesh.mfPi(x,2) ...
            + dphi3(:,i).*Mesh.mfPi(x,3) + dphi4(:,i).*Mesh.mfPi(x,4);
        Dxiadj = @(x,i) Mesh.mfPi(dphi1(:,i).*x,1) + Mesh.mfPi(dphi2(:,i).*x,2) ...
            + Mesh.mfPi(dphi3(:,i).*x,3) + Mesh.mfPi(dphi4(:,i).*x,4);

        if nargout==1
            % matrix free mode, return D*xc and operators
            Dx1.D    = @(x) getGradientMF(x,Dxi,dim,'Du');
            Dx1.Dadj = @(x) getGradientMF(x,Dxiadj,dim,'DTu');
        else
            Dx1.D    = @(x) Dxi(x,1);
            Dx1.Dadj = @(x) Dxiadj(x,1);
            Dx2.D    = @(x) Dxi(x,2);
            Dx2.Dadj = @(x) Dxiadj(x,2);
            Dx3.D    = @(x) Dxi(x,3);
            Dx3.Dadj = @(x) Dxiadj(x,3);
        end

    case '3D-mb'
        % compute edges
        e1   =  Mesh.P1*xn - Mesh.P4*xn;
        e2   =  Mesh.P2*xn - Mesh.P4*xn;
        e3   =  Mesh.P3*xn - Mesh.P4*xn;
        % compute inverse transformation to reference element
        cofA =  [
            e2(:,2).*e3(:,3)-e2(:,3).*e3(:,2), ...
            -(e1(:,2).*e3(:,3)-e1(:,3).*e3(:,2)), ...
            e1(:,2).*e2(:,3)-e1(:,3).*e2(:,2), ...
            -(e2(:,1).*e3(:,3)-e2(:,3).*e3(:,1)), ...
            e1(:,1).*e3(:,3)-e1(:,3).*e3(:,1), ...
            -(e1(:,1).*e2(:,3)-e1(:,3).*e2(:,1)), ...
            e2(:,1).*e3(:,2)-e2(:,2).*e3(:,1), ...
            -(e1(:,1).*e3(:,2)-e1(:,2).*e3(:,1)), ...
            e1(:,1).*e2(:,2)-e1(:,2).*e2(:,1), ...
            ];
        detA = e1(:,1).*cofA(:,1) + e2(:,1).*cofA(:,2)+ e3(:,1).*cofA(:,3);

        % compute gradients of basis functions
        dphi1 =   cofA(:,[1 4 7])./repmat(detA,[1 3]);
        dphi2 =   cofA(:,[2 5 8])./repmat(detA,[1 3]);
        dphi3 =   cofA(:,[3 6 9])./repmat(detA,[1 3]);
        dphi4 =   - dphi1 - dphi2 - dphi3;

        Dxi = @(i) sdiag(dphi1(:,i))*Mesh.P1 +  sdiag(dphi2(:,i))*Mesh.P2  ...
            + sdiag(dphi3(:,i))*Mesh.P3 +  sdiag(dphi4(:,i))*Mesh.P4;

        Dx1 = Dxi(1); Dx2 = Dxi(2); Dx3 = Dxi(3);

        if nargout==1, % return gradient operator for vector field
            B = [Dx1;Dx2;Dx3];
            Dx1 = blkdiag(B,B,B);
        end

    otherwise
        error('unknown flag: %s',flag);
end

function Du = getGradientMF(x,Dxi,dim,flag)
flag = [num2str(dim) 'D-' flag];
switch flag
    case '2D-Du'
        x = reshape(x,[],dim);
        Du  = [  Dxi(x(:,1),1);  Dxi(x(:,1),2);  Dxi(x(:,2),1);  Dxi(x(:,2),2); ];
    case '2D-DTu'
        x = reshape(x,[],4);
        Du  =  [Dxi(x(:,1),1)+Dxi(x(:,2),2); Dxi(x(:,3),1)+Dxi(x(:,4),2)];
    case '3D-Du'
        x = reshape(x,[],dim);
        Du  = [
                Dxi(x(:,1),1);  Dxi(x(:,1),2);  Dxi(x(:,1),3);...
                Dxi(x(:,2),1);  Dxi(x(:,2),2);  Dxi(x(:,2),3);...
                Dxi(x(:,3),1);  Dxi(x(:,3),2);  Dxi(x(:,3),3);...
               ];
    case '3D-DTu'
        x = reshape(x,[],9);
        Du  =  [
                    Dxi(x(:,1),1)+Dxi(x(:,2),2)+Dxi(x(:,3),3);...
                    Dxi(x(:,4),1)+Dxi(x(:,5),2)+Dxi(x(:,6),3);...
                    Dxi(x(:,7),1)+Dxi(x(:,8),2)+Dxi(x(:,9),3);...
                ];
    otherwise
        error('nyi');
end

% shortcut for sparse diagonal matrices
function D = sdiag(v)
D = spdiags(v(:),0,numel(v),numel(v));

function runMinimalExample
omega = [0 1 0 3]; m = [ 3 5];
Mesh = TriMesh1(omega,m);
xn = Mesh.xn + 6*1e-2 * randn(size(Mesh.xn));

D = feval(mfilename,Mesh,0);

figure(1); clf;
subplot(1,3,1);
triplot(Mesh.tri,xn(:,1),xn(:,2)); axis(omega); hold on;
plot(xn(:,1),xn(:,2),'or');
title('discretization');
subplot(1,3,2); spy(D);
title('gradient operator B');
subplot(1,3,3); spy(D'*D);
title('Laplace operator B''*B');

Dmf = feval(mfilename,Mesh,1);

err1 = norm(D*xn(:)-Dmf.D(xn))/norm(D*xn(:));
x = randn(size(xn)); y=randn(4*size(Mesh.tri,1),1);
err2 = abs(x(:)'*Dmf.Dadj(y)-y(:)'*Dmf.D(x));
fprintf('MFerror: %1.3e;  AdjointError: %1.3e OK? %d\n',err1,err2,max(err1,err2)<1e-12);
%========= 3D ==========
omega = [0 1 0 3 0 3]; m = [ 3 5 4];
Mesh = getTriangleMesh(omega,m,'matrixFree',0);
xn = Mesh.xn + 6*1e-2 * randn(size(Mesh.xn));

D = feval(mfilename,Mesh);

figure(1); clf;
subplot(1,3,1);
tetramesh(Mesh.tri,Mesh.xn,'FaceAlpha',.3); hold on;
plot3(xn(:,1),xn(:,2),xn(:,3),'or');
title('discretization');
subplot(1,3,2); spy(D);
title('gradient operator B');
subplot(1,3,3); spy(D'*D);
title('Laplace operator B''*B');

MeshMF = getTriangleMesh(omega,m,'matrixFree',1);
Dmf = feval(mfilename,MeshMF);

err1 = norm(D*xn(:)-Dmf.D(xn))/norm(D*xn(:));
x = randn(size(xn)); y=randn(9*size(Mesh.tri,1),1);
err2 = abs(x(:)'*Dmf.Dadj(y)-y(:)'*Dmf.D(x));
fprintf('MFerror: %1.3e AdjointError: %1.3e OK? %d\n',err1,err2,max(err1,err2)<1e-12);


text:
function dStates = Dynamics(States, Inputs, Params)

%This function is used to compute the dynamics of the tractor trailer truck
%system. The equations are derived by Derive_EoM.m

%States is a (4xN) matrix of states
%Inputs is a (2xN) matrix of inputs
%Params is a structure of parameters

%x = States(1,:);
%y = States(2,:);
th = States(3,:);
phi = States(4,:);

v = Inputs(1,:);
psi = Inputs(2,:)';

Lt = Params.Dyn.Lt;
Lc = Params.Dyn.Lc;

%[dStates, A, B] = Derive_EoM()

dStates =[...

                      -v.*cos(phi).*cos(psi).*sin(th);
                      v.*cos(phi).*cos(psi).*cos(th);
                         (v.*cos(psi).*sin(phi))/Lt;
 (v.*(Lt.*sin(psi) - Lc.*cos(psi).*sin(phi)))/(Lc*Lt)];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function quality2 = tetrahedron_quality2_3d ( tetra )

%*****************************************************************************80
%
%% TETRAHEDRON_QUALITY2_3D: "quality" of a tetrahedron in 3D.
%
%  Discussion:
%
%    The quality measure #2 of a tetrahedron is:
%
%      QUALITY2 = 2 * sqrt ( 6 ) * RIN / LMAX
%
%    where
%
%      RIN = radius of the inscribed sphere;
%      LMAX = length of longest side of the tetrahedron.
%
%    An equilateral tetrahredron achieves the maximum possible quality of 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du and Desheng Wang,
%    The Optimal Centroidal Voronoi Tesselations and the Gersho's
%    Conjecture in the Three-Dimensional Space,
%    Computers and Mathematics with Applications,
%    Volume 49, 2005, pages 1355-1373.
%
%  Parameters:
%
%    Input, real TETRA(3,4), the tetrahedron vertices.
%
%    Output, real QUALITY2, the quality of the tetrahedron.
%
  dim_num = 3;

  edge_length(1:6) = tetrahedron_edge_length_3d ( tetra );

  l_max = max ( edge_length(1:6) );

  [ r_in, pc ] = tetrahedron_insphere_3d ( tetra );

  quality2 = 2.0 * sqrt ( 6.0 ) * r_in / l_max;

  return
end


text:
function [y,x]=newlagmatrix(x,nlags,c)
% Construction of a matrix of lags (X) and a vector (Y) for use in an
% autoregression
%
% USAGE:
%     [Y,X]=newlagmatrix(Y,NLAGS,C)
%
% INPUTS:
%     Y     - The dependant variable(Tx1)
%     NLAGS - The number of lags(scalar)
%     C     - [OPTIONAL] 1 if you want to include a constant; 0 is default
%
% OUTPUTS:
%     Y     - (n-p) by 1 vector of contemporaneous values
%     X     - (n-p) by p (or p+1 if c=1) matrix of lags and possibly constant.  The
%               matrix is of the form [constant(if included) t-1 t-2 ... t-p]
%
% COMMENTS:
%     Original name, 'lagmatrix' conflicts with a Matlab file, so newlagmatrix
%     replaces original

% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 5    Date: 12/1/2005

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 2
        c=0;
    case 3
        %Nothing
    otherwise
        error('newlagmatrix requires either 2 or 3 arguments.')
end

[T,K]=size(x);
if min(T,K)~=1
    error('X must be a vector')
end

if T<K
    x=x'; %Transpose if a column vector is passed
end

if  length(c)~=1 || any(c~=0 && any(c~=1))
    error('C must be 1 or 0')
end

if length(nlags)~=1 || floor(nlags)~=nlags || any(nlags<0)
    error('NLAGS must be a positive integer')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nlags>0
    nlags=nlags+1;
    newX=[x;zeros(nlags,1)];
    lagmatrix=repmat(newX,nlags,1);
    lagmatrix=reshape(lagmatrix(1:size(lagmatrix,1)-nlags),T+nlags-1,nlags);
    lagmatrix=lagmatrix(nlags:T,:);
    y=lagmatrix(:,1);
    x=lagmatrix(:,2:nlags);
    if c==1
        x=[ones(size(x,1),1) x];
    end
else
    if c==1
        y=x;
        x=ones(T,1);
    else
        y=x;
        x=[];
    end
end

author: bashtage

repo: mfe-toolbox


text:
function i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d: ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
%  jSpectral:  Multitaper spectral analysis, and other time series tools
%
%  Multitaper spectral analysis
%   sleptap    - Calculate Slepian tapers.
%   mspec      - Multitaper power and cross spectra.
%   mconf      - Confidence intervals for the multitaper spectral estimate.
%
%  Multitaper polarization analysis
%   msvd       - Singular value decomposition for polarization analysis.
%   polparams  - Spectral matrix polarization parameters.
%   specdiag   - Diagonalize a 2 x 2 spectral matrix.
%
% Assorted other transforms
%   slidetrans  - Sliding-window ('moving-window') Fourier transform.
%   anatrans    - Analytic part of signal.
%   wigdist     - Wigner distribution (alias-free algorithm).
%
% Time series analysis utilities
%   doublen     - Interpolates a time series to double its length.
%   fourier     - The one-sided Fourier frequencies for a given length time series.
%   sampletimes - Computes mean sampling intervals and their statistics.
%
%  Plotting tools
%   twospecplot - Plots a pair of rotary or Cartesian spectra.
%
%  See also jWavelet, jEllipse, jMatern.

%   Low-level functions
%   timeseries_boundary - Apply boundary conditions to data before transform.

help jspectral

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


validInputSize = [3,3,18,1];
layer = yolov3Layer('yolov3', rand(3,2),1,1,[416,416],'default');
checkLayer(layer,validInputSize,'ObservationDimension' ,4)

%%
layer = mishLayer('mishACT');
validInputSize = [24 24 20];
checkLayer(layer,validInputSize,'ObservationDimension',4)


%% NISwGSP-03_SantaMaria-all
% %{
imfolder = 'images\NISwGSP-03_SantaMaria-all';
im_n = 8;
imfile = cell(im_n,1);
for ii = 1:im_n
    imfile{ii} = sprintf('%s\\%02d.jpg', imfolder, ii);
end

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, [], 0, 'equi', 0.02, imfolder );
%}


text:
%DATASETS Info on the dataset class construction for PRTools
%
% This is not a command, just an information file.
%
% Datasets in PRTools are in the MATLAB language defined as objects of the
% class PRDATASET. Below, the words 'object' and 'class' are used in the pattern
% recognition sense.
%
% A dataset is a set consisting of M objects, each described by K features.
% In PRTools, such a dataset is represented by a M x K matrix: M rows, each
% containing an object vector of K elements. Usually, a dataset is labeled.
% An example of a definition is:
%
%  DATA = [RAND(3,2) ; RAND(3,2)+0.5];
%  LABS = ['A';'A';'A';'B';'B';'B'];
%  A = PRDATASET(DATA,LABS)
%
% which defines a [6 x 2] dataset with 2 classes.
%
% The [6 x 2] data matrix (6 objects given by 2 features) is accompanied by
% labels, assigning each of the objects to one of the two classes A and B.
% Class labels can be numbers or strings and should always be given as rows
% in the label list. A lable may also have the value NaN or may be an empty
% string, indicating an ulabeled object. If the label list is not given,
% all objects are marked as unlabeled.
%
% Various other types of information can be stored in a dataset. The most
% simple way to get an overview is by typing:
%
%  STRUCT(A)
%
% which for the above example displays the following:
%
%         DATA: [6x2 double]
%      LABLIST: [2x1 double]
%         NLAB: [6x1 double]
%      LABTYPE: 'crisp'
%      TARGETS: []
%      FEATLAB: [2x1 double]
%      FEATDOM: {1x2 cell}
%        PRIOR: []
%         COST: []
%      OBJSIZE: 6
%     FEATSIZE: 2
%        IDENT: {6x1 cell}
%      VERSION: {1x2 cell}
%         NAME: []
%         USER: []
%
% These fields have the following meaning:
%
% DATA     : an array containing the objects (the rows) represented by
%            features (the columns). In the software and help-files, the number
%            of objects is usually denoted by M and the number of features is
%            denoted by K. So, DATA has the size of [M,K]. This is also defined
%            as the size of the entire dataset.
% LABLIST  : The names of the classes, stored row-wise. These class names
%            should be integers, strings or cells of strings. Mixtures of
%            these are not supported. LABLIST has as many rows as there are
%            classes. This number is usually denoted by C. LABLIST is
%            constructed from the set of LABELS given in the DATASET command
%            by determining the unique names while ordering them alphabetically.
% NLAB     : an [M x 1] vector of integers between 1 and C, defining for each
%            of the M objects its class. They are indexing LABLIST.
% LABTYPE  : 'CRISP', 'SOFT' or 'TARGETS' are the three possible label types.
%            In case of 'CRISP' labels, a unique class, defined by NLAB, is
%            assigned to each object, pointing to the class names given in
%            LABLIST.
%            For 'SOFT' labels, each object has a corresponding vector of C
%            numbers between 0 and 1 indicating its membership (or confidence
%            or posterior probability) of each of the C classes. These numbers
%            are stored in the array TARGETS of the size M x C. They don't
%            necessarily sum to one for individual row vectors.
%            Labels of type 'TARGETS' are in fact no labels, but merely target
%            vectors of length C. The values are again stored in TARGETS and
%            are not restricted in value.
% TARGETS  : [M,C] array storing the values of the soft labels or targets.
% FEATLAB  : A label list (like LABLIST) of K rows storing the names of the
%            features.
% FEATDOM  : A cell array describing for each feature its domain.
% PRIOR    : Vector of length C storing the class prior probabilities. They
%            should sum to one. If PRIOR is empty ([]) it is assumed that the
%            class prior probabilities correspond to the class frequencies.
% COST     : Classification cost matrix. COST(I,J) are the costs
%            of classifying an object from class I as class J. Column C+1
%            generates an alternative reject class and may be omitted,
%            yielding a size of [C,C]. An empty cost matrix, COST = []
%            (default) is interpreted as COST = ONES(C) - EYE(C) (identical
%            costs of misclassification).
% OBJSIZE  : The number of objects, M. In case the objects are related to a
%            n-dimensional structure, OBJSIZE is a vector of length n, storing
%            the size of this structure. For instance, if the objects are pixels
%            in a [20 x 16] image, then OBJSIZE = [20,16] and M = 320.
% FEATSIZE : The number of features, K. In case the features are related to
%            an n-dimensional structure, FEATSIZE is a vector of length n,
%            storing the size of this structure. For instance, if the features
%            are pixels in a [20 x 16] image, then FEATSIZE = [20,16] and
%            K = 320.
% IDENT    : A cell array of M elements storing indicators of the M objects.
%            They are initialized by integers 1:M.
% VERSION  : Some information related to the version of PRTools used for
%            defining the dataset.
% NAME     : A character string naming the dataset, possibly used to annotate
%            related graphics.
% USER     : Free field for the user, not used by PRTools.
%
%
% The fields can be set by commands like SETDATA, SETFEATLAB, SETLABELS,
% see below for a complete list.
% Note that there is no field LABELS in the DATASET definition. Labels are
% converted to NLAB and LABLIST. The command SETLABELS however exists and
% takes care of the conversion.
%
% The data and information stored in a dataset can be retrieved as follows:
%
% - By DOUBLE(A) and by +A, the content of A.DATA is returned.
% - [N,LABLIST] = CLASSSIZES(A);
%   It returns the numbers of objects per class and the class names stored
%   in LABLIST.
% - By DISPLAY(A), it writes the size of the dataset, the number of classes
%   and the label type on the terminal screen.
% - By SIZE(A), it returns the size of A.DATA: numbers of objects and features.
% - By SCATTERD(A), it makes a scatter plot of a dataset.
% - By SHOW(A), it may be used to display images that are stored as features
%   or as objects in a dataset.
% - By commands like: GETDATA, GETFEATLAB, etcetera, see below. With some
%   exceptions they point to a single dataset field. E.g. GETSIZE(A) returns
%   [M,K,C]. A aet of commands does not return data, but instead they return
%   indices to objects that have specific identifiers, labels or class indices:
%   FINDIDENT, FINDLABELS, FINDNLAB.
%
% Many standard MATLAB operations and a number of general MATLAB commands have
% been overloaded for variables of the DATASET type.
%
% SEE ALSO (PRTools Guide)
% PRDATASET, DATA2IM, OBJ2FEAT, FEAT2OBJ, IM2FEAT, IM2OBJ, DATAIM
% SETDATA, SETFEATLAB, SETFEATDOM, SETFEATSIZE, SETIDENT, SETLABELS,
% SETLABLIST, SETLABTYPE, SETNAME, SETNLAB, SETOBJSIZE, SETPRIOR, SETCOST,
% SETTARGETS, SETUSER, SETLABLISTNAMES, SETVERSION
% GETDATA, GETFEATLAB, GETFEATDOM, GETFEATSIZE, GETIDENT, GETLABELS,
% GETLABLIST, GETLABTYPE, GETNAME, GETNLAB, GETOBJSIZE, GETPRIOR, GETCOST,
% GETSIZE, GETTARGETS, GETUSER, GETVERSION, GETCLASSI, GETLABLISTNAMES,
% FINDIDENT, FINDLABELS, FINDNLAB

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

author: marianux

repo: ecg-kit

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


text:
% eegdrawg() - subroutine used by eegplotgold() to plot data.
%
% Author: Colin Humphries, CNL, Salk Institute, La Jolla, 7/96

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

% 4-4-97 shortened name to eegdrawq() -sm
% 4-7-97 allowed data names other than 'data' -ch
% 01-25-02 reformated help & license -ad

function y = eegdrawg(fighandle)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extract variables from figure and axes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

userdata = get(fighandle,'UserData');
samplerate = userdata(1);
PLOT_TIME = userdata(2);
spacing_var = userdata(3);
time = userdata(4);
maxtime = userdata(5);
axhandle = userdata(6);
plotcolor = userdata(7);
disp_scale = userdata(12);
colors = ['y','w'];
dataname = get(axhandle,'UserData');
eval(['global ',dataname])

cla % Clear figure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define internal variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eval(['[chans,frames] = size(',dataname,');']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Label x-axis
% This routine relabels the x-axis based on the new value of time.
% Labels are placed on one second intervals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Xlab = num2str(time);
for j = 1:1:PLOT_TIME
   Q = num2str(time+j);
   Xlab = str2mat(Xlab, Q);
end
set (gca, 'Ytick', 0:spacing_var:chans*spacing_var)
set (gca, 'XTickLabels', Xlab)
set (gca, 'XTick',(0:samplerate:PLOT_TIME*samplerate))

axis([0 PLOT_TIME*samplerate 0 (chans+1)*spacing_var]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plotting routine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:chans                 %repeat for each channel
   if (maxtime-time>PLOT_TIME)
      eval(['F = ',dataname,'(chans-i+1,(time*samplerate)+1:((time+PLOT_TIME)*samplerate));'])
   else
      eval(['F = ',dataname,'(chans-i+1,(time*samplerate)+1:(maxtime*samplerate));'])
   end
   F = F - mean(F) + i*spacing_var;
   plot (F,'clipping','off','color',colors(plotcolor))
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Scaling I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if disp_scale == 1
   line([(PLOT_TIME+.3)*samplerate,(PLOT_TIME+.3)*samplerate],[1.5*spacing_var 2.5*spacing_var],'clipping','off','color','w')
   line([(PLOT_TIME+.3)*samplerate-10,(PLOT_TIME+.3)*samplerate+10],[2.5*spacing_var,2.5*spacing_var],'clipping','off','color','w')
   line([(PLOT_TIME+.3)*samplerate-10,(PLOT_TIME+.3)*samplerate+10],[1.5*spacing_var,1.5*spacing_var],'clipping','off','color','w')
   text((PLOT_TIME+.5)*samplerate,2*spacing_var,num2str(round(spacing_var)),'clipping','off')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set slider=edit
% This routine resets the value of the user controls so that
%    they agree.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   H = findobj(fighandle,'style','slider');
   D = findobj(fighandle,'style','edit');
   set (D, 'string', num2str(time));

author: PatternRecognition

repo: OpenBMI


text:
function plot_vel(gnss_time,gnss_vel, matlab_time,matlab_vel, mcu_time,mcu_vel)
%速度估计曲线
    for i=2:2:nargin
        if i==2
%             gnss_time = gnss_time/60;
            figure('name', '速度估计曲线');
            subplot(3,1,1);
            plot(gnss_time, gnss_vel(:,1), 'r', 'Marker','.'); hold on; grid on;
            xlim([gnss_time(1) gnss_time(end)]);
            ylabel('东向速度(m/s)'); legend('GNSS');
            subplot(3,1,2);
            plot(gnss_time, gnss_vel(:,2), 'r', 'Marker','.'); hold on; grid on;
            xlim([gnss_time(1) gnss_time(end)]);
            ylabel('北向速度(m/s)'); legend('GNSS');
            subplot(3,1,3);
            plot(gnss_time, gnss_vel(:,3), 'r', 'Marker','.'); hold on; grid on;
            xlim([gnss_time(1) gnss_time(end)]);
            xlabel('时间(s)'); ylabel('天向速度(m/s)'); legend('GNSS');
        elseif i==4
%             matlab_time = matlab_time/60;
            subplot(3,1,1);
            plot(matlab_time, matlab_vel(:,1), 'b', 'Marker','.');
            legend('GNSS', 'MATLAB');
            subplot(3,1,2);
            plot(matlab_time, matlab_vel(:,2), 'b', 'Marker','.');
            legend('GNSS', 'MATLAB');
            subplot(3,1,3);
            plot(matlab_time, matlab_vel(:,3), 'b', 'Marker','.');
            legend('GNSS', 'MATLAB');
        elseif i==6
%             mcu_time = mcu_time/60;
            subplot(3,1,1);
            plot(mcu_time, mcu_vel(:,1), 'm', 'Marker','.');
            legend('GNSS', 'MATLAB', 'MCU');
            subplot(3,1,2);
            plot(mcu_time, mcu_vel(:,2), 'm', 'Marker','.');
            legend('GNSS', 'MATLAB', 'MCU');
            subplot(3,1,3);
            plot(mcu_time, mcu_vel(:,3), 'm', 'Marker','.');
            legend('GNSS', 'MATLAB', 'MCU');
        end
    end

    sgtitle('速度估计曲线');
    set(gcf, 'Units', 'normalized', 'Position', [0.025, 0.05, 0.95, 0.85]);
end

author: yandld


text:
function pass = test_coeffs2diskfun( )
% test diskfun coeff2diskfun() command

tol = 10*chebfunpref().cheb2Prefs.chebfun2eps;

f = diskfun.coeffs2diskfun(0);
pass(1) = iszero(f);

f = diskfun(@(t,r) r.^2, 'polar');
g = diskfun.coeffs2diskfun(coeffs2(f));
pass(2) = norm(f-g) < tol;

f = diskfun(@(t,r) r.*sin(t), 'polar');
g = diskfun.coeffs2diskfun(coeffs2(f));
pass(3) = norm(f-g) < tol;

f = diskfun(@(x,y) exp(-10*((x-0.5/sqrt(2)).^2 + (y-0.5/sqrt(2)).^2)));
g = diskfun.coeffs2diskfun(coeffs2(f));
pass(4) = norm(f-g);

f = diskfun(@(t,r) r.*sin(t), 'polar');
c = 1i/(2)*[0 0 0 ; 1 0 -1 ];
pass(5) = norm( f - diskfun.coeffs2diskfun(c) ) < tol;

f = diskfun(@(t,r) r.^3.*cos(3*t)+r.^2.*(sin(t)).^2, 'polar');
c = (1/8)*[0 -1 0 2 0 -1 0; 3 0 0 0 0 0 3; 0 -1 0 2 0 -1 0;1 0 0 0 0 0 1];
pass(6) = norm( f - diskfun.coeffs2diskfun(c) ) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function sys = dual(X)
%DUAL Extract dual variable
%
%   Z = DUAL(F)     Returns the dual variable for the constraint F
%
%   See also SOLVESDP, DUALIZE

sys = dual(lmi(X));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function fmri_loc2_GLMcontrasts

% LO Localizer scans processing code
%
% USER IMPORTANT!!!: Before running this script, make sure that all loloc
% scans have the correct annotation fields and scan groups (are grouped
% correctly). Otherwise, fmri_FFA_checkDataTypesForGLM will fail.
%
% 1. Set GLM parameters for all scans
% 2. Check LO localizer experiments for:
%    a. Datatype and scan numbers
%    b. Record parfiles
%    c. Make sure scan groupings are correct
% 3. Run GLM for LO localizer experiments
% 4. Run contrasts for the LO localizer experiments
%
% Calls the functions:
% * fmri_loc2_setGLMparams (subfunction)
% * fmri_loc2_checkDataTypesForGLM
% * fmri_loc2_computeContrastsLocalizer
%
%
% DY 06/05/2008
% Modified 06/18/2008: there are calls to er_setParams and saveSessions all
% over the place, which results in the params for the MotionComp_RefScan1
% being changed for the worse. I'm not sure how to efficiently prevent this
% from happening, so instead, I let it happen, and just fix the ANNOTATION
% and PARFILE fields at the end, since accuracy here is important.
%
% Code taken from fmri_FFA_motionCompFix
% AL 11/08 Modified for loc2 GLM

clear all

% Set directory and subject list
% Alternatively: subs = {'sess1','sess2','sess3',...};
if ispc
    fmriDir = 'W:\projects\Kids\fmri\localizer\';
else
    %fmriDir = '/biac1/kgs/projects/Kids/fmri/localizer/';
    fmriDir='/biac1/kgs/projects/Kids/fmri/code/Alina/Subject/';
end
cd(fmriDir); s=dir('*0*'); subs={s.name};

% Start a log text file to document successes and failures in preprocessing
dateAndTime=datestr(now); dateAndTime(12)='_'; dateAndTime(15)='h';dateAndTime(18)='m';
logDir=fullfile(fmriDir,'logs');
if ~isdir(logDir), mkdir(logDir), end;
logFile = fullfile(fmriDir,'logs',['LO2_Log_' dateAndTime '.txt']);
fid=fopen(logFile,'w');
[junk,logFileName]=fileparts(logFile)
scriptName = 'fmri_loc2_GLMcontrasts.m';
startTime = clock;

% 1. Set GLM parameters
params=fmri_loc2_setGLMparams(fid,logFileName,scriptName); % Subfunction within this script

for ii=1:length(subs)

    % Get rid of all mrVista globals in workspace. Put this at the top,
    % just in case we break from the rest of the loop.
    clear dataTYPES mrSESSION vANATOMYPATH HOMEDIR hI newDt newScan

    thisDir = fullfile(fmriDir,subs{ii}); cd(thisDir);
    fprintf(fid,'\n ------------------------------------------ \n');
    fprintf(fid,'Processing %s\n\n',thisDir);
    fprintf('Processing %s\n\n',thisDir);

    %% Check if the subject already has been processed by script (and thus
    %% should not be re-entered to this script). If they do, skip to the
    %% next subject by using CONTINUE.
    load mrSESSION.mat
    theglm=find(strcmp({dataTYPES.name},'GLMs'));
     if ~isempty(theglm)
        checkScans=length(dataTYPES(theglm).eventAnalysisParams);
        alreadyDone=[];
        for xx=1:checkScans
            tmp=(strfind(lower(dataTYPES(theglm).eventAnalysisParams(xx).annotation),'lo localizer run 2'));
            alreadyDone=[alreadyDone tmp];
        end
        if(~isempty(alreadyDone))
            fprintf(fid,'Error: Subject already processed (has GLM with ''LO localizer Run 2'' in annotation field), skipping... \n\n');
            fprintf('Error: Subject already processed (has GLM with ''LO localizer Run 2'' in annotation field), skipping... \n\n');
            continue
        end
    end

    %% 2a/b/c. Find datatype number for MotionComp_RefScan1, find scan
    %% numbers for lo localizer scans, record parfiles for these scans and
    %% make sure scan groupings are correct. If everything checks out,
    %% proceed with GLM. If not, SKIP THIS SUBJECT.
    [go,s2]=fmri_loc2_checkDataTypesForGLM(fid,dataTYPES);
    %clear dataTYPES mrSESSION vANATOMYPATH

    % If there was a problem in step 3, skip GLM and proceed to next
    % subject
    if go==0
        continue

    %% 3. Otherwise, run GLM (to debug: fid=1)
    elseif go==1
        try
            % Initialize a hidden view of the 'inplane' data (like invoking
            % mrVista without a GUI, dataTYPES/etc in workspace)
            hI = initHiddenInplane(s1.mc, s1.scan);
            params.parfiles={s1.parfile};
            er_setParams(hI, params, s2.scan, s2.mc); % write params as eventAnalysisParams fields
            tic; [hI, lo_scan] = applyGlm(hI,s2.mc, s2.scangroup,params);
            time=toc;
            fprintf(fid,'GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,lo_scan,floor(time/60),mod(time,60));
            fprintf('GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,lo_scan,floor(time/60),mod(time,60));
            glmOK=1;
        catch
            fprintf(fid,'FAILURE TO SUCCESSFULLY RUN GLM \n');
            fprintf('FAILURE TO SUCCESSFULLY RUN GLM \n');
            theerror=lasterror;
            fprintf(fid,'%s \n\n',theerror.message);
            fprintf('%s \n\n',theerror.message);
            glmOK=0;
        end

        %% 4. Compute contrast maps
        if glmOK

            try
                fmri_loc2_computeContrastsLocalizer(hI, lo_scan)
                fprintf(fid,'Computed contrasts using fmri_FFA_computeContrastsLocalizer \n');
                fprintf('Computed contrasts using fmri_FFA_computeContrastsLocalizer \n');
            catch
                fprintf(fid,'FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                fprintf('FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                theerror=lasterror;
                fprintf(fid,'%s \n\n',theerror.message);
                fprintf('%s \n\n',theerror.message);

            end
        end
    end

    % Fix the annotation and parfile fields for the MotionComp_RefScan
    % dataTYPE, since these are undesirably changed (collateral damage).
    dataTYPES(s2.mc).eventAnalysisParams(s2.scan).annotation='loloc_run2';
    par2=dataTYPES(s2.mc).scanParams(s2.scan).parfile;
    dataTYPES(s2.mc).eventAnalysisParams(s2.scan).parfiles=par2;
    save(fullfile(thisDir,'mrSESSION.mat'),'mrSESSION','dataTYPES', '-append');
    clear par1 par2

end

totalTime=etime(clock,startTime);

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'Total running time for script: %f minutes \n',totalTime/60);
fprintf('Total running time for script: %f minutes \n',totalTime/60);

fclose(fid); % close out the log file

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params=fmri_loc2_setGLMparams(fid,logFileName,scriptName)

% Set GLM parameters: We set the PARAMS struct fields to the values on the
% wiki, which are further explained there. Also we print the values to the
% log. http://vpnl.stanford.edu/internal/wiki/index.php/ROIanalysis#GLMs

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'\n GLM PARAMETERS \n');
fprintf(fid,' ------------------------------------------ \n');

params=er_defaultParams % Initialize param struct with defaults
p.al='params.alpha=0.0500'; eval(p.al); fprintf(fid,'%s \n',p.al);
p.at='params.ampType=''betas'''; eval(p.at); fprintf(fid,'%s \n',p.at);
p.an='params.annotation=''LO localizer Run 2'''; eval(p.an); fprintf(fid,'%s \n',p.an);
p.ap='params.assignParfiles=0'; eval(p.ap); fprintf(fid,'%s \n',p.ap);
p.bp='params.bslPeriod=[-6 -5 -4 -3 -2 -1 0]'; eval(p.bp); fprintf(fid,'%s \n',p.bp);
p.dt='params.detrend=1'; eval(p.dt); fprintf(fid,'%s \n',p.dt);
p.df='params.detrendFrames=20;'; eval(p.df); fprintf(fid,'%s \n',p.df);
p.ea='params.eventAnalysis=1'; eval(p.ea); fprintf(fid,'%s \n',p.ea);
p.eb='params.eventsPerBlock=6'; eval(p.eb); fprintf(fid,'%s \n',p.eb);
p.fp='params.framePeriod=2'; eval(p.fp); fprintf(fid,'%s \n',p.fp);
p.gt='params.glmHRF=3'; eval(p.gt); fprintf(fid,'%s \n',p.gt);
p.gp='params.glmHRF_params=[6 16 1 1 6 0 24]'; eval(p.gp); fprintf(fid,'%s \n',p.gp);
p.gw='params.glmWhiten=0'; eval(p.gw); fprintf(fid,'%s \n',p.gw); % temporal autocorrelation
p.ic='params.inhomoCorrect=1'; eval(p.ic); fprintf(fid,'%s \n',p.ic);
p.lp='params.lowPassFilter=0'; eval(p.lp); fprintf(fid,'%s \n',p.lp);
p.nb='params.normBsl=1'; eval(p.nb); fprintf(fid,'%s \n',p.nb);
p.od='params.onsetDelta=0'; eval(p.od); fprintf(fid,'%s \n',p.od);
p.pf='params.parfiles='''''; eval(p.pf); fprintf(fid,'%s \n',p.pf); % fill in after checkDataTypes...
p.pp='params.peakPeriod=[4 5 6 7 8 9 10 11 12 13 14]'; eval(p.pp); fprintf(fid,'%s \n',p.pp);
p.sh='params.setHRFParams=0'; eval(p.sh); fprintf(fid,'%s \n',p.sh);
p.sc='params.snrConds=1'; eval(p.sc); fprintf(fid,'%s \n',p.sc);
p.tn='params.temporalNormalization=0'; eval(p.tn); fprintf(fid,'%s \n',p.tn);
p.tw='params.timeWindow=[-6:1:24]'; eval(p.tw); fprintf(fid,'%s \n',p.tw);
p.log=['params.logFile=''' logFileName '''']; eval(p.log); fprintf(fid,'%s \n',p.log)
p.script=['params.scriptName=''' scriptName '''']; eval(p.script); fprintf(fid,'%s \n',p.script)

% params.framePeriod = dataTYPES(dt).scanParams(scan).framePeriod;


text:
% Test file for unbndfun/createMap.m.

function pass = test_createMap(pref)

tol = 1e1*eps;

% Doubly-unbounded domain:
dom = [-Inf Inf];
map = unbndfun.createMap(dom);

pass(1) = ( all(map.For([-1 1]) == [-Inf Inf]) ) && ...
    ( all(map.Inv([-1e100 1e100]) - [-1 1] < tol) );

% Singly-unbounded domain (left):
dom = [-Inf 3];
map = unbndfun.createMap(dom);
pass(2) = ( all(map.For([-1 1]) == [-Inf 3]) ) && ...
    ( all(map.Inv([-1e100 3]) == [-1 1]) );

% Singly-unbounded domain (right):
dom = [-100 Inf];
map = unbndfun.createMap(dom);
pass(3) = ( all(map.For([-1 1]) == [-100 Inf]) ) && ...
    ( all(map.Inv([-100 1e100]) == [-1 1]) );

end


text:
function images = loadMNISTImages(filename)
%loadMNISTImages returns a 28x28x[number of MNIST images] matrix containing
%the raw MNIST images

fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);

magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2051, ['Bad magic number in ', filename, '']);

numImages = fread(fp, 1, 'int32', 0, 'ieee-be');
numRows = fread(fp, 1, 'int32', 0, 'ieee-be');
numCols = fread(fp, 1, 'int32', 0, 'ieee-be');

images = fread(fp, inf, 'unsigned char=>unsigned char');
images = reshape(images, numCols, numRows, numImages);
images = permute(images,[2 1 3]);

fclose(fp);

% Reshape to #pixels x #examples
images = reshape(images, size(images, 1) * size(images, 2), size(images, 3));
% Convert to double and rescale to [0,1]
images = double(images) / 255;

end

author: xuzhenqi

repo: cnn

save_path: github-repos/MATLAB/xuzhenqi-cnn