function [ptrns,bmus,rmsErr,pfrq,sMap,sD] = som_make_njs(D,somRows,somCols,init,wt)
%SOM_MAKE_NJS Make SOMs according to particular specifications and output
%specific SOM needs and statistics
%   Function based on Nat Johnson's SOM code.
%
%   Inputs:
%       D[lat*lon,time] = input data with one space and one time dimensioni
%	varargin => includes options for only computing the data structure and not SOMs
%
%   Outputs:
%       ptrn[k,lat*lon] = k SOM patterns
%       bmus[time] = best matching units (pattern)
%       pfrq[k] = pattern frequency
%
%   Nathan Steiger -- Apr 2017

% ADD PATHS TO SOM FUNCTIONS
addpath('/home/nsteiger/MATLAB/SOM-Toolbox')
addpath('/home/nsteiger/MATLAB/SOM-Toolbox/som')
addpath('/home/nsteiger/MATLAB/SOM-Toolbox/dijkstra')

% INITIAL PARAMETER CHOICES

% See Liu et al.,2006, "Performance evaluation of the self-organizing map for feature
% extraction," JGR, for some analysis of the sensitivity wrt some of these
% parameter choices.  In general, the biggest choice is the number of rows
% and columns (i.e., number of patterns) in the SOM.)

%somRows = 2; % number of rows in SOM
%somCols = 3; % number of columns in SOM

%init = 1; disp('Linear initialization...')
%init = 0; disp('Random initialization...')
% 1 = linear initialization, 0 = random initialization.  (Note:
% for input matrices with many variables (i.e., grid points in this case),
% linear initialization will be the most time consuming part of the
% program, but the initial SOM ordering will be better than in random
% ordering.  If you choose random initialization, then you probably should
% repeat the analysis several times to make sure that results are not too
% sensitive to the initial ordering.)

lattice = 'rect'; % 'rect' for rectangular or 'hexa' for hexagonal
shape = 'sheet'; % 'sheet', 'cyl', or 'toroid'

neighborhood_fct = 'ep'; % choice of neighborhood function: 'gaussian',
% 'cutgauss', 'bubble', or 'ep' (Epanechikov function)

rad_ini = somRows*somCols; % initial neighborhood radius used in the training of the
% SOM (around the size of the minimum or maximum dimension should be
% fine for most maps)
rad_fin = 1; % final neighborhood radius at the end of SOM training (1
% usually provides a good fit and good final ordering, 0 provides the
% maximum fit, but possibly at the expense of ordering and with a
% greater risk of overfitting - depending on the purpose, that might
% not be a major drawback)
trainlen_rough = 50; % Number of training iterations during the rough
% initial phase to achieve a broad ordering
trainlen_finetune = 300; % Number of training iterations during the
% finetuning phase

%==================
% Pre-process data
%==================

[~,yrs]=size(D);

% Standardize the data using SOM function, along appropriate dimension
Dn=som_normalize(D','var');
if isempty(wt)
   Dnw=Dn;
else
   % Weight standardized data according to grid cell area
   %W=cosd(repmat(lat,[1 length(lon)]));
   %Dnw=bsxfun(@times,Dn,W(:)');
   Dnw=bsxfun(@times,Dn,wt(:)');
end

% Add a mask if desired
%if ~isempty(mask)
%   Dnw=bsxfun(@times,Dnw,mask(:)');
%end

%===================
% Compute the SOMs
%===================

% Establish the data struct, as described in the SOM toolbox manual
sD = som_data_struct(Dnw);


% Initialize the map
if init=='l'
    disp('Initializing the map using EOFs...')
    sMap = som_lininit(sD, 'msize', [somCols somRows], lattice, shape);
elseif init=='r'
    disp('Initializing the map using random initialization...')
    sMap = som_randinit(sD, 'msize', [somCols somRows], lattice, shape);
end

% Train the map with the batch version.  First do a broad ordering
% with relatively large neighborhood radius and a small number of
% training interations.
disp('Rough training phase...')
% Initializing things with recommended parameters into 'sTrain' plus some amendations of my own
%sTrain = som_train_struct('train', sD);
sTrain = som_train_struct('rough',sD);
sTrain.neigh = neighborhood_fct; % Set the neighborhood function
sTrain.radius_ini = rad_ini; % Set the initial neighborhood radius (default = 3)
sTrain.radius_fin = rad_fin; % Set the final neighborhood radius (default = 1)
sTrain.trainlen = trainlen_rough;

% Perform the initial rough training
sMap = som_batchtrain(sMap, sD, sTrain);

% Finetune in the second phase of training
disp('Finetuning traning phase...')
sTrain.radius_ini = rad_fin + 2; % this will almost always = 3
sTrain.trainlen = trainlen_finetune;
[sMap, sTopol] = som_batchtrain(sMap,sD, sTrain);

disp('Computing errors of maps and best-matching-units...')
% Get the average quantization (rms) error and topographic error
% as measures of quality
%[aqe, te] = som_quality(sMap,sD); % aqe is the average quantization
% error, and te is the topographic error

% Get the best-matching unit (BMU) time series.  The SOM patterns
% are numbered in a row-wise manner under this configuration, and
% so a BMU of 1 at time 1 corresponds to the occurrence of the SOM
% pattern at the top left corner of the SOM at time 1.   Pattern 2
% would be the second pattern in row 1, etc., and the last pattern
% would be on the bottom right corner of the SOM.  Also, we get the
% quantization (rms) error at each time, which measures the
% disparity between the input data field and the best-matching SOM
% pattern at each time.
[bmus,rmsErr] = som_bmus(sMap,sD); % bmus is the SOM pattern time
% series, and qerrs is the quantization error time series.

% Get the SOM pattern matrix from the codebook
ptrns0 = sMap.codebook;

if isempty(wt)
   ptrns=ptrns0;
else
   % Remove latitude weighting
   ptrns=bsxfun(@rdivide,ptrns0,wt(:)');
end


% Calculate the frequency of occurrence of each pattern
pfrqT = histcounts(bmus,somCols*somRows);
pfrq=pfrqT/yrs;


end

