% File Inputs for the Fortran inversion code
%
% Fortran inversion code requires three types of input binary files:
%
% (1) Dense vectors (measurement vector, measurement error-covariance stored as a diagonal, prior flux
%
% (2) A dense matrix (parameter and coordinate files)
%
% (3) A sparse matrix (Jacobians in Sparse-Coordinate format)
%
% We show how to write these files in Binary formats in Matlab and R
%
% To write a dense matrix in binary format in Matlab and R we would create some arbitrary data and then write in a binary form
%
% ALWAYS FOR PARAMETER FILE FIRST COLUMN is VARIANCE, SECOND COLUMN IS SPATIAL CORRELATION LENGTH, AND THIRD COLUMN IN TEMPORAL CORRELATION LENGTH
%
% ALWAYS FOR COORDINATE FILE FIRST COLUMN IS LATITUDE AND SECOND COLUMN IS LONGITUDE
%
% create random matrix with 10 rows and 5 columns and write in binary format
%
% numbers of rows in a file would always be a 4 byte integer
%
% numbers of columns in a file would always be a 4 byte integer
% ALL THE TEST FILES WOULD BE WRITTEN IN THE PWD (PATH OF WORKING
% DIRECTORY)
clc
clear all
% MATLAB [Write Dense Matrix]
z=rand(10,5);
fid = fopen( 'z.bin' , 'w' );
fwrite(fid, size(z,1), 'int32' );
fwrite(fid, size(z,2), 'int32' );
fwrite(fid, z, 'double');
fclose(fid);
% Read the above written file
%% MATLAB [Read Dense Matrix]
fid = fopen( 'z.bin' , 'r' );
rows=fread(fid, 1, 'int32' );
cols=fread(fid, 1, 'int32' );
z=fread(fid,[rows,cols], 'double' );
fclose(fid);
%% MATLAB [Write Dense Vector]
% create random matrix with 10 rows and 1 columns and write in binary format
% numbers of rows in a file would always be a 4 byte integer. The only difference between the two is that we never write columns in a vector
z=rand(10,1);
fid = fopen( 'z.bin' , 'w' );
fwrite(fid, size(z,1), 'int32' );
fwrite(fid, z, 'double' );
fclose(fid);
%% Read the above written file
MATLAB [Read Dense Vector]
fid = fopen( 'z.bin' , 'r' );
rows=fread(fid, 1, 'int32' );
z=fread(fid,rows, 'double' );
fclose(fid);
%% For creating H matrices in sparse format in Matlab
%MATLAB [Create SPARSE H]
% Lets say your original Jabobian is called as H
% We use coordinate form for sparse matrices
% These matrices are defined by there indices and non-zero entries
% however there are two other things required one knowledge on of
% non-zero entries in the full matrix and the entry of zero in last column
% if this has no other entry. This is required as we can convert back these
% matrices to full matrices if there is a requirement to do so
% Due to this reason a matrix always at least would have one non-zero entry i.e.
% of last row and column even everything is zero in a matrix.
% For writing H in strips we account for indices for a strip (matrix) not the
% full matrix from which the matrix is derived
precision= 'double' ;
k=1;
H=rand(10,5);
H=sparse(H);
% create some zeros in this matrix. In reality H would be extremely sparse
H(1,5)=0;H(10,5)=0;H(10,4)=0;H(6,4)=0;
if (nnz(H)~=0)
[i,j,val] = find(H); % find non-zero indexes and values
data_dump = [i,j,val];
clear i j val
% this if basically see that whether the last row and
% last column is non zero if it is then this if statement is skipped
% otherwise we add the row and column to data_dump with val=0
% this is important as it allows us to create a dense matrix from a
% sparse matrix if necessary
if ~(data_dump(end,1)==size(H,1)&& data_dump(end,2)==size(H,2))
data_dump = [data_dump;size(H,1),size(H,2),0];
end
fid = fopen(['H' ,num2str(k), '.bin' ], 'w' );
fwrite(fid, size(data_dump,1), 'int32' );
fwrite(fid, data_dump(:,1:2), 'int32' );
fwrite(fid, data_dump(:,3), precision);
fclose(fid);
clear data_dump
elseif (nnz(H)==0) % if all entries are zeros in a H
fid = fopen(['H' ,num2str(k), '.bin' ], 'w' );
fwrite(fid, 1, 'int32' );
fwrite(fid, [size(H,1), size(H,2)],'int32');
fwrite(fid, 0, precision);
end
%% Writing these H matrices in sparse format in a for loop as we store them on harddrive
%MATLAB [Create SPARSE H]
locations=5;
time_periods=5;
measurements=2;
precision= 'double' ;
H_big=rand(measurements,locations*time_periods);
for k=1:time_periods % no of H's strip
H=sparse(H_big(:,((k*locations)+1)-locations:k*locations)); % this is actual big H files
if (nnz(H)~=0)
[i,j,val] = find(H);
data_dump = [i,j,val];
clear i j val
if ~(data_dump(end,1)==size(H,1)&& data_dump(end,2)==size(H,2))
data_dump = [data_dump;size(H,1),size(H,2),0];
end
fid = fopen(['H' num2str(k) '.bin' ], 'w' ); % save in path of working directory pwd
fwrite(fid, size(data_dump,1), 'int32' );
fwrite(fid, data_dump(:,1:2), 'int32' );
fwrite(fid, data_dump(:,3), precision);
fclose(fid);
clear data_dump
elseif (nnz(H)==0)
fid = fopen(['H' num2str(k) '.bin' ], 'w' );
fwrite(fid, 1, 'int32' );
fwrite(fid, [size(H,1), size(H,2)],'int32');
fwrite(fid, 0, precision);
end
display(k)
end