function [new_I_t,new_S_t] = SIS_hypergraph_disease_iterate(beta,delta,I_t,S_t,I,dt,K,fun)
% beta: rate of infection
% delta: recovery rate
% I_t: infectious population at some time step
% (output) new_I_t: infectious population at next time step
% S_t: susceptible population at some time
% (output) new_S_t: susceptible population at next time step
% f: non-linear dependence on the number of infectious individuals of infection rate within hyperedges
% I: incidence matrix
% dt: small time step

n = size(I_t,1);
S_idx = find(S_t==1);
I_idx = find(I_t==1);
new_S_idx = S_idx;

if length(delta)==1
    delta = delta*ones(n,1);
end

%% infection trial
s = length(S_idx);
for i=1:s
    hyperedges=find(I(S_idx(i),:)==1); % hyperedges containing S_idx(i)
    if gets_infected(hyperedges,I_idx,I,beta,dt,K,fun)
            new_S_idx = setdiff(new_S_idx,S_idx(i));
    end
end

%% recovery trial
l = length(I_idx);
for i=1:l
    is_susceptible = 0;
    if rand(1)<1-exp(-delta(I_idx(i))*dt)
        is_susceptible = 1;
        %hyperedges = find(I(I_idx(i),:)==1);
        %h = length(hyperedges);
        %for e=1:h
        %    if gets_infected(find(I(:,hyperedges(e))==1),I_idx,beta) && rand(1)<1/2
        %        is_susceptible = 0;
        %        break;
        %    end
        %end
        if is_susceptible
            new_S_idx = union(new_S_idx,I_idx(i));
        end
    end
end

new_S_t = zeros(n,1);
nS = size(new_S_idx,1);
new_S_t(new_S_idx,1) = ones(nS,1);
new_I_t = ones(n,1)-new_S_t;

end