function [med,w,norms] = SpMed(t,x,dtyp)
%function [med,w,norms] = SpMed(t,x,dtyp)
%
%Functional Spatial Median
%
%INPUT: It is assumed that all curves are observed on the same time grid.
% The time grid can be irregular (doesn't need to be equispaced).
%
% Input variables:
%   t: Time grid (1 x m vector). Must be SORTED (in increasing order).
%   x: Discretized curves (n x m matrix). Each row represents a different
%      curve. Missing values NOT allowed. If your data set has missing
%      values, use either interpolation or smoothing to fill in the gaps.
%   dtyp: Data type (optional; character variable). Enter 's' if the data 
%         is "reasonably" smooth, or 'n' if the data is very noisy (in that
%         case a bias correction is necessary and the error variance is 
%         estimated nonparametrically using the GSJ estimator; this will
%         slow things down, so avoid if possible. It may also produce an 
%         inner-product matrix that is not nonnegative definite, in which
%         case the uncorrected matrix will be used). Default value: 's'.
%
% Output variables:
%   med: Spatial median (1 x m vector).
%   w: Weights (n x 1 vector). The spatial median is w'*x.
%   norms: Distances between each curve and the median (n x 1 vector).
%          Note that the median minimizes sum(norms) by definition, and
%          that w = norms.^(-1)/sum(norms.^(-1)).
%
% External functions called: gsj.m
%

%% Input check
if nargin<2
    error('Not enough input variables')
elseif nargin==2
    dtyp = 's';
end
[n,m] = size(x);
if size(t,1)>1
    t = t';
end
if size(t,2)~=m
    error('Dimensions of T and X not compatible')
end
if ~(strcmpi(dtyp,'n')||strcmpi(dtyp,'s'))
    error('Wrong input value for DTYP')
end

%% Inner-product matrix
A = 0.5*(x(:,1:m-1)*diag(t(2:m)-t(1:m-1))*x(:,1:m-1)' + ...
    x(:,2:m)*diag(t(2:m)-t(1:m-1))*x(:,2:m)');
if strcmpi(dtyp,'n')
    se2 = gsj(repmat(t',1,n),x');
    A1 = A - diag(se2)*(t(m)-t(1));
    if all(eig(A1)>-1e-6)
        A = A1;
    else
        disp('Corrected inner-product matrix is not nonnegative definite')
        disp('Using uncorrected matrix instead')
    end
end

%% Iterative minimization from sample mean
w = ones(n,1)/n;
norms = sqrt(diag(A) + w'*A*w - 2*A*w);
f = sum(norms);
err = 1;
iter = 0;
while err>1e-5 && iter<50
    iter = iter+1;
    f0 = f;
    if any(norms<eps)
        i0 = find(norms<eps);
        w = zeros(n,1);
        w(i0) = 1/length(i0);
    else
        w = 1./norms;
        w = w/sum(w);
    end
    norms = sqrt(diag(A) + w'*A*w - 2*A*w);
    f = sum(norms);
    err = abs(f/f0-1);
end

%% Output
med = w'*x;