Programming Assignment 1: Percolation

The Monte Carlo simulation method is used to estimate the seepage threshold.

Percolation.A random distribution of insulation and metal composite system is provided. For example, we want to know which parts must be metallic materials to make this composite system an electrical conductor. Or in a porous terrain, there is water or oil on the surface, where water or oil can penetrate from the surface to the bottom layer. Scientists call this process a model called percolation.

The model.In assignment, a nxn lattice is used to represent the percolation system. Each lattice is open or closed, and the white is closed and black. If a grid is full, it must first be the opening amount, and then it indicates that it can penetrate to this position from the top through the opened grid connected in 4 directions. When a system is percolates, it indicates that it can penetrate from the top layer to the bottom layer. That is to say, the bottom layer is full.

 

The problem.The researchers are interested in the problem. If each lattice is independent and the probability of being opened is P, what is the probability of the system percolates? P = 0, the probability of percolates is 0, P = 100, and the probability of percolates is 100. Is the probability P Distribution of 20 x20 and x100 grids:

When N is large enough, there is a threshold p *, so that when P <p *, any n * n mesh can hardly be penetrated, and when P> P *, it can basically be infiltrated. P * has no accurate numerical solutions. The task is to write an algorithm for calculating the estimated p.

Percolation data type.

public class Percolation {   public Percolation(int N)              // create N-by-N grid, with all sites blocked   public void open(int i, int j)         // open site (row i, column j) if it is not already   public boolean isOpen(int i, int j)    // is site (row i, column j) open?   public boolean isFull(int i, int j)    // is site (row i, column j) full?   public boolean percolates()            // does the system percolate?}

The upper left corner of the grid is (1, 1), and the lower right corner is (N, N ).

Use Boolean [] matrix to mark the opening and closing of the grid.

Isopen indicates whether the current grid is open. You can directly return the value of the current position of the matrix.

Think about how open, isfull, and percolates are implemented and complex.

What should open do? The idea is simple. You can open a closed grid and connect it with the opened grid in four directions.

Isfull determines whether the current vertex can penetrate. The simplest idea is to enumerate each opened grid in the first layer, and check whether the set query is consistent with the current grid in one set, if so, it is successful. The worst case of complexity enumeration is N, plus the need to judge one find () each time.

Percolates identifies whether the system can be infiltrated. The simplest idea is to enumerate the grids opened at the first and last layers and determine whether the system is in a collection. If the penetration succeeds, complexity n ^ 2, plus and query set find () to judge.

The simplest idea is too violent and complicated, so we begin to optimize it.

1. add a virtual node and two virtual nodes at the beginning and end. The first virtual node is associated with all open elements at the first layer. The last virtual node is associated with all open nodes at the last layer, for the isfull operation, if the first virtual node and the current node are in the same set, the operation succeeds. For the percolates operation, if the first and last virtual nodes are in the same set, the system penetration is successful. In this way, there will be no enumerated for loops in all these operations, but there will also be a problem, backwash, because some grids cannot penetrate from above, however, after connecting to the end virtual node, they are also in a set with the first virtual node.

2. How to Avoid backwash? I thought about a comparison method and added another query set. Now I have two query sets. One is only responsible for maintaining the first virtual node. When I perform isfull to determine whether to only consider this and check the set, maintenance of the first and last virtual nodes in another query set is used in percolates's judgment. This solves the backwash problem, but sacrifices the space.

The following is the implementation code:

public class Percolation {        private boolean[] matrix;    private int row, col;    private WeightedQuickUnionUF wquUF;    private WeightedQuickUnionUF wquUFTop;    private boolean alreadyPercolates;        public Percolation(int N) {        if (N < 1) throw new IllegalArgumentException("Illeagal Argument");        wquUF = new WeightedQuickUnionUF(N*N+2);        wquUFTop = new WeightedQuickUnionUF(N*N+1);        alreadyPercolates = false;        row = N;        col = N;        matrix = new boolean[N*N+1];    }        private void validate(int i, int j) {        if (i < 1 || i > row)             throw new IndexOutOfBoundsException("row index i out of bounds");        if (j < 1 || j > col)             throw new IndexOutOfBoundsException("col index j out of bounds");            }        public void open(int i, int j) {        validate(i, j);        int curIdx = (i-1)*col + j;         matrix[curIdx] = true;        if (i == 1) {            wquUF.union(curIdx, 0);            wquUFTop.union(curIdx, 0);        }        if (i == row) {            wquUF.union(curIdx, row*col+1);        }                int[] dx = {1, -1, 0, 0};        int[] dy = {0, 0, 1, -1};        for (int dir = 0; dir < 4; dir++) {            int posX = i + dx[dir];            int posY = j + dy[dir];            if (posX <= row && posX >= 1                     && posY <= row && posY >= 1                     && isOpen(posX, posY)) {                wquUF.union(curIdx, (posX-1)*col+posY);                wquUFTop.union(curIdx, (posX-1)*col+posY);            }        }    }        public boolean isOpen(int i, int j) {        validate(i, j);        return matrix[(i-1)*col + j];    }        public boolean isFull(int i, int j) {        validate(i, j);        int curIdx = (i-1)*col+j;        if (wquUFTop.find(curIdx) == wquUFTop.find(0)) return true;        return false;    }        public boolean percolates() {        if (alreadyPercolates) return true;        if (wquUF.find(0) == wquUF.find(row*col+1)) {            alreadyPercolates = true;            return true;        }         return false;    }        public static void main(String[] args) {        Percolation perc = new Percolation(2);        perc.open(1, 1);        perc.open(1, 2);        perc.open(2, 1);        System.out.println(perc.percolates());    }    }

Monte Carlo simulation.It is estimated that the threshold value of percolation is disabled during initialization. The system randomly finds a closed position to open until the system can penetrate. The total number of opened grids is the threshold value.

This part is a simple part of numerical calculation.

public class PercolationStats {    private double[] x;    private int expTimes;        public PercolationStats(int N, int T) {        if (N <= 0 || T <= 0)             throw new IllegalArgumentException("Illeagal Argument");        x = new double[T+1];        expTimes = T;        for (int i = 1; i <= T; ++i) {            Percolation perc = new Percolation(N);            boolean[] isEmptySiteLine = new boolean[N+1];            int numOfLine = 0;            while (true) {                    int posX, posY;                do {                    posX = StdRandom.uniform(N)+1;                    posY = StdRandom.uniform(N)+1;                } while (perc.isOpen(posX, posY));                perc.open(posX, posY);                x[i] += 1;                if (!isEmptySiteLine[posX]) {                    isEmptySiteLine[posX] = true;                    numOfLine++;                }                if (numOfLine == N) {                    if (perc.percolates())                        break;                }            }            x[i] = x[i] / (double) (N*N);        }    }        public double mean() {        double mu = 0.0;        for (int i = 1; i <= expTimes; ++i) {            mu += x[i];        }        return mu / (double) expTimes;    }        public double stddev() {        if (expTimes == 1)            return Double.NaN;        double sigma = 0.0;        double mu = mean();        for (int i = 1; i <= expTimes; ++i) {            sigma += (x[i]-mu)*(x[i]-mu);        }        return Math.sqrt(sigma / (double) (expTimes-1));    }        public double confidenceLo() {        double mu = mean();        double sigma = stddev();        return mu - 1.96*sigma / Math.sqrt(expTimes);    }        public double confidenceHi() {        double mu = mean();        double sigma = stddev();        return mu + 1.96*sigma / Math.sqrt(expTimes);    }        public static void main(String[] args) {        int N = Integer.parseInt(args[0]);        int T = Integer.parseInt(args[1]);        PercolationStats percStats = new PercolationStats(N, T);        StdOut.printf("mean = %f\n", percStats.mean());        StdOut.printf("stddev = %f\n", percStats.stddev());        StdOut.printf("95%% confidence interval = %f, %f\n",                       percStats.confidenceLo(), percStats.confidenceHi());            }        }