function hp_samples = gen_hp_samples(w,mu,A, num_samples,nc)
tic;
hp_samples = struct('timestamp', [],'event', []);
for n = 1:num_samples
    t=0;
    timestamp_and_event = [];
    lambdat = mu;                               % mu is the 1000 dim vector base intensity
    lambdat_sum = sum(lambdat);                 % Sum of all base intensities,used for gen next time-step
    while size(timestamp_and_event, 2) < nc
        rand_u = rand;
        dt = random('exp', 1/lambdat_sum);      % For generating the next time step
        lambda_ts = comp_lambda(t+dt, [], t, lambdat,w,mu,A);
        lambdats_sum = sum(lambda_ts);
        if rand_u > lambdats_sum / lambdat_sum  % Discard sample
            t = t+dt;
            lambdat = lambda_ts;
        else                                    % Keep sample
            u = rand * lambdats_sum;
            lambda_sum = 0;
            for d=1:length(mu)                  % d=1:1000
                lambda_sum = lambda_sum + lambda_ts(d);%sum over all lamda one by one
                if lambda_sum >= u              % Decide which dim event occurs
                    break;                      %end if
                end
            end
            lambdat = comp_lambda(t+dt, d, t, lambdat,w,mu,A);
            t = t+dt;
            timestamp_and_event = [timestamp_and_event,[t;d]];%append new sample
        end
        lambdat_sum = sum(lambdat);            % Time difference for next loop
    end                                        % While size(timestamp_and_event, 2)<nc
    hp_samples(n).timestamp = timestamp_and_event(1,:);
    hp_samples(n).event = timestamp_and_event(2,:);   
    if mod(n, 10)==0                           % Print out the generation process
        fprintf('#samples=%d/%d, time=%.2fsec\n',n, num_samples, toc);
    end
end                                            %for n = 1:num_samples