text:
function P = legendrepol(N, a, b)
%-------------------------------------------------------------------------
%
%       P = legendrepol(N)
%
% returns the Legendre polynomials up to order N, orthogonal on
% the interval [-1,1]. Each row contain the polynomial coefficients in
% descending order, i.e. the first row is P0 = [0 0 .... 1].
%
%       P = legendrepol(N),xmin,xmax)
%
% returns the Legendre polynomial orthogonal on [xmin, xmax].
%-------------------------------------------------------------------------

% T. Wik, April 2012

% L is the Legendre polynomials orthogonal on (-1,1):
if N<1
	L = 1;
else
	L 		   = zeros(N+1,N+1);
	L(1,N+1)   = 1;              % L0
	L(2,N:N+1) = [1 0];          % L1
   for n=1:N-1                   % L2 to LN
   	L(n+2,N-n:N+1) 	= 1/(n+1)* ...
                   	 ((2*n+1)*[L(n+1,N+1-n:N+1) 0]-n*[0 0 L(n,N+2-n:N+1)]);
   end
end

if nargin == 3
    if a > b
        error('legendrepol:wrongArguments','xmin > xmax !')
    end
% Transform such that P is orthogonal on (xmin,xmax):
    P = zeros(N+1);
    A = zeros(N+1);

    for n = 1:N+1
        A(end-n+1,end-n+1:end) = binomial(n-1,2/(b-a),-(a+b)/(b-a));
        Lambda  = diag(L(n,:));
        P(n,:) = sum(Lambda*A);
    end
elseif nargin == 1
        P = L;
else
    error('legendrepol:wrongArguments','Use one or three arguments!')
end

function P = binomial(n,a,b)
%------------------------------------------------------
%
% function P = binomial(n,a,b)
%
% Returns the polynomial coefficients p(i) for
%
% P(x) = (ax + b)^n
%      = p(1)x^n + p(2)x^(n-1) + ... + p(n)x + p(n+1)
%
%------------------------------------------------------

% First binomial coefficients for (x+1)^n

if n==0, P=1; return, end
p = zeros(n,n+1);

p(1,1:2) = [1 1];

for j = 2:n
    p(j,1) = 1;
    for k = 2:n
        p(j,k) = p(j-1,k-1) + p(j-1,k);
    end
    p(j,n+1) = 1;
end

% Adjust for a and b

for k = 1:n+1;
    P(k) = p(n,k)*a^(n-k+1)*b^(k-1);
end


text:
 function [PhaseLockingData] = bz_PhaseModulation(varargin)
% USAGE
%[PhaseLockingData] = bz_PhaseModulation(varargin)
%
% INPUTS
% spikes        -spike time cellinfo struct
%
% lfp           -lfp struct with a single channel from bz_GetLFP()
%
% passband      -frequency range for phase modulation [lowHz highHz] form
%
% intervals     -(optional) may specify timespans over which to calculate
%               phase modulation.  Formats accepted: tstoolbox intervalSet
%               or a 2column matrix of [starts stops] in seconds
%
% samplingRate  -specifies lfp sampling frequency default=1250
%
% method        -method selection for how to generate phase,
%               possibilties are: 'hilbert' (default) or 'wavelet'
%
% powerThresh   -integer power threshold to use as cut off,
%               measured in standard deviations (default = 2)
%
% plotting      -logical if you want to plot, false if not, default=true
%
% saveMat       -logical to save cellinfo .mat file with results, default=false
%
%
% OUTPUTS
%
% phasedistros  - Spike distribution perecentages for each cell in each bin
%               specified by phasebins
%
% phasebins     - 180 bins spanning from 0 to 2pi
%
% phasestats    - ncellsx1 structure array with following (via
%                 CircularDistribution.m from FMAToolbox)
%                    phasestats.m        mean angle
%                    phasestats.mode     distribution mode
%                    phasestats.k        concentration
%                    phasestats.p        p-value for Rayleigh test
%                    phasestats.r        mean resultant length
%
%
% Calculates distribution of spikes over various phases from a specified
% cycle of an lfp vector.   Phase 0 means peak of lfp wave.
%
% Brendon Watson 2015
% edited by david tingley, 2017

%% defaults
p = inputParser;
addRequired(p,'spikes',@bz_isCellInfo);
addRequired(p,'lfp',@bz_isLFP);
addRequired(p,'passband',@isnumeric)
addParameter(p,'intervals',[0 inf],@isnumeric)
addParameter(p,'samplingRate',1250,@isnumeric)
addParameter(p,'method','hilbert',@isstr)
addParameter(p,'plotting',true,@islogical)
addParameter(p,'numBins',180,@isnumeric)
addParameter(p,'powerThresh',2,@isnumeric)
addParameter(p,'saveMat',false,@islogical)

parse(p,varargin{:})

spikes = p.Results.spikes;
lfp = p.Results.lfp;
passband = p.Results.passband;

intervals = p.Results.intervals; % interval(s) over which to calculate
samplingRate = p.Results.samplingRate; % sampling rate of continuous signal (LFP)
method = p.Results.method;
plotting = p.Results.plotting;
numBins = p.Results.numBins;
powerThresh = p.Results.powerThresh;
saveMat = p.Results.saveMat;

%% Get phase for every time point in LFP
switch lower(method)
    case ('hilbert')
        [b a] = butter(3,[passband(1)/(samplingRate/2) passband(2)/(samplingRate/2)],'bandpass'); % order 3
%         [b a] = cheby2(4,20,passband/(samplingRate/2));
        filt = FiltFiltM(b,a,double(lfp.data(:,1)));
        power = fastrms(filt,ceil(samplingRate./passband(1)));  % approximate power is frequency band
        hilb = hilbert(filt);
        lfpphase = mod(angle(hilb),2*pi);
        clear fil
    case ('wavelet')% Use Wavelet transform to calulate the signal phases
%         nvoice = 12;
%         freqlist= 2.^(log2(passband(1)):1/nvoice:log2(passband(2)));
%         error('awt_freqlist, where did this come from?')
%         wt = awt_freqlist(double(lfp.data(:,1)), samplingRate, freqlist);
%         amp = (real(wt).^2 + imag(wt).^2).^.5;
%         phase = atan2(imag(wt),real(wt));
%         [~,mIdx] = max(amp'); %get index with max power for each timepiont
%         for i = 1:size(wt,1)
%             lfpphase(i) = phase(i,mIdx(i));
%         end
%         lfpphase = mod(lfpphase,2*pi);
        [wave,f,t,coh,wphases,raw,coi,scale,priod,scalef]=getWavelet(double(lfp.data(:,1)),samplingRate,passband(1),passband(2),8,0);
        [~,mIdx]=max(wave);%get index max power for each timepiont
        pIdx=mIdx'+[0;size(f,2).*cumsum(ones(size(t,1)-1,1))];%converting to indices that will pick off single maxamp index from each of the freq-based phases at eacht timepoint
        lfpphases=wphases(pIdx);%get phase of max amplitude wave at each timepoint
        lfpphases = mod(lfpphases,2*pi);%covert to 0-2pi rather than -pi:pi
% %     case ('peaks')
        % not yet coded
        % filter, smooth, diff = 0, diffdiff = negative
end

%% update intervals to remove sub-threshold power periods
disp('finding intervals below power threshold...')
thresh = mean(power) + std(power)*powerThresh;
minWidth = (samplingRate./passband(2)) * 2; % set the minimum width to two cycles

below=find(power<thresh);
if max(diff(diff(below))) == 0
    below_thresh = [below(1) below(end)];
elseif length(below)>0;
    ends=find(diff(below)~=1);
    ends(end+1)=length(below);
    ends=sort(ends);
    lengths=diff(ends);
    stops=below(ends)./samplingRate;
    starts=lengths./samplingRate;
    starts = [1; starts];
    below_thresh(:,2)=stops;
    below_thresh(:,1)=stops-starts;
else
    below_thresh=[];
end

% now merge interval sets from input and power threshold
intervals = SubtractIntervals(intervals,below_thresh);  % subtract out low power intervals

intervals = intervals(diff(intervals')>minWidth./samplingRate,:); % only keep min width epochs

%% Get phases for each spike for each cell
h = [];
% cum_spkphases = [];
phasebins=[];
spkphases = cell(1,length(spikes.times));
for a = 1:length(spikes.times)

    bools = InIntervals(spikes.times{a},intervals);
    s =spikes.times{a}(bools);
%     s = spikes{a};
    if isempty(s)
        phasedistros(:,a) = zeros(numBins,1);
        phasestats.m(a) = nan;
        phasestats.r(a) = nan;
        phasestats.k(a) = nan;
        phasestats.p(a) = nan;
        phasestats.mode(a) = nan;
        spkphases{a} = nan;
    else
        spkphases{a} = lfpphase(ceil(s*samplingRate));

%         cum_spkphases = vertcat(cum_spkphases, spkphases{a});

    %% Gather binned counts and stats (incl Rayleigh Test)
        [phasedistros(:,a),phasebins,ps]=CircularDistribution(spkphases{a},'nBins',numBins);
        phasestats.m(a) = mod(ps.m,2*pi);
        phasestats.r(a) = ps.r;
        phasestats.k(a) = ps.k;
        phasestats.p(a) = ps.p;
        phasestats.mode(a) = ps.mode;

    %% plotting
        if plotting
            if ~exist('PhaseModulationFig','dir')
                mkdir('PhaseModulationFig');
            end
            h(end+1) = figure;
            hax = subplot(1,2,1);
            rose(spkphases{a})
            title(hax,['Cell #' num2str(a) '. Rayleigh p = ' num2str(phasestats.p(a)) '.'])

            hax = subplot(1,2,2);
            bar(phasebins*180/pi,phasedistros(:,a))
            xlim([0 360])
            set(hax,'XTick',[0 90 180 270 360])
            hold on;
            plot([0:360],cos(pi/180*[0:360])*0.05*max(phasedistros(:,a))+0.95*max(phasedistros(:,a)),'color',[.7 .7 .7])
            set(h(end),'name',['PhaseModPlotsForCell' num2str(a)]);
            print(fullfile('PhaseModulationFig',['PhaseModPlotsForCell' num2str(a)]),'-dpng','-r0');
        end
    end
end
%% Cumulative effect across all spikes from all cells... not saving these stats for now
% phasebins=[];
% if length(cum_spkphases) > 10
%     [cpd,phasebins,cps]=CircularDistribution(cum_spkphases,'nBins',180);
%     cRp = cps.p;
%
%     if plotting
%         h(end+1) = figure;
%         hax = subplot(1,2,1);
%         rose(cum_spkphases)
%         title(hax,['All Spikes/Cells Accumulated. Rayleigh p = ' num2str(cps.p) '.'])
%
%         hax = subplot(1,2,2);
%         bar(phasebins*180/pi,cpd)
%         xlim([0 360])
%         set(hax,'XTick',[0 90 180 270 360])
%         hold on;
%         plot([0:360],cos(pi/180*[0:360])*0.05*max(cpd)+0.95*max(cpd),'color',[.7 .7 .7])
%         set(h(end),'name',['PhaseModPlotsForAllCells']);
%     end
% end

detectorName = 'bz_PhaseModulation';
channels = lfp.channels;
detectorParams = v2struct(intervals,samplingRate,method,plotting,numBins,...
    passband,powerThresh,channels);

PhaseLockingData = v2struct(phasedistros,phasebins,...
                            phasestats,spkphases,...
                            detectorName, detectorParams);
PhaseLockingData.region = spikes.region;
PhaseLockingData.UID = spikes.UID;
PhaseLockingData.sessionName = spikes.sessionName;

if saveMat
 save([lfp.Filename(1:end-4) '.PhaseLockingData.cellinfo.mat'],'PhaseLockingData');
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function gX = rbfardKernDiagGradX(kern, X)

% RBFARDKERNDIAGGRADX Gradient of RBFARD kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the automatic relevance determination radial basis function kernel matrix with
% respect to the elements of the design matrix given in X.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG X : the input data in the form of a design matrix.
% RETURN gX : the gradients of the diagonal with respect to each element
% of X. The returned matrix has the same dimensions as X.
%
% SEEALSO : rbfardKernParamInit, kernDiagGradX, rbfardkernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

gX = zeros(size(X));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function region = p10_region ( )

%*****************************************************************************80
%
%% P10_REGION returns the name of the integration region for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bF = const_bF_appA(PATH)

% number of segments
m = length(PATH)-1;
b = zeros(10*m,1);

% Start
b(1) = PATH(1);
% end
b(end-4) = PATH(end);
% waypoint positions
for i = 1:m-1
    b(10*(i-1)+6) = PATH(i+1);
    b(10*(i-1)+11) = PATH(i+1);
end

% Permutation
C = permut_mat_appA(m);
b_sorted = C*b;
% number of unknowns
unkns = 8*(m-1);
bF = b_sorted(1:end-unkns);

author: yorgoon

repo: minimum-snap-geometric-control


text:
% Usage: e = Faddeeva_erfi(z [, relerr])
%
% Compute erfi(z) = -i*erf(i*z), the imaginary error function,
% for an array or matrix of complex values z.
%
% relerr, if supplied, indicates a desired relative error tolerance in
% w; the default is 0, indicating that machine precision is requested (and
% a relative error < 1e-13 is usually achieved).  Specifying a larger
% relerr may improve performance for some z (at the expense of accuracy).
%
% S. G. Johnson, http://ab-initio.mit.edu/Faddeeva

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function p = legendre(n,var)
%LEGENDRE     Legendre polynomial
%
%   p = legendre(n,var)
%
%Computed by recursion
%  p(0,x) = 1
%  p(1,x) = x
%  p(n,x) = (2*n-1)/n * x * p(n-1,x) - (n-1)/n * p(n-2,x)    for n>1
%
%Parameter var for dependent variable optional (default 'x')
%

% written  07/30/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if nargin==1
    var = 'x';
  else
    if ~ischar(var)
      error('variable must be string')
    end
  end

  if n==0
    p = polynom(1,var);
  elseif n==1
    p = polynom([1 0],var);
  else
    t1 = [1 zeros(1,n)];
    t2 = [0 1 zeros(1,n-1)];
    for i=3:n+1
      t3 = (2*i-3)/(i-1)*[0 t2(1:n)] - (i-2)/(i-1)*t1;
      t1 = t2;
      t2 = t3;
    end
    p = polynom(fliplr(t3),var);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function c=comp_filterbank(f,g,a);
%COMP_FILTERBANK  Compute filtering
%
%   Function groups filters in g according to a presence of .h and .H
%   fields. If .H is present, it is further decided whether it is a full
%   frequency response or a band-limited freq. resp.
%

[L,W]=size(f);
M=numel(g);
c = cell(M,1);

% Divide filters into time domain and frequency domain groups
mFreq = 1:M;
mTime = mFreq(cellfun(@(gEl) isfield(gEl,'h') ,g)>0);
mFreq(mTime) = [];

if ~isempty(mTime)
   % Pick imp. resp.
   gtime = cellfun(@(gEl) gEl.h, g(mTime),'UniformOutput',0);

   % Call the routine
   gskip = cellfun(@(gEl) gEl.offset ,g(mTime));
   c(mTime) = comp_filterbank_td(f,gtime,a(mTime),gskip,'per');
end

if ~isempty(mFreq)
   % Filtering in the frequency domain
   F=fft(f);
   % Pick frequency domain filters
   gfreq = g(mFreq);
   % Divide filters into the full-length and band-limited groups
   mFreqFullL = 1:numel(gfreq);
   amFreqCell = mat2cell(a(mFreq,:).',size(a,2),ones(1,numel(mFreq)));
   mFreqBL = mFreqFullL(cellfun(@(gEl,aEl) numel(gEl.H)~=L || (numel(aEl)>1 && aEl(2) ~=1), gfreq(:),amFreqCell(:))>0);
   mFreqFullL(mFreqBL) = [];

   mFreqFullL = mFreq(mFreqFullL);
   mFreqBL = mFreq(mFreqBL);

   if ~isempty(mFreqFullL)
      G = cellfun(@(gEl) gEl.H, g(mFreqFullL),'UniformOutput',0);
      c(mFreqFullL) = comp_filterbank_fft(F,G,a(mFreqFullL));
   end

   if ~isempty(mFreqBL)
      G = cellfun(@(gEl) gEl.H, g(mFreqBL),'UniformOutput',0);
      foff = cellfun(@(gEl) gEl.foff, g(mFreqBL));
      % Cast from logical to double.
      realonly = cellfun(@(gEl) cast(isfield(gEl,'realonly') && gEl.realonly,'double'), g(mFreqBL));
      c(mFreqBL) = comp_filterbank_fftbl(F,G,foff,a(mFreqBL,:),realonly);
   end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%  rows(M)		equivalent to a=size(M); a=a(end-2);

function [a] = rows(M)

   a = size(M);

   if length(a) == 2,
      a = a(1);
   elseif length(a) > 2,
      a = a(end-2);
   else
      error('woops size wasn''t length 2 or bigger!');
   end;

   return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function retVec=ICRS2J2000F(vec,TT1,TT2)
%%ICRS2J2000F Convert vectors from the International Celestial Reference
%             System (ICRS) to the J2000.0 dynamical frame. This involves
%             removing a frame bias rotation.
%
%INPUTS: vec The 3XN matrix of N vectors that are to be rotated rom the
%            ICRS to the J2000.0 dynamical frame.
%   TT1, TT2 Two parts of a Julian date given in terrestrial time (TT).
%            The units of the date are days. The full date is the sum of
%            both terms. The date is broken into two parts to provide more
%            bits of precision. It does not matter how the date is split.
%
%OUTPUTS: retVec The 3XN set of vectors rotated into the J2000.0 dynamical
%                frame.
%
%As described in [1], the ICRS is ofset from the J2000 dynamical frame by
%a bias rotation.
%
%This is a wrapper for the function iauBp06 and some matrix operations in
%the International Astronomical Union's (IAU) Standard's of Fundamental
%Astronomy library.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%retVec=ICRS2J2000F(vec,TT1,TT2);
%
%REFERENCES:
%[1] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


function isBelow = isVersionBelow(versionA, versionB)
% Checks if versionA is smaller than versionB
    vA         = versionArray(versionA);
    vB         = versionArray(versionB);
    n          = min(length(vA), length(vB));
    deltaAB    = vA(1:n) - vB(1:n);
    difference = find(deltaAB, 1, 'first');
    if isempty(difference)
        isBelow = false; % equal versions
    else
        isBelow = (deltaAB(difference) < 0);
    end
end
% ==============================================================================
function arr = versionArray(str)
% Converts a version string to an array.
    if ischar(str)
        % Translate version string from '2.62.8.1' to [2; 62; 8; 1].
        switch getEnvironment
            case 'MATLAB'
                split = regexp(str, '\.', 'split'); % compatibility MATLAB < R2013a
            case  'Octave'
                split = strsplit(str, '.');
            otherwise
                errorUnknownEnvironment();
        end
        arr = str2num(char(split)); %#ok
    else
        arr = str;
    end
    arr = arr(:)';
end
% ==============================================================================
function errorUnknownEnvironment()
error('matlab2tikz:unknownEnvironment',...
      'Unknown environment "%s". Need MATLAB(R) or Octave.', getEnvironment);
end
% ==============================================================================


text:
function [llh, nx] = tapas_ti_llh_state(y, x, u, theta, ptheta)
%% Generate state space and compute the likelihood.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

nx = tapas_ti_gen_state(y, x, u, theta, ptheta);
llh = tapas_ti_llh(y, nx, u, theta, ptheta);

end

author: translationalneuromodeling

repo: tapas


function ap_multi = compute_AP_multiCam(good_image, junk_image, index, queryCam, testCam)
good_cam = testCam(good_image);
good_cam_uni = unique(good_cam);
ap_multi = zeros(1, 6);

% on the same camera
good_cam_now = queryCam;
ngood = length(junk_image)-1;
junk_image_now = [good_image; index(1)];
good_image_now = setdiff(junk_image, index(1));
old_recall = 0;
old_precision = 1.0;
ap = 0;
intersect_size = 0;
j = 0;
good_now = 0;
for n = 1:length(index)
    flag = 0;
    if ~isempty(find(good_image_now == index(n), 1))
        flag = 1; % good image
        good_now = good_now+1;
    end
    if ~isempty(find(junk_image_now == index(n), 1))
        continue; % junk image
    end

    if flag == 1%good
        intersect_size = intersect_size + 1;
    end
    if ngood == 0
        ap_multi(good_cam_now) = 0;
        break;
    end
    recall = intersect_size/ngood;
    precision = intersect_size/(j + 1);
    ap = ap + (recall - old_recall)*((old_precision+precision)/2);
    old_recall = recall;
    old_precision = precision;
    j = j+1;

    if good_now == ngood
        ap_multi(good_cam_now) = ap;
        break;
    end
end

for k = 1:length(good_cam_uni)
    good_cam_now = good_cam_uni(k);
    ngood = length(find(good_cam == good_cam_now));
    pos_junk = find(good_cam ~= good_cam_now);
    junk_image_now = [junk_image; good_image(pos_junk)];
    pos_good = find(good_cam == good_cam_now);
    good_image_now = good_image(pos_good);
    old_recall = 0;
    old_precision = 1.0;
    ap = 0;
    intersect_size = 0;
    j = 0;
    good_now = 0;
    for n = 1:length(index)
        flag = 0;
        if ~isempty(find(good_image_now == index(n), 1))
            flag = 1; % good image
            good_now = good_now+1;
        end
        if ~isempty(find(junk_image_now == index(n), 1))
            continue; % junk image
        end

        if flag == 1%good
            intersect_size = intersect_size + 1;
        end
        recall = intersect_size/ngood;
        precision = intersect_size/(j + 1);
        ap = ap + (recall - old_recall)*((old_precision+precision)/2);
        old_recall = recall;
        old_precision = precision;
        j = j+1;

        if good_now == ngood
            ap_multi(good_cam_now) = ap;
            break;
        end
    end
end

end


text:
function [Phase_3, Phase_2, Phase_1, UPhase_3] = learn_struct_bnpc(data,node_sizes,epsilon,star)
% G = learn_struct_bnpc(Data,node_sizes,epsilon,star)
%
% Data(i,m) is node i in case m.
% node_sizes and epsilon are optionnals.
% star = 0 to use try_to_separate_B instead of try_to_separate_B_star
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.
%
% Things to do : rewrite function orient_edges !
%                ! sometimes it causes crashes !
%
% V0.91 : 18 sept 2003 (olivier.francois@insa-rouen.fr)

verbose=1;
%if nargin < 5, mwst=0; end
if nargin < 4, star=1; end
if nargin < 3, epsilon=0.05; end
if nargin < 2, node_sizes=max(data'); end

if verbose
  fprintf('================== phase I : \n');
end
tmp1=cputime;
[Phase_1 II JJ score_mat score_mat2] = phaseI(data, node_sizes, epsilon);
tmp1=cputime-tmp1;

if verbose
  fprintf('Execution time : %2.5f\n',tmp1);
  fprintf('\n================== phase II : \n');
end
tmp1=cputime;
Phase_2 = phaseII(Phase_1, data, node_sizes, epsilon, II, JJ, score_mat);
tmp1=cputime-tmp1;

if verbose
  fprintf('Execution time : %2.5f\n',tmp1);
  fprintf('\n================== phase III : \n');
end
tmp1=cputime;
[Phase_3 UPhase_3]  = phaseIII(Phase_2, data, node_sizes, epsilon, score_mat2, star);
tmp1=cputime-tmp1;
if verbose
  fprintf('Execution time : %2.5f\n',tmp1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [G, II, JJ, score_mat, sc2] = phaseI(data,node_sizes,alpha)
% [G, II , JJ, score_mat, s2] = phaseI(data,node_sizes,epsilon)
%
% G is an acyclic graph
% [II JJ] is the list of important edges not processed in phase I (for phase II)
% score_mat is the mutual information score matrix
%
% data(i,m) is node i in case m.
% alpha is the significant level for CI tests ( default=0.05 ).
% node_sizes is the vector of sizes ( default=max(data') ).
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

% 0.
if nargin < 3, alpha=0.05; end
if nargin < 2, node_sizes=max(data'); end
[N m] = size(data);
score_mat = zeros(N);
edges=0;

% 1.
G = zeros(N);
L=[];

% 2. Use of Chi2 instead of MI ... allow using a confidence level alpha instead of an arbitrary epsilon
for i=1:(N-1)
  for j=(i+1):N
    [I score_mat(i,j)] = cond_indep_chisquare(i,j,[],data,'LRT',alpha,node_sizes);
  end
end
sc2=score_mat;

[tmp ordre]=sort(-score_mat(:));
ordre2=ordre(find(-tmp>alpha));
[II JJ]=ind2sub([N N],ordre2);

pointer=1 ;
fini=length(II);

% 3.
edges=2;
for pointer=1:min(2,fini),
  %fprintf('%d-%d\n',II(pointer),JJ(pointer));
  G(II(pointer),JJ(pointer))=1;
  G(JJ(pointer),II(pointer))=1;
  score_mat(II(pointer),JJ(pointer))=-inf;
end

pointer=min(2,fini);
arret=0;

while pointer<fini & ~arret
  % 4.
  pointer=pointer+1;
  C = ~reachability_graph(G);
  if C(II(pointer),JJ(pointer))
    %fprintf('%d-%d\n',II(pointer),JJ(pointer));
    G(II(pointer),JJ(pointer))=1;
    G(JJ(pointer),II(pointer))=1;
    score_mat(II(pointer),JJ(pointer))=-inf;
    edges=edges+1;
    if edges==N-1
      arret=1;
    end
  end

  % 5.
end

[tmp ordre]=sort(-score_mat(:));
ordre2=ordre(find(-tmp>alpha));
[II JJ]=ind2sub([N N],ordre2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G = phaseII(G1,data,node_sizes,alpha,II, JJ, score_mat)
% G = phaseII(G1,data,node_sizes,epsilon,II,JJ, score_mat)
%
% G1, II, JJ,score_mat are given by phaseI.
% data(i,m) is node i in case m.
% node_sizes is the vector of sizes ( default=max(data') ).
% alpha is the significant level for CI tests ( default=0.05 ).
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

st{1}='added';
st{2}='';
% 0.
[N m] = size(data);

G=G1;

% 6.
II=II(end:-1:1);
JJ=JJ(end:-1:1);

pointer=length(II);

while pointer>0
  % 7.
  trysep = try_to_separate_A(G,II(pointer),JJ(pointer),data,alpha,node_sizes);
  %fprintf('Try to separate %d and %d : %s\n',II(pointer),JJ(pointer),st{trysep+1});
  if ~trysep
    G(II(pointer),JJ(pointer))=1;
    G(JJ(pointer),II(pointer))=1;
    %fprintf('%d-%d\n',II(pointer),JJ(pointer));
  end
  % 8.
  pointer=pointer-1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [G,U] = phaseIII(G1,data,node_sizes,alpha,s2,star)
% [G] = phaseIII(G1,data,node_sizes,alpha,s2,star)
%
% data(i,m) is node i in case m.
% if star~=0, use try_to_separate_B_star instead of try_to_separate_B (default star=1).
% G is the non-oriented graph and G1 the oriented result.
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

% 0.
if nargin < 2, disp('Not enough arguments');return; end
if nargin < 3, node_sizes=max(data'); end
if nargin < 4, alpha = 0.05; end
if nargin < 5, star=1; end
G=G1;
N=length(G);
% reachability_matrix of G
M = expm(full(G)) - eye(length(G)); M = (M>0);

% 9.
fprintf('Thinning - separateA\n');

% Edges are examined in the inverse order of their Chi2 (or MI) score
s2(find(~G))=0;
[tmp ordre]=sort(s2(:));
ordre2=ordre(find(tmp>0));
[I J]=ind2sub([N N],ordre2);
ii=1:length(I);
for i=ii,
  %fprintf('%d-%d : ',I(i),J(i));
  G(I(i),J(i))=0;
  G(J(i),I(i))=0;
  trysep = try_to_separate_A(G,I(i),J(i),data,alpha,node_sizes);

  if ~trysep,
    G(I(i),J(i))=1;
    G(J(i),I(i))=1;
    %fprintf(' keep\n');
    %else
    %fprintf('delete\n');
  end
end

% 10.
fprintf('Thinning - separateB'); if star; fprintf('star'); end; fprintf('\n');
s2(find(~G))=0;
[tmp ordre]=sort(s2(:));
ordre2=ordre(find(tmp>0));
[I J]=ind2sub([N N],ordre2);
ii=1:length(I);
for i=ii,
  %fprintf('%d-%d : ',I(i),J(i));
  G(I(i),J(i))=0; G(J(i),I(i))=0;
  if star==0
    trysep = try_to_separate_B(G,I(i),J(i),data,node_sizes,alpha);
  else
    trysep = try_to_separate_B_star(G,I(i),J(i),data,node_sizes,alpha);
  end
  if ~trysep,
    G(I(i),J(i))=1; G(J(i),I(i))=1;
    %fprintf(' keep\n');
    %else
    %fprintf('delete\n');
  end
end

%11.
fprintf('Thinning - orient_edges\n');
U=G;
G=orient_edges(U,data,node_sizes,alpha);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [I, N1, N2] = try_to_separate_A(G,node1,node2,data,alpha,node_sizes)
% [I, N1, N2] = try_to_separate_A(G,node1,node2,data,alpha,node_sizess)
%
% G is the current partially directed graph.
% I is a boolean : I=1 <==> separated.
% N1 : neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2).
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

% 0.
if node1==node2 | length(G)<2
  disp('Error: Check your arguments in try_to_separate_A.');I=0; return
end
N=size(data,1);
if nargin==4, alpha=0.05; node_sizes=max(data'); end
if nargin==5, node_sizes=max(data'); end

% 1.
N1=find(G(node1,:)==1);N01=N1;
GG1=G(setdiff(1:N,node1),setdiff(1:N,node1));
node22=node2-(node2>node1);
% reachability_matrix of GG1
M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
% N1 is the neighbors of node1 that are on the adjacency between node1 and node2
for i=N1
  j=i-(i>node1);
  if M(j,node22)~=1, N01=setdiff(N01,i); N1=setdiff(N1,i); end
  % 2.
  if ~G(node1,i), N1=setdiff(N1,i); end
end

N2=find(G(node2,:)==1);N02=N2;
GG2=G(setdiff(1:N,node2),setdiff(1:N,node2));
node12=node1-(node2<node1);
% reachability_matrix of GG2
M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
% N2 is the neighbors of node2 that are on the adjacency between node1 and node2
for i=N2
  j=i-(i>node2);
  if M(node12,j)~=1, N02=setdiff(N02,i); N2=setdiff(N2,i); end
  % 2.
  if ~G(node2,i), N2=setdiff(N2,i); end
end

%fprintf('%d : N1=',node1); fprintf('%d',N1); fprintf('\n');
%fprintf('%d : N2=',node2); fprintf('%d',N2); fprintf('\n');
% 3.
if length(N1)>length(N2), tmp=N1; N1=N2; N2=tmp; clear tmp, end
% 4.
C=N1;
for test=1:2
  if test==2, C=N2; end
  % 5.
  [I v1] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes);
  if I, %fprintf('%d-%d separated (%2.5f) by C=',node1,node2,v1); fprintf('%d',C); fprintf('\n');
    return,
    %else
    %fprintf('%d-%d not separated (%2.5f) by C=',node1,node2,v1); fprintf('%d',C); fprintf('\n');
  end;

  % 6.
  step6=1;
  while step6
    step6=0;
    if length(C)>=1
      v=[];
      for i=C
	Ci = setdiff(C,i);
	[I(i) v(i)] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
      end
      [vm ind] = min(v);

      % 7.
      if I(ind)
	I = 1; return
      else
	if vm < v1
	  v1 = vm;
	  C = setdiff(C,ind);
	  % goto step 6.
	  step6 = 1;
	end
      end
    end

    % 8.
    if test==2, I=0; return, end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function I = try_to_separate_B(G,node1,node2,data,node_sizes,alpha,N1,N2)
% I = try_to_separate_B(G,node1,node2,data,node_sizes,epsilon,N1,N2)
%
% G is the current partially directed graph.
% data(i,m), node i in case m.
% node_sizes is the vector of size of the attributs in data ( default max(data') )
% N1 (optionnal) is the neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2).
% I is a boolean : I+1 <==> separated.
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

% 0.
if node1==node2 | length(G)<2
  disp('Error: Verify your arguments in try_to_separate_B.');I=0; return
end
if nargin < 4,
  disp('Error : not enougth arguments'); I=0; return;
end
N=size(data,1);

% 1.
if nargin < 8
  N1=find(G(node1,:)==1);
  GG1=G(setdiff(1:N,node1),setdiff(1:N,node1));
  node22=node2-(node2>node1);
  M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
  for i=N1
    j=i-(i>node1);
    if M(j,node22)~=1, N1=setdiff(N1,i); end
  end
  N2=find(G(node2,:)==1);
  GG2=G(setdiff(1:N,node2),setdiff(1:N,node2));
  node12=node1-(node2<node1);
  M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
  for i=N2
    j=i-(i>node2);
    if M(node12,j)~=1, N2=setdiff(N2,i); end
  end
end
if nargin < 6, alpha=0.05; end
if nargin < 5, node_sizes=max(data'); end
M = expm(full(G)) - eye(length(G)); M = (M>0);

% 2.
N1b=[];
for i=N1
  NN1=find(G(i,:)==1);
  for j=NN1
    if M(i,j)~=1 & ~ismember(j,N1), N1b=union(N1b,j); end
  end
end

% 3.
N2b=[];
for i=N2
  NN2=find(G(i,:)==1);
  for j=NN2
    if M(i,j)~=1 & ~ismember(j,N2), N2b=union(N2b,j); end
  end
end

% 4.
if length(union(N1,N1b)) < length(union(N2,N2b))
  C=union(N1,N1b);
else
  C=union(N2,N2b);
end

% 5.
continu=1;
while continu
  l=length(C);
  %fprintf('%d',continu);
  [I v] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes);
  if I==1; return, elseif l<2, I=0; return, end

  % 6.
  Cb=C;
  for i=1:l
    Ci=setdiff(C,C(i));
    [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
    e = (v+1)/3;    % e is a small value...
	if I==1,return, elseif vi<v+e, Cb=setdiff(Cb,C(i)); end
  end

  % 7.
  if length(Cb) < l, C=Cb; else continu==0; I=0; return; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function I = try_to_separate_B_star(G,node1,node2,data,node_sizes,alpha,N1,N2)
% I = try_to_separate_B_star(G,node1,node2,data,epsilon,node_sizess)
%
% G is the current partially directed graph.
% I is a boolean.
% N1 is the neighbors of node1 that are on an adjacency path between node1 and node2 (ditto for N2).
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

% 0.
if node1==node2 | length(G)<2
  disp('Error: Verify your arguments in try_to_separate_B_star.');I=0; return
end
if nargin < 4,
  disp('Error : not enougth arguments');I=0; return;
end
N=size(data,1);

% 1.
if nargin < 8
  N1=find(G(node1,:)==1);
  GG1=G(setdiff(1:N,node1),setdiff(1:N,node1));
  node22=node2-(node2>node1);
  M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
  for i=N1
    j=i-(i>node1);
    if M(j,node22)~=1, N1=setdiff(N1,i); end
    % 2.
    if ~G(node1,i), N1=setdiff(N1,i); end
  end
  N2=find(G(node2,:)==1);
  GG2=G(setdiff(1:N,node2),setdiff(1:N,node2));
  node12=node1-(node2<node1);
  M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
  for i=N2
    j=i-(i>node2);
    if M(node12,j)~=1, N2=setdiff(N2,i); end
    % 2.
    if ~G(node2,i), N2=setdiff(N2,i); end
  end
end
if nargin < 6, alpha=0.05; end
if nargin < 5, node_sizes=max(data'); end

% 3.
if length(N1)>length(N2), tmp=N1; N1=N2; N2=tmp; end

% 4.
C=N1;
l=length(C);
I=0;

% 5.
for test=1:2
  continu=1;
  %test
  if test==2 & ~isempty(N2)
    C=N2; l=length(C); IsInCi=zeros(1,l); IsInCi(l)=1;
  else IsInCi=zeros(1,l);
  end
  s=ones(1,l);
  % Pour tous les sous-ensemble Ci de C :
  while continu & ~isempty(C)
    Ci = setdiff(C.*IsInCi,0);
    [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
    if I, %fprintf('%d-%d separated (%2.5f) by C=',node1,node2,vi); fprintf('%d',Ci); fprintf('\n');
      return,
      %else
      %fprintf('%d-%d not separated (%2.5f) by C=',node1,node2,vi); fprintf('%d',Ci); fprintf('\n');
    end;
    %    if I, return, end

    if IsInCi==s, continu=0;
    else
      IsInCi(l)=IsInCi(l)+1;

      notOK=1; i=l;
      while notOK & i>1
	if IsInCi(i)>s(i), IsInCi(i)=0; IsInCi(i-1)=IsInCi(i-1)+1; else notOK=0; end
	i=i-1;
      end
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G1 = orient_edges(G,data,node_sizes,alpha)
% G1 = orient_edges(G,data,node_sizes)
%
% G is the partially directed graph.
% data(i,m), node i in case m.
% node_sizes is the vector of size of the attributs in data ( default max(data') )
%
% see "Learning bayesian Networks from Data: A Efficient Approach Based on Information Theorie"
%     Jie Cheng, David Bell and Weird Liu.

% 0.
if nargin < 4, alpha=0.05; end
if nargin < 3, node_sizes=max(data'); end
if nargin < 2, disp(' Require at least two arguments.'); return, end
N=length(G);
G1=G;

% 1.
[Lnode1 Lnode2]=find(triu(1-triu(G),1)); %[Lnode1 Lnode2]=ind2sub([N N],find(triu(1-triu(G),1)));
for ii = 1:length(Lnode1),
  node1=Lnode1(ii);
  node2=Lnode2(ii);
  %fprintf('%d %d\n',node1,node2);
  N1 = find(G(node1,:)==1);
  N2 = find(G(node2,:)==1);
  if ~isempty(intersect(N1,N2))
    %fprintf('V1='); fprintf('%d ',N1); fprintf('\n');
    %fprintf('V2='); fprintf('%d ',N2); fprintf('\n');
    GG1 = G(setdiff(1:N,node1),setdiff(1:N,node1));
    node22 = node2-(node2>node1);
    % reachability_matrix of GG1
    M = expm(full(GG1)) - eye(length(GG1)); M = (M>0);
    % N1 is the neighbors of node1 that are on the adjacency between node1 and node2
    for i = N1
      j = i-(i>node1);
      if M(j,node22)~=1, N1=setdiff(N1,i);
      end
    end
    GG2 = G(setdiff(1:N,node2),setdiff(1:N,node2));
    node12 = node1-(node2<node1);
    % reachability_matrix of GG2
    M = expm(full(GG2)) - eye(length(GG2)); M = (M>0);
    % N2 is the neighbors of node2 that are on the adjacency between node1 and node2
    for i=N2
      j = i-(i>node2);
      if M(node12,j)~=1, N2=setdiff(N2,i); end
    end
    %fprintf('%d %d\n',node1,node2);
    %fprintf('N1='); fprintf('%d ',N1); fprintf('\n');
    %fprintf('N2='); fprintf('%d ',N2); fprintf('\n');

    % 2.
    M = expm(full(G)) - eye(length(G)); M = (M>0);
    N1b=N1;
    for i=N1
      NN1 = find(G(i,:)==1);
      for j=NN1
	if M(i,j)~=1 & ~ismember(j,N1), N1b=union(N1b,j); end
      end
    end
    %fprintf('N1''='); fprintf('%d ',N1b); fprintf('\n');

    % 3.
    N2b=N2;
    for i=N2
      NN2 = find(G(i,:)==1);
      for j=NN2
	if M(i,j)~=1 & ~ismember(j,N2), N2b=union(N2b,j); end
      end
    end
    %fprintf('N2''='); fprintf('%d ',N2b); fprintf('\n');

    % 4.
    if length(N1b) < length(N2b)
      C = N1b;
    else
      C = N2b;
    end
    %l=length(C);
    %fprintf('C='); fprintf('%d ',C); fprintf('\n');

    % 7.
    step5=1;
    while step5
      step5=0;
      %fprintf('.');
      % 5.
      l=length(C);
      [I v] = cond_indep_chisquare(node1,node2,C,data,'LRT',alpha,node_sizes);
      %fprintf('C='); fprintf('%d ',C);
      %fprintf(': %d %2.5f\n',I,v);

      step8=0;
      if I==1 & v~=0 % v < epsilon
	step8=1;
      else
	if  l==1
	  G1(C,node1)=0; G1(C,node2)=0;
	  fprintf('%d -> %d <- %d\n',node1,C,node2);
	  step8=1;
	end
      end
      %fprintf('%d\n',step8);

      % 6.
      if ~step8
	Cb=C;
	for i=1:l
	  Ci=setdiff(C,C(i));
	  [I vi] = cond_indep_chisquare(node1,node2,Ci,data,'LRT',alpha,node_sizes);
	  % e = (v+1)/3;    % e is a small value...
	  if I==1 % vi < v+e
	    Cb=setdiff(Cb,C(i));
	    if ismember(C(i),N1) & ismember(C(i),N2)
	      G1(C(i),node1)=0; G1(C(i),node2)=0;
	      fprintf('%d -> %d <- %d\n',node1,C(i),node2);
	    end
	    if I==1 % vi < epsilon
	      step8=1;
	    end
	  end
	end % for
      end % if

      % 7.
      if ~step8
	if length(Cb) < length(C), C=Cb; end
	if length(C) > 0, step5==1; end
      end
    end % while step5
    % step8 : passer � la paire de noeud suivant
    %fprintf('\n');
    %else
    %fprintf(' No common neighbor\n');
  end % if
end % for

% 11.
step9=0;
fprintf('Infering directions ');
test = pdag_to_dag(G1);
while ~isdag(test) %& step9<N
  step9=step9+1;
  %if ~isempty(test)
    %fprintf('.');
    % 9.
    for a=1:N, for b=1:N, for c=1:N,
	  if a~=b & b~=c & c~=a
	    if G1(a,b)==1 & G1(b,a)==0
	      %fprintf('%d -> %d ... \n',a,b);
	      if G1(b,c)==1 & G1(c,b)==1
		if G1(a,c)+G1(c,a)==0,
		  G1(c,b)=0;
		  fprintf('%d -> %d (9)\n',b,c);
		end
	      end
	    end
	  end
    end, end, end

    % 10.
    for a=1:N-1, for b=a+1:N
	if G1(a,b)==1 & G1(b,a)==1
	  GGG1=xor(G1,G1'); % matrice des arcs orient�s de G1
	  M = expm(full(GGG1)) - eye(length(GGG1)); M = (M>0);
	  if M(a,b)==1,
	    G1(b,a)=0;
	    fprintf('%d -> %d (10)\n',a,b);
	  end
	end
    end, end

      test = pdag_to_dag(G1);
      if isempty(test), G1=return_one_edge(G1); end
  %else
  %  G1, return,
  %end
end % while 11.
G1 = pdag_to_dag(G1);
fprintf('%d boucles\n',step9);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function b = isdag(G)
b = sum(sum(G.*G'));        % How many undirected arcs ? (x2)
b=~b & ~isempty(G);
if b
  M = expm(full(G)) - eye(length(G)); M = (M>0);
  b = b & find(sum(sum(eye(length(G)).*M))); % is there no cycle ?
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G=return_one_edge(G)
N=length(G); fam=[]; node2=[];
permnode = randperm(5); i=0; %fprintf('rev\n');
while isempty(fam) | isempty(node2) | i==N
 i=i+1;
 node = permnode(i); %node = ceil(rand(1)*N);
 fam=find(G(node,:)==1);
 par=find(G(:,node)==1);
 fam = myunion(fam, par);
 if ~isempty(fam)
  node2 = fam(ceil(rand(1)*length(fam)));
  par2=find(G(:,node)==1);
  if ~isempty(intersect(par2, node)), node2=[]; end
 end
end
if isempty(myintersect(node, par2)),
 G(node, node2)=0;
 G(node2, node)=1;
 fprintf('%d -> %d (Rev)\n',node, node2);
else
 G(node, node2)=1;
 G(node2, node)=0;
 fprintf('%d -> %d (Rev)\n',node2, node);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function pass = test_splitting(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% Store current globlal state.
p = chebfunpref();
chebfunpref.setDefaults('factory')

warnState = warning('off', 'CHEBFUN:splitting:deprecated');

try
    % Testing takes place inside this TRY CATCH.

    pass(1) = strcmp(splitting(), 'off');

    pass(2) = strcmp(splitting('on'), 'off');

    pass(3) = strcmp(splitting(), 'on');

    F = @(x) cos(pi*x).*sign(x);
    f = chebfun(@(x) F(x), [-1 1]);
    x = linspace(-1, 1, 20);
    err = norm(f(x) - F(x), inf);
    pass(4) = numel(f.domain > 1) && err < 1e-5;

    pass(5) = strcmp(splitting('off'), 'on');

    pass(6) = strcmp(splitting(), 'off');

catch ME
    % Reset preferences and warning state
    chebfunpref.setDefaults(p);
    warning(warnState);

    % Rethrow error:
    rethrow(ME)
end

% Return to default settings:
chebfunpref.setDefaults(p);
warning(warnState);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function w=antelope2waveform_test()
    %dbpath = '/aerun/sum/db/archive/archive_2011/archive_2011_11_16';
    gismopath = fileparts(which('startup_GISMO'));
    dbpath = fullfile(gismopath, 'tests', 'test_data', 'antelope2waveform_testdb')
    sta = '.*';
    sta = 'AKT';
    sta = '.*';
    %sta = 'AKT|CRAP|VIB';
    chan = 'HHZ.*';
    %chan = 'HH.*';
    starttime = datenum2epoch(datenum('16-Nov-2011 16:29:00'));
    endtime = datenum2epoch(datenum('16-Nov-2011 16:43:24'));
%     dbexample_get_demodb_path;
%     dbpath = demodb_path;
%     sta = '.*';
%     chan = '.*';
%     starttime = [];
%     endtime = [];
     w=antelope.antelope2waveform(dbpath, sta, chan, starttime, endtime);
end

author: geoscience-community-codes

repo: GISMO


text:
function lattice_rule_test01 ( )

%*****************************************************************************80
%
%% LATTICE_RULE_TEST01 tests FIBONACCI_LATTICE_Q;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994, page 78-80, page 145.
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATTICE_RULE_TEST01\n' );
  fprintf ( 1, '  FIBONACCI_LATTICE_Q applies a Fibonacci lattice rule\n' );
  fprintf ( 1, '  to integrate a function over the unit square.\n' );
  fprintf ( 1, '  These Fibonacci rules are only available in 2D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM = %d\n', dim_num );

  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  exact = e_01_2d ( dim_num, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '         K         M    EXACT       ESTIMATE    ERROR\n' );
  fprintf ( 1, '\n' );

  for k = 3 : 18

    quad = fibonacci_lattice_q ( k, @f_01_2d );

    error = abs ( exact - quad );
    m = fibonacci ( k );

    fprintf ( 1, '  %8d  %8d  %10.6f  %10.6f  %10.6e\n', k, m, exact, quad, error );

  end

  return
end


text:
function recastMrImage = recast_as_MrImage(this)
% recasts an MrImageSpm4D as a 4D MrImage, changes available methods etc.
%
%   Y = MrImageSpm4D()
%   YasMrImage = Y.recast_as_MrImage()
%
% This is a method of class MrImageSpm4D.
%
% IN
%
% OUT
%
% EXAMPLE
%   recast_as_MrImage
%
%   See also MrImageSpm4D

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-05-03
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

recastMrImage = MrImage();
recastMrImage.update_properties_from(this);

% house keeping: rename, if default name was used before, add info
% about recast
if strcmp(recastMrImage.name, 'MrImageSpm4D')
    recastMrImage.name = 'MrImage';
end
recastMrImage.info{end+1,1} = 'recast_as_MrImage';


text:
function a=v_polygonarea(p)
%V_POLYGONAREA Calculate the area of a polygon
%
% Inputs:
%    P(n,2) is the polygon vertices
%
% Outputs:
%    A is teh area of teh polygon
%
% The area is positive if the vertices go anti-clockwise around the polygon.
%

%      Copyright (C) Mike Brookes 2009
%      Version: $Id: v_polygonarea.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p(end+1,:)=p(1,:);      % append an extra point
a=0.5*sum((p(1:end-1,1)-p(2:end,1)).*(p(1:end-1,2)+p(2:end,2)),1);
if ~nargout
    plot(p(:,1),p(:,2),'b-');
    mnx=[1.05 -0.05;-0.05 1.05]*[min(p);max(p)];
    set(gca,'xlim',mnx(:,1)','ylim',mnx(:,2)');
    title(sprintf('Area = %.2g',a));
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [map,timeBins] = SyncMap(synchronized,indices,varargin)

%SyncMap - Create a map from successive event-synchronized data.
%
% Using N repetitions of similar data centered around synchronizing
% events (e.g. N evoked potentials), create a map v = f(t,i) where
% t is the time relative to the synchronizing events and i is the
% occurrence (from 1 to N).
%
%  USAGE
%
%    [map,timeBins] = SyncMap(synchronized,indices,<options>)
%
%    synchronized   event-synchronized samples
%    indices        list of synchronizing event indices for each sample
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'durations'   durations before and after synchronizing events for each
%                   trial (in s) (default = [-0.5 0.5])
%     'nBins'       total number of time bins (default 100)
%     'smooth'      smoothing size (0 = no smoothing) (default = 0.01*nBins)
%    =========================================================================
%
%  SEE
%
%    See also Sync, SyncHist, PlotSync, PETHTransition.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Default values
durations = [-0.5 0.5];
nBins = 100;
smooth = [];

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help SyncMap'' for details).');
end

% Check parameter sizes
if size(synchronized,2) > 2,
	error('Parameter ''synchronized'' is not a Nx2 matrix (type ''help SyncMap'' for details).');
end
if size(indices,2) ~= 1,
	error('Parameter ''indices'' is not a vector (type ''help SyncMap'' for details).');
end
if size(indices,1) ~= size(synchronized,1),
	error('Parameters ''synchronized'' and ''indices'' have different lengths (type ''help SyncMap'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help SyncMap'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'durations',
			durations = varargin{i+1};
			if ~isdvector(durations,'#2','<'),
				error('Incorrect value for property ''durations'' (type ''help SyncMap'' for details).');
			end
		case 'nbins',
			nBins = varargin{i+1};
			if ~isiscalar(nBins,'>0'),
				error('Incorrect value for property ''nBins'' (type ''help SyncMap'' for details).');
			end
		case 'smooth',
			smooth = varargin{i+1};
			if ~isdvector(smooth,'>=0'),
				error('Incorrect value for property ''smooth'' (type ''help SyncMap'' for details).');
			end

		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SyncMap'' for details).']);
	end
end

if isempty(smooth),
	smooth = ceil(0.01*nBins);
end

start = durations(1);
stop = durations(2);
timeBinSize = (stop - start)/nBins;
timeBins = (start:timeBinSize:stop-timeBinSize)+timeBinSize/2;
binnedTime = Bin(synchronized(:,1),durations,nBins);
if size(synchronized,2) == 1,
	% Occurrences
	s = Accumulate([indices binnedTime],1);
	map = Smooth(s,smooth);
else
	% Values
	s = Accumulate([indices binnedTime],synchronized(:,2));
	n = Accumulate([indices binnedTime],1);
	n(n==0) = 1;
	map = Smooth(s,[smooth 0])./Smooth(n,[smooth 0]);
end


text:
function [r, c, hnew] = houghpeaks(h, numpeaks, threshold, nhood)
%HOUGHPEAKS Detect peaks in Hough transform.
%   [R, C, HNEW] = HOUGHPEAKS(H, NUMPEAKS, THRESHOLD, NHOOD) detects
%   peaks in the Hough transform matrix H.  NUMPEAKS specifies the
%   maximum number of peak locations to look for.  Values of H below
%   THRESHOLD will not be considered to be peaks.  NHOOD is a
%   two-element vector specifying the size of the suppression
%   neighborhood.  This is the neighborhood around each peak that is
%   set to zero after the peak is identified.  The elements of NHOOD
%   must be positive, odd integers.  R and C are the row and column
%   coordinates of the identified peaks.  HNEW is the Hough transform
%   with peak neighborhood suppressed.
%
%   If NHOOD is omitted, it defaults to the smallest odd values >=
%   size(H)/50.  If THRESHOLD is omitted, it defaults to
%   0.5*max(H(:)).  If NUMPEAKS is omitted, it defaults to 1.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 13:34:50 $

if nargin < 4
   nhood = size(h)/50;
   % Make sure the neighborhood size is odd.
   nhood = max(2*ceil(nhood/2) + 1, 1);
end
if nargin < 3
   threshold = 0.5 * max(h(:));
end
if nargin < 2
   numpeaks = 1;
end

done = false;
hnew = h; r = []; c = [];
while ~done
   [p, q] = find(hnew == max(hnew(:)));
   p = p(1); q = q(1);
   if hnew(p, q) >= threshold
      r(end + 1) = p; c(end + 1) = q;

      % Suppress this maximum and its close neighbors.
      p1 = p - (nhood(1) - 1)/2; p2 = p + (nhood(1) - 1)/2;
      q1 = q - (nhood(2) - 1)/2; q2 = q + (nhood(2) - 1)/2;
      [pp, qq] = ndgrid(p1:p2,q1:q2);
      pp = pp(:); qq = qq(:);

      % Throw away neighbor coordinates that are out of bounds in
      % the rho direction.
      badrho = find((pp < 1) | (pp > size(h, 1)));
      pp(badrho) = []; qq(badrho) = [];

      % For coordinates that are out of bounds in the theta
      % direction, we want to consider that H is antisymmetric
      % along the rho axis for theta = +/- 90 degrees.
      theta_too_low = find(qq < 1);
      qq(theta_too_low) = size(h, 2) + qq(theta_too_low);
      pp(theta_too_low) = size(h, 1) - pp(theta_too_low) + 1;
      theta_too_high = find(qq > size(h, 2));
      qq(theta_too_high) = qq(theta_too_high) - size(h, 2);
      pp(theta_too_high) = size(h, 1) - pp(theta_too_high) + 1;

      % Convert to linear indices to zero out all the values.
      hnew(sub2ind(size(hnew), pp, qq)) = 0;

      done = length(r) == numpeaks;
   else
      done = true;
   end
end

repo: weiyanmin


text:
function hpp_test03 ( )

%*****************************************************************************80
%
%% HPP_TEST03 tests HEPP_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HPP_TEST03:\n' );
  fprintf ( 1, '  HePP_VALUE evaluates a Hermite product polynomial.\n' );
  fprintf ( 1, '  POLYNOMIAL_VALUE evaluates a polynomial.\n' );

  m = 3;
  n = 1;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( m, -1.0, +1.0, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Evaluate at X = ( %g, %g, %g )\n', x(1:3) )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rank  I1  I2  I3:  He(I1,X1)*He(I2,X2)*He(I3,X3)  P(X1,X2,X3)\n' );
  fprintf ( 1, '\n' );

  for rank = 1 : 20

    l = comp_unrank_grlex ( m, rank );
%
%  Evaluate the HePP directly.
%
    v1 = hepp_value ( m, n, l, x );
%
%  Convert the HePP to a polynomial, and reevaluate.
%
    o_max = prod ( floor ( ( l(1:m) + 2 ) / 2 ) );

    [ o, c, e ] = hepp_to_polynomial ( m, l, o_max );

    v2 = polynomial_value ( m, o, c, e, n, x );
%
%  Compare results.
%
    fprintf ( 1, '  %4d  %2d  %2d  %2d  %12g  %12g\n', rank, l(1:m), v1, v2 );

  end

  return
end


text:
function sVF = grad(sF, varargin)
% calculates the gradient of a spherical harmonic
%
% Syntax
%   sVF = grad(sF) % returns the gradient as a spherical vector field
%   g = grad(sF, v) % return the gradient in point v as vector3d
%
% Input
%  sF - @S2FunHarmonic
%  v - @vector3d
%
% Output
%  sVF - @sphericalVectorFieldHarmonic
%    g - @vector3d
%

if nargin > 1
  sF = [sF.drho; sF.dthetasin];
  v = varargin{1};
  y = eval(sF, v);
  sVF = ...
    y(:, 1)./sin(v.theta).^2.*S2VectorField.rho(v)+ ...
    y(:, 2) .* S2VectorField.theta(v);

  sVF(isnan(sVF)) = vector3d([0 0 0]);

else
  sF = [sF.drho; sF.dtheta];
  sVF = S2VectorFieldHarmonic(sF);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [network, mappedX, reconX] = train_autoencoder(X, layers, noise, max_iter)
%TRAIN_AUTOENCODER Trains an simple autoencoder
%
%   [network, mappedX, reconX] = train_encoder(X, layers, noise, max_iter)
%
% Trains up an autoencoder with the structure that is specified in layers.
% The low-dimensional data is returned in mappedX, and the network in
% network. The variable noise indicates how much of the input data should be
% blacked out (default = 0).
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if nargin < 2
        error('Not enough inputs.');
    end
    if isempty(layers)
        error('There should be at least one hidden layer.');
    end
    if ~exist('noise', 'var') || isempty(noise)
        noise = 0;
    end
    if ~exist('max_iter', 'var') || isempty(max_iter)
        max_iter = 50;
    end

    % Initialize the network
    D = size(X, 2);
    no_layers = length(layers) + 1;
    network = cell(no_layers, 1);
    network{1}.W = randn(D, layers(1)) * .0001;
    network{1}.bias_upW = zeros(1, layers(1));
    for i=2:no_layers - 1
        network{i}.W = randn(layers(i - 1), layers(i)) * .0001;
        network{i}.bias_upW = zeros(1, layers(i));
    end
    network{no_layers}.W = randn(layers(end), D) * .0001;
    network{no_layers}.bias_upW = zeros(1, D);
    reconX = run_data_through_autoenc(network, X);
    disp(['Initial MSE of reconstructions: ' num2str(mean((X(:) - reconX(:)) .^ 2))]);

    % Perform backpropagation to minimize reconstruction error
    network = backprop(network, X, X, max_iter, noise);

    % Get representation from hidden layer
    [reconX, mappedX] = run_data_through_autoenc(network, X);
    disp(['Final MSE of reconstructions: ' num2str(mean((X(:) - reconX(:)) .^ 2))]);

repo: cluster


text:
function msm_to_st ( a, output_filename )

%*****************************************************************************80
%
%% MSM_TO_ST writes a MATLAB Sparse Matrix to an ST Sparse Triplet file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real sparse A(), the MATLAB Sparse Matrix.
%
%    Input, string OUTPUT_FILENAME, the name of the ST file.
%
  [ nrow, ncol ] = size ( a );
  nnzero = nnz ( a );

  [ rows, cols, vals ] = find ( a );

  base0 = 0;
  base1 = 1;
  rows = st_rebase ( base1, base0, nnzero, rows );
  cols = st_rebase ( base1, base0, nnzero, cols );

  st_header_write ( output_filename, nrow, ncol, nnzero );

  st_data_write ( output_filename, nrow, ncol, nnzero, rows, cols, vals );

  return
end


text:
function cspy(S,epsilon)
  % cspy Visualize sparsity pattern. coloring positive and negative entries of S
  % in red and blue respectively
  %
  % cspy(S)
  % cspy(S,epsilon)
  %
  % Inputs:
  %   S  m by n (sparse) matrix
  %   epsilon  zero value {0}
  %
  % See also: spy
  %

  if nargin < 2
    epsilon = 0;
  end
  spy(S>epsilon,'r');
  hold on;
  spy(S<-epsilon,'b');
  hold off;
end

author: alecjacobson

repo: gptoolbox


text:
function [joints, parts] = util_get_parts24()

parts = [];
joints = [];

parts.id    = 0;   parts.pos    = [0 0];     parts.xaxis    = [1 0];
parts(2).id = 1;   parts(2).pos = [0 1];     parts(2).xaxis = [1 0];
parts(3).id = 2;   parts(3).pos = [1 1];     parts(3).xaxis = [1 0];
parts(4).id = 3;   parts(4).pos = [1 2];     parts(4).xaxis = [2 1];
parts(5).id = 4;   parts(5).pos = [2 2];     parts(5).xaxis = [2 1];
parts(6).id = 5;   parts(6).pos = [3 3];     parts(6).xaxis = [3 4];
parts(7).id = 6;   parts(7).pos = [3 4];     parts(7).xaxis = [3 4];
parts(8).id = 7;   parts(8).pos = [4 4];     parts(8).xaxis = [4 5];
parts(9).id = 8;   parts(9).pos = [4 5];     parts(9).xaxis = [4 5];
parts(10).id = 9;  parts(10).pos = [5 5];    parts(10).xaxis = [4 5];
parts(11).id = 10; parts(11).pos = [6 7];    parts(11).xaxis = [7 6];
parts(12).id = 11; parts(12).pos = [8 9];    parts(12).xaxis = [9 8];
parts(13).id = 12; parts(13).pos = [10 10];  parts(13).xaxis = [11 10];
parts(14).id = 13; parts(14).pos = [11 10];  parts(14).xaxis = [11 10];
parts(15).id = 14; parts(15).pos = [11 11];  parts(15).xaxis = [11 10];
parts(16).id = 15; parts(16).pos = [12 11];  parts(16).xaxis = [12 11];
parts(17).id = 16; parts(17).pos = [12 12];  parts(17).xaxis = [12 11];
parts(18).id = 17; parts(18).pos = [13 13];  parts(18).xaxis = [13 14];
parts(19).id = 18; parts(19).pos = [13 14];  parts(19).xaxis = [13 14];
parts(20).id = 19; parts(20).pos = [14 14];  parts(20).xaxis = [14 15];
parts(21).id = 20; parts(21).pos = [14 15];  parts(21).xaxis = [14 15];
parts(22).id = 21; parts(22).pos = [15 15];  parts(22).xaxis = [14 15];
parts(23).id = 22; parts(23).pos = [8 8];  parts(23).xaxis = [9 8];
parts(24).id = 23; parts(24).pos = [9 9];  parts(24).xaxis = [9 8];

joints.id    = 0;   joints.child    = parts(1);   joints.parent    = parts(2);   joints.pos    = 0;
joints(2).id = 1;   joints(2).child = parts(2);   joints(2).parent = parts(3);   joints(2).pos = 1;
joints(3).id = 2;   joints(3).child = parts(3);   joints(3).parent = parts(4);  joints(3).pos = 1;
joints(4).id = 3;   joints(4).child = parts(4);   joints(4).parent = parts(5);   joints(4).pos = 2;
joints(5).id = 4;   joints(5).child = parts(5);   joints(5).parent = parts(11);  joints(5).pos = 2;

joints(6).id = 5;   joints(6).child = parts(6);   joints(6).parent = parts(11);  joints(6).pos = 3;
joints(7).id = 6;   joints(7).child = parts(7);   joints(7).parent = parts(6);   joints(7).pos = 3;
joints(8).id = 7;   joints(8).child = parts(8);   joints(8).parent = parts(7);   joints(8).pos = 4;
joints(9).id = 8;   joints(9).child = parts(9);   joints(9).parent = parts(8);  joints(9).pos = 4;
joints(10).id = 9;  joints(10).child = parts(10);  joints(10).parent = parts(9);  joints(10).pos = 4;

joints(11).id = 10; joints(11).child = parts(12); joints(11).parent = parts(11);  joints(11).pos = 7;

joints(12).id = 11; joints(12).child = parts(13); joints(12).parent = parts(14);  joints(12).pos = 10;
joints(13).id = 12; joints(13).child = parts(14); joints(13).parent = parts(15); joints(13).pos = 11;
joints(14).id = 13; joints(14).child = parts(15); joints(14).parent = parts(16); joints(14).pos = 11;
joints(15).id = 14; joints(15).child = parts(16); joints(15).parent = parts(17); joints(15).pos = 12;
joints(16).id = 15; joints(16).child = parts(17); joints(16).parent = parts(11); joints(16).pos = 12;

joints(17).id = 16; joints(17).child = parts(18); joints(17).parent = parts(11);  joints(17).pos = 13;
joints(18).id = 17; joints(18).child = parts(19); joints(18).parent = parts(18); joints(18).pos = 13;
joints(19).id = 18; joints(19).child = parts(20); joints(19).parent = parts(19); joints(19).pos = 14;
joints(20).id = 19; joints(20).child = parts(21); joints(20).parent = parts(20); joints(20).pos = 14;
joints(21).id = 20; joints(21).child = parts(22); joints(21).parent = parts(21); joints(21).pos = 15;
joints(22).id = 21; joints(22).child = parts(24); joints(22).parent = parts(23); joints(22).pos = 8;

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [near,in,on,B,L] = near_mesh(V,F,Q,epsilon)
  % NEAR_MESH test whether a list of points are in or near a given mesh
  %
  % [near] = near_mesh(V,F,Q)
  % [near,in,on,B,L] = near_mesh(V,F,Q)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of face indices
  %   Q  #Q by dim list of query points
  %   epsilon  minium distance allowed after collapses are complete, default is to
  %     use fraction of maximum edge length
  % Outputs:
  %   near #Q list of flags revealing whether queries are near (V,F)
  %   in #Q list of flags revealing whether queries are in (V,F)
  %   on #Q list of flags revealing whether queries are on boundary of (V,F)
  %   B  #B by 1 list of mesh outline edges
  %   L  #loops+1 by 1 list of boundary loop start indices into B, the last
  %     entries is (by tradition) always the numel of B + 1
  %
  % See in_mesh, inpolygon
  %

  if ~exist('epsilon','var') || isempty(epsilon)
    EE = edges(F);
    % maximum edge length
    %maxD = max(sqrt(sum((V(EE(:,1),:) - V(EE(:,2),:)).^2,2)));
    minD = min(sqrt(sum((V(EE(:,1),:) - V(EE(:,2),:)).^2,2)));
    epsilon = minD/2;
  end

  dim = size(V,2);
  % only works in 2D
  assert(dim == 2);

  % first determine points strictly in or on mesh
  [in,on,B,L] = in_mesh(V,F,Q);

  % avoid sqrts
  sqr_eps = epsilon.^2;

  % boundary edges
  BE = [B; B(2:end) B(1)]';
  % compute projection of each point to each boundary line segment
  [T,sqrD] = project_to_lines(Q,V(BE(:,1),:),V(BE(:,2),:));
  % each vertex seen by each edge
  QBE = repmat(Q,[1 1 size(BE,1)]);
  % edge start positions
  S = V(BE(:,1),:);
  % edge destination positions
  D = V(BE(:,2),:);
  % distance of each point to each edge start
  sqrDS = ...
    squeeze(sum((QBE - permute(repmat(S,[1 1 size(Q,1)]),[3 2 1])).^2,2));
  % distance of each point to each edge dest
  sqrDD = ...
    squeeze(sum((QBE - permute(repmat(D,[1 1 size(Q,1)]),[3 2 1])).^2,2));
  % replace distances to edges when point is closest to start or dest endpoints
  % respectively
  sqrD(T<0) = sqrDS(T<0);
  sqrD(T>1) = sqrDD(T>1);
  % compute minimum distance to boundary
  [minD] = min(sqrD,[],2);
  % mask telling whether closest point for each edge is close enough
  near = minD<sqr_eps;
  % all strictly in points are also close
  near = in | near;
end


text:
% SCRIPT TO TEST THE GRAPHIC CAPABILITIES OF THE TOOLBOX: ROBOT IN A MANUFACTURING CELL
% MORE MANUFACTURING DEMOS ARE FOUND UNDER:
%   arte/RAPID/programs
%   The scripts found in this directory can also be translated to RAPID
%   language using the matlab2RAPID function.
%
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
close all
echo on
 q = [0 0 0 0 0 0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% THE DEMO PRESENTS THE ROBOT IN A MANUFACTURING CELL AND SEVERAL BASIC COMMANDS
%
%  START BY LOADING LOAD A ROBOT BY MEANS OF THE LOAD_ROBOT FUNCTION'
%
% >> robot=load_robot('ABB', 'IRB140');
%
%  LOTS OF OTHER ROBOTS CAN BE LOADED, PLEASE TRY:
%
%   robot=load_robot('ABB', 'IRB6620');
%   robot=load_robot('KUKA', 'KR5_arc');
%   robot=load_robot('KUKA', 'KR5_sixx_R650');
%   MORE ROBOTS AVAILABLE UNDER THE arte_lib3.x/robots DIRECTORY
%
%   Example:
%
% >> robot=load_robot('ABB', 'IRB140');
%
% >> adjust_view(robot);
%
% fprintf('\nNOW LOAD AUXILIAR EQUIPMENT')
% >> robot.equipment{1}=load_robot('equipment', 'bumper_cutting');
%
% >> drawrobot3d(robot, q)
%
% >> fprintf('\nNOW LOAD A GRIPPER')
% >> robot.tool=load_robot('equipment/end_tools', 'parallel_gripper_0');
%
% >> drawrobot3d(robot, q)
% >> fprintf('\nNOW THE TEACHING PENDANT')
% >> teach
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo off

fprintf('\nLOAD A ROBOT:\n')

robot=load_robot('ABB', 'IRB140');

fprintf('\nNOW DRAW THE ROBOT')
drawrobot3d(robot, q)

adjust_view(robot);

fprintf('\nNOW LOAD AUXILIAR EQUIPMENT')
robot.equipment{1}=load_robot('equipment', 'bumper_cutting');

%draw the robot
drawrobot3d(robot, q)

fprintf('\nNOW LOAD A GRIPPER')
robot.tool=load_robot('equipment/end_tools', 'parallel_gripper_0');

%draw the robot
drawrobot3d(robot, q)

adjust_view(robot);
robot.tool_activated=1;

drawrobot3d(robot,q)

fprintf('\nNOW LOAD A DIFFERENT TOOL')
robot.tool=load_robot('equipment/end_tools', 'vacuum_2');

drawrobot3d(robot,q)

adjust_view(robot);
%NOW MOVE THE ROBOT WITH THE TEACHING PENDANT

teach

repo: ARTE