#!/usr/bin/env perl 

#
# Simula una ciudad con densidad tiempo de inyección y distancia de
# recorrido dados 
#

=encoding UTF-8

=head1 SYNOPSIS

./td.pl -size N -cuadra C -tiempo t -distancia d -pvuelta p -tmax tmax -tomaturnos -despliega -imprime -total

    simula una ciudad de tamaño N con cuadras de tamaño C donde se
    inyectan coches con constante de tiempo t y se remueven con constante
    de distancia d que dan vuelta con probabilidad p. La simulación dura
    un tiempo tmax (en unidades de 1/velocidad). Despliega los resultados en
    una simulación y/o imprime una tabla de densidad vs. tiempo o sólo
    los resultados totales   

=head1 AUTHOR

    W. Luis Mochán

=head1 LICENCE

 This work is licensed under a Creative Commons Attribution 4.0
 International License. 

=cut

use PDL;
use PDL::NiceSlice;
use PDL::Graphics::TriD;
use Time::HiRes qw(usleep);
use Getopt::Long;
use Pod::Usage;
use strict;
use warnings;

# parametros
my $size=200;
my $cuadra=10; # tamaño retícula de cuadras
my $tau=200; #tiempo caracteristico (p.ej 2 viajes al día=> medio día)
my $recorrido =140; #distancia de recorrido caracteristica
my $pvuelta=.1; # probabilidad de dar vuelta. 
my $Tmax=10000; # tiempo maximo

#banderas
my $despliegaflag=0;
my $imprimeflag=0;
my $totalflag=0;
my $tomaturnoflag=0;
my $helpflag=0;

pod2usage unless GetOptions(
    "size=i"=>\$size,
    "cuadra=i"=>\$cuadra,
    "tiempo=f"=>\$tau,
    "distancia=f"=>\$recorrido,
    "pvuelta=f"=>\$pvuelta,
    "tmax=i"=>\$Tmax,
    "despliega"=>\$despliegaflag,
    "imprime"=>\$imprimeflag,
    "total"=>\$totalflag,
    "tomaturnos"=>\$tomaturnoflag,
    "help|?"=>\$helpflag,
    ) and !$helpflag;

die "Nada que hacer: indica si quieres desplegar o imprimir o totales" unless
    $despliegaflag or $imprimeflag or $totalflag; 


#constantes
#coches que caben en la calle.
my $capacidad=2*$size*$size/$cuadra-($size*$size)/($cuadra*$cuadra); 

# medidores
my $tiempo=0; #tiempo transcurrido
my $porsalir=0; #coches esperando salir de su casa.
my $presentes=0; #coches en la ciudad
my $densidad=0; #coches/capacidad
my $tperdido=0;
my $tperdidoacum=0;

# escenario
my $pantalla;
my $ciudad;
if($despliegaflag){
    # la pantalla es para que coches vecinos no se toquen.
    $pantalla=zeros(2*$size,2*$size);
    $ciudad=$pantalla(0:-1:2,0:-1:2);
} else { #save some memory
    $ciudad=zeros($size,$size);
}

# prepara calles: incrementar primer indice avanza sobre las calles
my $ns=$ciudad(0:-1:2*$cuadra, -1:0)->transpose; # calles norte sur
my $sn=$ciudad($cuadra:-1:2*$cuadra, 0:-1)->transpose; # calles sur norte
my $ew=$ciudad(-1:0, 0:-1:2*$cuadra); # calles este oeste
my $we=$ciudad(0:-1,$cuadra:-1:2*$cuadra); # calles oeste este
my @calles=($ns,$sn,$ew,$we);


#srand(0); #temporal para pruebas

while($tiempo<$Tmax){
    if($tomaturnoflag){ #Toman turno todas las direcciones
	map 
	{
	    pueblayavanza($_); #avanza coches de una calle
	    print "$tiempo $densidad $tperdido $porsalir\n"
		if $imprimeflag and $tiempo % 10 == 0; 
	    despliega() if $despliegaflag;
	} 
	(0..3);
    }else{ # avanzan unos u otros al azar
	pueblayavanza(floor(rand(4)));
	    print "$tiempo $densidad $tperdido $porsalir\n"
		if $imprimeflag and $tiempo % 10 == 0; 
	despliega() if $despliegaflag;
    }
}
print $tperdido/$tiempo, "\n";

sub pueblayavanza { # 
    my $idx=shift; #identifica la calle.
    my $calle=$calles[$idx];
    # añade coches a la pila por salir con constante de tiempo $tau
    $porsalir += $capacidad/4/$tau;
    my $disponibles=!$calle; #calcula donde puede añadir coches
    my $cuantosdisponibles=$disponibles->sum;
    $disponibles*=random($calle)<($porsalir/$cuantosdisponibles)
	if $cuantosdisponibles;
    # $disponibles es una máscara de donde añadir coches
    # añádelos
    my $cuantosanhadidos = $disponibles->sum;
    $porsalir -= $cuantosanhadidos;
    $presentes += $cuantosanhadidos;
    $calle += $disponibles*($idx+1);

    # identifica qué coches darán vuelta hacia esta calle
    my $r=($calle!=0) & (random($calle) < $pvuelta);
    $calle.=($idx+1)*$r+$calle*!$r;

    # identifica quienes pueden avanzar
    my $estos=$calle==$idx+1; # los que van en esta calle
    my $atravezados=($calle!=0) & ($calle!=$idx+1); # los que van en otra dirección
    my $rezagados=($estos*$calle->rotate(-1))!=0; #chocarían de avanzar
    my $avanzados=($estos & !$rezagados)->rotate(1); 

    #identifica cuales ya llegaron a su destino
    my $arribados = $avanzados*(random($calle) < 1/$recorrido);
    #y eliminalos
    $avanzados *= !$arribados;
    $presentes -= $arribados->sum;

    #actualiza
    $tiempo+=0.25;
    $calle*=$atravezados+$rezagados;
    $calle+=$avanzados*($idx+1);
    $densidad = ($porsalir+$presentes)/$capacidad;
    $tperdido = $densidad*$tau-$recorrido;
    $tperdidoacum += $tperdido;
}

sub despliega {
    	  usleep(50000); 
	  nokeeptwiddling3d(); 
	  imagrgb [$pantalla];
}

__END__
  
