package ductsim;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.concurrent.BrokenBarrierException;
import java.util.concurrent.CyclicBarrier;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.ThreadLocalRandom;
import java.util.concurrent.locks.ReentrantLock;

public class Duct implements Runnable {

    public class Coordinate
    {
        public int x;
        public int y;
        public int z;

        public Coordinate(int x, int y, int z)
            { this.x = x; this.y = y; this.z = z; }
        public boolean equals(Coordinate c)
            { return (this.x == c.x) && (this.y == c.y) &&(this.z == c.z); }
    }

    public static final int BASAL          = 1;
    public static final int PROG_STEM      = 2;
    public static final int PROG_BIPOTENT  = 3;
    public static final int PROG_MYOEPI    = 4;
    public static final int PROG_LUMEPI    = 5;
    public static final int MYOEPI         = 6;
    public static final int LUMEPI         = 7;

    public static final int gene_length = 32;

    private static int[][][] duct;

    private static int     length;
    private static int     radius;
    private static int     diameter;
    private static boolean HGP;           // Indica presencia o no de HGP
    private static int     capacity;      // Capacidad total del ducto
    private static int     stem_capacity; // Capacidad de las dos capas principales del ducto

    private static final int mutation_rate = 5; // Representado como mutation_rate/1000

    private static int     generation_number;
    private static boolean stem_initialized;  // Indica si se han colocado las células madre
    private static boolean initialized;       // Indica si se ha inicializado el ducto
    private boolean        local_initialized;

    private static int total_population;            // Contadores de población
    private static int dysfunctional_population;
    private static int total_reproductions;         // Contadores de reproducción
    private static int dysfunctional_reproductions;
    private static int stem_counter;                // Contadores de mutación
    private static int bipotent_counter;
    private static int proglumepi_counter;
    private static int lumepi_counter;
    private static int mut_stem_counter;
    private static int mut_bipotent_counter;
    private static int mut_proglumepi_counter;
    private static int mut_lumepi_counter;

    private static boolean DCIS_achieved;   // Indica si se ha alcanzado DCIS o no
    private static int     DCIS_generation; // Generación en que se generó DCIS

    private static int n_threads; // nº de hilos
    private static ExecutorService threadPool;
    private static CyclicBarrier   threadBarrier;
    private static Duct[]          threadArray;
    private static ReentrantLock[] thread_locks;
    private static final ReentrantLock   counter_lock = new ReentrantLock();

    private static int total_generations;
    private Integer[]  index_array;
    private int        thread_index;
    private Coordinate boundaries;


    /**
     * Constructor del ducto
     *
     * @param l Longitud del ducto (número de secciones)
     * @param r Radio del ducto
     * @param h Presencia o no de HGP
     * @param n Número de hilos
     */
    public Duct(int l, int r, boolean h, int n)
    {
        if (threadPool != null)
        {
            threadPool.shutdown();
            while(!threadPool.isTerminated()){}
        }

        threadPool    = null;
        threadBarrier = null;
        threadArray   = null;
        n_threads     = n;
        duct = new int[l][r*2][r*2];
        length   = l;
        radius   = r;
        HGP      = h;
        diameter = r*2;
        capacity = 0;
        stem_capacity = 0;

        generation_number = 0;
        stem_initialized  = false;
        initialized       = false;

        DCIS_achieved   = false;
        DCIS_generation = 0;
        index_array = new Integer[length];

        for(int z=0 ; z<length ; z++)
        {
            index_array[z] = z;
            for(int x=0 ; x<=(r/2) ; x++)
            {
                for(int y=x ; y<=r ; y++)
                {
                    if (get_point_radius(x,y) == r)
                    {
                        generate_cells_8way(z, y, x, BASAL);
                    }
                }
            }
        }

        for(int i=0 ; i<diameter ; i++)
        {
            for(int j=0 ; j<diameter ; j++)
            {
                if(get_point_radius(i,j) < radius)
                    capacity++;
                if(get_point_radius(i,j) == radius-1 || get_point_radius(i,j) == radius-2)
                    stem_capacity++;
            }
        }
        capacity *= length;
        stem_capacity *= length;
        threadBarrier = new CyclicBarrier(n_threads+1); // Hilos que lanzamos mas éste, el principal
        threadArray = new Duct[n_threads];
        thread_locks  = new ReentrantLock[n_threads-1];

        int step = length/n_threads;
        for(int i=0 ; i<n_threads ; i++)
        {
            int begin =     i*step;
            int end   = (i+1)*step-1;
            if ((i+1) == n_threads)
                end = length-1;

            threadArray[i] = new Duct(begin, end, i);
            if (i < n_threads-1)
                thread_locks[i]  = new ReentrantLock();
        }
    }

    /* Constructor para crear hilos (begin y end inclusive en rango)*/
    public Duct(int begin, int end, int index)
    {
        thread_index = index;
        index_array = new Integer[end-begin+1];
        // Aquí Coordinate se usa como un par de números, no como una coordenada
        boundaries = new Coordinate(begin, end, 0);
        // [begin, end] ambos inclusive
        for(int z=begin, i=0 ; z<=end ; z++, i++)
        {
            index_array[i] = z;
        }
    }

    // [OPT] Cambiar las llamadas a función por acceso directo a la matriz
    // duct.[c.z][c.y][c.x] o algo así
    // en realidad para lo que se usa da igual
    public int  get_cell (Coordinate c)        { return duct[c.z][c.y][c.x]; }
    public void set_cell(Coordinate co, int c) { this.duct[co.z][co.y][co.x] = c; }
    public void erase_cell(Coordinate co)      { this.duct[co.z][co.y][co.x] = 0; }

    public void generate_cells_8way(int z, int y, int x, int type)
    {
        int r = radius;
        this.duct[z][x][y]         = Cell(type);
        this.duct[z][x][2*r-y]     = Cell(type);
        this.duct[z][2*r-x][y]     = Cell(type);
        this.duct[z][2*r-x][2*r-y] = Cell(type);
        // simetría diagonal
        this.duct[z][y][x]         = Cell(type);
        this.duct[z][y][2*r-x]     = Cell(type);
        this.duct[z][2*r-y][x]     = Cell(type);
        this.duct[z][2*r-y][2*r-x] = Cell(type);

    }

    public static int get_point_radius(int x, int y)
    {
        // [ALT]: Fórmula que se ajuste al círculo que se ve en el paper
        return (int)(Math.sqrt((x-radius)*(x-radius)+(y-radius)*(y-radius))+1);
    }

    public int Cell (int type)
    {
        return type * (int)1e8;
    }

    public int get_cell_type(int c)
    {
        return c/(int)1e8;
    }

    public int get_gene_housekeeping(int c)
    {
        return (c % (int)1e8) / (int)1e6;
    }

    public int get_gene_protooncogene(int c)
    {
        return (c % (int)1e6) / (int)1e4;
    }

    public int get_gene_supressor(int c)
    {
        return (c % (int)1e4) / (int)1e2;
    }

    public int get_gene_apoptosis(int c)
    {
        return c % (int)1e2;
    }

    public int set_cell_type(int c, int type)
    {
        return (c % (int)1e8) + type * (int)1e8;
    }

    // Muta 10% del genoma
    public int mutate_HGP(int c)
    {
        int counter = (int)((((float)gene_length) / 100f) * 10f);
        return c += counter * (int)1e6 + counter * (int)1e4 + counter * (int)1e2 + counter;
    }

    // Mutation Rate representado como por mil, no por cien
    public int mutate(int c, int mutation_rate)
    {
        int hk_c = gene_length - get_gene_housekeeping(c);
        int po_c = gene_length - get_gene_protooncogene(c);
        int sp_c = gene_length - get_gene_supressor(c);
        int ap_c = gene_length - get_gene_apoptosis(c);

        for(int i=0 ; i < hk_c ; i++)
        {
            if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
                c += (int)1e6;
        }
        for(int i=0 ; i < po_c ; i++)
        {
            if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
                c += (int)1e4;
        }
        for(int i=0 ; i < sp_c ; i++)
        {
            if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
                c += (int)1e2;
        }
        for(int i=0 ; i < ap_c ; i++)
        {
            if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
                c += 1;
        }

        return c;
    }

    // Método que usamos para ver si la célula es cancerígena según el paper
    public boolean is_considered_cancerous(int c)
    {
        return //get_gene_housekeeping(c)() >= 16 ||
               get_gene_protooncogene(c) >= 16 ||
               //get_gene_supressor(c)() >= 16 ||
               get_gene_apoptosis(c)     >= 16;
    }

    public int get_mutations(int c)
    {
        return get_gene_housekeeping(c)  +
               get_gene_protooncogene(c) +
               get_gene_supressor(c)     +
               get_gene_apoptosis(c);
    }

    public boolean is_dysfunctional(int c)
    {
        return get_gene_housekeeping(c)  >= 16 ||
               get_gene_protooncogene(c) >= 16 ||
               get_gene_supressor(c)     >= 16 ||
               get_gene_apoptosis(c)     >= 16;
    }

    public boolean housekeeping_is_dysfunctional(int c)
    {
        return get_gene_housekeeping(c) >= 16;
    }

    public boolean protooncogene_is_dysfunctional(int c)
    {
        return get_gene_protooncogene(c) >= 16;
    }

    public boolean supressor_is_dysfunctional(int c)
    {
        return get_gene_supressor(c)>= 16;
    }

    public boolean apoptosis_is_dysfunctional(int c)
    {
        return get_gene_apoptosis(c) >= 16;
    }

    public List<Float> mutations_stats()
    {
        List<Float> result = new ArrayList<>();
        result.add((float)mut_stem_counter/(float)stem_counter);
        result.add((float)mut_bipotent_counter/(float)bipotent_counter);
        result.add((float)mut_proglumepi_counter/(float)proglumepi_counter);
        result.add((float)mut_lumepi_counter/(float)lumepi_counter);
        return result;
    }

    public void update_counters(Coordinate c)
    {
        switch(get_cell_type(get_cell(c)))
        {
            case PROG_STEM:
                stem_counter++;
                mut_stem_counter += get_mutations(get_cell(c));
                break;
            case PROG_BIPOTENT:
                bipotent_counter++;
                mut_bipotent_counter += get_mutations(get_cell(c));
                break;
            case PROG_LUMEPI:
                proglumepi_counter++;
                mut_proglumepi_counter += get_mutations(get_cell(c));
                break;
            case LUMEPI:
                lumepi_counter++;
                mut_lumepi_counter += get_mutations(get_cell(c));
                break;
        }
    }

    protected void execute(int nGens)
    {
        if(generation_number == 0)
        {
            init_routine_stem(); // Hacemos la primera generación (colocación de stem cells) secuencialmente
        }

        total_generations = nGens;

        try {
        threadPool = Executors.newFixedThreadPool(n_threads);
        for (int i=0; i<threadArray.length; i++)
            threadPool.execute(threadArray[i]);
        } catch(Exception e){System.out.println(e);}

        for (int i=0; i<total_generations; i++)
        {
            total_population = 0;
            dysfunctional_population = 0;
            total_reproductions = 0;
            dysfunctional_reproductions = 0;

            stem_counter = 0;
            bipotent_counter = 0;
            proglumepi_counter = 0;
            lumepi_counter = 0;
            mut_stem_counter = 0;
            mut_bipotent_counter = 0;
            mut_proglumepi_counter = 0;
            mut_lumepi_counter = 0;

            try
            {
                threadBarrier.await();
                generation_number++;

                if(!initialized)
                {
                    initialized = threadArray[0].get_local_initialized();
                    for (int t=1; t<threadArray.length; t++)
                        initialized = initialized && threadArray[t].get_local_initialized();
                }
            }
            catch (BrokenBarrierException | InterruptedException e){}
        }
        threadPool.shutdown();
        while(!threadPool.isTerminated()){}
    }

    public void run()
    {
        for (int i=0; i<total_generations; i++)
        {
            try
            {
                next_generation();
                threadBarrier.await();
            }
            catch (BrokenBarrierException | InterruptedException e){}
        }
    }

    public int get_generation_number()           { return generation_number; }
    public int get_length()                      { return length; }
    public int get_radius()                      { return radius; }
    public int get_diameter()                    { return diameter; }
    public int get_capacity()                    { return capacity; }
    public boolean get_local_initialized()       { return local_initialized; }

    public int get_dysfunctional_population()    { return dysfunctional_population; }
    public int get_total_population()            { return total_population; }
    public int get_dysfunctional_reproductions() { return dysfunctional_reproductions; }
    public int get_total_reproductions()         { return total_reproductions; }
    public int[][] get_slice(int z)              { return duct[z]; }

    public int DCIS_achieved()
    {
        if(!DCIS_achieved)
            return 0;
        else
            return DCIS_generation;
    }

    /**
     * @param c Coordenada
     * @return true si la coordenada c se encuentra dentro del array ducto
     *         (no necesariamente dentro de la membrana basal)
     *         false en caso contrario
     */
    public boolean coordinate_is_inbounds(Coordinate c)
    {
        return (c.x >= 0) && (c.x < diameter) &&
               (c.y >= 0) && (c.y < diameter) &&
               (c.z >= 0) && (c.z < length);
    }

    /**
     * @param c Coordenada
     * @return ArrayList de Coordenadas con todos los vecinos (vecindad Moore 3D) de la cordenada c
     */
    private List<Coordinate> get_neighbors(Coordinate c)
    {
        List<Coordinate> neighbors = new ArrayList<>();
        int x = c.x;
        int y = c.y;
        int z = c.z;
        // Las 6 coordenadas "principales" (vecindad directa)
        neighbors.add(new Coordinate(x,   y,   z-1));
        neighbors.add(new Coordinate(x,   y,   z+1));
        neighbors.add(new Coordinate(x,   y-1, z));
        neighbors.add(new Coordinate(x-1, y,   z));
        neighbors.add(new Coordinate(x+1, y,   z));
        neighbors.add(new Coordinate(x,   y+1, z));
        // Vecinos diagonales
        neighbors.add(new Coordinate(x-1, y-1, z-1)); // diagonal
        neighbors.add(new Coordinate(x,   y-1, z-1)); // diagonal
        neighbors.add(new Coordinate(x+1, y-1, z-1)); // diagonal
        neighbors.add(new Coordinate(x-1, y,   z-1)); // diagonal
        neighbors.add(new Coordinate(x+1, y,   z-1)); // diagonal
        neighbors.add(new Coordinate(x-1, y+1, z-1)); // diagonal
        neighbors.add(new Coordinate(x,   y+1, z-1)); // diagonal
        neighbors.add(new Coordinate(x+1, y+1, z-1)); // diagonal
        neighbors.add(new Coordinate(x-1, y-1, z+1)); // diagonal
        neighbors.add(new Coordinate(x,   y-1, z+1)); // diagonal
        neighbors.add(new Coordinate(x+1, y-1, z+1)); // diagonal
        neighbors.add(new Coordinate(x-1, y,   z+1)); // diagonal
        neighbors.add(new Coordinate(x+1, y,   z+1)); // diagonal
        neighbors.add(new Coordinate(x-1, y+1, z+1)); // diagonal
        neighbors.add(new Coordinate(x,   y+1, z+1)); // diagonal
        neighbors.add(new Coordinate(x+1, y+1, z+1)); // diagonal
        neighbors.add(new Coordinate(x-1, y-1, z));   // diagonal
        neighbors.add(new Coordinate(x+1, y-1, z));   // diagonal
        neighbors.add(new Coordinate(x-1, y+1, z));   // diagonal
        neighbors.add(new Coordinate(x+1, y+1, z));   // diagonal

        return neighbors;
    }

    /**
     * @param cell_coord
     * @return Devuelve una coordenada aleatoria que cumple las siguientes condiciones:
     *          - Es vecina de cell_coord
     *          - Se encuentra dentro de la membrana basal
     *          - Está vacía (no está ocupada por una célula)
     *         Si tal coordenada no existe, devuelve null
     */
    private Coordinate get_vacant_neighbor(Coordinate cell_coord)
    {
        List<Coordinate> neighbors = get_neighbors(cell_coord);
        Collections.shuffle(neighbors);
        for(int i=0 ; i < neighbors.size() ; i++)
        {
            Coordinate aux_coord = neighbors.get(i);
            if(coordinate_is_inbounds(aux_coord) == true &&
               get_cell(aux_coord) == 0 &&
               get_point_radius(aux_coord.x, aux_coord.y) < radius)
            {
                return neighbors.get(i);
            }
        }
        return null;
    }

    /**
     * @param cell_coord Coordenada a la que se va a migrar
     * @param old_coord  Coordenada desde la que se migra
     * @return true si cell_coord tiene alguna célula vecina aparte de old_coord
     *         false en caso contrario
     */
    private boolean has_adjadcent_neighbor(Coordinate cell_coord, Coordinate old_coord)
    {
        List<Coordinate> neighbors = get_neighbors(cell_coord);
        for(int i=0 ; i < neighbors.size() ; i++)
        {
            Coordinate aux_coord = neighbors.get(i);
            if(coordinate_is_inbounds(aux_coord) == true &&
               get_cell(aux_coord) != 0 &&
               !aux_coord.equals(old_coord))
            {
                return true;
            }
        }
        return false;
    }

    private boolean reproduce(Coordinate current_cell_coord)
    {
        List<Coordinate> neighbors = get_neighbors(current_cell_coord);
        Collections.shuffle(neighbors);
        int current_cell_type = get_cell_type(get_cell(current_cell_coord));

        /*
        INIT:   Allow progenitor cells to reproduce into vacant neighboring lattice points
                without cells according to progenitor hierarchy
        NORMAL: if an adjacent neighboring point is vacant and does not break the bi-layer
                ductal structure (Fig. 3C), the cell reproduces
        */

        /* [ALT]: (a) Permitir que una célula mioepitelial se reproduzca en el sitio que sea, y luego
                      se corrija ésto en la migración, o bien
                  (b) Permitir que sólo se genere una célula en un sitio que le 'corresponde' [actual]
        */

        int max_radius = radius-1;
        int min_radius = radius-2;

        if(current_cell_type == PROG_LUMEPI)
            max_radius--;
        if(current_cell_type == PROG_MYOEPI)
            min_radius++;

        Coordinate new_cell_coord = null;
        boolean end_flag = false;
        for(int i=0 ; i<neighbors.size() ; i++)
        {
            // Condiciones: la célula está inbounds, vacía, y entre min_radius y max_radius (inclusive)
            new_cell_coord = neighbors.get(i);
            if(coordinate_is_inbounds(new_cell_coord) == true &&
               get_cell(new_cell_coord) == 0 &&
               // (b)
               get_point_radius(new_cell_coord.x, new_cell_coord.y) <= max_radius &&
               get_point_radius(new_cell_coord.x, new_cell_coord.y) >= min_radius)
               // (a)
               //get_point_radius(new_cell_coord.x, new_cell_coord.y) <= radius &&
               //get_point_radius(new_cell_coord.x, new_cell_coord.y) >= radius-2)
            {
                end_flag = true;
                break;
            }
        }

        if(end_flag == false)
        {
            return false;
        }

        int parent = get_cell(current_cell_coord);
        int new_cell = parent;
        switch (get_cell_type(parent))
        {
            // [OPT] Creo que los switch se compilan como una cadena de if-else;
            // por lo que es mejor poner la opción más común primero
            case PROG_LUMEPI:
                new_cell = set_cell_type(new_cell, LUMEPI);
                break;
            case PROG_BIPOTENT:
                // [ALT] Decidimos qué tipo de célula va a generar una bipotente en base a
                // la posición en que se va a generar
                if (get_point_radius(new_cell_coord.x, new_cell_coord.y) == radius-1)
                    new_cell = set_cell_type(new_cell, PROG_MYOEPI);
                else
                    new_cell = set_cell_type(new_cell, PROG_LUMEPI);
                break;
            case PROG_STEM:
                new_cell = set_cell_type(new_cell, PROG_BIPOTENT);
                break;
            case PROG_MYOEPI:
                new_cell = set_cell_type(new_cell, MYOEPI);
                break;
            default:
                break;
        }
        // Mutamos tanto madre como hija
        parent = mutate(parent, mutation_rate);
        new_cell = mutate(new_cell, mutation_rate);
        set_cell(current_cell_coord, parent);
        set_cell(new_cell_coord, new_cell);

        return true;
    }

    private boolean cancerous_reproduce(Coordinate current_cell_coord)
    {
        Coordinate new_cell_coord = null;
        Coordinate pushed_position = null;

        List<Coordinate> neighbors = get_neighbors(current_cell_coord);
        Collections.shuffle(neighbors);
        boolean end_flag = false;
        for (int i=0 ; i<neighbors.size() ; i++)
        {
            pushed_position = null;
            new_cell_coord = neighbors.get(i);

            if (coordinate_is_inbounds(new_cell_coord) == true)
            {
                if (get_cell(new_cell_coord) == 0 &&
                    get_point_radius(new_cell_coord.x, new_cell_coord.y) < radius)
                {
                    end_flag = true;
                    break;
                }
                // [ALT] Sólo empujamos células LUMEPI (?)
                if (get_point_radius(new_cell_coord.x, new_cell_coord.y) < radius &&
                    (pushed_position = get_vacant_neighbor(new_cell_coord)) != null)
                {
                    end_flag = true;
                    break;
                }
            }
        }

        if(end_flag == false)
        {
            return false;
        }

        if(pushed_position != null)
        {
            int pushed_cell = get_cell(new_cell_coord);
            set_cell(pushed_position, pushed_cell);
            erase_cell(new_cell_coord);
        }

        int parent = get_cell(current_cell_coord);
        int new_cell = parent;
        switch (get_cell_type(parent))
        {
            case PROG_LUMEPI:
                new_cell = set_cell_type(new_cell, LUMEPI);
                break;
            case PROG_BIPOTENT:
                if (get_point_radius(new_cell_coord.x, new_cell_coord.y) == radius-1)
                    new_cell = set_cell_type(new_cell, PROG_MYOEPI);
                else
                    new_cell = set_cell_type(new_cell, PROG_LUMEPI);
                break;
            case PROG_STEM:
                new_cell = set_cell_type(new_cell, PROG_BIPOTENT);
                break;
            case PROG_MYOEPI:
                new_cell = set_cell_type(new_cell, MYOEPI);
                break;

            default:
                break;
        }
        parent =   mutate(parent, mutation_rate);
        new_cell = mutate(new_cell, mutation_rate);
        set_cell(current_cell_coord, parent);
        set_cell(new_cell_coord, new_cell);

        return true;
    }

    private boolean hormonal_reproduce(Coordinate current_cell_coord)
    {
        List<Coordinate> neighbors = get_neighbors(current_cell_coord);
        Collections.shuffle(neighbors);

        int current_cell_type = get_cell_type(get_cell(current_cell_coord));
        Coordinate new_cell_coord = null;

        boolean end_flag = false;
        Coordinate pushed_position = null;
        for (int i=0 ; i<neighbors.size() ; i++)
        {
            new_cell_coord = neighbors.get(i);
            if (coordinate_is_inbounds(new_cell_coord) != false &&
                get_point_radius(new_cell_coord.x, new_cell_coord.y) < radius &&
                get_cell(new_cell_coord) != 0 &&
                (pushed_position = get_vacant_neighbor(new_cell_coord)) != null
               )
            {
                end_flag = true;
                break;
            }
        }

        if(end_flag == false)
        {
            return false;
        }

        int pushed_cell = get_cell(new_cell_coord);
        set_cell(pushed_position, pushed_cell);
        erase_cell(new_cell_coord);

        int parent = get_cell(current_cell_coord);
        int new_cell = parent;
        switch (get_cell_type(parent))
        {
            case PROG_LUMEPI:
                new_cell = set_cell_type(new_cell, LUMEPI);
                break;
            case PROG_BIPOTENT:
                if (get_point_radius(new_cell_coord.x, new_cell_coord.y) == radius-1)
                    new_cell = set_cell_type(new_cell, PROG_MYOEPI);
                else
                    new_cell = set_cell_type(new_cell, PROG_LUMEPI);
                break;
            case PROG_STEM:
                new_cell = set_cell_type(new_cell, PROG_BIPOTENT);
                break;
            case PROG_MYOEPI:
                new_cell = set_cell_type(new_cell, MYOEPI);
                break;
            default:
                break;
        }
        parent   = mutate(parent, mutation_rate);
        new_cell = mutate(new_cell, mutation_rate);
        set_cell(current_cell_coord, parent);
        set_cell(new_cell_coord, new_cell);

        return true;
    }

    private boolean migrate(Coordinate current_cell_coord)
    {
        List<Coordinate> neighbors = get_neighbors(current_cell_coord);
        Collections.shuffle(neighbors);

        int current_cell = get_cell(current_cell_coord);
        int current_cell_type = get_cell_type(current_cell);

        int max_radius = radius-1;
        // [ALT]: No permitimos que una célula migre hacia "dentro"; siempre hacia un punto en el mismo radio
        // que el punto en el que está, o mayor
        int min_radius = get_point_radius(current_cell_coord.x, current_cell_coord.y);

        /*
        if(get_cell(current_cell_coord).apoptosis_is_dysfunctional())
        {
            // WARN: Arreglar la migración de células cancerígenas (asegurar que migran a un lugar conectado
            // con el nucleo principal de poblacion)
            //min_radius = 0;
        }
        */
        if(current_cell_type == PROG_LUMEPI || current_cell_type == LUMEPI)
            max_radius = radius-2;
            //if(current_cell_type == PROG_MYOEPI || current_cell_type == MYOEPI)
                //min_radius = radius-1;
        Coordinate new_cell_coord = null;
        boolean end_flag = false;
        for (int i=0 ; i<neighbors.size() ; i++)
        {
            new_cell_coord = neighbors.get(i);
            if (coordinate_is_inbounds(new_cell_coord) == true &&
                get_cell(new_cell_coord) == 0 &&
                get_point_radius(new_cell_coord.x, new_cell_coord.y) <= max_radius &&
                get_point_radius(new_cell_coord.x, new_cell_coord.y) >= min_radius &&
                has_adjadcent_neighbor(new_cell_coord, current_cell_coord))
            {
                end_flag = true;
                break;
            }

        }

        if(end_flag == false)
        {
            return false;
        }

        erase_cell(current_cell_coord);
        set_cell(new_cell_coord, current_cell);

        return true;
    }

    private void init_routine_stem()
    {
        if(radius == 1)
            return;

        // 1. Draw the BASAL membrane by designating outer boundary points in a cylindrical orientation

        /* Hecho en constructor de SliceSec */

        // 2. Place STEM CELLS within the basal membrane
        //    (less than 5% of the population - 200 for a 24,000 cell pop. => [ ~ 0.84% ])

        // [ALT]: Con nuestro modelo de radio, las células que forman las dos capas iniciales son 20,800 contra las
        //        24,000 del paper.
        float stem_cells_percentage = 1f;  // Porcentaje del total de células que queremos que sean células madre
        int n_stem_cells = (int)((((float)stem_capacity) / 100) * stem_cells_percentage);

        if(n_stem_cells == 0)
            n_stem_cells = 1;

        int chosen_x, chosen_y, chosen_z;
        // Elegimos posición de células madre por polling de puntos aleatorios
        for(int counter = 1 ; counter <= n_stem_cells ; counter++)
        {
            // Elegir SliceSec aleatorio
            chosen_z = ThreadLocalRandom.current().nextInt(0, length);
            // Recorrer slice, elegir punto aleatorio hasta encontrar uno de radio r-1
            do {
                chosen_x = ThreadLocalRandom.current().nextInt(0, diameter);
                chosen_y = ThreadLocalRandom.current().nextInt(0, diameter);
            } while (get_point_radius(chosen_x, chosen_y) != radius-1 ||
                     get_cell(new Coordinate(chosen_x, chosen_y, chosen_z)) != 0);

            // Establecer nueva célula basal
            Coordinate new_coordinate = new Coordinate(chosen_x, chosen_y, chosen_z);
            int new_stem_cell = Cell(PROG_STEM);

            // 3. If HGP, mutate all genes in all stem cells by 10%
            // (ver método en Cell.java)
            if(HGP)
            {
                // [ALT]: Mutación determinista (mutar 10% del genoma) o estocástica (mutar con 10% de probabilidad) ?
                new_stem_cell = mutate_HGP(new_stem_cell);
            }
            set_cell(new_coordinate, new_stem_cell);
        }
        total_population += n_stem_cells;
        stem_initialized = true;
    }

    private void init_routine()
    {
        if(!stem_initialized)
        {
            init_routine_stem();
            return;
        }

        boolean end_flag = true; // Indica si la inicialización del ducto ha acabado

        Integer[] index_array_aux = new Integer[index_array.length];
        System.arraycopy(index_array, 0, index_array_aux, 0, index_array.length);
        Collections.shuffle(Arrays.asList(index_array_aux));
        // Iteramos sobre todas las coordenadas posibles (abierto a optimizaciones)
        for(int z_i = 0 ; z_i < index_array_aux.length ; z_i++) {
            int z = index_array_aux[z_i];
            for(int y=0 ; y<diameter ; y++) {
                for(int x=0 ; x<diameter ; x++) {
                    Coordinate current_cell_coord = new Coordinate(x,y,z);

                    // Si la coordenada no contiene ninguna célula o una célula basal, pasamos a la siguiente
                    if(get_cell(current_cell_coord) == 0 || get_cell_type(get_cell(current_cell_coord)) == BASAL)
                        continue;

                    total_population++;
                    if(is_dysfunctional(get_cell(current_cell_coord)))
                        dysfunctional_population++;

                    int current_cell_type = get_cell_type(get_cell(current_cell_coord));
                    // 4. Allow PROGENITOR cells to REPRODUCE into VACANT NEIGHBORING lattice points WITHOUT CELLS
                    //    according to progenitor hierarchy
                    if (current_cell_type == PROG_STEM     ||
                        current_cell_type == PROG_BIPOTENT ||
                        current_cell_type == PROG_MYOEPI   ||
                        current_cell_type == PROG_LUMEPI)
                    {
                        if (reproduce(current_cell_coord) == true)
                        {
                            end_flag = false;
                            total_reproductions++;
                        }
                    }
                    update_counters(current_cell_coord);
                    // 5. Allow NON-STEM cells to MIGRATE 1 point to a VACANT NEIGHBORING lattice point consistent
                    //    with the bi-layer ductal structure
                    if(current_cell_type != PROG_STEM)
                    {
                        if (migrate(current_cell_coord) == true)
                            end_flag = false;
                    }
                }
            }
        }
        local_initialized = end_flag;
        if(local_initialized)
            System.out.println(Thread.currentThread().getName() + " inicializado en " + generation_number);
    }

    private void normal_routine()
    {
        for(int z_i = 0 ; z_i < index_array.length ; z_i++) {
            int z = index_array[z_i];
            for(int y=0 ; y<diameter ; y++) {
                for(int x=0 ; x<diameter ; x++) {

                    Coordinate current_cell_coord = new Coordinate(x,y,z);
                    //System.out.println(Thread.currentThread().getName() + " accede a " + x + "," + y + "," + z + "[gen" + generation_number + "]");
                    boolean left_boundary_flag  = false;
                    boolean right_boundary_flag = false;
                    boolean reproduction_flag = false;
                    boolean cancerous_reproduction_flag = false;

                    // Barrera izquierda
                    if(boundaries.x != 0 && z <= boundaries.x+1)
                    {
                        thread_locks[thread_index-1].lock();
                        left_boundary_flag = true;
                    }
                    // Barrera derecha
                    if(boundaries.y != length-1 && z >= boundaries.y-1)
                    {
                        thread_locks[thread_index].lock();
                        right_boundary_flag = true;
                    }

                    if(get_cell(current_cell_coord) == 0 || get_cell_type(get_cell(current_cell_coord)) == BASAL)
                    {
                        if(left_boundary_flag)
                            thread_locks[thread_index-1].unlock();
                        if(right_boundary_flag)
                            thread_locks[thread_index].unlock();
                        continue;
                    }

                    // Si la célula se encuentra en zona de conflicto, debe tomar el cerrojo de la partición contigua
                    // Cuando esto suceda, para evitar deadlock, la partición más a la izquierda tiene prioridad sobre el
                    // cerrojo de la partición a su derecha.

                    // Partición derecha quiere acceder a partición izquierda: partición derecha no tiene prioridad.
                    // Ceder cerrojo, y luego adquirir los dos cerrojos, tras lo que se re-comproará la condición inicial


                    int current_cell = get_cell(current_cell_coord);

                    // 1. If housekeeping or apoptosis genes are dysfunctional (50% mutated), mark the cell for death
                    // 3. If dead, remove from lattice
                    if(housekeeping_is_dysfunctional(current_cell) ||
                          (protooncogene_is_dysfunctional(current_cell) &&
                           !supressor_is_dysfunctional(current_cell)    &&
                           !apoptosis_is_dysfunctional(current_cell)))
                    {
                        erase_cell(current_cell_coord);

                        if(left_boundary_flag)
                            thread_locks[thread_index-1].unlock();
                        if(right_boundary_flag)
                            thread_locks[thread_index].unlock();

                        continue;
                    }

                    // 2. If apoptosis genes are functional and cell is not adjacent to the basal membrane
                    //    or a myoepithelial cell, mark the cell for death
                    // 3. If dead, remove from lattice
                    if(!apoptosis_is_dysfunctional(current_cell) &&
                       get_point_radius(current_cell_coord.x, current_cell_coord.y) < radius-2)
                    {
                        erase_cell(current_cell_coord);
                        if(left_boundary_flag)
                            thread_locks[thread_index-1].unlock();
                        if(right_boundary_flag)
                            thread_locks[thread_index].unlock();
                        continue;
                    }

                    // 4. Progenitor cells may reproduce into neighboring positions according to
                    //    progenitor hierarchy for three reasons:
                    int current_cell_type = get_cell_type(get_cell(current_cell_coord));
                    if (current_cell_type == PROG_STEM     ||
                        current_cell_type == PROG_BIPOTENT ||
                        current_cell_type == PROG_MYOEPI   ||
                        current_cell_type == PROG_LUMEPI)
                    {
                        // Cancerous reproduction—if the cell has a dysfunctional (50%) protooncogene
                        // and tumor suppressor gene, cell reproduces into any vacant neighboring lattice
                        // point or by pushing a neighboring cell to a vacant lattice point and reproducing
                        if(protooncogene_is_dysfunctional(current_cell) &&
                           supressor_is_dysfunctional(current_cell))
                        {
                            if (cancerous_reproduce(current_cell_coord));
                            {
                                reproduction_flag = true;
                                cancerous_reproduction_flag = true;
                            }
                        }
                        else
                        {
                            // Hormonal reproduction—if cell is stochastically chosen (1% chance), cell
                            // pushes a adjcance neighboring cell to a vacant lattice point and reproduces
                            // [ALT]: No permitimos rep. hormonal de MYOEPI para que no haya MYOEPIs en la capa LUMEPI
                            if(current_cell_type != PROG_MYOEPI && ThreadLocalRandom.current().nextInt(1, 101) == 1)
                            {
                                if (hormonal_reproduce(current_cell_coord))
                                {
                                    reproduction_flag = true;
                                }
                            }
                            else
                            {
                                // Normal reproduction—if an adjacent neighboring point is vacant and does
                                // not break the bi-layer ductal structure (Fig. 3C), the cell reproduces
                                if (reproduce(current_cell_coord))
                                {
                                    reproduction_flag = true;
                                    if (is_considered_cancerous(current_cell))
                                        cancerous_reproduction_flag = true;
                                }
                            }
                        }
                    }

                    // 5. Allow non-stem cells to migrate one point to a vacant neighboring lattice point
                    //    consistent with the bi-layer ductal structure (Fig. 3C)
                    if(current_cell_type != PROG_STEM)
                    {
                        migrate(current_cell_coord);
                    }

                    if(left_boundary_flag)
                        thread_locks[thread_index-1].unlock();
                    if(right_boundary_flag)
                        thread_locks[thread_index].unlock();

                    counter_lock.lock();
                    total_population++;
                    if(is_dysfunctional(current_cell))
                        dysfunctional_population++;
                    if (reproduction_flag)
                    {
                        total_reproductions++;
                        if (cancerous_reproduction_flag)
                            dysfunctional_reproductions++;
                    }
                    // Para Graph1
                    if(!DCIS_achieved && is_considered_cancerous(current_cell))
                    {
                        DCIS_achieved = true;
                        DCIS_generation = generation_number;
                    }
                    update_counters(current_cell_coord);
                    counter_lock.unlock();
                }
            }
        }
    }

    protected void next_generation()
    {
        if (initialized)
            normal_routine();
        else
        {
            if(!local_initialized)
                init_routine();
        }
    }
}