text:
function [P,exitCode]=RiccatiPostNoClutter(H,F,R,Q,PD,RelTol,AbsTol,maxIter)
%%RICCATIPOSTNOCLUTTER Solve the posterior Riccati equation that has been
%                      modified for a probability of detection less
%                      than or equal to 1. This provides the average
%                      asymptotic covariance of a linear Kalman filter with
%                      linear measurements just after a (possible with
%                      probability PD) measurement update.
%
%INPUTS: H The zDim X xDim measurement matrix such that H*x+w is the
%          measurement, where x is the state and w is zero-mean Gaussian
%          noise with covariance matrix R.
%        F The xDim X xDim state transition matrix The state at discrete-
%          time k+1 is modeled as F times the state at time k plus zero-
%          mean Gaussian process noise with covariance matrix Q.
%        R The zDim X zDim measurement covariance matrix.
%        Q The xDim X xDim process noise covariance matrix. This can not be
%          a singular matrix.
%       PD The optional detection probability of the target at each scan.
%          If omitted, PD is assumed to be one.
%   RelTol The maximum relative error tolerance allowed, a positive scalar.
%          This tolerance applies to all elements in the matrix and is only
%          needed if PD<1. If omitted or an empty matrix is passed, the
%          default value of 1e-13 is used.
%   AbsTol The absolute error tolerance allowed, a positive scalar. This
%          tolerance applies to all elements in the matrix and is only
%          needed if PD<1. If omitted or an empty matrix is passed, the
%          default value of 1e-10 is used.
%  maxIter An optional integer specifying the maximum number of iterations.
%          If omitted or an empty matrix is passed, the default is 5000.
%
%OUTPUTS: P The average asymptotic error covariance matrix of the linear
%           Kalman filter with a probability of detection of PD after a
%           (possible) measurement update.
%  exitCode This is zero if convergence is achieved and 1 if it is not
%           achieved.
%
%The notion of a Riccati equation with a detection probability less than
%one is discussed in [1]. However, they do not explicitly write out the
%posterior form. The posterior recursion goes from P_{k|k} to P_{k+1|k+1}.
%This stands in contrast to the prior Riccati equation, which expresses the
%transition from P_{k|k-1} to P_{k+1|k}.
%
%The posterior Riccati equation is
%P=F*P*F'+Q-PD*(F*P*F'*H'+Q*H')*inv(H*F*P*F'*H'+H*Q*H'+R)*(F*P*F'*H'+Q*H')'
%and is derived by substituting the prediction step into the measurement
%update step of the linear Kalman filter with additive process noise and
%additive measurement noise.
%
%EXAMPLE:
% T=1;
% F=FPolyKal(T,6,1);
% q0=1;
% Q=QPolyKal(T,6,1,q0);
% R=diag([10;10;10]);
% H=[1,0,0,0,0,0;
%    0,1,0,0,0,0;
%    0,0,1,0,0,0];
% PD=0.5;
% P=RiccatiPostNoClutter(H,F,R,Q,PD)
%
%REFERENCES:
%[1] Y. Boers and H. Driessen, "Modified Riccati equation and its
%    application to target tracking," IEE Proceedings Radar, Sonar and
%    Navigation, vol. 153, no. 1, pp. 7-12, Feb. 2006.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(PD))
   PD=1;
end

if(nargin<6||isempty(AbsTol))
    AbsTol=1e-13;
end

if(nargin<7||isempty(RelTol))
    RelTol=1e-10;
end

if(nargin<8||isempty(maxIter))
   maxIter=5000;
end

%Get an initial estimate by solving the Riccati equation with PD=1.
AD=F';
BD=F'*H';
SD=Q*H';
RD=H*Q*H'+R;
QD=Q;
PPrev=RiccatiSolveD(AD,BD,QD,RD,SD);
exitCode=0;
if(PD~=1)
    %Now, iterate the Riccati equation with PD<1 until convergence of the
    %Frobenis norm or until a maximum number of iterations has occurred.
    curIter=0;
    while(curIter<maxIter)
        P=F*PPrev*F'+Q-PD*(F*PPrev*F'*H'+Q*H')*inv(H*F*PPrev*F'*H'+H*Q*H'+R)*(F*PPrev*F'*H'+Q*H')';
        %Ensure symmetry is preserved.
        P=(P+P')/2;

        diffMag=abs(P-PPrev);
        if(all((diffMag(:)<=RelTol*abs(P(:)))|(diffMag(:)<=AbsTol)))
            return;
        end

        PPrev=P;
        curIter=curIter+1;
    end
    exitCode=1;%Convergence not achieved.
else
    P=PPrev;
end

end

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

author: USNavalResearchLaboratory


text:
function DetectionList=OSCFAR(data2D,numGuardCells,numAvgCells,k,PFA)
%%OSCFAR Perform order-statistic constant false alarm rate (OS-CFAR)
%        detection on a two-dimensional grid. This can be, for example, a
%        range-Doppler plot. This function just performs the detection, it
%        does not centroid the detections.
%
%INPUTS: data2D The delay-Doppler plot or a set of delay-Doppler plots.
%               This is a numRowsXnumColsXnumPlots set of numPlots 2D
%               range-Doppler maps (or similar matrices on which CFAR
%               should be performed). This can contain complex values.
% numGuardCells The integer number of cells in each dimension around each
%               test cell that are not considered in the average used for
%               determining the changing detection threshold. This is a 2X1
%               vector. If the same value is used in both directions, then
%               a scalar can be passed.
%   numAvgCells The width of the region in cells in each dimension after
%               the guard cell region that define the average used for
%               determinig the threshold. This is a 2X1 vector. If the same
%               value is used in both directions, then a scalar can be
%               passed.
%             k The integer order to use k. That is, the kth largest sample
%               in the test region is used as the test statistic.
%           PFA The scalar probability of false alarm between 0 and 1 that
%               determines the threshold for detection.
%
%OUTPUTS: DetectionList A numPlotsX1 collection of structures.
%               DetectionList(i).Index provides a 2XnumDetect set of the
%               row and column indices of each detection in the ith plot.
%               The values from data2D of the detections are given in
%               DetectionList(i).Value
%
%This implements the algorithm described in Section 4 of [1].
%
%EXAMPLE:
% fc=1e9;%1GHz carrier frequency.
% B=2e6;%2Mhz bandwidth.
% %Baseband start and end frequencies.
% fStart=-B/2;
% fEnd=B/2;
% %Sampling rate is two times the Nyquist rate.
% T0=1/(2*2*fEnd);%Sampling period in seconds.
% T=2e-5;%Chirp duration in seconds.
%
% PRF=2000;%Pulse repetition frequency (Hertz)
% TB=1/PRF;%The pulse repetition period.
% %The number of samples per PRI. The above parameters were chosen so that
% %this is an integer. Fix just deals with finite precision errors.
% Ns=fix(TB/T0);
%
% %Generate the reference signal. This is an up-chirp.
% x=LFMChirp(T,fStart,fEnd,T0);
% x=x(:);
%
% %We will use 64 pulse repetition intervals.
% NB=64;
%
% %True target parameters.
% c=Constants.speedOfLight;
% rTrue=[40e3;50e3;60e3];
% tau=rTrue/c;%The true delay (s).
%
% %The true range rate (m/s).
% rrTrue=[100;60;-100];
% a=rrTrue/c;
%
% %Complex amplitudes
% A=[48;512*exp(1j*2*pi*rand(1));32*exp(1j*2*pi*rand(1))];
%
% numTargets=length(rTrue);
%
% %Allocate space for the received signal. The first dimensions is "fast
% %time"; the second dimension if "slow time".
% y=zeros(Ns,NB);
%
% %Create the received signal, properly delayed and Doppler shifted for
% %each PRI. The same waveform is used in each PRI.
% t=0:T0:((Ns-1)*T0);%Sample times
% for i=0:(NB-1)
%
%     for curTar=1:numTargets
%         tCur=t-a(curTar)*t-tau(curTar)-i*TB;
%         %The signal simulated with range migration.
%         y(:,i+1)=y(:,i+1)+(A(curTar)*exp(-1j*2*pi*fc*(tau(curTar)+a(curTar)*t)).*LFMChirp(T,fStart,fEnd,tCur)).';
%     end
%
%     y(:,i+1)=y(:,i+1)+ComplexGaussianD.rand(Ns).';
%
%     t=t+TB;%Increment to the next time step.
% end
%
% M1=1;
% M2=1;
%
% %Windowing in range and Doppler to lower sidelobes.
% Nx=length(x);
% wRange=windowFunSym(Nx,'Blackman',0);
% x=wRange.*x;
% wDoppler=windowFunSym(NB,'Nuttall',1);
%
% [delayDopPlot,Doppler,delay]=delayDopplerPlotNBPulseDop(x,y,M1,M2,wDoppler,T0);
%
% numGuardCells=[2;4];
% numAvgCells=[5;3];
% PFA=1e-7;
%
% range=c*delay;
% rangeRate=Doppler*(c/fc);
%
% %Display the range-Doppler plot
% figure(1)
% clf
% imagesc([rangeRate(1),rangeRate(end)],[range(1), range(end)]/1e3,10*log10(abs(delayDopPlot)));
% set(gca,'YDir','normal')
% caxis([-30 27])
% colormap(jet(256))
% h1=xlabel('Range Rate (m/s)');
% h2=ylabel('Range (km)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
% k=165;
% DetectionList=OSCFAR(delayDopPlot,numGuardCells,numAvgCells,k,PFA);
% idx=DetectionList(1).Index;
%
% rVals=range(idx(1,:));
% DopVals=rangeRate(idx(2,:));
%
% %Display the CFAR detections.
% figure(2)
% clf
% hold on
% scatter(DopVals,rVals)
% h1=xlabel('Range Rate (m/s)');
% h2=ylabel('Range (km)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% axis([rangeRate(1), rangeRate(end), range(1), range(end)])
%
%REFERENCES
%[1] P. P. Gandhi and S. A. Kassam, "Analysis of CFAR processors in
%    nonhomogeneous background," IEEE Transactions on Aerospace and
%    Electronic Systems, vol. 24, no. 4, pp. 427-445, Jul. 1988.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(length(numGuardCells)==1)
    numGuardCells=[numGuardCells;numGuardCells];
end

if(length(numAvgCells)==1)
    numAvgCells=[numAvgCells;numAvgCells];
end

NFilter=numGuardCells+numAvgCells;

numRows=size(data2D,1);
numCols=size(data2D,2);
numEls=numRows*numCols;

numPlots=size(data2D,3);

%The number of elements in the mask.
NMask=2*numGuardCells+2*numAvgCells+1;

%Create the mask. The ones in the mask select which cells below the mask
%are considered  when performing the sorting and then averaging.
Mask1=ones(NMask(:)');
Mask2=zeros(NMask(:)');
Mask2((numAvgCells(1)+1):(numAvgCells(1)+2*numGuardCells(1)+1),(numAvgCells(2)+1):(numAvgCells(2)+2*numGuardCells(2)+1))=1;
Mask=logical(Mask1-Mask2);

%The number of cells in the mask
NCFAR=prod(2*NFilter+1)-prod(2*numGuardCells+1);

%The threshold for detection.
T=OSCFARThreshold4PFA(PFA,NCFAR,k);

%Allocate space for the return variables.
DetectionList(numPlots).Index=[];
DetectionList(numPlots).Value=[];

%For each of the range-Doppler plots that was passed.
for curPlot=1:numPlots
    plotMag2=abs(data2D(:,:,curPlot)).^2;

    %Preallocate the maximum amount of space for detections. It can get
    %shrunk on return.
    indexR=zeros(numEls,1);
    indexD=zeros(numEls,1);
    numDet=0;

    %Extend
    plotMag2=addAliasedPadding(plotMag2,NFilter);

    %The initial span of the rows covered by the mask.
    rowSpan=1:NMask(1);
    for curRow=1:numRows
        %The initial span of the columns covered by the mask.
        colSpan=1:NMask(2);
        for curCol=1:numCols
            curVals=plotMag2(rowSpan,colSpan);
            sortedMaskedVals=sort(curVals(Mask),'ascend');

            if(sortedMaskedVals(k)*T<plotMag2(NFilter(1)+curRow,NFilter(2)+curCol))
                numDet=numDet+1;
                indexR(numDet)=curRow;
                indexD(numDet)=curCol;
            end
            colSpan=colSpan+1;
        end
        rowSpan=rowSpan+1;
    end

    %Shink to fit the actual detections.
    indexR=indexR(1:numDet);
    indexD=indexD(1:numDet);

    DetectionList(curPlot).Index=[indexR.';indexD.'];
    idx=sub2ind([numRows,numCols],indexR,indexD)+(curPlot-1)*numRows*numCols;

    DetectionList(curPlot).Value=data2D(idx);
end
end

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

repo: TrackerComponentLibrary


text:
function ncomb = multinomial_coef2 ( nfactor, factor )

%*****************************************************************************80
%
%% MULTINOMIAL_COEF2 computes a Multinomial coefficient.
%
%  Discussion:
%
%    The multinomial coefficient is a generalization of the binomial
%    coefficient.  It may be interpreted as the number of combinations of
%    N objects, where FACTOR(1) objects are indistinguishable of type 1,
%    ... and FACTOR(NFACTOR) are indistinguishable of type NFACTOR,
%    and N is the sum of FACTOR(1) through FACTOR(NFACTOR).
%
%    NCOMB = N! / ( FACTOR(1)! FACTOR(2)! ... FACTOR(NFACTOR)! )
%
%    A direct method is used, which should be exact.  However, there
%    is a possibility of intermediate overflow of the result.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NFACTOR, the number of factors.
%    1 <= NFACTOR.
%
%    Input, integer FACTOR(NFACTOR), contains the factors.
%    0.0 <= FACTOR(I).
%
%    Output, integer NCOMB, the value of the multinomial coefficient.
%
  ncomb = 1;
  k = 0;

  for i = 1 : nfactor

    for j = 1 : factor(i)
      k = k + 1;
      ncomb = ( ncomb * k ) / j;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% dtiEpilepsy
%
% This script is applied to compare a single subject with epilepsy (an
% adult) with the template of comparable adults without epilepsy.  We are
% using the SIRL20 template for the moment.
%
% To run this script you will need the statistics toolbox on your path.
% This toolbox is part of our normal lab Matlab and it is part of the Unix
% side.  We don't have it yet for license 11723 version.  We should get it.
%

% Specify and load the subject data and the atlas data
%
%  Does the subject need to be spatially aligned with the template? Usually
%  this is done once for each dti scan.  If you haven't done it, you can
%  get it done here.
doSpatialNorm = true;

% Is there a warping needed for the B0 field?  Radiation necrosis needed
% this, but normally we don't.  See paragraph below.  The b0Norm is only
% invoked during the doSpatialNorm operation.
b0Norm = false;

% There are several possible templates.  They should be listed here.
templateName = 'SIRL20adult';

% Set up the files and directories.  We have had trouble reading across the
% network from time to time.  So on PCs, we copy the key files locally.  On
% Linux boxes, everything works in place.
if(ispc)
    dataDir = '\\white.stanford.edu\biac2-wandell2\data\Epilepsy';
    templateDir = '\\white.stanford.edu\biac2-wandell2\data\templates\adult';
else
    dataDir = '/biac2/wandell2/data/Epilepsy';
    templateDir = '/biac2/wandell2/data/templates/adult';
end
avgdir = fullfile(templateDir,[templateName 'warp3_averageDataset']);
% Full path to the file
%subjectDt6 = epiSubjectFile(dataDir, 1,'postictal');
subjectDt6 = epiSubjectFile(dataDir, 1,'control');

% This normalization operation should be a integrated script, we think.  It
% requires SPM2.  It requires comments.
if(doSpatialNorm)
    if(ispc)
        % To avoid a memory paging error in windows, we copy the templates to a
        % local directory.
        copyfile(fullfile(templateDir,[templateName '_brain.img']), fullfile(tempdir,[templateName '_brain.img']));
        copyfile(fullfile(templateDir,[templateName '_brain.hdr']), fullfile(tempdir,[templateName '_brain.hdr']));
        copyfile(fullfile(templateDir,[templateName '_EPI.img']), fullfile(tempdir,[templateName '_EPI.img']));
        copyfile(fullfile(templateDir,[templateName '_EPI.hdr']), fullfile(tempdir,[templateName '_EPI.hdr']));
        copyfile(fullfile(templateDir,[templateName 'warp3_averageDataset'],'tensorSummary.mat'), fullfile(tempdir,['tensorSummary.mat']));
        copyfile(fullfile(templateDir,[templateName 'warp3_averageDataset'],'average_dt6.mat'), fullfile(tempdir,['average_dt6.mat']));

        % On the PC, we put all the files in one directory, the temp
        % directory for this user.
        templateDir = tempdir;
        avgdir = templateDir;
    end
    dtiSpatialNormalize;
else
    % No normalization needed
    subjectDt6 = [subjectDt6 '_' templateName '.mat'];
end

% Needs further definition.
tensorSumFile = fullfile(avgdir,'tensorSummary.mat');

% VOXELWISE ANALYSIS
% We are not sure about the difference between the template, the
% average and the tensor summary.
% We think that the
if(exist(tensorSumFile,'file'))
    disp(['Loading control group tensor summary (' tensorSumFile ')']);
    % This loads various summary structures such as md, fa, logTensor,
    % meanB0 and notes.
    load(tensorSumFile);
else
    % Compute the tensor summary.
    dtiTensorSummary;
end

% Load the atlas information
template = load(fullfile(avgdir,'average_dt6'));
xformDtToAcpc = template.xformToAcPc;

% Directory containing this subject's data
subDir = fileparts(subjectDt6);

% Single subject DTI data, spatially normalized
disp(['Loading dt6 file (' subjectDt6 ')']);
ssDt_sn = load(subjectDt6);
%mask = mask&ssDt_sn.dt6(:,:,:,1)>0;

%figure; imagesc(makeMontage(template.b0)); axis image; colormap gray;
%figure; imagesc(makeMontage(ssDt_sn.b0)); axis image; colormap gray;

%------------------------
% At this point the atlas data and the subject are loaded.
% We begin statistical analysis.  We begin with FA and MD analyses in this
% section. Then we move on to the full tensor analyses in the following
% section.

% The mask is the valid data mask computed  within dtiTensorSummary.
ssDt_ind = dtiImgToInd(ssDt_sn.dt6, mask);
[eigVec, eigVal] = dtiEig(ssDt_ind);

% Compute a test comparing FA in this subject with the atlas. This is done
% by a t-test implemented in dtiTTest.
[ss_fa, ss_md]   = dtiComputeFA(eigVal);

% Create the T-test images and related data for FA
[Tfa, DISTR, df] = dtiTTest(fa.mean, fa.stdev, fa.n, ss_fa);
[TfaImg, tFDRfa, n_signif,index_signif, pvals] = dtiTTestImage(Tfa, DISTR, df, mask);

% showSlices = [15:62];
% figure; imagesc(makeMontage(TfaImg,showSlices)); axis image; colormap cool; colorbar;
% set(gcf,'Name','FA test'); title(sprintf('tthresh, no FDR (p<10^-^4) = %0.1f',tThresh));

% Save out images that can be viewed in dtiFiberUI
dtiWriteNiftiWrapper(TfaImg, xformDtToAcpc, fullfile(subDir,['fa_' DISTR '-test_' num2str(df(1)) 'df.nii.gz']));
dtiWriteNiftiWrapper(-1*TfaImg, xformDtToAcpc, fullfile(subDir,['negfa_' DISTR '-test_' num2str(df(1)) 'df.nii.gz']));
dtiWriteNiftiWrapper(dtiIndToImg(fa.stdev,mask), xformDtToAcpc, fullfile(subDir,'fa_variance.nii.gz'));

% Repeat the process for mean diffusivity
[Tmd, DISTR, df] = dtiTTest(md.mean, md.stdev, md.n, ss_md);
[TmdImg, tFDRmd, n_signif,index_signif, pvals] = dtiTTestImage(Tmd, DISTR, df, mask);

% showSlices = [15:62];
% figure; imagesc(makeMontage(TmdImg,showSlices)); axis image; colormap cool; colorbar;
% set(gcf,'Name','MD test'); title(sprintf('tthresh, no FDR (p<10^-^4) = %0.1f',tThresh));

% Save out images that can be viewed in dtiFiberUI.  Until we fix the
% dtiFiberUI, we can't really select out the negative values.  So we write
% the data out both positive and negative for the moment.  We need to
% amend dtiFiberUI so that it
dtiWriteNiftiWrapper(TmdImg, xformDtToAcpc, fullfile(subDir,['md_' DISTR '-test_' num2str(df(1)) 'df.nii.gz']));
dtiWriteNiftiWrapper(-1*TmdImg, xformDtToAcpc, fullfile(subDir,['negmd_' DISTR '-test_' num2str(df(1)) 'df.nii.gz']));
dtiWriteNiftiWrapper(dtiIndToImg(md.stdev,mask), xformDtToAcpc, fullfile(subDir,'md_variance.nii.gz'));

% To view the locations of the significantly different voxels, use
% dtiFiberUI.  Load the file, say, fa_t-test_19df.nii.gz as a NIFTI image
% overlay into the program.  You can adjust the cutoff in the user
% interface.  You can use, say, the fdr threshold (tFDRxx) as the p < 0.05
% cutoff.

% clear TmdImg TfaImg Tmd Tfa

%--------------------------------------
% The next two analyses require log-space tensors.

% Log-transform the single subject data
eigVal(eigVal<0) = 0;
eigVal = log(eigVal);
ssDt_ind = dtiEigComp(eigVec, eigVal);

% Test for VECTOR differences
[T, DISTR, df] = dtiLogTensorTest('vec', logTensor.mean, logTensor.stdev, logTensor.n, ssDt_ind);

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(:));
Simg = dtiIndToImg(logTensor.stdev,mask);

% figure; imagesc(makeMontage(Timg,showSlices)); axis image; colormap hot; colorbar;
% set(gcf,'Name','Vec test'); title(sprintf('fthresh (p<10^-^4) = %0.1f',fThresh));
%
% figure; imagesc(makeMontage(Simg,showSlices)); axis image; colormap hot; colorbar;
% set(gcf,'Name','Vec test variance');

%
dtiWriteNiftiWrapper(Timg, xformDtToAcpc, fullfile(subDir,['vec_' DISTR '-test_' num2str(df(1)) ',' num2str(df(2)) 'df.nii.gz']));
dtiWriteNiftiWrapper(Simg, xformDtToAcpc, fullfile(subDir,'vec_variance.nii.gz'));
disp(['dtiFiberUI threshold: ' num2str(fThresh/fMax)]);

% Test for VALUE differences
[T, DISTR, df] = dtiLogTensorTest('val', logTensor.mean, logTensor.stdev, logTensor.n, ssDt_ind);

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(:));
Simg = dtiIndToImg(logTensor.stdev,mask);

% figure; imagesc(makeMontage(Timg,showSlices)); axis image; colormap hot; colorbar;
% set(gcf,'Name','Values test'); title(sprintf('fthresh (p<10^-^4) = %0.1f',fThresh));
%
% figure; imagesc(makeMontage(Simg,showSlices)); axis image; colormap hot; colorbar;
% set(gcf,'Name','Values test variance');

%
dtiWriteNiftiWrapper(Timg, xformDtToAcpc, fullfile(subDir,['val_' DISTR '-test_' num2str(df(1)) ',' num2str(df(2)) 'df.nii.gz']));
dtiWriteNiftiWrapper(Simg, xformDtToAcpc, fullfile(subDir,'val_variance.nii.gz'));
disp(['dtiFiberUI threshold: ' num2str(fThresh/fMax)]);

% figure; imagesc(makeMontage(sqrt(Timg),[20:55])); axis image; colormap hot; colorbar;
% set(gcf,'Name','Vec Standardized Distance');

% Possible to test for 'full' also.  Need more comments in dtiLogTensorTest

% FDR analysis for eigenvector differences
%
% Sqrt(F) is the standardized distance between the groups.  Ask about this.
%
fdrVal = 0.05; fdrType = 'general';
T(isnan(T)) = 0;
pvals = 1-fcdf(T, df(1), df(2));
[n_signif,index_signif] = fdr(pvals,fdrVal,fdrType,'mean');
disp(n_signif);max(pvals(index_signif))

% Convert back to an fThreshold
pThreshFDR = max(pvals(index_signif));
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));

%---- Make locally viewable versions of the images ----
%
logPimg = dtiIndToImg(-log10(pvals), mask);
cmap = autumn(256); maxLogP = 10; minLogP = -log10(pThreshFDR);

anatRgb = repmat(mrAnatHistogramClip(double(ssDt_sn.anat.img),0.4,0.98),[1,1,1,3]);
tmp = mrAnatResliceSpm(logPimg, inv(ssDt_sn.xformToAcPc), [], ssDt_sn.anat.mmPerVox, [1 1 1 0 0 0]);
tmp(tmp>maxLogP) = maxLogP;
tmp = (tmp-minLogP)./(maxLogP-minLogP);
overlayMask = (tmp>=0);
tmp(~overlayMask) = 0;
overlayMask = repmat(overlayMask,[1 1 1 3]);
overlayRgb = reshape(cmap(round(tmp*255+1),:),[size(tmp) 3]);
anatRgb(overlayMask) = overlayRgb(overlayMask);

% reorient so that the eyes point up
anatRgb = flipdim(permute(anatRgb,[2 1 3 4]),1);
%sl = [2:2:40];
sl = [-36:2:60];
for(ii=1:length(sl)) slLabel{ii} = sprintf('Z = %d',sl(ii)); end
slImg = inv(ssDt_sn.anat.xformToAcPc)*[zeros(length(sl),2) sl' ones(length(sl),1)]';
slImg = round(slImg(3,:));
anatOverlay = makeMontage3(anatRgb, slImg, ssDt_sn.anat.mmPerVox(1), 0, slLabel);
mrUtilPrintFigure(fullfile(subDir,'ss_t1_vecSPM'));
legendLabels = explode(',',sprintf('%0.1f,',[minLogP:1:maxLogP]));
legendLabels{end} = ['>=' num2str(maxLogP)];
mrUtilMakeColorbar(cmap, legendLabels, '-log10(p)', fullfile(subDir,'vecSPM_legend'));

templateBrain = template.anat.img;
templateBrain(template.anat.brainMask<0.25) = 0;
templateBrain(template.anat.brainMask<0.5) = templateBrain(template.anat.brainMask<0.5)*.5;
avgRgb = repmat(templateBrain,[1,1,1,3]);
avgRgb(overlayMask) = overlayRgb(overlayMask);

% reorient so that the eyes point up
avgRgb = flipdim(permute(avgRgb,[2 1 3 4]),1);
avgOverlay = makeMontage3(avgRgb, slImg, ssDt_sn.anat.mmPerVox(1), 0, slLabel);

slImg = inv(ssDt_sn.xformToAcPc)*[zeros(length(sl),2) sl' ones(length(sl),1)]';
slImg = round(slImg(3,:));
overlayRgb = flipdim(permute(overlayRgb,[2 1 3 4]),1);
overlay = makeMontage3(overlayRgb, slImg, [], 2);

% figure; imagesc(makeMontage(logPimg,showSlices)); axis image; colormap hot; colorbar;
threshMask = zeros(size(pvals));
threshMask(index_signif) = 1;
threshMask = dtiIndToImg(threshMask, mask);
img = logPimg; img(threshMask<1) = 0;
%  figure; imagesc(makeMontage(img,showSlices)); axis image; colormap hot; colorbar;

%%%END HERE%%%

%--------------------------------------------------------------------
% FDR Analysis

% Quantile transformation
Tchisq = chi2inv(cdf(DISTR, T, df(1), df(2)), df(1));
DISTR = 'chi2';
% Histograms
H = fdrHist(Tchisq,0.2,1);

% Empirical null
Tmax = prctile(Tchisq,90);
w = (H.x < Tmax);
[params, paramsCov, H0hat] = fdrEmpNull(H, w, DISTR, {}, df);
[p0, s0, df0] = deal(params(1), params(2), params(3));
paramsConf = [[log(p0), s0, df0] + 1.95*sqrt(diag(paramsCov))';
    [log(p0), s0, df0] - 1.95*sqrt(diag(paramsCov))'];

% p0 adjustment for theoretical null
[params, paramsCov, H0] = fdrEmpNull(H, w, DISTR, {'df','s'}, df);
p0H0 = params(1);

% FDR curves
[fdrH0, t] = fdrCurveHist('FDR', H0, 1);
[fdrH0hat, t] = fdrCurveHist('FDR', H0hat, 1);

% Threshold
thr = fdrThresh(1, fdrH0(:,1), t, level);
thr = fdrThresh(1, fdrH0hat(:,1), t, level);

%-----------------------------------------------------------------------
% FDR Plots

% Histogram of test stats
figure, set(gcf, 'name', 'Histograms'), hold on
h = bar(H.x, H.hist, 1, 'w');
h0 = plot(H0.x, H0.yhat, 'b');
h1 = plot(H0hat.x, H0hat.yhat, 'r');
hold off, legend([h0 h1], 'theo null','emp null',1)
xlabel('T'); ylabel('voxel count');
a=axis; axis([0 prctile(T,99) a(3:4)]);

% FDR curves
figure,	set(gcf, 'name', 'FDR'), hold on
plot(t, fdrH0(:,1), 'b')
plot(t, fdrH0hat(:,1), 'r')
legend('theo null','emp null')
plot(t, fdrH0(:,2), 'b:', t, fdrH0(:,3), 'b:')
plot(t, fdrH0hat(:,2), 'r:', t, fdrH0hat(:,3), 'r:')
hold off, axis([0 prctile(Tchisq,99.99) 0 1]);
xlabel('threshold'); ylabel('FDR');

pthresh = 1-chi2cdf(t,df(1));
figure;plot(log10(pthresh), fdrH0(:,1), 'b')
xlabel('log10(p-value)'); ylabel('FDR');

% Test for eigenvalue differences
[T, M, S, DISTR, df] = dtiLogTensorTest(1, [2:size(allDt6_ind,3)], allDt6_ind, 'val');
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,[20:55])); axis image; colormap hot; colorbar;
set(gcf,'Name','Val test'); title(sprintf('fthresh (p<10^-^4) = %0.1f',fThresh));
dtiWriteNiftiWrapper(Timg, xformDtToAcpc, fullfile(subDir,['val_' DISTR '-test_' num2str(df(1)) ',' num2str(df(2)) 'df.nii.gz']));
dtiWriteNiftiWrapper(dtiIndToImg(S,mask), xformDtToAcpc, fullfile(subDir,'val_variance.nii.gz'));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function cdf = binomial_cdf ( x, a, b )

%*****************************************************************************80
%
%% BINOMIAL_CDF evaluates the Binomial CDF.
%
%  Discussion:
%
%    CDF(X)(A,B) is the probability of at most X successes in A trials,
%    given that the probability of success on a single trial is B.
%
%    A sequence of trials with fixed probability of success on
%    any trial is known as a sequence of Bernoulli trials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the desired number of successes.
%    0 <= X <= A.
%
%    Input, integer A, the number of trials.
%    1 <= A.
%
%    Input, real B, the probability of success on one trial.
%    0.0D+00 <= B <= 1.0.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x < 0 )

    cdf = 0.0;

  elseif ( a <= x )

    cdf = 1.0;

  elseif ( b == 0.0 )

    cdf = 1.0;

  elseif ( b == 1.0 )

    cdf = 0.0;

  else

    cdf = 0.0;

    for j = 0 : x

      cnk = binomial_coef ( a, j );

      pr = cnk * b^j * ( 1.0 - b )^( a - j );

      cdf = cdf + pr;

    end

  end

  return
end


text:
function  [tscounts, wfcounts, evcounts, contcounts] = plx_info(filename, fullread)
% plx_info(filename, fullread) -- read and display .plx or .pl2 file info
%
% [tscounts, wfcounts, evcounts, contcounts] = plx_info(filename, fullread)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   fullread - if 0, reads only the file header
%              if 1, reads the entire file
%               for .pl2 files, this parameter is ignored
%
% OUTPUT:
%   tscounts - 2-dimensional array of timestamp counts for each unit
%      tscounts(i, j) is the number of timestamps for channel j-1, unit i
%                                (see comment below)
%   wfcounts - 2-dimensional array of waveform counts for each unit
%     wfcounts(i, j) is the number of waveforms for channel j-1, unit i
%                                (see comment below)
%   evcounts - 1-dimensional array of external event counts
%     evcounts(i) is the number of events for event channel i
%
%   contcounts - 1-dimensional array of sample counts for continuous channels
%     contcounts(i) is the number of continuous for slow channel i-1
%
% Note that for tscounts, wfcounts, the unit,channel indices i,j are off by one.
% That is, for channels, the count for channel n is at index n+1, and for units,
%  index 1 is unsorted, 2 = unit a, 3 = unit b, etc
% The dimensions of the tscounts and wfcounts arrays are
%   (NChan+1) x (MaxUnits+1)
% where NChan is the number of spike channel headers in the plx file, and
% MaxUnits is 4 if fullread is 0, or 26 if fullread is 1. This is because
% the header of a .plx file can only accomodate 4 units, but doing a
% fullread on the file may show that there are actually up to 26 units
% present in the file. Likewise, NChan will have a maximum of 128 channels
% if fullread is 0.
% The dimension of the evcounts and contcounts arrays is the number of event
% and continuous (slow) channels.
% The counts for slow channel 0 is at contcounts(1)

tscounts = [];
wfcounts = [];
evcounts = [];
contcounts = [];

if nargin ~= 2
    error 'expected 2 input arguments';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2 = PL2GetFileIndex(filename);
    numSpikeChannels = numel(pl2.SpikeChannels);
    % pl2 files support up to 256 units, but we limit to
    % 26 sorted plus 1 unsorted to be compatible with plx_info
    tscounts = zeros(27,numSpikeChannels+1);
    for i=1:numSpikeChannels
        tscounts(:,i+1) = pl2.SpikeChannels{i}.UnitCounts(1:27);
    end
    wfcounts = tscounts;

    numAnalogChannels = numel(pl2.AnalogChannels);
    contcounts = zeros(1,numAnalogChannels);
    for i=1:numAnalogChannels
        contcounts(1,i) = pl2.AnalogChannels{i}.NumValues;
    end

    evcounts = pl2.EventCounts;
    return
end

[tscounts, wfcounts, evcounts, contcounts] = mexPlex(4, filename, fullread);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function  [bnde,bnd] = extdom_edges2( t , p )
% bnde -- boundary edges
% bnd  -- boundary points
% Given an element table, extract boundary edges in no order.
edge = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])];                               % Non-unique edges
edge = sortrows(sort(edge,2));                                             % Put shared edges next to each other
idx = all(diff(edge)==0,2);                                                % Find shared edges
idx = [idx;false]|[false;idx];                                             % True for all shared edges
bnde = edge(~idx,:);                                                       % Boundary edges
bnde  = unique(bnde,'rows');                                               % Boundary nodes
bnd=p(unique(bnde(:)),:);
end

author: CHLNDDEV


text:
%This file will generate the library we prepared for the MMK example.
% Last Updated: 2019/04/22
% Coded By: K

function [Data,Sym_Struct]=SINDyLib(X,dX,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order)
%% First get the size of the X matrix, determin the data length and the number of variables we have.
[Data_Length,Variable_Number]=size(X);
[~,Variable_Number_dX]=size(dX);
[~,Variable_Number_u]=size(u);
%Also create the symbolic variable
Symbol=sym('z',[Variable_Number,1]);
Symbol_dX=sym('dz',[Variable_Number,1]);
Symbol_u=sym('u',[Variable_Number_u,1]);

%Now according the Highest Polynomial Order entered, we will calculate the data matrix.
Data=[];
Index=1;

%% First calculate the polynomial term
%Order zero:
Index=1;
Data(:,Index)=ones(Data_Length,1);
Sym_Struct{1,Index}=1;

%Order One:
if Highest_Poly_Order>=1
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=X(:,i);
        Sym_Struct{1,Index}=Symbol(i,1);
    end
end

%Order Two:
if Highest_Poly_Order>=2
    for i=1:Variable_Number
        for j=i:Variable_Number
            Index=Index+1;
            Data(:,Index)=X(:,i).*X(:,j);
            Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1);
        end
    end
end

%Order Three:
if Highest_Poly_Order>=3
    for i=1:Variable_Number
        for j=i:Variable_Number
            for k=j:Variable_Number
                Index=Index+1;
                Data(:,Index)=X(:,i).*X(:,j).*X(:,k);
                Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1)*Symbol(k,1);
            end
        end
    end
end

%Order Four:
if Highest_Poly_Order>=4
    for i=1:Variable_Number
        for j=i:Variable_Number
            for k=j:Variable_Number
                for pi=k:Variable_Number
                    Index=Index+1;
                    Data(:,Index)=X(:,i).*X(:,j).*X(:,k).*X(:,pi);
                    Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1)*Symbol(k,1)*Symbol(pi,1);
                end
            end
        end
    end
end

%% Then add the Trigonometric Function in the output data:

%Order One:
if Highest_Trig_Order>=1
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=sin(X(:,i));
        Sym_Struct{1,Index}=sin(Symbol(i,1));
    end
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=cos(X(:,i));
        Sym_Struct{1,Index}=cos(Symbol(i,1));
    end
end

%Order Two:
if Highest_Trig_Order>=2
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)+X(:,j));
            Sym_Struct{1,Index}=sin(Symbol(i,1)+Symbol(j,1));
        end
    end
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)+X(:,j));
            Sym_Struct{1,Index}=cos(Symbol(i,1)+Symbol(j,1));
        end
    end
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-X(:,j));
            Sym_Struct{1,Index}=sin(Symbol(i,1)-Symbol(j,1));
        end
    end
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-X(:,j));
            Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1));
        end
    end
    %
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-2*X(:,j));
            Sym_Struct{1,Index}=sin(Symbol(i,1)-2*Symbol(j,1));
        end
    end
    %
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-2*X(:,j));
            Sym_Struct{1,Index}=cos(Symbol(i,1)-2*Symbol(j,1));
        end
    end
end

%Order Three:
if Highest_Trig_Order>=3
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-X(:,j)).^2;
            Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1))^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-X(:,j)).^2;
            Sym_Struct{1,Index}=sin(Symbol(i,1)-Symbol(j,1))^2;
        end
    end
end

%Order Four
if Highest_Trig_Order>=4
    %
    for i=1:Variable_Number
        for j=i:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(2*X(:,i)-2*X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=sin(2*Symbol(i,1)-2*Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=i:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(2*X(:,i)-2*X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=cos(2*Symbol(i,1)-2*Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=sin(X(:,i)-X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=sin(Symbol(i,1)-Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        for j=i+1:Variable_Number
            Index=Index+1;
            Data(:,Index)=cos(X(:,i)-X(:,j)).*X(:,i).^2;
            Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1))*Symbol(i,1)^2;
        end
    end
    %
    for i=1:Variable_Number
        j=2;
        Index=Index+1;
        Data(:,Index)=cos(X(:,i)-X(:,j)).*sin(X(:,i));
        Sym_Struct{1,Index}=cos(Symbol(i,1)-Symbol(j,1))*sin(Symbol(i,1));
    end
end

pin=Index;

%% From here, we add the dX*Theta elements in our data.
if Highest_dPoly_Order>=1
    for j=1:Variable_Number_dX
    for k=1:pin
        Index=Index+1;
        Data(:,Index)=dX(:,j).*Data(:,k);
        Sym_Struct{1,Index}=Symbol_dX(j,1)*(Sym_Struct{1,k});
    end
    end
end

%% Frome here, we add the u*Theta elements in our data
if Highest_U_Order>=1
    for j=1:Variable_Number_u
        for k=1:pin
            Index=Index+1;
            Data(:,Index)=u(:,j).*Data(:,k);
            Sym_Struct{1,Index}=Symbol_u(j,1)*(Sym_Struct{1,k});
        end
    end
end

author: dynamicslab

repo: SINDy-PI


text:
%PLOT_ELLIPSE_INV Draw an ellipse or ellipsoid
%
% PLOT_ELLIPSE_INV(A, OPTIONS) draws an ellipse defined by X'.inv(A).X = 0 on the
% current plot, centred at the origin.
%
% PLOT_ELLIPSE_INV(A, C, OPTIONS) as above but centred at C=[X,Y].  If
% C=[X,Y,Z] the ellipse is parallel to the XY plane but at height Z.
%
% H = PLOT_ELLIPSE_INV(A, C, OPTIONS) as above but return graphic handle.
%
% Options::
% 'edgecolor'   the color of the circle's edge, Matlab color spec
% 'fillcolor'   the color of the circle's interior, Matlab color spec
% 'alpha'       transparency of the filled circle: 0=transparent, 1=solid
% 'alter',H     alter existing circles with handle H
%
% Notes::
% - For the case where the inverse of ellipse parameters are known, perhaps
%   an inverse covariance matrix.
% - If A (2x2) draw an ellipse, else if A(3x3) draw an ellipsoid.
% - The ellipse is added to the current plot.
%
% See also PLOT_ELLIPSE, PLOT_CIRCLE, PLOT_BOX, PLOT_POLY.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
% See also PLOT_ELLIPSE, PLOT_CIRCLE, PLOT_BOX, PLOT_POLY.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function h = plot_ellipse_inv(A, xc, varargin)

    if nargin == 1
        h = plot_ellipse(inv(A));
    elseif nargin == 2
        h = plot_ellipse(inv(A), xc);
    else
        h = plot_ellipse(inv(A), xc, varargin{:});
    end

author: Allopart

repo: rbpf-gmapping


text:
% This is the main FastSLAM loop. This script calls all the required
% functions in the correct order.
%
% You can disable the plotting or change the number of steps the filter
% runs for to ease the debugging. You should however not change the order
% or calls of any of the other lines, as it might break the framework.
%
% If you are unsure about the input and return values of functions you
% should read their documentation which tells you the expected dimensions.

% Turn off pagination:
more off;

clear all;
close all;

% Make tools available
addpath('tools');

% Read world data, i.e. landmarks. The true landmark positions are not given to the robot
landmarks = read_world('../data/world.dat');
% Read sensor readings, i.e. odometry and range-bearing sensor
data = read_data('../data/sensor_data.dat');

% Get the number of landmarks in the map
N = size(landmarks,2);

noise = [0.005, 0.01, 0.005]';

% how many particles
numParticles = 100;

% initialize the particles array
particles = struct;
for i = 1:numParticles
  particles(i).weight = 1. / numParticles;
  particles(i).pose = zeros(3,1);
  particles(i).history = cell();
  for l = 1:N % initialize the landmarks aka the map
    particles(i).landmarks(l).observed = false;
    particles(i).landmarks(l).mu = zeros(2,1);    % 2D position of the landmark
    particles(i).landmarks(l).sigma = zeros(2,2); % covariance of the landmark
  end
end

% toogle the visualization type
%showGui = true;  % show a window while the algorithm runs
showGui = false; % plot to files instead

% Perform filter update for each odometry-observation pair read from the
% data file.
for t = 1:size(data.timestep, 2)
%for t = 1:50
    printf('timestep = %d\n', t);

    % Perform the prediction step of the particle filter
    particles = prediction_step(particles, data.timestep(t).odometry, noise);

    % Perform the correction step of the particle filter
    particles = correction_step(particles, data.timestep(t).sensor);

    % Generate visualization plots of the current state of the filter
    plot_state(particles, landmarks, t, data.timestep(t).sensor, showGui);

    % Resample the particle set
    particles = resample(particles);
end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function rgb = ind2color(ind,varargin)
% convert ind to rgb values

cmap = getMTEXpref('colors');

if check_option(varargin,'ordered')
  order = [1,6,19,12,20,22,3,11,18,16,...
    23,8,7,2,14,15,9,5,21,4,24,13,10,17];

  iorder(order) = 1:24;

  cmap = cmap(order,:);

elseif nargin > 1 && isa(varargin{1},'double')

  cmap = varargin{1};

end

isInd = ind > 0 & ~isnan(ind);
rgb = nan(numel(ind),3);

rgb(isInd,:) = cmap(ind(isInd),:);

author: mtex-toolbox

repo: mtex

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


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha
    [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters);
end

author: SaveTheRbtz

repo: ml-class


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2011 Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function PI = b2nst__ (space, knots, degree, msh)

  U = knots{1};
  V = knots{2};
  dU = degree(1);
  dV = degree(2);

  %% FIXME only odd degree supported at the moment
  if (mod (dU, 2) == 0 || mod (dV, 2) == 0)
   error ('b2nst__: only odd degree supported at the moment')
  end

  mcp = numel (U) - dU - 1;
  ncp = numel (V) - dV - 1;
  UT = U;
  VT = V;
  UT([dU+2 end-dV-1]) = [];
  VT([dV+2 end-dV-1]) = [];

  %% get numbers of shape functions to replace
  ucenter = (dU+1)/2+1;
  vcenter = (dV+1)/2+1;
  aux = setdiff (1:dV+2, vcenter);
  ns = 0;
  for jj = aux
    ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], ucenter, jj);
    ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], mcp+1-ucenter, jj);
  end

  for jj = ncp+1-aux;
    ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], ucenter, jj);
    ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], mcp+1-ucenter, jj);
  end

  %% get numbers of shape functions to remove
  ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], ucenter, vcenter);
  ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], ucenter, ncp+1-vcenter);
  ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], mcp+1-ucenter, ncp+1-vcenter);
  ns=ns+1; rmshp(ns) = sub2ind ([mcp, ncp], mcp+1-ucenter, vcenter);

  nrmshp = numel (rmshp);

  %% construct projection matrix
  P = speye (space.ndof);
  %% remove basis functions to be removed
  P(:, rmshp) = [];

 % uv = [reshape(msh.quad_nodes(1, :, :),[],1), reshape(msh.quad_nodes(2, :, :),[],1)]';
  ns = 0;
  for jj = 1:dV+1
    ns=ns+1; fun{ns} = @(pts) tbasisfun (pts, [dU, dV], ...
                 {U(ucenter+[0:dU+1]), VT(jj+[0:dV+1])});
    ns=ns+1; fun{ns} = @(pts) tbasisfun (pts, [dU, dV], ...
                 {U(ucenter+[0:dU+1]), VT(end+1-jj-[dV+1:-1:0])});
    ns=ns+1; fun{ns} = @(pts) tbasisfun (pts, [dU, dV], ...
                 {U(end-ucenter+1-[dU+1:-1:0]), VT(jj+[0:dV+1])});
    ns=ns+1; fun{ns} = @(pts) tbasisfun (pts, [dU, dV], ...
                 {U(end-ucenter+1-[dU+1:-1:0]), VT(end+1-jj-[dV+1:-1:0])});
  end

% Computation of the projection
%  We cannot use the op_f_v_tp operator, because our function handles fun{jj}
%  are for points in the parametric domain, not in the physical domain
  M = op_u_v_tp (space, space, msh);

  rhs = zeros (space.ndof, ns);
  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    sp_col  = sp_evaluate_col (space, msh_col, 'gradient', false);
    uv = [reshape(msh_col.quad_nodes(1,:,:),[],1), reshape(msh_col.quad_nodes(2,:,:),[],1)]';
    for jj = 1:ns
      rhs(:,jj) = rhs(:,jj) + op_f_v (sp_col, msh_col, ...
                               reshape (fun{jj}(uv), msh_col.nqn, msh_col.nel));
    end
  end

  for jj = 1:ns
    PI(:, jj) = M \ rhs(:,jj);
%    PI(:, jj) = M \ op_f_v_tp (space, msh, reshape (fun(:,jj), size (msh.quad_weights)));
  end

  PI = [P PI];

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [KM] = calc(kern,dat1,dat2,ind1,ind2)

% [KM] = calc(k,dat1,[d2],[ind1],[ind2])
%
% Calculate the kernel with specific indexed data from
% datasets d1 and (optionally) d2.
% If d2 is not given it is assumed to make the matrix between d1 and itself

  global caching_kernel;

  if nargin<3,
      dat2=[];
  end;

  if isempty(dat2)
      dat2=dat1;
  end;

  if nargin<4
    ind1=[1:length(get_index(dat1))];
  end

  if nargin<5
    ind2=[1:length(get_index(dat2))];
  end

  if kern.calc_on_output
      dat1.X=dat1.Y;
      dat2.X=dat2.Y;
  end %% output kernels

  index1 = get_index(dat1);
  index2 = get_index(dat2);

  if strcmp(kern.ker,'custom')
   % the kernel matrix is stored in param{1} of kernel
   KM = kern.kerparam;
   KM = KM(index2(ind2),index1(ind1));
   return;
  end

  if strcmp(kern.ker,'custom_fast')
   global K;
   KM = K(index2(ind2),index1(ind1));
   return;
 end

 if kern.kercaching==1 & (~isa(dat1,'data_global')|~isa(dat2,'data_global')),
   error('Caching kernel only for data_global');
 else
   if (kern.kercaching==1)&~ismyX(dat1)&~ismyX(dat2),
     if caching_kernel==kern,
       global K;
       index1 = get_index(dat1);
       index2 = get_index(dat2);
       KM = K(index2(ind2),index1(ind1));
       return;
     else
       global K;
       %global X;global Y;
       dtmp=data_global('get_kernel tmp',[],[]);
       [numEx vDim oDim]  = get_dim(dtmp);
       ind1=[1:numEx]; ind2=[1:numEx];
       caching_kernel=kern;
     end
   end
 end

  paramtmp=[];
  if iscell(kern.kerparam),
    for i=1:length(kern.kerparam),
      paramtmp{i} = kern.kerparam{i};
    end
  else
    paramtmp = kern.kerparam;
  end
  if kern.kercaching~=1|ismyX(dat1)|ismyX(dat2),

    %% --- normal Kernel calculation ----------------------
    KM = feval(kern.ker,kern,dat1,dat2,ind1,ind2,paramtmp);

  else
    K = feval(kern.ker,kern,dtmp,dtmp,ind1,ind2,paramtmp);
    KM = K(get_index(dat2),get_index(dat1));
  end;

  if kern.output_distance  %% ------------ convert kernel to a distance
    kern.output_distance=0;
    Kdn = get_norm(kern,dat1,ind1).^2;
    Kn = get_norm(kern,dat2,ind2).^2;
    KM = ones(length(Kn),1)*Kdn' + Kn*ones(1,length(Kdn)) - 2*KM;
    KM = sqrt(KM);
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ber = get_ber(V_equal, W_equal)

global Ns Nsym Nr Vn H hMod hDemod;
%one channel realization for Nsym times data streams transmission
% 2 is for real and imaginary
data = randi([0 1],Nsym*Ns*2,1);

%QPSK modulation generate all original signals s in Nsym times
s = reshape(step(hMod,data),Ns, Nsym);

%generate noise vector u
u = sqrt(Vn/2).*(randn(Nr,Nsym)+1i*randn(Nr,Nsym));

%get the receive vector
%colloct Nsym receive vector in the r matrix, where nth column is a
%receive vector at n time
r = W_equal' * H * V_equal * s + W_equal' * u;

%QPSK demodulation
de_data = step(hDemod, r(:));
%get the number of error bits
ber = biterr(data,de_data);

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This files deals the parameters for the Testing environment
%
% options:
% FAIRtestPara('disp')          displays the current settings
% FAIRtestPara('set',list{:})	sets defaults and overwrites with user supplied values
% FAIRtestPara('get',name)		gets value of variable name
% FAIRtestPara('clear')			clears all variables
%
%==============================================================================

function varargout = FAIRtestPara(varargin)

caller = dbstack; % identify the name of the calling function

if nargin == 0,
  help(mfilename);
  FAIRtestPara('disp');
  return;
end;

if nargin > 0 && isstr(varargin{1}),
  task = varargin{1};
  varargin(1) = [];

  switch task,
    case 'set',
      setFAIRtestpara(varargin{:});

    case 'get',
      for j = 1:length(varargin)
        varargout{j} = getappdata(0,varargin{2*j-1});
      end;

%    case 'value',
%      varargout{1} = strcmp(getappdata(0,varargin{1}),'on');

    case 'disp',
      dispFAIRtestpara;

    case 'clear',
      FAIRpara = getappdata(0);
      Fnames   = fieldnames(FAIRpara);
      K        = reshape(find(cellfun(@isempty,strfind(Fnames,'FAIR'))==0),1,[]);
      for k = K
        rmappdata(0,Fnames{k});
      end;
    otherwise, error('1');
  end;
else
  error(1)
end;

%------------------------------------------------------------------------------

function setFAIRtestpara(varargin)

para = {
  'FAIRcaller',         ''
  'FAIRfolder',         ''
  'FAIRedit',           'on'
  'FAIRcompile',        'on'
  'FAIRignoreCompiles', 'on'
  'FAIRrun',            'on'
  'FAIRrange',          'full'
  'FAIRextension',      'all'
  'FAIRerror',          'on'
  'FAIRkeyboard',       'on'
  'FAIRfilter',         'off'
  'FAIRpause',          'off'
  'FAIRinput',          'off'
  'FAIRclearTemp',      'on'
  'FAIRcheckList',      []
};

for j=1:size(para,1)
  value = getappdata(0,para{j,1});
  if isempty(value)
    setappdata(0,para{j,1},para{j,2})
  end;
end;

for k=1:length(varargin)/2, % overwrites default parameter
    setappdata(0,varargin{2*k-1},varargin{2*k});
end;

%------------------------------------------------------------------------------

function dispFAIRtestpara;
FAIRpara = getappdata(0);
Fnames   = fieldnames(FAIRpara);
K        = reshape(find(cellfun(@isempty,strfind(Fnames,'FAIR'))==0),1,[]);

if isempty(K)
return;
end;

FAIRmessage('.');
for k = K
  name  = Fnames{k};
  value = getappdata(0,name);
  if ischar(value)
    fprintf('%-20s: %s\n',name,value)
  elseif isnumeric(value) && (length(value) == 1)
    fprintf('%-20s= %s\n',name,num2str(value))
  elseif isnumeric(value),
    fprintf('%-20s: %d files\n',name,length(value))
  elseif isstruct(value),
    fprintf('%-20s: %d files\n',['struct ',name],length(value))
  else
    keyboard
  end;
end;
FAIRmessage('.');

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


text:
function [ xy, w ] = triangle_nco_rule ( rule, order_num )

%*****************************************************************************80
%
%% TRIANGLE_NCO_RULE returns the points and weights of an NCO rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Input, integer ORDER_NUM, the order (number of points) of the rule.
%
%    Output, real XY(2,ORDER_NUM), the points of the rule.
%
%    Output, real W(ORDER_NUM), the weights of the rule.
%

%
%  Get the suborder information.
%
  suborder_num = triangle_nco_suborder_num ( rule );

  suborder = triangle_nco_suborder ( rule, suborder_num );

  [ suborder_xyz, suborder_w ] = triangle_nco_subrule ( rule, suborder_num );
%
%  Expand the suborder information to a full order rule.
%
  o = 0;

  for s = 1 : suborder_num

    if ( suborder(s) == 1 )

      o = o + 1;
      xy(1:2,o) = suborder_xyz(1:2,s);
      w(o) = suborder_w(s);

    elseif ( suborder(s) == 3 )

      for k = 1 : 3
        o = o + 1;
        xy(1,o) = suborder_xyz ( i4_wrap(k,  1,3), s );
        xy(2,o) = suborder_xyz ( i4_wrap(k+1,1,3), s );
        w(o) = suborder_w(s);
      end

    elseif ( suborder(s) == 6 )

      for k = 1 : 3
        o = o + 1;
        xy(1,o) = suborder_xyz ( i4_wrap(k,  1,3), s );
        xy(2,o) = suborder_xyz ( i4_wrap(k+1,1,3), s );
        w(o) = suborder_w(s);
      end

      for k = 1 : 3
        o = o + 1;
        xy(1,o) = suborder_xyz ( i4_wrap(k+1,1,3), s );
        xy(2,o) = suborder_xyz ( i4_wrap(k,  1,3), s );
        w(o) = suborder_w(s);
      end

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'TRIANGLE_NCO_RULE - Fatal error!\n' );
      fprintf ( 1, '  Illegal SUBORDER(%d) = %d\n', s, suborder(s) );
      error ( 'TRIANGLE_NCO_RULE - Fatal error!' );

    end

  end

  return
end


text:
function [data_out, indd] = array_padd(data_in, padsize, paddvalue, direction, paddmode)

% function to pad data array with various border conditions
% INPUTS:
%   DATA_IN - input data array
%   PADSIZE - same as in padarray.m, [rowpad, colpad] number of samples to
%             pad in row and in column direction
%   PADDVALUE - numerical value used to pad (ignored within some pad modes)
%   DIRECTION - same as in padarray.m {'both' 'post' and 'pre'}
%   PADDMODE  - 'circular', 'replicate' and 'symmetric' are the  same as in
%             padarray.m.  New options are:
%    'barthannwin', 'bartlett', 'blackman','blackmanharris', 'bohmanwin',
%    'flattopwin','gausswin','hamming' ,'hann', 'nuttallwin','parzenwin', 'triang'
%     In it, 'symmetric' padded values are multiplied by the 1/2 of corresponding
%     window to taper off the values to zero. Symmetrical padding option in
%     it can be replaced by 'replicate' by uncommenting line 70 and
%     commenting line 71 in function body. With these windowing options,
%     direction and paddvalue is ignored (but must be present for consistency)
%     and internaly 'symmetric' and  'both' are used
%
%
% OUTPUTS:  DATA_OUT    - output padded data array
%           INDD        - indeces of padded array used to recover the original array
%
%EXAMPLES:
%       data_in = [1 1 1 1 1; 1 2 3 2 1; 1 2 3 2 1; 1 1 1 1 1]
%       [data_out, indd] = array_padd(data_in, [3, 5])
%       [data_out, indd] = array_padd(data_in, [3, 5], 5)
%       [data_out, indd] = array_padd(data_in, [3, 5], 0, 'both')
%       [data_out, indd] = array_padd(data_in, [3, 5], 0, 'both', 'replicate')
%       [data_out, indd] = array_padd(data_in, [3, 5], 0, 'both', 'symmetric')
%       [data_out, indd] = array_padd(data_in, [3, 5], 0, 'both', 'hamming')
%       imagesc(data_out); colorbar
% original array size and position within padded array can be recovered as
%       data_out = data_out(indd(1):indd(2),indd(3):indd(4));

% Other m-files required: none
% Subfunctions: none
% MAT-files required:
% window.m and padarray.m from Signal and Image Processing Toolboxes

%____________________________________________
% 	Sergei Koptenko, Resonant Medical Inc.,
%           Montreal, Qc., Canada
%	sergei.koptenko@resonantmedical.com
%   Website: http://www.resonantmedical.com
%____________Feb/30/2005_____________________

if nargin <2,
	disp('Not enough arguments')
    elseif  nargin <3, paddvalue =0; direction = 'both'; paddmode = 'simple';
        elseif  nargin <4, direction = 'both'; paddmode = 'simple';
            elseif nargin <5, paddmode = 'simple';
end

[rrow,ccol] = size(data_in);
%Find indices of the original array within the padded array
switch direction
    case {'both','pre' }
       indd = [padsize(1)+1,(padsize(1)+rrow), padsize(2)+1,(padsize(2)+ccol)];
    case 'post'
        indd = [1,rrow, 1,ccol];
end

% Create the padded array
switch paddmode
    case 'simple'
          data_out = padarray(data_in, padsize, paddvalue, direction);

   case {'circular', 'replicate' , 'symmetric' }
        data_out = padarray(data_in,padsize, paddmode, direction);

    case  {'barthannwin', 'bartlett', 'blackman', 'blackmanharris',...
            'bohmanwin', 'flattopwin','gausswin', 'hann' , 'nuttallwin',...
            'parzenwin' , 'triang','hamming'} % This option forces direction == 'both'
        eval(['rowwind =window(@' paddmode ', ' num2str(2*padsize(1)) ');']);
        eval(['colwind =window(@' paddmode ', ' num2str(2*padsize(2)) ');']);

%       data_out = padarray(data_in, padsize, 'replicate', 'both');
        data_out = padarray(data_in, padsize, 'symmetric', 'both');
        [mrow, mcol] = size(data_out);
%________Create a Column mask______________
        tc =ones(mrow,1); % single column mask
        tc(1:indd(1)-1) = tc(1:indd(1)-1) .* rowwind(1:padsize(1));  % mask for the row start
        tc(indd(2)+1:end) = tc(indd(2)+1:end) .* rowwind(padsize(1)+1:end);%mask for the row end
        mc = repmat(tc, 1, mcol); %column mask array
%________Create a ROW mask______________
        tr =ones(1, mcol); % single ROW mask
        tr(1:indd(3)-1) = tr(1:indd(3)-1) .* colwind(1:padsize(2))';    % mask for the col start
        tr(indd(4)+1:end) = tr(indd(4)+1:end) .* colwind(padsize(2)+1:end)'; % mask for the col end
        mr = repmat(tr, mrow, 1); %  row mask array
%________Create FULL mask______________
        data_out = data_out .*  mr .* mc;

   otherwise
      disp('Unknown method.')
end
return

author: Sable

repo: mcbench-benchmarks


text:
function [sbatch2 sbatch3] = select_neighbors(sbatch1, data)

l = double(data.Ltraining(sbatch1)) + 1;
nnn1 = size(data.indnlabels,1);
sbatch2 = double(data.indnlabels(ceil(rand(numel(sbatch1),1).*(data.nnlabels(l))) + (l-1)*nnn1));

nnn2 = size(data.nnTraining,1);
sbatch3 = double(data.nnTraining(ceil(rand(numel(sbatch1),1).*(data.nns(sbatch1))) + (sbatch1-1)*nnn2));

author: norouzi

repo: hdml

save_path: github-repos/MATLAB/norouzi-hdml


text:
function vistformfwd(tform, wdata, zdata, N)
%VISTFORMFWD Visualize forward geometric transform.
%   VISTFORMFWD(TFORM, WRANGE, ZRANGE, N) shows two plots: an N-by-N
%   grid in the W-Z coordinate system, and the spatially transformed
%   grid in the X-Y coordinate system.  WRANGE and ZRANGE are
%   two-element vectors specifying the desired range for the grid.  N
%   can be omitted, in which case the default value is 10.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/04/18 05:07:34 $

if nargin < 4
   N = 10;
end

% Create the w-z grid and transform it.
[w, z] = meshgrid(linspace(wdata(1), zdata(2), N), ...
                  linspace(wdata(1), zdata(2), N));
wz = [w(:) z(:)];
xy = tformfwd([w(:) z(:)], tform);

% Calculate the minimum and maximum values of w and x,
% as well as z and y. These are used so the two plots can be
% displayed using the same scale.
x = reshape(xy(:, 1), size(w)); % reshape is discussed in Sec. 8.2.2.
y = reshape(xy(:, 2), size(z));
wx = [w(:); x(:)];
wxlimits = [min(wx) max(wx)];
zy = [z(:); y(:)];
zylimits = [min(zy) max(zy)];

% Create the w-z plot.
subplot(1,2,1) % See Section 7.2.1 for a discussion of this function.
plot(w, z, 'b'), axis equal, axis ij
hold on
plot(w', z', 'b')
hold off
xlim(wxlimits)
ylim(zylimits)
set(gca, 'XAxisLocation', 'top')
xlabel('w'), ylabel('z')

% Create the x-y plot.
subplot(1, 2, 2)
plot(x, y, 'b'), axis equal, axis ij
hold on
plot(x', y', 'b')
hold off
xlim(wxlimits)
ylim(zylimits)
set(gca, 'XAxisLocation', 'top')
xlabel('x'), ylabel('y')

repo: weiyanmin


text:
% Steady state transport solutions for one Pe=0 and various values of Da2
x = [0:0.01:1];
Da2 = 8; mu1 = sqrt(Da2); mu2 = -mu1;
s = mu2*exp(mu2)-mu1*exp(mu1);
c = (mu2*exp(mu2)*exp(mu1*x)-mu1*exp(mu1)*exp(mu2*x))./s
for Da2 = [4 2 1 0.5 0.25 0.125];
    s = sqrt(Da2); mu1 = s; mu2 = -s;
    s = mu2*exp(mu2)-mu1*exp(mu1);
    c = [c;(mu2*exp(mu2)*exp(mu1*x)-mu1*exp(mu1)*exp(mu2*x))./s];
end
plot (x,c);
legend('Da_2=8','Da_2=4','Da_2=2','Da_2=1','Da_2=0.5','Da_2=0.25','Da_2=0.125');
xlabel ('x/L [-]'); ylabel ('c/c_{in} [-]');


text:
function M1 = image_resize(M,p1,q1,r1)

% image_resize - resize an image using bicubic interpolation
%
%   M1 = image_resize(M,nx,ny,nz);
% or
%   M1 = image_resize(M,newsize);
%
%   Works for 2D, 2D 2 or 3 channels, 3D images.
%
%   Copyright (c) 2004 Gabriel Peyre

if nargin==2
    % size specified as an array
    q1 = p1(2);
    if length(p1)>2
        r1 = p1(3);
    else
        r1 = size(M,3);
    end
    p1 = p1(1);
end

if nargin<4
    r1 = size(M,3);
end

if ndims(M)>=3
    % RVB image
    M1 = zeros(p1,q1, size(M,3), size(M,4));
    for i=1:size(M,3)
        for j=1:size(M,4)
            M1(:,:,i,j) = image_resize(M(:,:,i,j), p1, q1);
        end
    end
    return;
% elseif ndims(M)==3
%     p = size(M,1);
%     q = size(M,2);
%     r = size(M,3);
%     [Y,X,Z] = meshgrid( (0:q-1)/(q-1), (0:p-1)/(p-1), (0:r-1)/(r-1)  );
%     [YI,XI,ZI] = meshgrid( (0:q1-1)/(q1-1), (0:p1-1)/(p1-1), (0:r1-1)/(r1-1) );
%     M1 = interp3( Y,X,Z, M, YI,XI,ZI ,'cubic');
%     return;
end

p = size(M,1);
q = size(M,2);
[Y,X] = meshgrid( (0:q-1)/(q-1), (0:p-1)/(p-1) );
[YI,XI] = meshgrid( (0:q1-1)/(q1-1), (0:p1-1)/(p1-1) );
M1 = interp2( Y,X, M, YI,XI ,'cubic');

author: optimaltransport


text:
function [lat,lon] = maptrimp(lat0,lon0,latlim,lonlim)

    %MAPTRIMP  Trims a patch map to a specified region
    %
    %  [lat,lon] = MAPTRIMP(lat0,lon0,latlim,lonlim) trims a patch map
    %  to a region specified by latlim and lonlim.  Latlim and lonlim
    %  are two element vectors, defining the latitude and longitude limits
    %  respectively.  The inputs lat0 and lon0 must be vectors representing
    %  patch map vector data.
    %
    %  See also MAPTRIMS, MAPTRIML

    %  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
    %  Written by:  E. Byrns, E. Brown
    %   $Revision: 1399 $    $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    report_this_filefun(mfilename('fullpath'));

    if nargin < 4;   error('Incorrect number of arguments');   end

    %  Test the inputs

    if  ~isequal(sort(size(latlim)),sort(size(lonlim)),[1 2])
        error('Lat and lon limit inputs must be 2 element vectors')
    end

    %  Test for real inputs

    if any([~isreal(lat0) ~isreal(lon0) ~isreal(latlim) ~isreal(lonlim)])
        warning('Imaginary parts of complex arguments ignored')
        lat0 = real(lat0);       lon0 = real(lon0);
        latlim = real(latlim);   lonlim = real(lonlim);
    end

    %  Get the corners of the submap region

    up    = max(latlim);   low  = min(latlim);
    right = max(lonlim);   left = min(lonlim);

    %  Copy the input data and ensure column vectors.

    lat = lat0(:);   lon = lon0(:);

    %  Get the vector of patch items and remove any NaN padding
    %  at the beginning or end of the column.  This eliminates potential
    %  multiple NaNs at the beginning and end of the patches.

    while isnan(lat(1)) | isnan(lon(1))
        lat(1) = [];   lon(1) = [];
    end
    while isnan(lat(length(lat))) | isnan(lon(length(lon)));
        lat(length(lat)) = [];   lon(length(lon)) = [];
    end

    %  Add a NaN to the end of the data vector.  Necessary for processing
    %  of multiple patches.

    lat(length(lat)+1) = NaN;   lon(length(lon)+1) = NaN;

    %  Find the individual patches and then trim the data

    indx = find(isnan(lon) | isnan(lat));
    if isempty(indx);   indx = length(lon)+1;   end

    for i = 1:length(indx)

        if i == 1;      startloc = 1;
        else;       startloc = indx(i-1)+1;
        end
        endloc   = indx(i)-1;

        indices = (startloc:endloc)';   %  Indices will be empty if NaNs are
        %  neighboring in the vector data.
        if ~isempty(indices)            %  Should not happen, but test just in case

            %  Patches which lie completely outside the trim window.  Replace
            %  with NaNs and then eliminate it entirely later.  Replacing with
            %  NaNs is useful so that the indexing with indices is not messed
            %  up if an entire patch is eliminated at this point.

            %  If at least one point of the patch edge does not lie with the
            %  specified window limits, then the entire patch is trimmed.

            if ~any(lon(indices) >= left & lon(indices) <= right  &&  ...
                    lat(indices) >= low  & lat(indices) <= up)
                lon(indices) = NaN;    lat(indices) = NaN;
            end

            %  Need to only test along edge since patch must lie somehow within the window.

            %  Points off the bottom

            loctn = find( lon(indices) < left );
            if ~isempty(loctn);  lon(indices(loctn)) = left;     end

            %  Points off the top

            loctn = find( lon(indices) > right );
            if ~isempty(loctn);   lon(indices(loctn)) = right;   end

            %  Points off the left

            loctn = find( lat(indices) < low );
            if ~isempty(loctn);   lat(indices(loctn)) = low;     end

            %  Points off the right

            loctn = find( lat(indices) > up );
            if ~isempty(loctn);   lat(indices(loctn)) = up;      end
        end
    end

    %  Eliminate multiple NaNs in the vector.  Will occur if a patch
    %  lies entirely outside the window of interest.

    if ~isempty(lat)
        nanloc = isnan(lat);	[r,c] = size(nanloc);
        nanloc = find(nanloc(1:r-1,:) & nanloc(2:r,:));
        lat(nanloc) = [];  lon(nanloc) = [];
    end

author: CelsoReyes


text:
function [c, ceq] = collectConstraints(t,x,u,defects, pathCst, bndCst)
% [c, ceq] = collectConstraints(t,x,u,defects, pathCst, bndCst)
%
% OptimTraj utility function.
%
% Collects the defects, calls user-defined constraints, and then packs
% everything up into a form that is good for fmincon.
%
% INPUTS:
%   t = time vector
%   x = state matrix
%   u = control matrix
%   defects = defects matrix
%   pathCst = user-defined path constraint function
%   bndCst = user-defined boundary constraint function
%
% OUTPUTS:
%   c = inequality constraint for fmincon
%   ceq = equality constraint for fmincon
%

ceq_dyn = reshape(defects,numel(defects),1);

%%%% Compute the user-defined constraints:
if isempty(pathCst)
    c_path = [];
    ceq_path = [];
else
    [c_pathRaw, ceq_pathRaw] = pathCst(t,x,u);
    c_path = reshape(c_pathRaw,numel(c_pathRaw),1);
    ceq_path = reshape(ceq_pathRaw,numel(ceq_pathRaw),1);
end
if isempty(bndCst)
    c_bnd = [];
    ceq_bnd = [];
else
    t0 = t(1);
    tF = t(end);
    x0 = x(:,1);
    xF = x(:,end);
    [c_bnd, ceq_bnd] = bndCst(t0,x0,tF,xF);
end

%%%% Pack everything up:
c = [c_path;c_bnd];
ceq = [ceq_dyn; ceq_path; ceq_bnd];

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [filtered_signal,filtb,filta]=hipass_butterworth(inputsignal,cutoff_freq,Fs,order)
% High-pass Butterworth filter
% [filtered_signal,filtb,filta] = hipass_butterworth(inputsignal,cutoff_freq,Fs,order)
%
% This is simply a set of built-in Matlab functions, repackaged for ease of
% use by Chad Greene, October 2012.
%
% INPUTS:
% inputsignal = input time series
% cutoff_freq = filter corner frequency
% Fs = data sampling frequency
% order = order of Butterworth filter
%
% OUTPUTS:
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% y_filt = hipass_butterworth(y,900,Fs,4); % cut off at 900 Hz
% figure
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% legend('unfiltered','filtered')
% sound(y,Fs)      % play original time series
% pause(2)         % pause two seconds
% sound(y_filt,Fs) % play filtered time series
%
%
% EXAMPLE 2:
% load train
% t = (1:length(y))/Fs;
% [y_filt,filtb,filta] = hipass_butterworth(y,900,Fs,4); % cut off at 900 Hz
% [h1,f1] = freqz(filtb,filta,256,Fs);
%
% figure
% subplot(3,1,1)
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% text(0,.1,' time series','units','normalized')
%
% subplot(3,1,2)
% AX = plotyy(f1,10*log10(abs(h1)),f1,angle(h1),'semilogx');
% set(get(AX(1),'ylabel'),'string','gain (dB)')
% set(get(AX(2),'ylabel'),'string','phase (rad)')
% xlim(AX(1),[min(f1) max(f1)])
% xlim(AX(2),[min(f1) max(f1)])
% text(0,.1,' filter response','units','normalized')
% box off
%
% [Pxx,f] = pwelch(y,512,256,[],Fs,'onesided');
% [Pxxf,f_f]= pwelch(y_filt,512,256,[],Fs,'onesided');
% subplot(3,1,3)
% semilogx(f,10*log10(Pxx))
% hold on
% semilogx(f_f,10*log10(Pxxf),'r')
% xlabel('frequency (Hz)')
% ylabel('PSD (dB)')
% xlim([min(f1) max(f1)])
% box off
% legend('unfiltered','filtered','location','northwest')
% legend boxoff

nyquist_freq = Fs/2;  % Nyquist frequency
Wn=cutoff_freq/nyquist_freq;    % non-dimensional frequency
[filtb,filta]=butter(order,Wn,'high'); % construct the filter
filtered_signal=filtfilt(filtb,filta,inputsignal); % filter the data with zero phase


text:
% op_addrcvrs.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,fids_presum,specs_presum,ph,sig]=op_addrcvrs(in,point,mode,coilcombos);
%
% DESCRIPTION:
% Perform weighted coil recombination for MRS data acquired with a reciever
% coil array.
%
% INPUTS:
% in            = input spectrum in matlab structure format.
% point         = point of fid to use for phase estimation (optional. Default = 1);
% mode          = Method for estimating the coil weights and phases (optional.  Default = 'w').
%                 -'w' performs amplitude weighting of channels based on the
%                 maximum signal of each coil channel.
%                 -'h' performs amplitude weighting of channels based on the
%                 maximum signal of each coil channel divided by the square of
%                 the noise in each coil channel (as described by Hall et al.
%                 Neuroimage 2014).
% coilcombos    = The predetermined coil weights and phases and amplitudes as
%                 generated by the op_getcoilcombos.m function.  If this
%                 argument is provided, the 'point', and 'mode', arguments
%                 will be ignored.  (optional.  Default = []).
%
% OUTPUTS:
% out           = Output dataset with coil channels combined.
% fids_presum   = Input data with coil channels in phase (time domain).
% specs_presum  = Input data with coil channels in phase (frequency domain).
% ph            = Vector of applied coil phases (in degrees).
% sig           = Vector of coil weights.

function [out,fids_presum,specs_presum,ph,sig]=op_addrcvrs(in,point,mode,coilcombos);

if in.flags.addedrcvrs || ~in.dims.coils
    disp('WARNING:  Only one receiver channel found!  Returning input without modification!');
    out=in;
    out.flags.addedrcvrs=1;
    fids_presum=in.fids;
    specs_presum=in.specs;
    ph=0;
    sig=1;
else

    %To get best possible SNR, add the averages together (if it hasn't already been done):
    if ~in.flags.averaged
        av=op_averaging(in);
    else
        av=in;
    end

    %also, for best results, we will combine all subspectra:
    if nargin<4
        if in.flags.isFourSteps
            av=op_fourStepCombine(av);
        end
        if in.dims.subSpecs>0
            av=op_combinesubspecs(av,'summ');
        end
        if nargin<3
            mode = 'w';
            if nargin<2
                point=1;
            end
        end
    end
    avfids=av.fids;
    avspecs=av.specs;

    %initialize phase matrix and the amplitude maxtrix that are the size of nPoints x Coils
    %ph=ones(in.sz(in.dims.t),in.sz(in.dims.coils));
    %sig=ones(in.sz(in.dims.t),in.sz(in.dims.coils));

    %now start finding the relative phases between the channels and populate
    %the ph matrix
    for n=1:in.sz(in.dims.coils)
        if nargin<4
            %ph(:,n)=phase(avfids(point,n,1,1))*ph(:,n);
            phs(n)=phase(avfids(point,n,1,1));
            switch mode
                case 'w'
                    %sig(:,n)=abs(avfids(point,n,1,1))*sig(:,n);
                    sigs(n)=abs(avfids(point,n,1,1));
                case 'h'
                    S=max(abs(avfids(:,n,1,1)));
                    N=std(avfids(end-100:end,n,1,1));
                    %sig(:,n)=(S/(N.^2))*sig(:,n);
                    sigs(n)=(S/(N.^2));
            end
        else
            %ph(:,n)=coilcombos.ph(n)*ph(:,n);
            phs(n)=coilcombos.ph(n);
            sigs(n)=coilcombos.sig(n);
        end
    end

    %now replicate the phase matrix to equal the size of the original matrix:
    % replicate=in.sz;
    % replicate(1)=1;
    % replicate(in.dims.coils)=1;
    % ph=repmat(ph,replicate);
    % sig=repmat(sig,replicate);
    sigs=sigs/norm(sigs(:));

    ph=ones(in.sz);
    sig=ones(in.sz);

    if in.dims.coils==1
        for n=1:in.sz(1)
            ph(n,:)=phs(n)*ph(n,:);
            sig(n,:)=sigs(n)*sig(n,:);
        end
    elseif in.dims.coils==2
        for n=1:in.sz(2)
            ph(:,n,:)=phs(n)*ph(:,n,:);
            sig(:,n,:)=sigs(n)*sig(:,n,:);
        end
    elseif in.dims.coils==3
        for n=1:in.sz(3)
            ph(:,:,n,:)=phs(n)*ph(:,:,n,:);
            sig(:,:,n,:)=sigs(n)*sig(:,:,n,:);
        end
    elseif in.dims.coils==4
        for n=1:in.sz(4)
            ph(:,:,:,n,:)=phs(n)*ph(:,:,:,n,:);
            sig(:,:,:,n,:)=sigs(n)*sig(:,:,:,n,:);
        end
    elseif in.dims.coils==5
        for n=1:in.sz(5)
            ph(:,:,:,:,n)=phs(n)*ph(:,:,:,:,n);
            sig(:,:,:,:,n)=sigs(n)*sig(:,:,:,:,n);
        end
    end

    %now apply the phases by multiplying the data by exp(-i*ph);
    fids=in.fids.*exp(-i*ph);
    fids_presum=fids;
    specs_presum=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

    %Apply the amplitude factors by multiplying the data by amp;
    if mode=='w' || mode=='h'
        fids=fids.*sig;
    end

    %now sum along coils dimension
    fids=sum(fids,in.dims.coils);
    fids=squeeze(fids);

    %re-calculate Specs using fft
    specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

    %change the dims variables
    if in.dims.t>in.dims.coils
        dims.t=in.dims.t-1;
    else
        dims.t=in.dims.t;
    end
    dims.coils=0;
    if in.dims.averages>in.dims.coils
        dims.averages=in.dims.averages-1;
    else
        dims.averages=in.dims.averages;
    end
    if in.dims.subSpecs>in.dims.coils
        dims.subSpecs=in.dims.subSpecs-1;
    else
        dims.subSpecs=in.dims.subSpecs;
    end
    if in.dims.extras>in.dims.coils
        dims.extras=in.dims.extras-1;
    else
        dims.extras=in.dims.extras;
    end

    %re-calculate the sz variable
    sz=size(fids);

    %FILLING IN DATA STRUCTURE
    out=in;
    out.fids=fids;
    out.specs=specs;
    out.sz=sz;
    out.dims=dims;

    %FILLING IN THE FLAGS
    out.flags=in.flags;
    out.flags.writtentostruct=1;
    out.flags.addedrcvrs=1;

end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function F = cscd(F, varargin)
%CSCD   Cosecant of a CHEBFUN, result in degrees.
%   CSCD(F) computes the cosecant (in degrees) of the CHEBFUN F.
%
%   CSCD(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also ACSCD, CSC.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Call the compose method:
F = compose(F, @cscd, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
% function varargout = viewIP(I,omega,m,varargin)
%
% visualizes the intensity projections of a 3D image
%
% Input:
%   I           discretized image
%   omega		describing the domain
%	m 			number of discretization points
%   varargin    optional parameters like {'axis','off'}
%
% Output:
%  ih			image handle
%  I			[I1,I2;I3',0]
%==============================================================================

function varargout = viewIP(I,omega,m,varargin)

if nargin==0
help(mfilename);
    runMinimalExample;
	return;
end

mode = 'max';
% reshape and form the average images
I  = reshape(I,m);
switch mode,
  case 'average',
    I1 = squeeze(sum(I,3))/m(3);
    I2 = squeeze(sum(I,2))/m(2);
    I3 = squeeze(sum(I,1))/m(1);
  case 'max',
    I1 = squeeze(max(I,[],3));
    I2 = squeeze(max(I,[],2));
    I3 = squeeze(max(I,[],1));
  case 'min',
    I1 = squeeze(min(I,[],3));
    I2 = squeeze(min(I,[],2));
    I3 = squeeze(min(I,[],1));
  otherwise, eror('nyi');
end;

I   = [I1,I2;I3',255+zeros(m(3),m(3))]; cla;
ih  = imagesc(I); axis image; hold on
ph1 = plot(0.5+m(2)*[1,1],0.5+[0,m(1)+m(3)],'b-','linewidth',2);
ph2 = plot(0.5+[0,m(2)+m(3)],0.5+m(1)*[1,1],'b-','linewidth',2);
ih  = [ih;ph1;ph2];

% the following lines add some nice stuff to the code.
% if varargin = {'title','FAIR','xlabel','x'}
% the code evaluates "title('FAIR');xlabel('x');"
for k=1:2:length(varargin),
  if not(isempty(varargin{k})),
    feval(varargin{k},varargin{k+1});
  end;
end;

if nargout == 1,
  varargout = {ih};
elseif nargout ==2,
  varargout = {ih,I};
end;

%------------------------------------------------------------------------------

function runMinimalExample
load mice3D;
viewIP(dataT,omega,m);
%==============================================================================


text:
function F = minus (F,w)
%MINUS downdate a dense Cholesky factorization
%
% Example
%   % F becomes the Cholesky factorization of A-w*w'
%   F = factorize (A) ;
%   w = rand (size (A,1),1) ;
%   F = F - w ;
%   x = F\b ;               % computes x = (A-w*w')\b
%
% See also factorize, cholupdate.

% Copyright 2009, Timothy A. Davis, University of Florida

if (F.kind ~= 6)
    error ('Only dense Cholesky factorization downdate supported.') ;
end

F.R = cholupdate (F.R, w, '-') ;
F.A = F.A - w*w' ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function q = wrap2 ( m, q )

%*****************************************************************************80
%
%% WRAP2 is a circular wrap of the pointer offset Q.
%
%  Discussion:
%
%    Input values of Q between 0 and M are 'legal'.
%    Values of Q below 0 are incremented by M + 1 until they are legal.
%    Values of Q above M are decremented by M + 1 until they become legal.
%    The legal value is the output value of the function.
%
%  Example:
%
%    M  Qin  Qout
%
%    3  -5   3
%    3  -4   0
%    3  -3   1
%    3  -2   2
%    3  -1   3
%    3   0   0
%    3   1   1
%    3   2   2
%    3   3   3
%    3   4   0
%    3   5   1
%    3   6   2
%    3   7   3
%    3   8   0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2010
%
%  Author:
%
%    Original C version by Sophocles Orfanidis.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sophocles Orfanidis,
%    Introduction to Signal Processing,
%    Prentice-Hall, 1995,
%    ISBN: 0-13-209172-0,
%    LC: TK5102.5.O246.
%
%  Parameters:
%
%    Input, integer M, the maximum acceptable value for outputs.
%    M must be at least 0.
%
%    Input, integer Q, the value to be wrapped.
%
%    Output, integer Q, the wrapped value.
%
  if ( m < 0 )
    fprintf ( 2, '\n' );
    fprintf ( 2, 'WRAP2 - Fatal error!\n' );
    fprintf ( 2, '  M < 0.\n' );
    error ( 'WRAP2 - Fatal error!' );
  end
%
%  When Q = M + 1, it wraps to Q = 0.
%
  while ( m < q )
    q = q - m - 1;
  end
%
%  When Q = - 1, it wraps to Q = M.
%
  while ( q < 0 )
    q = q + m + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = bab_eigen_right ( n, alpha, beta )

%*****************************************************************************80
%
%% BAB_EIGEN_RIGHT returns the right eigenvectors of the BAB matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real ALPHA, BETA, the parameters.
%
%    Output, real A(N,N), the right eigenvector matrix.
%
  a = zeros ( n, n );

  for i = 1: n
    for j = 1 : n
      angle = ( i * j ) * pi / ( n + 1 );
      a(i,j) = sqrt ( 2.0 / ( n + 1 ) ) * sin ( angle );
    end
  end

  return
end


text:
function [Y] = mnbinpdf(X, W, K)
%MNBINPDF       Maszle's negative binomial probability density function
%
%  Y = MNBINPDF(X, W, K)
%
%  returns the probabilities of the values X using a negative binomial
%  distribution with mean W and clumping parameter K.
%
%  This is an alternate form of the negative binomial commonly
%  employed in biology to describe aggregated count data.
%  W and K are related to the traditional parametrization by the
%  relationships:
%
%      P = K/(K + W),  and K = R,
%
%  with the distinction that K is any positive real number, not just
%  positive integers.
%
%  X should only be integers and is rounded if not.  Calculation is
%  optimized by taking the exponential of the GAMMALN function.
%
%  See also MNBINRND

%----------------------------------------------------------------------
% Copyright (c) 1997.  Don R. Maszle.  All rights reserved.
%
%   -- Revisions -----
%      Author:  Don R. Maszle
%      E-mail:  maze@sparky.berkeley.edu
%   -- SCCS  ---------
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% We use the definition from Bradley and May, Trans. R.Soc.Trop.Med.Hyg.
% v72(3), 1978, p262.

% A notational convenience
a = W/(W + K);
X = round(X);

Y = ((1 - a).^K) ...
    .* (exp( gammaln(X + K) + X.*log(a) - gammaln(K) - gammaln(X + 1)));


text:
function a = upshift ( n )

%*****************************************************************************80
%
%% UPSHIFT returns the UPSHIFT matrix.
%
%  Formula:
%
%    if ( J-I == 1 mod ( n ) )
%      A(I,J) = 1
%    else
%      A(I,J) = 0
%
%  Example:
%
%    N = 4
%
%    0 1 0 0
%    0 0 1 0
%    0 0 0 1
%    1 0 0 0
%
%  Rectangular properties:
%
%    A is integral: int ( A ) = A.
%
%    A is a zero/one matrix.
%
%  Square Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is nonsingular.
%
%    A is a permutation matrix.
%
%    If N is even, det ( A ) = -1.
%    If N is odd,  det ( A ) = +1.
%
%    A is unimodular.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is a Hankel matrix: constant along anti-diagonals.
%
%    A is an N-th root of the identity matrix.
%
%    The inverse of A is the downshift matrix.
%
%    A circulant matrix C, whose first row is (c1, c2, ..., cn), can be
%    written as a polynomial in A:
%
%      C = c1 * I + c2 * A + c3 * A**2 + ... + cn * A**n-1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows and columns
%    of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      if ( i4_modp ( j - i, n ) == 1 )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end
    end
  end

  return
end


text:
function [fe] = tapas_ti_trapezoidal(llh, T)
%% Computes the free energy using the trapezoidal rule
%
% Input
%       llh     -- Samples from the log likelihod
%       T       -- Temperatures at which the samples were drawn
% Output
%       fe      -- Free energy
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

ellh = mean(llh, 1);
fe = trapz(T, llh);

end

author: translationalneuromodeling

repo: tapas


text:
function n = p15_n ( )

%*****************************************************************************80
%
%% P15_N returns the number of variables for problem 15.
%
%  Discussion:
%
%    The number of variables may be any multiple of 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = - 4;

  return
end


text:
function roidb = fast_rcnn_generate_sliding_windows(conf, imdb, roidb, roipool_in_size)
% [pred_boxes, scores] = fast_rcnn_conv_feat_detect(conf, im, conv_feat, boxes, max_rois_num_in_gpu, net_idx)
% --------------------------------------------------------
% Fast R-CNN
% Reimplementation based on Python Fast R-CNN (https://github.com/rbgirshick/fast-rcnn)
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    regions.images = imdb.image_ids;

    im_sizes = imdb.sizes;
    regions.boxes = cellfun(@(x) generate_sliding_windows_one_image(conf, x, roipool_in_size), num2cell(im_sizes, 2), 'UniformOutput', false);

    roidb = roidb_from_proposal(imdb, roidb, regions);
end

function boxes = generate_sliding_windows_one_image(conf, im_size, roipool_in_size)
    im_scale = prep_im_for_blob_size(im_size, conf.scales, conf.max_size);
    im_size = round(im_size * im_scale);

    x1 = 1:conf.feat_stride:im_size(2);
    y1 = 1:conf.feat_stride:im_size(1);
    [x1, y1] = meshgrid(x1, y1);
    x1 = x1(:);
    y1 = y1(:);
    x2 = x1 + roipool_in_size * conf.feat_stride - 1;
    y2 = y1 + roipool_in_size * conf.feat_stride - 1;

    boxes = [x1, y1, x2, y2];
    boxes = filter_boxes(im_size, boxes);

    boxes = bsxfun(@times, boxes-1, 1/im_scale) + 1;
end

function boxes = filter_boxes(im_size, boxes)
    valid_ind = boxes(:, 1) >= 1 & boxes(:, 1) <= im_size(2) & ...
                boxes(:, 2) >= 1 & boxes(:, 2) <= im_size(1) & ...
                boxes(:, 3) >= 1 & boxes(:, 3) <= im_size(2) & ...
                boxes(:, 4) >= 1 & boxes(:, 4) <= im_size(1);

    boxes = boxes(valid_ind, :);
end

author: ShaoqingRen


text:
function [align] = flydra(in,config)

Cst = in.Cst;
Rot = in.Rot;

Octave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if ~Octave,
  drawscene(in.Xe,Cst',Rot,41,'cloud','Graphical Output Validation: View from top or bottom (no sRt)',config.cal.cams2use);
end

% definition of the absolute world frame (in cm)

%% for hummingbird cage
%
%cam(3).C = [168, 173, 185]';
%cam(4).C = [142, 76, 89]';
%cam(1).C = [81, 72, 165]';
%cam(2).C = [22,  15, 176]';
%cam(5).C = [53.5, 175, 78]';

% definition of the absolute world frame (in mm) (Doug's cube)

% cam(1).C = [10, 390, 600]';
% cam(2).C = [460,  380, 600]';
% cam(3).C = [430, 115, 580]';
% cam(4).C = [0, 0, 610]';
% cam(5).C = [240, 150, 610]';

% definition of the absolute world frame (in mm) (windtunnel June 2005)

%cam(1).C = [1510, 180, 500+322]';
%cam(2).C = [1140,  260, 450+322]';
%cam(3).C = [1040, 25, 480+322]';
%cam(4).C = [490, 160, 390+322]';
%cam(5).C = [900, 25, 400+322]';

% the above numbers (from June 2005) transformed to acheive better
% agreement (July 3, 2005)

%cam(1).C = [1585.0266, 157.1808, 725.3687]';
%cam(2).C = [1221.0045, 230.7819, 717.1524]';
%cam(3).C = [1123.2774, -2.1424, 731.8758]';
%cam(4).C = [580.8019, 124.5745, 705.3737]';
%cam(5).C = [978.8653, 3.3878, 663.1311]';

% cal July 5, 2005

%cam(1).C = [1710, 375, 820]';
%cam(2).C = [1110, 316, 800]';
%cam(3).C = [1110, 100, 466]';
%cam(4).C = [-380, 290, 875]';
%cam(5).C = [440, 15, 680]';

% touch-up July 6, 2005 (5c)

%cam(1).C = [1853.5325, 228.1160, 971.5365]';
%cam(2).C = [1202.0386, 199.5437, 953.7752]';
%cam(3).C = [1183.3920, 41.6797, 537.2034]';
%cam(4).C = [-406.8365, 229.6378, 1073.0449]';
%cam(5).C = [460.2840, -60.0974, 778.1218]';

% touch-up Aug 22, 2005 (touch-up from A applied to B)

%cam(1).C = cam(1).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(2).C = cam(2).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(3).C = cam(3).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(4).C = cam(4).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(5).C = cam(5).C - [130.0*0.83, 0., 110.0*0.83]';

% touch-up Jan 27, 2006 ( hacked-in estimate )

%cam(1).C = cam(1).C - [130.0*0.83+20.0, 0., 110.0*0.83]';
%cam(2).C = cam(2).C - [130.0*0.83+10.0, 0., 110.0*0.83]';
%cam(3).C = cam(3).C - [130.0*0.83+10.0, 0., 110.0*0.83]';
%cam(4).C = cam(4).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(5).C = cam(5).C - [130.0*0.83, 0., 110.0*0.83]';

% new Feb 14, 2006

%cam(1).C = [1760, 190, 850]';
%cam(2).C = [1060, 260, 785]';
%cam(3).C = [1070, 70, 470 ]';
%cam(4).C = [-380, 230, 873]';
%cam(5).C = [380, -10, 670]';

% moved cam 2006 03 07 19:45

%cam(1).C = [1760, 190, 850]';
%cam(2).C = [1060, 260, 785]';
%cam(3).C = [1070, 190, 510 ]';
%cam(4).C = [-380, 230, 873]';
%cam(5).C = [380, 120, 670]';

if 1,

% 2006 03 31

%cam(1).C = [640, -105, 170]';
%cam(2).C = [1220, 172, 912]';
%cam(3).C = [475, 172, 580 ]';
%cam(4).C = [-350, 172, 932]';
%cam(5).C = [320, -105, 170]';

% 2006 04 03e

%cam(3).C = [380, 240, 530 ]';

% 2006 04 04a

%cam(4).C = [-115, 220, 800]';
end

% 2006 09 14

%cam(1).C = [685, 220, 775]';
%cam(2).C = [1220, 172, 912]';
%cam(3).C = [380, 240, 530 ]';
%cam(4).C = [-115, 220, 800]';
%cam(5).C = [180, 240, 740]';

%% 2006 10 23 (from DLT)
%cam(1).C = [ 651.7963301   169.95615265  749.51356735]';
%cam(2).C = [ 1090.40810501   190.82014771   866.8980674 ]';
%cam(3).C = [ 327.41874025  237.43022637  539.80942371]';
%cam(4).C = [-241.364468    194.27934314  767.35249515]';
%cam(5).C = [ 147.20028399  171.45774441  725.08522715]';

% 2006 12 19 riverside hummingbirds
% estimates from DLT
%cam(1).C = [976.22314128 -2289.85157108   489.78104927]';
%cam(2).C = [-17278.94752856   1794.85211148  15855.51116479]';
%cam(3).C = [ 2266.45466953  2590.2964507    -74.53159994]';
%cam(4).C = [-329.01069449  796.46308431  624.04191004]';

% hand measurements
%ft2mm = 12*25.4;
%cam(1).C = [ 1350 -1580 540]';
%cam(2).C = [ -450 770 (750+5*ft2mm) ]';
%cam(3).C = [ 1500 -1600 2250]';
%cam(4).C = [-280 820 570]';

% 2006 12 01 DLT in flydra WT
cam(1).C = [867.37446261  198.97798549  832.51824736]';
cam(2).C = [1194.97461953 206.718389208 898.66723201]';
cam(3).C = [459.669418616 168.068949056 713.890827]';
cam(4).C = [-186.0833207 207.892635698 775.298571621]';
cam(5).C = [115.230367611 202.584140956 904.723766629]';

% 2007 10 03 DLT in flydra WT
cam(1).C = [ 985.827871     182.49267738  1134.97921258]';
cam(2).C = [362.27817754   189.63204145  1134.23310277]';
cam(3).C = [721.10852967  167.31210242  932.32564817]';
cam(4).C = [-330.47764957   168.38518854  1014.81801646]';
cam(5).C = [1774.53765317   144.79725152   928.10954671]';

% of the similarity computation

[align.simT.s, align.simT.R, align.simT.t]  = estsimt([Cst'],[cam(:).C]);
[align.P, align.X] = align3d(in.Pe,in.Xe,align.simT);
% save aligned data
if 1 % SAVE_STEPHI | SAVE_PGUHA
	[align.Cst,align.Rot] = savecalpar(align.P,config);
end

if ~Octave,
  drawscene(align.X,align.Cst',align.Rot,61,'cloud','Graphical Output Validation: Aligned data',config.cal.cams2use);

  set(gca,'CameraTarget',[0,0,0]);
  set(gca,'CameraPosition',[0,0,1]);

  figure(61),
  % print -depsc graphevalaligned.eps
  eval(['print -depsc ', config.paths.data, 'topview.eps'])

  drawscene(align.X,align.Cst',align.Rot,62,'cloud','Graphical Output Validation: Aligned data',config.cal.cams2use);

  set(gca,'CameraTarget',[0,0,0.9]);
  set(gca,'CameraPosition',[2,0,0.9]);

  %figure(62),
  % print -depsc graphevalaligned.eps
  %eval(['print -depsc ', config.paths.data, 'sideview.eps'])
end

return

author: strawlab

repo: MultiCamSelfCal


text:
function gcm = fmri_scm2gcm(X,Nnnc,TR,tPreStim,delta,tau)
%
% gcm = fmri_scm2gcm(X,Nnnc,TR,tPreStim,delta,tau)
%
% Produces a Gamma Convolution Matrix from a stimulus
% convolution matrix (X) and parameters of the gamma
% function (delta, tau).  The gamma functions
% are interpreted as basis vectors.
%
% $Id: fmri_scm2gcm.m,v 1.3 2005/06/01 01:04:52 sayres Exp $

[Ntp Nch Nr] = size(X);
Nh = Nch/Nnnc;
Ng = length(delta);

t = TR*[0:Nh-1] - tPreStim;
h = fmri_hemodyn(t,delta,tau);
h = h./(repmat(max(h),[Nh 1]));

h_all = zeros(Nch,Nnnc*Ng);
h0 = zeros(Nh,Nnnc*Ng);
h0(1:Nh,1:Ng) = h;
for c = 1:Nnnc,
    r1 = Nh*(c-1)+1;
    r2 = r1 + Nh - 1;
    h_all(r1:r2,:) = fmri_shiftcol(h0,Ng*(c-1));
end

gcm = zeros(Ntp,Nnnc*Ng,Nr);

for r = 1:Nr,
    gcm(:,:,r) = X(:,:,r)*h_all;
end

return;

author: vistalab

repo: vistasoft