log tnt.log;
macro=;
taxname+1000;
taxname=;
mxram 10240;
nstates 32;
nstates NOGAPS;

/*Arguments*/
if ( argnumber != 1 )
  silent -console;
  quote -   /-----------------------------------------------\;
  quote -   |        GUOYI TNT SCRIPT 2022-2023 MIT         |;
  quote -   |        You need to give your filename         |;
  quote -   |           tnt*> guoyi filename                |;
  quote -   \-----------------------------------------------/;
  proc/; 
end

/*Report what will be done*/
quote -   /-----------------------------------------------\;
quote -   |    Implied weighting will be estimated.        |;
quote -   |    TBR Mult and Xmult will be performed.       |;
quote -   |    Strict consensus will be used.              |;
quote -   |    Relative bremer support, jackknifing and    |;
quote -   |        boot will be shown on the resample.svg. |;
quote -   |    Apomorphic characters mapping will be shown |;
quote -   |        on the apo.svg and saved to apo*.tre.   |;
quote -   |    TL, CI and RI will be calculated finally.   |;
quote -   \------------------------------------------------/;


/*Set K*/
procedure %1;
var:
 actK
 minK
 maxK
 gVal
 a
 b
 Nval
 maxIts
;
if (ntax<0)
  quote There is no data;
  proc/;
end
/** find the value of g **/
loop 0 nchar
  if (!isinfo[#1]) continue; end
  if (!isact[#1]) continue; end
  set a=maxsteps[#1];
  set b=minsteps[#1];
  if (('a'-'b')>'gVal')
    set gVal 'a'-'b';
  end
stop
quote g value: 'gVal';
set actK 10;
set minK 0;
set maxK 500; /** Assume 500 as the maximum possible value of K **/
set maxIts 0;
/** find the best value of K **/
loop 0 1
  quote actual k value: 'actK';
  set a 1-('actK'/('actK'+1));
  quote 'a';
  set b ('actK'/('actK'+'gVal'-1))-('actK'/('actK'+'gVal'));
  quote 'b';
  set Nval 'a'/'b';
  quote N value: 'Nval';
  if (('Nval'>14.8)&&('Nval'<15.2)) endloop; end /** N is in range **/
  if ('Nval'>15.2)
    set minK 'actK';
  else
    set maxK 'actK';
  end
  set a ('maxK'-'minK')/2;
  set actK 'minK' + 'a';
  set maxIts ++;
  if ('maxIts'==100) endloop; end /** emergency exit **/
  setloop 0;
stop
/** If out of emergency **/
if ('maxIts'==100)
  quote Iterations are not finished;
  quote Best K is 'actK';
  proc /;
end
quote K value: 'actK' (N='Nval');
piwe='actK';

/*Reopen tnt*/
procedure %1;
hold 10000;

/*Search trees*/
xpiwe(*;
mult=replic 1000 tbr hold 1000;
sect: slack 1000; 
xmult=hit 1000 noupdate replications 20 drift 10 ratchet 10 fuse 10 hold 1 keepall;
piwe&;

/*Export trees*/
export= trees.tre;
taxname-;
export= trees_no.tre;
taxname=;
tsave = trees.ctf;
tsave *= trees.tnt.tre;

/*Consensus tree*/
bbreak=tbr;
nelsen *;
tchoose/;

/*Export consensus tree*/
export - original.tre;
taxname-;
export - original_no.tre;
taxname=;
tsave = original.ctf;
tsave *= original.tnt.tre;

/*Get Bremer/jak/boot support*/
tchoose/;
ttags=;
sub 1; hold +1000; bbreak=fill;
sub 2; hold +2000; bbreak=fill;
sub 3; hold +3000; bbreak=fill;
sub 4; hold +4000; bbreak=fill;
sub 5; hold +5000; bbreak=fill;
sub 6; hold +6000; bbreak=fill;
sub 7; hold +7000; bbreak=fill;
sub 8; hold +8000; bbreak=fill;
bsupport[;
resample jak replications 1000;
resample boot replications 1000;

/*Export consensus tree with supports*/
ttags & resample.svg thickness 7 italics fontsize 15;
quote -         /----------------------------------------------\;
quote -                       resample tags start               ;
ttags/;
quote -                       resample tags stop                ;
quote -         \----------------------------------------------/;
export < resample.tre;
taxname-;
export - resample_no.tre;

/*Apomorphic characters*/
export = winclada.tre;
taxname =;
ttags-;
ttags=;
apo >0;
quote -         /----------------------------------------------\;
quote -                     apomorphy tags start                ;
ttags/;
quote -                     apomorphy tags stop                 ;
quote -         \----------------------------------------------/;
ttags & apo.svg thickness 7 italics fontsize 15;
export < apo.tre;
taxname-;
export - apo_no.tre;
ttags-;

/*Caulculate TL/CI/RI score*/
length *;
report-;
var-;
var =
 0 themin
 1 themax
 + CIs [ ( ntrees + 1 ) ]
 + RIs [ ( ntrees + 1 ) ]
 + this
;
set themin minsteps;
set themax maxsteps;
loop 0 ntrees
   progress #1 (ntrees+1) Calculating indices...;
   set this length[#1];
   set CIs[#1] 'themin'/'this'; /*CI=1 means no homoplasy*/
   set RIs[#1] ('themax'-'this')/('themax'-'themin'); /*RI=1 character fits perfetcly*/
   stop
progress/;
report=;
macfloat 3;
maketable CIs Consistency index;
maketable RIs Retention index;

/*Report*/
quote -         /----------------------------------------------\;
quote -         |      The analysis has been finished.         |;
quote -         |      The file `tnt.log` contains             |;
quote -         |          K, TL, CI and RI                    |;
quote -         |      The file `trees*.tre` contain           |;
quote -         |          trees found by mult and xmult       |;
quote -         |      The file `original*.tre` contain        |;
quote -         |          strict consensus tree without label |;
quote -         |      The file `resample*.tre` contain        |;
quote -         |          consensus tree with support         |;
quote -         |      The file `apo*.tre` contain             |;
quote -         |          tree with apomorphic character      |;
quote -         |      The file `*_no.tre` contain             |;
quote -         |          tree with `taxaname-`               |;
quote -         |      The file `*.ctf` tree file is           |;
quote -         |          only readable for TNT               |;
quote -         |      The file `*.tnt.tre` contain            |;
quote -         |          is the ctf file with taxname        |;
quote -         |      The file `resample.svg` contain         |;
quote -         |          strict consensus tree with          |;
quote -         |          relative bremer, jak and boot       |;
quote -         |          support on the tree                 |;
quote -         |      The file `apo.svg` contain the          |;
quote -         |          tree with apomorphy mapping         |;
quote -         |      The file `winclada.tre` can be          |;
quote -         |          converted by tnt2winclada           |;
quote -         \----------------------------------------------/;

/*Quit*/
zzz;