Qmatrix[D_] :=
  Module[{n, P, i, j, Q},
   n = Length[D[[1]]];
   P = D - DiagonalMatrix[Diagonal[D]];
   P = Table[(n - 2)*P[[i]][[j]] - Total[P[[i]]] - Total[P[[j]]], {i, 
      n}, {j, n}];
   Q = P - DiagonalMatrix[Diagonal[P]];
   Q = Q + DiagonalMatrix[ConstantArray[10^15, n]];
   Q
   ];

SeedRandom[165];
BaggageStatisticsChoice[Qmins_, CurrntPart_] := 
  Module[{baggageStats, minQIndices, baggageSum, baggMax, ext},
   baggageStats = Map[Length, CurrntPart];
   If[Length[Qmins] > 2,
    (*THEN*)
    minQIndices = DeleteDuplicates[Map[Sort, Qmins]];
    baggageSum = Map[Total[baggageStats[[#]]] &, minQIndices];
    baggMax = Position[baggageSum, Max[baggageSum]];
    If[Length[baggMax] > 1,
     ext = minQIndices[[RandomChoice[baggMax]]];,
     ext = minQIndices[[First[baggMax]]] ;
     ];,
    (*ELSE*)
    ext = Sort[{Sort[Qmins][[1, 1]], Sort[Qmins][[1, 2]]}];
    ];
   Flatten[ext]
   ];

NJ[dstart_] := 
  Module[{Q, d, n, i, CurrentPartition, PartitionList, minQ, 
    minQIndices, baggageSum, baggMax, Index1, Index2, Part1, Part2, 
    NewVec},
   d = N[dstart];(* Abr: 
   seems to me that this N function is doing nothing *)

   Q = Qmatrix[d];
   n = Length[d];
   CurrentPartition = {};
   Do[
    CurrentPartition = Append[CurrentPartition, {i}]
    , {i, 1, n}]; 
   PartitionList = {CurrentPartition};
   Do[
    minQ = Position[Q, Min[Q]];(* 
    minQ stores all the indices in Q where the minimum occurs *)
    {Index1, Index2} = BaggageStatisticsChoice[minQ, CurrentPartition];
    Part1 = CurrentPartition[[Index1]];
    Part2 = CurrentPartition[[Index2]];
    CurrentPartition =   
     Append[Drop[Drop[CurrentPartition, {Index2}], {Index1}],
      Sort[Flatten[{Part1, Part2}]]]; 
    NewVec = (1/2)*(d[[Index1]] + 
        d[[Index2]] - ConstantArray[d[[Index1]][[Index2]], Length[d]]);
    d = Append[d, NewVec];
    d = Append[Transpose[d], Append[NewVec, 10^15]];
    d = Drop[Drop[d, {Index2}] , {Index1}];
    d = Drop[Drop[Transpose[d], {Index2}] , {Index1}];
    PartitionList = Append[PartitionList, CurrentPartition];
    Q = Qmatrix[d];,
    {i, 1, n - 3}];
   PartitionList
   ];

RandomEpsilonConeVolumeDiscover[numsamples_, numtaxa_, epsilon_] := 
  Module[{origin, infmatrix, treecountvector, tempmat, newtree, 
    treepos, newtempmat, treelist, record1, NormalizationConstant, i, 
    j, l, m},
   origin = ConstantArray[0, {numtaxa, numtaxa}];
   infmatrix = 
    DiagonalMatrix[
     ConstantArray[10^15, 
      numtaxa]];(*This assures that NJ will never choose a diagonal entry as the smallest*)
   treelist = {};
   treecountvector = {};
   record1 = {};
   Do[
    tempmat = RandomReal[NormalDistribution[0, 1], {numtaxa, numtaxa}];
    tempmat = UpperTriangularize[tempmat, 1];
    NormalizationConstant = Norm[Flatten[tempmat]];
    tempmat = tempmat*(epsilon/NormalizationConstant);
    newtempmat = 
     origin + infmatrix + Abs[tempmat + Transpose[tempmat]];
    newtree = NJ[newtempmat];
    If [Check[NJ[newtempmat], err] === err, 
     record1 = Append[record1, newtempmat]];
    treepos = Flatten[Position[treelist, newtree]];
    If[treepos === {},
     treelist = Append[treelist, newtree];
     treecountvector = Append[treecountvector, 1],
     treecountvector[[treepos[[1]]]] = 
      treecountvector[[treepos[[1]]]] + 1]
    , {i, 1, numsamples}];
   {treelist, treecountvector, record1}
   ];


(*This works: the convention is TimingAlgorithmNumsamplesNumtaxa*)


Data4taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 4, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_4.dat", Data4taxa, "CSV"] 

Data5taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 5, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_5.dat", Data5taxa, "CSV"] 

Data6taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 6, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_6.dat", Data6taxa, "CSV"] 

Data7taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 7, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_7.dat", Data7taxa, "CSV"] 

Data8taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 8, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_8.dat", Data8taxa, "CSV"] 

Data9taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 9, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_9.dat", Data9taxa, "CSV"] 

Data10taxa = Timing[RandomEpsilonConeVolumeDiscover[1000000, 10, 1]];

Export["/home/abraham.mc/Ruth/BaggageNJ/TimingBaggNJ1000000_10.dat", Data10taxa, "CSV"] 


Quit[]