Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add a function to save a slice in MDH format #149

Merged
merged 11 commits into from
Sep 12, 2023
55 changes: 55 additions & 0 deletions +sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
classdef unittest_spinw_spec2MDHisto < sw_tests.unit_tests.unittest_super
properties
swModel = [];
tmpdir = '';
testfilename = '';
nsteps = {100};
end

properties (TestParameter)
testpars = struct(...
'test_1_0_0', struct('q0', [0 0 0], 'qdir', [1 0 0], 'proj', [[1 0 0]' [0 1 0]' [0 0 1]'], 'nxs', 'test100mdh.nxs'), ...
'test_1_1_0', struct('q0', [0 0 0], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test110mdh.nxs'), ...
'test_1_1_1', struct('q0', [0 0 0], 'qdir', [1 1 1], 'proj', [[1 1 1]' [1 -1 0]' [1 1 -2]'], 'nxs', 'test111mdh.nxs'), ...
'test_1_1_2', struct('q0', [0 0 2], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test112mdh.nxs'), ...
'test_1_1_2_2', struct('q0', [0 0 2], 'qdir', [1 1 0], 'proj', [[1 -1 0]' [1 1 0]' [0 0 1]'], 'nxs', 'test112_2mdh.nxs'), ...
'test_2_2_2', struct('q0', [2 2 2], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test222mdh.nxs'));
end

methods (TestClassSetup)
function setup_model(testCase)
testCase.swModel = sw_model('triAF', 1);
end
function setup_tempdir(testCase)
testCase.tmpdir = tempdir;
end
end

methods (TestMethodTeardown)
function remove_tmpdir(testCase)
delete(testCase.testfilename);
end
end

methods (Test)
function test_qdirs(testCase, testpars)
q0 = testpars.q0;
qdir = testpars.qdir;
proj = testpars.proj;
testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian');
spec = sw_egrid(spinwave(testCase.swModel, {q0 q0+qdir testCase.nsteps{1}}));
dproj = [norm((qdir-q0))/testCase.nsteps{1}, 1e-6, 1e-6];
testCase.testfilename = fullfile(testCase.tmpdir, testpars.nxs);
sw_spec2MDHisto(spec, proj, dproj, testCase.testfilename);
end

function test_non_ortho(testCase)
q0 = [0 0 2];
qdir = [1 1 0];
spec = sw_egrid(spinwave(testCase.swModel, {q0 q0+qdir testCase.nsteps{1}}));
proj = [qdir(:) [1 0 0]' [0 0 1]'];
dproj = [1e-6, norm((qdir-q0))/testCase.nsteps{1}, 1e-6];
verifyError(testCase,@() sw_spec2MDHisto(spec, proj, dproj, 'tmp/test_blank.nxs'), "read_struct:nonorthogonal")
end
end
end
286 changes: 286 additions & 0 deletions swfiles/sw_spec2MDHisto.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,286 @@
function sw_spec2MDHisto(spectra,proj, dproj, filename)
% saves spectrum to MDHisto
%
% ### Syntax
%
% sw_spec2MDHisto(spectra,proj,dproj,filename)`
%
% ### Description
%
% `sw_spec2MDHisto(spectra,proj,dproj,filename)` saves a
% spectrum that is calculated by sw_egrid
%
% ### Input Arguments
% spectra: a structure calculated by sw_egrid
%
% proj: a 3x3 matrix defining an orthogonal coordinate system
% where each column is a vector defining the orientation
% of the view. One of the vectors must be identical to the Q axis
% defined by the direction of the calculation.
%
% dproj: is a 3 vector that is the bin size in each of the
% directions defined in proj. For the direction of the
% calculation, the value used is internally calcualted from the spectrum.
% It is wise to enter the step size for clarity.
%
% filename: is the name of the nexus file. It will overwrite the existing
% file if one already exists
%
% Example:
% q0 = [0 0 0];
% qdir = [1 0 0];
% nsteps = 100;
% spec = sw_egrid(spinwave(sw_model('triAF', 1), {q0 q0+qdir nsteps}))
% proj = [qdir(:) [0 1 0]' [0 0 1]'];
% dproj = [(qdir(1)-q0(1)/steps, 1e-6, 1e-6];
% sw_spec2MDHisto(spec, proj, dproj, 'testmdh.nxs');
% Note that:
% (1) In the call to `spinwave`, only one q-direction may be specified
% e.g. the HKL specifier must be of the form {q0 q0+qdir nsteps}
% (2) one column in the `proj` matrix must be the q-direction used in
% `spinwave` (e.g. `qdir`).


if nargin==0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We typically add a unit test for this like so

function test_no_input_calls_help(testCase)
% Tests that if call addmatrix with no input, it calls the help
help_function = sw_tests.utilities.mock_function('swhelp');
testCase.swobj.addatom();
testCase.assertEqual(help_function.n_calls, 1);
testCase.assertEqual(help_function.arguments, ...
{{'spinw.addatom'}});
end

swhelp sw_spec2MDHisto
return
end
[unit_cell,Bmat,proj_out,D,dat,proj,name] = read_struct(spectra,proj,dproj);
%check if hdf file exists and delete if it does.
if exist(filename,'file')
delete(filename)
end

h5createnwrite(filename,'/MDHistoWorkspace/coordinate_system',3); % None = 0, QLab = 1, QSample = 2, HKL = 3
h5createnwrite(filename,'/MDHistoWorkspace/visual_normalization',0);
h5writeatt(filename,'/MDHistoWorkspace','NX_class','NXentry');
h5writeatt(filename,'/MDHistoWorkspace','Qconvention','Inelastic');
h5writeatt(filename,'/MDHistoWorkspace','SaveMDVersion', 2);
% write data
rtpth = NXScreategroup(filename,'/MDHistoWorkspace','data','NXdata');
% write D dimensions
Dszs=zeros(1,4);
Dstrcell={};
for idx=1:length(D)
szd=size(D{idx});
Dszs(idx)=szd(2)-1;
Dstrcell{idx} = strcat('D',num2str(idx-1));
Dpth = strcat(rtpth,'/',Dstrcell{idx});
h5createnwritevec(filename,rtpth,Dstrcell{idx},D{idx});
if idx<length(D)
h5writeatt(filename,Dpth,'units','r.l.u');
h5writeatt(filename,Dpth,'frame',mat2str(transpose(proj(:,idx))));
h5writeatt(filename,Dpth,'long_name',mat2str(transpose(proj(:,idx))));
else
h5writeatt(filename,Dpth,'long_name','DeltaE');
h5writeatt(filename,Dpth,'frame','General Frame');
h5writeatt(filename,Dpth,'units','DeltaE');
end
end
%write signal
signal = reshape(dat,Dszs); % change signal array dimensions to match the number of changing dimensions
h5createnwrite(filename,strcat(rtpth,'/signal'),signal);
h5createnwrite(filename,strcat(rtpth,'/errors_squared'),zeros(size(signal)));
h5createnwrite(filename,strcat(rtpth,'/num_events'),zeros(size(signal)));
h5createnwrite(filename,strcat(rtpth,'/mask'),zeros(size(signal),'int8'));
axesstr='';
for idx=1:length(Dstrcell)
if idx<length(Dstrcell)
axesstr=strcat(axesstr,Dstrcell{idx},':');
else
axesstr=strcat(axesstr,Dstrcell{idx});
end
end
h5writeatt(filename,strcat(rtpth,'/signal'),'axes',axesstr);
h5writeatt(filename,strcat(rtpth,'/signal'),'signal',1);

%write experiment
exppth = NXScreategroup(filename,'/MDHistoWorkspace','experiment0','NXgroup');
h5writeatt(filename,exppth,'version', 1)
% write logs
log_pth = NXScreategroup(filename,exppth,'logs','NXgroup');
h5writeatt(filename,log_pth,'version',1)
writeNXlog(filename,log_pth,'W_MATRIX',proj_out,' ')
writeNXlog(filename,log_pth,'RUBW_MATRIX',proj_out,' ')
%write sample
smplpth = NXScreategroup(filename,exppth,'sample','NXsample');
h5writeatt(filename,smplpth,'version',1);
h5writeatt(filename,smplpth,'name',name);
h5writeatt(filename,smplpth,'shape_xml','<type name="userShape"> </type>');
h5createnwritevec(filename,smplpth,'num_oriented_lattice',int32(1))
h5createnwritevec(filename,smplpth,'num_other_samples',int32(0))
h5createnwritevec(filename,smplpth,'geom_height', 0)
h5createnwritevec(filename,smplpth,'geom_id', int32(0))
h5createnwritevec(filename,smplpth,'geom_thickness', 0)
h5createnwritevec(filename,smplpth,'geom_width' ,0)
%write material
mtl_pth = NXScreategroup(filename,smplpth,'material','NXdata');
h5writeatt(filename,mtl_pth,'formulaStyle','empty')
h5writeatt(filename,mtl_pth,'name',' ')
h5writeatt(filename,mtl_pth,'version',int32(2))
h5createnwritevec(filename,mtl_pth,'packing_fraction',1)
h5createnwritevec(filename,mtl_pth,'number_density',0)
h5createnwritevec(filename,mtl_pth,'pressure',0)
h5createnwritevec(filename,mtl_pth,'temperature',0)
%write crystal lattice
OL_pth = NXScreategroup(filename,smplpth,'oriented_lattice','NXcrystal');
%write lattice parameters
u_parm_names={'a','b','c','alpha','beta','gamma'};
for idx=1:length(u_parm_names)
parm_name = strcat('unit_cell_',u_parm_names{idx});
h5createnwritevec(filename,OL_pth,parm_name,unit_cell(idx))
end
%write orientation matrix
om_path =strcat(OL_pth,'/orientation_matrix');
h5createnwrite(filename,om_path,Bmat);
%write instrument
instr_pth = NXScreategroup(filename,exppth,'instrument','NXinstrument');
h5writeatt(filename,instr_pth,'version',int32(1))
h5createnwritevec(filename,instr_pth,'name','SEQUOIA');
end

function h5createnwrite(filename,path,val)
dtyp = class(val);
h5create(filename,path,size(val),'Datatype',dtyp);
h5write(filename,path,val);
end

function h5createnwritevec(filename,group,ds,val)
fh = H5F.open(filename,'H5F_ACC_RDWR','H5P_DEFAULT');
sph = H5S.create_simple(1,length(val),length(val));
gsh = H5G.open(fh,group);
memtype = 'H5ML_DEFAULT';
switch class(val)
case 'int32'
H5type = 'H5T_NATIVE_INT';

case 'int8'
H5type = 'H5T_NATIVE_CHAR';
case 'char'
dims= size(val);
sph = H5S.create_simple(1,fliplr(dims(1)),[]);
H5type = H5T.copy('H5T_FORTRAN_S1');
H5T.set_size(H5type,(dims(2)+1))
memtype = H5T.copy ('H5T_C_S1');
H5T.set_size(memtype,dims(2));
otherwise
H5type = 'H5T_NATIVE_DOUBLE';
end
dsh = H5D.create(gsh,ds,H5type,sph,'H5P_DEFAULT');
H5D.write(dsh,memtype,'H5S_ALL','H5S_ALL','H5P_DEFAULT',val)
H5S.close(sph)
H5D.close(dsh)
H5G.close(gsh)
H5F.close(fh)
end

function pthout = NXScreategroup(filename,pth,group,NX_class)
% ### Syntax

% pthout = NXScreategroup(filename,pth,group,NX_class)

% ### Description
%
% create a group with attributes of a nexus class
%
% ### Input Arguments
% filename, name of hdf5 file
% pth, path to where the group should be created
% group the group name
% NX_class a string containing a valid NX_class definition
%
% ### Output Arguments
% returns a path to the group

fh = H5F.open(filename,'H5F_ACC_RDWR','H5P_DEFAULT');
pthh = H5G.open(fh,pth);
gh = H5G.create(pthh,group,'H5P_DEFAULT','H5P_DEFAULT','H5P_DEFAULT');
pthout =strcat([pth,'/',group]);
H5G.close(gh)
H5G.close(pthh)
H5F.close(fh)
h5writeatt(filename,pthout,'NX_class',NX_class)
end

function writeNXlog(filename,log_pth,log_nm,value,units)
% ### Description
%
% create and write a Nexus log
%
% ### Input Arguments
%filename, name of hdf5 file
%log_pth, path to where the log should be created
%log_nm the name of the log
% value the vlaue of the log
% the units if any
pth = NXScreategroup(filename,log_pth,log_nm,'NXlog' );
h5createnwritevec(filename,pth,'value',value)
h5writeatt(filename,pth,'units',units)

end
function [latt_parms,Bmat,proj_out,D,signal,proj,name] = read_struct(dstruct,proj,dproj)

% ### Input Arguments
% dstruct: is the spinw structure
% proj: is the viewing projection matrix each column is a different axis.
% One axis must be parallel to the drection of propogation in the
% dstruct
% drproj: is the size of the bin in each direction (ignored for the direction
% of the cut).
%
% ### Output Arguments
% latt_parms: the lattice parameters from the spinw structure
% Bmat: a matrix for converting from inverse angstroms to rlu
% proj_out:
% D: a cell array of the number of steps in each direction
% signal: the signal array from the spinw spec strcuture
% name : the chemical formula from the spinW file
check_ortho(proj)
fnames=fieldnames(dstruct);
objnum = find(strcmp(fnames,'obj'));
swobj = dstruct.(subsref(fnames,substruct('{}',{objnum})));
latt_parms = abc(swobj);
M = basisvector(swobj);
Bmat = inv(M);
name =formula(swobj).chemform;
%proj_out = proj(:);
hkls = dstruct.hkl;
hkls_sz = size(hkls);
dir_vec = hkls(:,hkls_sz(2))-hkls(:,1);
%qout = hkls'/dir_vec';
D={};
for idx=1:3
procjv = proj(:,idx)/norm(proj(:,idx));

if abs(norm(cross(dir_vec,procjv)))> 1e-6
dtmp = dot(hkls(:,2),procjv);
D{idx} = dtmp+dproj(idx)/2.*[-1 1];
else
hkl_proj = hkls'/dir_vec';
dhkl = hkl_proj(2)-hkl_proj(1); % get the spacing along the q axis
%hkl_proj = dot(hkls(:,1),procjv);
D{idx} = zeros([1,length(hkl_proj)+1]);
D{idx}(1:length(hkl_proj)) = hkl_proj-dhkl/2;
D{idx}(length(D{idx})) = hkl_proj(length(hkl_proj))+dhkl/2;% change to bin boundaries
proj(:,idx) = dir_vec; %set varying projection vector to spectra object
end
end
D{4}=dstruct.Evect;

signal = dstruct.swConv;
proj_out = proj(:);
end

function check_ortho(mat)
% check if the three column vectors in proj are orthogonal to each other
sm = size(mat);
for idx = 1:sm(2)
for idx2 = 1:sm(2)
if idx~=idx2
if norm(dot(mat(:,idx),mat(:,idx2))) >1e-6
error("read_struct:nonorthogonal","the 3 vectors in proj must form an orthogonal basis set")
end
end
end
end
end