#!/usr/local/bin/perl
#
# Read in the data file
# Print out a HTML formatted lines
#

$nn=4; #nn = tamanho da amostra no grafico de control
$yinicial=1000;
$y1=40;
$c2=2;
$c3=3;

$file = 'espectro.txt' ;                 
open(INFO, "<$file") ;                
@espectro = <INFO> ;                   
close(INFO) ;                        

$k=0;
$kk=0;

# print ("1: @espectro \n");
   
#lee Z y demas hiervas

foreach $linea  (@espectro)        # assign @espectro to $line, one at a time
{                                    
$k++;
 @l_sep= split ( /=/ , $linea ) ;

  
 if($k<=15){

     if(@l_sep[0] eq 'xmin '){$xmin =@l_sep[1];}
    if(@l_sep[0] eq 'xmax '){$xmax =@l_sep[1];}
    if(@l_sep[0] eq 'nx '){$nx =@l_sep[1];}
    if(@l_sep[0] eq 'dx '){$dx =@l_sep[1];}
    if(@l_sep[0] eq 'x1 '){$x1 =@l_sep[1];}
    if(@l_sep[0] eq 'ymin '){$ymin =@l_sep[1];}
    if(@l_sep[0] eq 'ymax '){$ymax =@l_sep[1];}
    if(@l_sep[0] eq 'ny '){$ny =@l_sep[1];}
    if(@l_sep[0] eq 'dy '){$dy =@l_sep[1];}
    if(@l_sep[0] eq 'y1 '){$y1 =@l_sep[1];}
 }
#print("xmin:$xmin, xmax:$xmax \n
#ymin:$ymin, ymax:$ymax \n
#nx:$nx,   ny:$ny \n
#dx:$dx,   dy:$dy
#x1:$x1,   y1:$y1");
 if($k>=16){
        if(@l_sep[1]ne''){
	    $mediaz=$mediaz+@l_sep[1];
	    $yy=int($kk/$nx)+1;
	    $xx=$kk-($yy-1)*$nx+1;
	    @z[$kk]=@l_sep[1];
					$kk2=$xx-1+($yy-1)*$nx;
# print ("$kk=$kk2: $xx,$yy,@z[$kk], @l_sep[0], 1:@l_sep[1] \n");
	    $kk++;
     }
 }
}

$y0=int($yinicial/$dy);

 print ("$y0 ");

#Media de todo...
$mediaz=$mediaz/$kk;



#saca las medias xb y los desvios s

for($k=1;$k<=$nx;$k++)
{
    for($l=$y0;$l<=$ny;$l++){
	@xb[$k]=@xb[$k]+@z[$k-1+($l-1)*$nx];
	@s[$k]=@s[$k]+(@z[$k-1+($l-1)*$nx])**2;
    }
    @xtotal[$k]=@xb[$k];
    @xb[$k]=@xb[$k]/$ny;
    @s[$k]=@s[$k]-(@xb[$k])**2;

 #  print("k=$k,xb=$xb[$k],s=$s[$k]\n");
}


#calcula las frecuencias para cada tripla posible para cada columna
for($k=1;$k<=$nx;$k++)
{
    $mediac=0;
    $mediac1;
    for($l=$y0;$l<=($ny);$l++){
	$mediac=$mediac+@z[$k-1+($l)*$nx];
	$mediac1=$mediac1+@z[$k+($l)*$nx];
	$media01=$mediac+$mediac1;
	if($l>=$40){@hf[$k]=@hf[$k]+@z[$k-1+($l)*$nx];}
    }
    for($l=$y0;$l<=($ny-3);$l++){
	@disp[$k]=@disp[$k]+(@z[$k-1+($l)*$nx]/$mediac-@z[$k+($l)*$nx]/$mediac1)**2;
	@loc[$k]=@loc[$k]+abs(@z[$k-1+($l)*$nx]-@z[$k+($l)*$nx])/$media01;
	if(@z[$k-1+($l-1)*$nx]<=$mediac){$k1=0;}
	else{$k1=1};
 	if(@z[$k-1+($l)*$nx]<=$mediac){$k2=0;}
	else{$k2=1};
 	if(@z[$k-1+($l+1)*$nx]<=$mediac){$k3=0;}
	else{$k3=1};
	@emp[$k-1+($nx-1)*($k1+2*$k2+4*$k3)]++;
	#print("$k1,$k2,$k3--@emp[$k-1+($nx-1)*($k1+2*$k2+4*$k3)]\n");
	#print("$k,$mediaz,--$k1,$k2,$k3\n");
    }
    @disp[$k]=@disp[$k]**(.5);
    @disp3[$k]=(@disp[$k]+@disp[$k-1]+@disp[$k-2]);
}
# Entropias, entropias relativas y relativas a dospasos  usando los pesos de cada frecuencia (normalizados) como probabilidades
for($k=1;$k<=($nx-3);$k++){
    for($l=$y0;$l<=$ny;$l++){
	$pk=@z[$k-1+($l-1)*$nx]/@xtotal[$k];
	$pkp=@z[$k+($l-1)*$nx]/@xtotal[$k+1];
	$pk2p=@z[$k+1+($l-1)*$nx]/@xtotal[$k+2];
	$pk3p=@z[$k+2+($l-1)*$nx]/@xtotal[$k+3];
	@entro2[$k]=@entro2[$k]-$pk*log($pk);
	@entro2r[$k]=@entro2r[$k]+$pk*log($pk/$pkp);
	@entro2r2p[$k]=@entro2r2p[$k]+$pk*log($pk/$pk2p);
	@entro2r3p[$k]=@entro2r3p[$k]+$pk*log($pk/$pk3p);
	#print("entro2parcial $k,$l,($k+$ll),@entro2[$k]\n");
    }
#print("entro2 $k,$l,($k+$ll),@entro2[$k]\n");
}


#calcula la entropia relativa de cada columna respecto a la siguienbnte

for($k=1;$k<=($nx-1);$k++){
    for($l=0;$l<=7;$l++){
	$ll=($nx-1)*$l;
	#print("entro   $k,$l,($k+$ll),  @emp[$k+$ll] \n");
	if(@emp[$k+$ll]==0){@emp[$k+$ll]=.01}
	if(@emp[$k-1+$ll]==0){@emp[$k-1+$ll]=.01}
	#print("entro   $k,$l,($k+$ll),  @emp[$k+$ll] \n");
	@entropia[$k]=@entropia[$k]+(@emp[$k-1+$ll]/($ny-3-$y0))*log(@emp[$k-1+$ll]/@emp[$k+$ll]);
	#if(@emp[$k+$ll]!=0.001){print("@emp[$k+$ll],---,@entropia[$k]\n");}
    }
} 


open ( REP, ">entropias.txt" ) ;
for($k=1;$k<=$nx;$k++)
{
printf (REP  "$k||@entro2[$k],@entro2r[$k],@entro2r2p[$k],@entro2r3p[$k],@entropia[$k], @disp[$k], @loc[$k],@hf[$k],@disp3[$k]  \n");
#print( "entr[$k]=@entropia[$k] \n");
}
close(REP);


open ( REP, ">gentropias.txt" ) ;
for($k=3;$k<=$nx-3;$k++)
{
    $e=int(10*@entro2[$k]);
    $e1=int(10*@entro2r[$k]);
    $e2=int(10*@entro2r2p[$k]);
    $e3=int(10*@entro2r3p[$k]);
    $se123=int(10*(@entro2r[$k]+@entro2r2p[$k-1]+@entro2r3p[$k-2])/3);
    $rse123=5*(@entro2r[$k]+@entro2r2p[$k-1]+@entro2r3p[$k-2])/3;
    @suave[$k]=int($seviejo+$rse123);
    if(@disp3[$k-1]>1){$c2=4;
			$c3=7;
		    }
    if(@disp3[$k-1]<=1){$c2=2;
			$c3=3;
		    }
    if(@clase[$k-2]==0){
	if (@suave[$k-1]<$c3){@clase[$k-1]=0;}
	if (@suave[$k-1]>=$c3){
	    if  (@suave[$k]>=$c3){ @clase[$k-1]=1;}
	    if  (@suave[$k]<$c3){ @clase[$k-1]=0;}
}}
  if(@clase[$k-2]==1){
	if (@suave[$k-1]>$c2){@clase[$k-1]=1;}
	if (@suave[$k-1]<=$c2){
	    if  (@suave[$k]<=$c2){ @clase[$k-1]=0;}
	    if  (@suave[$k]>$c2){ @clase[$k-1]=1;}
	}}
 if(@hf[$k-1]<.001){@clase[$k-1]=1; }
    $seviejo=$rse123;
    $kd=2*$k/10000;
    $me=int(1000*@xb[$k]);
printf (REP "$kd | $e|| $e1 | $e2 | $e3 | $se123 | @suave[$k]|  @clase[$k-1]\n");
#print( "entr[$k]=@entropia[$k] \n");
}
close(REP);
 
$ca=@clase[3];
for($k=3;$k<=$nx-3;$k++)
{
if(@clase[$k]==$ca){$consec=$consec+1;}
if(@clase[$k] ne $ca)
{
    if($ca==1){
	if($consec<=4){
	    for($ind=$k-$consec;$ind<=$k-1;$ind++){@clase[$ind]=0;}
	}
    }
#    if($ca==0){
#	if($consec<=2){
#	    for($ind=$k-$consec;$ind<=$k-1;$ind++){@clase[$ind]=1;}
#	}
#    }
    $ca=@clase[$k];
    $consec=1;
}
}

$ca=@clase[3];
for($k=3;$k<=$nx-3;$k++)
{
if(@clase[$k]==$ca){$consec=$consec+1;}
if(@clase[$k] ne $ca)
{
#    if($ca==1){
#	if($consec<=3){
#	    for($ind=$k-$consec;$ind<=$k-1;$ind++){@clase[$ind]=0;}
#	}
#    }
#    if($ca==1){
#	if($consec<=3){
#	    for($ind=$k-$consec;$ind<=$k-1;$ind++){@clase[$ind]=0;}
#	}
#    }
    if($ca==0){
	if($consec<=2){
	    for($ind=$k-$consec;$ind<=$k-1;$ind++){@clase[$ind]=1;}
	}
    }
    $ca=@clase[$k];
    $consec=1;
}
}



$ca=@clase[3];
$cp=$ca;            #OJO!!!!!!!!
for($k=3;$k<=$nx-3;$k++)
{
if(@clase[$k]==$ca){$consec=$consec+1;}
if(@clase[$k] ne $ca)
{

    if($ca==0){
	if($consecp>=7){
	    if($consec<=7){
	    for($ind=$k-$consec;$ind<=$k-1;$ind++){@clase[$ind]=1;}
	
	}
	}
    }
    $cp=$ca;
    $consecp=$consec;
    $ca=@clase[$k];
    $consec=1;
}
}



open ( REP, ">clases.txt" ) ;
for($k=3;$k<=$nx-3;$k++)
{
$kd=2*$k/10000;
$su=int( (@suave[$k]+ @suave[$k-1]+ @suave[$k-2])/3);
printf (REP "$kd || @suave[$k],  $su, @clase[$k]\n");
}
close(REP);

for($k=1;$k<=$nx;$k++){$tvocal=$tvocal+@clase[$k];}
$pvocal=$tvocal/$nx;

$nc=-1;
open ( REP, ">>clas.txt" ) ;
open ( REV, ">>dvog.txt" ) ;
open ( REC, ">>dcon.txt" ) ;
$ca=@clase[3];
for($k=3;$k<=$nx-10;$k++)
{
$kd=2*$k/1000;
if(@clase[$k]!=$ca){
    if (@clase[$k]==0){	
	$nc++;
	@con[$nc]=$kd-$ncon;
	$nvog=$kd;
	printf (REP "$kd $nc @con[$nc] c \n");
	if($nc>=1){printf (REC "@con[$nc] \n");}
    }
    if (@clase[$k]==1){	
	$nv++;
	@vog[$nv]=$kd-$nvog;
	$ncon=$kd;
	printf (REP "$kd $nv @vog[$nv] v \n");
	printf (REV "@vog[$nv] \n");}
    $ca=@clase[$k];
}
}
close(REC);
close(REV);
close(REP);

for($k=1;$k<=$nv;$k++){$tvog=$tvog+@vog[$k];}
for($k=1;$k<=$nc;$k++){$tcon=$tcon+@con[$k];}
$pvog=$tvog/($tvog+$tcon);
$mcon=$tcon/$nc;
for($k=1;$k<=$nc;$k++){$dcon=$dcon+(1-@con[$k]/$mcon)**2;}
$dcon=($dcon/($nc-1))**(.5);
for($k=1;$k<=$nc;$k++){$dconb=$dconb+($mcon-@con[$k])**2;}
$dconb=($dconb/($nc-1))**(.5);
$cvar=$dconb/$mcon;
open ( REP, ">>pv.txt" ) ;
printf (REP "$pvog\n");
close(REP);
open ( REP, ">>deltac.txt" ) ;
printf (REP "$dcon\n");
close(REP);
open ( REP, ">>sigmac.txt" ) ;
printf (REP "$dconb\n");
close(REP);
open ( REP, ">>coefvar.txt" ) ;
printf (REP "$cvar\n");
close(REP);