function [A_intersection_B,A_minus_B,B_minus_A,A_union_B] = ioami_jaccard(A,B,type)
% note this gives A-B,B-A and [A intersection B (or IOAMI)] all normalized by A Union B
AB=[A(:) B(:)];
AB(any(isnan(AB), 2), :) = [];
A_intersection_B=sum(min(AB,[],2)); % total shared intensity vector i.e., A intersection B
A_plus_B=sum(sum(AB)); % A+B = total intensity of the signal sum(A union B) but intersected area counted twice as part of A as well as B
A_union_B=A_plus_B-A_intersection_B; % sum(A union B) by removing intersected area that has been counted as part of A as well as B
%% Non Normalized
A_minus_B=A_union_B-nansum(AB(:,2));% This sum(A-B)
B_minus_A=A_union_B-nansum(AB(:,1));% % This sum(B-A)
if strcmpi(type,'normalized')
    %% Normalization by A union B
    A_minus_B=A_minus_B/A_union_B;% This sum(A-B)/sum(A union B) in terms of sets
    B_minus_A=B_minus_A/A_union_B;% This sum(B-A)/sum(A union B)
    %IOAMI AS DEFINED IN PAPER%
    A_intersection_B=A_intersection_B/A_union_B; %This is sum(A intersection B)/sum( A union B)
end
% Non nomalized ioami_jaccard has units ppm by micromoles -2 sec-1
% Normalized in correlation form is unitless
% Note both are positive definite