%-*-octave-*-
% shotGun Simulation 
%    G = Laenge der DNA-Sequenz
%    Lmin = Minimale Laenge eines Fragments
%    Lmax = Maximale Laenge eines Fragments
%    N = Anzahl Fragmente ( > 1)
%    M = Anzahl Wiederholungen
%    vis = !=0: Ausgaben an     0: Ausgaben aus
%    [ fc, fl ] = [ Mittlere Contig Anzahl, Mittlere Contig Laenge ]
function [FC,FL] = shotGun(G,Lmin,Lmax,N,M,vis)

  Lm = (Lmax+Lmin)/2;

  %%%%% -- A Fragmentanfang -- L Fragmentlänge
  r = rand( N,M );
  A = sort(floor( r.*(G-Lmin) ) .+ 1);
  
  r = rand( N,M );
  L = floor( r.*(Lmax-Lmin+1) ) .+ Lmin;
  
  % falls Fragment hinten rausragt: stutzen
  L = ((A + L) > G) .* (G-(A+L)) .+ L;
  
  %%%%% -- NC Anzahl contigs
  AD = [A(2:N,:);G+2 .* ones(1,M)] .- A; % Differenz zum nächstem Fragment
  CEI = (AD .- L) > 0;    % 1 falls das Fragment ein contig Ende ist (Contig Ende Index in der Liste)
  NC = sum(CEI);

  %%%%% -- Laenge der contigs
  CL = zeros(N,M);
  MCL = zeros(M,1); % Mittelwert der contig Laengen pro Versuch
  for j = 1:M
    % find gibt die != 0 indices einer Matrix
    i = find(CEI(:,j));
    % Contig Ende
    CE = A(i,j) .+ L(i,j);
    i = [0;i];
    A2=[A(:,j);0];
    % Contig Anfang
    CA = A2(i+1); CA = CA(1:length(CA)-1);
    CL(1:NC(j),j) = CE .- CA;
    MCL(j,1) = mean( CL(1:NC(j),j) );
  end

  FC = mean(NC);
  FL = mean(MCL);

  if (vis != 0)
    amin=Lmin*N/G; amax=Lmax*N/G; a=Lm*N/G;
    printf("G = %d,  [ Lmin,Lmax ] = [ %d, %d ] (Lmean=%d),  N = %d\n",G,Lmin,Lmax,Lm,N);
    printf("Coverage = [ %f, %f ] (L=%f)\n", Lmin*N/G, Lmax*N/G, Lm*N/G );
    printf("Contigs: mittlere Anzahl = %f   \t\t( erwartet: [%f,%f](%f) )\n", mean(NC), N*exp(-amin), N*exp(-amax), N*exp(-a));
    printf("         mittlere Laenge = %f   \t( erwarted: [%f,%f](%f) )\n", mean(MCL), G/N*(exp(amin)-1), G/N*(exp(amax)-1), G/N*(exp(a)-1) );
  end
end
