function [Bss,Pnl]=MakBss(X,Prm)
%**************************************************************************
%Illustration of non-stationary extreme value modelling in MATLAB accompanying 
%Threshold modeling of nonstationary extremes
%Paul Northrop, Philip Jonathan and David Randell
%in 
%Extreme Value Modeling and Risk Analysis: Methods and Applications
%Dipak Dey and Jun Yan (Editors)
%2015
%**************************************************************************

%**************************************************************************
%Generate spline basis
%
% X        n x 1 data
% Xl       1 x 1 left limit of X data
% Xr       1 x 1 right limit of X
% ndX      1 x 1 number of segments wanted
% SplDgr   1 x 1 degree of spline
% PnlOrd   1 x 1 difference penalty order
% IsPrd    1 x 1 optional periodic b-spline
%
% Bss      n x q basis matrix
% Pnl      q x q penatly matrix   Pnl=(D' * D)
%**************************************************************************

Xl=Prm.BssDmn(1);  %left limit of data
Xr=Prm.BssDmn(2);  %left limit of data
ndX=Prm.BssnSgm;
SplDgr=Prm.BssSplDgr;
PnlOrd=Prm.BssPnlOrd;
IsPrd=Prm.BssPnlOrd;

%Standardise for stability
X=(X-Xl)/(Xr-Xl);
Xl=0;
Xr=1;

%Estimate truncated power functions
dX = (Xr - Xl) / ndX ;
knots = (Xl - SplDgr * dX :dX: Xr + SplDgr * dX)';
Pw=NaN(length(X),length(knots));
for i=1:length(X)
	for j=1:length(knots)
		Pw(i,j)=tpower(X(i),knots(j),SplDgr); %grid of values of the truncated power function for all possible (X,knots)
	end
end

%Convert truncated power functions into b-spline basis set
n = length(knots);
if ndX>SplDgr
    tD = diff(eye(n),SplDgr + 1) ./ (gamma(SplDgr + 1) * dX ^ SplDgr);
    Bss = (-1) ^ (SplDgr + 1) * Pw * tD';
    m = length(Bss(1,:));
else
    error('not enough knots for spline degree')    
end

%Make periodic b-spline if required
if IsPrd==1; %periodic version
	Bss(:,1:SplDgr) = Bss(:,1:SplDgr)+Bss(:,m-(SplDgr-1):m);
	Bss=Bss(:,1:m-SplDgr);
end;

%Compute Penalty Matrix
if IsPrd==0; %non-periodic
	if PnlOrd==0;
		D = eye(m);
    else
        if m>PnlOrd
            D = diff(eye(m), PnlOrd);
        else
            error('not enough knots for penalty degree')
        end
	end;
	Pnl = D' * D;
else; %periodic
  %20130115: this was previously using PnlOrd+1 instead of SplDgr, which works when
  %SplDgr=3 and PnlOrd=2 (our default) but was not correct in general
  if (m-SplDgr)>PnlOrd
        d1=diff(eye(m-SplDgr),PnlOrd);
        d2=d1(1,:);
        d3=NaN(m-SplDgr,m-SplDgr);
        for i=1:m-SplDgr
            d3(i,:)=circshift(d2,[0,i-1]);
        end        
        Pnl = d3' * d3;
    else
        error('not enough knots for penalty degree')
    end
end;

return;

%*************************************************************
function Pnl = tpower(x, knots, SplDgr)

% TRUNCATED POWER FUNCTION - USING TO CONSTRUCT Bss-SPLINE BASIS MATRIX
%See "Splines, Knots & Penalties" (Eilers & Marx 2004), and "The Craft of Smoothing"
%presentation - slide 23-24

[X,T] = meshgrid(x,knots);
Pnl=((X - T) .^ SplDgr .* (X > T))';

%Normal completion
return;
%*************************************************************