function varargout = cms_read_sol(filename,varargin)
% sol = cms_read_sol(filename,'name1',value1,'name2',value2,...)
%
% Description:
%   Reads a CMS solution file and creates
%   a structure variable containing the solution
%   datasets values and times
%
% Input:
%   filename - input file name including full path
%   grd - Grid structure containing the following fields
%     nx - Grid size in x direction
%     ny - Grid size in y direction
%   varargin - dataset names
%
% Output:
%   sol - structure variable containing dataset values and times
%   varargout - variable length datasets corresponding to varargin
%
% Usage:
%   filename = 'test_sol.h5';
%   sol = cms_read_sol(filename,grd);
%   wse = cms_read_sol(filename,grd,'Water_Elevation');
%   [wse,uv] = cms_read_sol(filename,'Water_Elevation','Current_Velocity');
%   sol = cms_read_sol(filename,'grid',grd,'datasetnames',dnames,...
%     'points',pts,'cells',ids,'polygon',poly,'steps',nstep,'times',t)
%
% Author: Alex Sanchez, USACE-ERDC-CHL

grd = [];
pts = [];
% poly = [];
steps = [];
cells = [];
t = [];
dnames = [];

if nargin==0
    demo_cms_read_sol
    return
end

if mod(length(varargin),2)~=0
    error('Invalid number of input options')
end

%Get options
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'grid'
            grd = varargin{i+1};
        case {'datasets','dataset','datasetnames','datasetname',...
                'names','name'}
            dnames = varargin{i+1};
            if ~iscell(dnames)
                dnames = {dnames};
            end
        case {'point','points'}
            pts = varargin{i+1};
            %         case {'poly','polygon'}
            %             poly = varargin{i+1};
        case {'times','time'}
            t = varargin{i+1};
        case {'step','steps'}
            steps = varargin{i+1};
        otherwise
            warning(['Unidentified Property: ',varargin{i}])
    end
end

if ~isempty(pts)
    if isempty(grd)
        error('Cannot extract points without grid information.')
    end
    [x,y] = meshgrid(grd.x,grd.y);
    x = x'; y = y';
    x = x(:); y = y(:);
    [ii,jj] = meshgrid(1:grd.nx,1:grd.ny);
    ii = ii'; jj = jj';
    ii = ii(:); jj = jj(:);
    %Find nearhest cells to points
    npts = size(pts,1);
    cells = zeros(npts,1);
    icol = cells; jrow = cells;
    for k=1:npts
        dist = sqrt((x-pts(k,1)).^2+(y-pts(k,2)).^2);
        [~,cells(k)] = min(dist);
        icol(k) = ii(cells(k));
        jrow(k) = jj(cells(k));
    end
%     figure
%     contourf(grd.x,grd.y,grd.depth,20)
%     hold on
% %     plot(x(cells),y(cells),'ok','markersize',4)
%     plot(grd.x(icol),grd.y(jrow),'ok','markersize',4)
%     hold off
end

info = hdf5info(filename);
s = info.GroupHierarchy.Groups.Groups;
ns = length(s);

nd = length(dnames);

%Determine what datasets to read
if nd>0 %try selecting only few datasets
    ilist = zeros(nd,1);
    j = 0;
    for k=1:nd %Search for each dataset
        for i=1:ns
            field = s(i).Name;
            ind = strfind(field,'/');
            field = field(ind(end)+1:end);
            if strcmpi(field,dnames{k})
                j = j + 1;
                ilist(j) = i;
                break
            end
        end
    end
    ilist(ilist==0) = []; %Remove datasets not found
else
    ilist = 1:ns;
end
nout = length(ilist);

if nd>nargout && nargout>1
    error('Cannot have more output variables than input datasets')
end

%Initialize output variable
varargout = cell(nout,1);

for k=1:nout
    i = ilist(k);
    %Read dataset times
    varargout{k}.times = cms_read_sol_times(filename,s(i).Name,steps);
    %If times specified, find time steps
    if ~isempty(t)
        steps = zeros(1,length(t));
        for j=1:length(t)
            ii = find(t(j)==varargout{k}.times);
            if ~isempty(ii)
                steps(j) = ii;
            end
        end
        steps(steps==0) = [];
        varargout{k}.times = varargout{k}.times(steps);
    end
    %Read dataset values
    varargout{k}.values = cms_read_sol_values(filename,...
        s(i).Name,steps,cells);
    if ~isempty(steps)
        varargout{k}.steps = steps;
    end
    if ~isempty(cells)
        varargout{k}.cells = cells;
        varargout{k}.icol = icol;
        varargout{k}.jrow = jrow;
    end
    if ~isempty(pts)
        varargout{k}.x = x(cells);
        varargout{k}.y = y(cells);
    end
    %Reshape if grid data is specified and not point output
    if ~isempty(grd) && isempty(cells)
        nt = length(varargout{k}.times);
        varargout{k}.values = reshape(varargout{k}.values,[grd.nx,grd.ny,nt]);
        varargout{k}.values = permute(varargout{k}.values,[2,1,3]);
    end
end

%Combine all into single variable
if nargout==1 && nout>1
    sol = [];
    for k=1:nout
        i = ilist(k);
        field = s(i).Name;
        ind = strfind(field,'/');
        field = field(ind(end)+1:end);
        sol.(field) = varargout{k};
    end
    varargout{1} = sol;
end

%Simulation label
simlabel = info.GroupHierarchy.Groups.Name(2:end);
for k=1:nargout
    varargout{k}.label = simlabel;
end

return
end

%------------------------------------------------------------------
function times = cms_read_sol_times(fname,dpath,steps)
% CMS_READ_SOL_TIMES  Reads a time step from a CMS solution file
%
% Input:
%   fname - XMDF file name
%   dpath - Dataset path 
%   steps - Time steps
%
% Output:
%    times - Dataset time [hrs]
%
% Author: Alex Sanchez, USACE-CHL


fid = H5F.open(fname); %Open file
dset_id = H5D.open(fid,[dpath,'/Times']); %Open times dataset
if nargin==2
    steps = [];
end
if isempty(steps) %read all times
    times = H5D.read(dset_id)';
else
    file_space_id = H5D.get_space(dset_id);
    [~,dims] = H5S.get_simple_extent_dims(file_space_id);
    if steps(end)>dims(1)
        warning('Specified time step larger than steps in file')
        steps(end) = dims(1);
    end
    stepmax = max(steps);
    stepmin = min(steps);
    block = stepmax - stepmin + 1;
    mem_space_id = H5S.create_simple(1,block,[]);
    offset = stepmin - 1;
    H5S.select_hyperslab(file_space_id,'H5S_SELECT_SET',offset,[],[],block);
    times = H5D.read(dset_id,'H5ML_DEFAULT',mem_space_id,file_space_id,'H5P_DEFAULT')';
end
H5D.close(dset_id); %close dataset
H5F.close(fid); %close file

return
end

%----------------------------------------------------------------
function values = cms_read_sol_values(fname,dpath,steps,cells)
% CMS_READ_SOL_VALUES  Reads a time step from a CMS solution file
%
% Input:
%   fname - XMDF file name
%   dpath - Dataset path
%   steps - Time steps
%   cells - Cell ID's
%
% Output:
%   values - Dataset values
%
% Author: Alex Sanchez, USACE-CHL

%Open file
fid = H5F.open(fname);
%Values
dset_id = H5D.open(fid,[dpath,'/Values']);

if nargin==2
    steps = [];
end
if nargin==3
    cells = [];
end

file_space_id = H5D.get_space(dset_id);
[~,dims] = H5S.get_simple_extent_dims(file_space_id);
% nt = dims(1);
% ncells = dims(2);
% ndim = dims(3);
rank = length(dims);
offset = zeros(1,rank);
block = dims;
if ~isempty(steps)
    stepmin = min(steps);
    stepmax = max(steps);
    offset(1) = stepmin-1;
    block(1) = stepmax-stepmin+1; %for each time step
end
if ~isempty(cells)
    cellmin = min(cells);
    cellmax = max(cells);
    offset(2) = cellmin-1;
    block(2) = cellmax-cellmin+1;
end
if rank==2 %Time-varying scalar
    mem_space_id = H5S.create_simple(rank,block,[]);
    H5S.select_hyperslab(file_space_id,'H5S_SELECT_SET',offset,[],[],block);
    values = H5D.read(dset_id,'H5ML_DEFAULT',mem_space_id,file_space_id,'H5P_DEFAULT');
elseif rank==3 %time-varying vector
    block(3) = 1; %for each component
    mem_space_id = H5S.create_simple(rank,block,[]);
    H5S.select_hyperslab(file_space_id,'H5S_SELECT_SET',offset,[],[],block);
    values = H5D.read(dset_id,'H5ML_DEFAULT',mem_space_id,file_space_id,'H5P_DEFAULT');
    offset(3) = 1;
    H5S.select_hyperslab(file_space_id,'H5S_SELECT_SET',offset,[],[],block)
    vec2 = H5D.read(dset_id,'H5ML_DEFAULT',mem_space_id,file_space_id,'H5P_DEFAULT');
    values = values + 1i*vec2;
    values = permute(values,[2,3,1]);
end
if ~isempty(cells)
    ind = cells-cellmin+1;
    values = values(ind,:);
end
if ~isempty(steps)
    ind = steps-stepmin+1;
    values = values(:,ind);
end
H5D.close(dset_id); %close dataset
H5F.close(fid); %close file

values(real(values)==-999) = NaN;
values = double(values);

return
end

%------------------------------------------------------------------
function demo_cms_read_sol

fgrid = 'D:\CMS\test_cases\hypothetical\hdf5\Flow_small_grid.h5';
grd = cms_read_grid(fgrid);
fname = 'D:\CMS\test_cases\hypothetical\hdf5\Flow_small_vel.h5';
%Option 1: Read whole solution file with no grid information
sol = cms_read_sol(fname);
disp(sol)
%Option 2: Read whole solution file with grid information
sol2 = cms_read_sol(fname,'grid',grd);
disp(sol2)
%Option 3: Read single dataset with grid information
sol3 = cms_read_sol(fname,'grid',grd,...
    'datasets','Current_Velocity');
disp(sol3)
%Option 4: Read multiple datasets with grid information
sol4 = cms_read_sol(fname,'grid',grd,...
    'datasets',{'Current_Velocity','Current_Magnitude'});
disp(sol4)
%Option 5: Read multiple datasets with grid information and separate output
[vel,mag] = cms_read_sol(fname,'grid',grd,...
    'datasets',{'Current_Velocity','Current_Magnitude'});
disp(vel)
disp(mag)
%Option 5: Read multiple datasets with grid information and separate output
%at a point
pt = [1500,1500];
solpt = cms_read_sol(fname,'grid',grd,...
    'datasets','Current_Magnitude','point',pt);
disp(solpt)
%Option 6: Read snap shot of single dataset
sol6 = cms_read_sol(fname,'grid',grd,...
    'datasets','Current_Velocity','times',15);
disp(sol6)
%Option 7: Read snap shot of single dataset
sol7 = cms_read_sol(fname,'grid',grd,...
    'datasets','Current_Velocity','times',15:18);
disp(sol7)

return
end

