Guest User

Untitled

a guest
Jan 20th, 2018
90
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. package ductsim;
  2.  
  3. import java.util.ArrayList;
  4. import java.util.Arrays;
  5. import java.util.Collections;
  6. import java.util.List;
  7. import java.util.concurrent.BrokenBarrierException;
  8. import java.util.concurrent.CyclicBarrier;
  9. import java.util.concurrent.ExecutorService;
  10. import java.util.concurrent.Executors;
  11. import java.util.concurrent.ThreadLocalRandom;
  12. import java.util.concurrent.locks.ReentrantLock;
  13.  
  14. public class Duct implements Runnable {
  15.    
  16.     public class Coordinate
  17.     {
  18.         public int x;
  19.         public int y;
  20.         public int z;
  21.        
  22.         public Coordinate(int x, int y, int z)
  23.             { this.x = x; this.y = y; this.z = z; }
  24.         public boolean equals(Coordinate c)
  25.             { return (this.x == c.x) && (this.y == c.y) &&(this.z == c.z); }
  26.     }
  27.    
  28.     public static final int BASAL          = 1;
  29.     public static final int PROG_STEM      = 2;
  30.     public static final int PROG_BIPOTENT  = 3;
  31.     public static final int PROG_MYOEPI    = 4;
  32.     public static final int PROG_LUMEPI    = 5;
  33.     public static final int MYOEPI         = 6;
  34.     public static final int LUMEPI         = 7;
  35.    
  36.     public static final int gene_length = 32;
  37.    
  38.     private static int[][][] duct;
  39.    
  40.     private static int     length;
  41.     private static int     radius;
  42.     private static int     diameter;
  43.     private static boolean HGP;           // Indica presencia o no de HGP
  44.     private static int     capacity;      // Capacidad total del ducto
  45.     private static int     stem_capacity; // Capacidad de las dos capas principales del ducto
  46.  
  47.     private static final int mutation_rate = 5; // Representado como mutation_rate/1000
  48.  
  49.     private static int     generation_number;
  50.     private static boolean stem_initialized;  // Indica si se han colocado las células madre
  51.     private static boolean initialized;       // Indica si se ha inicializado el ducto
  52.     private boolean        local_initialized;
  53.  
  54.     private static int total_population;            // Contadores de población
  55.     private static int dysfunctional_population;
  56.     private static int total_reproductions;         // Contadores de reproducción
  57.     private static int dysfunctional_reproductions;
  58.     private static int stem_counter;                // Contadores de mutación
  59.     private static int bipotent_counter;
  60.     private static int proglumepi_counter;
  61.     private static int lumepi_counter;
  62.     private static int mut_stem_counter;
  63.     private static int mut_bipotent_counter;
  64.     private static int mut_proglumepi_counter;
  65.     private static int mut_lumepi_counter;
  66.    
  67.     private static boolean DCIS_achieved;   // Indica si se ha alcanzado DCIS o no
  68.     private static int     DCIS_generation; // Generación en que se generó DCIS
  69.  
  70.     private static int n_threads; // nº de hilos
  71.     private static ExecutorService threadPool;
  72.     private static CyclicBarrier   threadBarrier;
  73.     private static Duct[]          threadArray;
  74.     private static ReentrantLock[] thread_locks;
  75.     private static final ReentrantLock   counter_lock = new ReentrantLock();
  76.  
  77.     private static int total_generations;
  78.     private Integer[]  index_array;
  79.     private int        thread_index;
  80.     private Coordinate boundaries;
  81.    
  82.    
  83.     /**
  84.      * Constructor del ducto
  85.      *
  86.      * @param l Longitud del ducto (número de secciones)
  87.      * @param r Radio del ducto
  88.      * @param h Presencia o no de HGP
  89.      * @param n Número de hilos
  90.      */
  91.     public Duct(int l, int r, boolean h, int n)
  92.     {
  93.         if (threadPool != null)
  94.         {
  95.             threadPool.shutdown();
  96.             while(!threadPool.isTerminated()){}
  97.         }  
  98.  
  99.         threadPool    = null;
  100.         threadBarrier = null;
  101.         threadArray   = null;
  102.         n_threads     = n;
  103.         duct = new int[l][r*2][r*2];
  104.         length   = l;
  105.         radius   = r;
  106.         HGP      = h;
  107.         diameter = r*2;
  108.         capacity = 0;
  109.         stem_capacity = 0;
  110.        
  111.         generation_number = 0;
  112.         stem_initialized  = false;
  113.         initialized       = false;
  114.    
  115.         DCIS_achieved   = false;
  116.         DCIS_generation = 0;
  117.         index_array = new Integer[length];
  118.        
  119.         for(int z=0 ; z<length ; z++)
  120.         {
  121.             index_array[z] = z;
  122.             for(int x=0 ; x<=(r/2) ; x++)
  123.             {
  124.                 for(int y=x ; y<=r ; y++)
  125.                 {
  126.                     if (get_point_radius(x,y) == r)
  127.                     {
  128.                         generate_cells_8way(z, y, x, BASAL);
  129.                     }
  130.                 }
  131.             }
  132.         }
  133.  
  134.         for(int i=0 ; i<diameter ; i++)
  135.         {
  136.             for(int j=0 ; j<diameter ; j++)
  137.             {
  138.                 if(get_point_radius(i,j) < radius)
  139.                     capacity++;
  140.                 if(get_point_radius(i,j) == radius-1 || get_point_radius(i,j) == radius-2)
  141.                     stem_capacity++;
  142.             }
  143.         }
  144.         capacity *= length;
  145.         stem_capacity *= length;
  146.         threadBarrier = new CyclicBarrier(n_threads+1); // Hilos que lanzamos mas éste, el principal
  147.         threadArray = new Duct[n_threads];
  148.         thread_locks  = new ReentrantLock[n_threads-1];
  149.        
  150.         int step = length/n_threads;
  151.         for(int i=0 ; i<n_threads ; i++)
  152.         {
  153.             int begin =     i*step;
  154.             int end   = (i+1)*step-1;
  155.             if ((i+1) == n_threads)
  156.                 end = length-1;
  157.            
  158.             threadArray[i] = new Duct(begin, end, i);
  159.             if (i < n_threads-1)
  160.                 thread_locks[i]  = new ReentrantLock();
  161.         }
  162.     }
  163.  
  164.     /* Constructor para crear hilos (begin y end inclusive en rango)*/
  165.     public Duct(int begin, int end, int index)
  166.     {
  167.         thread_index = index;
  168.         index_array = new Integer[end-begin+1];
  169.         // Aquí Coordinate se usa como un par de números, no como una coordenada
  170.         boundaries = new Coordinate(begin, end, 0);
  171.         // [begin, end] ambos inclusive
  172.         for(int z=begin, i=0 ; z<=end ; z++, i++)
  173.         {
  174.             index_array[i] = z;
  175.         }
  176.     }
  177.    
  178.     // [OPT] Cambiar las llamadas a función por acceso directo a la matriz
  179.     // duct.[c.z][c.y][c.x] o algo así
  180.     // en realidad para lo que se usa da igual
  181.     public int  get_cell (Coordinate c)        { return duct[c.z][c.y][c.x]; }
  182.     public void set_cell(Coordinate co, int c) { this.duct[co.z][co.y][co.x] = c; }
  183.     public void erase_cell(Coordinate co)      { this.duct[co.z][co.y][co.x] = 0; }
  184.    
  185.     public void generate_cells_8way(int z, int y, int x, int type)
  186.     {
  187.         int r = radius;
  188.         this.duct[z][x][y]         = Cell(type);
  189.         this.duct[z][x][2*r-y]     = Cell(type);
  190.         this.duct[z][2*r-x][y]     = Cell(type);
  191.         this.duct[z][2*r-x][2*r-y] = Cell(type);
  192.         // simetría diagonal
  193.         this.duct[z][y][x]         = Cell(type);
  194.         this.duct[z][y][2*r-x]     = Cell(type);
  195.         this.duct[z][2*r-y][x]     = Cell(type);
  196.         this.duct[z][2*r-y][2*r-x] = Cell(type);
  197.        
  198.     }
  199.    
  200.     public static int get_point_radius(int x, int y)
  201.     {
  202.         // [ALT]: Fórmula que se ajuste al círculo que se ve en el paper
  203.         return (int)(Math.sqrt((x-radius)*(x-radius)+(y-radius)*(y-radius))+1);
  204.     }
  205.  
  206.     public int Cell (int type)
  207.     {
  208.         return type * (int)1e8;
  209.     }
  210.    
  211.     public int get_cell_type(int c)
  212.     {
  213.         return c/(int)1e8;
  214.     }
  215.  
  216.     public int get_gene_housekeeping(int c)
  217.     {
  218.         return (c % (int)1e8) / (int)1e6;
  219.     }
  220.  
  221.     public int get_gene_protooncogene(int c)
  222.     {
  223.         return (c % (int)1e6) / (int)1e4;
  224.     }
  225.  
  226.     public int get_gene_supressor(int c)
  227.     {
  228.         return (c % (int)1e4) / (int)1e2;
  229.     }
  230.  
  231.     public int get_gene_apoptosis(int c)
  232.     {
  233.         return c % (int)1e2;
  234.     }
  235.    
  236.     public int set_cell_type(int c, int type)
  237.     {
  238.         return (c % (int)1e8) + type * (int)1e8;
  239.     }
  240.    
  241.     // Muta 10% del genoma
  242.     public int mutate_HGP(int c)
  243.     {
  244.         int counter = (int)((((float)gene_length) / 100f) * 10f);
  245.         return c += counter * (int)1e6 + counter * (int)1e4 + counter * (int)1e2 + counter;
  246.     }
  247.  
  248.     // Mutation Rate representado como por mil, no por cien
  249.     public int mutate(int c, int mutation_rate)
  250.     {
  251.         int hk_c = gene_length - get_gene_housekeeping(c);
  252.         int po_c = gene_length - get_gene_protooncogene(c);
  253.         int sp_c = gene_length - get_gene_supressor(c);
  254.         int ap_c = gene_length - get_gene_apoptosis(c);
  255.        
  256.         for(int i=0 ; i < hk_c ; i++)
  257.         {
  258.             if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
  259.                 c += (int)1e6;
  260.         }
  261.         for(int i=0 ; i < po_c ; i++)
  262.         {
  263.             if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
  264.                 c += (int)1e4;
  265.         }
  266.         for(int i=0 ; i < sp_c ; i++)
  267.         {
  268.             if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
  269.                 c += (int)1e2;
  270.         }
  271.         for(int i=0 ; i < ap_c ; i++)
  272.         {
  273.             if(ThreadLocalRandom.current().nextInt(1, 1001) <= mutation_rate)
  274.                 c += 1;
  275.         }
  276.        
  277.         return c;
  278.     }
  279.    
  280.     // Método que usamos para ver si la célula es cancerígena según el paper
  281.     public boolean is_considered_cancerous(int c)
  282.     {
  283.         return //get_gene_housekeeping(c)() >= 16 ||
  284.                get_gene_protooncogene(c) >= 16 ||
  285.                //get_gene_supressor(c)() >= 16 ||
  286.                get_gene_apoptosis(c)     >= 16;
  287.     }
  288.    
  289.     public int get_mutations(int c)
  290.     {
  291.         return get_gene_housekeeping(c)  +
  292.                get_gene_protooncogene(c) +
  293.                get_gene_supressor(c)     +
  294.                get_gene_apoptosis(c);
  295.     }
  296.    
  297.     public boolean is_dysfunctional(int c)
  298.     {
  299.         return get_gene_housekeeping(c)  >= 16 ||
  300.                get_gene_protooncogene(c) >= 16 ||
  301.                get_gene_supressor(c)     >= 16 ||
  302.                get_gene_apoptosis(c)     >= 16;
  303.     }
  304.  
  305.     public boolean housekeeping_is_dysfunctional(int c)
  306.     {
  307.         return get_gene_housekeeping(c) >= 16;
  308.     }
  309.    
  310.     public boolean protooncogene_is_dysfunctional(int c)
  311.     {
  312.         return get_gene_protooncogene(c) >= 16;
  313.     }
  314.        
  315.     public boolean supressor_is_dysfunctional(int c)
  316.     {
  317.         return get_gene_supressor(c)>= 16;
  318.     }
  319.        
  320.     public boolean apoptosis_is_dysfunctional(int c)
  321.     {
  322.         return get_gene_apoptosis(c) >= 16;
  323.     }
  324.  
  325.     public List<Float> mutations_stats()
  326.     {
  327.         List<Float> result = new ArrayList<>();
  328.         result.add((float)mut_stem_counter/(float)stem_counter);
  329.         result.add((float)mut_bipotent_counter/(float)bipotent_counter);
  330.         result.add((float)mut_proglumepi_counter/(float)proglumepi_counter);
  331.         result.add((float)mut_lumepi_counter/(float)lumepi_counter);
  332.         return result;
  333.     }
  334.  
  335.     public void update_counters(Coordinate c)
  336.     {  
  337.         switch(get_cell_type(get_cell(c)))
  338.         {
  339.             case PROG_STEM:
  340.                 stem_counter++;
  341.                 mut_stem_counter += get_mutations(get_cell(c));
  342.                 break;
  343.             case PROG_BIPOTENT:
  344.                 bipotent_counter++;
  345.                 mut_bipotent_counter += get_mutations(get_cell(c));
  346.                 break;
  347.             case PROG_LUMEPI:
  348.                 proglumepi_counter++;
  349.                 mut_proglumepi_counter += get_mutations(get_cell(c));
  350.                 break;
  351.             case LUMEPI:
  352.                 lumepi_counter++;
  353.                 mut_lumepi_counter += get_mutations(get_cell(c));
  354.                 break;
  355.         }
  356.     }
  357.  
  358.     protected void execute(int nGens)
  359.     {
  360.         if(generation_number == 0)
  361.         {
  362.             init_routine_stem(); // Hacemos la primera generación (colocación de stem cells) secuencialmente
  363.         }
  364.        
  365.         total_generations = nGens;
  366.        
  367.         try {
  368.         threadPool = Executors.newFixedThreadPool(n_threads);
  369.         for (int i=0; i<threadArray.length; i++)
  370.             threadPool.execute(threadArray[i]);
  371.         } catch(Exception e){System.out.println(e);}
  372.        
  373.         for (int i=0; i<total_generations; i++)
  374.         {
  375.             total_population = 0;
  376.             dysfunctional_population = 0;
  377.             total_reproductions = 0;
  378.             dysfunctional_reproductions = 0;
  379.  
  380.             stem_counter = 0;
  381.             bipotent_counter = 0;
  382.             proglumepi_counter = 0;
  383.             lumepi_counter = 0;
  384.             mut_stem_counter = 0;
  385.             mut_bipotent_counter = 0;
  386.             mut_proglumepi_counter = 0;
  387.             mut_lumepi_counter = 0;
  388.            
  389.             try
  390.             {  
  391.                 threadBarrier.await();
  392.                 generation_number++;
  393.                
  394.                 if(!initialized)
  395.                 {
  396.                     initialized = threadArray[0].get_local_initialized();
  397.                     for (int t=1; t<threadArray.length; t++)
  398.                         initialized = initialized && threadArray[t].get_local_initialized();
  399.                 }
  400.             }      
  401.             catch (BrokenBarrierException | InterruptedException e){}
  402.         }
  403.         threadPool.shutdown();
  404.         while(!threadPool.isTerminated()){}
  405.     }
  406.    
  407.     public void run()
  408.     {
  409.         for (int i=0; i<total_generations; i++)
  410.         {
  411.             try
  412.             {
  413.                 next_generation();
  414.                 threadBarrier.await();
  415.             }
  416.             catch (BrokenBarrierException | InterruptedException e){}
  417.         }
  418.     }
  419.  
  420.     public int get_generation_number()           { return generation_number; }
  421.     public int get_length()                      { return length; }
  422.     public int get_radius()                      { return radius; }
  423.     public int get_diameter()                    { return diameter; }
  424.     public int get_capacity()                    { return capacity; }
  425.     public boolean get_local_initialized()       { return local_initialized; }
  426.    
  427.     public int get_dysfunctional_population()    { return dysfunctional_population; }
  428.     public int get_total_population()            { return total_population; }
  429.     public int get_dysfunctional_reproductions() { return dysfunctional_reproductions; }
  430.     public int get_total_reproductions()         { return total_reproductions; }
  431.     public int[][] get_slice(int z)              { return duct[z]; }
  432.    
  433.     public int DCIS_achieved()
  434.     {
  435.         if(!DCIS_achieved)
  436.             return 0;
  437.         else
  438.             return DCIS_generation;
  439.     }
  440.  
  441.     /**
  442.      * @param c Coordenada
  443.      * @return true si la coordenada c se encuentra dentro del array ducto
  444.      *         (no necesariamente dentro de la membrana basal)
  445.      *         false en caso contrario
  446.      */
  447.     public boolean coordinate_is_inbounds(Coordinate c)
  448.     {
  449.         return (c.x >= 0) && (c.x < diameter) &&
  450.                (c.y >= 0) && (c.y < diameter) &&
  451.                (c.z >= 0) && (c.z < length);
  452.     }
  453.    
  454.     /**
  455.      * @param c Coordenada
  456.      * @return ArrayList de Coordenadas con todos los vecinos (vecindad Moore 3D) de la cordenada c
  457.      */
  458.     private List<Coordinate> get_neighbors(Coordinate c)
  459.     {
  460.         List<Coordinate> neighbors = new ArrayList<>();
  461.         int x = c.x;
  462.         int y = c.y;
  463.         int z = c.z;
  464.         // Las 6 coordenadas "principales" (vecindad directa)
  465.         neighbors.add(new Coordinate(x,   y,   z-1));
  466.         neighbors.add(new Coordinate(x,   y,   z+1));
  467.         neighbors.add(new Coordinate(x,   y-1, z));
  468.         neighbors.add(new Coordinate(x-1, y,   z));
  469.         neighbors.add(new Coordinate(x+1, y,   z));
  470.         neighbors.add(new Coordinate(x,   y+1, z));
  471.         // Vecinos diagonales
  472.         neighbors.add(new Coordinate(x-1, y-1, z-1)); // diagonal
  473.         neighbors.add(new Coordinate(x,   y-1, z-1)); // diagonal
  474.         neighbors.add(new Coordinate(x+1, y-1, z-1)); // diagonal
  475.         neighbors.add(new Coordinate(x-1, y,   z-1)); // diagonal
  476.         neighbors.add(new Coordinate(x+1, y,   z-1)); // diagonal
  477.         neighbors.add(new Coordinate(x-1, y+1, z-1)); // diagonal
  478.         neighbors.add(new Coordinate(x,   y+1, z-1)); // diagonal
  479.         neighbors.add(new Coordinate(x+1, y+1, z-1)); // diagonal
  480.         neighbors.add(new Coordinate(x-1, y-1, z+1)); // diagonal
  481.         neighbors.add(new Coordinate(x,   y-1, z+1)); // diagonal
  482.         neighbors.add(new Coordinate(x+1, y-1, z+1)); // diagonal
  483.         neighbors.add(new Coordinate(x-1, y,   z+1)); // diagonal
  484.         neighbors.add(new Coordinate(x+1, y,   z+1)); // diagonal
  485.         neighbors.add(new Coordinate(x-1, y+1, z+1)); // diagonal
  486.         neighbors.add(new Coordinate(x,   y+1, z+1)); // diagonal
  487.         neighbors.add(new Coordinate(x+1, y+1, z+1)); // diagonal
  488.         neighbors.add(new Coordinate(x-1, y-1, z));   // diagonal      
  489.         neighbors.add(new Coordinate(x+1, y-1, z));   // diagonal      
  490.         neighbors.add(new Coordinate(x-1, y+1, z));   // diagonal  
  491.         neighbors.add(new Coordinate(x+1, y+1, z));   // diagonal
  492.  
  493.         return neighbors;
  494.     }
  495.    
  496.     /**
  497.      * @param cell_coord
  498.      * @return Devuelve una coordenada aleatoria que cumple las siguientes condiciones:
  499.      *          - Es vecina de cell_coord
  500.      *          - Se encuentra dentro de la membrana basal
  501.      *          - Está vacía (no está ocupada por una célula)
  502.      *         Si tal coordenada no existe, devuelve null
  503.      */
  504.     private Coordinate get_vacant_neighbor(Coordinate cell_coord)
  505.     {
  506.         List<Coordinate> neighbors = get_neighbors(cell_coord);
  507.         Collections.shuffle(neighbors);
  508.         for(int i=0 ; i < neighbors.size() ; i++)
  509.         {
  510.             Coordinate aux_coord = neighbors.get(i);
  511.             if(coordinate_is_inbounds(aux_coord) == true &&
  512.                get_cell(aux_coord) == 0 &&
  513.                get_point_radius(aux_coord.x, aux_coord.y) < radius)
  514.             {
  515.                 return neighbors.get(i);
  516.             }
  517.         }
  518.         return null;
  519.     }
  520.  
  521.     /**
  522.      * @param cell_coord Coordenada a la que se va a migrar
  523.      * @param old_coord  Coordenada desde la que se migra
  524.      * @return true si cell_coord tiene alguna célula vecina aparte de old_coord
  525.      *         false en caso contrario
  526.      */
  527.     private boolean has_adjadcent_neighbor(Coordinate cell_coord, Coordinate old_coord)
  528.     {
  529.         List<Coordinate> neighbors = get_neighbors(cell_coord);
  530.         for(int i=0 ; i < neighbors.size() ; i++)
  531.         {
  532.             Coordinate aux_coord = neighbors.get(i);
  533.             if(coordinate_is_inbounds(aux_coord) == true &&
  534.                get_cell(aux_coord) != 0 &&
  535.                !aux_coord.equals(old_coord))
  536.             {
  537.                 return true;
  538.             }
  539.         }
  540.         return false;
  541.     }
  542.    
  543.     private boolean reproduce(Coordinate current_cell_coord)
  544.     {  
  545.         List<Coordinate> neighbors = get_neighbors(current_cell_coord);
  546.         Collections.shuffle(neighbors);
  547.         int current_cell_type = get_cell_type(get_cell(current_cell_coord));
  548.        
  549.         /*
  550.         INIT:   Allow progenitor cells to reproduce into vacant neighboring lattice points
  551.                 without cells according to progenitor hierarchy
  552.         NORMAL: if an adjacent neighboring point is vacant and does not break the bi-layer
  553.                 ductal structure (Fig. 3C), the cell reproduces
  554.         */
  555.  
  556.         /* [ALT]: (a) Permitir que una célula mioepitelial se reproduzca en el sitio que sea, y luego
  557.                       se corrija ésto en la migración, o bien
  558.                   (b) Permitir que sólo se genere una célula en un sitio que le 'corresponde' [actual]
  559.         */
  560.  
  561.         int max_radius = radius-1;
  562.         int min_radius = radius-2;
  563.        
  564.         if(current_cell_type == PROG_LUMEPI)
  565.             max_radius--;
  566.         if(current_cell_type == PROG_MYOEPI)
  567.             min_radius++;
  568.        
  569.         Coordinate new_cell_coord = null;
  570.         boolean end_flag = false;
  571.         for(int i=0 ; i<neighbors.size() ; i++)
  572.         {
  573.             // Condiciones: la célula está inbounds, vacía, y entre min_radius y max_radius (inclusive)
  574.             new_cell_coord = neighbors.get(i);
  575.             if(coordinate_is_inbounds(new_cell_coord) == true &&
  576.                get_cell(new_cell_coord) == 0 &&
  577.                // (b)
  578.                get_point_radius(new_cell_coord.x, new_cell_coord.y) <= max_radius &&
  579.                get_point_radius(new_cell_coord.x, new_cell_coord.y) >= min_radius)
  580.                // (a)
  581.                //get_point_radius(new_cell_coord.x, new_cell_coord.y) <= radius &&
  582.                //get_point_radius(new_cell_coord.x, new_cell_coord.y) >= radius-2)                 
  583.             {
  584.                 end_flag = true;
  585.                 break;
  586.             }
  587.         }
  588.        
  589.         if(end_flag == false)
  590.         {
  591.             return false;
  592.         }
  593.  
  594.         int parent = get_cell(current_cell_coord);
  595.         int new_cell = parent;
  596.         switch (get_cell_type(parent))
  597.         {
  598.             // [OPT] Creo que los switch se compilan como una cadena de if-else;
  599.             // por lo que es mejor poner la opción más común primero
  600.             case PROG_LUMEPI:
  601.                 new_cell = set_cell_type(new_cell, LUMEPI);
  602.                 break;
  603.             case PROG_BIPOTENT:
  604.                 // [ALT] Decidimos qué tipo de célula va a generar una bipotente en base a
  605.                 // la posición en que se va a generar
  606.                 if (get_point_radius(new_cell_coord.x, new_cell_coord.y) == radius-1)
  607.                     new_cell = set_cell_type(new_cell, PROG_MYOEPI);
  608.                 else
  609.                     new_cell = set_cell_type(new_cell, PROG_LUMEPI);
  610.                 break;
  611.             case PROG_STEM:
  612.                 new_cell = set_cell_type(new_cell, PROG_BIPOTENT);
  613.                 break;
  614.             case PROG_MYOEPI:
  615.                 new_cell = set_cell_type(new_cell, MYOEPI);
  616.                 break;
  617.             default:
  618.                 break;
  619.         }
  620.         // Mutamos tanto madre como hija
  621.         parent = mutate(parent, mutation_rate);
  622.         new_cell = mutate(new_cell, mutation_rate);
  623.         set_cell(current_cell_coord, parent);
  624.         set_cell(new_cell_coord, new_cell);
  625.        
  626.         return true;
  627.     }
  628.    
  629.     private boolean cancerous_reproduce(Coordinate current_cell_coord)
  630.     {  
  631.         Coordinate new_cell_coord = null;
  632.         Coordinate pushed_position = null;
  633.  
  634.         List<Coordinate> neighbors = get_neighbors(current_cell_coord);
  635.         Collections.shuffle(neighbors);
  636.         boolean end_flag = false;
  637.         for (int i=0 ; i<neighbors.size() ; i++)
  638.         {
  639.             pushed_position = null;
  640.             new_cell_coord = neighbors.get(i);
  641.  
  642.             if (coordinate_is_inbounds(new_cell_coord) == true)
  643.             {
  644.                 if (get_cell(new_cell_coord) == 0 &&
  645.                     get_point_radius(new_cell_coord.x, new_cell_coord.y) < radius)
  646.                 {
  647.                     end_flag = true;
  648.                     break;
  649.                 }
  650.                 // [ALT] Sólo empujamos células LUMEPI (?)
  651.                 if (get_point_radius(new_cell_coord.x, new_cell_coord.y) < radius &&
  652.                     (pushed_position = get_vacant_neighbor(new_cell_coord)) != null)
  653.                 {
  654.                     end_flag = true;
  655.                     break;
  656.                 }
  657.             }  
  658.         }
  659.  
  660.         if(end_flag == false)
  661.         {
  662.             return false;
  663.         }
  664.        
  665.         if(pushed_position != null)
  666.         {
  667.             int pushed_cell = get_cell(new_cell_coord);
  668.             set_cell(pushed_position, pushed_cell);
  669.             erase_cell(new_cell_coord);
  670.         }
  671.        
  672.         int parent = get_cell(current_cell_coord);
  673.         int new_cell = parent;
  674.         switch (get_cell_type(parent))
  675.         {
  676.             case PROG_LUMEPI:
  677.                 new_cell = set_cell_type(new_cell, LUMEPI);
  678.                 break;
  679.             case PROG_BIPOTENT:
  680.                 if (get_point_radius(new_cell_coord.x, new_cell_coord.y) == radius-1)
  681.                     new_cell = set_cell_type(new_cell, PROG_MYOEPI);
  682.                 else
  683.                     new_cell = set_cell_type(new_cell, PROG_LUMEPI);
  684.                 break;
  685.             case PROG_STEM:
  686.                 new_cell = set_cell_type(new_cell, PROG_BIPOTENT);
  687.                 break;
  688.             case PROG_MYOEPI:
  689.                 new_cell = set_cell_type(new_cell, MYOEPI);
  690.                 break;
  691.  
  692.             default:
  693.                 break;
  694.         }
  695.         parent =   mutate(parent, mutation_rate);
  696.         new_cell = mutate(new_cell, mutation_rate);
  697.         set_cell(current_cell_coord, parent);
  698.         set_cell(new_cell_coord, new_cell);
  699.        
  700.         return true;
  701.     }
  702.    
  703.     private boolean hormonal_reproduce(Coordinate current_cell_coord)
  704.     {
  705.         List<Coordinate> neighbors = get_neighbors(current_cell_coord);
  706.         Collections.shuffle(neighbors);
  707.        
  708.         int current_cell_type = get_cell_type(get_cell(current_cell_coord));
  709.         Coordinate new_cell_coord = null;
  710.        
  711.         boolean end_flag = false;
  712.         Coordinate pushed_position = null;
  713.         for (int i=0 ; i<neighbors.size() ; i++)
  714.         {
  715.             new_cell_coord = neighbors.get(i);
  716.             if (coordinate_is_inbounds(new_cell_coord) != false &&
  717.                 get_point_radius(new_cell_coord.x, new_cell_coord.y) < radius &&
  718.                 get_cell(new_cell_coord) != 0 &&
  719.                 (pushed_position = get_vacant_neighbor(new_cell_coord)) != null
  720.                )
  721.             {
  722.                 end_flag = true;
  723.                 break;
  724.             }
  725.         }
  726.  
  727.         if(end_flag == false)
  728.         {
  729.             return false;
  730.         }
  731.  
  732.         int pushed_cell = get_cell(new_cell_coord);
  733.         set_cell(pushed_position, pushed_cell);
  734.         erase_cell(new_cell_coord);
  735.        
  736.         int parent = get_cell(current_cell_coord);
  737.         int new_cell = parent;
  738.         switch (get_cell_type(parent))
  739.         {
  740.             case PROG_LUMEPI:
  741.                 new_cell = set_cell_type(new_cell, LUMEPI);
  742.                 break;
  743.             case PROG_BIPOTENT:
  744.                 if (get_point_radius(new_cell_coord.x, new_cell_coord.y) == radius-1)
  745.                     new_cell = set_cell_type(new_cell, PROG_MYOEPI);
  746.                 else
  747.                     new_cell = set_cell_type(new_cell, PROG_LUMEPI);
  748.                 break;
  749.             case PROG_STEM:
  750.                 new_cell = set_cell_type(new_cell, PROG_BIPOTENT);
  751.                 break;
  752.             case PROG_MYOEPI:
  753.                 new_cell = set_cell_type(new_cell, MYOEPI);
  754.                 break;
  755.             default:
  756.                 break;
  757.         }
  758.         parent   = mutate(parent, mutation_rate);
  759.         new_cell = mutate(new_cell, mutation_rate);
  760.         set_cell(current_cell_coord, parent);
  761.         set_cell(new_cell_coord, new_cell);
  762.        
  763.         return true;
  764.     }
  765.    
  766.     private boolean migrate(Coordinate current_cell_coord)
  767.     {
  768.         List<Coordinate> neighbors = get_neighbors(current_cell_coord);
  769.         Collections.shuffle(neighbors);
  770.        
  771.         int current_cell = get_cell(current_cell_coord);
  772.         int current_cell_type = get_cell_type(current_cell);
  773.        
  774.         int max_radius = radius-1;
  775.         // [ALT]: No permitimos que una célula migre hacia "dentro"; siempre hacia un punto en el mismo radio
  776.         // que el punto en el que está, o mayor
  777.         int min_radius = get_point_radius(current_cell_coord.x, current_cell_coord.y);
  778.        
  779.         /*
  780.         if(get_cell(current_cell_coord).apoptosis_is_dysfunctional())
  781.         {
  782.             // WARN: Arreglar la migración de células cancerígenas (asegurar que migran a un lugar conectado
  783.             // con el nucleo principal de poblacion)
  784.             //min_radius = 0;
  785.         }
  786.         */
  787.         if(current_cell_type == PROG_LUMEPI || current_cell_type == LUMEPI)
  788.             max_radius = radius-2;
  789.             //if(current_cell_type == PROG_MYOEPI || current_cell_type == MYOEPI)
  790.                 //min_radius = radius-1;
  791.         Coordinate new_cell_coord = null;
  792.         boolean end_flag = false;
  793.         for (int i=0 ; i<neighbors.size() ; i++)
  794.         {
  795.             new_cell_coord = neighbors.get(i);
  796.             if (coordinate_is_inbounds(new_cell_coord) == true &&
  797.                 get_cell(new_cell_coord) == 0 &&
  798.                 get_point_radius(new_cell_coord.x, new_cell_coord.y) <= max_radius &&
  799.                 get_point_radius(new_cell_coord.x, new_cell_coord.y) >= min_radius &&
  800.                 has_adjadcent_neighbor(new_cell_coord, current_cell_coord))
  801.             {
  802.                 end_flag = true;
  803.                 break;
  804.             }
  805.        
  806.         }
  807.        
  808.         if(end_flag == false)
  809.         {
  810.             return false;
  811.         }
  812.        
  813.         erase_cell(current_cell_coord);
  814.         set_cell(new_cell_coord, current_cell);
  815.        
  816.         return true;
  817.     }
  818.  
  819.     private void init_routine_stem()
  820.     {
  821.         if(radius == 1)
  822.             return;
  823.        
  824.         // 1. Draw the BASAL membrane by designating outer boundary points in a cylindrical orientation
  825.        
  826.         /* Hecho en constructor de SliceSec */
  827.        
  828.         // 2. Place STEM CELLS within the basal membrane
  829.         //    (less than 5% of the population - 200 for a 24,000 cell pop. => [ ~ 0.84% ])
  830.        
  831.         // [ALT]: Con nuestro modelo de radio, las células que forman las dos capas iniciales son 20,800 contra las
  832.         //        24,000 del paper.
  833.         float stem_cells_percentage = 1f;  // Porcentaje del total de células que queremos que sean células madre
  834.         int n_stem_cells = (int)((((float)stem_capacity) / 100) * stem_cells_percentage);
  835.        
  836.         if(n_stem_cells == 0)
  837.             n_stem_cells = 1;
  838.        
  839.         int chosen_x, chosen_y, chosen_z;
  840.         // Elegimos posición de células madre por polling de puntos aleatorios
  841.         for(int counter = 1 ; counter <= n_stem_cells ; counter++)
  842.         {
  843.             // Elegir SliceSec aleatorio
  844.             chosen_z = ThreadLocalRandom.current().nextInt(0, length);
  845.             // Recorrer slice, elegir punto aleatorio hasta encontrar uno de radio r-1
  846.             do {
  847.                 chosen_x = ThreadLocalRandom.current().nextInt(0, diameter);
  848.                 chosen_y = ThreadLocalRandom.current().nextInt(0, diameter);
  849.             } while (get_point_radius(chosen_x, chosen_y) != radius-1 ||
  850.                      get_cell(new Coordinate(chosen_x, chosen_y, chosen_z)) != 0);
  851.            
  852.             // Establecer nueva célula basal
  853.             Coordinate new_coordinate = new Coordinate(chosen_x, chosen_y, chosen_z);
  854.             int new_stem_cell = Cell(PROG_STEM);
  855.            
  856.             // 3. If HGP, mutate all genes in all stem cells by 10%
  857.             // (ver método en Cell.java)
  858.             if(HGP)
  859.             {
  860.                 // [ALT]: Mutación determinista (mutar 10% del genoma) o estocástica (mutar con 10% de probabilidad) ?
  861.                 new_stem_cell = mutate_HGP(new_stem_cell);
  862.             }
  863.             set_cell(new_coordinate, new_stem_cell);
  864.         }
  865.         total_population += n_stem_cells;
  866.         stem_initialized = true;
  867.     }
  868.    
  869.     private void init_routine()
  870.     {
  871.         if(!stem_initialized)
  872.         {
  873.             init_routine_stem();
  874.             return;
  875.         }
  876.        
  877.         boolean end_flag = true; // Indica si la inicialización del ducto ha acabado
  878.        
  879.         Integer[] index_array_aux = new Integer[index_array.length];
  880.         System.arraycopy(index_array, 0, index_array_aux, 0, index_array.length);
  881.         Collections.shuffle(Arrays.asList(index_array_aux));
  882.         // Iteramos sobre todas las coordenadas posibles (abierto a optimizaciones)
  883.         for(int z_i = 0 ; z_i < index_array_aux.length ; z_i++) {
  884.             int z = index_array_aux[z_i];
  885.             for(int y=0 ; y<diameter ; y++) {
  886.                 for(int x=0 ; x<diameter ; x++) {
  887.                     Coordinate current_cell_coord = new Coordinate(x,y,z);
  888.  
  889.                     // Si la coordenada no contiene ninguna célula o una célula basal, pasamos a la siguiente
  890.                     if(get_cell(current_cell_coord) == 0 || get_cell_type(get_cell(current_cell_coord)) == BASAL)
  891.                         continue;
  892.  
  893.                     total_population++;
  894.                     if(is_dysfunctional(get_cell(current_cell_coord)))
  895.                         dysfunctional_population++;
  896.  
  897.                     int current_cell_type = get_cell_type(get_cell(current_cell_coord));
  898.                     // 4. Allow PROGENITOR cells to REPRODUCE into VACANT NEIGHBORING lattice points WITHOUT CELLS
  899.                     //    according to progenitor hierarchy
  900.                     if (current_cell_type == PROG_STEM     ||
  901.                         current_cell_type == PROG_BIPOTENT ||
  902.                         current_cell_type == PROG_MYOEPI   ||
  903.                         current_cell_type == PROG_LUMEPI)
  904.                     {
  905.                         if (reproduce(current_cell_coord) == true)
  906.                         {
  907.                             end_flag = false;
  908.                             total_reproductions++;
  909.                         }
  910.                     }
  911.                     update_counters(current_cell_coord);
  912.                     // 5. Allow NON-STEM cells to MIGRATE 1 point to a VACANT NEIGHBORING lattice point consistent
  913.                     //    with the bi-layer ductal structure
  914.                     if(current_cell_type != PROG_STEM)
  915.                     {
  916.                         if (migrate(current_cell_coord) == true)
  917.                             end_flag = false;
  918.                     }
  919.                 }
  920.             }
  921.         }  
  922.         local_initialized = end_flag;
  923.         if(local_initialized)
  924.             System.out.println(Thread.currentThread().getName() + " inicializado en " + generation_number);
  925.     }
  926.    
  927.     private void normal_routine()
  928.     {
  929.         for(int z_i = 0 ; z_i < index_array.length ; z_i++) {
  930.             int z = index_array[z_i];
  931.             for(int y=0 ; y<diameter ; y++) {
  932.                 for(int x=0 ; x<diameter ; x++) {
  933.        
  934.                     Coordinate current_cell_coord = new Coordinate(x,y,z);
  935.                     //System.out.println(Thread.currentThread().getName() + " accede a " + x + "," + y + "," + z + "[gen" + generation_number + "]");
  936.                     boolean left_boundary_flag  = false;
  937.                     boolean right_boundary_flag = false;
  938.                     boolean reproduction_flag = false;
  939.                     boolean cancerous_reproduction_flag = false;
  940.                    
  941.                     // Barrera izquierda
  942.                     if(boundaries.x != 0 && z <= boundaries.x+1)
  943.                     {
  944.                         thread_locks[thread_index-1].lock();
  945.                         left_boundary_flag = true;
  946.                     }
  947.                     // Barrera derecha
  948.                     if(boundaries.y != length-1 && z >= boundaries.y-1)
  949.                     {
  950.                         thread_locks[thread_index].lock();
  951.                         right_boundary_flag = true;
  952.                     }
  953.                    
  954.                     if(get_cell(current_cell_coord) == 0 || get_cell_type(get_cell(current_cell_coord)) == BASAL)
  955.                     {
  956.                         if(left_boundary_flag)
  957.                             thread_locks[thread_index-1].unlock();
  958.                         if(right_boundary_flag)
  959.                             thread_locks[thread_index].unlock();
  960.                         continue;
  961.                     }
  962.                    
  963.                     // Si la célula se encuentra en zona de conflicto, debe tomar el cerrojo de la partición contigua
  964.                     // Cuando esto suceda, para evitar deadlock, la partición más a la izquierda tiene prioridad sobre el
  965.                     // cerrojo de la partición a su derecha.
  966.                    
  967.                     // Partición derecha quiere acceder a partición izquierda: partición derecha no tiene prioridad.
  968.                     // Ceder cerrojo, y luego adquirir los dos cerrojos, tras lo que se re-comproará la condición inicial
  969.                    
  970.                    
  971.                     int current_cell = get_cell(current_cell_coord);
  972.  
  973.                     // 1. If housekeeping or apoptosis genes are dysfunctional (50% mutated), mark the cell for death
  974.                     // 3. If dead, remove from lattice
  975.                     if(housekeeping_is_dysfunctional(current_cell) ||
  976.                           (protooncogene_is_dysfunctional(current_cell) &&
  977.                            !supressor_is_dysfunctional(current_cell)    &&
  978.                            !apoptosis_is_dysfunctional(current_cell)))
  979.                     {
  980.                         erase_cell(current_cell_coord);
  981.                        
  982.                         if(left_boundary_flag)
  983.                             thread_locks[thread_index-1].unlock();
  984.                         if(right_boundary_flag)
  985.                             thread_locks[thread_index].unlock();
  986.                            
  987.                         continue;
  988.                     }
  989.  
  990.                     // 2. If apoptosis genes are functional and cell is not adjacent to the basal membrane
  991.                     //    or a myoepithelial cell, mark the cell for death
  992.                     // 3. If dead, remove from lattice
  993.                     if(!apoptosis_is_dysfunctional(current_cell) &&
  994.                        get_point_radius(current_cell_coord.x, current_cell_coord.y) < radius-2)
  995.                     {
  996.                         erase_cell(current_cell_coord);
  997.                         if(left_boundary_flag)
  998.                             thread_locks[thread_index-1].unlock();
  999.                         if(right_boundary_flag)
  1000.                             thread_locks[thread_index].unlock();
  1001.                         continue;
  1002.                     }
  1003.  
  1004.                     // 4. Progenitor cells may reproduce into neighboring positions according to
  1005.                     //    progenitor hierarchy for three reasons:
  1006.                     int current_cell_type = get_cell_type(get_cell(current_cell_coord));
  1007.                     if (current_cell_type == PROG_STEM     ||
  1008.                         current_cell_type == PROG_BIPOTENT ||
  1009.                         current_cell_type == PROG_MYOEPI   ||
  1010.                         current_cell_type == PROG_LUMEPI)
  1011.                     {
  1012.                         // Cancerous reproduction—if the cell has a dysfunctional (50%) protooncogene
  1013.                         // and tumor suppressor gene, cell reproduces into any vacant neighboring lattice
  1014.                         // point or by pushing a neighboring cell to a vacant lattice point and reproducing
  1015.                         if(protooncogene_is_dysfunctional(current_cell) &&
  1016.                            supressor_is_dysfunctional(current_cell))
  1017.                         {
  1018.                             if (cancerous_reproduce(current_cell_coord));
  1019.                             {
  1020.                                 reproduction_flag = true;
  1021.                                 cancerous_reproduction_flag = true;
  1022.                             }
  1023.                         }
  1024.                         else
  1025.                         {
  1026.                             // Hormonal reproduction—if cell is stochastically chosen (1% chance), cell
  1027.                             // pushes a adjcance neighboring cell to a vacant lattice point and reproduces
  1028.                             // [ALT]: No permitimos rep. hormonal de MYOEPI para que no haya MYOEPIs en la capa LUMEPI
  1029.                             if(current_cell_type != PROG_MYOEPI && ThreadLocalRandom.current().nextInt(1, 101) == 1)
  1030.                             {
  1031.                                 if (hormonal_reproduce(current_cell_coord))
  1032.                                 {
  1033.                                     reproduction_flag = true;
  1034.                                 }
  1035.                             }
  1036.                             else
  1037.                             {
  1038.                                 // Normal reproduction—if an adjacent neighboring point is vacant and does
  1039.                                 // not break the bi-layer ductal structure (Fig. 3C), the cell reproduces
  1040.                                 if (reproduce(current_cell_coord))
  1041.                                 {
  1042.                                     reproduction_flag = true;
  1043.                                     if (is_considered_cancerous(current_cell))
  1044.                                         cancerous_reproduction_flag = true;
  1045.                                 }
  1046.                             }
  1047.                         }
  1048.                     }
  1049.                    
  1050.                     // 5. Allow non-stem cells to migrate one point to a vacant neighboring lattice point
  1051.                     //    consistent with the bi-layer ductal structure (Fig. 3C)
  1052.                     if(current_cell_type != PROG_STEM)
  1053.                     {
  1054.                         migrate(current_cell_coord);
  1055.                     }
  1056.                    
  1057.                     if(left_boundary_flag)
  1058.                         thread_locks[thread_index-1].unlock();
  1059.                     if(right_boundary_flag)
  1060.                         thread_locks[thread_index].unlock();
  1061.                    
  1062.                     counter_lock.lock();
  1063.                     total_population++;
  1064.                     if(is_dysfunctional(current_cell))
  1065.                         dysfunctional_population++;
  1066.                     if (reproduction_flag)
  1067.                     {
  1068.                         total_reproductions++;
  1069.                         if (cancerous_reproduction_flag)
  1070.                             dysfunctional_reproductions++;
  1071.                     }
  1072.                     // Para Graph1
  1073.                     if(!DCIS_achieved && is_considered_cancerous(current_cell))
  1074.                     {
  1075.                         DCIS_achieved = true;
  1076.                         DCIS_generation = generation_number;
  1077.                     }
  1078.                     update_counters(current_cell_coord);
  1079.                     counter_lock.unlock();
  1080.                 }
  1081.             }
  1082.         }
  1083.     }
  1084.    
  1085.     protected void next_generation()
  1086.     {  
  1087.         if (initialized)
  1088.             normal_routine();
  1089.         else
  1090.         {
  1091.             if(!local_initialized)
  1092.                 init_routine();
  1093.         }
  1094.     }
  1095. }
RAW Paste Data