text:
# ------------------------------------
# calcucation of the viscosity term by central difference
# ------------------------------------
function central_diff(E_vis_hat, F_vis_hat, G_vis_hat, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
                    QconU, QconD, QconL, QconR, QconF, QconB, cellxmax, cellymax, cellzmax, mu, mut, mut_bd, lambda,
                    vecAx, vecAy, vecAz, specific_heat_ratio, volume, Rd, nval, yplus, swith_wall, icell)

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax+1 -icell
                rho_av = 0.5*(QbaseR[i,j,k,1] + QbaseL[i,j,k,1])
                u_av   = 0.5*(QbaseR[i,j,k,2] + QbaseL[i,j,k,2])
                v_av   = 0.5*(QbaseR[i,j,k,3] + QbaseL[i,j,k,3])
                w_av   = 0.5*(QbaseR[i,j,k,4] + QbaseL[i,j,k,4])

                mu_av     = 0.5*(mu[i-1,j,k] + mu[i,j,k])
                lambda_av = (0.5*(1/lambda[i-1,j,k] + 1/lambda[i,j,k]))^(-1)

                volume_av = 0.5*(volume[i,j,k] + volume[i-1,j,k])

                dudxi = QbaseR[i,j,k,2] - QbaseL[i,j,k,2]
                dvdxi = QbaseR[i,j,k,3] - QbaseL[i,j,k,3]
                dwdxi = QbaseR[i,j,k,4] - QbaseL[i,j,k,4]
                dTdxi = QbaseR[i,j,k,5]/(QbaseR[i,j,k,1]*Rd) - QbaseL[i,j,k,5]/(QbaseL[i,j,k,1]*Rd)

                dudeta = 0.25 * (QbaseU[i,  j+1,k,2] - QbaseD[i,  j+1,k,2] + QbaseU[i,  j,k,2] - QbaseD[i,  j,k,2]
                               + QbaseU[i-1,j+1,k,2] - QbaseD[i-1,j+1,k,2] + QbaseU[i-1,j,k,2] - QbaseD[i-1,j,k,2])
                dvdeta = 0.25 * (QbaseU[i,  j+1,k,3] - QbaseD[i,  j+1,k,3] + QbaseU[i,  j,k,3] - QbaseD[i,  j,k,3]
                               + QbaseU[i-1,j+1,k,3] - QbaseD[i-1,j+1,k,3] + QbaseU[i-1,j,k,3] - QbaseD[i-1,j,k,3])
                dwdeta = 0.25 * (QbaseU[i,  j+1,k,4] - QbaseD[i,  j+1,k,4] + QbaseU[i,  j,k,4] - QbaseD[i,  j,k,4]
                               + QbaseU[i-1,j+1,k,4] - QbaseD[i-1,j+1,k,4] + QbaseU[i-1,j,k,4] - QbaseD[i-1,j,k,4])

                dT1 = QbaseU[  i,j+1,k,5]/(QbaseU[  i,j+1,k,1]*Rd) - QbaseD[  i,j+1,k,5]/(QbaseD[  i,j+1,k,1]*Rd)
                dT2 = QbaseU[  i,  j,k,5]/(QbaseU[  i,  j,k,1]*Rd) - QbaseD[  i,  j,k,5]/(QbaseD[  i,  j,k,1]*Rd)
                dT3 = QbaseU[i-1,j+1,k,5]/(QbaseU[i-1,j+1,k,1]*Rd) - QbaseD[i-1,j+1,k,5]/(QbaseD[i-1,j+1,k,1]*Rd)
                dT4 = QbaseU[i-1,  j,k,5]/(QbaseU[i-1,  j,k,1]*Rd) - QbaseD[i-1,  j,k,5]/(QbaseD[i-1,  j,k,1]*Rd)
                dTdeta = 0.25 * (dT1 + dT2 + dT3 + dT4)

                dudzeta = 0.25 * (QbaseF[i,  j,k+1,2] - QbaseB[i,  j,k+1,2] + QbaseF[i,  j,k,2] - QbaseB[i,  j,k,2]
                                + QbaseF[i-1,j,k+1,2] - QbaseB[i-1,j,k+1,2] + QbaseF[i-1,j,k,2] - QbaseB[i-1,j,k,2])
                dvdzeta = 0.25 * (QbaseF[i,  j,k+1,3] - QbaseB[i,  j,k+1,3] + QbaseF[i,  j,k,3] - QbaseB[i,  j,k,3]
                                + QbaseF[i-1,j,k+1,3] - QbaseB[i-1,j,k+1,3] + QbaseF[i-1,j,k,3] - QbaseB[i-1,j,k,3])
                dwdzeta = 0.25 * (QbaseF[i,  j,k+1,4] - QbaseB[i,  j,k+1,4] + QbaseF[i,  j,k,4] - QbaseB[i,  j,k,4]
                                + QbaseF[i-1,j,k+1,4] - QbaseB[i-1,j,k+1,4] + QbaseF[i-1,j,k,4] - QbaseB[i-1,j,k,4])

                dT1 = QbaseF[  i,  j,k+1,5]/(QbaseF[  i,  j,k+1,1]*Rd) - QbaseB[  i,  j,k+1,5]/(QbaseB[  i,  j,k+1,1]*Rd)
                dT2 = QbaseF[  i,  j,  k,5]/(QbaseF[  i,  j,  k,1]*Rd) - QbaseB[  i,  j,  k,5]/(QbaseB[  i,  j,  k,1]*Rd)
                dT3 = QbaseF[i-1,  j,k+1,5]/(QbaseF[i-1,  j,k+1,1]*Rd) - QbaseB[i-1,  j,k+1,5]/(QbaseB[i-1,  j,k+1,1]*Rd)
                dT4 = QbaseF[i-1,  j,  k,5]/(QbaseF[i-1,  j,  k,1]*Rd) - QbaseB[i-1,  j,  k,5]/(QbaseB[i-1,  j,  k,1]*Rd)
                dTdzeta = 0.25 * (dT1 + dT2 + dT3 + dT4)

                vecAy_xav    = 0.25*( vecAy[i,j,k,1] + vecAy[i,j+1,k,1] + vecAy[i-1,j,k,1] + vecAy[i-1,j+1,k,1] )
                vecAy_yav    = 0.25*( vecAy[i,j,k,2] + vecAy[i,j+1,k,2] + vecAy[i-1,j,k,2] + vecAy[i-1,j+1,k,2] )
                vecAy_zav    = 0.25*( vecAy[i,j,k,3] + vecAy[i,j+1,k,3] + vecAy[i-1,j,k,3] + vecAy[i-1,j+1,k,3] )

                vecAz_xav    = 0.25*( vecAz[i,j,k,1] + vecAz[i,j,k+1,1] + vecAz[i-1,j,k,1] + vecAz[i-1,j,k+1,1] )
                vecAz_yav    = 0.25*( vecAz[i,j,k,2] + vecAz[i,j,k+1,2] + vecAz[i-1,j,k,2] + vecAz[i-1,j,k+1,2] )
                vecAz_zav    = 0.25*( vecAz[i,j,k,3] + vecAz[i,j,k+1,3] + vecAz[i-1,j,k,3] + vecAz[i-1,j,k+1,3] )

                dudx = vecAx[i,j,k,1] * dudxi + vecAy_xav * dudeta + vecAz_xav * dudzeta
                dudy = vecAx[i,j,k,2] * dudxi + vecAy_yav * dudeta + vecAz_yav * dudzeta
                dudz = vecAx[i,j,k,3] * dudxi + vecAy_zav * dudeta + vecAz_zav * dudzeta

                dvdx = vecAx[i,j,k,1] * dvdxi + vecAy_xav * dvdeta + vecAz_xav * dvdzeta
                dvdy = vecAx[i,j,k,2] * dvdxi + vecAy_yav * dvdeta + vecAz_yav * dvdzeta
                dvdz = vecAx[i,j,k,3] * dvdxi + vecAy_zav * dvdeta + vecAz_zav * dvdzeta

                dwdx = vecAx[i,j,k,1] * dwdxi + vecAy_xav * dwdeta + vecAz_xav * dwdzeta
                dwdy = vecAx[i,j,k,2] * dwdxi + vecAy_yav * dwdeta + vecAz_yav * dwdzeta
                dwdz = vecAx[i,j,k,3] * dwdxi + vecAy_zav * dwdeta + vecAz_zav * dwdzeta

                sigma_xx = 2/3*mu_av*( 2*dudx -   dvdy -   dwdz )
                sigma_yy = 2/3*mu_av*( - dudx + 2*dvdy -   dwdz )
                sigma_zz = 2/3*mu_av*( - dudx -   dvdy + 2*dwdz )

                sigma_xy = mu_av*( dudy + dvdx )
                sigma_xz = mu_av*( dvdz + dwdy )
                sigma_yz = mu_av*( dwdx + dudz )

                dTdx = vecAx[i,j,k,1]*dTdxi + vecAy_xav*dTdeta + vecAz_xav*dTdzeta
                dTdy = vecAx[i,j,k,2]*dTdxi + vecAy_yav*dTdeta + vecAz_yav*dTdzeta
                dTdz = vecAx[i,j,k,3]*dTdxi + vecAy_zav*dTdeta + vecAz_zav*dTdzeta

                betax = u_av*sigma_xx + v_av*sigma_xy + w_av*sigma_xz + lambda_av * dTdx
                betay = u_av*sigma_xy + v_av*sigma_yy + w_av*sigma_yz + lambda_av * dTdy
                betaz = u_av*sigma_xz + v_av*sigma_yz + w_av*sigma_zz + lambda_av * dTdz

                yplus_av = 0.5 *(yplus[i-1,j,k] + yplus[i,j,k])

                if swith_wall[1] == 1 && i == 1+icell
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,1] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
                elseif swith_wall[2] == 1 && i == cellxmax - icell
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,1] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
                else
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,1] =
                    Smagorinsky_model(dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, dTdx, dTdy, dTdz,
                                    rho_av, u_av, v_av, w_av, volume_av, yplus_av, i, j)

                    #tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    #e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,1] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

                #=
                    if i == 51 && j ==3
                        println(" tur tau ")
                        println(tau_xx)
                        println(tau_xy)
                        println(tau_yy)
                        println(e_sgs_x)
                        println(e_sgs_y)
                        println(" kiyo ")
                        println(vecAx[i,j,1]*sigma_xx + vecAx[i,j,2]*sigma_xy)
                        println(vecAx[i,j,1]*tau_xx + vecAx[i,j,2]*tau_xy)
                        println(vecAx[i,j,1]*tau_xx)
                        println(vecAx[i,j,2]*tau_xy)
                    end
                    =#
                end
                #=
                if i == 20 && j == 20
                    println("(20,20,du)")
                    println(dudxi)
                    println(dudeta)
                    println(vecAx[i,j,1])
                    println(vecAy_xav)
                    println(volume_av)
                    println((vecAx[i,j,1] * dudxi + vecAy_xav * dudeta)/volume_av)
                    println(tau_xx)
                    println(tau_xy)
                    println(tau_yy)
                end
                if i == 150 && j == 150
                    println("(150,150,du)")
                    println(dudxi)
                    println(dudeta)
                    println(vecAx[i,j,1])
                    println(vecAy_xav)
                    println(volume_av)
                    println((vecAx[i,j,1] * dudxi + vecAy_xav * dudeta)/volume_av)
                    println(tau_xx)
                    println(tau_xy)
                    println(tau_yy)
                end
                =#

                E_vis_hat[i,j,k,1] = 0.0
                E_vis_hat[i,j,k,2] = ((vecAx[i,j,k,1]*sigma_xx + vecAx[i,j,k,2]*sigma_xy + vecAx[i,j,k,3]*sigma_xz) -
                                      (vecAx[i,j,k,1]*tau_xx   + vecAx[i,j,k,2]*tau_xy   + vecAx[i,j,k,3]*tau_xz)) / volume_av
                E_vis_hat[i,j,k,3] = ((vecAx[i,j,k,1]*sigma_xy + vecAx[i,j,k,2]*sigma_yy + vecAx[i,j,k,3]*sigma_yz) -
                                      (vecAx[i,j,k,1]*tau_xy   + vecAx[i,j,k,2]*tau_yy   + vecAx[i,j,k,3]*tau_yz)) / volume_av
                E_vis_hat[i,j,k,4] = ((vecAx[i,j,k,1]*sigma_xz + vecAx[i,j,k,2]*sigma_yz + vecAx[i,j,k,3]*sigma_zz) -
                                      (vecAx[i,j,k,1]*tau_xz   + vecAx[i,j,k,2]*tau_yz   + vecAx[i,j,k,3]*tau_zz)) / volume_av
                E_vis_hat[i,j,k,5] = ((vecAx[i,j,k,1]*betax    + vecAx[i,j,k,2]*betay    + vecAx[i,j,k,3]*betaz) -
                                      (vecAx[i,j,k,1]*e_sgs_x  + vecAx[i,j,k,2]*e_sgs_y  + vecAx[i,j,k,3]*e_sgs_z))/ volume_av
            end
        end
    end

    for k in 1+icell:cellzmax -icell
        for j in 1+icell:cellymax+1 -icell
            for i in 1+icell:cellxmax -icell
                rho_av = 0.5*(QbaseU[i,j,k,1] + QbaseD[i,j,k,1])
                u_av   = 0.5*(QbaseU[i,j,k,2] + QbaseD[i,j,k,2])
                v_av   = 0.5*(QbaseU[i,j,k,3] + QbaseD[i,j,k,3])
                w_av   = 0.5*(QbaseU[i,j,k,4] + QbaseD[i,j,k,4])

                mu_av     = 0.5*(mu[i,j-1,k] + mu[i,j,k])
                lambda_av = (0.5*(1/lambda[i,j-1,k] + 1/lambda[i,j,k]))^(-1)

                volume_av = 0.5*(volume[i,j,k] + volume[i,j-1,k])

                dudxi = 0.25 * (QbaseR[i+1,  j,k,2] - QbaseL[i+1,  j,k,2] + QbaseR[i,  j,k,2] - QbaseL[i,  j,k,2]
                              + QbaseR[i+1,j-1,k,2] - QbaseL[i+1,j-1,k,2] + QbaseR[i,j-1,k,2] - QbaseL[i,j-1,k,2])
                dvdxi = 0.25 * (QbaseR[i+1,  j,k,3] - QbaseL[i+1,  j,k,3] + QbaseR[i,  j,k,3] - QbaseL[i,  j,k,3]
                              + QbaseR[i+1,j-1,k,3] - QbaseL[i+1,j-1,k,3] + QbaseR[i,j-1,k,3] - QbaseL[i,j-1,k,3])
                dwdxi = 0.25 * (QbaseR[i+1,  j,k,4] - QbaseL[i+1,  j,k,4] + QbaseR[i,  j,k,4] - QbaseL[i,  j,k,4]
                              + QbaseR[i+1,j-1,k,4] - QbaseL[i+1,j-1,k,4] + QbaseR[i,j-1,k,4] - QbaseL[i,j-1,k,4])

                dT1 = QbaseR[i+1,  j,k,5]/(QbaseR[i+1,  j,k,1]*Rd) - QbaseL[i+1,  j,k,5]/(QbaseL[i+1,  j,k,1]*Rd)
                dT2 = QbaseR[  i,  j,k,5]/(QbaseR[  i,  j,k,1]*Rd) - QbaseL[  i,  j,k,5]/(QbaseL[  i,  j,k,1]*Rd)
                dT3 = QbaseR[i+1,j-1,k,5]/(QbaseR[i+1,j-1,k,1]*Rd) - QbaseL[i+1,j-1,k,5]/(QbaseL[i+1,j-1,k,1]*Rd)
                dT4 = QbaseR[  i,j-1,k,5]/(QbaseR[  i,j-1,k,1]*Rd) - QbaseL[  i,j-1,k,5]/(QbaseL[  i,j-1,k,1]*Rd)
                dTdxi = 0.25 * (dT1 + dT2 + dT3 + dT4)

                dudeta = QbaseU[i,j,k,2] - QbaseD[i,j,k,2]
                dvdeta = QbaseU[i,j,k,3] - QbaseD[i,j,k,3]
                dwdeta = QbaseU[i,j,k,4] - QbaseD[i,j,k,4]
                dTdeta = (QbaseU[i,j,k,5]/(QbaseU[i,j,k,1]*Rd) - QbaseD[i,j,k,5]/(QbaseD[i,j,k,1]*Rd))

                dudzeta = 0.25 * (QbaseF[i,  j,k+1,2] - QbaseB[i,  j,k+1,2] + QbaseF[i,  j,k,2] - QbaseB[i,  j,k,2]
                                + QbaseF[i,j-1,k+1,2] - QbaseB[i,j-1,k+1,2] + QbaseF[i,j-1,k,2] - QbaseB[i,j-1,k,2])
                dvdzeta = 0.25 * (QbaseF[i,  j,k+1,3] - QbaseB[i,  j,k+1,3] + QbaseF[i,  j,k,3] - QbaseB[i,  j,k,3]
                                + QbaseF[i,j-1,k+1,3] - QbaseB[i,j-1,k+1,3] + QbaseF[i,j-1,k,3] - QbaseB[i,j-1,k,3])
                dwdzeta = 0.25 * (QbaseF[i,  j,k+1,4] - QbaseB[i,  j,k+1,4] + QbaseF[i,  j,k,4] - QbaseB[i,  j,k,4]
                                + QbaseF[i,j-1,k+1,4] - QbaseB[i,j-1,k+1,4] + QbaseF[i,j-1,k,4] - QbaseB[i,j-1,k,4])

                dT1 = QbaseF[  i,  j,k+1,5]/(QbaseF[  i,  j,k+1,1]*Rd) - QbaseB[  i,  j,k+1,5]/(QbaseB[  i,  j,k+1,1]*Rd)
                dT2 = QbaseF[  i,  j,  k,5]/(QbaseF[  i,  j,  k,1]*Rd) - QbaseB[  i,  j,  k,5]/(QbaseB[  i,  j,  k,1]*Rd)
                dT3 = QbaseF[  i,j-1,k+1,5]/(QbaseF[  i,j-1,k+1,1]*Rd) - QbaseB[  i,j-1,k+1,5]/(QbaseB[  i,j-1,k+1,1]*Rd)
                dT4 = QbaseF[  i,j-1,  k,5]/(QbaseF[  i,j-1,  k,1]*Rd) - QbaseB[  i,j-1,  k,5]/(QbaseB[  i,j-1,  k,1]*Rd)
                dTdzeta = 0.25 * (dT1 + dT2 + dT3 + dT4)

                vecAx_xav    = 0.25*( vecAx[i,j,k,1] + vecAx[i+1,j,k,1] + vecAx[i,j-1,k,1] + vecAx[i+1,j-1,k,1] )
                vecAx_yav    = 0.25*( vecAx[i,j,k,2] + vecAx[i+1,j,k,2] + vecAx[i,j-1,k,2] + vecAx[i+1,j-1,k,2] )
                vecAx_zav    = 0.25*( vecAx[i,j,k,3] + vecAx[i+1,j,k,3] + vecAx[i,j-1,k,3] + vecAx[i+1,j-1,k,3] )

                vecAz_xav    = 0.25*( vecAz[i,j,k,1] + vecAz[i,j,k+1,1] + vecAz[i,j-1,k,1] + vecAz[i,j-1,k+1,1] )
                vecAz_yav    = 0.25*( vecAz[i,j,k,2] + vecAz[i,j,k+1,2] + vecAz[i,j-1,k,2] + vecAz[i,j-1,k+1,2] )
                vecAz_zav    = 0.25*( vecAz[i,j,k,3] + vecAz[i,j,k+1,3] + vecAz[i,j-1,k,3] + vecAz[i,j-1,k+1,3] )

                dudx = vecAx_xav * dudxi + vecAy[i,j,k,1] * dudeta + vecAz_xav * dudzeta
                dudy = vecAx_yav * dudxi + vecAy[i,j,k,2] * dudeta + vecAz_yav * dudzeta
                dudz = vecAx_zav * dudxi + vecAy[i,j,k,3] * dudeta + vecAz_zav * dudzeta

                dvdx = vecAx_xav * dvdxi + vecAy[i,j,k,1] * dvdeta + vecAz_xav * dvdzeta
                dvdy = vecAx_yav * dvdxi + vecAy[i,j,k,2] * dvdeta + vecAz_yav * dvdzeta
                dvdz = vecAx_zav * dvdxi + vecAy[i,j,k,3] * dvdeta + vecAz_zav * dvdzeta

                dwdx = vecAx_xav * dwdxi + vecAy[i,j,k,1] * dwdeta + vecAz_xav * dwdzeta
                dwdy = vecAx_yav * dwdxi + vecAy[i,j,k,2] * dwdeta + vecAz_yav * dwdzeta
                dwdz = vecAx_zav * dwdxi + vecAy[i,j,k,3] * dwdeta + vecAz_zav * dwdzeta

                sigma_xx = 2/3*mu_av*( 2*dudx -   dvdy -   dwdz )
                sigma_yy = 2/3*mu_av*( - dudx + 2*dvdy -   dwdz )
                sigma_zz = 2/3*mu_av*( - dudx -   dvdy + 2*dwdz )

                sigma_xy = mu_av*( dudy + dvdx )
                sigma_xz = mu_av*( dvdz + dwdy )
                sigma_yz = mu_av*( dwdx + dudz )

                dTdx = vecAx_xav*dTdxi + vecAy[i,j,k,1]*dTdeta + vecAz_xav*dTdzeta
                dTdy = vecAx_yav*dTdxi + vecAy[i,j,k,2]*dTdeta + vecAz_yav*dTdzeta
                dTdz = vecAx_zav*dTdxi + vecAy[i,j,k,3]*dTdeta + vecAz_zav*dTdzeta

                betax = u_av*sigma_xx + v_av*sigma_xy + w_av*sigma_xz + lambda_av * dTdx
                betay = u_av*sigma_xy + v_av*sigma_yy + w_av*sigma_yz + lambda_av * dTdy
                betaz = u_av*sigma_xz + v_av*sigma_yz + w_av*sigma_zz + lambda_av * dTdz

                yplus_av = 0.5 *(yplus[i,j-1,k] + yplus[i,j,k])

                if swith_wall[3] == 1 && j == 1+icell
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,2] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
                elseif swith_wall[3] == 1 && j == cellymax - icell
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,2] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
                else
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,2] =
                    Smagorinsky_model(dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, dTdx, dTdy, dTdz,
                                    rho_av, u_av, v_av, w_av, volume_av, yplus_av, i, j)
                end

                F_vis_hat[i,j,k,1] = 0.0
                F_vis_hat[i,j,k,2] = ((vecAy[i,j,k,1]*sigma_xx + vecAy[i,j,k,2]*sigma_xy + vecAy[i,j,k,3]*sigma_xz) -
                                      (vecAy[i,j,k,1]*tau_xx   + vecAy[i,j,k,2]*tau_xy   + vecAy[i,j,k,3]*tau_xz)) / volume_av
                F_vis_hat[i,j,k,3] = ((vecAy[i,j,k,1]*sigma_xy + vecAy[i,j,k,2]*sigma_yy + vecAy[i,j,k,3]*sigma_yz) -
                                      (vecAy[i,j,k,1]*tau_xy   + vecAy[i,j,k,2]*tau_yy   + vecAy[i,j,k,3]*tau_yz)) / volume_av
                F_vis_hat[i,j,k,4] = ((vecAy[i,j,k,1]*sigma_xz + vecAy[i,j,k,2]*sigma_yz + vecAy[i,j,k,3]*sigma_zz) -
                                      (vecAy[i,j,k,1]*tau_xz   + vecAy[i,j,k,2]*tau_yz   + vecAy[i,j,k,3]*tau_zz)) / volume_av
                F_vis_hat[i,j,k,5] = ((vecAy[i,j,k,1]*betax    + vecAy[i,j,k,2]*betay    + vecAy[i,j,k,3]*betaz) -
                                      (vecAy[i,j,k,1]*e_sgs_x  + vecAy[i,j,k,2]*e_sgs_y  + vecAy[i,j,k,3]*e_sgs_z))/ volume_av
            end
        end
    end

    for k in 1+icell:cellzmax+1 -icell
        for j in 1+icell:cellymax -icell
            for i in 1+icell:cellxmax -icell
                rho_av = 0.5*(QbaseF[i,j,k,1] + QbaseB[i,j,k,1])
                u_av   = 0.5*(QbaseF[i,j,k,2] + QbaseB[i,j,k,2])
                v_av   = 0.5*(QbaseF[i,j,k,3] + QbaseB[i,j,k,3])
                w_av   = 0.5*(QbaseF[i,j,k,4] + QbaseB[i,j,k,4])

                mu_av     = 0.5*(mu[i,j,k-1] + mu[i,j,k])
                lambda_av = (0.5*(1/lambda[i,j,k-1] + 1/lambda[i,j,k]))^(-1)

                volume_av = 0.5*(volume[i,j,k] + volume[i,j,k-1])

                dudxi = 0.25 * (QbaseR[i+1,  j,  k,2] - QbaseL[i+1,  j,  k,2] + QbaseR[i,  j,  k,2] - QbaseL[i,  j,  k,2]
                              + QbaseR[i+1,  j,k-1,2] - QbaseL[i+1,  j,k-1,2] + QbaseR[i,  j,k-1,2] - QbaseL[i,  j,k-1,2])
                dvdxi = 0.25 * (QbaseR[i+1,  j,  k,3] - QbaseL[i+1,  j,  k,3] + QbaseR[i,  j,  k,3] - QbaseL[i,  j,  k,3]
                              + QbaseR[i+1,  j,k-1,3] - QbaseL[i+1,  j,k-1,3] + QbaseR[i,  j,k-1,3] - QbaseL[i,  j,k-1,3])
                dwdxi = 0.25 * (QbaseR[i+1,  j,  k,4] - QbaseL[i+1,  j,  k,4] + QbaseR[i,  j,  k,4] - QbaseL[i,  j,  k,4]
                              + QbaseR[i+1,  j,k-1,4] - QbaseL[i+1,  j,k-1,4] + QbaseR[i,  j,k-1,4] - QbaseL[i,  j,k-1,4])

                dT1 = QbaseR[i+1,  j,  k,5]/(QbaseR[i+1,  j,  k,1]*Rd) - QbaseL[i+1,  j,  k,5]/(QbaseL[i+1,  j,  k,1]*Rd)
                dT2 = QbaseR[  i,  j,  k,5]/(QbaseR[  i,  j,  k,1]*Rd) - QbaseL[  i,  j,  k,5]/(QbaseL[  i,  j,  k,1]*Rd)
                dT3 = QbaseR[i+1,  j,k-1,5]/(QbaseR[i+1,  j,k-1,1]*Rd) - QbaseL[i+1,  j,k-1,5]/(QbaseL[i+1,  j,k-1,1]*Rd)
                dT4 = QbaseR[  i,  j,k-1,5]/(QbaseR[  i,  j,k-1,1]*Rd) - QbaseL[  i,  j,k-1,5]/(QbaseL[  i,  j,k-1,1]*Rd)
                dTdxi = 0.25 * (dT1 + dT2 + dT3 + dT4)

                dudeta = 0.25 * (QbaseU[  i,j+1,  k,2] - QbaseD[  i,j+1,  k,2] + QbaseU[i,  j,  k,2] - QbaseD[i,  j,  k,2]
                               + QbaseU[  i,j+1,k-1,2] - QbaseD[  i,j+1,k-1,2] + QbaseU[i,  j,k-1,2] - QbaseD[i,  j,k-1,2])
                dvdeta = 0.25 * (QbaseU[  i,j+1,  k,3] - QbaseD[  i,j+1,  k,3] + QbaseU[i,  j,  k,3] - QbaseD[i,  j,  k,3]
                               + QbaseU[  i,j+1,k-1,3] - QbaseD[  i,j+1,k-1,3] + QbaseU[i,  j,k-1,3] - QbaseD[i,  j,k-1,3])
                dwdeta = 0.25 * (QbaseU[  i,j+1,  k,4] - QbaseD[  i,j+1,  k,4] + QbaseU[i,  j,  k,4] - QbaseD[i,  j,  k,4]
                               + QbaseU[  i,j+1,k-1,4] - QbaseD[  i,j+1,k-1,4] + QbaseU[i,  j,k-1,4] - QbaseD[i,  j,k-1,4])

                dT1 = QbaseU[  i,j+1,  k,5]/(QbaseU[  i,j+1,  k,1]*Rd) - QbaseD[  i,j+1,  k,5]/(QbaseD[  i,j+1,  k,1]*Rd)
                dT2 = QbaseU[  i,  j,  k,5]/(QbaseU[  i,  j,  k,1]*Rd) - QbaseD[  i,  j,  k,5]/(QbaseD[  i,  j,  k,1]*Rd)
                dT3 = QbaseU[  i,j+1,k-1,5]/(QbaseU[  i,j+1,k-1,1]*Rd) - QbaseD[  i,j+1,k-1,5]/(QbaseD[  i,j+1,k-1,1]*Rd)
                dT4 = QbaseU[  i,  j,k-1,5]/(QbaseU[  i,  j,k-1,1]*Rd) - QbaseD[  i,  j,k-1,5]/(QbaseD[  i,  j,k-1,1]*Rd)
                dTdeta = 0.25 * (dT1 + dT2 + dT3 + dT4)

                dudzeta = QbaseF[i,j,k,2] - QbaseB[i,j,k,2]
                dvdzeta = QbaseF[i,j,k,3] - QbaseB[i,j,k,3]
                dwdzeta = QbaseF[i,j,k,4] - QbaseB[i,j,k,4]
                dTdzeta = (QbaseF[i,j,k,5]/(QbaseF[i,j,k,1]*Rd) - QbaseB[i,j,k,5]/(QbaseB[i,j,k,1]*Rd))

                vecAx_xav    = 0.25*( vecAx[i,j,k,1] + vecAx[i+1,j,k,1] + vecAx[i,j,k-1,1] + vecAx[i+1,j,k-1,1] )
                vecAx_yav    = 0.25*( vecAx[i,j,k,2] + vecAx[i+1,j,k,2] + vecAx[i,j,k-1,2] + vecAx[i+1,j,k-1,2] )
                vecAx_zav    = 0.25*( vecAx[i,j,k,3] + vecAx[i+1,j,k,3] + vecAx[i,j,k-1,3] + vecAx[i+1,j,k-1,3] )

                vecAy_xav    = 0.25*( vecAy[i,j,k,1] + vecAy[i,j+1,k,1] + vecAy[i,j,k-1,1] + vecAy[i,j+1,k-1,1] )
                vecAy_yav    = 0.25*( vecAy[i,j,k,2] + vecAy[i,j+1,k,2] + vecAy[i,j,k-1,2] + vecAy[i,j+1,k-1,2] )
                vecAy_zav    = 0.25*( vecAy[i,j,k,3] + vecAy[i,j+1,k,3] + vecAy[i,j,k-1,3] + vecAy[i,j+1,k-1,3] )

                dudx = vecAx_xav * dudxi + vecAy_xav * dudeta + vecAz[i,j,k,1] * dudzeta
                dudy = vecAx_yav * dudxi + vecAy_yav * dudeta + vecAz[i,j,k,2] * dudzeta
                dudz = vecAx_zav * dudxi + vecAy_zav * dudeta + vecAz[i,j,k,3] * dudzeta

                dvdx = vecAx_xav * dvdxi + vecAy_xav * dvdeta + vecAz[i,j,k,1] * dvdzeta
                dvdy = vecAx_yav * dvdxi + vecAy_yav * dvdeta + vecAz[i,j,k,2] * dvdzeta
                dvdz = vecAx_zav * dvdxi + vecAy_zav * dvdeta + vecAz[i,j,k,3] * dvdzeta

                dwdx = vecAx_xav * dwdxi + vecAy_xav * dwdeta + vecAz[i,j,k,1] * dwdzeta
                dwdy = vecAx_yav * dwdxi + vecAy_yav * dwdeta + vecAz[i,j,k,2] * dwdzeta
                dwdz = vecAx_zav * dwdxi + vecAy_zav * dwdeta + vecAz[i,j,k,3] * dwdzeta

                sigma_xx = 2/3*mu_av*( 2*dudx -   dvdy -   dwdz )
                sigma_yy = 2/3*mu_av*( - dudx + 2*dvdy -   dwdz )
                sigma_zz = 2/3*mu_av*( - dudx -   dvdy + 2*dwdz )

                sigma_xy = mu_av*( dudy + dvdx )
                sigma_xz = mu_av*( dvdz + dwdy )
                sigma_yz = mu_av*( dwdx + dudz )

                dTdx = vecAx_xav*dTdxi + vecAy_xav*dTdeta + vecAz[i,j,k,1]*dTdzeta
                dTdy = vecAx_yav*dTdxi + vecAy_yav*dTdeta + vecAz[i,j,k,2]*dTdzeta
                dTdz = vecAx_zav*dTdxi + vecAy_zav*dTdeta + vecAz[i,j,k,3]*dTdzeta

                betax = u_av*sigma_xx + v_av*sigma_xy + w_av*sigma_xz + lambda_av * dTdx
                betay = u_av*sigma_xy + v_av*sigma_yy + w_av*sigma_yz + lambda_av * dTdy
                betaz = u_av*sigma_xz + v_av*sigma_yz + w_av*sigma_zz + lambda_av * dTdz

                yplus_av = 0.5 *(yplus[i,j,k] + yplus[i,j,k-1])

                if swith_wall[5] == 1 && i == 1+icell
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,3] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
                elseif swith_wall[6] == 1 && i == cellzmax - icell
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,3] = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
                else
                    tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz,
                    e_sgs_x, e_sgs_y, e_sgs_z, mut_bd[i,j,k,3] =
                    Smagorinsky_model(dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, dTdx, dTdy, dTdz,
                                    rho_av, u_av, v_av, w_av, volume_av, yplus_av, i, j)
                end

                G_vis_hat[i,j,k,1] = 0.0
                G_vis_hat[i,j,k,2] = ((vecAz[i,j,k,1]*sigma_xx + vecAz[i,j,k,2]*sigma_xy + vecAz[i,j,k,3]*sigma_xz) -
                                      (vecAz[i,j,k,1]*tau_xx   + vecAz[i,j,k,2]*tau_xy   + vecAz[i,j,k,3]*tau_xz)) / volume_av
                G_vis_hat[i,j,k,3] = ((vecAz[i,j,k,1]*sigma_xy + vecAz[i,j,k,2]*sigma_yy + vecAz[i,j,k,3]*sigma_yz) -
                                      (vecAz[i,j,k,1]*tau_xy   + vecAz[i,j,k,2]*tau_yy   + vecAz[i,j,k,3]*tau_yz)) / volume_av
                G_vis_hat[i,j,k,4] = ((vecAz[i,j,k,1]*sigma_xz + vecAz[i,j,k,2]*sigma_yz + vecAz[i,j,k,3]*sigma_zz) -
                                      (vecAz[i,j,k,1]*tau_xz   + vecAz[i,j,k,2]*tau_yz   + vecAz[i,j,k,3]*tau_zz)) / volume_av
                G_vis_hat[i,j,k,5] = ((vecAz[i,j,k,1]*betax    + vecAz[i,j,k,2]*betay    + vecAz[i,j,k,3]*betaz) -
                                      (vecAz[i,j,k,1]*e_sgs_x  + vecAz[i,j,k,2]*e_sgs_y  + vecAz[i,j,k,3]*e_sgs_z))/ volume_av
            end
        end
    end

    # 乱流粘性の代入
    for k in 1:cellzmax
        for j in 1:cellymax
            for i in 1:cellxmax
                mut[i,j,k] = (mut_bd[i+1,j,k,1] + mut_bd[i,j,k,1] + mut_bd[i,j+1,k,2]
                            + mut_bd[i,j,k,2] + mut_bd[i,j,k+1,3] + mut_bd[i,j,k,3]) / 6.0
            end
        end
    end

    #=
    println("(u,v ,20)")
    println(Qbase[20,20,2])
    println(Qbase[19,20,2])
    println(Qbase[20,19,2])
    println(Qbase[19,19,2])
    println("(u,v,150)")
    println(Qbase[150,150,2])
    println(Qbase[149,150,2])
    println(Qbase[150,149,2])
    println(Qbase[149,149,2])
    =#

    return E_vis_hat, F_vis_hat, G_vis_hat,mut
end

function Smagorinsky_model(dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, dTdx, dTdy, dTdz,
                            rho_av, u_av, v_av, w_av, volume_av, yplus_av, i, j)
    Cs = 0.18
    Cl = 0.09
    Pr_sgs = 0.6
    k = 0.41

    S11 = 0.5*(dudx + dudx)
    S12 = 0.5*(dudy + dvdx)
    S13 = 0.5*(dudz + dwdx)

    S21 = 0.5*(dvdx + dudy)
    S22 = 0.5*(dvdy + dvdy)
    S23 = 0.5*(dvdz + dwdy)

    S31 = 0.5*(dwdx + dudz)
    S32 = 0.5*(dwdy + dvdz)
    S33 = 0.5*(dwdz + dwdz)

    absS = ( 2 * (S11^2 + S12^2 + S13^2 + S21^2 + S22^2 + S23^2 + S31^2 + S32^2 + S33^2) )^0.5
    Delta = (volume_av) ^ (1/3) * wallf_Van_Driest(yplus_av)

    #lsgs = 0.5*((Cs*Delta - k*y) - abs(Cs*Delta - k*y)) + k*y # min([ky, Cs*Delta])
    #nu_sgs = (lsgs)^2 * absS
    nu_sgs = (Cs*Delta)^2 * absS
    mu_sgs = nu_sgs * rho_av

    #=
    if i == 20 && j == 20
        println("(20,20)")
        println(S11)
        println(S12)
        println(S22)
        println(Delta)
        println(absS)
        println(nu_sgs/volume_av)
    end
    if i == 150 && j == 150
        println("(150,150)")
        println(S11)
        println(S12)
        println(S22)
        println(Delta)
        println(absS)
        println(nu_sgs/volume_av)
    end
    =#

    # SGS stress
    #tau_xx = -2 * rho_av * nu_sgs * 2/3 * dudx +  1/3 * Cl * 2 * rho_av * Delta^2 * absS^2
    #tau_xy = - rho_av * nu_sgs * (dvdx + dudy)
    #tau_yy = -2 * rho_av * nu_sgs * 2/3 * dvdy +  1/3 * Cl * 2 * rho_av * Delta^2 * absS^2
    tau_xx = -2 * rho_av * nu_sgs * S11
    tau_xy = -2 * rho_av * nu_sgs * S12
    tau_xz = -2 * rho_av * nu_sgs * S13

    tau_yy = -2 * rho_av * nu_sgs * S22
    tau_yz = -2 * rho_av * nu_sgs * S23

    tau_zz = -2 * rho_av * nu_sgs * S33

    # SGS energy
    e_sgs_x = -rho_av * nu_sgs / Pr_sgs * dTdx + tau_xx*u_av + tau_xy*v_av + tau_xz*w_av
    e_sgs_y = -rho_av * nu_sgs / Pr_sgs * dTdy + tau_xy*u_av + tau_yy*v_av + tau_yz*w_av
    e_sgs_z = -rho_av * nu_sgs / Pr_sgs * dTdz + tau_xz*u_av + tau_yz*v_av + tau_zz*w_av

    #
    tau_xx = tau_xx / volume_av   *1.0e0
    tau_xy = tau_xy / volume_av   *1.0e0
    tau_xz = tau_xz / volume_av   *1.0e0
    tau_yy = tau_yy / volume_av   *1.0e0
    tau_yz = tau_yz / volume_av   *1.0e0
    tau_zz = tau_zz / volume_av   *1.0e0
    e_sgs_x = e_sgs_x / volume_av *1.0e0
    e_sgs_y = e_sgs_y / volume_av *1.0e0
    e_sgs_z = e_sgs_z / volume_av *1.0e0

    #=
    if i == 51 && j ==3
        println(" iroiro ")
        println(absS)
        println(nu_sgs)
        println(-2 * rho_av * nu_sgs * 2/3 * dudx)
        println(1/3 * Cl * 2 * rho_av * Delta^2 * absS^2)
    end
    =#

    #tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz, e_sgs_x, e_sgs_y, e_sgs_z, mu_sgs = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    return tau_xx, tau_xy, tau_xz, tau_yy, tau_yz, tau_zz, e_sgs_x, e_sgs_y, e_sgs_z, mu_sgs
end


text:
function read_binary_vec(io::IO)
    vec = Float32[]
    while !eof(io)
        push!(vec, read(io, Float32))
    end
    vec
end

function infer_shape{F<:AbstractFloat}(vec::Vector{F}, shape::Tuple{Vararg{Int}})

    #=
    Setting a dimension to -1 will have us infer the other one
    =#

    l = length(vec)

    tot = 1
    dim_to_infer = 0
    for (i,v) in enumerate(shape)
        if v == -1
            @assert(dim_to_infer == 0)
            dim_to_infer = i
        else
            tot *= v
        end
    end

    if dim_to_infer != 0
        @assert(mod(l, tot) == 0)
        dims = collect(shape)
        dims[dim_to_infer] = div(l, tot)
        shape = tuple(dims...)
    end

    shape
end

function convert_to_column_major_array{F<:AbstractFloat}(vec::Vector{F}, shape::Tuple{Int})
    @assert(shape[1] == length(vec))
    deepcopy(vec)
end
function convert_to_column_major_array{F<:AbstractFloat}(vec::Vector{F}, shape::Tuple{Int, Int})
    n, m = shape
    @assert(n*m == length(vec))

    retval = Array(F, n, m)

    count = 0
    for i in 1 : n
        for j in 1: m
            count += 1
            retval[i,j] = vec[count]
        end
    end

    retval
end
function convert_to_column_major_array{F<:AbstractFloat}(vec::Vector{F}, shape::Tuple{Int, Int, Int})
    m, n, o = shape
    @assert(n*m*o == length(vec))

    retval = Array(F, m, n, o)

    count = 0
    for i in 1 : m
        for j in 1 : n
            for k in 1 : o
                count += 1
                retval[i,j, k] = vec[count]
            end
        end
    end

    retval
end
function convert_to_column_major_array{F<:AbstractFloat}(vec::Vector{F}, shape::Tuple{Int, Int, Int, Int})
    m, n, o, p = shape
    @assert(n*m*o*p == length(vec))

    retval = Array(F, m, n, o, p)

    count = 0
    for i in 1 : m
        for j in 1 : n
            for k in 1 : o
                for l in 1 : p
                    count += 1
                    retval[i,j, k, l] = vec[count]
                end
            end
        end
    end

    retval
end

ext: jl


text:
"""
    write_netcdf2d(fname,var,lon,lat,varname,filenc)

Write field `var`to netcdf file `fname`, with coordinates `lon` and `lat` and variable name `varname`.
File `filenc` is used to copy metadata from.
"""
function write_netcdf2d(fname,var,lon,lat,varname,filenc)

if(length(size(lon))==1)
  (nsx,)=size(lon)
  (nsy,)=size(lat)
else
  (nsx,nsy)=size(lon)
end
if(length(size(var))==2)
      nt=1
      fnotime=1
      tim=[0]
else
      (ns,ns,nt)=size(var)
      fnotime=0
end

mode=NC_NETCDF4

#println("Reading ",filenc)
ncin= NetCDF.open(filenc)
if( haskey(ncin.vars,"lon") )
   lonatt= ncin.vars["lon"].atts
   latatt= ncin.vars["lat"].atts
elseif (haskey(ncin.vars,"longitude"))
   lonatt= ncin.vars["longitude"].atts
   latatt= ncin.vars["latitude"].atts
elseif (haskey(ncin.vars,"x"))
   lonatt= ncin.vars["x"].atts
   latatt= ncin.vars["y"].atts
elseif (haskey(ncin.vars,"X"))
   lonatt= ncin.vars["X"].atts
   latatt= ncin.vars["Y"].atts
else
   println("Input reference file does not contain lon or longitude dimensional variables")
   quit(1)
end

if(fnotime==0)
tim=NetCDF.readvar(ncin,"time");
timeatt= ncin.vars["time"].atts
end
varatt= ncin.vars[varname].atts
NetCDF.close(ncin)

#ncgetatt(filenc, "global", gatts)

varattr=Dict()
for k in collect(keys(varatt))
   if((k!="_FillValue")&&(k!="add_offset")&&(k!="scale_factor"))
       varattr[k]=varatt[k]
   end
end

isfile(fname) ? rm(fname) : nothing

if(length(size(lon))==1)
   if(fnotime==0)
      nccreate( fname, varname , "lon" , nsx, "lat", nsy, "time", tim[1:nt] ,timeatt, atts=varattr,mode=mode,t=NC_FLOAT);
   else
      nccreate( fname, varname , "lon" , nsx, "lat", nsy, atts=varattr,mode=mode,t=NC_FLOAT);
   end
   nccreate( fname, "lon" , "lon" , nsx,  atts=lonatt,mode=mode,t=NC_FLOAT);
   nccreate( fname, "lat" , "lat" , nsy,  atts=latatt,mode=mode,t=NC_FLOAT);
else
   if(fnotime==0)
      nccreate( fname, varname , "x" , nsx, "y", nsy, "time", tim[1:nt] ,timeatt, atts=varattr,mode=mode,t=NC_FLOAT);
   else
      nccreate( fname, varname , "x" , nsx, "y", nsy, atts=varattr,mode=mode,t=NC_FLOAT);
   end
   nccreate( fname, "lon" , "x" , nsx, "y", nsy,  atts=lonatt,mode=mode,t=NC_FLOAT);
   nccreate( fname, "lat" , "x" , nsx, "y", nsy,  atts=latatt,mode=mode,t=NC_FLOAT);
end

ncwrite(var,fname,varname)
ncwrite(lon,fname,"lon")
ncwrite(lat,fname,"lat")

ncclose(fname)
end


text:
"""
    OPFBlocks(
        opfdata::OPFData,
        rawdata::RawData;
        modelinfo::ModelInfo = ModelInfo(),
        backend=JuMPBlockBackend,
        algparams::AlgParams = AlgParams()
    )

Create a structure `OPFBlocks` to decompose the original OPF problem
specified in `opfdata` timestep by timestep, by dualizing the ramping
constraint. One block corresponds
to one optimization subproblem (and hence, to a particular timestep),
and the attribute `blkCount` enumerates the total number of subproblems.
The subproblems are specified using `AbstractBlockModel` objects,
allowing to define them either with JuMP
(if `backend=JuMPBlockBackend` is chosen) or with ExaPF (`backend=ExaBlockBackend`).

### Decomposition by contingencies

By default, OPFBlocks decomposes the problem only timestep by
timestep (single-period multiple-contingency scheme),
leading to a total of `T` subproblems.
However, if the option `algparams.decompCtgs` is set to `true`,
the original problem is also decomposed contingency by contingency
(single-period single-contingency scheme).
In this case the total number of subproblems is `T * K` (with `K` the
total number of contingencies).

### Deporting the resolution on the GPU

When the backend is set to `ExaBlockBackend` (and a CUDA GPU is available), the user
could chose to deport the resolution of each subproblem directly on
the GPU simply by setting `algparams.device=CUDADevice`. However, note that
we could not instantiate more subproblems on the GPU than the number of GPU
available.

"""
mutable struct OPFBlocks <: AbstractBlocks
    blkCount::Int64
    blkIndex::CartesianIndices
    blkModel::Vector{AbstractBlockModel}
    colCount::Int64
    colValue::Array{Float64,2}
end

function OPFBlocks(
    opfdata::OPFData,
    rawdata::RawData;
    modelinfo::ModelInfo = ModelInfo(),
    backend=JuMPBlockBackend,
    algparams::AlgParams = AlgParams(),
    comm::Union{MPI.Comm,Nothing}
)
    ngen  = length(opfdata.generators)
    nbus  = length(opfdata.buses)
    nline = length(opfdata.lines)
    T = modelinfo.num_time_periods
    K = modelinfo.num_ctgs + 1 # base case counted separately

    blkCount = algparams.decompCtgs ? (T*K) : T
    blkIndex = algparams.decompCtgs ? CartesianIndices((1:K,1:T)) : CartesianIndices((1:1,1:T))
    colCount = ((algparams.decompCtgs ? 1 : K)*(
                    2*nbus + 4*ngen + 1)) + 2*ngen
    colValue = zeros(colCount,blkCount)
    blkModel = AbstractBlockModel[]

    function local_copy(modelinfo)
        info = deepcopy(modelinfo)
        info.num_time_periods = 1
        return info
    end

    function load_local_data(rawdata, opfdata, modelinfo, t, k; decompCtgs=false)
        lineOff = Line()
        if length(rawdata.ctgs_arr) < k - 1
            error("Not enough contingencies in .ctg file while trying to read contingency $(k-1).")
        end
        if decompCtgs
            if k > 1
                lineOff = opfdata.lines[rawdata.ctgs_arr[k - 1]]
            end
        end
        data = opf_loaddata(
            rawdata;
            time_horizon_start=t,
            time_horizon_end=t,
            load_scale=modelinfo.load_scale,
            ramp_scale=modelinfo.ramp_scale,
            corr_scale=modelinfo.corr_scale,
            lineOff=lineOff,
        )
        return data
    end

    for blk in LinearIndices(blkIndex)
        k = blkIndex[blk][1]
        t = blkIndex[blk][2]
        if is_my_work(blk, comm)

            # Local info
            localinfo = local_copy(modelinfo)
            if algparams.decompCtgs
                @assert k > 0
                localinfo.num_ctgs = 0
            end
            localdata = load_local_data(rawdata, opfdata, localinfo, t, k;
                                        decompCtgs=algparams.decompCtgs)
            # Create block model
            localmodel = backend(blk, localdata, rawdata, algparams, localinfo, t, k, T)
        else
            localmodel = EmptyBlockModel()
        end
        push!(blkModel, localmodel)
    end

    return OPFBlocks(blkCount, blkIndex, blkModel, colCount, colValue)
end

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
# The configuration file format was updated in v0.6.0.
# However, this configuration file still seems to be in the old format.
#
# To update your configuration file to the new format (v0.6.0 and newer),
# open a new Julia session and load this file through
#
#     include("<path_to_SolidStateDetectors.jl>/test/update_config_files.jl")
#
# Afterwards, run
#
#     update_config_file("<path_to_configuration_file>")
#
# This method returns the file name of the updated configuration file.
# Please close the Julia session after updating the configuration files, as some
# parsing methods are overridden with old methods.

include("../src/IO/Update.jl")


text:
module SHA2

##############################################################################
##
## TODO
##
##############################################################################

# - Refactor RIPEMD-160 and SHA-256 to share boilerplate (E.g., read/write, padding)
# - Digest -> init, update, finalize

##############################################################################
##
## Refereces
##
##############################################################################

# - http://en.wikipedia.org/wiki/SHA-2

##############################################################################
##
## Exports
##
##############################################################################

export sha256

##############################################################################
##
## Notes
##
##############################################################################

# WARNING: This has only been tested on OS X 10.0.1 on a 64-bit machine.
#          The implementation below may die on a 32-bit machine.
# WARNING: This implementation assumes string inputs

# - All variables are 32-bit unsigned integers. From the Julia documentation
#
#     http://docs.julialang.org/en/release-0.3/manual/integers-and-floating
#     -point-numbers/
#
#   we see that
#
#   Unsigned integers are input and output using the 0x prefix and hexadecimal
#   (base 16) digits 0-9a-f (the capitalized digits A-F also work for input).
#   The size of the unsigned value is determined by the number of hex digits.
#
# - Addition is calculated modulo 2^32
#
# - For each round, there is one round constant k[i] and one entry in the
#   message schedule array w[i], 0 ≤ i ≤ 63
#
# - The compression function uses 8 working variables, a through h
#
# - Big-endian convention is used when expressing the constants in this
#   pseudocode, and when parsing message block data from bytes to chunks,
#   for example, the first chunk of the input message "abc" after padding is
#   0x61626380

##############################################################################
##
## Algorithm parameters
##
##############################################################################

const DIGEST_SIZE = 32 # SHA256 outputs a 32 byte (256 bit) digest
const BLOCK_SIZE = 64  # SHA256 operates on 64 (512 bit) blocks
const CHUNKS_PER_BLOCK = 16
const CHUNKS_PER_SCHEDULE = 64

##############################################################################
##
## Initialize array of round constants
##
##############################################################################

# First 32 bits of the fractional parts of the cube roots of the first 64
# primes, 2 through 311. A sample generation function below:
#
# function initial_array_of_round_constraints(n)
#   fractional_cuberoot = cbrt(n) % 1
#   first_32_bits = floor(fractional_cube_root * 2^32)
#   @printf "%x" first_32_bits
#   return first_32_bits
# end

const k = [0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1,
           0x923f82a4, 0xab1c5ed5, 0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3,
           0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174, 0xe49b69c1, 0xefbe4786,
           0x0fc19dc6, 0x240ca1cc, 0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
           0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7, 0xc6e00bf3, 0xd5a79147,
           0x06ca6351, 0x14292967, 0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13,
           0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85, 0xa2bfe8a1, 0xa81a664b,
           0xc24b8b70, 0xc76c51a3, 0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
           0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5, 0x391c0cb3, 0x4ed8aa4a,
           0x5b9cca4f, 0x682e6ff3, 0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208,
           0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2]

##############################################################################
##
## Initialize helper functions
##
##############################################################################

# Hard-coding for 32-bit unsigned ints
# Technically, we could use the '>>' operator for right shift because
# we're only using unsigned ints, but we explicitly use logical right shift
ROTRIGHT(num, shift) = (num >>> shift) | (num << (32 - shift))
ROTLEFT(num, shift) = (num << shift) | (num >>> (32 - shift))

CH(x, y, z) = z $ (x & (y $ z))
MA(x, y, z) = ((x | y) & z) | (x & y)

E0(x) = (ROTRIGHT(x, 2) $ ROTRIGHT(x, 13) $ ROTRIGHT(x, 22))
E1(x) = (ROTRIGHT(x, 6) $ ROTRIGHT(x, 11) $ ROTRIGHT(x, 25))

S0(x) = (ROTRIGHT(x, 7) $ ROTRIGHT(x, 18) $ (x >>> 3))
S1(x) = (ROTRIGHT(x, 17) $ ROTRIGHT(x, 19) $ (x >>> 10))

##############################################################################
##
## Function definitions
##
##############################################################################

function transform!(state, block)

  # Pre-allocate the message schedule array (64 8-bit words)
  m = zeros(Uint32, CHUNKS_PER_SCHEDULE)

  for i = 1:CHUNKS_PER_BLOCK
    m[i] = uint32(block[4 * (i - 1) + 1]) << 24 +
           uint32(block[4 * (i - 1) + 2]) << 16 +
           uint32(block[4 * (i - 1) + 3]) << 8 +
           uint32(block[4 * (i - 1) + 4])
  end
  for i = CHUNKS_PER_BLOCK + 1:CHUNKS_PER_SCHEDULE
    s0 = S0(m[i-15])
    s1 = S1(m[i-2])

    m[i] = m[i-16] + s0 + m[i-7] + s1
  end

  a, b, c, d, e, f, g, h = state

  for i = 1:CHUNKS_PER_SCHEDULE
    t1 = h + E1(e) + CH(e, f, g) + k[i] + m[i]
    t2 = E0(a) + MA(a, b, c)

    h = g;
    g = f;
    f = e;
    e = uint32(d + t1)
    d = c;
    c = b;
    b = a;
    a = uint32(t1 + t2)
  end

  state[1] += a
  state[2] += b
  state[3] += c
  state[4] += d
  state[5] += e
  state[6] += f
  state[7] += g
  state[8] += h

  return
end

function sha256(msg::String; is_hex=false)

  if is_hex
    msg = [uint8(parseint(msg[2*i-1:2*i], 16)) for i in 1:length(msg)/2]
  else
    # We only want byte array literal (i.e., character array)
    msg = msg.data
  end

  # Get original length and bit lengths
  len = length(msg)
  bitlen = len * 8

  # Append the bit '1' to the message.
  append!(msg, [0x80])

  # First 32 bits of the fractional parts of the square roots of the first 8
  # primes, 2 through 19. A sample generation function below:
  #
  # function initial_hash_value(n)
  #   fractional_square_root = sqrt(n) % 1
  #   first_32_bits = floor(fractional_square_root * 2^32)
  #   @printf "%x" first_32_bits
  #   return first_32_bits
  # end

  state = [0x6a09e667, 0xbb67ae85, 0x3c6ef372, 0xa54ff53a,
           0x510e527f, 0x9b05688c, 0x1f83d9ab, 0x5be0cd19]

  # Divide up message into blocks of BLOCK_SIZE = 512 bits
  # and run through transformation
  while length(msg) >= BLOCK_SIZE
    transform!(state, msg[1:BLOCK_SIZE])
    msg = msg[BLOCK_SIZE + 1:end]
  end

  # Get the number of characters untransformed
  rem = length(msg)

  # If there are any characters remaining
  if rem > 0

    # Append k bits '0', where k is the minimum number >= 0 such that the
    # resulting message length (modulo 512 in bits) is 448.
    if length(msg) > BLOCK_SIZE - 8
      msg = append!(msg, zeros(Uint8, BLOCK_SIZE - rem))
      transform!(state, msg)
      msg = zeros(Uint8, BLOCK_SIZE)
    else
      msg = append!(msg, zeros(Uint8, BLOCK_SIZE - rem))
    end

    # Append length of message (without the '1' bit or padding), in bits, as
    # 64-bit big-endian integer (this will make the entire post-processed
    # length a multiple of 512 bits)
    msg[57] = (bitlen >>> 56) & 0xff
    msg[58] = (bitlen >>> 48) & 0xff
    msg[59] = (bitlen >>> 40) & 0xff
    msg[60] = (bitlen >>> 32) & 0xff
    msg[61] = (bitlen >>> 24) & 0xff
    msg[62] = (bitlen >>> 16) & 0xff
    msg[63] = (bitlen >>> 8) & 0xff
    msg[64] = bitlen & 0xff

    # Process the last block
    transform!(state, msg)
  end

  # Assemble digest and return
  return reduce((x, y) -> string(x, y), map(x -> lpad(hex(uint32(x)), 8, "0"), state))
end

end # module SHA256


text:
# An example of semi-supervised node classification

using Flux
using Flux: onecold, onehotbatch
using Flux.Losses: logitcrossentropy
using GraphNeuralNetworks
using MLDatasets: Cora
using Statistics, Random
using CUDA
CUDA.allowscalar(false)

function eval_loss_accuracy(X, y, ids, model, g)
    ŷ = model(g, X)
    l = logitcrossentropy(ŷ[:,ids], y[:,ids])
    acc = mean(onecold(ŷ[:,ids] |> cpu) .== onecold(y[:,ids] |> cpu))
    return (loss = round(l, digits=4), acc = round(acc*100, digits=2))
end

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 1f-3             # learning rate
    epochs = 100          # number of epochs
    seed = 17             # set seed > 0 for reproducibility
    usecuda = true      # if true use cuda (if available)
    nhidden = 128        # dimension of hidden features
    infotime = 10 	     # report every `infotime` epochs
end

function train(; kws...)
    args = Args(; kws...)

    args.seed > 0 && Random.seed!(args.seed)

    if args.usecuda && CUDA.functional()
        device = gpu
        args.seed > 0 && CUDA.seed!(args.seed)
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    # LOAD DATA
    data = Cora.dataset()
    g = GNNGraph(data.adjacency_list) |> device
    X = data.node_features |> device
    y = onehotbatch(data.node_labels, 1:data.num_classes) |> device
    train_ids = data.train_indices |> device
    val_ids = data.val_indices |> device
    test_ids = data.test_indices |> device
    ytrain = y[:,train_ids]

    nin, nhidden, nout = size(X,1), args.nhidden, data.num_classes

    ## DEFINE MODEL
    model = GNNChain(GCNConv(nin => nhidden, relu),
                     Dropout(0.5),
                     GCNConv(nhidden => nhidden, relu),
                     Dense(nhidden, nout))  |> device

    ps = Flux.params(model)
    opt = ADAM(args.η)

    @info g

    ## LOGGING FUNCTION
    function report(epoch)
        train = eval_loss_accuracy(X, y, train_ids, model, g)
        test = eval_loss_accuracy(X, y, test_ids, model, g)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
    end

    ## TRAINING
    report(0)
    for epoch in 1:args.epochs
        gs = Flux.gradient(ps) do
            ŷ = model(g, X)
            logitcrossentropy(ŷ[:,train_ids], ytrain)
        end

        Flux.Optimise.update!(opt, ps, gs)

        epoch % args.infotime == 0 && report(epoch)
    end
end

train()


text:
#=
Test every function once in sequence.
Sequence could matter ?
=#
#using Test
using Spotify
using Albums, Artists, Browse, Episodes, Follow, Library
using Markets, Personalization, Player, Playlists, Search
using  Shows, Tracks, UsersProfile, Objects

#= Albums =#
album_get(album_id, market="US")
album_get_tracks(album_id, limit=20, offset=0, market="US")

#= Artists =#
artist_get(artist_id)
artist_get_albums(artist_id, include_groups="None", country="US", limit=20, offset=0)
artist_top_tracks(artist_id, country="US")
artist_get_related_artists(artist_id)

#= Browse =#
category_get(category_id, country="US", locale="en")
category_get_playlist(category_id, country="US", limit = 20, offset = 0)
category_get_several(country="US", locale="en", limit=20, offset=0)
category_get_new_releases(country="US", locale="en", limit=20, offset=0)
category_get_featured_playlist(country="US",locale="en", limit=50, offset=0, timestamp=string(now())[1:19])
recommendations_get(
    seeds,
    track_attributes,
    limit = 50,
    market = "US"
)
recommendations_dict_parser(track_attributes::Dict)

#= Episodes =#
episodes_get(episode_id, market="US")

#= Follow =#
follow_check(type, ids)
follow_check_playlist(playlist_id, ids)
follow_artists(type="artist", limit=20)
follow_artists_users(type, ids)
follow_playlist(playlist_id)
unfollow_artists_users(type, ids)
unfollow_playlist(playlist_id)

#= Library =#
library_get_saved_tracks(limit=20, offset=0, market="US")
library_get_saved_shows(limit=20, offset=0)
library_get_saved_albums(limit=20, offset=0, market="US")
library_check_saved_tracks(track_ids)
library_check_saved_shows(show_ids)
library_check_saved_albums(album_ids)
library_remove_albums(album_ids)
library_remove_shows(show_ids)
library_remove_tracks(track_ids)
library_save_album(album_ids)
library_save_show(shows_ids)
library_save_track(track_ids)

#= Personalization =#
top_tracks(offset=0, limit=20, time_range="medium")
top_artists(offset=0, limit=20, time_range="medium")

#= Player =#
show_get(show_id, market="US")
show_get_episodes(show_id, market="US", limit=20, offset=0)

#= Playlists =#
tracks_get_audio_analysis(track_id)
tracks_get_audio_features(track_id)
tracks_get(track_id, market="US")

#= Search =#
#= Shows =#
#= Tracks =#
#= UsersProfile =#
#= Objects =#


text:
Base.@kwdef struct PokemonBattleResult
    battles::Vector{Battle}
    parsed_battles::Vector{ParsedBattle}
    time_created::DateTime = now()
end


text:
struct NanoDate <: Dates.AbstractDateTime
    datetime::DateTime      # DateTime(UTM(microsec))
    nanosecs::Nanosecond    # Nanosecond(microsecond * 1_000 + nanosecond)

    function NanoDate(datetime::DateTime, nanosecs::Nanosecond)
        nanosec = value(nanosecs)
        0 <= nanosec < 1_000_000 && return new(datetime, Nanosecond(nanosec))

        datetime2, nanosecs2 = canonical(millis=value(datetime), nanos=nanosec)
        new(datetime2, nanosecs2)
    end
end

function canonical_mn(millis::T, nanos::T) where {T<:Integer}
    micros, nanos = fldmod_1000(nanos)
    micromillis, micros = fldmod_1000(micros)
    millis += micromillis
    nanos = micros * NanosecondsPerMicrosecond + nanos
    return Millisecond(millis), Nanosecond(nanos)
end

function canonical_cn(micros::T, nanos::T) where {T<:Integer}
    nanomicros, nanos = fldmod_1000(nanos)
    micros += nanomicros
    micromillis, micros = fldmod_1000(micros)
    nanos = micros * NanosecondsPerMicrosecond + nanos
    return Millisecond(micromillis), Nanosecond(nanos)
end

function canonical_mc(millis::T, micros::T) where {T<:Integer}
    micromills, micros = fldmod_1000(micros)
    millis += micromillis
    nanos = micros * NanosecondsPerMicrosecond
    return Millisecond(millis), Nanosecond(nanos)
end

function canonical_mcn(millis::T, micros::T, nanos::T) where {T<:Integer}
    micronanos, nanos = fldmod_1000(nanos)
    micros += micronanos
    millimicros, micros = fldmod_1000(micros)
    millis += millimicros
    nanos = micros * NanosecondsPerMicrosecond + nanos
    return Millisecond(millis), Nanosecond(nanos)
end

canonical(micros::Microsecond, nanos::Nanosecond) =
    canonical_cn(value(micros), value(nanos))

canonical(millis::Millisecond, nanos::Nanosecond) =
    canonical_mn(value(millis), value(nanos))

canonical(millis::Millisecond, micros::Microsecond) =
    canonical_mc(value(millis), value(micros))

canonical(millis::Millisecond, micros::Microsecond, nanos::Nanosecond) =
    canonical_mcn(value(millis), value(micros), value(nanos))

@inline function canonical(;millis, nanos)
    micros, nanos = fldmod_1000(nanos)
    micromillis, micros = fldmod_1000(micros)
    millis += micromillis
    nanos += micros * NanosecondsPerMicrosecond
    datetime = DateTime(Dates.UTM(millis))
    nanosecs = Nanosecond(nanos)
    datetime, nanosecs
end

datetime(x::NanoDate) = x.datetime
nanosecs(x::NanoDate) = x.nanosecs

NanoDate(x::NanoDate) = x

#### Additional Constructors

function NanoDate(date::Date, time::Time)
    slowtime, fasttime = fldmod(value(time), NanosecondsPerMillisecond)
    datetime = DateTime(date) + Millisecond(slowtime)
    NanoDate(datetime, Nanosecond(fasttime))
end

@inline NanoDate(time::Time, date::Date) = NanoDate(date, time)

NanoDate(dt::DateTime, cs::Microsecond) = NanoDate(dt, Nanosecond(value(cs) * 1_000))
NanoDate(dt::DateTime, cs::Microsecond, ns::Nanosecond) =
    NanoDate(dt, nanosecs(cs, ns))

NanoDate(dt::Date, cs::Microsecond, ns::Nanosecond) = NanoDate(DateTime(dt), cs, ns)
NanoDate(dt::Date, cs::Microsecond) = NanoDate(DateTime(dt), cs)
NanoDate(dt::Date, ns::Nanosecond) = NanoDate(DateTime(dt), ns)

NanoDate(yr::Year, mn::Month=Month(1), dy::Day=Day(1)) = NanoDate(Date(yr, mn, dy))
NanoDate(yr::Year, mn::Month, dy::Day, hr::Hour, mi::Minute=Minute(0), sc::Second=Second(0),
    ms::Millisecond=Millisecond(0), cs::Microsecond=Microsecond(0), ns::Nanosecond=Nanosecond(0)) =
    NanoDate(DateTime(yr, mn, dy, hr, mi, sc, ms), nanosecs(cs, ns))

NanoDate(yr, mn=1, dy=1, hr=0, mi=0, sc=0, ms=0, cs=0, ns=0) =
    NanoDate(DateTime(yr, mn, dy, hr, mi, sc, ms), nanosecs(cs, ns))

NanoDate(cs::Microsecond, ns::Nanosecond) = NanoDate(DateTime0, nanosecs(cs,ns))

# for internal cse

csns(ns) = ns
@inline csns(cs, ns) = cs * 1_000 + ns

nanosecs(ns) = Nanosecond(ns)
nanosecs(ns::Nanosecond) = ns
nanosecs(cs::Microsecond) = Nanosecond(1_000 * value(cs))

@inline nanosecs(cs, ns) = Nanosecond(csns(cs, ns))
@inline nanosecs(cs::Microsecond, ns::Nanosecond) = nanosecs(value(cs), value(ns))

NanoDate(dt::DateTime, x) = NanoDate(dt, nanosecs(x))
NanoDate(dt::DateTime, cs, ns) = NanoDate(dt, nanosecs(cs, ns))
NanoDate(d::Date, x) = NanoDate(d, nanosecs(x))
NanoDate(d::Date, cs, ns) = NanoDate(d, nanosecs(cs, ns))

NanoDate(dy::Day) = NanoDate(Date(UTD(dy)), Nanosecond0)
NanoDate(dy::Day, ns::Nanosecond) = NanoDate(Date(UTD(dy)), ns)
NanoDate(dy::Day, cs::Microsecond) = NanoDate(Date(UTD(dy)), nanosecs(cs))
NanoDate(dy::Day, cs::Microsecond, ns::Nanosecond) = NanoDate(Date(UTD(dy)), nanosecs(cs, ns))

NanoDate(ms::Millisecond) = NanoDate(DateTime(UTM(ms)), Nanosecond0)
NanoDate(ms::Millisecond, ns::Nanosecond) = NanoDate(DateTime(UTM(ms)), ns)
NanoDate(ms::Millisecond, cs::Microsecond) = NanoDate(DateTime(UTM(ms)), nanosecs(cs))
NanoDate(ms::Millisecond, cs::Microsecond, ns::Nanosecond) = NanoDate(DateTime(UTM(ms)), nanosecs(cs, ns))

NanoDate() = NanoDate(now())
NanoDate(::Type{UTC}) = NanoDate(now(UTC))

for T in (:Year, :Quarter, :Month, :Week, :Day,
            :Hour, :Minute, :Second, :Millisecond)
    @eval NanoDate(nd::NanoDate, x::$T) = NanoDate(nd.datetime - $T(nd) + x, nd.nanosecs)
end

for T in (:Microsecond, :Nanosecond)
    @eval NanoDates.NanoDate(nd::NanoDate, x::$T) = nd - $T(nd) + x
end

NanoDate(nd::NanoDate, cs::Microsecond, ns::Nanosecond) =
    nd - Microsecond(nd) - Nanosecond(nd) + ns + cs

NanoDate(nd::NanoDate, dtm::DateTime) =
    NanoDate(dtm, nd.nanosecs)
NanoDate(nd::NanoDate, tm::Time) =
    NanoDate(Date(nd.datetime),tm)
NanoDate(nd::NanoDate, dt::Date) =
    NanoDate(DateTime(dt, trunc(Time(nd.datetime), Millisecond)), nd.nanosecs)

const NanoDate0 = NanoDate(0,1,1,0,0,0,0,0,0)


text:
function extract_dload!(
    dloads::Dict{String,TF},
    facesets::Dict{String,Vector{Tuple{TI,TI}}},
    file,
    ::Type{Val{dim}},
    offset::TI,
) where {TI,TF,dim}
    pattern = r"(\d+)\s*,\s*P(\d+)\s*,\s*(\-?\d+\.\d*)"
    dload_heading_pattern = r"\*DLOAD"

    faceset_name = "DLOAD_SET_$(length(dloads)+1)"
    facesets[faceset_name] = Tuple{TI,TI}[]

    first = true
    prevload = zero(TF)
    load = zero(TF)

    line = readline(file)
    m = match(stopping_pattern, line)
    while m isa Nothing
        m = match(dload_heading_pattern, line)
        if m != nothing
            dloads[faceset_name] = load
            first = true
            faceset_name = "DLOAD_SET_$(length(dloads)+1)"
            facesets[faceset_name] = Tuple{TI,TI}[]
        end
        m = match(pattern, line)
        if m != nothing
            cellidx = parse(TI, m[1]) - offset
            faceidx = parse(TI, m[2])
            load = parse(TF, m[3])
            if !first && prevload != load
                throw("Loads in the same DLOAD set are not equal.")
            end
            prevload = load
            if first
                first = false
            end
            push!(facesets[faceset_name], (cellidx, faceidx))
        end
        line = readline(file)
        m = match(stopping_pattern, line)
    end
    dloads[faceset_name] = load

    return line
end


text:
module MyExampleILIC

include("extra_file.jl")

export my_f, greetings

end


text:
module P4est

using CBinding: @cstruct
using Reexport: @reexport

include("libp4est.jl")

"""
    uses_mpi()

Returns true if the p4est library was compiled with MPI enabled.
"""
uses_mpi() = isdefined(@__MODULE__, :P4EST_ENABLE_MPI)

end


text:
"""
$(SIGNATURES)

Find the closest point to `p` within `hull`. If `p` is inside `hull`,
`p` itself is returned.
"""
function closest_point(p::PointLike, hull::ConvexHull)
    op = orientation_comparator(hull)
    vertices = hull.vertices
    n = length(vertices)
    @inbounds begin
        if n === 0
            throw(ArgumentError())
        elseif n === 1
            return vertices[1]
        else
            closest = p
            i = 1
            v1 = vertices[n]
            v2 = vertices[1]
            while true
                δ = v2 - v1
                p′ = p - v1
                strictly_outside_edge = op(cross2(p′, δ), 0)
                if strictly_outside_edge
                    # Find closest point to current edge.
                    λ = clamp(p′ ⋅ δ / (δ ⋅ δ), false, true)
                    closest_to_edge = v1 + λ * δ

                    # Accept if previous candidate was p itself or if it's closer.
                    if closest == p || normsquared(closest_to_edge - p) < normsquared(closest - p)
                        closest = closest_to_edge
                    end
                end
                i == n && break
                v1 = vertices[i]
                v2 = vertices[i + 1]
                i += 1
            end
            return closest
        end
    end
end


text:
# Fusion trees:
#==============================================================================#
struct FusionTree{G<:Sector,N,M,L,T}
    outgoing::NTuple{N,G}
    incoming::G
    innerlines::NTuple{M,G} # M = N-2
    vertices::NTuple{L,T} # L = N-1
end
FusionTree{G,N}(outgoing::NTuple{N}, incoming, innerlines::NTuple{M}, vertices::NTuple{L}) where {G<:Sector,N,M,L} =
    FusionTree{G,N,M,L,vertex_labeltype(G)}(outgoing, incoming, innerlines, vertices)
FusionTree{G,N}(outgoing, incoming, innerlines) where {G<:Sector,N} =
    FusionTree{G,N}(outgoing, incoming, innerlines, ntuple(n->nothing, valsub(Val(N),Val(1))))

FusionTree{G}(outgoing::NTuple{N}, incoming, innerlines::NTuple{M}, vertices::NTuple{L}) where {G<:Sector,N,M,L} =
    FusionTree{G,N,M,L,vertex_labeltype(G)}(outgoing, incoming, innerlines, vertices)
FusionTree{G}(outgoing::NTuple{N}, incoming, innerlines) where {G<:Sector,N} =
    FusionTree{G,N}(outgoing, incoming, innerlines, ntuple(n->nothing, valsub(Val(N),Val(1))))

FusionTree(outgoing::NTuple{N,G}, incoming::G, innerlines::NTuple{M,G}) where {G<:Sector, N, M} = FusionTree{G,N}(outgoing, incoming, innerlines, ntuple(n->nothing,valsub(Val(N),Val(1))))

function FusionTree(outgoing::NTuple{N,G}, incoming::G = first(⊗(outgoing...))) where {G<:Sector, N}
    @assert fusiontype(G) == Abelian
    FusionTree{G,N}(outgoing, incoming, _abelianinner((outgoing...,dual(incoming))))
end

# Properties
sectortype(::Type{<:FusionTree{G}}) where {G<:Sector} = G
fusiontype(::Type{<:FusionTree{G}}) where {G<:Sector} = fusiontype(G)
Base.length(::Type{<:FusionTree{<:Sector,N}}) where {N} = N

sectortype(t::FusionTree) = sectortype(typeof(t))
fusiontype(t::FusionTree) = fusiontype(typeof(t))
Base.length(t::FusionTree) = length(typeof(t))

# Fusion tree methods

# permute fusion tree
"""
    function permute(t::FusionTree{<:Sector,N}, i) where {N} -> (Immutable)Dict{typeof(t),<:Number}

Performs a permutation of the outgoing indices of the fusion tree `t` and returns the result
as a `Dict` (or `ImmutableDict`) of output trees and corresponding coefficients.
"""
function Base.permute(t::FusionTree{G,N}, p::NTuple{N,Int}) where {G<:Sector, N}
    @assert braidingtype(G) <: SymmetricBraiding

    if fusiontype(G) == Abelian
        coeff = Rsymbol(one(G), one(G), one(G))
        for i = 1:N
            for j = 1:i-1
                if p[j] > p[i]
                    a, b = t.outgoing[j], t.outgoing[i]
                    coeff *= Rsymbol(a, b, first(a ⊗ b))
                end
            end
        end
        t = FusionTree(tpermute(t.outgoing, p), t.incoming)
        return ImmutableDict(t=>coeff)
    else
        coeff = Rsymbol(one(G), one(G), one(G))
        trees = Dict(t=>coeff)
        newtrees = similar(trees)
        for s in permutation2swaps(p)
            for (t, c) in trees
                for (t′,c′) in braid(t, s)
                    newtrees[t′] = get(newtrees, t′, zero(coeff)) + c*c′
                end
            end
            trees, newtrees = newtrees, trees
            empty!(newtrees)
        end
        return trees
    end
end

"""
    function braid(t::FusionTree{<:Sector,N}, i) where {N} -> ImmutableDict{typeof(t),<:Number}

Performs a braid of neighbouring outgoing indices `i` and `i+1` on a fusion tree `t`,
and returns the result as a linked list of output trees and corresponding coefficients.
"""
function braid(t::FusionTree{<:Sector,N}, i) where {N}
    1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
    outer = t.outgoing
    inner = t.innerlines
    if i == 1
        a, b = outer[1], outer[2]
        c = N > 2 ? inner[1] : t.incoming
        outer = setindex(outer, b, 1)
        outer = setindex(outer, a, 2)
        return ImmutableDict(FusionTree(outer, t.incoming, inner, t.vertices) => Rsymbol(a, b, c))
    end
    # case i > 1:
    b = outer[i]
    d = outer[i+1]
    a = i == 2 ? outer[1] : inner[i-2]
    c = inner[i-1]
    e = i == N-1 ? t.incoming : inner[i]
    outer′ = setindex(outer, d, i)
    outer′ = setindex(outer′, b, i+1)
    if fusiontype(t) == Abelian
        inner′ = setindex(inner, first(a ⊗ d), i-1)
        return ImmutableDict(FusionTree(outer′, t.incoming, inner′) => Rsymbol(b, d, first(b ⊗ d)))
    elseif fusiontype(t) == SimpleNonAbelian
        iter = a ⊗ d
        s = start(iter)
        c′, s = next(iter, s)
        while iszero(Nsymbol(b, c′, e))
            c′, s = next(iter, s)
        end
        coeff = conj(Rsymbol(b,a,c))*Fsymbol(b,a,d,e,c,c′)*Rsymbol(b,c′,e)
        inner′ = setindex(inner, c′, i-1)
        output = ImmutableDict(FusionTree(outer′, t.incoming, inner′, t.vertices)=>coeff)
        while !done(iter, s)
            c′, s = next(iter, s)
            iszero(Nsymbol(b, c′, e)) && continue
            coeff = conj(Rsymbol(b,a,c))*Fsymbol(b,a,d,e,c,c′)*Rsymbol(b,c′,e)
            inner′ = setindex(inner, c′, i-1)
            output = ImmutableDict(output, FusionTree(outer′, t.incoming, inner′, t.vertices) => coeff)
        end
        return output
    else
        # TODO: DegenerateNonAbelian case
    end
end

# Fusion tree iterators
include("iterator.jl")

# Show methods
function Base.show(io::IO, t::FusionTree{G,0}) where {G<:Sector}
    if get(io, :compact, false)
        println(io, "FusionTree{", G, ", 0}((), (), (), ())")
    else
        println(io, "FusionTree{", G, ", 0}: (empty fusion tree)")
    end
end
function Base.show(io::IO, t::FusionTree{G,N,M,K,Void}) where {G<:Sector,N,M,K}
    if get(io, :compact, false)
        print(io, "FusionTree{", G, ", ", N, "}(", t.outgoing, ", ", t.incoming, ", ", t.innerlines, ")")
        return
    end
    println(io, "FusionTree{", G, ", ", N, "}:")
    cprint = x->sprint(showcompact, x)
    outlabels = map(cprint, t.outgoing)
    L = max(map(strwidth, outlabels)...)
    up = lpad(rpad(" ∧ ", L>>1), L)
    down = lpad(rpad(" ∨ ", L>>1), L)
    L = strwidth(up)
    f = lpad(rpad("⋅+-",L>>1,"-"), L, " ")
    m = lpad(rpad("-+-",L>>1,"-"), L, "-")
    a = "-<-"
    pad = " "^L
    line1 = ""
    line2 = ""
    line3 = ""
    line4 = ""
    n = 1
    while n <= N
        label = lpad(rpad(outlabels[n], 2), L)
        line4 *= label
        line3 *= down
        if n == 1
            line2 *= f
        else
            line2 *= m
        end
        line1 *= pad
        # print inner label
        if n == 1
            label = cprint(t.outgoing[1])
            label = lpad(rpad(label, 2), strwidth(a))
            line1 *= label
            ℓ = strwidth(label)
            line2 *= lpad(rpad(a, ℓ>>1, "-"), ℓ, "-")
            line3 *= " "^ℓ
            line4 *= " "^ℓ
        elseif n < N
            label = cprint(t.innerlines[n-1])
            label = lpad(rpad(label, 2), strwidth(a))
            line1 *= label
            ℓ = strwidth(label)
            line2 *= lpad(rpad(a, ℓ>>1, "-"), ℓ, "-")
            line3 *= " "^ℓ
            line4 *= " "^ℓ
        else
            line2 *= a
            line2 *= cprint(t.incoming)
        end
        n += 1
    end
    println(io, line1)
    println(io, line2)
    println(io, line3)
    println(io, line4)
end

# auxiliary routines

# _abelianinner: generate the inner indices for given outer indices in the abelian case
_abelianinner(outer::Tuple{}) = ()
_abelianinner(outer::Tuple{G}) where {G<:Sector} = ()
_abelianinner(outer::Tuple{G,G}) where {G<:Sector} = ()
_abelianinner(outer::Tuple{G,G,G}) where {G<:Sector} = ()
function _abelianinner(outer::NTuple{N,G}) where {G<:Sector,N}
    c = first(outer[1] ⊗ outer[2])
    return (c, _abelianinner((c, tail2(outer)...))...)
end


text:
#*****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Experiment 3
test combined time series and binary logical operations
"""
module NNGrammarExpt3

export circuit_fgandor, restarts

import Compat.ASCIIString
using TFTools
using Datasets
using TensorFlow
using TensorFlow.CoreTypes
import TensorFlow.API: l2_loss, AdamOptimizer, cast, round_, reshape_,
    reduce_max, reduce_min, pack, minimum_, maximum_, transpose_, reduce_sum,
    moments, gradients
using StatsBase
using RLESUtils, Confusion

using PyCall
@pyimport tensorflow as tf

function restarts(f::Function, N::Int64; kwargs...)
   [f(; kwargs...) for i = 1:N]
end

function circuit_fgandor(;
    dataname::AbstractString="bin_ts_synth",
    labelfield::AbstractString="F_x1_and_x3",
    learning_rate::Float64=0.001,
    max_training_epochs::Int64=500,
    target_cost::Float64=0.001,
    batch_size::Int64=5000,
    embed_hidden_units::Vector{Int64}=Int64[50, 30, 15],
    mux_hidden_units::Vector{Int64}=Int64[50,30],
    display_step::Int64=1,
    b_debug::Bool=false,
    nshow::Int64=20)

    Ds = dataset(dataname) #DFSet
    data_set = TFDataset(Ds, getmeta(Ds)[Symbol(labelfield)])

    # Construct model
    (n_examples, n_steps, n_feats) = size(Ds)
    n_featsflat = n_steps * n_feats

    # inputs
    feats = Placeholder(DT_FLOAT32, [-1, n_steps, n_feats])
    inputs = Tensor(feats)
    feats_flat = reshape_(feats, Tensor([-1, n_featsflat]))

    # common (embedding) layer
    #embed_in = feats_flat
    #embed_blk = ReluStack(embed_in, embed_hidden_units)
    #embed_out = out(embed_blk)

    # mux select input
    muxselect = feats_flat #embed layer not helping, disable for now...

    #softness parameter
    softness = collect(linspace(0.1, 1.0, max_training_epochs))
    softness_pl = Placeholder(DT_FLOAT32, [1])

    #toggle hard or soft output
    harden_pl = Placeholder(DT_BOOL, [1])

    # x mux
    x_muxin = inputs
    x_mux = SoftMux(n_feats, mux_hidden_units, x_muxin, muxselect, Tensor(harden_pl), Tensor(softness_pl))
    x_muxout = out(x_mux)

    # y mux
    y_muxin = inputs
    y_mux = SoftMux(n_feats, mux_hidden_units, y_muxin, muxselect, Tensor(harden_pl), Tensor(softness_pl))
    y_muxout = out(y_mux)

    # logical op block
    ops1_in = (x_muxout, y_muxout)
    ops1_list = [ops1_And, ops1_Or]
    ops1_blk = SoftOpsMux(ops1_in, ops1_list, mux_hidden_units, muxselect, Tensor(harden_pl), Tensor(softness_pl))
    ops1_out = out(ops1_blk)

    # temporal op block
    ops2_in = (ops1_out,)
    ops2_list = [ops2_F, ops2_G]
    ops2_blk = SoftOpsMux(ops2_in, ops2_list, mux_hidden_units, muxselect, Tensor(harden_pl), Tensor(softness_pl))
    ops2_out = out(ops2_blk)

    # outputs
    pred = ops2_out
    labels = Placeholder(DT_FLOAT32, [-1])

    #take nnout of batch, compute moments, take variance component and sum
    x_var = reduce_sum(moments(x_mux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    y_var = reduce_sum(moments(y_mux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    ops1_var = reduce_sum(moments(ops1_blk.softmux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    ops2_var = reduce_sum(moments(ops2_blk.softmux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    sum_var = x_var + y_var + constant(10.0) .* ops1_var + constant(10.0) .* ops2_var

    # Define loss and optimizer
    #cost = l2_loss(pred - labels) # Squared loss
    cost = l2_loss(pred - labels) + constant(50.0) .* sum_var

    #optimizer
    #optimizer = minimize(AdamOptimizer(learning_rate), cost) # Adam Optimizer
    opt = Optimizer(tf.train[:GradientDescentOptimizer](learning_rate))
    gvs = opt.x[:compute_gradients](cost.x)
    capped_gvs = [(tf.nn[:l2_normalize](tf.clip_by_value(grad, -1.0, 1.0), 0), var) for (grad, var) in gvs]
    optimizer = Operation(opt.x[:apply_gradients](capped_gvs))

    #compiled hardselect
    ckt = Circuit([
        x_mux, y_mux,
        ops1_blk,
        ops2_blk,
        ],
        Vector{ASCIIString}[
        recordcolnames(Ds),
        recordcolnames(Ds),
        ["and", "or"],
        ["F", "G"]])

    # Initializing the variables
    init = initialize_all_variables()

    # Rock and roll
    sess = Session()
    run(sess, init)

    # Training cycle
    for epoch in 1:max_training_epochs
        avg_cost = 0.0
        total_batch = div(num_examples(data_set), batch_size)

        # Loop over all batches
        for i in 1:total_batch
            batch_xs, batch_ys = next_batch(data_set, batch_size)
            fd = FeedDict(feats => batch_xs, labels => batch_ys, softness_pl => [softness[epoch]], harden_pl => [false])
            # Fit training using batch data
            run(sess, optimizer, fd)

            #debug
            #@show run(sess, pred - labels, fd)
            #/debug

            # Compute average loss
            batch_average_cost = run(sess, cost, fd)
            avg_cost += batch_average_cost / (total_batch * batch_size)
        end

        # Display logs per epoch step
        if epoch % display_step == 0
            fd = FeedDict(feats => data_set.X[1,:,:], labels => [data_set.Y[1]], softness_pl => [softness[epoch]], harden_pl => [false])
            softsel = softselect_by_example(sess, ckt, fd)
            x_grad = run(sess, Tensor(gradients(cost, x_mux.nnout)), fd)
            y_grad = run(sess, Tensor(gradients(cost, y_mux.nnout)), fd)
            a1_grad = run(sess, Tensor(gradients(cost, ops1_blk.softmux.nnout)), fd)
            t1_grad = run(sess, Tensor(gradients(cost, ops2_blk.softmux.nnout)), fd)
            println("x_sel=", softsel[1][1])
            println("x_grad=",x_grad[1])
            println("y_sel=", softsel[1][2])
            println("y_grad=",y_grad[1])
            println("a1_sel=", softsel[1][3])
            println("a1_grad=",a1_grad[1])
            println("t1_sel=", softsel[1][4])
            println("t1_grad=",t1_grad[1])
            println("Epoch $(epoch)  cost=$(avg_cost)")
            if avg_cost < target_cost
                break;
            end
        end
    end
    println("Optimization Finished")

    # Test model
    correct_prediction = (round_(pred) == labels)
    accuracy = mean(cast(correct_prediction, DT_FLOAT32))

    #reload data_set to recover original order
    data_set = TFDataset(Ds, getmeta(Ds)[Symbol(labelfield)])
    X = data_set.X
    Y = data_set.Y

    #soft metrics
    fd = FeedDict(feats => data_set.X, labels => data_set.Y, softness_pl => [softness[end]], harden_pl => [false])
    Ypred_soft = run(sess, round_(pred), fd)
    acc_soft = run(sess, accuracy, fd)
    stringout = simplestring(sess, ckt, fd; order=[4,1,3,2])
    top5 = topstrings(stringout, 5)
    softsel = softselect_by_example(sess, ckt, fd)
    #grads = run(sess, Tensor([f1_grad, v1_grad, a1_grad, t1_grad]), fd)
    conf = confusion(Ypred_soft.==1.0, Y.==1.0)
    conf_indices = confusion_indices(Ypred_soft.==1.0, Y.==1.0)

    #hard metrics
    fd = FeedDict(feats => data_set.X, labels => data_set.Y, softness_pl => [softness[end]], harden_pl => [true])
    Ypred_hard = run(sess, round_(pred), fd)
    acc_hard = run(sess, accuracy, fd)

    println("Soft Accuracy:", acc_soft)
    println("Hard Accuracy:", acc_hard)
    println(top5)
    d=Dict{ASCIIString,Any}()

    top5, acc_hard
end

function ops2_F(x::Tensor)
   reduce_max(x, Tensor(1))
end

function ops2_G(x::Tensor)
   reduce_min(x, Tensor(1))
end

function ops1_And(x::Tensor, y::Tensor)
    minimum_(x, y) #element-wise
end

function ops1_Or(x::Tensor, y::Tensor)
    maximum_(x, y) #element-wise
end

end #module


text:
import DataStructures.OrderedSet
import DataStructures.OrderedDict

const EXEC = ["exec"]

baremodule Aggregate
    const NotSet = ""
    const Sum = "sum"
    const Min = "min"
    const Max = "max"
end

# Key commands
@redisfunction "del" Integer key...
@redisfunction "dump" AbstractString key
@redisfunction "exists" Bool key
@redisfunction "expire" Bool key seconds
@redisfunction "expireat" Bool key timestamp
@redisfunction "keys" Set{AbstractString} pattern
@redisfunction "migrate" Bool host port key destinationdb timeout
@redisfunction "move" Bool key db
@redisfunction "persist" Bool key
@redisfunction "pexpire" Bool key milliseconds
@redisfunction "pexpireat" Bool key millisecondstimestamp
@redisfunction "pttl" Integer key
@redisfunction "randomkey" Union{AbstractString, Nothing}
@redisfunction "rename" AbstractString key newkey
@redisfunction "renamenx" Bool key newkey
@redisfunction "restore" Bool key ttl serializedvalue
@redisfunction "scan" Tuple{Integer, Array{AbstractString, 1}} cursor::Integer options...
@redisfunction "sort" Array{AbstractString, 1} key options...
@redisfunction "ttl" Integer key
function keytype(conn::RedisConnection, key)
    response = execute_command(conn, flatten_command("type", key))
    convert_response(AbstractString, response)
end
function keytype(conn::TransactionConnection, key)
    execute_command(conn, flatten_command("type", key))
end

# String commands
@redisfunction "append" Integer key value
@redisfunction "bitcount" Integer key options...
@redisfunction "bitop" Integer operation destkey key keys...
@redisfunction "bitpos" Integer key bit options...
@redisfunction "decr" Integer key
@redisfunction "decrby" Integer key decrement
@redisfunction "get" Union{String, Nothing} key
@redisfunction "getbit" Integer key offset
@redisfunction "getrange" AbstractString key start finish
@redisfunction "getset" AbstractString key value
@redisfunction "incr" Integer key
@redisfunction "incrby" Integer key increment::Integer

# Bulk string reply: the value of key after the increment,
# as per http://redis.io/commands/incrbyfloat
@redisfunction "incrbyfloat" AbstractString key increment::Float64
@redisfunction "mget" Array{Union{AbstractString, Nothing}, 1} key keys...
@redisfunction "mset" Bool keyvalues
@redisfunction "msetnx" Bool keyvalues
@redisfunction "psetex" AbstractString key milliseconds value
@redisfunction "set" Bool key value options...
@redisfunction "setbit" Integer key offset value
@redisfunction "setex" AbstractString key seconds value
@redisfunction "setnx" Bool key value
@redisfunction "setrange" Integer key offset value
@redisfunction "strlen" Integer key

# Hash commands
@redisfunction "hdel" Integer key field fields...
@redisfunction "hexists" Bool key field
@redisfunction "hget" Union{AbstractString, Nothing} key field
@redisfunction "hgetall" Dict{AbstractString, AbstractString} key
@redisfunction "hincrby" Integer key field increment::Integer

# Bulk string reply: the value of key after the increment,
# as per http://redis.io/commands/hincrbyfloat
@redisfunction "hincrbyfloat" AbstractString key field increment::Float64

@redisfunction "hkeys" Array{AbstractString, 1} key
@redisfunction "hlen" Integer key
@redisfunction "hmget" Array{Union{AbstractString, Nothing}, 1} key field fields...
@redisfunction "hmset" Bool key value
@redisfunction "hset" Bool key field value...
@redisfunction "hsetnx" Bool key field value
@redisfunction "hvals" Array{AbstractString, 1} key
@redisfunction "hscan" Tuple{Integer, Dict{AbstractString, AbstractString}} key cursor::Integer options...

# List commands
@redisfunction "blpop" Array{AbstractString, 1} keys timeout
@redisfunction "brpop" Array{AbstractString, 1} keys timeout
@redisfunction "brpoplpush" AbstractString source destination timeout
@redisfunction "lindex" Union{AbstractString, Nothing} key index
@redisfunction "linsert" Integer key place pivot value
@redisfunction "llen" Integer key
@redisfunction "lpop" Union{AbstractString, Nothing} key
@redisfunction "lpush" Integer key value values...
@redisfunction "lpushx" Integer key value
@redisfunction "lrange" Array{AbstractString, 1} key start finish
@redisfunction "lrem" Integer key count value
@redisfunction "lset" AbstractString key index value
@redisfunction "ltrim" AbstractString key start finish
@redisfunction "rpop" Union{AbstractString, Nothing} key
@redisfunction "rpoplpush" Union{AbstractString, Nothing} source destination
@redisfunction "rpush" Integer key value values...
@redisfunction "rpushx" Integer key value

# Set commands
@redisfunction "sadd" Integer key member members...
@redisfunction "scard" Integer key
@redisfunction "sdiff" Set{AbstractString} key keys...
@redisfunction "sdiffstore" Integer destination key keys...
@redisfunction "sinter" Set{AbstractString} key keys...
@redisfunction "sinterstore" Integer destination key keys...
@redisfunction "sismember" Bool key member
@redisfunction "smembers" Set{AbstractString} key
@redisfunction "smove" Bool source destination member
@redisfunction "spop" Union{AbstractString, Nothing} key
@redisfunction "srandmember" Union{AbstractString, Nothing} key
@redisfunction "srandmember" Set{AbstractString} key count
@redisfunction "srem" Integer key member members...
@redisfunction "sunion" Set{AbstractString} key keys...
@redisfunction "sunionstore" Integer destination key keys...
@redisfunction "sscan" Set{AbstractString} key cursor::Integer options...

# Sorted set commands
#=
merl-dev: a number of methods were added to take AbstractString for score value
to enable score ranges like '(1 2,' or "-inf", "+inf",
as per docs http://redis.io/commands/zrangebyscore
=#

@redisfunction "zadd" Integer key score::Number member::AbstractString

# NOTE:  using ZADD with Dicts could introduce bugs if some scores are identical
@redisfunction "zadd" Integer key scorememberdict

#=
This following version of ZADD enables adding new members using `Tuple{Int64, AbstractString}` or
`Tuple{Float64, AbstractString}` for single or multiple additions to the sorted set without
resorting to the use of `Dict`, which cannot be used in the case where all entries have the same score.
=#
@redisfunction "zadd" Integer key scoremembertup scorememberstup...

@redisfunction "zcard" Integer key
@redisfunction "zcount" Integer key min max

# Bulk string reply: the new score of member (a double precision floating point number),
# represented as string, as per http://redis.io/commands/zincrby
@redisfunction "zincrby" AbstractString key increment member

@redisfunction "zlexcount" Integer key min max
@redisfunction "zrange" OrderedSet{AbstractString} key start finish options...
@redisfunction "zrangebylex" OrderedSet{AbstractString} key min max options...
@redisfunction "zrangebyscore" OrderedSet{AbstractString} key min max options...
@redisfunction "zrank" Union{Integer, Nothing} key member
@redisfunction "zrem" Integer key member members...
@redisfunction "zremrangebylex" Integer key min max
@redisfunction "zremrangebyrank" Integer key start finish
@redisfunction "zremrangebyscore" Integer key start finish
@redisfunction "zrevrange" OrderedSet{AbstractString} key start finish options...
@redisfunction "zrevrangebyscore" OrderedSet{AbstractString} key start finish options...
@redisfunction "zrevrank" Union{Integer, Nothing} key member
# ZCORE returns a Bulk string reply: the score of member (a double precision floating point
# number), represented as string.
@redisfunction "zscore" Union{AbstractString, Nothing} key member
@redisfunction "zscan" Set{AbstractString} key cursor::Integer options...

# Stream commands
@redisfunction "xadd" Integer key ID fvpairs...
@redisfunction "xrange" Dict{AbstractString, Dict{AbstractString,AbstractString}} key start end_ options...
@redisfunction "xrevrange" OrderedDict{AbstractString, Dict{AbstractString,AbstractString}} key end_ start options...
@redisfunction "xlen" Integer key
@redisfunction "xread" Dict{AbstractString, OrderedDict{AbstractString, Dict{AbstractString, AbstractString}}} options...

function _build_store_internal(destination, numkeys, keys, weights, aggregate, command)
    length(keys) > 0 || throw(ClientException("Must supply at least one key"))
    suffix = []
    if length(weights) > 0
        suffix = map(string, weights)
        pushfirst!(suffix, "weights")
    end
    if aggregate != Aggregate.NotSet
        push!(suffix, "aggregate")
        push!(suffix, aggregate)
    end
    vcat([command, destination, numkeys], keys, suffix)
end

# TODO: PipelineConnection and TransactionConnection
function zinterstore(conn::RedisConnectionBase, destination, numkeys,
    keys::Array, weights=[]; aggregate=Aggregate.NotSet)
    command = _build_store_internal(destination, numkeys, keys, weights, aggregate, "zinterstore")
    execute_command(conn, command)
end

function zunionstore(conn::RedisConnectionBase, destination, numkeys::Integer,
    keys::Array, weights=[]; aggregate=Aggregate.NotSet)
    command = _build_store_internal(destination, numkeys, keys, weights, aggregate, "zunionstore")
    execute_command(conn, command)
end

# HyperLogLog commands
@redisfunction "pfadd" Bool key element elements...
@redisfunction "pfcount" Integer key keys...
@redisfunction "pfmerge" Bool destkey sourcekey sourcekeys...

# Connection commands
@redisfunction "auth" AbstractString password
@redisfunction "echo" AbstractString message
@redisfunction "ping" AbstractString
@redisfunction "quit" Bool
@redisfunction "select" AbstractString index

# Transaction commands
@redisfunction "discard" Bool
@redisfunction "exec" Array{Bool} # only one element ever in this array?
@redisfunction "multi" Bool
@redisfunction "unwatch" Bool
@redisfunction "watch" Bool key keys...

# Scripting commands
# TODO: PipelineConnection and TransactionConnection
function evalscript(conn::RedisConnection, script, numkeys::Integer, args)
    response = execute_command(conn, flatten_command("eval", script, numkeys, args))
    convert_eval_response(Any, response)
end

#################################################################
# TODO: NEED TO TEST BEYOND THIS POINT
@redisfunction "evalsha" Any sha1 numkeys keys args
@redisfunction "script_exists" Array script scripts...
@redisfunction "script_flush" AbstractString
@redisfunction "script_kill" AbstractString
@redisfunction "script_load" AbstractString script

# Server commands
@redisfunction "bgrewriteaof" Bool
@redisfunction "bgsave" AbstractString
@redisfunction "client_getname" AbstractString
@redisfunction "client_id" Integer
@redisfunction "client_list" AbstractString
@redisfunction "client_pause" Bool timeout
@redisfunction "client_setname" Bool name
@redisfunction "cluster_slots" Array
@redisfunction "command" Array
@redisfunction "command_count" Integer
@redisfunction "command_info" Array command commands...
@redisfunction "config_get" Array parameter
@redisfunction "config_resetstat" Bool
@redisfunction "config_rewrite" Bool
@redisfunction "config_set" Bool parameter value
@redisfunction "dbsize" Integer
@redisfunction "debug_object" AbstractString key
@redisfunction "debug_segfault" Any
@redisfunction "flushall" AbstractString
@redisfunction "flushdb" AbstractString Integer
@redisfunction "info" AbstractString
@redisfunction "info" AbstractString section
@redisfunction "lastsave" Integer
@redisfunction "role" Array
@redisfunction "save" Bool
@redisfunction "shutdown" AbstractString
@redisfunction "shutdown" AbstractString option
@redisfunction "slaveof" AbstractString host port
@redisfunction "_time" Array{AbstractString, 1}

# Sentinel commands
@sentinelfunction "master" Dict{AbstractString, AbstractString} mastername
@sentinelfunction "reset" Integer pattern
@sentinelfunction "failover" Any mastername
@sentinelfunction "monitor" Bool name ip port quorum
@sentinelfunction "remove" Bool name
@sentinelfunction "set" Bool name option value

function sentinel_masters(conn::SentinelConnection)
    response = execute_command(conn, flatten_command("sentinel", "masters"))
    [convert_response(Dict, master) for master in response]
end

function sentinel_slaves(conn::SentinelConnection, mastername)
    response = execute_command(conn, flatten_command("sentinel", "slaves", mastername))
    [convert_response(Dict, slave) for slave in response]
end

function sentinel_getmasteraddrbyname(conn::SentinelConnection, mastername)
    execute_command(conn, flatten_command("sentinel", "get-master-addr-by-name", mastername))
end

# Custom commands (PubSub/Transaction)
@redisfunction "publish" Integer channel message

function _subscribe(conn::SubscriptionConnection, channels::Array)
    execute_command_without_reply(conn, pushfirst!(channels, "subscribe"))
end

function subscribe(conn::SubscriptionConnection, channel::AbstractString, callback::Function)
    conn.callbacks[channel] = callback
    _subscribe(conn, [channel])
end

function subscribe(conn::SubscriptionConnection, subs::Dict{AbstractString, Function})
    for (channel, callback) in subs
        conn.callbacks[channel] = callback
    end
    _subscribe(conn, collect(keys(subs)))
end

# REF: https://github.com/stejin/Redis.jl/commit/eace3a1ace7b464ba1e5ea074b27e6e8f652820e
function unsubscribe(conn::SubscriptionConnection, channels...)
    for channel in channels
        delete!(conn.callbacks, channel)
    end
    for channel in channels
        execute_command_without_reply(conn, unshift!([channel], "unsubscribe"))
    end
end
# function unsubscribe(conn::SubscriptionConnection, channels...)
#     for channel in channels
#         delete!(conn.callbacks, channel)
#     end
#     execute_command(conn, pushfirst!(channels, "unsubscribe"))
# end

function _psubscribe(conn::SubscriptionConnection, patterns::Array)
    execute_command_without_reply(conn, pushfirst!(patterns, "psubscribe"))
end

function psubscribe(conn::SubscriptionConnection, pattern::AbstractString, callback::Function)
    conn.callbacks[pattern] = callback
    _psubscribe(conn, [pattern])
end

function psubscribe(conn::SubscriptionConnection, subs::Dict{AbstractString, Function})
    for (pattern, callback) in subs
        conn.callbacks[pattern] = callback
    end
    _psubscribe(conn, collect(values(subs)))
end

function punsubscribe(conn::SubscriptionConnection, patterns...)
    for pattern in patterns
        delete!(conn.pcallbacks, pattern)
    end
    execute_command(conn, pushfirst!(patterns, "punsubscribe"))
end

#Need a specialized version of execute to keep the connection in the transaction state
function exec(conn::TransactionConnection)
    response = execute_command(conn, EXEC)
    multi(conn)
    response
end

###############################################################
# The following Redis commands can be typecast to Julia structs
###############################################################

function time(c::RedisConnection)
    t = _time(c)
    s = parse(Int,t[1])
    ms = parse(Float64, t[2])
    s += (ms / 1e6)
    return unix2datetime(s)
end


text:
module TensorNetworkAD
using Zygote, BackwardsLinalg
using OMEinsum

export trg, num_grad
export ctmrg
export optimiseipeps
export hamiltonian, model_tensor, mag_tensor
export Ising, TFIsing, Heisenberg

include("hamiltonianmodels.jl")

include("trg.jl")
include("fixedpoint.jl")
include("ctmrg.jl")
include("ipeps.jl")
include("autodiff.jl")
include("variationalipeps.jl")
include("exampletensors.jl")

end # module


text:
testpath = joinpath(dirname(dirname(pathof(PressureDrop))), "test/testdata")

@testset "Mayes 3-16 errors" begin
    default_GOR = 500
    max_depth = 8500

    GLpath = joinpath(testpath, "Mayes_3_16/Test_GLdesign_Mayes_3_16.csv")
    surveypath = joinpath(testpath, "Mayes_3_16/Test_survey_Mayes_3_16.csv")

    well = read_survey(path = surveypath, skiplines = 4, maxdepth = max_depth) #include IDs so that tailpipe and liner are properly accounted
    valves = read_valves(path = GLpath, skiplines = 1)

    model = WellModel(wellbore = well, roughness = 0.001,
                    valves = valves, pressurecorrelation = HagedornAndBrown,
                    WHP = 0, CHP = 0, dp_est = 25, temperature_method = "Shiu",
                    BHT = 160, geothermal_gradient = 0.8,
                    q_o = 0, q_w = 500,
                    GLR = 0, naturalGLR = 0,
                    APIoil = 38.2, sg_water = 1.05, sg_gas = 0.65)

    ## 2018-02-14 (domain error in H&B flow pattern)
    #duplicating the exact leadup from the loop where problem originally occurred, to avoid masking any problematic calculation
    oil = 41.3966666
    gas =  581.0
    water = 136.8298755
    TP = 138.0
    CP = 473.0
    gasinj = 533.1899861698156

    model.WHP = TP
    model.CHP = CP
    model.q_o = oil
    model.q_w = water
    model.naturalGLR = model.q_o + model.q_w <= 0 ? 0 : max( gas * 1000 / (model.q_o + model.q_w) , default_GOR * model.q_o / (model.q_o + model.q_w) )
    model.GLR = model.q_o + model.q_w <= 0 ? 0 : max( (gas + gasinj) * 1000 / (model.q_o + model.q_w) , model.naturalGLR)

    @test (gaslift_model!(model, find_injectionpoint = true, dp_min = 100) |> x-> x[1][end]) ≈ 678 atol = 1

    #2018-02-12 (no output)
    oil = 45.0
    gas =  584.0
    water = 149.0972222
    TP = 137.0
    CP = 474.0
    gasinj = 518.08798656154

    model.WHP = TP
    model.CHP = CP
    model.q_o = oil
    model.q_w = water
    model.naturalGLR = model.q_o + model.q_w <= 0 ? 0 : max( gas * 1000 / (model.q_o + model.q_w) , default_GOR * model.q_o / (model.q_o + model.q_w) )
    model.GLR = model.q_o + model.q_w <= 0 ? 0 : max( (gas + gasinj) * 1000 / (model.q_o + model.q_w) , model.naturalGLR)

    @test (gaslift_model!(model, find_injectionpoint = true, dp_min = 100) |> x-> x[1][end]) ≈ 697 atol = 1

    # nonconverging segment at 6395'
    oil = 2
    gas =  584.0
    water = 149.0972222
    TP = 137.0
    CP = 474.0
    gasinj = 518.08798656154

    model.WHP = TP
    model.CHP = CP
    model.q_o = oil
    model.q_w = water
    model.naturalGLR = model.q_o + model.q_w <= 0 ? 0 : max( gas * 1000 / (model.q_o + model.q_w) , default_GOR * model.q_o / (model.q_o + model.q_w) )
    model.GLR = model.q_o + model.q_w <= 0 ? 0 : max( (gas + gasinj) * 1000 / (model.q_o + model.q_w) , model.naturalGLR)

    @test (gaslift_model!(model, find_injectionpoint = true, dp_min = 100) |> x-> x[1][end]) ≈ 674 atol = 1
end