text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returs the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

for i = 1:K
    centroids(i,:) = mean(X(idx==i,:),1);
end

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

end

repo: MachineLearningEx


text:
function [error_avg,error_A] = calc_abundance_error(A1,A2)
%CALC_ABUNDANCE_ERROR Calculate RMSE for each column of A1 and A2
if ~isfloat(A1)
    A1 = double(A1);
end

if length(unique(A1(:))) == 2
    is_real_dataset = 1;
else
    is_real_dataset = 0;
end

P = permute_abundances(A1,A2);
if nanmean(nanmean(abs(eye(size(A1,2)) - P))) ~= 0
    A2_1 = (P*A2')';
else
    A2_1 = A2;
end

if is_real_dataset
    inds = any(A1==1,2);
%     error_A = nanmean(abs(A1(inds,:) - A2_1(inds,:)));
    error_A = calc_abundance_rmse(A1(inds,:), A2_1(inds,:));
else
%     error_A = nanmean(abs(A1-A2_1));
    error_A = calc_abundance_rmse(A1, A2_1);
end

error_avg = nanmean(error_A);

function rmse = calc_abundance_rmse(A1,A2)
rmse = sqrt(nanmean(abs(A1 - A2).^2, 1));

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [refinedModel,summary] = performDataDrivenRefinement(model, microbeID, biomassReaction, database, inputDataFolder,summary)
% This function is part of the DEMETER pipeline and performs data-driven
% refinement of a genome-scale reconstruction based on available
% species-specific experimental data.
%
% USAGE
%      [model,summary] = performDataDrivenRefinement(model, microbeID, biomassReaction, database, inputDataFolder,summary)
%
% INPUTS
% model             COBRA model structure to refine
% microbeID         ID of the reconstructed microbe that serves as the
%                   reconstruction name and to identify it in input tables
% inputDataFolder   Folder with input tables with experimental data and
%                   databases that inform the refinement process
% summary           Structure with information on performed refinement
%
% OUTPUT
% refinedModel      Refined COBRA model structure
% summary           Structure with information on performed refinement
%
% .. Authors:
%       - Almut Heinken and Stefania Magnusdottir, 2016-2021

%% Fermentation pathways
% Based on the fermentation pathway data for the microbe (table prepared above),
% add and remove reactions as defined in the following script.
% Perform fermentation pathway gap fill:
[model, addedRxns, removedRxns] = fermentationPathwayGapfill(model, microbeID, database, inputDataFolder);
summary.('addedRxns_fermentation') = addedRxns;
summary.('removedRxns_fermentation') = removedRxns;
%% Carbon sources
% Based on the carbon source data for the microbe (table prepared above), add
% and remove reactions as defined in the following script.
% Perform carbon source pathway gap fill:
[model, addedRxns, removedRxns] = carbonSourceGapfill(model, microbeID, database, inputDataFolder);
summary.('addedRxns_carbonSources') = addedRxns;
summary.('removedRxns_carbonSources') = removedRxns;
%% Putrefaction pathways
[model,putrefactionRxnsAdded]=putrefactionPathwaysGapfilling(model,microbeID,database);
summary.('putrefactionRxnsAdded') = putrefactionRxnsAdded;
%% Known secretion products
[model,secretionRxnsAdded] = secretionProductGapfill(model,microbeID,database,inputDataFolder);
summary.('secretionRxnsAdded') = secretionRxnsAdded;
%% Known consumed metabolites
[model,uptakeRxnsAdded] = uptakeMetaboliteGapfill(model,microbeID,database, inputDataFolder);
summary.('uptakeRxnsAdded') = uptakeRxnsAdded;

%% test pathways to make sure they work

model = rebuildModel(model,database,biomassReaction);

FNs = {};
% Carbon sources
[TruePositives, FalseNegatives] = testCarbonSources(model, microbeID, biomassReaction, database, inputDataFolder);
FNs=union(FNs,FalseNegatives);

% Metabolite uptake
[TruePositives, FalseNegatives] = testMetaboliteUptake(model, microbeID, biomassReaction, database, inputDataFolder);
FNs=union(FNs,FalseNegatives);

%% gapfill if there are any false negatives
osenseStr='min';

dataDrivenGapfill={};
if ~isempty(FNs)
    for j=1:length(FNs)
        metExch=['EX_' database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j})),1} '(e)'];
        % find reactions that could be gap-filled to enable flux
        % seems to try to fix exchanges that are not part of the model,
        % leading it to crash:
        %         Objective reactions not found in model!
        %
        % Error in runGapfillingFunctions (line 41)
        % model = changeObjective(model, objectiveFunction);
     %   add exchange reaction
     if isempty(find(strcmp(model.rxns,metExch)))
         met = [database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j}))} '(e)'];
         model = addExchangeRxn(model,met,-1,1000);
     end
     [model,gapfilledRxns] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
     dataDrivenGapfill=union(dataDrivenGapfill,gapfilledRxns);
     %  end
    end
    if ~isempty(dataDrivenGapfill)
        summary.('DataDrivenGapfill')=dataDrivenGapfill;
    end
end

% Fermentation products
[TruePositives, FalseNegatives] = testFermentationProducts(model, microbeID, biomassReaction, database, inputDataFolder);
FNs=union(FNs,FalseNegatives);

% Putrefaction pathways
[TruePositives, FalseNegatives] = testPutrefactionPathways(model, microbeID, biomassReaction, database);
FNs=union(FNs,FalseNegatives);

% Secretion products
[TruePositives, FalseNegatives] = testSecretionProducts(model, microbeID, biomassReaction, database, inputDataFolder);
FNs=union(FNs,FalseNegatives);

% gapfill if there are any false negatives
osenseStr='max';

if ~isempty(FNs)
    for j=1:length(FNs)
        metInd=find(strcmp(database.metabolites(:,2),FNs{j}));
        if length(metInd)>1
            metInd=metInd(1);
        end
        metExch=['EX_' database.metabolites{metInd,1} '(e)'];
        if isempty(find(ismember(model.rxns,metExch)))
            % reaction ID itself provided
            metExch = FNs{j};
        end
        % find reactions that could be gap-filled to enable flux
        try
        [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
        catch
        end
        summary.('conditionSpecificGapfill') = union(summary.('conditionSpecificGapfill'),condGF);
        summary.('targetedGapfill') = union(summary.('targetedGapfill'),targetGF);
        summary.('relaxFBAGapfill') = union(summary.('relaxFBAGapfill'),relaxGF);
    end
end

%% change back to biomass reaction
refinedModel=changeObjective(model,biomassReaction);

end

author: opencobra

repo: cobratoolbox


text:
function F=scarsplslda(X,y,A,K,method,num,OPT)
%+++ This is the simplified version of CARS, only retaining the EDF element.
%+++ num: the number of Monte Carlo Sampling.
%+++ A:   the maximal number of PLS components to extract.
%+++ fold: number of folds for cross validation.
%+++ method: pretreat method, 'center' or 'autoscaling'.
%+++ OPT: 1: Plot. 0: No plot.
%+++ Advisor: Yizeng Liang, yizeng_liang@263.net.
%+++ Hongdong Li, Jan.3, 2009.
%+++ Reference: Hongdong Li, Yizeng Liang, Qingsong Xu and Dongsheng Cao,
%         Key wavelengths screening using competitive adaptive reweighted sampling method
%         for multivariate calibration J?. Anal. Chim. Acta, 2009, 648 (1): 77-84

%+++ Initial settings.
if nargin<7;OPT=1;end;  %+++ OPT==1: then figure output.
if nargin<6;num=50;end
if nargin<5;method='autoscaling';end
if nargin<4;K=5;end
if nargin<3;A=2;end

[Mx,Nx]=size(X);
A=min([Mx Nx A]);
index=1:Nx;

r0=1;
r1=2/Nx;
Vsel=1:Nx;

W=zeros(num,Nx);
Ratio=zeros(1,num);

%+++ Parameter of exponentially decreasing function.
b=log(r0/r1)/(num-1);  a=r0*exp(b);

%+++ Main Loop
for iter=1:num

     LDA=plslda(X(:,Vsel),y,A,method);    %+++ PLS model
     w=zeros(Nx,1);coef=LDA.coef_lda_origin(1:end-1);
     w(Vsel)=coef;W(iter,:)=w;
     w=abs(w);                                  %+++ weights
     [ws,indexw]=sort(-w);                      %+++ sort weights

     ratio=a*exp(-b*(iter+1));                      %+++ Ratio of retained variables.
     Ratio(iter)=ratio;
     K=round(Nx*ratio);

     w(indexw(K+1:end))=0;                      %+++ Eliminate some variables with small coefficients.

     Vsel=find(w~=0);
     fprintf('The %dth variable sampling finished.\n',iter);    %+++ Screen output.
 end

%+++  Cross-Validation to choose an optimal subset;
RMSEP=zeros(1,num);
Rpc=zeros(1,num);
for i=1:num
   vsel=find(W(i,:)~=0);
   CV=plsldacv(X(:,vsel),y,A,10,method,0);
   RMSEP(i)=min(CV.cv);
   Rpc(i)=CV.optPC;
   fprintf('The %dth subset finished.\n',i);
end
Rmin=min(RMSEP);
indexOPT=find(RMSEP==Rmin);
indexOPT=indexOPT(end);

%+++ output
F.method=method;
F.nPC=A;
F.W=W;
F.vim=sum(W,2);
F.cv=RMSEP;
F.minCV=Rmin;
F.iterOPT=indexOPT;
F.optPC=Rpc(indexOPT);
F.ratio=Ratio;
F.vsel=find(W(indexOPT,:)~=0)';
F.vim=abs(sum(W));

%+++ Plot
 if OPT==1;plotcars(F);end;
%+++ END

author: viggin

repo: domain-adaptation-toolbox


text:
%% This function applies Canny Operator to an input
%% image and can be used for edge detection.

%%% input arguments
%% image = String representing input file.
%% si = Size of Gaussian Kernel used for antialiasing the input image and calculating
%%      gaussian derivatives.
%% sigma = Standard deviation of Gaussian function.
%% thresh = upper threshold used for thresholding in Canny method.

%%% output Arguments
%% gnh = output Image.

%%% References :- DIGITAL IMAGE PROCESSING Third Edition,
%%                Rafael C. Gonzalez
%%                Richard E. Woods

function gnh= canny(image,si,sigma,thresh)
%% input argument check.

error(nargchk(4,4,nargin));
if(thresh>1)
 error('Thresh value should be between 0 and 1');
end
I=imread(image);
if(size(I,3)==3)
I=rgb2gray(I);
end
figure,imagesc(I),impixelinfo,title('Original Image'),colormap('gray');
lthresh=0.4*thresh;
uthresh=thresh;
if(uthresh>1.0)
    uthresh=1;
end
%% Antialias the image by convolving with gaussian filter

h=fspecial('gaussian',si,sigma);
I=im2double(I);
I=imfilter(I,h,'conv');
figure,imagesc(I),impixelinfo,title('Original Image after Convolving with gaussian'),colormap('gray');
%% Compute Gaussian derivatives

 x=-si:1:si;
 y=-si:1:si;
gaussx=-(x/(sigma*sigma)).*exp(-(x.*x+y.*y)/(2*sigma*sigma));
gaussy=gaussx';
Ix=imfilter(I,gaussx,'conv');
Iy=imfilter(I,gaussy,'conv');
%% Compute magnitude and orientation of gradient vector.

Mag=sqrt(Ix .^ 2 + Iy .^ 2);
Magmax=max(Mag(:));
Mag=Mag/Magmax;
figure,imagesc(Mag),impixelinfo,title('Magnitude of gradient Vectors'),colormap('gray');
gnl=zeros(size(I,1),size(I,2));
gnh=zeros(size(I,1),size(I,2));
angle=atand(Iy./Ix);
figure,quiver(20*Ix,20*Iy,0),title('Orientation of gradient vectors'),colormap('gray');
%% Apply non maxima suppression and thresholding. The orientation of gradient vecctor is
%% resolved into four directions.

for y=2:1:size(I,1)-1
    for x=2:1:size(I,2)-1
        if(((angle(y,x)>=-22.5)&&(angle(y,x)<22.5))||((angle(y,x)<-157.5))||(angle(y,x)>157.5))
            if((Mag(y,x)>Mag(y,x+1)) && (Mag(y,x)>Mag(y,x-1)))
                if(Mag(y,x)>=lthresh)
                    if(Mag(y,x)>=uthresh)
                gnh(y,x)=1;
                    else
                    gnl(y,x)=1;
                    end
                end
            end
        end
        if(((angle(y,x)>=-112.5)&&(angle(y,x)<-67.5))||((angle(y,x)>=67.5)&&(angle(y,x)<112.5)))
            if((Mag(y,x)>Mag(y+1,x)) && (Mag(y,x)>Mag(y-1,x)))
                if(Mag(y,x)>=lthresh)
                    if(Mag(y,x)>=uthresh)
                gnh(y,x)=1;
                    else
                    gnl(y,x)=1;
                    end
                end
            end
        end
        if(((angle(y,x)>=-67.5)&&(angle(y,x)<-22.5))||((angle(y,x)>=112.5)&&(angle(y,x)<157.5)))
            if((Mag(y,x)>Mag(y-1,x+1)) && (Mag(y,x)>Mag(y+1,x-1)))
                if(Mag(y,x)>=lthresh)
                    if(Mag(y,x)>=uthresh)
                gnh(y,x)=1;
                    else
                    gnl(y,x)=1;
                    end
                end
            end
        end
        if(((angle(y,x)>=-157.5)&&(angle(y,x)<-112.5))||((angle(y,x)>=22.5)&&(angle(y,x)<67.5)))
            if((Mag(y,x)>Mag(y+1,x+1)) && (Mag(y,x)>Mag(y-1,x-1)))
                if(Mag(y,x)>=lthresh)
                    if(Mag(y,x)>=uthresh)
                gnh(y,x)=1;
                    else
                    gnl(y,x)=1;
                    end
                end
            end
        end
    end
end
figure,imagesc(gnh),impixelinfo,title('Image after thresholding'),colormap('gray');
%% Apply Connectivity Analysis.

[row,col]=find(gnh>0);
    for t=1:1:size(col)
        if(gnl(row(t)+1,col(t))>0)
            gnh(row(t)+1,col(t))=1;
            gnl(row(t)+1,col(t))=0;
        end
        if(gnl(row(t),col(t)+1)>0)
            gnh(row(t),col(t)+1)=1;
            gnl(row(t),col(t)+1)=0;
        end
        if(gnl(row(t)+1,col(t)+1)>0)
            gnh(row(t)+1,col(t)+1)=1;
            gnl(row(t)+1,col(t)+1)=0;
        end
        if(gnl(row(t)-1,col(t)-1)>0)
            gnh(row(t)-1,col(t)-1)=1;
            gnl(row(t)-1,col(t)-1)=0;
        end
        if(gnl(row(t)+1,col(t)-1)>0)
            gnh(row(t)+1,col(t)-1)=1;
            gnl(row(t)+1,col(t)-1)=0;
        end
        if(gnl(row(t)-1,col(t)+1)>0)
            gnh(row(t)-1,col(t)+1)=1;
            gnl(row(t)-1,col(t)+1)=0;
        end
        if(gnl(row(t)-1,col(t))>0)
            gnh(row(t)-1,col(t))=1;
            gnl(row(t)-1,col(t))=0;
        end
        if(gnl(row(t),col(t)-1)>0)
            gnh(row(t),col(t)-1)=1;
            gnl(row(t),col(t)-1)=0;
        end
    end
    %% Display output Image.

figure,imagesc(gnh),impixelinfo,title('Image after Applying Connectivity Analysis'),colormap('gray');
end

author: Sable

repo: mcbench-benchmarks


text:
function [ o, x, w ] = epn_glg_02_xiu ( n, alpha )

%*****************************************************************************80
%
%% EPN_GLG_02_XIU implements the Xiu rule for region EPN_GLG.
%
%  Discussion:
%
%    The rule has order
%
%      O = N + 1.
%
%    The rule has precision P = 2.
%
%    EPN_GLG is the N-dimensional positive space [0,+oo)^N with generalized
%    Laguerre weight function:
%
%      w(alpha;x) = product ( 1 <= i <= n ) x(i)^alpha exp ( - x(i) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongbin Xiu,
%    Numerical integration formulas of degree two,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 1515-1520.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, the exponent of X in the weight function.
%    -1.0 < ALPHA.
%
%    Input, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EPN_GLG_02_XIU - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= -1.0\n' );
    error ( 'EPN_GLG_02_XIU - Fatal error!' );
  end

  o = n + 1;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  for j = 1 : o

    i = 0;
    for r = 1 : floor ( n / 2 )
      arg = 2 * r * ( j - 1 ) * pi / ( n + 1 );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * cos ( arg );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * sin ( arg );
    end

    if ( i < n )
      i = i + 1;
      x(i,j) = r8_mop ( j - 1 );
    end

  end

  gamma0 = - 1.0;
  delta0 = alpha + 1.0;
  c1 = - alpha - 1.0;

  x(1:n,1:o) = ( sqrt ( gamma0 * c1 ) * x(1:n,1:o) - delta0 ) / gamma0;

  expon = 0;
  volume_1d = ep1_glg_monomial_integral ( expon, alpha );
  volume = volume_1d ^ n;

  w(1:o) = volume / o;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function makecmap(bins,fighandle)
    %  MAKECMAP Interactive colormap maker.
    %
    %     usage:    makecmap
    %               makecmap(map)
    %               makecmap(map,fighandle)
    %               makecmap(bins)
    %    		makecmap(bins,fighandle)
    %		makecmap(fighandle,'fig')
    %
    %       bins - size of colormap (default - 64)
    %       map - rgb colormap (matrix)
    %       fighandle - figure(s); to apply colormap (optional)
    %	fighandle,'fig' - reads colormap from figure (fighandle)
    %
    %
    %     Nodes are displayed as circles. To move a node use mouse
    %       button #1. To add a node use mouse button #2. To delete use
    %       button #3.
    %
    %    Mouse buttons     button 1 - move    (3 button mouse)
    %                      button 2 - add
    %                      button 3 - erase
    %
    %    There are four graphs. The first three corresond to the relative
    %       rgb values. And the fourth graph is a luminance factor
    %       which the rgb values are multiplied by.
    %
    %     colormap = [red,green,blue] * luminance
    %
    %    Reset button - resets the colormap to the original values.
    %    Export button - assigns the colormap to a desktop variable.
    %    Close button - closes the figure window.
    %
    %    Runs under Matlab 5
    %

    % Note: the readcolormap option is not perfect. It works well if the
    % colormap contains straight lines like most matlab colormaps. You can play
    % with the tolerance. Also, the more nodes in an axes the slower the mouse
    % handler routines will run.

    % Written by Colin Humphries, Salk Institute, March 1997
    %             colin@salk.edu
    %

    % TODO untangle the confusing global situation -CGR

    %report_this_filefun();

    nargs = nargin;

    if nargs == 0
        bins = 64;
        nargs = 1;
    end

    if ~ischar(bins)

        if size(bins,2) == 3
            map = bins;
            bins = size(map,1);
            if nargs == 1
                fighandle = [];
            end
        else
            map = [];
            if nargs == 2
                if ischar(fighandle)
                    fighandle = bins;
                    map = get(fighandle,'Colormap');
                    bins = size(map,1);
                end
            else
                fighandle = [];
            end
        end

        if isempty(map)
            red = [1 0;bins 1];       % default colormap
            green = [1 0;bins 1];
            blue = [1 0;bins 1];
        else
            % Note: this routine cycles through all
            % the data points and assigns nodes to
            % the points where straight lines end.
            % This does not work well in every
            % occasion. If the colormap is nonlinear
            % in many places, a lot of nodes will be
            % assigned.

            for j = 1:3           % cycle through each column
                for i = 1:bins      % cycle through each row
                    if i == 1         % assign the first point as the first node
                        nodes(1,:) = [1 map(1,j)];
                        lastnode = nodes;   % last node assigned
                    elseif i == 2
                        oldtempnode = [i map(2,j)]; % assigns point 2 as the test node
                    else
                        tempnode = [i map(i,j)];    % look at current point
                        if (abs((oldtempnode(2)-lastnode(2))/(oldtempnode(1)...
                                -lastnode(1))) <= abs((tempnode(2)-lastnode(2))...
                                /(tempnode(1)-lastnode(1)))*1.01) & ...
                                (abs((oldtempnode(2)-lastnode(2))/...
                                (oldtempnode(1)-lastnode(1))) >= ...
                                abs((tempnode(2)-lastnode(2))/...
                                (tempnode(1)-lastnode(1)))*.99) % Check whether the line
                            % from the last node to
                            % tempnode and
                            % oldtempnode have the
                            % same slope.
                            % Tolerance is +/- 1%

                            oldtempnode = tempnode;    % if so then move oldtempnode
                        else
                            nodes = [nodes;oldtempnode]; % if not then assign a new
                            lastnode = oldtempnode;      % node.
                            oldtempnode = tempnode;
                        end
                    end
                end
                nodes = [nodes;bins map(bins,j)];  % the last node is the last point
                if j == 1
                    red = nodes;  % get red nodes
                    nodes = [];
                elseif j == 2
                    green = nodes;% get green nodes
                    nodes = [];
                elseif j == 3
                    blue = nodes; % get blue nodes
                    nodes = [];
                end
            end
        end

        mcmfighandle = figure_w_normalized_uicontrolunits('color','k');

        if isempty(fighandle)
            fighandle = mcmfighandle;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Set up Colormap from nodes
        %   Converts each set of nodes into a discrete vector

        for i = 1:size(red,1)-1
            W = linspace(red(i,2),red(i+1,2),red(i+1,1)-red(i,1)+1);
            mapred(red(i,1):red(i+1,1)-1) = W(1:red(i+1,1)-red(i,1))';
        end
        mapred(red(i+1,1)) = red(i+1,2);
        for i = 1:size(green,1)-1
            W = linspace(green(i,2),green(i+1,2),green(i+1,1)-green(i,1)+1);
            mapgreen(green(i,1):green(i+1,1)-1) = W(1:green(i+1,1)-green(i,1))';
        end
        mapgreen(green(i+1,1)) = green(i+1,2);
        for i = 1:size(blue,1)-1
            W = linspace(blue(i,2),blue(i+1,2),blue(i+1,1)-blue(i,1)+1);
            mapblue(blue(i,1):blue(i+1,1)-1) = W(1:blue(i+1,1)-blue(i,1))';
        end
        mapblue(blue(i+1,1)) = blue(i+1,2);
        map = [mapred(:),mapgreen(:),mapblue(:)]; % get new colormap

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Set up Graphs

        nodenum = size(red,1);
        subplot(3,2,1)                 % Red Graph
        set(gca,'NextPlot','add')
        for i = 1:nodenum-1
            plot(red(i,1),red(i,2),'ow')
            line([red(i,1),red(i+1)],[red(i,2),red(i+1,2)],'color','r')
        end
        plot(red(nodenum,1),red(nodenum,2),'ow')
        axis([1 red(nodenum,1) 0 1])
        redax = gca;
        set(redax,'UserData',red,'tag','r','box','on','color','k',...
            'xcolor','w','ycolor','w','xgrid','on','ygrid','on',...
            'Ytick',(0:.2:1));
        ylabel('Red','FontSize',14)

        subplot(3,2,3)                % Green Graph
        nodenum = size(green,1);
        set(gca,'NextPlot','add')
        for i = 1:nodenum-1
            plot(green(i,1),green(i,2),'ow')
            line([green(i,1),green(i+1)],[green(i,2),green(i+1,2)],'color','g')
        end
        plot(green(nodenum,1),green(nodenum,2),'ow')
        axis([1 green(nodenum,1) 0 1])
        greenax = gca;
        set(greenax,'UserData',green,'tag','g','box','on','color','k',...
            'xcolor','w','ycolor','w','xgrid','on','ygrid','on',...
            'YTick',(0:.2:1));
        ylabel('Green','FontSize',14)

        subplot(3,2,5)                % Blue Graph
        nodenum = size(blue,1);
        set(gca,'NextPlot','add')
        for i = 1:nodenum-1
            plot(blue(i,1),blue(i,2),'ow')
            line([blue(i,1),blue(i+1)],[blue(i,2),blue(i+1,2)],'color','b')
        end
        plot(blue(nodenum,1),blue(nodenum,2),'ow')
        axis([1 blue(nodenum,1) 0 1])
        blueax = gca;
        set(blueax,'UserData',blue,'tag','b','box','on','color','k',...
            'xcolor','w','ycolor','w','xgrid','on','ygrid','on',...
            'YTick',(0:.2:1));
        ylabel('Blue','FontSIze',14)

        colormap(map)                % Apply Colormap
        set(fighandle,'Colormap',map)

        transax = axes('Position',[.578 .4056 .327 .2238]); % Transfer Graph
        set(gca,'NextPlot','add')
        trans = [1 1;bins 1];
        nodenum = 2;
        cla
        for i = 1:nodenum-1
            plot(trans(i,1),trans(i,2),'ow')
            line([trans(i,1),trans(i+1)],[trans(i,2),trans(i+1,2)],'color','y')
        end
        plot(trans(nodenum,1),trans(nodenum,2),'ow')
        axis([1 trans(nodenum,1) 0 1])
        trans = ones(1,bins);
        title('luminance','color','w','FontSize',14)
        set(transax,'UserData',[1 1;bins 1],'tag','y','box','on','color','k',...
            'xcolor','w','ycolor','w','xgrid','on','ygrid','on',...
            'YTick',(0:.2:1))
        set(mcmfighandle,'UserData',[map,trans(:)],'tag',int2str(fighandle))

        cbax = axes('Position',[.578 .75 .327 .1072]);  % Colorbar Graph
        colorbar(cbax)
        set(cbax,'tag','cb','xcolor','w','ycolor','w',...
            'UserData',[size(red,1),0;...
            size(green,1),0;...
            size(blue,1),0;...
            red;...
            green;...
            blue])

        % Title axis

        hdaxes = axes('Position',[.57 .93 .327 .05],'Visible','off',...
            'tag','hd','UserData',map);
        text(.5,0,'Colormap Editor','FontSize',16,'Color','w',...
            'HorizontalAlignment','center')

        % Uicontrols

        uicontrol('style','pushbutton','string','Reset',...
            'Units','Normalized',...
            'Position',[.55 .0929 .1071 .0595],...
            'callback',@callbackfun_001)

        uicontrol('style','pushbutton','string','Export',...
            'Units','Normalized',...
            'Position',[.675 .0929 .1071 .0595],...
            'callback',@callbackfun_002)

        uicontrol('style','pushbutton','string','Close',...
            'Units','Normalized',...
            'Position',[.80 .0929 .1071 .0595],...
            'callback',@callbackfun_003)

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Set up callbacks

        set(mcmfighandle,'WindowButtonDownFcn','makecmap(''down'');')
        set(mcmfighandle,'WindowButtonUpFcn','makecmap(''up'');')
        set(mcmfighandle,'WindowButtonMotionFcn','makecmap(''motion'');')
        set(mcmfighandle,'CloseRequestFcn','makecmap(''close'');')

        % End of main function
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    else
        if strcmp(bins,'down')           % Mouse Button Down Handler
            figh = gcbf;
            axhandle = gca;
            lcolor = get(axhandle,'tag');  % find out what axis we are in
            if strcmp(lcolor,'cb') | strcmp(lcolor,'hd')
                return           % if colorbar or title axis then return
            end
            xlimits = get(axhandle,'Xlim');
            ylimits = get(axhandle,'Ylim');
            nodes = get(axhandle,'UserData');
            nodenum = size(nodes,1);      % get nodes
            P = get(axhandle,'CurrentPoint'); % get mouse position
            X = P(1,1);
            Y = P(1,2);
            button = get(gcbf,'SelectionType'); % get button
            if X <= xlimits(1)-xlimits(2)*.05 || X >= xlimits(2)*1.05 || ...
                    Y <= ylimits(1)-ylimits(2)*.07 || Y >= ylimits(2)*1.07
                return    % if outside of axis limits then return
            end

            if strcmp(button,'normal')    % left mouse button
                global M_BUTTON_DOWN M_PNTS M_LINES M_TEXT
                ii = find(X > nodes(:,1));  % find node above point
                minnode = max(ii);
                ii = find(X < nodes(:,1));  % find node below point
                maxnode = min(ii);
                if (X - nodes(minnode,1)) > (nodes(maxnode,1)-X) % if closer to maxnode
                    if maxnode == nodenum;     % if maxnode is the endnode
                        nodes(maxnode,2) = max(min(Y,1),0);  % only move Y-position
                    else
                        nodes(maxnode,1) = round(X);         % move node
                        nodes(maxnode,2) = max(min(Y,1),0);
                    end
                    M_BUTTON_DOWN = maxnode;               % save what node was selected
                elseif (X - nodes(minnode,1)) < (nodes(maxnode,1)-X)% if closer to minnode
                    if minnode == 1     %if firstnode
                        nodes(minnode,2) = max(min(Y,1),0); % only move Y-position
                    else
                        nodes(minnode,1) = round(X);        % move node
                        nodes(minnode,2) = max(min(Y,1),0);
                    end
                    M_BUTTON_DOWN = minnode;              % save what node was selected
                end
                cla
                for i = 1:nodenum-1                 % Redraw axis
                    M_PNTS(i) = plot(nodes(i,1),nodes(i,2),'ow');
                    M_LINES(i) = line([nodes(i,1),nodes(i+1)],...
                        [nodes(i,2),nodes(i+1,2)],'color',lcolor);
                end
                M_PNTS(i+1) = plot(nodes(nodenum,1),nodes(nodenum,2),'ow');

                M_TEXT = text(0,1.05,...
                    [' ',int2str(nodes(M_BUTTON_DOWN,1)),', ',...
                    num2str(nodes(M_BUTTON_DOWN,2),2)],...
                    'Color','w','clipping','off',...
                    'VerticalAlignment','bottom',...
                    'FontSize',10);

                axis([1 nodes(nodenum,1) 0 1])

            elseif strcmp(button,'extend')     % Middle Mouse Button
                %global M_BUTTON_DOWN M_PNTS M_LINES M_TEXT

                ii = find(X > nodes(:,1));       % find node below
                minnode = max(ii);
                ii = find(X < nodes(:,1));       % find node above
                maxnode = min(ii);
                nodes = [nodes(1:minnode,:);[round(X),max(min(Y,1),0)];...
                    nodes(maxnode:nodenum,:)];   % add node shifting those above
                nodenum = size(nodes,1);
                M_BUTTON_DOWN = minnode+1;
                cla
                for i = 1:nodenum-1              % Redraw axis
                    M_PNTS(i) = plot(nodes(i,1),nodes(i,2),'ow');
                    M_LINES(i) = line([nodes(i,1),nodes(i+1)],...
                        [nodes(i,2),nodes(i+1,2)],'color',lcolor);
                end
                M_PNTS(i+1) = plot(nodes(nodenum,1),nodes(nodenum,2),'ow');
                M_TEXT = text(0,1.05,...
                    [' ',int2str(nodes(M_BUTTON_DOWN,1)),', ',...
                    num2str(nodes(M_BUTTON_DOWN,2),2)],...
                    'Color','w','clipping','off',...
                    'VerticalAlignment','bottom',...
                    'FontSize',10);

                axis([1 nodes(nodenum,1) 0 1])

            elseif strcmp(button,'alt')          % Right Mouse Button
                ii = find(X > nodes(:,1));         % find node below
                minnode = max(ii);
                ii = find(X < nodes(:,1));         % find node above
                maxnode = min(ii);
                if (X - nodes(minnode,1)) > (nodes(maxnode,1)-X)    % reassign nodes
                    if maxnode < nodenum
                        nodes = [nodes(1:maxnode-1,:);nodes(maxnode+1:nodenum,:)];
                    end
                elseif (X - nodes(minnode,1)) < (nodes(maxnode,1)-X)
                    if minnode > 1
                        nodes = [nodes(1:minnode-1,:);nodes(minnode+1:nodenum,:)];
                    end
                end
                nodenum = size(nodes,1);
                cla
                for i = 1:nodenum-1       %redraw axis
                    plot(nodes(i,1),nodes(i,2),'ow')
                    line([nodes(i,1),nodes(i+1)],[nodes(i,2),nodes(i+1,2)],'color',lcolor)
                end
                plot(nodes(nodenum,1),nodes(nodenum,2),'ow')
                axis([1 nodes(nodenum,1) 0 1])
            end
            set(axhandle,'UserData',nodes)  % update nodes

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        elseif strcmp(bins,'up')            % Mouse Button Up Handler
            figh = gcbf;
            %global M_TEXT
            delete(M_TEXT)
            clear global M_BUTTON_DOWN M_LINES M_PNTS M_TEXT % clear variables
            UserData = get(figh,'UserData');
            fighandle = str2num(get(figh,'tag')); % get figure to apply colormap
            axhandle = gca;
            nodes = get(axhandle,'UserData');    % get nodes
            tagval = get(axhandle,'tag');
            if strcmp(tagval,'cb') | strcmp(tagval,'hd')
                return   % if colorbar or title axis then return
            end
            for i = 1:size(nodes,1)-1
                W = linspace(nodes(i,2),nodes(i+1,2),nodes(i+1,1)-nodes(i,1)+1);
                map(nodes(i,1):nodes(i+1,1)-1) = W(1:nodes(i+1,1)-nodes(i,1))';
            end
            map(nodes(i+1,1)) = nodes(i+1,2);
            if strcmp(tagval,'r')
                UserData(:,1) = map';
            elseif strcmp(tagval,'g')
                UserData(:,2) = map';
            elseif strcmp(tagval,'b')
                UserData(:,3) = map';
            elseif strcmp(tagval,'y')
                UserData(:,4) = map';
            end
            set(figh,'UserData',UserData,'Colormap',...      % assign new colormap
                UserData(:,1:3).*(UserData(:,4)*ones(1,3)))
            if fighandle ~= figh
                set(fighandle,'Colormap',UserData(:,1:3).*(UserData(:,4)*ones(1,3)))
            end

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif strcmp(bins,'motion')      % Mouse Motion Handler
            %global M_BUTTON_DOWN M_LINES M_PNTS M_TEXT
            if isempty(M_BUTTON_DOWN)
                return     % If mouse button is not down then return
            end
            axhandle = gca;
            nodes = get(gca,'UserData');
            tagval = get(gca,'tag');
            P = get(axhandle,'CurrentPoint');
            X = P(1,1);
            Y = P(1,2);
            nodenum = size(nodes,1);
            if M_BUTTON_DOWN == 1 | M_BUTTON_DOWN == nodenum % move selected node to
                nodes(M_BUTTON_DOWN,2) = max(min(Y,1),0);       % new mouse position
            else
                nodes(M_BUTTON_DOWN,1) = max(nodes(M_BUTTON_DOWN-1,1), ...
                    min(nodes(M_BUTTON_DOWN+1,1),round(X)));
                nodes(M_BUTTON_DOWN,2) = max(min(Y,1),0);
            end

            delete(M_PNTS(M_BUTTON_DOWN))
            M_PNTS(M_BUTTON_DOWN) = plot(nodes(M_BUTTON_DOWN,1),...
                nodes(M_BUTTON_DOWN,2),'ow');
            if M_BUTTON_DOWN > 1
                delete(M_LINES(M_BUTTON_DOWN-1))
                M_LINES(M_BUTTON_DOWN-1) = line([nodes(M_BUTTON_DOWN-1,1),...
                    nodes(M_BUTTON_DOWN,1)],...
                    [nodes(M_BUTTON_DOWN-1,2),nodes(M_BUTTON_DOWN,2)],'color',tagval);
            end
            if M_BUTTON_DOWN < nodenum
                delete(M_LINES(M_BUTTON_DOWN))
                M_LINES(M_BUTTON_DOWN) = line([nodes(M_BUTTON_DOWN,1),...
                    nodes(M_BUTTON_DOWN+1,1)],...
                    [nodes(M_BUTTON_DOWN,2),nodes(M_BUTTON_DOWN+1,2)],'color',tagval);
            end

            delete(M_TEXT)
            M_TEXT = text(0,1.05,...
                [' ',int2str(nodes(M_BUTTON_DOWN,1)),', ',...
                num2str(nodes(M_BUTTON_DOWN,2),2)],...
                'Color','w','clipping','off',...
                'VerticalAlignment','bottom',...
                'FontSize',10);

            axis([1 nodes(nodenum,1) 0 1])
            set(axhandle,'UserData',nodes)

        elseif strcmp(bins,'export')   % export colormap to desktop
            fig = gcbf;
            pos = get(fig,'Position');
            figx = 400;                  % set up dialog figure
            figy = 200;
            figh = figure_w_normalized_uicontrolunits('Units','pixels',...
                'Position',...
                [pos(1)+pos(3)/2-figx/2 pos(2)+pos(4)/2-figy/2 figx figy],...
                'Resize','off','CloseRequestFcn','',...
                'numbertitle','off');
            uicolor = get(figh,'Color');

            % text uicontrol
            uicontrol('Style','Text','Units','Pixels',...
                'String','Output Colormap to Desktop Variable:',...
                'Position',[20 figy-40 300 22],'HorizontalAlignment','left',...
                'FontSize',14,'BackGroundColor',uicolor)

            % edit uicontrol
            ui1 = uicontrol('Style','Edit','Units','Pixels',...
                'String','cmap','FontSize',12,...
                'Position',[120 figy-100 150 30]);

            % OK Button
            uicontrol('Style','PushButton','Units','Pixels',...
                'String','OK','FontSize',14,...
                'Position',[figx/4-20 10 65 30],...
                'UserData',[ui1 fig],'Callback',@timestringA) %does not get data back proerly, I'd bet

            % Cancel Button
            uicontrol('Style','PushButton','Units','Pixels',...
                'String','Cancel','FontSize',14,...
                'Position',[3*figx/4-20 10 65 30],...
                'Callback',@timestringB) %does not get data back proerly, I'd bet

        elseif strcmp(bins,'close')   % Close Request function
            figh = gcbf;
            clear global M_BUTTON_DOWN M_LINES M_PNTS M_TEXT % get rid of global vars
            delete(figh)

        elseif strcmp(bins,'reset')   % reset Colormap to original values
            fighandle = gcbf;
            obj = findobj('tag','hd','parent',fighandle);
            map = get(obj,'UserData');  % get original colormap stored in hd axes
            set(fighandle,'Colormap',map,...
                'UserData',[map , ones(size(map,1),1)])

            obj = findobj('tag','cb','parent',fighandle);
            nodes = get(obj,'UserData');  % get node values stored in cb axes
            red = nodes(4:nodes(1,1)+3,:);
            green = nodes(nodes(1,1)+4:nodes(1,1)+nodes(2,1)+3,:);
            blue = nodes(nodes(1,1)+nodes(2,1)+4:nodes(1,1)+nodes(2,1)+...
                nodes(3,1)+3,:);

            obj = findobj('tag','r','parent',fighandle);
            axes(obj)                         % reset colormap graphs
            cla
            nodenum = nodes(1,1);
            for i = 1:nodenum-1
                plot(red(i,1),red(i,2),'ow')
                line([red(i,1),red(i+1)],[red(i,2),red(i+1,2)],'color','r')
            end
            plot(red(nodenum,1),red(nodenum,2),'ow')
            set(obj,'UserData',red)

            obj = findobj('tag','g','parent',fighandle);
            axes(obj)
            cla
            nodenum = nodes(2,1);
            for i = 1:nodenum-1
                plot(green(i,1),green(i,2),'ow')
                line([green(i,1),green(i+1)],[green(i,2),green(i+1,2)],'color','g')
            end
            plot(green(nodenum,1),green(nodenum,2),'ow')
            set(obj,'UserData',green)

            obj = findobj('tag','b','parent',fighandle);
            axes(obj)
            cla
            nodenum = nodes(3,1);
            for i = 1:nodenum-1
                plot(blue(i,1),blue(i,2),'ow')
                line([blue(i,1),blue(i+1)],[blue(i,2),blue(i+1,2)],'color','b')
            end
            plot(blue(nodenum,1),blue(nodenum,2),'ow')
            set(obj,'UserData',blue)

            obj = findobj('tag','y','parent',fighandle);
            axes(obj)
            cla
            trans = [1 1;size(map,1) 1];
            nodenum = 2;
            for i = 1:nodenum-1
                plot(trans(i,1),trans(i,2),'ow')
                line([trans(i,1),trans(i+1)],[trans(i,2),trans(i+1,2)],'color','y')
            end
            plot(trans(nodenum,1),trans(nodenum,2),'ow')
            set(obj,'UserData',trans)
            figh = str2double(get(fighandle,'tag'));
            if figh~=fighandle
                set(figh,'Colormap',map)
            end
        end
    end

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        makecmap('reset');
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        makecmap('export');
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        makecmap('close');
    end

    function timestringA(src,ev)
        [OBJ1,FIGH1] = gcbo;
        OBJHAN = get(OBJ1,'UserData');
        LAB1 = get(OBJHAN(1),'string');
        DATA1 = get(OBJHAN(2),'Colormap');
        eval([LAB1,' = DATA1;']); % YUCK. Yuck. yuck.
        delete(FIGH1);
    end

    function timestringB(src,ev)
        [OBJ1,FIGH1] = gcbo;
        delete(FIGH1);
    end
end

author: CelsoReyes


text:
function parent = tree_rb_to_parent ( n, a )

%*****************************************************************************80
%
%% TREE_RB_TO_PARENT converts rooted binary tree to parent node representation.
%
%  Discussion:
%
%    Parent node representation of a tree assigns to each node a "parent" node,
%    which represents the first link of the path between the node and the
%    root node.  The root node itself is assigned a parent of 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of nodes in the tree.
%
%    Input, integer A(N), the preorder traversal form for the
%    rooted binary tree.
%
%    Output, integer PARENT(N), the parent node representation
%    of the tree.
%
  node = 0;
  node_num = 0;

  for k = 1 : n

    dad = node;
    node_num = node_num + 1;
    node = node_num;
    parent(node) = dad;

    if ( a(k) == 1 )

      use(node) = 0;

    else

      use(node) = 2;

      while ( use(node) == 2 )
        node = dad;
        if ( node == 0 )
          break
        end
        use(node) = use(node) + 1;
        dad = parent(node);
      end

    end

  end

  return
end


text:
function [range_h,range_w,reg_window]=init_regwindow(sz,target_sz,params)
        reg_scale =target_sz;
        use_sz = sz;
        reg_window = ones(use_sz) * params.reg_window_max;
        range = zeros(numel(reg_scale), 2);

        % determine the target center and range in the regularization windows
        for j = 1:numel(reg_scale)
            range(j,:) = [0, reg_scale(j) - 1] - floor(reg_scale(j) / 2);
        end
        center = floor((use_sz + 1)/ 2) + mod(use_sz + 1,2);
        range_h = (center(1)+ range(1,1)) : (center(1) + range(1,2));
        range_w = (center(2)+ range(2,1)) : (center(2) + range(2,2));
        reg_window(range_h, range_w) = params.reg_window_min;

repo: AutoTrack


text:
% h = drawellipse(trk,varargin)
% h = drawellipse(x,y,theta,semimaj,semimin,varargin)
% extra arguments will be fed to plot
function varargout = updateellipse(h,x,varargin)

% draw an isosceles triangle with center (x,y)
% with rotation theta
% with height maj*4
% with base min*4

if isstruct(x),
  fly = x;
  if nargin == 1 || ~isnaturalnumber(varargin{1}),
    x = [fly.x];
    y = [fly.y];
    theta = [fly.theta];
    a = [fly.a];
    b = [fly.b];
  else
    t = varargin{1};
    varargin = varargin(2:end);
    x = fly.x(t);
    y = fly.y(t);
    theta = fly.theta(t);
    a = fly.a(t);
    b = fly.b(t);
  end
else
  if nargin < 5,
    error('not enough arguments; usage: drawflyo(x,y,theta,a,b,...)');
  end
  y = varargin{1};
  theta = varargin{2};
  a = varargin{3};
  b = varargin{4};
  varargin = varargin(5:end);
end

phi = -0.03:0.01:2*pi;

nx = numel(x); ny = numel(y);
na = numel(a); nb = numel(b);
ntheta = numel(theta);
n = max([nx,ny,na,nb,ntheta]);

for i = 1:n,
  X1 = a(min(i,na))*cos(phi);
  Y1 = b(min(i,nb))*sin(phi);
  costheta = cos(theta(min(i,ntheta)));
  sintheta = sin(theta(min(i,ntheta)));
  X = costheta*X1 - sintheta*Y1 + x(min(i,nx));
  Y = sintheta*X1 + costheta*Y1 + y(min(i,ny));
  set(h(i),'XData',X,'YData',Y);
end

if nargout > 0,
  varargout{1} = h;
end;

author: kristinbranson

repo: JAABA


text:
function out=sum(x,dim)

precision=x(1).precision;
% can only handle up to 2 dimensions
s=size(x);
if nargin==1
 if any(s==1)
  out=mp(0,precision);
  for ii=1:numel(x)
   out=out+x(ii);
  end
 else
  out=mp(zeros(1,s(2)),precision);
  for j=1:s(2)
   for i=1:s(1)
    out(1,j)=out(1,j)+x(i,j);
   end
  end
 end
elseif nargin==2
 if dim==1
  out=mp(zeros(1,s(2)),precision);
  for j=1:s(2)
   for i=1:s(1)
    out(1,j)=out(1,j)+x(i,j);
   end
  end
 elseif dim==2
  out=mp(zeros(s(1),1),precision);
  for j=1:s(2)
   for i=1:s(1)
    out(i,1)=out(i,1)+x(i,j);
   end
  end
 else
  error('mp sum can only handle up to 2-D arrays (11/04)');
 end
end


text:
function K = compute_cubic_spline_kernel(x,y);
nx = size(x,1);
ny = size(y,1);
Kmax = max(repmat(abs(x),1,ny),repmat(abs(y)',nx,1) );
Kmin = min(repmat(abs(x),1,ny),repmat(abs(y)',nx,1) );
K = Kmin .* Kmin .* ( 3 * Kmax - Kmin );
K = K .* (  x*y' > 0 );

repo: REST


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function image = inRead(fileName)

image = im2double(imread(fileName));

end

author: mahmoudnafifi


text:
function [r,v] = keplerUniversal(r0,v0,t,mu)
%Purpose:
%Most effecient way to propagate any type of two body orbit using kepler's
%equations.
%-------------------------------------------------------------------------%
%                                                                         %
% Inputs:                                                                 %
%--------
%r_ECI                  [3 x N]                         Position Vector in
%                                                       ECI coordinate
%                                                       frame of reference
%
%v_ECI                  [3 x N]                         Velocity vector in
%                                                       ECI coordinate
%                                                       frame of reference
%
%t                      [1 x N]                         time vector in
%                                                       seconds
%
%mu                     double                          Gravitational Constant
%                                                       Defaults to Earth if
%                                                       not specified
% Outputs:
%---------                                                                %
%r_ECI                  [3 x N]                         Final position
%                                                       vector in ECI
%
%v_ECI                  [3 x N]                         Final velocity
%                                                       vector in ECI
%--------------------------------------------------------------------------
% Programmed by Darin Koblick 03-04-2012                                  %
%--------------------------------------------------------------------------
if ~exist('mu','var'); mu = 398600.4418; end
tol = 1e-9;
v0Mag = sqrt(sum(v0.^2,1));  r0Mag = sqrt(sum(r0.^2,1));
alpha = -(v0Mag.^2)./mu + 2./r0Mag;
%% Compute initial guess (X0) for Newton's Method
X0 = NaN(size(t));
%Check if there are any Eliptic/Circular orbits
idx = alpha > 0.000001;
if any(idx)
    X0(idx) = sqrt(mu).*t(idx).*alpha(idx);
end
%Check if there are any Parabolic orbits
idx = abs(alpha) < 0.000001;
if any(idx)
   h = cross(r0(:,idx),v0(:,idx)); hMag = sqrt(sum(h.^2,1));
   p = (hMag.^2)./mu; s = acot(3.*sqrt(mu./(p.^3)).*t(idx))./2;
   w = atan(tan(s).^(1/3)); X0(idx) = sqrt(p).*2.*cot(2.*w);
end
%Check if there are any Hyperbolic orbits
idx = alpha < -0.000001;
if any(idx)
   a = 1./alpha(idx);
   X0(idx) = sign(t(idx)).*sqrt(-a).*...
       log(-2.*mu.*alpha(idx).*t(idx)./ ...
       (dot(r0(:,idx),v0(:,idx))+sign(t(idx)).*sqrt(-mu.*a).*...
       (1-r0Mag(idx).*alpha(idx))));
end
%% Newton's Method to converge on solution
% Declare Constants that do not need to be computed within the while loop
err = Inf;
dr0v0Smu = dot(r0,v0)./sqrt(mu);
Smut = sqrt(mu).*t;
while any(abs(err) > tol)
    X02 = X0.^2;
    X03 = X02.*X0;
    psi = X02.*alpha;
    [c2,c3] = c2c3(psi);
    X0tOmPsiC3 = X0.*(1-psi.*c3);
    X02tC2 = X02.*c2;
    r = X02tC2 + dr0v0Smu.*X0tOmPsiC3 + r0Mag.*(1-psi.*c2);
    Xn = X0 + (Smut-X03.*c3-dr0v0Smu.*X02tC2-r0Mag.*X0tOmPsiC3)./r;
    err = Xn-X0; X0 = Xn;
end
f = 1 - (Xn.^2).*c2./r0Mag; g = t - (Xn.^3).*c3./sqrt(mu);
gdot = 1 - c2.*(Xn.^2)./r; fdot = Xn.*(psi.*c3-1).*sqrt(mu)./(r.*r0Mag);
r = bsxfun(@times,f,r0) + bsxfun(@times,g,v0);
v = bsxfun(@times,fdot,r0) + bsxfun(@times,gdot,v0);
%% Ensure Solution Integrity
%idx = round((f.*gdot - fdot.*g)./tol).*tol ~= 1; r(:,idx) = NaN; v(:,idx) = NaN;
end

function [c2,c3] = c2c3(psi)
%Vallado pg. 71 Algorithm 1
c2 = NaN(size(psi));
c3 = NaN(size(psi));
idx = psi > 1e-6;
if any(idx)
    c2(idx) = (1-cos(sqrt(psi(idx))))./psi(idx);
    c3(idx) = (sqrt(psi(idx))-sin(sqrt(psi(idx))))./sqrt(psi(idx).^3);
end
idx = psi < -1e-6;
if any(idx)
    c2(idx) = (1 - cosh(sqrt(-psi(idx))))./psi(idx);
    c3(idx) = (sinh(sqrt(-psi(idx)))-sqrt(-psi(idx)))./sqrt(-psi(idx).^3);
end
idx = abs(psi) <= 1e-6;
if any(idx)
    c2(idx) = 0.5;
    c3(idx) = 1/6;
end
end


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: This is where you can perform data analysis of a simulation
%           from your stored viz_IB2d information
%
%     Note:
%           (1) USER-DEFINED functions should be made made to analyze specific
%               data sets, regions of interest, etc. (see
%               EXAMPLE_FOR_DATA_ANALYSIS for an example of this)
%           (2) MUST make sure to 'addpath' to where DA_Blackbox is, i.e.,
%               line 56
%           (3) MUST make sure to set path to desired dataset, i.e., in line
%               59
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Insect_Analysis()

% TEMPORAL INFO FROM input2d %
dt = 2.5e-5;      % Time-step
Tfinal = 0.05;   % Final time in simulation
pDump=40;        % Note: 'print_dump' should match from input2d

% DATA ANALYSIS INFO %
start=1;                             % 1ST interval # included in data analysis
finish=25;                           % LAST interval # included in data analysis
dump_Times = (start:1:finish)*pDump; % Time vector when data was printed in analysis

% SET PATH TO DESIRED viz_IB2d DATA %
pathMain = '/Users/nick_battista/Desktop/IB2d/matIB2d/Examples/Example_Insect_Flight/Insect_1024x1024';

% SET PATH TO DA_BLACKBOX %
addpath('/Users/nick_battista/Desktop/IB2d/data_analysis/analysis_in_matlab/DA_Blackbox');

pathViz = [pathMain '/viz_IB2d/'];
pathForce = [pathMain '/hier_IB2d_data'];

%simVec = {'32x32','64x64','96x96','128x128','256x256','512x512','768x768','1024x1024'};

for i=start:1:finish

    % Points to desired data viz_IB2d data file
    if i<10
       numSim = ['000', num2str(i) ];
    elseif i<100
       numSim = ['00', num2str(i) ];
    elseif i<1000
       numSim = ['0', num2str(i)];
    else
       numSim = num2str(i);
    end

    % Imports immersed boundary positions %
    %[xLag,yLag] = give_Lag_Positions(pathViz,numSim);

    % Imports (x,y) grid values and ALL Eulerian Data %
    %                      DEFINITIONS
    %          x: x-grid                y: y-grid
    %       Omega: vorticity           P: pressure
    %    uMag: mag. of velocity
    %    uX: mag. of x-Velocity   uY: mag. of y-Velocity
    %    U: x-directed velocity   V: y-directed velocity
    %    Fx: x-directed Force     Fy: y-directed Force
    %
    %  Note: U(j,i): j-corresponds to y-index, i to the x-index
    %
    %[x,y,Omega,P,uMag,uX,uY,U,V,Fx,Fy] = import_Eulerian_Data(pathViz,numSim);

    % Imports Lagrangian Pt. FORCE (magnitude) DATA %
    %                      DEFINITIONS
    %
    %      fX_Lag: forces in x-direction on boundary
    %      fY_Lag: forces in y-direction on boundary
    %       fLagMag: magnitude of force at boundary
    %   fLagNorm: magnitude of NORMAL force at boundary
    %   fLagTan: magnitude of TANGENT force at boundary
    %
    [fX_Lag,fY_Lag] = import_Lagrangian_Force_Data_Insect(pathForce,numSim);

    fx_Mat(i,1) = mean( fX_Lag );
    fy_Mat(i,1) = mean( fY_Lag );

end

cd ..;

strName = 'fX_1024';
print_Matrix_To_Txt_File(fx_Mat,strName)

strName = 'fY_1024';
print_Matrix_To_Txt_File(fy_Mat,strName)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION prints matrix to file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Matrix_To_Txt_File(a,strName)

nameTxt = [strName '.txt'];

fid = fopen(nameTxt, 'wt'); % Open for writing
for i=1:size(a,1)
   fprintf(fid, '%d ', a(i,:));
   fprintf(fid, '\n');
end
fclose(fid);

author: nickabattista


text:
function [g,a]=ar1nv(x)
% AR1NV - Estimate the parameters for an AR(1) model
% Syntax: [g,a]=ar1nv(x);
%
% Input: x - a time series.
%
% Output: g - estimate of the lag-one autocorrelation.
%         a - estimate of the noise variance.

% (c) Eric Breitenberger

x=x(:);
N=length(x);
m=mean(x);
x=x-m;

% Lag zero and one covariance estimates:
c0=x'*x/N;
c1=x(1:N-1)'*x(2:N)/(N-1);

g=c1/c0;
a=sqrt((1-g^2)*c0);

author: grinsted

repo: wavelet-coherence


text:
% The COBRAToolbox: testSimpleOptKnock.m
%
% Purpose:
%     - test the simpleOptKnock function
%
% Authors:
%     - Jacek Wachowiak
global CBTDIR
% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testSimpleOptKnock'));
cd(fileDir);

% test variables
model = getDistributedModel('ecoli_core_model.mat');
targetRxn = model.rxns(20); % 'EX_ac(e)'

% function outputs
[wtRes, delRes] = simpleOptKnock(model, targetRxn);
[wtRes2, delRes2] = simpleOptKnock(model, targetRxn, model.genes(1), 1, 0.5, 1);

% tests
assert(abs(12.1585 - delRes.maxProd(23)) < 1e-4); % 'EX_co2(e)'
assert(abs(14.3123 - delRes.maxProd(12)) < 1e-4);% 'ATPS4r'
assert(abs(12.1584 - delRes.minProd(14)) < 1e-4);% 'CO2t'
assert(abs(8.5036 - delRes.minProd(16)) < 1e-4);% 'CYTBD'
assert(abs(0.8739 - wtRes.growth) < 1e-4);% always the same
%for gene deletions
assert(abs(0.8739 - wtRes2.growth) < 1e-4);% always the same
assert(abs(3e-07 - delRes2.maxProd) < 1e-4);
assert(abs(delRes2.minProd) < 1e-4);

% change to old directory
cd(currentDir);


text:
function [W,H] = nmfsc( V, rdim, sW, sH, fname, showflag )
% nmfsc - non-negative matrix factorization with sparseness constraints
%
% SYNTAX:
% [W,H] = nmfsc( V, rdim, sW, sH, fname, showflag );
%
% INPUTS:
% V          - data matrix
% rdim       - number of components (inner dimension of factorization)
% sW         - sparseness of W, in [0,1]. (give [] if no constraint)
% sH         - sparseness of H, in [0,1]. (give [] if no constraint)
% fname      - name of file to write results into
% showflag   - binary flag. if set then graphically show progress
%
% Note: Sparseness is measured on the scale [0,1] where 0 means
% completely distributed and 1 means ultimate sparseness.
%
% NOTE: There is NO CONVERGENCE CRITERION. The estimation never ends,
% but rather has to be terminated manually. See README file of code
% package for details.
%

% Check that we have non-negative data
if min(V(:))<0, error('Negative values in data!'); end

% Globally rescale data to avoid potential overflow/underflow
V = V/max(V(:));

% Data dimensions
vdim = size(V,1);
samples = size(V,2);

% Create initial matrices
W = abs(randn(vdim,rdim));
H = abs(randn(rdim,samples));
H = H./(sqrt(sum(H.^2,2))*ones(1,samples));

% Make initial matrices have correct sparseness
if ~isempty(sW),
    L1a = sqrt(vdim)-(sqrt(vdim)-1)*sW;
    for i=1:rdim, W(:,i) = projfunc(W(:,i),L1a,1,1); end
end
if ~isempty(sH),
    L1s = sqrt(samples)-(sqrt(samples)-1)*sH;
    for i=1:rdim, H(i,:) = (projfunc(H(i,:)',L1s,1,1))'; end
end

% Initialize displays
if showflag,
   figure(1); clf; % this will show the energies and sparsenesses
   figure(2); clf; % this will show the objective function
   drawnow;
end

% Calculate initial objective
objhistory = 0.5*sum(sum((V-W*H).^2));

% Initial stepsizes
stepsizeW = 1;
stepsizeH = 1;

timestarted = clock;

% Start iteration
iter = 0;
while 1,

    % Show progress
    fprintf('[%d]: %.5f\n',iter,objhistory(end));

    % Save every once in a while
    if rem(iter,5)==0,
	elapsed = etime(clock,timestarted);
	fprintf('Saving...');
	save(fname,'W','H','sW','sH','iter','objhistory','elapsed');
	fprintf('Done!\n');
    end

    % Show stats
    if showflag & (rem(iter,5)==0),
	figure(1);
	subplot(3,1,1); bar(sqrt(sum(W.^2)).*sqrt(sum(H'.^2)));
	cursW = (sqrt(vdim)-(sum(abs(W))./sqrt(sum(W.^2))))/(sqrt(vdim)-1);
	subplot(3,1,2); bar(cursW);
	cursH = (sqrt(samples)-(sum(abs(H'))./sqrt(sum(H'.^2)))) ...
		/(sqrt(samples)-1);
	subplot(3,1,3); bar(cursH);
	if iter>1,
	    figure(2);
	    plot(objhistory(2:end));
    end
    % added for now
    figure(100); imstiled(reshape(W,15,15,43),[],'gray')
	drawnow;
    end

    % Update iteration count
    iter = iter+1;

    % Save old values
    Wold = W;
    Hold = H;

    % ----- Update H ---------------------------------------

    if ~isempty(sH),

	% Gradient for H
	dH = W'*(W*H-V);
	begobj = objhistory(end);

	% Make sure we decrease the objective!
	while 1,

	    % Take step in direction of negative gradient, and project
	    Hnew = H - stepsizeH*dH;
	    for i=1:rdim, Hnew(i,:) = (projfunc(Hnew(i,:)',L1s,1,1))'; end

	    % Calculate new objective
	    newobj = 0.5*sum(sum((V-W*Hnew).^2));

	    % If the objective decreased, we can continue...
	    if newobj<=begobj,
		break;
	    end

	    % ...else decrease stepsize and try again
	    stepsizeH = stepsizeH/2;
	    fprintf('.');
	    if stepsizeH<1e-200,
		fprintf('Algorithm converged.\n');
		return;
	    end

	end

	% Slightly increase the stepsize
	stepsizeH = stepsizeH*1.2;
	H = Hnew;

    else

	% Update using standard NMF multiplicative update rule
	H = H.*(W'*V)./(W'*W*H + 1e-9);

	% Renormalize so rows of H have constant energy
	norms = sqrt(sum(H'.^2));
	H = H./(norms'*ones(1,samples));
	W = W.*(ones(vdim,1)*norms);

    end

    % ----- Update W ---------------------------------------

    if ~isempty(sW),

	% Gradient for W
	dW = (W*H-V)*H';
	begobj = 0.5*sum(sum((V-W*H).^2));

	% Make sure we decrease the objective!
	while 1,

	    % Take step in direction of negative gradient, and project
	    Wnew = W - stepsizeW*dW;
	    norms = sqrt(sum(Wnew.^2));
	    for i=1:rdim,
		Wnew(:,i) = projfunc(Wnew(:,i),L1a*norms(i),(norms(i)^2),1);
	    end

	    % Calculate new objective
	    newobj = 0.5*sum(sum((V-Wnew*H).^2));

	    % If the objective decreased, we can continue...
	    if newobj<=begobj,
		break;
	    end

	    % ...else decrease stepsize and try again
	    stepsizeW = stepsizeW/2;
	    fprintf(',');
	    if stepsizeW<1e-200,
		fprintf('Algorithm converged.\n');
		return;
	    end

	end

	% Slightly increase the stepsize
	stepsizeW = stepsizeW*1.2;
	W = Wnew;

    else

	% Update using standard NMF multiplicative update rule
	W = W.*(V*H')./(W*H*H' + 1e-9);

    end

    % Calculate objective
    newobj = 0.5*sum(sum((V-W*H).^2));
    objhistory = [objhistory newobj];

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes,varargin)
% gridfit: estimates a surface on a 2d grid, based on scattered data
%          Replicates are allowed. All methods extrapolate to the grid
%          boundaries. Gridfit uses a modified ridge estimator to
%          generate the surface, where the bias is toward smoothness.
%
%          Gridfit is not an interpolant. Its goal is a smooth surface
%          that approximates your data, but allows you to control the
%          amount of smoothing.
%
% usage #1: zgrid = gridfit(x,y,z,xnodes,ynodes);
% usage #2: [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes);
% usage #3: zgrid = gridfit(x,y,z,xnodes,ynodes,prop,val,prop,val,...);
%
% Arguments: (input)
%  x,y,z - vectors of equal lengths, containing arbitrary scattered data
%          The only constraint on x and y is they cannot ALL fall on a
%          single line in the x-y plane. Replicate points will be treated
%          in a least squares sense.
%
%          ANY points containing a NaN are ignored in the estimation
%
%  xnodes - vector defining the nodes in the grid in the independent
%          variable (x). xnodes need not be equally spaced. xnodes
%          must completely span the data. If they do not, then the
%          'extend' property is applied, adjusting the first and last
%          nodes to be extended as necessary. See below for a complete
%          description of the 'extend' property.
%
%          If xnodes is a scalar integer, then it specifies the number
%          of equally spaced nodes between the min and max of the data.
%
%  ynodes - vector defining the nodes in the grid in the independent
%          variable (y). ynodes need not be equally spaced.
%
%          If ynodes is a scalar integer, then it specifies the number
%          of equally spaced nodes between the min and max of the data.
%
%          Also see the extend property.
%
%  Additional arguments follow in the form of property/value pairs.
%  Valid properties are:
%    'smoothness', 'interp', 'regularizer', 'solver', 'maxiter'
%    'extend', 'tilesize', 'overlap'
%
%  Any UNAMBIGUOUS shortening (even down to a single letter) is
%  valid for property names. All properties have default values,
%  chosen (I hope) to give a reasonable result out of the box.
%
%   'smoothness' - scalar or vector of length 2 - determines the
%          eventual smoothness of the estimated surface. A larger
%          value here means the surface will be smoother. Smoothness
%          must be a non-negative real number.
%
%          If this parameter is a vector of length 2, then it defines
%          the relative smoothing to be associated with the x and y
%          variables. This allows the user to apply a different amount
%          of smoothing in the x dimension compared to the y dimension.
%
%          Note: the problem is normalized in advance so that a
%          smoothness of 1 MAY generate reasonable results. If you
%          find the result is too smooth, then use a smaller value
%          for this parameter. Likewise, bumpy surfaces suggest use
%          of a larger value. (Sometimes, use of an iterative solver
%          with too small a limit on the maximum number of iterations
%          will result in non-convergence.)
%
%          DEFAULT: 1
%
%
%   'interp' - character, denotes the interpolation scheme used
%          to interpolate the data.
%
%          DEFAULT: 'triangle'
%
%          'bilinear' - use bilinear interpolation within the grid
%                     (also known as tensor product linear interpolation)
%
%          'triangle' - split each cell in the grid into a triangle,
%                     then linear interpolation inside each triangle
%
%          'nearest' - nearest neighbor interpolation. This will
%                     rarely be a good choice, but I included it
%                     as an option for completeness.
%
%
%   'regularizer' - character flag, denotes the regularization
%          paradignm to be used. There are currently three options.
%
%          DEFAULT: 'gradient'
%
%          'diffusion' or 'laplacian' - uses a finite difference
%              approximation to the Laplacian operator (i.e, del^2).
%
%              We can think of the surface as a plate, wherein the
%              bending rigidity of the plate is specified by the user
%              as a number relative to the importance of fidelity to
%              the data. A stiffer plate will result in a smoother
%              surface overall, but fit the data less well. I've
%              modeled a simple plate using the Laplacian, del^2. (A
%              projected enhancement is to do a better job with the
%              plate equations.)
%
%              We can also view the regularizer as a diffusion problem,
%              where the relative thermal conductivity is supplied.
%              Here interpolation is seen as a problem of finding the
%              steady temperature profile in an object, given a set of
%              points held at a fixed temperature. Extrapolation will
%              be linear. Both paradigms are appropriate for a Laplacian
%              regularizer.
%
%          'gradient' - attempts to ensure the gradient is as smooth
%              as possible everywhere. Its subtly different from the
%              'diffusion' option, in that here the directional
%              derivatives are biased to be smooth across cell
%              boundaries in the grid.
%
%              The gradient option uncouples the terms in the Laplacian.
%              Think of it as two coupled PDEs instead of one PDE. Why
%              are they different at all? The terms in the Laplacian
%              can balance each other.
%
%          'springs' - uses a spring model connecting nodes to each
%              other, as well as connecting data points to the nodes
%              in the grid. This choice will cause any extrapolation
%              to be as constant as possible.
%
%              Here the smoothing parameter is the relative stiffness
%              of the springs connecting the nodes to each other compared
%              to the stiffness of a spting connecting the lattice to
%              each data point. Since all springs have a rest length
%              (length at which the spring has zero potential energy)
%              of zero, any extrapolation will be minimized.
%
%          Note: The 'springs' regularizer tends to drag the surface
%          towards the mean of all the data, so too large a smoothing
%          parameter may be a problem.
%
%
%   'solver' - character flag - denotes the solver used for the
%          resulting linear system. Different solvers will have
%          different solution times depending upon the specific
%          problem to be solved. Up to a certain size grid, the
%          direct \ solver will often be speedy, until memory
%          swaps causes problems.
%
%          What solver should you use? Problems with a significant
%          amount of extrapolation should avoid lsqr. \ may be
%          best numerically for small smoothnesss parameters and
%          high extents of extrapolation.
%
%          Large numbers of points will slow down the direct
%          \, but when applied to the normal equations, \ can be
%          quite fast. Since the equations generated by these
%          methods will tend to be well conditioned, the normal
%          equations are not a bad choice of method to use. Beware
%          when a small smoothing parameter is used, since this will
%          make the equations less well conditioned.
%
%          DEFAULT: 'normal'
%
%          '\' - uses matlab's backslash operator to solve the sparse
%                     system. 'backslash' is an alternate name.
%
%          'symmlq' - uses matlab's iterative symmlq solver
%
%          'lsqr' - uses matlab's iterative lsqr solver
%
%          'normal' - uses \ to solve the normal equations.
%
%
%   'maxiter' - only applies to iterative solvers - defines the
%          maximum number of iterations for an iterative solver
%
%          DEFAULT: min(10000,length(xnodes)*length(ynodes))
%
%
%   'extend' - character flag - controls whether the first and last
%          nodes in each dimension are allowed to be adjusted to
%          bound the data, and whether the user will be warned if
%          this was deemed necessary to happen.
%
%          DEFAULT: 'warning'
%
%          'warning' - Adjust the first and/or last node in
%                     x or y if the nodes do not FULLY contain
%                     the data. Issue a warning message to this
%                     effect, telling the amount of adjustment
%                     applied.
%
%          'never'  - Issue an error message when the nodes do
%                     not absolutely contain the data.
%
%          'always' - automatically adjust the first and last
%                     nodes in each dimension if necessary.
%                     No warning is given when this option is set.
%
%
%   'tilesize' - grids which are simply too large to solve for
%          in one single estimation step can be built as a set
%          of tiles. For example, a 1000x1000 grid will require
%          the estimation of 1e6 unknowns. This is likely to
%          require more memory (and time) than you have available.
%          But if your data is dense enough, then you can model
%          it locally using smaller tiles of the grid.
%
%          My recommendation for a reasonable tilesize is
%          roughly 100 to 200. Tiles of this size take only
%          a few seconds to solve normally, so the entire grid
%          can be modeled in a finite amount of time. The minimum
%          tilesize can never be less than 3, although even this
%          size tile is so small as to be ridiculous.
%
%          If your data is so sparse than some tiles contain
%          insufficient data to model, then those tiles will
%          be left as NaNs.
%
%          DEFAULT: inf
%
%
%   'overlap' - Tiles in a grid have some overlap, so they
%          can minimize any problems along the edge of a tile.
%          In this overlapped region, the grid is built using a
%          bi-linear combination of the overlapping tiles.
%
%          The overlap is specified as a fraction of the tile
%          size, so an overlap of 0.20 means there will be a 20%
%          overlap of successive tiles. I do allow a zero overlap,
%          but it must be no more than 1/2.
%
%          0 <= overlap <= 0.5
%
%          Overlap is ignored if the tilesize is greater than the
%          number of nodes in both directions.
%
%          DEFAULT: 0.20
%
%
%   'autoscale' - Some data may have widely different scales on
%          the respective x and y axes. If this happens, then
%          the regularization may experience difficulties.
%
%          autoscale = 'on' will cause gridfit to scale the x
%          and y node intervals to a unit length. This should
%          improve the regularization procedure. The scaling is
%          purely internal.
%
%          autoscale = 'off' will disable automatic scaling
%
%          DEFAULT: 'on'
%
%
% Arguments: (output)
%  zgrid   - (nx,ny) array containing the fitted surface
%
%  xgrid, ygrid - as returned by meshgrid(xnodes,ynodes)
%
%
% Speed considerations:
%  Remember that gridfit must solve a LARGE system of linear
%  equations. There will be as many unknowns as the total
%  number of nodes in the final lattice. While these equations
%  may be sparse, solving a system of 10000 equations may take
%  a second or so. Very large problems may benefit from the
%  iterative solvers or from tiling.
%
%
% Example usage:
%
%  x = rand(100,1);
%  y = rand(100,1);
%  z = exp(x+2*y);
%  xnodes = 0:.1:1;
%  ynodes = 0:.1:1;
%
%  g = gridfit(x,y,z,xnodes,ynodes);
%
% Note: this is equivalent to the following call:
%
%  g = gridfit(x,y,z,xnodes,ynodes, ...
%              'smooth',1, ...
%              'interp','triangle', ...
%              'solver','normal', ...
%              'regularizer','gradient', ...
%              'extend','warning', ...
%              'tilesize',inf);
%
%
% Author: John D'Errico
% e-mail address: woodchips@rochester.rr.com
% Release: 2.0
% Release date: 5/23/06

% set defaults
params.smoothness = 1;
params.interp = 'triangle';
params.regularizer = 'gradient';
params.solver = 'backslash';
params.maxiter = [];
params.extend = 'warning';
params.tilesize = inf;
params.overlap = 0.20;
params.mask = [];
params.autoscale = 'on';
params.xscale = 1;
params.yscale = 1;

% was the params struct supplied?
if ~isempty(varargin)
  if isstruct(varargin{1})
    % params is only supplied if its a call from tiled_gridfit
    params = varargin{1};
    if length(varargin)>1
      % check for any overrides
      params = parse_pv_pairs(params,varargin{2:end});
    end
  else
    % check for any overrides of the defaults
    params = parse_pv_pairs(params,varargin);

  end
end

% check the parameters for acceptability
params = check_params(params);

% ensure all of x,y,z,xnodes,ynodes are column vectors,
% also drop any NaN data
x=x(:);
y=y(:);
z=z(:);
k = isnan(x) | isnan(y) | isnan(z);
if any(k)
  x(k)=[];
  y(k)=[];
  z(k)=[];
end
xmin = min(x);
xmax = max(x);
ymin = min(y);
ymax = max(y);

% did they supply a scalar for the nodes?
if length(xnodes)==1
  xnodes = linspace(xmin,xmax,xnodes)';
  xnodes(end) = xmax; % make sure it hits the max
end
if length(ynodes)==1
  ynodes = linspace(ymin,ymax,ynodes)';
  ynodes(end) = ymax; % make sure it hits the max
end

xnodes=xnodes(:);
ynodes=ynodes(:);
dx = diff(xnodes);
dy = diff(ynodes);
nx = length(xnodes);
ny = length(ynodes);
ngrid = nx*ny;

% set the scaling if autoscale was on
if strcmpi(params.autoscale,'on')
  params.xscale = mean(dx);
  params.yscale = mean(dy);
  params.autoscale = 'off';
end

% check to see if any tiling is necessary
if (params.tilesize < max(nx,ny))
  % split it into smaller tiles. compute zgrid and ygrid
  % at the very end if requested
  zgrid = tiled_gridfit(x,y,z,xnodes,ynodes,params);
else
  % its a single tile.

  % mask must be either an empty array, or a boolean
  % aray of the same size as the final grid.
  nmask = size(params.mask);
  if ~isempty(params.mask) && ((nmask(2)~=nx) || (nmask(1)~=ny))
    if ((nmask(2)==ny) || (nmask(1)==nx))
      error 'Mask array is probably transposed from proper orientation.'
    else
      error 'Mask array must be the same size as the final grid.'
    end
  end
  if ~isempty(params.mask)
    params.maskflag = 1;
  else
    params.maskflag = 0;
  end

  % default for maxiter?
  if isempty(params.maxiter)
    params.maxiter = min(10000,nx*ny);
  end

  % check lengths of the data
  n = length(x);
  if (length(y)~=n) || (length(z)~=n)
    error 'Data vectors are incompatible in size.'
  end
  if n<3
    error 'Insufficient data for surface estimation.'
  end

  % verify the nodes are distinct
  if any(diff(xnodes)<=0) || any(diff(ynodes)<=0)
    error 'xnodes and ynodes must be monotone increasing'
  end

  % do we need to tweak the first or last node in x or y?
  if xmin<xnodes(1)
    switch params.extend
      case 'always'
        xnodes(1) = xmin;
      case 'warning'
        warning('GRIDFIT:extend',['xnodes(1) was decreased by: ',num2str(xnodes(1)-xmin),', new node = ',num2str(xmin)])
        xnodes(1) = xmin;
      case 'never'
        error(['Some x (',num2str(xmin),') falls below xnodes(1) by: ',num2str(xnodes(1)-xmin)])
    end
  end
  if xmax>xnodes(end)
    switch params.extend
      case 'always'
        xnodes(end) = xmax;
      case 'warning'
        warning('GRIDFIT:extend',['xnodes(end) was increased by: ',num2str(xmax-xnodes(end)),', new node = ',num2str(xmax)])
        xnodes(end) = xmax;
      case 'never'
        error(['Some x (',num2str(xmax),') falls above xnodes(end) by: ',num2str(xmax-xnodes(end))])
    end
  end
  if ymin<ynodes(1)
    switch params.extend
      case 'always'
        ynodes(1) = ymin;
      case 'warning'
        warning('GRIDFIT:extend',['ynodes(1) was decreased by: ',num2str(ynodes(1)-ymin),', new node = ',num2str(ymin)])
        ynodes(1) = ymin;
      case 'never'
        error(['Some y (',num2str(ymin),') falls below ynodes(1) by: ',num2str(ynodes(1)-ymin)])
    end
  end
  if ymax>ynodes(end)
    switch params.extend
      case 'always'
        ynodes(end) = ymax;
      case 'warning'
        warning('GRIDFIT:extend',['ynodes(end) was increased by: ',num2str(ymax-ynodes(end)),', new node = ',num2str(ymax)])
        ynodes(end) = ymax;
      case 'never'
        error(['Some y (',num2str(ymax),') falls above ynodes(end) by: ',num2str(ymax-ynodes(end))])
    end
  end

  % determine which cell in the array each point lies in
  [junk,indx] = histc(x,xnodes); %#ok
  [junk,indy] = histc(y,ynodes); %#ok
  % any point falling at the last node is taken to be
  % inside the last cell in x or y.
  k=(indx==nx);
  indx(k)=indx(k)-1;
  k=(indy==ny);
  indy(k)=indy(k)-1;
  ind = indy + ny*(indx-1);

  % Do we have a mask to apply?
  if params.maskflag
    % if we do, then we need to ensure that every
    % cell with at least one data point also has at
    % least all of its corners unmasked.
    params.mask(ind) = 1;
    params.mask(ind+1) = 1;
    params.mask(ind+ny) = 1;
    params.mask(ind+ny+1) = 1;
  end

  % interpolation equations for each point
  tx = min(1,max(0,(x - xnodes(indx))./dx(indx)));
  ty = min(1,max(0,(y - ynodes(indy))./dy(indy)));
  % Future enhancement: add cubic interpolant
  switch params.interp
    case 'triangle'
      % linear interpolation inside each triangle
      k = (tx > ty);
      L = ones(n,1);
      L(k) = ny;

      t1 = min(tx,ty);
      t2 = max(tx,ty);
      A = sparse(repmat((1:n)',1,3),[ind,ind+ny+1,ind+L], ...
        [1-t2,t1,t2-t1],n,ngrid);

    case 'nearest'
      % nearest neighbor interpolation in a cell
      k = round(1-ty) + round(1-tx)*ny;
      A = sparse((1:n)',ind+k,ones(n,1),n,ngrid);

    case 'bilinear'
      % bilinear interpolation in a cell
      A = sparse(repmat((1:n)',1,4),[ind,ind+1,ind+ny,ind+ny+1], ...
        [(1-tx).*(1-ty), (1-tx).*ty, tx.*(1-ty), tx.*ty], ...
        n,ngrid);

  end
  rhs = z;

  % do we have relative smoothing parameters?
  if numel(params.smoothness) == 1
    % it was scalar, so treat both dimensions equally
    smoothparam = params.smoothness;
    xyRelativeStiffness = [1;1];
  else
    % It was a vector, so anisotropy reigns.
    % I've already checked that the vector was of length 2
    smoothparam = sqrt(prod(params.smoothness));
    xyRelativeStiffness = params.smoothness(:)./smoothparam;
  end

  % Build regularizer. Add del^4 regularizer one day.
  switch params.regularizer
    case 'springs'
      % zero "rest length" springs
      [i,j] = meshgrid(1:nx,1:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      m = nx*(ny-1);
      stiffness = 1./(dy/params.yscale);
      Areg = sparse(repmat((1:m)',1,2),[ind,ind+1], ...
        xyRelativeStiffness(2)*stiffness(j(:))*[-1 1], ...
        m,ngrid);

      [i,j] = meshgrid(1:(nx-1),1:ny);
      ind = j(:) + ny*(i(:)-1);
      m = (nx-1)*ny;
      stiffness = 1./(dx/params.xscale);
      Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny], ...
        xyRelativeStiffness(1)*stiffness(i(:))*[-1 1],m,ngrid)];

      [i,j] = meshgrid(1:(nx-1),1:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      m = (nx-1)*(ny-1);
      stiffness = 1./sqrt((dx(i(:))/params.xscale/xyRelativeStiffness(1)).^2 + ...
        (dy(j(:))/params.yscale/xyRelativeStiffness(2)).^2);

      Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny+1], ...
        stiffness*[-1 1],m,ngrid)];

      Areg = [Areg;sparse(repmat((1:m)',1,2),[ind+1,ind+ny], ...
        stiffness*[-1 1],m,ngrid)];

    case {'diffusion' 'laplacian'}
      % thermal diffusion using Laplacian (del^2)
      [i,j] = meshgrid(1:nx,2:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      dy1 = dy(j(:)-1)/params.yscale;
      dy2 = dy(j(:))/params.yscale;

      Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
        xyRelativeStiffness(2)*[-2./(dy1.*(dy1+dy2)), ...
        2./(dy1.*dy2), -2./(dy2.*(dy1+dy2))],ngrid,ngrid);

      [i,j] = meshgrid(2:(nx-1),1:ny);
      ind = j(:) + ny*(i(:)-1);
      dx1 = dx(i(:)-1)/params.xscale;
      dx2 = dx(i(:))/params.xscale;

      Areg = Areg + sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], ...
        xyRelativeStiffness(1)*[-2./(dx1.*(dx1+dx2)), ...
        2./(dx1.*dx2), -2./(dx2.*(dx1+dx2))],ngrid,ngrid);

    case 'gradient'
      % Subtly different from the Laplacian. A point for future
      % enhancement is to do it better for the triangle interpolation
      % case.
      [i,j] = meshgrid(1:nx,2:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      dy1 = dy(j(:)-1)/params.yscale;
      dy2 = dy(j(:))/params.yscale;

      Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
        xyRelativeStiffness(2)*[-2./(dy1.*(dy1+dy2)), ...
        2./(dy1.*dy2), -2./(dy2.*(dy1+dy2))],ngrid,ngrid);

      [i,j] = meshgrid(2:(nx-1),1:ny);
      ind = j(:) + ny*(i(:)-1);
      dx1 = dx(i(:)-1)/params.xscale;
      dx2 = dx(i(:))/params.xscale;

      Areg = [Areg;sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], ...
        xyRelativeStiffness(1)*[-2./(dx1.*(dx1+dx2)), ...
        2./(dx1.*dx2), -2./(dx2.*(dx1+dx2))],ngrid,ngrid)];

  end
  nreg = size(Areg,1);

  % Append the regularizer to the interpolation equations,
  % scaling the problem first. Use the 1-norm for speed.
  NA = norm(A,1);
  NR = norm(Areg,1);
  A = [A;Areg*(smoothparam*NA/NR)];
  rhs = [rhs;zeros(nreg,1)];
  % do we have a mask to apply?
  if params.maskflag
    unmasked = find(params.mask);
  end
  % solve the full system, with regularizer attached
  switch params.solver
    case {'\' 'backslash'}
      if params.maskflag
        % there is a mask to use
        zgrid=nan(ny,nx);
        zgrid(unmasked) = A(:,unmasked)\rhs;
      else
        % no mask
        zgrid = reshape(A\rhs,ny,nx);
      end

    case 'normal'
      % The normal equations, solved with \. Can be faster
      % for huge numbers of data points, but reasonably
      % sized grids. The regularizer makes A well conditioned
      % so the normal equations are not a terribly bad thing
      % here.
      if params.maskflag
        % there is a mask to use
        Aunmasked = A(:,unmasked);
        zgrid=nan(ny,nx);
        zgrid(unmasked) = (Aunmasked'*Aunmasked)\(Aunmasked'*rhs);
      else
        zgrid = reshape((A'*A)\(A'*rhs),ny,nx);
      end

    case 'symmlq'
      % iterative solver - symmlq - requires a symmetric matrix,
      % so use it to solve the normal equations. No preconditioner.
      tol = abs(max(z)-min(z))*1.e-13;
      if params.maskflag
        % there is a mask to use
        zgrid=nan(ny,nx);
        [zgrid(unmasked),flag] = symmlq(A(:,unmasked)'*A(:,unmasked), ...
          A(:,unmasked)'*rhs,tol,params.maxiter);
      else
        [zgrid,flag] = symmlq(A'*A,A'*rhs,tol,params.maxiter);
        zgrid = reshape(zgrid,ny,nx);
      end
      % display a warning if convergence problems
      switch flag
        case 0
          % no problems with convergence
        case 1
          % SYMMLQ iterated MAXIT times but did not converge.
          warning('GRIDFIT:solver',['Symmlq performed ',num2str(params.maxiter), ...
            ' iterations but did not converge.'])
        case 3
          % SYMMLQ stagnated, successive iterates were the same
          warning('GRIDFIT:solver','Symmlq stagnated without apparent convergence.')
        otherwise
          warning('GRIDFIT:solver',['One of the scalar quantities calculated in',...
            ' symmlq was too small or too large to continue computing.'])
      end

    case 'lsqr'
      % iterative solver - lsqr. No preconditioner here.
      tol = abs(max(z)-min(z))*1.e-13;
      if params.maskflag
        % there is a mask to use
        zgrid=nan(ny,nx);
        [zgrid(unmasked),flag] = lsqr(A(:,unmasked),rhs,tol,params.maxiter);
      else
        [zgrid,flag] = lsqr(A,rhs,tol,params.maxiter);
        zgrid = reshape(zgrid,ny,nx);
      end

      % display a warning if convergence problems
      switch flag
        case 0
          % no problems with convergence
        case 1
          % lsqr iterated MAXIT times but did not converge.
          warning('GRIDFIT:solver',['Lsqr performed ', ...
            num2str(params.maxiter),' iterations but did not converge.'])
        case 3
          % lsqr stagnated, successive iterates were the same
          warning('GRIDFIT:solver','Lsqr stagnated without apparent convergence.')
        case 4
          warning('GRIDFIT:solver',['One of the scalar quantities calculated in',...
            ' LSQR was too small or too large to continue computing.'])
      end

  end  % switch params.solver

end  % if params.tilesize...

% only generate xgrid and ygrid if requested.
if nargout>1
  [xgrid,ygrid]=meshgrid(xnodes,ynodes);
end

% ============================================
% End of main function - gridfit
% ============================================

% ============================================
% subfunction - parse_pv_pairs
% ============================================
function params=parse_pv_pairs(params,pv_pairs)
% parse_pv_pairs: parses sets of property value pairs, allows defaults
% usage: params=parse_pv_pairs(default_params,pv_pairs)
%
% arguments: (input)
%  default_params - structure, with one field for every potential
%             property/value pair. Each field will contain the default
%             value for that property. If no default is supplied for a
%             given property, then that field must be empty.
%
%  pv_array - cell array of property/value pairs.
%             Case is ignored when comparing properties to the list
%             of field names. Also, any unambiguous shortening of a
%             field/property name is allowed.
%
% arguments: (output)
%  params   - parameter struct that reflects any updated property/value
%             pairs in the pv_array.
%
% Example usage:
% First, set default values for the parameters. Assume we
% have four parameters that we wish to use optionally in
% the function examplefun.
%
%  - 'viscosity', which will have a default value of 1
%  - 'volume', which will default to 1
%  - 'pie' - which will have default value 3.141592653589793
%  - 'description' - a text field, left empty by default
%
% The first argument to examplefun is one which will always be
% supplied.
%
%   function examplefun(dummyarg1,varargin)
%   params.Viscosity = 1;
%   params.Volume = 1;
%   params.Pie = 3.141592653589793
%
%   params.Description = '';
%   params=parse_pv_pairs(params,varargin);
%   params
%
% Use examplefun, overriding the defaults for 'pie', 'viscosity'
% and 'description'. The 'volume' parameter is left at its default.
%
%   examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world')
%
% params =
%     Viscosity: 10
%        Volume: 1
%           Pie: 3
%   Description: 'Hello world'
%
% Note that capitalization was ignored, and the property 'viscosity'
% was truncated as supplied. Also note that the order the pairs were
% supplied was arbitrary.

npv = length(pv_pairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  return
end

if ~isstruct(params)
  error 'No structure for defaults was supplied'
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
  p_i = lower(pv_pairs{2*i-1});
  v_i = pv_pairs{2*i};

  ind = strmatch(p_i,lpropnames,'exact');
  if isempty(ind)
    ind = find(strncmp(p_i,lpropnames,length(p_i)));
    if isempty(ind)
      error(['No matching property found for: ',pv_pairs{2*i-1}])
    elseif length(ind)>1
      error(['Ambiguous property name: ',pv_pairs{2*i-1}])
    end
  end
  p_i = propnames{ind};

  % override the corresponding default in params
  params = setfield(params,p_i,v_i); %#ok

end

% ============================================
% subfunction - check_params
% ============================================
function params = check_params(params)

% check the parameters for acceptability
% smoothness == 1 by default
if isempty(params.smoothness)
  params.smoothness = 1;
else
  if (numel(params.smoothness)>2) || any(params.smoothness<=0)
    error 'Smoothness must be scalar (or length 2 vector), real, finite, and positive.'
  end
end

% regularizer  - must be one of 4 options - the second and
% third are actually synonyms.
valid = {'springs', 'diffusion', 'laplacian', 'gradient'};
if isempty(params.regularizer)
  params.regularizer = 'diffusion';
end
ind = find(strncmpi(params.regularizer,valid,length(params.regularizer)));
if (length(ind)==1)
  params.regularizer = valid{ind};
else
  error(['Invalid regularization method: ',params.regularizer])
end

% interp must be one of:
%    'bilinear', 'nearest', or 'triangle'
% but accept any shortening thereof.
valid = {'bilinear', 'nearest', 'triangle'};
if isempty(params.interp)
  params.interp = 'triangle';
end
ind = find(strncmpi(params.interp,valid,length(params.interp)));
if (length(ind)==1)
  params.interp = valid{ind};
else
  error(['Invalid interpolation method: ',params.interp])
end

% solver must be one of:
%    'backslash', '\', 'symmlq', 'lsqr', or 'normal'
% but accept any shortening thereof.
valid = {'backslash', '\', 'symmlq', 'lsqr', 'normal'};
if isempty(params.solver)
  params.solver = '\';
end
ind = find(strncmpi(params.solver,valid,length(params.solver)));
if (length(ind)==1)
  params.solver = valid{ind};
else
  error(['Invalid solver option: ',params.solver])
end

% extend must be one of:
%    'never', 'warning', 'always'
% but accept any shortening thereof.
valid = {'never', 'warning', 'always'};
if isempty(params.extend)
  params.extend = 'warning';
end
ind = find(strncmpi(params.extend,valid,length(params.extend)));
if (length(ind)==1)
  params.extend = valid{ind};
else
  error(['Invalid extend option: ',params.extend])
end

% tilesize == inf by default
if isempty(params.tilesize)
  params.tilesize = inf;
elseif (length(params.tilesize)>1) || (params.tilesize<3)
  error 'Tilesize must be scalar and > 0.'
end

% overlap == 0.20 by default
if isempty(params.overlap)
  params.overlap = 0.20;
elseif (length(params.overlap)>1) || (params.overlap<0) || (params.overlap>0.5)
  error 'Overlap must be scalar and 0 < overlap < 1.'
end

% ============================================
% subfunction - tiled_gridfit
% ============================================
function zgrid=tiled_gridfit(x,y,z,xnodes,ynodes,params)
% tiled_gridfit: a tiled version of gridfit, continuous across tile boundaries
% usage: [zgrid,xgrid,ygrid]=tiled_gridfit(x,y,z,xnodes,ynodes,params)
%
% Tiled_gridfit is used when the total grid is far too large
% to model using a single call to gridfit. While gridfit may take
% only a second or so to build a 100x100 grid, a 2000x2000 grid
% will probably not run at all due to memory problems.
%
% Tiles in the grid with insufficient data (<4 points) will be
% filled with NaNs. Avoid use of too small tiles, especially
% if your data has holes in it that may encompass an entire tile.
%
% A mask may also be applied, in which case tiled_gridfit will
% subdivide the mask into tiles. Note that any boolean mask
% provided is assumed to be the size of the complete grid.
%
% Tiled_gridfit may not be fast on huge grids, but it should run
% as long as you use a reasonable tilesize. 8-)

% Note that we have already verified all parameters in check_params

% Matrix elements in a square tile
tilesize = params.tilesize;
% Size of overlap in terms of matrix elements. Overlaps
% of purely zero cause problems, so force at least two
% elements to overlap.
overlap = max(2,floor(tilesize*params.overlap));

% reset the tilesize for each particular tile to be inf, so
% we will never see a recursive call to tiled_gridfit
Tparams = params;
Tparams.tilesize = inf;

nx = length(xnodes);
ny = length(ynodes);
zgrid = zeros(ny,nx);

% linear ramp for the bilinear interpolation
rampfun = inline('(t-t(1))/(t(end)-t(1))','t');

% loop over each tile in the grid
h = waitbar(0,'Relax and have a cup of JAVA. Its my treat.');
warncount = 0;
xtind = 1:min(nx,tilesize);
while ~isempty(xtind) && (xtind(1)<=nx)

  xinterp = ones(1,length(xtind));
  if (xtind(1) ~= 1)
    xinterp(1:overlap) = rampfun(xnodes(xtind(1:overlap)));
  end
  if (xtind(end) ~= nx)
    xinterp((end-overlap+1):end) = 1-rampfun(xnodes(xtind((end-overlap+1):end)));
  end

  ytind = 1:min(ny,tilesize);
  while ~isempty(ytind) && (ytind(1)<=ny)
    % update the waitbar
    waitbar((xtind(end)-tilesize)/nx + tilesize*ytind(end)/ny/nx)

    yinterp = ones(length(ytind),1);
    if (ytind(1) ~= 1)
      yinterp(1:overlap) = rampfun(ynodes(ytind(1:overlap)));
    end
    if (ytind(end) ~= ny)
      yinterp((end-overlap+1):end) = 1-rampfun(ynodes(ytind((end-overlap+1):end)));
    end

    % was a mask supplied?
    if ~isempty(params.mask)
      submask = params.mask(ytind,xtind);
      Tparams.mask = submask;
    end

    % extract data that lies in this grid tile
    k = (x>=xnodes(xtind(1))) & (x<=xnodes(xtind(end))) & ...
        (y>=ynodes(ytind(1))) & (y<=ynodes(ytind(end)));
    k = find(k);

    if length(k)<4
      if warncount == 0
        warning('GRIDFIT:tiling','A tile was too underpopulated to model. Filled with NaNs.')
      end
      warncount = warncount + 1;

      % fill this part of the grid with NaNs
      zgrid(ytind,xtind) = NaN;

    else
      % build this tile
      zgtile = gridfit(x(k),y(k),z(k),xnodes(xtind),ynodes(ytind),Tparams);

      % bilinear interpolation (using an outer product)
      interp_coef = yinterp*xinterp;

      % accumulate the tile into the complete grid
      zgrid(ytind,xtind) = zgrid(ytind,xtind) + zgtile.*interp_coef;

    end

    % step to the next tile in y
    if ytind(end)<ny
      ytind = ytind + tilesize - overlap;
      % are we within overlap elements of the edge of the grid?
      if (ytind(end)+max(3,overlap))>=ny
        % extend this tile to the edge
        ytind = ytind(1):ny;
      end
    else
      ytind = ny+1;
    end

  end % while loop over y

  % step to the next tile in x
  if xtind(end)<nx
    xtind = xtind + tilesize - overlap;
    % are we within overlap elements of the edge of the grid?
    if (xtind(end)+max(3,overlap))>=nx
      % extend this tile to the edge
      xtind = xtind(1):nx;
    end
  else
    xtind = nx+1;
  end

end % while loop over x

% close down the waitbar
close(h)

if warncount>0
  warning('GRIDFIT:tiling',[num2str(warncount),' tiles were underpopulated & filled with NaNs'])
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function hmm = hmmrotatepca(hmm,Gamma,XXGXX)

for k = 1:length(hmm.state)

    % Orthogonalize W to the standard PCA subspace
    W = hmm.state(k).W.Mu_W;
    %[W,~] = svd(hmm.state(k).W.Mu_W,'econ');

%     % Enforce a sign convention on the coefficients:
%     % the largest element in each column will have a positive sign.
    [~,maxind] = max(abs(W), [], 1);
    [d1, d2] = size(W);
    colsign = sign(W(maxind + (0:d1:(d2-1)*d1)));
    W = bsxfun(@times, W, colsign);
    hmm.state(k).W.Mu_W = W;

    % recompute covariance of W
    v = hmm.Omega.Gam_rate / hmm.Omega.Gam_shape;
    SW = XXGXX{k} * W / sum(Gamma(:,k));
    M = W'*W+v*eye(size(W,2));
    iS_W = v*eye(size(W,2))+M\W'*SW; S_W = inv(iS_W);
    for n = 1:size(hmm.state(k).W.iS_W,1)
        hmm.state(k).W.iS_W(n,:,:) = iS_W;
        hmm.state(k).W.S_W(n,:,:) = S_W;
    end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function rs = ccf(s1, s2, len)

%tstoolbox/@signal/ccf
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,3);
if (ndim(s1) > 2) | (~isreal(data(s1)))
	help(mfilename)
	return
end

if (nargin < 3)
	if dlens(s1,1) > 256
		len = 128;
	else
		len = nextpow2(dlens(s1,1)/4);
	end
end

c = xcorr(data(s1), data(s2),'coeff');
rs = signal(c, s1);	% special constructor calling syntax for working routines
a = getaxis(s1, 1);
dl = delta(a);
a = setfirst(a, 0);
rs = setaxis(rs, 1, a);
rs = addhistory(rs, 'Cross correlation function');
rs = addcommandlines(rs, 's = ccf(s, s2', len);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ mK ] = CreateConvMtx( vK, numElements, operationMode, convShape )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here

OPERATION_MODE_CONVOLUTION = 1;
OPERATION_MODE_CORRELATION = 2;

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

switch(operationMode)
    case(OPERATION_MODE_CONVOLUTION)
        vK = vK(end:-1:1);
    case(OPERATION_MODE_CORRELATION)
        % vK = vK; %<! Default Code is correlation
end

kernelLength    = length(vK);
mK = zeros([numElements + kernelLength - 1, numElements]);

for ii = 1:numElements + kernelLength - 1
    kernelLastIdx     = min(kernelLength, kernelLength + numElements - ii);
    kernelFirstIdx    = max(kernelLastIdx - ii + 1, 1);

    kernelEffLength = kernelLastIdx - kernelFirstIdx + 1;

    colLastIdx   = min(ii, numElements);
    colFirstIdx  = colLastIdx - kernelEffLength + 1;

    mK(ii, colFirstIdx:colLastIdx)   = vK(kernelFirstIdx:kernelLastIdx);
end

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        % mK = mK;
    case(CONVOLUTION_SHAPE_SAME)
        rowIdxFirst = 1 + floor(kernelLength / 2);
        rowIdxLast  = rowIdxFirst + numElements - 1;
        mK = mK(rowIdxFirst:rowIdxLast, :);
    case(CONVOLUTION_SHAPE_VALID)
        mK = mK(kernelLength:end - kernelLength + 1, :);
end

end

author: RoyiAvital

repo: StackExchangeCodes