# HW 3 sim_droughtstat.R JLCY, 19 Dec 2013

############################################################ 

# simulation statistics

{

    # Set the threshold level
    frac = 1
    th = frac * mean(array(t(LFmonflow[1:nyrs, 1:12])))

    # drought statistics for simulation

    # Longest Surplus
    A = array(t(simdismon))
    xx = A
    Y = A
    xx[xx <= th] = 0
    xx[xx > th] = 1  # surplus
    n = length(xx)
    x1 = xx
    if (xx[1] == 1) 
        x1 = c(0, x1)
    if (xx[1] == 1) 
        Y = c(0, Y)
    if (xx[n] == 1) 
        x1 = c(x1, 0)
    if (xx[n] == 1) 
        Y = c(Y, 0)
    n = length(x1[x1 == 1])
    xy = order(x1)[1:(length(x1) - n)]
    y3 = xy
    x3 = diff(xy) - 1
    x3 = x3[x3 > 0]
    # mxsp[k1] = max(x3) #Longest Surplus LS
    sim_LS[k] = max(x3)

    # Maximum surplus
    nsp = length(x3)
    surp = 1:nsp
    nsp = length(y3)

    isp = 0
    for (i in 1:(nsp - 1)) {
        i1 = y3[i] + 1
        i2 = y3[i + 1] - 1
        if (i2 >= i1) {
            isp = isp + 1
            surp[isp] = sum(Y[i1:i2])
        }
    }

    # print(c(k1,isp,nsp,length(x3))) print(c(k1,mxsp[k1],mxdef[k1],i1,i2,nsp))

    # maxs[k1] = max(surp) #Maximum Surplus MS
    sim_MS[k] = max(surp)

    # for longest drought
    xx = A
    Y = A
    xx[xx <= th] = 1  # drought
    xx[xx > th] = 0
    n = length(xx)
    x1 = xx
    if (xx[1] == 1) 
        x1 = c(0, x1)
    if (xx[1] == 1) 
        Y = c(0, Y)
    if (xx[n] == 1) 
        x1 = c(x1, 0)
    if (xx[n] == 1) 
        Y = c(Y, 0)
    n = length(x1[x1 == 1])
    xy = order(x1)[1:(length(x1) - n)]
    y3 = xy
    x2 = diff(xy) - 1
    x2 = x2[x2 > 0]
    # mxdef[k1] = max(x2) #Longest Drought LD
    sim_LD[k] = max(x2)

    # Maximum Deficit
    nsp = length(x2)
    surp = 1:nsp
    nsp = length(y3)
    isp = 0
    for (i in 1:(nsp - 1)) {
        i1 = y3[i] + 1
        i2 = y3[i + 1] - 1
        if (i2 >= i1) {
            isp = isp + 1
            surp[isp] = sum(Y[i1:i2])
        }
    }
    # print(c(k1,isp,nsp,length(x3)))
    sim_MD[k] = max(surp)  # Maximum Drought      MD      

}

## Error: 找不到物件 'LFmonflow'