text:
module Zint

using StaticArrays: SVector, MVector
using ..Sheets: SV2, RWGSheet
using ..Layers: Layer
using ..RWG: RWGData
using ..PGF: jksums
using LinearAlgebra: ⋅, norm
using Statistics: mean

export filljk!, zint

zhatcross(t) = [-t[2], t[1]]
zhatcross(t::SV2) = SV2(-t[2], t[1])

"""
    filljk!(metal::RWGSheet, rwgdat::RWGData, closed::Bool)

Compute matrices of frequency-independent integrals needed in filling
the generalized impedance matrix.

## Arguments:

- `metal`: A variable of `RWGSheet` containing the face
  information, unit cell incremental phase shifts, and
  Green's function smoothing parameter.
- `rwgdat`:  A variable of type `RWGData` from which we make
  use of the arrays `ufpm` and `ufp2fp`.
- `closed`:    A logical flag which, if true, instructs this function
  to perform singularity extraction and analytic integration
  for all source triangles. If false then analytic singularity
  integration is performed ONLY for the case when source and
  observation triangles are the same.

##  Outputs:

- `metal`:  The following fields of `metal` are modified:
  `J`, `J_ξ`, `J_η`, `K`, `K_ξ`, `K_η`.  These are matrices
  of face-pair scalar integrals defined in Equations
  (7.22) through (7.27) of the theory documentation. Note these
  integrals are unitless.  Also modified are the fields `ρ_r` (the
  vector of singular vector integrals defined in Equation (B.3b),
  divided by twice the source triangle area to make them unitless,
  and `rinv`, a vector of singular scalar integrals defined in Eq.
  (B.3a), divided by u times twice the source triangle
  area to make them unitless.

"""
function filljk!(metal::RWGSheet, rwgdat::RWGData, closed::Bool)
    ξ = SVector(0.33333333333333330, 0.10128650732345633, 0.79742698535308730,
                0.10128650732345633, 0.47014206410511505, 0.05971587178976989,
                0.47014206410511505)
    η = SVector(0.33333333333333333, 0.79742698535308730, 0.10128650732345633,
                0.10128650732345633, 0.05971587178976989, 0.47014206410511505,
                0.47014206410511505)
    wght = SVector(0.1125, 0.06296959027241358, 0.06296959027241358, 0.06296959027241358,
                   0.06619707639425308, 0.06619707639425308, 0.06619707639425308)

    next = (2, 3, 1)

    nface = size(metal.fe, 2) # Number of faces in triangulated sheet.
    nedge = length(metal.e1) # Number of edges in triangulated sheet.
    nufp = rwgdat.nufp # Number of unique face pairs
    J = zeros(ComplexF64, nufp)
    J_ξ = zeros(ComplexF64, nufp)
    J_η = zeros(ComplexF64, nufp)
    K = zeros(ComplexF64, nufp)
    K_ξ = zeros(ComplexF64, nufp)
    K_η = zeros(ComplexF64, nufp)
    ρ_r = zeros(typeof(SV2(0.,0.)), nufp)
    rinv = zeros(Float64, nufp)

    i2s = CartesianIndices((nface,nface))

    ulocal = metal.u  # Obtain smoothing parameter in units that
    #                 # are consistent with metal's length units.
    us1 = ulocal * metal.s₁
    us2 = ulocal * metal.s₂

    Threads.@threads for iufp ∈ 1:rwgdat.nufp  # Loop over each unique face pair
        ifmifs = rwgdat.ufp2fp[iufp][1]  # Obtain index into face/face matrix
        rowcol = i2s[ifmifs]
        ifm, ifs = rowcol[1], rowcol[2] # indices of match and source triangles
        is = @view metal.fe[:,ifs]     # Obtain the three edges of the source triangle.

        #vtxcrd!(rs, ifs, metal) # Obtain coordinates of source tri's vertices
        rs = vtxcrd(ifs, metal) # Obtain coordinates of source tri's vertices
        # Calculate twice the signed area of the source triangle
        rs32 = rs[3] - rs[2]
        rs12 = rs[1] - rs[2]
        area2 = zhatcross(rs32) ⋅ rs12
        unz = 1.0   # Sign of unit surface normal (in z direction)
        if area2 < 0
            area2 = abs(area2)
            unz = -1.0
        end

        rs31 = rs[3] - rs[1]
        rs21 = -rs12

        self_tri = ifm == ifs
        clsflg = self_tri || closed  # Extract if self or if always
        im = @view metal.fe[:,ifm]   # Obtain the three edges of the observation triangle
        rm = vtxcrd(ifm, metal) # observation triangle vertex coordinates
        rmc = mean(rm) # observation face centroid (meters)

        for i ∈ 1:length(ξ)   # Loop for numerical integration
            rt = rs[1] + rs21 * ξ[i] + rs31 * η[i] # Source point
            uρ00 = ulocal * (rmc - rt)
            Jsum, Ksum = jksums(uρ00, metal.ψ₁, metal.ψ₂, us1, us2, clsflg) # spatial sums
            # Accumulate results of numerical integration:
            Jsumwght = Jsum * wght[i]
            J_ξ[iufp] += Jsumwght * ξ[i]
            J_η[iufp] += Jsumwght * η[i]
            J[iufp] += Jsumwght
            Ksumwght = Ksum * wght[i]
            K_ξ[iufp] += Ksumwght * ξ[i]
            K_η[iufp] += Ksumwght * η[i]
            K[iufp] += Ksumwght
        end
        if clsflg
            # compute closed-form of singular integrals from Eqs (B.3)
            for i ∈ 1:3  # Loop over sides of source triangle
                ip1 = next[i]  # Next edge in cyclic list.
                # Find components of l, unit tangent vector to edge #i
                lvec = rs[ip1] - rs[i]
                lvec /= norm(lvec)  # Make into unit vector.

                # Find unit outward normal to edge i, lying in plane z = 0
                uvec = -unz * zhatcross(lvec) # so (ux,uy) = (ly * unz, -lx * unz)

                # Compute signed perp. distance from observation point to edge #i
                p0 = uvec ⋅ (rs[i] - rmc)
                p0sq = p0 * p0

                # Compute lplus and lminus as defined in the reference
                lp = lvec ⋅ (rs[ip1] - rmc)
                lm = lvec ⋅ (rs[i] - rmc)
                ledge = lp - lm   # Length of the edge

                pplus = sqrt(p0sq + lp * lp)
                pminus = sqrt(p0sq + lm * lm)

                # Check for special case of observation point on extension of edge #i
                if abs(p0/ledge) < 1e-4
                    p0 = 0.0
                    factor = lp * pplus  -  lm * pminus
                else
                    factor = p0 * log((pplus + lp) / (pminus + lm))
                    rinv[iufp] += factor
                    factor = factor * p0 + pplus * lp - pminus * lm
                end
                ρ_r[iufp] += uvec * factor
            end
            ρ_r[iufp] *= 0.5
            # rinv is now the quantity defined in (B.3a)
            # ρ_r is now the quantity defined in (B.3b)
            ρ_r[iufp] /= area2  # Make it unitless
            rinv[iufp] /= (area2 * ulocal) # Make it unitless
        end

    end

    metal.J = J
    metal.J_ξ = J_ξ
    metal.J_η = J_η
    metal.K = K
    metal.K_ξ = K_ξ
    metal.K_η = K_η
    metal.ρ_r = ρ_r
    metal.rinv = rinv

    return nothing
end

#=
"""
 Return the coordinates (in local units) of the triangle vertices for face iface.
"""
@inline function vtxcrd!(ρvec, iface, metal)
    vi = @view metal.fv[:,iface] # Vertex indices
    ρvec[:] = metal.ρ[vi]
end
=#

"""
 Return the coordinates (in local units) of the triangle vertices for face iface.
"""
@inline function vtxcrd(iface, metal)
    vi = @view metal.fv[:,iface] # Vertex indices
    @view metal.ρ[vi]
end

"""
    zint(Σm1_func, Σm2_func, rs, rmc)  --> (I1, I1_ξ, I1_η, I2)

Compute frequency-dependent integrals needed to fill the generalized impedance matrix.

## Arguments

- `Σm1_func`, `Σm2_func`:  Functions that evaluate the modal series.
- `rs`:  Vector of length three containing the source triangle vertices in meters.
         Each vertex is a 2-vector such as a StaticArrays SV2.
- `rmc`  A 2-vector containing the field observation (match) point in meters.

## Return Values:

- `I1`, `I1_ξ`, `I1_η`, `I2`: Complex, frequency-dependent, spectral integrals defined
                              by Eqs (7.22a), (7.27), and (7.32a).  Their units are (1/m).
"""
@inline function zint(Σm1_func, Σm2_func, rs, rmc)
    ξ = SVector(0.33333333333333330, 0.10128650732345633, 0.79742698535308730,
                0.10128650732345633, 0.47014206410511505, 0.05971587178976989,
                0.47014206410511505)
    η = SVector(0.33333333333333333, 0.79742698535308730, 0.10128650732345633,
                0.10128650732345633, 0.05971587178976989, 0.47014206410511505,
                0.47014206410511505)
    wght = SVector(0.1125, 0.06296959027241358, 0.06296959027241358, 0.06296959027241358,
                   0.06619707639425308, 0.06619707639425308, 0.06619707639425308)

    I1_ξ = zero(ComplexF64)
    I1_η = zero(ComplexF64)
    I1 = zero(ComplexF64)
    I2 = zero(ComplexF64)

    for i ∈ 1:length(ξ)
        rt = rs[1] + (rs[2] - rs[1]) * ξ[i]  +  (rs[3] - rs[1]) * η[i] # Source point
        ρdif = rmc - rt
        sig1 = Σm1_func(ρdif)::ComplexF64
        sig2 = Σm2_func(ρdif)::ComplexF64
        sig1wght = sig1 * wght[i]
        I1_ξ += sig1wght * ξ[i]
        I1_η += sig1wght * η[i]
        I1 += sig1wght
        I2 += sig2 * wght[i]
    end
    return   (I1, I1_ξ, I1_η, I2)
end  # function

end # module


text:
type SnapShot
     nz   :: Int64
     nx   :: Int64
     ext  :: Int64
     iflag:: Int64
     dt   :: Float64
     it   :: Int64
     Vxx  :: Array{Float64,1}
     Vxz  :: Array{Float64,1}
     Vzx  :: Array{Float64,1}
     Vzz  :: Array{Float64,1}
     Txxx :: Array{Float64,1}
     Txxz :: Array{Float64,1}
     Tzzx :: Array{Float64,1}
     Tzzz :: Array{Float64,1}
     Txzx :: Array{Float64,1}
     Txzz :: Array{Float64,1}
end

function InitSnapShot(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dt::Float64, it::Int64)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    spt= SnapShot(nz, nx, ext, iflag, dt, it,
                  zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx),
                  zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx))
    return spt
end

type Wfd
     nz   :: Int64
     nx   :: Int64
     ext  :: Int64
     iflag:: Int64
     dt   :: Float64
     it   :: Int64
     Vx   :: Array{Float64,1}
     Vz   :: Array{Float64,1}
     Txx  :: Array{Float64,1}
     Tzz  :: Array{Float64,1}
     Txz  :: Array{Float64,1}
end

function InitWfd(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dt::Float64, it::Int64)
    wfd= Wfd(nz, nx, ext, iflag, dt, it,
             zeros(nz*nx), zeros(nz*nx), zeros(nz*nx), zeros(nz*nx), zeros(nz*nx))
    return wfd
end

function CopySnapShot!(snapShot1, snapShot2)
    if snapShot2.nz != snapShot1.nz || snapShot2.nx != snapShot1.nx || snapShot2.ext != snapShot1.ext || snapShot2.iflag != snapShot1.iflag
       error("the two snapShot are different")
    end
    snapShot1.it   = snapShot2.it
    n = length(snapShot1.Vxx)
    @inbounds for i = 1 : n
        snapShot1.Vxx[i]  = snapShot2.Vxx[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Vxz[i]  = snapShot2.Vxz[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Vzx[i]  = snapShot2.Vzx[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Vzz[i]  = snapShot2.Vzz[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Txxx[i]  = snapShot2.Txxx[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Txxz[i]  = snapShot2.Txxz[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Tzzx[i]  = snapShot2.Tzzx[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Tzzz[i]  = snapShot2.Tzzz[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Txzx[i]  = snapShot2.Txzx[i]
    end
    @inbounds for i = 1 : n
        snapShot1.Txzz[i]  = snapShot2.Txzz[i]
    end
    return nothing
end

function CopySnapShot(snapShot::SnapShot)
    spt = SnapShot(snapShot.nz, snapShot.nx, snapShot.ext, snapShot.iflag, snapShot.dt, snapShot.it,
                   snapShot.Vxx , snapShot.Vxz , snapShot.Vzx , snapShot.Vzz,
                   snapShot.Txxx, snapShot.Txxz, snapShot.Tzzx, snapShot.Tzzz, snapShot.Txzx, snapShot.Txzz)
    return spt
end

function Add2SnapShots!(snapShot2::SnapShot, snapShot1::SnapShot)
    snapShot2.Vxx = snapShot2.Vxx + snapShot1.Vxx
    snapShot2.Vxz = snapShot2.Vxz + snapShot1.Vxz
    snapShot2.Vzx = snapShot2.Vzx + snapShot1.Vzx
    snapShot2.Vzz = snapShot2.Vzz + snapShot1.Vzz
    snapShot2.Txxx = snapShot2.Txxx + snapShot1.Txxx
    snapShot2.Txxz = snapShot2.Txxz + snapShot1.Txxz
    snapShot2.Tzzx = snapShot2.Tzzx + snapShot1.Tzzx
    snapShot2.Tzzz = snapShot2.Tzzz + snapShot1.Tzzz
    snapShot2.Txzx = snapShot2.Txzx + snapShot1.Txzx
    snapShot2.Txzz = snapShot2.Txzz + snapShot1.Txzz
    return nothing
end

function WriteSnapShot(pathout::String, snapShot::SnapShot)
    fid = open(pathout, "w")
    write(fid, convert(Int32  , snapShot.nz))
    write(fid, convert(Int32  , snapShot.nx))
    write(fid, convert(Int32  , snapShot.ext))
    write(fid, convert(Int32  , snapShot.iflag))
    write(fid, convert(Float32, snapShot.dt))
    write(fid, convert(Array{Float32}, snapShot.Vxx))
    write(fid, convert(Array{Float32}, snapShot.Vxz))
    write(fid, convert(Array{Float32}, snapShot.Vzx))
    write(fid, convert(Array{Float32}, snapShot.Vzz))
    write(fid, convert(Array{Float32}, snapShot.Txxx))
    write(fid, convert(Array{Float32}, snapShot.Txxz))
    write(fid, convert(Array{Float32}, snapShot.Tzzx))
    write(fid, convert(Array{Float32}, snapShot.Tzzz))
    write(fid, convert(Array{Float32}, snapShot.Txzx))
    write(fid, convert(Array{Float32}, snapShot.Txzz))
    return fid
end

function WriteSnapShot(fid::IOStream, snapShot::SnapShot)
    write(fid, convert(Array{Float32}, snapShot.Vxx))
    write(fid, convert(Array{Float32}, snapShot.Vxz))
    write(fid, convert(Array{Float32}, snapShot.Vzx))
    write(fid, convert(Array{Float32}, snapShot.Vzz))
    write(fid, convert(Array{Float32}, snapShot.Txxx))
    write(fid, convert(Array{Float32}, snapShot.Txxz))
    write(fid, convert(Array{Float32}, snapShot.Tzzx))
    write(fid, convert(Array{Float32}, snapShot.Tzzz))
    write(fid, convert(Array{Float32}, snapShot.Txzx))
    write(fid, convert(Array{Float32}, snapShot.Txzz))
    flush(fid)
    return nothing
end

function ReverseOrderSnapShots(pathin::String)
    fid = open(pathin, "r")
    pathin_tmp = join([pathin "_tmp"])
    fid1  = open(pathin_tmp, "w")
    nz    = read(fid, Int32)
    nx    = read(fid, Int32)
    ext   = read(fid, Int32)
    iflag = read(fid, Int32)
    dt    = read(fid, Float32)
    write(fid1, nz, nx, ext, iflag, dt);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    snapSize = sizeof(Float32)*Nz*Nx*10
    nt = round(Int64, ((filesize(fid)-sizeof(Float32)*5)/snapSize))
    for it = nt : -1 : 1
        position = sizeof(Float32)*5 + (it-1)*snapSize
        seek(fid, position)
        tmp = read(fid, Float32, Nz*Nx*10)
        write(fid1, tmp)
    end
    close(fid); close(fid1)
    mv(pathin_tmp, pathin, remove_destination=true)
    return nothing
end

function WriteWfd(pathout::String, snapShot::SnapShot)
    fid = open(pathout, "w")
    nz = snapShot.nz;  nx = snapShot.nx;
    ext= snapShot.ext; iflag = snapShot.iflag;
    write(fid, convert(Int32,  nz))
    write(fid, convert(Int32,  nx))
    write(fid, convert(Int32,  ext))
    write(fid, convert(Int32,  iflag))
    write(fid, convert(Float32, snapShot.dt))
    if iflag == 1
       Nz = nz + 2*ext
       zupper = ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    tmp = snapShot.Vxx+snapShot.Vxz;   tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Vzx+snapShot.Vzz;   tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Txxx+snapShot.Txxz; tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Tzzx+snapShot.Tzzz; tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Txzx+snapShot.Txzz; tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    return fid
end

function WriteWfd(fid::IOStream, snapShot::SnapShot)
    nz = snapShot.nz;  nx = snapShot.nx;
    ext= snapShot.ext; iflag = snapShot.iflag;
    if iflag == 1
       Nz = nz + 2*ext
       zupper = ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    tmp = snapShot.Vxx+snapShot.Vxz;   tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Vzx+snapShot.Vzz;   tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Txxx+snapShot.Txxz; tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Tzzx+snapShot.Tzzz; tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    tmp = snapShot.Txzx+snapShot.Txzz; tmp = reshape(tmp, Nz, Nx); write(fid, convert(Array{Float32}, vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])))
    flush(fid)
    return nothing
end

function extractBoundary!(bv::Array{Float64,1}, spt::SnapShot)
    nz = spt.nz;  nx = spt.nx;
    ext= spt.ext; iflag = spt.iflag;
    if iflag == 1
       Nz = nz + 2*ext
       zu = ext
    elseif iflag == 2
       Nz = nz +   ext
       zu = 0
    end
    Nx = nx + 2*ext
    # save vz
    tmp=spt.Vzx+spt.Vzz; tmp=reshape(tmp,Nz,Nx);
    tmp=tmp[zu+1:zu+nz,ext+1:ext+nx];
    il=   1; iu=   2*nz    ; bv[il:iu]=vec(tmp[:,1:2]);
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,end-1:end]);
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[1:2,3:nx-2])');
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[end-1:end,3:nx-2])');
    # save vx
    tmp=spt.Vxz+spt.Vxx; tmp=reshape(tmp,Nz,Nx);
    tmp=tmp[zu+1:zu+nz,ext+1:ext+nx];
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,1:2]);
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,end-1:end]);
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[1:2,3:nx-2])');
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[end-1:end,3:nx-2])');
    # save Tzz
    tmp=spt.Tzzz+spt.Tzzx; tmp=reshape(tmp,Nz,Nx);
    tmp=tmp[zu+1:zu+nz,ext+1:ext+nx];
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,1:2]);
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,end-1:end]);
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[1:2,3:nx-2])');
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[end-1:end,3:nx-2])');
    # save Txx
    tmp=spt.Txxz+spt.Txxx; tmp=reshape(tmp,Nz,Nx);
    tmp=tmp[zu+1:zu+nz,ext+1:ext+nx];
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,1:2]);
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,end-1:end]);
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[1:2,3:nx-2])');
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[end-1:end,3:nx-2])');
    # save Txz
    tmp=spt.Txzx+spt.Txzz; tmp=reshape(tmp,Nz,Nx);
    tmp=tmp[zu+1:zu+nz,ext+1:ext+nx];
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,1:2]);
    il=iu+1; iu=iu+2*nz    ; bv[il:iu]=vec(tmp[:,end-1:end]);
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[1:2,3:nx-2])');
    il=iu+1; iu=iu+2*(nx-4); bv[il:iu]=vec((tmp[end-1:end,3:nx-2])');
    return nothing
end

function WriteWfdBoundary(path::String, spt::SnapShot, bv::Array{Float64,1})
    fid = open(path, "w")
    nz = spt.nz;  nx = spt.nx;
    ext= spt.ext; iflag = spt.iflag;
    write(fid, nz, nx)
    extractBoundary!(bv, spt)
    write(fid, bv)
    flush(fid)
    return fid
end

function WriteWfdBoundary(fid::IOStream, spt::SnapShot, bv::Array{Float64,1})
    extractBoundary!(bv, spt)
    write(fid, bv)
    flush(fid)
    return fid
end

function ReadWfdBoundary(path::String)
    fid = open(path, "r")
    nz=read(fid,Int64); nx=read(fid,Int64);
    l = (4*nz+(nx-4)*4)*5
    nt= floor(Int64, (filesize(fid)-sizeof(Float64)*2)/(l*sizeof(Float64)))
    bv= read(fid, Float64, nt*l)
    bv= reshape(bv, l, nt)
    return bv
end

function ReadWfdBoundary(path::String, it::Int64)
    fid = open(path, "r")
    nz=read(fid,Int64); nx=read(fid,Int64);
    l = (4*nz+(nx-4)*4)*5
    position = sizeof(Int64) + (it-1)*l*sizeof(Float64)
    seek(fid, position)
    bv = read(fid, Float64, l)
    return bv
end

function ReadWfd(pathin::String, it::Int64)
    fid = open(pathin, "r")
    nz  = convert(Int64, read(fid, Int32))
    nx  = convert(Int64, read(fid, Int32))
    ext = convert(Int64, read(fid, Int32))
    iflag = convert(Int64, read(fid, Int32))
    dt  = convert(Float64, read(fid, Float32))
    wfd = InitWfd(nz, nx, ext, iflag, dt, it)
    wfdSize = sizeof(Float32)*nz*nx*5
    position = sizeof(Float32)*5 + (it-1)*wfdSize
    seek(fid, position)
    wfd.Vx  = convert(Array{Float64}, read(fid,Float32, nz*nx))
    wfd.Vz  = convert(Array{Float64}, read(fid,Float32, nz*nx))
    wfd.Txx = convert(Array{Float64}, read(fid,Float32, nz*nx))
    wfd.Tzz = convert(Array{Float64}, read(fid,Float32, nz*nx))
    wfd.Txz = convert(Array{Float64}, read(fid,Float32, nz*nx))
    close(fid)
    return wfd
end

function InfoWfd(pathin::String; flag=false)
    fid = open(pathin, "r")
    nz  = convert(Int64, read(fid, Int32))
    nx  = convert(Int64, read(fid, Int32))
    ext = convert(Int64, read(fid, Int32))
    iflag = convert(Int64, read(fid, Int32))
    dt  = convert(Float64, read(fid, Float32))
    wfdSize = sizeof(Float32)*nz*nx*5
    nt = round(Int64, (filesize(fid)-sizeof(Int32)*5)/wfdSize)
    if flag == true
       println("nz: $nz")
       println("nx: $nx")
       println("ext: $ext")
       println("iflag: $iflag")
       println("dt: $dt")
       println("nt: $nt")
       tmax = (nt-1)*dt
       println("tmax: $tmax")
       fSize = filesize(fid)/1024/1024/1024
       println("file size: $fSize G")
    end
    close(fid)
    return nz, nx, ext, iflag, dt, nt
end

function InfoSnapShots(pathin::String; flag=false)
    fid = open(pathin, "r")
    nz  = convert(Int64, read(fid, Int32))
    nx  = convert(Int64, read(fid, Int32))
    ext = convert(Int64, read(fid, Int32))
    iflag = convert(Int64, read(fid, Int32))
    dt  = convert(Float64, read(fid, Float32))
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    snapSize = sizeof(Float32)*Nz*Nx*10
    nt = round(Int64, (filesize(fid)-sizeof(Int32)*5)/snapSize)
    if flag == true
       println("nz: $nz")
       println("nx: $nx")
       println("ext: $ext")
       println("iflag: $iflag")
       println("dt: $dt")
       println("nt: $nt")
       tmax = (nt-1)*dt
       println("tmax: $tmax")
       fSize = filesize(fid)/1024/1024/1024
       println("file size: $fSize G")
    end
    close(fid)
    return nz, nx, ext, iflag, dt, nt
end

function ReadSnapShot(pathin::String, it::Int64)
    fid = open(pathin, "r")
    nz  = convert(Int64, read(fid, Int32))
    nx  = convert(Int64, read(fid, Int32))
    ext = convert(Int64, read(fid, Int32))
    iflag = convert(Int64, read(fid, Int32))
    dt  = convert(Float64, read(fid, Float32))
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    snapShot = InitSnapShot(nz, nx, ext, iflag, dt, it)
    snapSize = sizeof(Float32)*Nz*Nx*10
    position = sizeof(Float32)*5 + (it-1)*snapSize
    seek(fid, position)
    snapShot.Vxx  = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Vxz  = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Vzx  = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Vzz  = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Txxx = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Txxz = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Tzzx = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Tzzz = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Txzx = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    snapShot.Txzz = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    close(fid)
    return snapShot
end

function ExtractCpt(path::String, pathin::String; itype="wfd", cpt="Vx")
    fid = open(pathin, "r")
    nz    = read(fid, Int32)
    nx    = read(fid, Int32)
    ext   = read(fid, Int32)
    iflag = read(fid, Int32)
    dt    = read(fid, Float32)
    fido = open(path, "w")
    write(fido, nz, nx, ext, iflag, dt)
    if itype == "spt"
       if iflag == 1
          Nz = nz + 2*ext
          upper = ext
       elseif iflag == 2
          Nz = nz +   ext
          upper = 0
       end
       Nx = nx + 2*ext
       sptSize = sizeof(Float32)*Nz*Nx*10
       nt = round(Int64, ((filesize(fid)-sizeof(Float32)*5)/sptSize))
       inds = 0
       if cpt == "Vx"
          inds = sizeof(Float32)*5
       elseif cpt == "Vz"
          inds = sizeof(Float32)*5 + 2*Nz*Nx*sizeof(Float32)
       elseif cpt == "Txx"
          inds = sizeof(Float32)*5 + 4*Nz*Nx*sizeof(Float32)
       elseif cpt == "Txx"
          inds = sizeof(Float32)*5 + 6*Nz*Nx*sizeof(Float32)
       elseif cpt == "Txx"
          inds = sizeof(Float32)*5 + 8*Nz*Nx*sizeof(Float32)
       end
       for it = 1: nt
           position = (it-1)*sptSize + inds
           seek(fid, position)
           tmp = read(fid,Float32,Nz*Nx) + read(fid,Float32,Nz*Nx)
           tmp = reshape(tmp, Int64(Nz), Int64(Nx))
           write(fido, vec(tmp[upper+1:upper+nz, ext+1:ext+nx]))
       end
    elseif itype == "wfd"
       sptSize = sizeof(Float32)*nz*nx*5
       nt = round(Int64, ((filesize(fid)-sizeof(Float32)*5)/sptSize))
       inds = 0
       if cpt == "Vx"
         inds = sizeof(Float32)*5
       elseif cpt == "Vz"
         inds = sizeof(Float32)*5 +   nz*nx*sizeof(Float32)
       elseif cpt == "Txx"
         inds = sizeof(Float32)*5 + 2*nz*nx*sizeof(Float32)
       elseif cpt == "Txx"
         inds = sizeof(Float32)*5 + 2*nz*nx*sizeof(Float32)
       elseif cpt == "Txx"
         inds = sizeof(Float32)*5 + 2*nz*nx*sizeof(Float32)
       end
       for it = 1 : nt
           position = (it-1)*sptSize + inds
           seek(fid, position)
           tmp = read(fid, Float32, nz*nx)
           write(fido, tmp)
       end
    end
    close(fido)
    return nothing
end

function NormOfSnapShots(pathin::String, itype="l2")
    fid = open(pathin, "r")
    nz  = read(fid, Int32)
    nx  = read(fid, Int32)
    ext = read(fid, Int32)
    iflag = read(fid, Int32)
    dt    = read(fid, Float32)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    snapSize = sizeof(Float32)*Nz*Nx*10
    nt = round(Int64, ((filesize(fid)-sizeof(Float32)*5)/snapSize))
    if itype == "l2"
       tmp = 0.0
       for it = 1 : nt
           val = convert(Array{Float64}, read(fid, Float32, Nz*Nx*10))
           tmp = tmp + dot(val, val)
       end
       tmp = sqrt(tmp)
    elseif itype == "l1"
       tmp = 0.0
       for it = 1 : nt
           val = convert(Array{Float64}, read(fid, Float32, Nz*Nx*10))
           tmp = tmp + sum(abs(val))
       end
    end
    close(fid)
    return tmp
end

function NormOfSpt(snapShot::SnapShot, itype="l2")
  tmp = 0.0
  if itype == "l2"
     tmp = tmp + dot(snapShot.Vxx, snapShot.Vxx)
     tmp = tmp + dot(snapShot.Vxz, snapShot.Vxz)
     tmp = tmp + dot(snapShot.Vzx, snapShot.Vzx)
     tmp = tmp + dot(snapShot.Vzz, snapShot.Vzz)
     tmp = tmp + dot(snapShot.Txxx, snapShot.Txxx)
     tmp = tmp + dot(snapShot.Txxz, snapShot.Txxz)
     tmp = tmp + dot(snapShot.Tzzx, snapShot.Tzzx)
     tmp = tmp + dot(snapShot.Tzzz, snapShot.Tzzz)
     tmp = tmp + dot(snapShot.Txzx, snapShot.Txzx)
     tmp = tmp + dot(snapShot.Txzz, snapShot.Txzz)
     tmp = sqrt(tmp)
  elseif itype == "l1"
     tmp = tmp + sum(abs(snapShot.Vxx))
     tmp = tmp + sum(abs(snapShot.Vxz))
     tmp = tmp + sum(abs(snapShot.Vzx))
     tmp = tmp + sum(abs(snapShot.Vzz))
     tmp = tmp + sum(abs(snapShot.Txxx))
     tmp = tmp + sum(abs(snapShot.Txxz))
     tmp = tmp + sum(abs(snapShot.Tzzx))
     tmp = tmp + sum(abs(snapShot.Tzzz))
     tmp = tmp + sum(abs(snapShot.Txzx))
     tmp = tmp + sum(abs(snapShot.Txzz))
  end
  return tmp
end

function ScaleWfd!(pathin::String, alpha::Float64)
    alpha = convert(Float32, alpha)
    fid = open(pathin, "r+")
    nz    = read(fid, Int32)
    nx    = read(fid, Int32)
    snapSize = sizeof(Float32)*nz*nx*5
    nt = round(Int64, ((filesize(fid)-sizeof(Float32)*5)/snapSize))
    for it =1: nt
        position = sizeof(Float32)*5 + (it-1)*snapSize
        seek(fid, position)
        tmp = read(fid, Float32, nz*nx*5) .* alpha
        seek(fid, position)
        write(fid, tmp)
    end
    close(fid)
    return nothing
end

function ReverseOrderWfd(pathin::String)
    fid = open(pathin, "r")
    pathin_tmp = join([pathin "_tmp"])
    fid1  = open(pathin_tmp, "w")
    nz    = read(fid, Int32)
    nx    = read(fid, Int32)
    ext   = read(fid, Int32)
    iflag = read(fid, Int32)
    dt    = read(fid, Float32)
    write(fid1, nz, nx, ext, iflag, dt);
    snapSize = sizeof(Float32)*nz*nx*5
    nt = round(Int32, ((filesize(fid)-sizeof(Float32)*5)/snapSize))
    for it = nt : -1 : 1
        position = sizeof(Float32)*5 + (it-1)*snapSize
        seek(fid, position)
        tmp = read(fid, Float32, nz*nx*5)
        write(fid1, tmp)
    end
    close(fid); close(fid1)
    mv(pathin_tmp, pathin, remove_destination=true)
    return nothing
end

function IptWfd(path::String, path1::String)
    fid  = open(path , "r")
    fid1 = open(path1, "r")
    nz =read(fid ,Int32); nx =read(fid ,Int32); ext =read(fid ,Int32); iflag =read(fid ,Int32);
    nz1=read(fid1,Int32); nx1=read(fid1,Int32); ext1=read(fid1,Int32); iflag1=read(fid1,Int32);
    if nz1!=nz || nx1!=nx || ext1!=ext || iflag1!=iflag
       error("size dismatch")
    end
    snapSize = sizeof(Float32)*nz*nx*5
    nt = round(Int32, ((filesize(fid)-sizeof(Float32)*5)/snapSize))
    seek(fid , sizeof(Float32)*5)
    seek(fid1, sizeof(Float32)*5)
    tmp = 0.0
    for it = 1 : nt
        tmp1 = read(fid , Float32, nz*nx*5)
        tmp2 = read(fid1, Float32, nz*nx*5)
        tmp  = tmp + dot(tmp, tmp1)
    end
    close(fid); clsoe(fid1)
    return tmp
end

function partialV!(vxx::Array{Float64,2}, vxz::Array{Float64,2}, vzx::Array{Float64,2}, vzz::Array{Float64,2}, spt::SnapShot, dz::Float64, dx::Float64, vx::Array{Float64,1}, vz::Array{Float64,1})
    a1 = 9/8; a2 = -1/24;
    nz = spt.nz; nx = spt.nx; ext = spt.ext; iflag = spt.iflag
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    mysum1!(vx, spt.Vxx, spt.Vxz)
    mysum1!(vz, spt.Vzx, spt.Vzz)
    for ix = 1 : Nx
        for iz = 1 : Nz
          # compute Dvx/Dx, centered on (i,j)
            indp1h = (ix-3)*Nz + iz
            indph  = (ix-2)*Nz + iz
            indh   = (ix-1)*Nz + iz
            ind1h  = (ix  )*Nz + iz
            if ix == 1
               vxx[iz,ix] = a1/dx*(vx[indh]          ) + a2/dx*(vx[ind1h]           )
            elseif ix == 2
               vxx[iz,ix] = a1/dx*(vx[indh]-vx[indph]) + a2/dx*(vx[ind1h]           )
            elseif ix == Nx
               vxx[iz,ix] = a1/dx*(vx[indh]-vx[indph]) + a2/dx*(         -vx[indp1h])
            else
               vxx[iz,ix] = a1/dx*(vx[indh]-vx[indph]) + a2/dx*(vx[ind1h]-vx[indp1h])
            end

          # compute Dvx/Dz, centered on (i+1/2, j+1/2)
            indp1h = (ix-1)*Nz + iz-1
            indph  = (ix-1)*Nz + iz
            indh   = (ix-1)*Nz + iz+1
            ind1h  = (ix-1)*Nz + iz+2
            if iz == 1
               vxz[iz,ix] = a1/dz*(vx[indh]-vx[indph]) + a2/dz*(vx[ind1h]           )
            elseif iz == Nz-1
               vxz[iz,ix] = a1/dz*(vx[indh]-vx[indph]) + a2/dz*(         -vx[indp1h])
            elseif iz == Nz
               vxz[iz,ix] = a1/dz*(        -vx[indph]) + a2/dz*(         -vx[indp1h])
            else
               vxz[iz,ix] = a1/dz*(vx[indh]-vx[indph]) + a2/dz*(vx[ind1h]-vx[indp1h])
            end

            # compute Dvz/Dx
            indp1h = (ix-2)*Nz + iz
            indph  = (ix-1)*Nz + iz
            indh   = (ix  )*Nz + iz
            ind1h  = (ix+1)*Nz + iz
            if ix == 1
               vzx[iz,ix] = a1/dx*(vz[indh]-vz[indph]) + a2/dx*(vz[ind1h]           )
            elseif ix == Nx-1
               vzx[iz,ix] = a1/dx*(vz[indh]-vz[indph]) + a2/dx*(         -vz[indp1h])
            elseif ix == Nx
               vzx[iz,ix] = a1/dx*(        -vz[indph]) + a2/dx*(         -vz[indp1h])
            else
               vzx[iz,ix] = a1/dx*(vz[indh]-vz[indph]) + a2/dx*(vz[ind1h]-vz[indp1h])
            end

          # compute Dvz/Dz, centered on (i+1/2, j+1/2)
            indp1h = (ix-1)*Nz + iz-2
            indph  = (ix-1)*Nz + iz-1
            indh   = (ix-1)*Nz + iz
            ind1h  = (ix-1)*Nz + iz+1
            if iz == 1
               vzz[iz,ix] = a1/dz*(vz[indh]          ) + a2/dz*(vz[ind1h]           )
            elseif iz == 2
               vzz[iz,ix] = a1/dz*(vz[indh]-vz[indph]) + a2/dz*(vz[ind1h]           )
            elseif iz == Nz
               vzz[iz,ix] = a1/dz*(vz[indh]-vz[indph]) + a2/dz*(         -vz[indp1h])
            else
               vzz[iz,ix] = a1/dz*(vz[indh]-vz[indph]) + a2/dz*(vz[ind1h]-vz[indp1h])
            end
        end
    end
    return nothing
end

# function partialV!(vxx::Array{Float64,2}, vxz::Array{Float64,2}, vzx::Array{Float64,2}, vzz::Array{Float64,2}, spt::SnapShot, dz::Float64, dx::Float64, vx::Array{Float64,1}, vz::Array{Float64,1})
#     a1 = 9/8; a2 = -1/24;
#     nz = spt.nz; nx = spt.nx; ext = spt.ext; iflag = spt.iflag
#     if iflag == 1
#        Nz = nz + 2*ext
#        upper = ext
#     elseif iflag == 2
#        Nz = nz +   ext
#        upper = 0
#     end
#     mysum1!(vx, spt.Vxx, spt.Vxz)
#     mysum1!(vz, spt.Vzx, spt.Vzz)
#     for iz = 1 : nz
#         for ix = 1 : nx
#           # compute Dvx/Dx, centered on (i,j)
#             ind1 = (ix+ext-3)*Nz + (iz+upper)
#             ind2 = (ix+ext-2)*Nz + (iz+upper)
#             ind3 = (ix+ext-1)*Nz + (iz+upper)
#             ind4 = (ix+ext  )*Nz + (iz+upper)
#             vxx[iz,ix] =  a1/dx*(vx[ind3]-vx[ind2]) + a2/dx*(vx[ind4]-vx[ind1])
#           # compute Dvx/Dz
#             if iflag == 2 && iz == 1
#                ind1 = (ix+ext-1)*Nz + iz+upper
#                ind2 = (ix+ext-1)*Nz + iz+upper+1
#                ind3 = (ix+ext-1)*Nz + iz+upper+2
#                vxz[iz,ix] =  a1/dz*(vx[ind2]-vx[ind1]) + a2/dx*(vx[ind3])
#             else
#                ind1 = (ix+ext-1)*Nz + iz+upper-1
#                ind2 = (ix+ext-1)*Nz + iz+upper
#                ind3 = (ix+ext-1)*Nz + iz+upper+1
#                ind4 = (ix+ext-1)*Nz + iz+upper+2
#                vxz[iz,ix] =  a1/dz*(vx[ind3]-vx[ind2]) + a2/dz*(vx[ind4]-vx[ind1])
#             end
#             # compute Dvz/Dx
#                ind1 = (ix+ext-2)*Nz + (iz+upper)
#                ind2 = (ix+ext-1)*Nz + (iz+upper)
#                ind3 = (ix+ext  )*Nz + (iz+upper)
#                ind4 = (ix+ext+1)*Nz + (iz+upper)
#                vzx[iz,ix] =  a1/dx*(vz[ind3]-vz[ind2]) + a2/dx*(vz[ind4]-vz[ind1])
#           # compute Dvz/Dz
#             if iflag == 2 && iz == 1
#                ind1 = (ix+ext-1)*Nz + iz+upper
#                ind2 = (ix+ext-1)*Nz + iz+upper+1
#                vzz[iz,ix] =  a1/dz*vz[ind1] + a2/dz*vz[ind2]
#             elseif iflag == 2 && iz == 2
#                ind1 = (ix+ext-1)*Nz + iz+upper-1
#                ind2 = (ix+ext-1)*Nz + iz+upper
#                ind3 = (ix+ext-1)*Nz + iz+upper+1
#                vzz[iz,ix] =  a1/dz*(vz[ind2]-vz[ind1]) + a2/dz*vz[ind3]
#             else
#                ind1 = (ix+ext-1)*Nz + iz+upper-2
#                ind2 = (ix+ext-1)*Nz + iz+upper-1
#                ind3 = (ix+ext-1)*Nz + iz+upper
#                ind4 = (ix+ext-1)*Nz + iz+upper+1
#                vzz[iz,ix] =  a1/dz*(vz[ind3]-vz[ind2]) + a2/dz*(vz[ind4]-vz[ind1])
#             end
#         end
#     end
#     return nothing
# end

function WritePv(path::String, vxx::Array{Float64,2}, vxz::Array{Float64,2}, vzx::Array{Float64,2}, vzz::Array{Float64,2}, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dt::Float64)
    fid = open(path, "w")
    write(fid, Int32(nz), Int32(nx), Int32(ext), Int32(iflag), Float32(dt))
    write(fid, convert(Array{Float32}, vec(vxx)))
    write(fid, convert(Array{Float32}, vec(vxz)))
    write(fid, convert(Array{Float32}, vec(vzx)))
    write(fid, convert(Array{Float32}, vec(vzz)))
    flush(fid)
    return fid
end

function WritePv(fid::IOStream, vxx::Array{Float64,2}, vxz::Array{Float64,2}, vzx::Array{Float64,2}, vzz::Array{Float64,2})
    write(fid, convert(Array{Float32}, vec(vxx)))
    write(fid, convert(Array{Float32}, vec(vxz)))
    write(fid, convert(Array{Float32}, vec(vzx)))
    write(fid, convert(Array{Float32}, vec(vzz)))
    flush(fid)
    return nothing
end

function ReadPv(path::String, it::Int64)
    fid = open(path, "r")
    nz  = Int64(read(fid,Int32)); nx = Int64(read(fid,Int32));
    ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
    dt  = Float64(read(fid,Float32));
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    sliceSize = sizeof(Float32) * Nz * Nx * 4
    position  = sizeof(Int32)*5 + (it-1)*sliceSize
    seek(fid, position)
    vxx = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    vxz = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    vzx = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    vzz = convert(Array{Float64}, read(fid, Float32, Nz*Nx))
    close(fid)
    return vxx, vxz, vzx, vzz
end

# function ReadPv(path::String, it::Int64)
#     fid = open(path, "r")
#     nz  = Int64(read(fid,Int32)); nx = Int64(read(fid,Int32));
#     ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
#     dt  = Float64(read(fid,Float32));
#     sliceSize = sizeof(Float32) * nz * nx * 4
#     position  = sizeof(Int32)*5 + (it-1)*sliceSize
#     seek(fid, position)
#     vxx = convert(Array{Float64}, read(fid, Float32, nz*nx))
#     vxz = convert(Array{Float64}, read(fid, Float32, nz*nx))
#     vzx = convert(Array{Float64}, read(fid, Float32, nz*nx))
#     vzz = convert(Array{Float64}, read(fid, Float32, nz*nx))
#     vxx = reshape(vxx, nz, nx)
#     vxz = reshape(vxz, nz, nx)
#     vzx = reshape(vzx, nz, nx)
#     vzz = reshape(vzz, nz, nx)
#     close(fid)
#     return vxx, vxz, vzx, vzz
# end

function InfoPv(path::String; print_flag=false)
    fid = open(path, "r")
    nz  = Int64(read(fid,Int32)); nx = Int64(read(fid,Int32));
    ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
    dt = Float64(read(fid, Float32));
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    sliceSize = sizeof(Float32) * Nz * Nx * 4
    ns = round(Int64, (filesize(fid)-sizeof(Int32)*4)/sliceSize)
    if print_flag
       println("nz   : $nz"   )
       println("nx   : $nx"   )
       println("ext  : $ext"  )
       println("iflag: $iflag")
       println("dt   : $dt"   )
       println("ns   : $ns"   )
    end
    close(fid)
    return nz, nx, ext, iflag, dt, ns
end
# function InfoPv(path::String; print_flag=false)
#     fid = open(path, "r")
#     nz  = Int64(read(fid,Int32)); nx = Int64(read(fid,Int32));
#     ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
#     dt = Float64(read(fid, Float32));
#     sliceSize = sizeof(Float32) * nz * nx * 4
#     ns = round(Int64, (filesize(fid)-sizeof(Int32)*4)/sliceSize)
#     if print_flag
#        println("nz   : $nz"   )
#        println("nx   : $nx"   )
#        println("ext  : $ext"  )
#        println("iflag: $iflag")
#        println("dt   : $dt"   )
#        println("ns   : $ns"   )
#     end
#     close(fid)
#     return nz, nx, ext, iflag, dt, ns
# end

function checkPv(path, path1)
    (nz, nx, ext, iflag, dt, nt) = InfoPv(path)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    vxx = zeros(Nz,Nx); vxz = zeros(Nz,Nx);
    vzx = zeros(Nz,Nx); vzz = zeros(Nz,Nx);
    vxx1 = zeros(Nz,Nx); vxz1 = zeros(Nz,Nx);
    vzx1 = zeros(Nz,Nx); vzz1 = zeros(Nz,Nx);
    d=0.
    for it = 1 : nt
        (vxx , vxz , vzx , vzz ) = ReadPv(path, it)
        (vxx1, vxz1, vzx1, vzz1) = ReadPv(path1, it)
        d = d + norm(vxx-vxx1)
        d = d + norm(vxz-vxz1)
        d = d + norm(vzx-vzx1)
        d = d + norm(vzz-vzz1)
    end
    return d
end

# function checkPv(path, path1)
#     (nz, nx, ext, iflag, dt, nt) = InfoPv(path)
#     vxx = zeros(nz,nx); vxz = zeros(nz,nx);
#     vzx = zeros(nz,nx); vzz = zeros(nz,nx);
#     vxx1 = zeros(nz,nx); vxz1 = zeros(nz,nx);
#     vzx1 = zeros(nz,nx); vzz1 = zeros(nz,nx);
#     d=0.
#     for it = 1 : nt
#         (vxx , vxz , vzx , vzz ) = ReadPv(path, it)
#         (vxx1, vxz1, vzx1, vzz1) = ReadPv(path1, it)
#         d = d + vecnorm(vxx-vxx1)
#         d = d + vecnorm(vxz-vxz1)
#         d = d + vecnorm(vzx-vzx1)
#         d = d + vecnorm(vzz-vzz1)
#     end
#     return d
# end

# function IptSpts(path::String, path1::String)
#     (nz, nx, ext, iflag, dt, nt) = InfoSnapShots(path1)
#     tmp = 0.0
#     for it = 1 : nt
#         spt = ReadSnapShot(path , it)
#         spt1= ReadSnapShot(path1, it)
#         tmp = tmp + dot(spt.Vxx , spt1.Vxx )
#         tmp = tmp + dot(spt.Vxz , spt1.Vxz )
#         tmp = tmp + dot(spt.Vzx , spt1.Vzx )
#         tmp = tmp + dot(spt.Vzz , spt1.Vzz )
#         tmp = tmp + dot(spt.Txxx, spt1.Txxx)
#         tmp = tmp + dot(spt.Txxz, spt1.Txxz)
#         tmp = tmp + dot(spt.Tzzx, spt1.Tzzx)
#         tmp = tmp + dot(spt.Tzzz, spt1.Tzzz)
#         tmp = tmp + dot(spt.Txzx, spt1.Txzx)
#         tmp = tmp + dot(spt.Txzz, spt1.Txzz)
#     end
#     return tmp
# end

function IptSpts(spt1::SnapShot, spt2::SnapShot)
    d = 0.0
    d = d + dot(spt1.Vxx, spt2.Vxx)
    d = d + dot(spt1.Vxz, spt2.Vxz)
    d = d + dot(spt1.Vzx, spt2.Vzx)
    d = d + dot(spt1.Vzz, spt2.Vzz)
    d = d + dot(spt1.Txxx, spt2.Txxx)
    d = d + dot(spt1.Txxz, spt2.Txxz)
    d = d + dot(spt1.Tzzx, spt2.Tzzx)
    d = d + dot(spt1.Tzzz, spt2.Tzzz)
    d = d + dot(spt1.Txzx, spt2.Txzx)
    d = d + dot(spt1.Txzz, spt2.Txzz)
    return d
end

function dif2spt(spt1::SnapShot, spt2::SnapShot)
    d = 0.0
    d = d + norm(spt1.Vxx-spt2.Vxx)
    d = d + norm(spt1.Vxz-spt2.Vxz)
    d = d + norm(spt1.Vzx-spt2.Vzx)
    d = d + norm(spt1.Vzz-spt2.Vzz)
    d = d + norm(spt1.Txxx-spt2.Txxx)
    d = d + norm(spt1.Txxz-spt2.Txxz)
    d = d + norm(spt1.Tzzx-spt2.Tzzx)
    d = d + norm(spt1.Tzzz-spt2.Tzzz)
    d = d + norm(spt1.Txzx-spt2.Txzx)
    d = d + norm(spt1.Txzz-spt2.Txzz)
    return d
end

ext: jl

max_stars_repo_name: WneleiGao/ElasticWave


text:
# TODO: allow passing a rng as first parameter

"""
    randobs(data, [n])
Pick a random observation or a batch of `n` random observations
from `data`.
For this function to work, the type of `data` must implement
[`numobs`](@ref) and [`getobs`](@ref).
"""
randobs(data) = getobs(data, rand(1:numobs(data)))

randobs(data, n) = getobs(data, rand(1:numobs(data), n))


text:
using Test
using Base64, JSON

import IJulia
import IJulia: helpmode, error_content, docdict

@testset "errors" begin
    content = error_content(UndefVarError(:a))
    @test "UndefVarError" == content["ename"]
end

@testset "docdict" begin
    @test haskey(docdict("import"), "text/plain")
    @test haskey(docdict("sum"), "text/plain")
end

struct FriendlyData
    name::AbstractString
end

@testset "Custom MIME types" begin
    friend = FriendlyData("world")

    FRIENDLY_MIME_TYPE = MIME"application/vnd.ijulia.friendly-text"
    FRIENDLY_MIME = FRIENDLY_MIME_TYPE()
    Base.Multimedia.istextmime(::FRIENDLY_MIME_TYPE) = true
    Base.show(io, ::FRIENDLY_MIME_TYPE, x::FriendlyData) = write(io, "Hello, $(x.name)!")
    IJulia.register_mime(FRIENDLY_MIME)

    BINARY_MIME_TYPE = MIME"application/vnd.ijulia.friendly-binary"
    BINARY_MIME = BINARY_MIME_TYPE()
    Base.Multimedia.istextmime(::BINARY_MIME_TYPE) = false
    Base.show(io, ::BINARY_MIME_TYPE, x::FriendlyData) = write(io, "Hello, $(x.name)!")
    IJulia.register_mime(BINARY_MIME)

    JSON_MIME_TYPE = MIME"application/vnd.ijulia.friendly-json"
    JSON_MIME = JSON_MIME_TYPE()
    Base.Multimedia.istextmime(::JSON_MIME_TYPE) = true
    Base.show(io, ::JSON_MIME_TYPE, x::FriendlyData) = write(io, JSON.json(Dict("name" => x.name)))
    IJulia.register_jsonmime(JSON_MIME)

    FRIENDLY_MIME_TYPE_1 = MIME"application/vnd.ijulia.friendly-text-1"
    FRIENDLY_MIME_TYPE_2 = MIME"application/vnd.ijulia.friendly-text-2"
    FRIENDLY_MIME_1 = FRIENDLY_MIME_TYPE_1()
    FRIENDLY_MIME_2 = FRIENDLY_MIME_TYPE_2()
    FRIENDLY_MIME_TYPE_UNION = Union{FRIENDLY_MIME_TYPE_1, FRIENDLY_MIME_TYPE_2}
    Base.Multimedia.istextmime(::FRIENDLY_MIME_TYPE_UNION) = true
    Base.show(io, ::FRIENDLY_MIME_TYPE_UNION, x::FriendlyData) = write(io, "Hello, $(x.name)!")
    IJulia.register_mime([FRIENDLY_MIME_1, FRIENDLY_MIME_2])

    # We stringify then re-parse the dict so that JSONText's are parsed as
    # actual JSON objects and we can index into them.
    data = JSON.parse(JSON.json(IJulia.display_dict(friend)))
    @test data[string(FRIENDLY_MIME)] == "Hello, world!"
    @test data[string(BINARY_MIME)] == base64encode("Hello, world!")
    @test data[string(JSON_MIME)]["name"] == "world"
    @test data[string(FRIENDLY_MIME_1)] == "Hello, world!"
    @test !haskey(data, string(FRIENDLY_MIME_2))

end

struct AngryData
    thing::AbstractString
end

@testset "Render 1st available MIME in MIME-vector." begin
    ANGRY_MIME_TYPE_1 = MIME"application/vnd.ijulia.angry-1"
    ANGRY_MIME_1 = ANGRY_MIME_TYPE_1()
    ANGRY_MIME_TYPE_2 = MIME"application/vnd.ijulia.angry-2"
    ANGRY_MIME_2 = ANGRY_MIME_TYPE_2()
    ANGRY_MIME_TYPE_3 = MIME"application/vnd.ijulia.angry-3"
    ANGRY_MIME_3 = ANGRY_MIME_TYPE_3()
    ANGRY_MIME_VECTOR = [ANGRY_MIME_1, ANGRY_MIME_2, ANGRY_MIME_3]

    Base.Multimedia.istextmime(::Union{ANGRY_MIME_TYPE_1, ANGRY_MIME_TYPE_2, ANGRY_MIME_TYPE_3}) = true
    Base.show(io, ::Union{ANGRY_MIME_TYPE_2, ANGRY_MIME_TYPE_3}, x::AngryData) = write(io, "I hate $(x.thing)!")
    IJulia.register_mime(ANGRY_MIME_VECTOR)

    broccoli = AngryData("broccoli")
    @test IJulia._showable(ANGRY_MIME_VECTOR, broccoli)
    @test !IJulia._showable(ANGRY_MIME_1, broccoli)

    data = IJulia.display_dict(broccoli)
    @test data[string(ANGRY_MIME_2)] == "I hate broccoli!"
    @test !haskey(data, ANGRY_MIME_1)
    @test !haskey(data, ANGRY_MIME_3)
end


text:
using LowRankArithmetic, LinearAlgebra
using Test

include("basic_arithmetic.jl")
include("orthonormalization.jl")
include("rounding.jl")


text:
VERSION >= v"0.6.0" && __precompile__()

module ExcelTables

################################################################################
##
## Dependencies
##
################################################################################

using DataFrames, PyCall, Distributions, GLM, StatsBase, DataStructures,
    HypothesisTests, NamedArrays, FreqTables, Stella, Labels, ROCAnalysis, Dates

##############################################################################
##
## Exported methods and types (in addition to everything reexported above)
##
##############################################################################

export  univariatexls, # output univariate statistics in an excel worksheet
        bivariatexls,  # output bivariate statistics in an excel worksheet
        glmxls,        # output GLM models to an excel worksheet
        mglmxls,       # output multiple GLM regression models to an excel spreadsheet
        dfxls,         # output dataframe in an excel file
        mcfadden,
        nagelkerke,
        hltest

##############################################################################
##
## Load files
##
##############################################################################
include("xlsout.jl")
include("mglmxls.jl")
include("formats.jl")

end # module


text:
# ------------------------------------------------------------------------------------------
# # OnlineStats Integration
#
# - OnlineStats is a Julia package for running statistical algorithms both online (one
# observation at a time) and in parallel.
# - In this notebook we will go over how to run these calculations through JuliaDB.
#
# First we'll load OnlineStats, JuliaDB and Plots (with GR backend)
# ------------------------------------------------------------------------------------------

using OnlineStats, JuliaDB, Plots
gr()

# print table rather than column summary
IndexedTables.set_show_compact!(false);

# ------------------------------------------------------------------------------------------
# Now we'll see:
# 1. An Intro to OnlineStats
# 1. An Example Using OnlineStats with `partitionplot`, `reduce`, and `groupreduce` (all
# available for out-of-core processing!)
# 1. Mosaic Plots
# 1. Linear Regression
# 1. Approximate Solutions to Statistical Learning Models with `StatLearn`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Intro to OnlineStats
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Every stat is a type
# ------------------------------------------------------------------------------------------

m = Mean()

# ------------------------------------------------------------------------------------------
# ### Stats can be `fit!`-ted with more data
# ------------------------------------------------------------------------------------------

fit!(m, randn(100))

# ------------------------------------------------------------------------------------------
# ### Stats can be `merge!`-ed together
# ------------------------------------------------------------------------------------------

m2 = fit!(Mean(), randn(100))

merge!(m, m2)

# ------------------------------------------------------------------------------------------
# ### `fit!`-ting and `merge!`-ing works quite well with JuliaDB
#
# - JuliaDB can send stats to the worker processes and then merge them at the end.
#
# <img src="https://user-
# images.githubusercontent.com/8075494/32748459-519986e8-c88a-11e7-89b3-80dedf7f261b.png"
# width=400>
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Example: Diamonds Dataset
#
# - OnlineStats operations are available through `partitionplot`, `reduce` and
# `groupreduce`.
# - Here we will take a look at a dataset about diamond quality.
#
#
# - First let's examine the first few rows of the data:
# ------------------------------------------------------------------------------------------

;head diamonds.csv

# ------------------------------------------------------------------------------------------
# - The first column is the row number, so we'll only load the other columns.  Here we sort
# the data by column 2, `:carat`.
# ------------------------------------------------------------------------------------------

using JuliaDB
t = loadtable("diamonds.csv", indexcols = [2], datacols = 3:11)

# ------------------------------------------------------------------------------------------
# # `partitionplot`
#
# `partitionplot(table[, x], y; by, stat, nparts, dropmissing)`
#
# - We introduced `partitionplot` in the previous notebook.
# - Here we'll go through some more advanced examples.
# - If only one variable is provided, the x-axis will be the number of observations.
# - Note that `:carat` is our index variable and thus it is sorted:
# ------------------------------------------------------------------------------------------

partitionplot(t, :carat)

# ------------------------------------------------------------------------------------------
# - When we provide two variables, the summary of the `y` variable is plotted in each
# "section" of the `x` variable.
#
# - In this example, we'll use the Interact package to choose which statistic to summarize
# `y` with in the of plot `y = :price` vs. `x = :carat`:
# ------------------------------------------------------------------------------------------

import Interact

Interact.@manipulate for summarizer in [Mean(), Extrema(), Hist(10)]
    partitionplot(t, :carat, :price, stat = summarizer)
end

# ------------------------------------------------------------------------------------------
# - Any selector is a valid input to the `x` and `y` arguments of `partitionplot`:
# ------------------------------------------------------------------------------------------

partitionplot(t, :carat => x -> -x, :price, stat = Hist(10))

# ------------------------------------------------------------------------------------------
# - Any valid Plots keyword argument can also be included, such as changing the type of plot
# series to a barplot.
#
# - Here we are essentially plotting `groupby(mean, t, :cut, select = :carat)`:
# ------------------------------------------------------------------------------------------

partitionplot(t, :cut, :carat; stat = Mean(), seriestype = :bar)

# ------------------------------------------------------------------------------------------
# # `reduce` and `groupreduce`
#
# - In the following examples, we'll see that the reducer for `reduce` and `groupreduce` can
# come from OnlineStats.  The OnlineStats reducers can be:
#   1. A `Series`
#   1. An `OnlineStat`
#   1. A tuple of `OnlineStat`s
#
#
# - Here we get the mean of `:carat` for each level of `:cut`:
# ------------------------------------------------------------------------------------------

groupreduce(Mean(), t, :cut; select = :carat)

# ------------------------------------------------------------------------------------------
# - We could also calculate the same information with `groupby`
# - However `[group]reduce` with an OnlineStat is often more efficient.
# ------------------------------------------------------------------------------------------

@time groupreduce(Mean(), t, :cut; select = :carat)
@time groupby(mean, t, :cut; select = :carat)

# ------------------------------------------------------------------------------------------
# ## Size of Selections
#
# - One thing to note is that the selection in `reduce`/`groupreduce` must match the "input
# size" of the `OnlineStat` object passed as an argument.
# - Different stats can be applied to different columns with the `Group` type.
#     - Multiple stats of the same type can be created via integer multiplication (see
# example below)
# - For example, if you want to calculate the means for 5 different columns, you would use
# as reducer:
# ------------------------------------------------------------------------------------------

5 * Mean()

# ------------------------------------------------------------------------------------------
# - In this example we create 3 histograms of 50 bins each to use the columns `:x`, `:y`,
# and `:z` as input:
# ------------------------------------------------------------------------------------------

hists = reduce(3Hist(50), t; select = (:x, :y, :z))

# ------------------------------------------------------------------------------------------
# - The printed histograms do not provide much information.
# - Fortunately, OnlineStats implements many Plots recipes.
# - Plotting a `Group` will then plot the stats it contains:
# ------------------------------------------------------------------------------------------

plot(hists, layout = 3)

# ------------------------------------------------------------------------------------------
# - You may also want to calculate different statistics for different columns.
# ------------------------------------------------------------------------------------------

reducer = Group(Mean(), CountMap(String))

plot(reduce(reducer, t; select = (:carat, :cut)), layout = 2)

# ------------------------------------------------------------------------------------------
# ## Visualizing Continuous Distributions
#
# - The `Hist` type is particularly useful for visualizing continuous distributions, as it
# adaptively finds the "best" bin locations.
#
# - In the following cell, try increasing the number of bins and see how the histogram
# changes:
# ------------------------------------------------------------------------------------------

# try increasing the number of bins
plot(reduce(Hist(50), t; select = :carat))

# ------------------------------------------------------------------------------------------
# ## Visualizing Categorical Distributions
#
# - The `CountMap` type tracks the number of occurrences for each unique value in a column.
# - `CountMap` accepts the column type as its argument.
# - The `:cut` column is of type `String`, so here we use `CountMap(String)`:
# ------------------------------------------------------------------------------------------

plot(reduce(CountMap(String), t; select = :cut))

# ------------------------------------------------------------------------------------------
# ## Mosaic Plot
#
# - Mosaic plots are extremely useful in visualizing the association between two categorical
# variables, as it shows class probabilities of the `x` variables and conditional class
# probabilities of the `y` variable.
# - `Mosaic` accepts the two column types as its input.
#
#
# - The `:cut` and `:color` variables are both of type `String`, so here we use
# `Mosiac(String, String)`.
# - We can see most common `:cut` is `"Ideal"` and the least common is `"Fair"` (which we
# already know from above).
# ------------------------------------------------------------------------------------------

plot(reduce(Mosaic(String, String), t; select = (:cut, :color)))

# ------------------------------------------------------------------------------------------
# # Statistical Models
#
# - We can build linear (and ridge) regression models using the `LinRegBuilder` and `LinReg`
# types.
#
# ## `LinRegBuilder(p)`
#
# - `LinRegBuilder` builds a data structure that allows you to regress any of the observed
# columns on any subset of other columns.
# ------------------------------------------------------------------------------------------

x = (:carat,:depth,:table,:price,:x,:y,:z)

o = reduce(LinRegBuilder(7), t; select = x)

# ------------------------------------------------------------------------------------------
# - If we wish to fit a model with formula:
# ```
# price ~ carat + depth + table + x + y + z + 1
# ```
#   We can specify the `y` variable to be the 4-th column, which was `:price`.  By default,
# an intercept (bias) term is included as the last coefficient.
# ------------------------------------------------------------------------------------------

coef(o; y=4, verbose=true)

# ------------------------------------------------------------------------------------------
# - We can then create a different regression (without revisiting data) based on the formula
#
#   ```
#   carat ~ x + y + z
#   ```
#
#   by specifying the `y` and `x` variables and removing the `bias`:
# ------------------------------------------------------------------------------------------

coef(o; y=1, x = 5:7, bias = false)

# ------------------------------------------------------------------------------------------
# ## `LinReg()`
#
# - The `LinReg` type allows you to fit linear regression with an optional L2 (ridge)
# penalty.
# - One thing to note is that besides `LinRegBuilder`, models in OnlineStats expect data as
# a tuple: `(x, y)`.
#   - This requires selections to take the form: `((xvars...), yvar)`.
#
#
# - In this example, we'll create the model:
# ```
# carat ~ x + y + z
# ```
# - Note that we found the same result as we did with `LinRegBuilder` above.
# ------------------------------------------------------------------------------------------

o = reduce(LinReg(), t; select = ((:x, :y, :z), :carat))

coef(o, 0.1)

# ------------------------------------------------------------------------------------------
# ## Other Models
#
# - **OnlineStats** has a variety of methods that are more advanced than linear regression.
# - They're beyond the scope of this introduction, but interested readers can investigate
# the following:
#
# ### `StatLearn` (Stochastic approximation for linear statistical learning models)
#
# ### `NBClassifier` (Naive Bayes Classifier)
#
# ### `FastTree`/ `FastForest` (experimental online decision trees/random forests)
#
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
##
# LowRankFun represents f(x,y) by r column and row slices stored as Funs.
##

export LowRankFun

"""
`LowRankFun` gives an approximation to a bivariate function in low rank form.
"""
type LowRankFun{S<:Space,M<:Space,SS<:AbstractProductSpace,T<:Number} <: BivariateFun{T}
    A::Vector{Fun{S,T}}
    B::Vector{Fun{M,T}}
    space::SS

    function LowRankFun(A::Vector{Fun{S,T}},B::Vector{Fun{M,T}},space::SS)
        @assert length(A) == length(B)
        @assert length(A) > 0
        new(A,B,space)
    end
end

LowRankFun{S,M,SS,T}(A::Vector{Fun{S,T}},B::Vector{Fun{M,T}},space::SS)=LowRankFun{S,M,SS,T}(A,B,space)
LowRankFun{S,M,T}(A::Vector{Fun{S,T}},B::Vector{Fun{M,T}})=LowRankFun(A,B,space(first(A))⊗space(first(B)))
LowRankFun{S,M,T,V}(A::Vector{Fun{S,T}},B::Vector{Fun{M,V}})=LowRankFun(convert(Vector{Fun{S,promote_type(T,V)}},A),convert(Vector{Fun{M,promote_type(T,V)}},B),space(first(A))⊗space(first(B)))

Base.rank(f::LowRankFun)=length(f.A)
Base.size(f::LowRankFun,k::Integer)=k==1?mapreduce(length,max,f.A):mapreduce(length,max,f.B)
Base.size(f::LowRankFun)=size(f,1),size(f,2)

## Construction via a Matrix of coefficients

function LowRankFun{S<:Space,M<:Space,T<:Number}(X::Array{T},dx::S,dy::M)
    U,Σ,V=svd(X)
    m=max(1,count(s->s>10eps(T),Σ))

    A=Fun{S,T}[Fun(U[:,k].*sqrt(Σ[k]),dx) for k=1:m]
    B=Fun{M,T}[Fun(conj(V[:,k]).*sqrt(Σ[k]),dy) for k=1:m]

    LowRankFun(A,B)
end

## Construction in a TensorSpace via a Vector of Funs

function LowRankFun{S,T,SV}(X::Vector{Fun{S,T}},d::TensorSpace{SV,T,2})
    @assert d[1] == space(X[1])
    LowRankFun(X,d[2])
end

function LowRankFun{S,T}(X::Vector{Fun{S,T}},dy::Space)
    m=mapreduce(ncoefficients,max,X)
    M=zeros(T,m,length(X))
    for k=1:length(X)
        M[1:ncoefficients(X[k]),k]=X[k].coefficients
    end

    LowRankFun(M,space(X[1]),dy)
end

## Adaptive constructor selector

function LowRankFun(f::Function,dx::Space,dy::Space;method::Symbol=:standard,tolerance::Union{Symbol,Tuple{Symbol,Number}}=:relative,retmax::Bool=false,gridx::Integer=64,gridy::Integer=64,maxrank::Integer=100)
    if method == :standard
        F,maxabsf=standardLowRankFun(f,dx,dy;tolerance=tolerance,gridx=gridx,gridy=gridy,maxrank=maxrank)
    elseif method == :Cholesky
        @assert domain(dx) == domain(dy)
        if dx == dy
            F,maxabsf=CholeskyLowRankFun(f,dx;tolerance=tolerance,grid=max(gridx,gridy),maxrank=maxrank)
        else
            G,maxabsf=CholeskyLowRankFun(f,dx;tolerance=tolerance,grid=max(gridx,gridy),maxrank=maxrank)
            F=LowRankFun(copy(G.A),map(b->Fun(b,dy),G.B))
        end
    end
    retmax ? (F,maxabsf) : F
end

## Standard adaptive construction

function standardLowRankFun(f::Function,dx::Space,dy::Space;tolerance::Union{Symbol,Tuple{Symbol,Number}}=:relative,gridx::Integer=64,gridy::Integer=64,maxrank::Integer=100)
    xy = checkpoints(dx⊗dy)
    T = promote_type(eltype(f(first(xy)...)),eltype(dx),eltype(domain(dx)),eltype(dy),eltype(domain(dy)))

    # We start by sampling on the given grid, find the approximate maximum and create the first rank-one approximation.
    ptsx,ptsy=points(dx,gridx),points(dy,gridy)
    X = zeros(T,gridx,gridy)
    maxabsf,r=findapproxmax!(f,X,ptsx,ptsy,gridx,gridy)
    if maxabsf < eps(zero(T))/eps(T) return LowRankFun([Fun([zero(T)],dx)],[Fun([zero(T)],dy)]),maxabsf end
    a,b=Fun(x->f(x,r[2]),dx),Fun(y->f(r[1],y),dy)

    # If necessary, we resize the grid to be at least as large as the
    # lengths of the first row and column Funs and we recompute the values of X.
    if gridx < ncoefficients(a) || gridy < ncoefficients(b)
        gridx,gridy = max(gridx,ncoefficients(a)),max(gridy,ncoefficients(b))
        ptsx,ptsy=points(dx,gridx),points(dy,gridy)
        X = zeros(T,gridx,gridy)
        maxabsf,r=findapproxmax!(f,X,ptsx,ptsy,gridx,gridy)
        a,b=Fun(x->f(x,r[2]),dx),Fun(y->f(r[1],y),dy)
    end

    A,B=typeof(a)[],typeof(b)[]
    if tolerance == :relative
        tol = 100maxabsf*eps(T)
    elseif tolerance[1] == :absolute
        tol = 100*tolerance[2]*eps(T)
    end
    tol10 = tol/10
    Avals,Bvals = zeros(T,gridx),zeros(T,gridy)
    p₁,p₂ = plan_transform(dx,Avals),plan_transform(dy,Bvals)

    # Eat, drink, subtract rank-one, repeat.
    for k=1:maxrank
        if norm(a.coefficients,Inf) < tol || norm(b.coefficients,Inf) < tol return LowRankFun(A,B),maxabsf end
        A,B =[A;a/sqrt(abs(a(r[1])))],[B;b/(sqrt(abs(b(r[2])))*sign(b(r[2])))]
        r=findapproxmax!(A[k],B[k],X,ptsx,ptsy,gridx,gridy)
        Ar,Br=evaluate(A,r[1]),evaluate(B,r[2])
        for i=1:gridx
            @inbounds Avals[i] = f(ptsx[i],r[2])
        end
        for j=1:gridy
            @inbounds Bvals[j] = f(r[1],ptsy[j])
        end
        a,b = Fun(transform(dx,Avals,p₁),dx) - dotu(Br,A),Fun(transform(dy,Bvals,p₂),dy) - dotu(Ar,B)
        chop!(a,tol10),chop!(b,tol10)
    end
    warn("Maximum rank of " * string(maxrank) * " reached")
    return LowRankFun(A,B),maxabsf
end

## Adaptive Cholesky decomposition, when f is Hermitian positive (negative) definite

function CholeskyLowRankFun(f::Function,dx::Space;tolerance::Union{Symbol,Tuple{Symbol,Number}}=:relative,grid::Integer=64,maxrank::Integer=100)
    xy = checkpoints(dx⊗dx)
    T = promote_type(eltype(f(first(xy)...)),eltype(dx),eltype(domain(dx)))

    # We start by sampling on the given grid, find the approximate maximum and create the first rank-one approximation.
    pts=points(dx,grid)
    X = zeros(T,grid)
    maxabsf,r=findcholeskyapproxmax!(f,X,pts,grid)
    if maxabsf < eps(zero(T))/eps(T) return LowRankFun([Fun([zero(T)],dx)],[Fun([zero(T)],dx)]),maxabsf end
    a=Fun(x->f(x,r),dx)

    # If necessary, we resize the grid to be at least as large as the
    # ncoefficients of the first row/column Fun and we recompute the values of X.
    if grid < ncoefficients(a)
        grid = max(grid,ncoefficients(a))
        pts=points(dx,grid)
        X = zeros(T,grid)
        maxabsf,r=findcholeskyapproxmax!(f,X,pts,grid)
        a=Fun(x->f(x,r),dx)
    end

    A,B=typeof(a)[],typeof(a)[]
    if tolerance == :relative
        tol = 100maxabsf*eps(T)
    elseif tolerance[1] == :absolute
        tol = 100*tolerance[2]*eps(T)
    end
    tol10 = tol/10
    Avals = zeros(T,grid)
    p₁ = plan_transform(dx,Avals)

    # Eat, drink, subtract rank-one, repeat.
    for k=1:maxrank
        if norm(a.coefficients,Inf) < tol return LowRankFun(A,B),maxabsf end
        A,B = [A;a/sqrt(abs(a(r)))],[B;a/(sqrt(abs(a(r)))*sign(a(r)))]
        r=findcholeskyapproxmax!(A[k],B[k],X,pts,grid)
        Br=evaluate(B,r)
        for i=1:grid
            @inbounds Avals[i] = f(pts[i],r)
        end
        a = Fun(transform(dx,Avals,p₁),dx) - dotu(Br,A)
        chop!(a,tol10)
    end
    warn("Maximum rank of " * string(maxrank) * " reached")
    return LowRankFun(A,B),maxabsf
end

## Construction via TensorSpaces and ProductDomains

LowRankFun{SV,T}(f::Function,S::TensorSpace{SV,T,2};kwds...)=LowRankFun(f,S[1],S[2];kwds...)
LowRankFun(f::Function,dx::Domain,dy::Domain;kwds...)=LowRankFun(f,Space(dx),Space(dy);kwds...)
LowRankFun{D,T}(f::Function,d::ProductDomain{D,T,2};kwds...)=LowRankFun(f,d[1],d[2];kwds...)

LowRankFun(f::Function,d1::Vector,d2::Vector;kwds...)=LowRankFun(f,convert(Domain,d1),convert(Domain,d2);kwds...)
LowRankFun(f::Function;kwds...)=LowRankFun(f,Interval(),Interval();kwds...)

## Construction from values

LowRankFun{T<:Number}(A::Array{T})=LowRankFun(A,Interval{T}(),Interval{T}())
LowRankFun(c::Number,etc...)=LowRankFun((x,y)->c,etc...)

## Construction from other LowRankFuns

LowRankFun(f::LowRankFun,d1::IntervalDomain,d2::IntervalDomain)=LowRankFun(map(g->Fun(g.coefficients,d1),f.A),map(g->Fun(g.coefficients,d2),f.B))
LowRankFun(f::LowRankFun)=LowRankFun(f,Interval(),Interval())

## Utilities

function findapproxmax!(f::Function,X::Matrix,ptsx::Vector,ptsy::Vector,gridx,gridy)
    for j=1:gridy
        ptsyj = ptsy[j]
        @simd for k=1:gridx
            @inbounds X[k,j]+=f(ptsx[k],ptsyj)
        end
    end
    maxabsf,impt = findmaxabs(X)
    imptple = ind2sub((gridx,gridy),impt)
    maxabsf,[ptsx[imptple[1]],ptsy[imptple[2]]]
end

function findapproxmax!(A::Fun,B::Fun,X::Matrix,ptsx::Vector,ptsy::Vector,gridx,gridy)
    Ax,By = A(ptsx),B(ptsy)
    subtractrankone!(Ax,By,X,gridx,gridy)
    maxabsf,impt = findmaxabs(X)
    imptple = ind2sub((gridx,gridy),impt)
    [ptsx[imptple[1]],ptsy[imptple[2]]]
end

function findcholeskyapproxmax!(f::Function,X::Vector,pts::Vector,grid)
    @simd for k=1:grid
        @inbounds X[k]+=f(pts[k],pts[k])
    end
    maxabsf,impt = findmaxabs(X)
    maxabsf,pts[impt]
end

function findcholeskyapproxmax!(A::Fun,B::Fun,X::Vector,pts::Vector,grid)
    Ax,By = A(pts),B(pts)
    subtractrankone!(Ax,By,X,grid)
    maxabsf,impt = findmaxabs(X)
    pts[impt]
end

function subtractrankone!(A::AbstractVector,B::AbstractVector,X::AbstractMatrix,gridx::Int,gridy::Int)
    for j=1:gridy
        @inbounds Bj = B[j]
        @simd for k=1:gridx
            @inbounds X[k,j] -= A[k]*Bj
        end
    end
end

function subtractrankone!(A::AbstractVector,B::AbstractVector,X::AbstractVector,grid::Int)
    @simd for k=1:grid
        @inbounds X[k] -= A[k]*B[k]
    end
end

## TODO: in Julia base?
function findmaxabs(a)
    if isempty(a)
        throw(ArgumentError("collection must be non-empty"))
    end
    m = abs(a[1])
    mi = 1
    for i in eachindex(a)
        ai = abs(a[i])
        if ai > m || m!=m
            m = ai
            mi = i
        end
    end
    return (m, mi)
end

@compat (f::LowRankFun)(x,y)=evaluate(f,x,y)

domain(f::LowRankFun,k::Integer)=k==1? domain(first(f.A)) : domain(first(f.B))
space(f::LowRankFun,k::Integer)=k==1? space(first(f.A)) : space(first(f.B))
space(f::LowRankFun)=f.space

Base.transpose{S,M,SS,T}(f::LowRankFun{S,M,SS,T})=LowRankFun(f.B,f.A,transpose(space(f)))

function values(f::LowRankFun)
    xm=mapreduce(ncoefficients,max,f.A)
    ym=mapreduce(ncoefficients,max,f.B)
    ret=zeros(xm,ym)
    for k=1:length(f.A)
        ret+=values(pad(f.A[k],xm))*values(pad(f.B[k],ym)).'
    end
    ret
end

#TODO: this is inconsistent with 1D where it does canonical
function coefficients(f::LowRankFun)
    xm=mapreduce(ncoefficients,max,f.A)
    ym=mapreduce(ncoefficients,max,f.B)
    ret=zeros(xm,ym)
    for k=1:length(f.A)
        ret+=pad(f.A[k].coefficients,xm)*pad(f.B[k].coefficients,ym).'
    end
    ret
end

function coefficients(f::LowRankFun,n::Space,m::Space)
    xm=mapreduce(ncoefficients,max,f.A)
    ym=mapreduce(ncoefficients,max,f.B)
    ret=zeros(xm,ym)
    for k=1:length(f.A)
        ret+=pad(coefficients(f.A[k],n),xm)*pad(coefficients(f.B[k],m),ym).'
    end
    ret
end

function vecpoints(f::LowRankFun,k::Integer)
    if k==1
        xm=mapreduce(ncoefficients,max,f.A)
        points(space(first(f.A)),xm)
    else
        ym=mapreduce(ncoefficients,max,f.B)
        points(space(first(f.B)),ym)
    end
end

evaluate{T<:Fun,M<:Fun}(A::Vector{T},B::Vector{M},x,y)=dotu(evaluate(A,x),evaluate(B,y))
evaluate{T<:Fun,M<:Fun}(A::Vector{T},B::Vector{M},x::AbstractVector,y::AbstractVector)=evaluate(A,x).'*evaluate(B,y)

evaluate(f::LowRankFun,x,y)=evaluate(f.A,f.B,x,y)
evaluate(f::LowRankFun,::Colon,::Colon)=f
evaluate(f::LowRankFun,x::Number,::Colon)=dotu(f.B,evaluate(f.A,x))
evaluate{T<:Number}(f::LowRankFun,x::Vector{T},::Colon)=f.B.'*evaluate(f.A,x)
function evaluate(f::LowRankFun,::Colon,y::Number)
    m = maximum(map(ncoefficients,f.A))
    r=rank(f)
    ret = zeros(m)

    for k=1:r
        for j=1:ncoefficients(f.A[k])
            @inbounds ret[j] += f.A[k].coefficients[j]*f.B[k](y)
        end
    end

    Fun(ret,first(f.A).space)
end

## Truncate
#TODO: should reduce rank if needed
Base.chop(f::LowRankFun,tol)=LowRankFun(map(g->chop(g,tol),f.A),map(g->chop(g,tol),f.B),f.space)
function pad(f::LowRankFun,m::Integer,n::Integer)
    A,B = deepcopy(f.A),deepcopy(f.B)

    for k=1:rank(f)
        pad!(A[k],m);pad!(B[k],n)
    end

    LowRankFun(A,B,f.space)
end

## Algebra

for op = (:*,:.*,:./,:/)
    @eval ($op){T<:Fun}(A::Array{T,1},c::Number)=map(f->($op)(f,c),A)
    @eval ($op)(f::LowRankFun,c::Number) = LowRankFun(($op)(f.A,c),f.B)
    @eval ($op)(c::Number,f::LowRankFun) = LowRankFun(($op)(c,f.A),f.B)
end

# Let K be a LowRankFun and f be a Fun.
# op(f,K) acts as operating in the x variable, and
# op(K,f) acts as operating in the y variable.

for op = (:*,:.*,:./,:/)
    @eval ($op){S,T,U,V}(f::Fun{S,T},A::Vector{Fun{U,V}})=map(a->($op)(f,a),A)
    @eval ($op)(f::Fun,K::LowRankFun) = LowRankFun(($op)(f,K.A),K.B)
    @eval ($op){S,T,U,V}(B::Vector{Fun{U,V}},f::Fun{S,T})=map(b->($op)(b,f),B)
    @eval ($op)(K::LowRankFun,f::Fun) = LowRankFun(K.A,($op)(K.B,f))
end

+(f::LowRankFun,g::LowRankFun)=LowRankFun([f.A;g.A],[f.B;g.B])
-(f::LowRankFun)=LowRankFun(-f.A,f.B)
-(f::LowRankFun,g::LowRankFun)=f+(-g)

## QR factorization of a LowRankFun

function Base.qr(f::LowRankFun)
    sp,r = space(f),rank(f)
    Q,R = qr(coefficients(f.A))
    BR = coefficients(f.B)*R.'
    LowRankFun(map(i->Fun(Q[:,i],sp[1]),1:r),map(i->Fun(BR[:,i],sp[2]),1:r),sp)
end

## Special functions

Base.real(u::LowRankFun)=LowRankFun([map(real,u.A),map(imag,u.A)],[map(real,u.B),-map(imag,u.B)])
Base.imag(u::LowRankFun)=LowRankFun([map(real,u.A),map(imag,u.A)],[map(imag,u.B),map(real,u.B)])

## Calculus

Base.sum(g::LowRankFun)=dotu(map(sum,g.A),map(sum,g.B))
Base.sum(g::LowRankFun,n::Integer)=(n==1)?dotu(g.B,map(sum,g.A)):dotu(g.A,map(sum,g.B))
Base.cumsum(g::LowRankFun,n::Integer)=(n==1)?LowRankFun(map(cumsum,g.A),copy(g.B)):LowRankFun(copy(g.A),map(cumsum,g.B))
differentiate(g::LowRankFun,n::Integer)=(n==1)?LowRankFun(map(differentiate,g.A),copy(g.B)):LowRankFun(copy(g.A),map(differentiate,g.B))
integrate(g::LowRankFun,n::Integer)=(n==1)?LowRankFun(map(integrate,g.A),copy(g.B)):LowRankFun(copy(g.A),map(integrate,g.B))


text:
# Variation of Information

"""
    varinfo(a, b) -> Float64

Compute the *variation of information* between the two clusterings of the same
data points.

`a` and `b` can be either [`ClusteringResult`](@ref) instances or
assignments vectors (`AbstractVector{<:Integer}`).

# References
> Meila, Marina (2003). *Comparing Clusterings by the Variation of
> Information.* Learning Theory and Kernel Machines: 173–187.
"""
function varinfo(a, b)
    C = counts(a, b)
    isempty(C) && return 0.0
    countsA = a isa ClusteringResult ? counts(a) : sum(C, dims=2)
    countsB = b isa ClusteringResult ? counts(b) : sum(C, dims=1)
    I = 0.0
    @inbounds for (i, ci) in enumerate(countsA), (j, cj) in enumerate(countsB)
        cij = C[i, j]
        if cij > 0.0
            I += cij * log(cij*cij / (ci*cj))
        end
    end
    return -I/sum(countsA)
end


text:
using SensorFeatureTracking
using Base: Test

@testset "SensorFeatureTracking" begin
    # Common
    include("sensorgeometry.jl")
    include("blockflow.jl")

end

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
using LinearAlgebra, Rotations

function nullspace_method(jac_joints, delta; regularization_strength=0.0)
    hess_approx = transpose(jac_joints)*jac_joints
    joint_delta = transpose(jac_joints)*delta

    if regularization_strength > 0.0
        hess_approx += I(size(hess_approx)[1])*regularization_strength
        return hess_approx\joint_delta
    else
        return hess_approx\joint_delta
    end
end

function inverse_kinematics!(sim::MJSim, site_name::String,
    target_pos::Union{Vector{Float64},Nothing},
    target_quat::Union{Vector{Float64},Nothing}=nothing,
    joint_names=nothing,
    tol::Float64=1E-14,
    rot_weight::Float64=1.0,
    regularization_threshold::Float64=0.1,
    regularization_strength::Float64=3E-2,
    max_update_norm::Float64=2.0,
    progress_thresh::Float64=20.0,
    max_steps::Int64=100)

    if isnothing(target_quat)
        jac = Array{Float64}(undef, 3, sim.m.nv)#zeros((3, sim.m.nv))
        err = Array{Float64}(undef, 3)#zeros(3)
        jac_pos, jac_rot = jac, nothing
        err_pos, err_rot = err, nothing

    else
        jac = Array{Float64}(undef, 6, sim.m.nv)#zeros((6, sim.m.nv))
        err = Array{Float64}(undef, 6)#zeros(6)
        jac_pos, jac_rot = jac[1:4], jac[4:end]
        err_pos, err_rot = err[1:4], err[4:end]
    end

    update_nv = zeros(sim.m.nv)

    if !isnothing(target_quat)
        site_xquat = Array{Float64}(undef, 4)#zeros(4)
        neg_site_xquat = Array{Float64}(undef, 4)#zeros(4)
        err_rot_quat = Array{Float64}(undef, 4)
    end

    mj_fwdPosition(sim.m, sim.d)

    site_id = jl_name2id(sim.m, MJCore.mjOBJ_SITE, site_name)

    site_xpos = sim.d.site_xpos[:,site_id]
    site_xmat = sim.d.site_xmat[:,site_id]
    println("site pos: ",site_xpos)

    if isnothing(joint_names)
        dof_indices = sim.m.dof_jntid
    else
        dof_indices = []
        for joint_name in joint_names
            jid = jl_name2id(sim.m, MJCore.mjOBJ_JOINT, joint_name)
            push!(dof_indices, jid)
        end
        println("dof inds: ",dof_indices)
    end

    steps = 0
    success = false

    for steps = 1:max_steps

        err_norm = 0.0

        if !isnothing(target_pos)
            err_pos = target_pos - site_xpos
            err_norm += norm(err_pos)
        end

        if !isnothing(target_quat)
            # q = UnitQuaternion(SMatrix{3,3}(transpose(reshape(site_xmat, (3,3)))))
            MJCore.mju_mat2Quat(site_xquat, site_xmat)
            MJCore.mju_negQuat(neg_site_xquat, site_xquat)
            MJCore.mju_mulQuat(err_rot_quat, target_quat, neg_site_xquat)
            println("error rot: ",err_rot)
            MJCore.mju_quat2Vel(err_rot, err_rot_quat, 1.)
            err_norm += norm(err_rot)*rot_weight
        end

        if err_norm < tol
            success = true
            break
        else
            mj_jacSite(sim.m, sim.d, jac_pos, jac_rot, site_id)
            jac_joints = Array{Float64}(undef, 6, length(dof_indices))
            for i = 1:length(dof_indices)
                jac_joints[:,i] = jac[:,dof_indices[i]]
            end

            reg_strength = err_norm > regularization_threshold ?
                            regularization_strength : 0.0

            update_joints = nullspace_method(jac_joints, err,
                                regularization_strength=reg_strength)
            update_norm = norm(update_joints)
            println(update_joints)
            progress_criterion = err_norm/update_norm
            if progress_criterion > progress_thresh
                break
            end
            if update_norm > max_update_norm
                update_joints *= max_update_norm/update_norm
            end

            for i in dof_indices
                update_nv[i] = update_joints[i]
            end

            mj_integratePos(sim.m, sim.d.qpos, update_nv, 1.)
            mj_fwdPosition(sim.m, sim.d)
        end
    end

    if !success && steps == max_steps
        println("Failed to converge")
    end

    qpos = sim.d.qpos
    return qpos
end


text:
using CSV
using DataFrames
using Plots

data_r = CSV.read(homedir() * "/memory_usage_rust.csv", DataFrame)
data_f = CSV.read(homedir() * "/memory_usage_fortran.csv", DataFrame)

timestamp_r = data_r[:, 1]
allocated_r = data_r[:, 2] ./ (1024^3)

timestamp_f = data_f[:, 1]
allocated_f = data_f[:, 2] ./ (1024^3)

common = min(size(timestamp_r)[1], size(timestamp_f)[1])

plot(timestamp_f[1:common], [allocated_r[1:common], allocated_f[1:common]], label = ["Rust fits_web_ql v4" "C / FORTRAN FITSWEBQLSE v5"], xlabel = "elapsed time [s]", ylabel = "allocated memory [GB]", legend = :bottomright, title = "memory consumption comparison")
savefig(homedir() * "/mem_comparison.pdf")


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -6.103907844236108 <= q <= 8.378675476959124)

                     add_NL_constraint(m, :((-0.8297607126794757 + -0.5064858067335005*(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))*(1 + erf((-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))/sqrt(2)))/2 + -0.8619658571743893*(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))*(1 + erf((-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))/sqrt(2)))/2 + 0.8415580445654567*(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))*(1 + erf((-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.8297607126794757 + -0.5064858067335005*(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))*(1 + erf((-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))/sqrt(2)))/2 + -0.8619658571743893*(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))*(1 + erf((-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))/sqrt(2)))/2 + 0.8415580445654567*(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))*(1 + erf((-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.0809505831099564 + 0.6185835236357913*(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))*(1 + erf((-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))/sqrt(2)))/2 + -0.5288114683363174*(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))*(1 + erf((-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))/sqrt(2)))/2 + -0.42070009030061817*(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))*(1 + erf((-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.0809505831099564 + 0.6185835236357913*(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))*(1 + erf((-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))/sqrt(2)))/2 + -0.5288114683363174*(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))*(1 + erf((-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))/sqrt(2)))/2 + -0.42070009030061817*(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))*(1 + erf((-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.16011867292674786 + -0.14898607305663658*(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))*(1 + erf((-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))/sqrt(2)))/2 + -0.8220978033005912*(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))*(1 + erf((-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))/sqrt(2)))/2 + -0.8518875556726564*(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))*(1 + erf((-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))/sqrt(2)))/2)*(1 + erf((0.16011867292674786 + -0.14898607305663658*(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))*(1 + erf((-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2]))/sqrt(2)))/2 + -0.8220978033005912*(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))*(1 + erf((-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2]))/sqrt(2)))/2 + -0.8518875556726564*(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))*(1 + erf((-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
###########
# Image preparation
#
# TODO
# - scale
# - co-registration

"""
    crop(img, topleft, bottomright)

Crops an image.
"""
function crop(img, topleft, bottomright)
    i1, i2 = indices(img)
    img[max(topleft[1],i1[1]):min(bottomright[1],i1[end]),
        max(topleft[2],i2[1]):min(bottomright[2],i2[end])]
end
crop(img, left, right, down, up) = img[1+left:end-right,1+down:end-up]
"""
    thin(img,step)

Thin an image
"""
function thin(img,step, alg=[:imresize, :restrict, :subsampling][3])
    if alg==:imresize
        ## using filtering: a bit slow
        sz = (length(1:step:size(img,1)), length(1:step:size(img,2)))
        σ = maximum(map((o,n)->0.25*o/n, size(img), sz))
        kern = KernelFactors.gaussian((σ,σ))   # from ImageFiltering
        return imresize(imfilter(img, kern, NA()), sz)
    elseif alg==:restrict
        ### using restrict, slightly faster
        l2 = log2(step)
        @assert isinteger(l2) "Can only thin by power of 2."
        for i=1:l2
            s = size(img)
            img = ImageFiltering.padarray(img, ImageFiltering.Pad(:reflect,1,1)) # does not work
            img = ImageTransformations.restrict(img)[1:fld1(s[1],2),1:fld1(s[2],2)]
        end
        return img
    elseif alg==:subsampling
        ### just sub-sampling: causes aliasing
        return img[1:step:end,1:step:end]
    end
    error()
end

"""
    rotate_n_crop(img, exi::ExpImgs; verbose=false)
    rotate_n_crop(img, p1::Tuple{Int,Int}, p2::Tuple{Int,Int},
                  halfheight::Integer; verbose=false)

Rotates the image such that the line through p1 and p2 is horizontal.
It then crops it such that the width spans p1 and p2 and the
height is 2x halfheight.

Besides the processed image it also returns the distance of the (p1,p2) line to
the center of the original image.  This is needed to calculate real distances:
`top+center_dist` and `bottom-center_dist` gives the number of pixels from a line
through the center-point.

Notes: needs to be in sync with ExpImgs defaults definition.
"""
rotate_n_crop(img, exi::ExpImgs; verbose=false) =
    rotate_n_crop(img, exi.p1, exi.p2, exi.halfheight_crop, verbose=verbose)
function rotate_n_crop(img, p1::Tuple{Int,Int}, p2::Tuple{Int,Int}, halfheight::Integer; verbose=false)
    mid_point = map(x->fld1(x,2), size(img))
    if verbose
        guidict = imshow(img);
        annotate!(guidict, AnnotationPoint(p1[2], p1[1], shape='x', size=20, linewidth=2));
        annotate!(guidict, AnnotationPoint(p2[2], p2[1], shape='x', size=20, linewidth=2));
    end
    @assert p1[2]<=p2[2] "point p1 must be on the left of p2"
    angle = atan2(p2[1]-p1[1], p2[2]-p1[2])
    len = round(Int, sqrt( (p2[1]-p1[1])^2 + (p2[2]-p1[2])^2 ))
    tfm = recenter(RotMatrix(-angle), [p1...])
    img2 = warp(img, tfm);
    topleft = (p1[1]-halfheight, p1[2])
    bottomright = (p1[1]+halfheight, p1[2]+len)
    if verbose
        guidict = imshow(img2);
        annotate!(guidict, AnnotationPoint(p1[2], p1[1], shape='x', size=20, linewidth=2));
        annotate!(guidict, AnnotationBox(topleft[2], topleft[1], bottomright[2], bottomright[1], linewidth=2, color=RGB(0,0,1)))
    end
    center_dist = cross([(mid_point.-p1)...,0], [(p2.-p1)...,0])[3]/norm([(p1.-p2)...])
    crop(img2, topleft, bottomright), round(Int, center_dist)
end

"""
    prep_img(path::String, exi::ExpImgs; verbose=false)
    prep_img(img_color, exi; verbose=false)

Load and prepare image by:
- rotate and crop
- colordiff it

Also returns:
- center_dist -- distance to previous center
- capture_time -- time of image capture
"""
function prep_img(path::String, exi::ExpImgs; verbose=false)
    verbose && println(path)
    img_color = load(path);

    dt = ImageMagick.magickinfo(path, "exif:DateTimeOriginal")["exif:DateTimeOriginal"]
    _prep_img(img_color, exi; verbose=verbose)..., DateTime(dt, "yyyy:mm:dd HH:MM:SS") - exi.time_correction
end
function _prep_img(img_color::AbstractArray, exi::ExpImgs; verbose=false)::Tuple{Matrix{Float64}, Int}
    @unpack p1, p2, halfheight_crop, thin_num = exi
    img_color, center_dist = rotate_n_crop(img_color, p1, p2, halfheight_crop, verbose=verbose)
    img_color = thin(img_color, thin_num);
    @assert size(img_color)==size(exi) "size(img_color)=$(size(img_color)) not equal size(exi)=$(size(exi))"
    # calculate the difference in color
    img = colordiffit(img_color, exi, verbose=verbose);
    # find EXIF time and correct it
    return img, center_dist
end

"""
    colordiffit(img_color, exi::ExpImgs; verbose=false)

Use `colordiff` make an image of perceived color difference
"""
function colordiffit(img_color, exi::ExpImgs; verbose=false)
    loc1, loc2 = round.(Int, (size(img_color,1)*exi.color_loc[1], size(img_color,2)*exi.color_loc[2]))
    c0 = img_color[loc1,loc2]
    if verbose
        guidict = imshow(img_color)
        dp = size(img_color,1)÷50
        annotate!(guidict, AnnotationBox(loc2+dp, loc1+dp, loc2-dp, loc1-dp, linewidth=2, color=RGB(0,0,1)))
    end
    colordiff.(img_color, c0)
end


text:
function f(x::Integer, y::Integer)
    return x, y
end

a = zeros(3, 3)
b = zeros(3, 3)

function g()
    for i = 1 : 3
        for j = 1 : 3
            a[i, j], b[i, j] = f(i, j)
        end
    end
end
println(a, b)

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
#####################################
#     Basic interface functions     #
#####################################

"""
    cost,i1,i2 = dtw(seq1, seq2, [dist=SqEuclidean])

Find a set of indices (`i1`,`i2`) that align two time series (`seq1`,`seq2`) by
dynamic time warping. Also returns the distance (after warping) according to
the SemiMetric `dist`, which defaults to squared Euclidean distance (see
Distances.jl). If `seq1` and `seq2` are matrices, each column is considered
an observation.
"""
function dtw(
        seq1::Sequence,
        seq2::Sequence,
        dist::SemiMetric=SqEuclidean()
    )

    D = dtw_cost_matrix(seq1, seq2, dist)
    return trackback(D)
end

"""
    cost,i1,i2 = dtw(seq1,seq2,i2min,i2max,[dist=SqEuclidean])

Do DTW to align `seq1` and `seq2` confined to a window. Vectors `i2min` and
`i2max` specify (inclusive) lower and upper bounds for `seq2` for each index in
`seq1`. Thus, `i2min` and `i2max` are required to be the same length as `seq1`.
"""
function dtw(
        seq1::Sequence,
        seq2::Sequence,
        i2min::AbstractVector,
        i2max::AbstractVector,
        dist::SemiMetric = SqEuclidean()
    )

    D = dtw_cost_matrix(seq1, seq2, i2min, i2max, dist)
    return trackback(D)
end

##############################
#  Cost matrix computations  #
##############################
function dtw_cost_matrix(
        seq1::Sequence{N,T},
        seq2::Sequence{N,T},
        dist::SemiMetric=SqEuclidean()
    ) where {N,T}
    # Build the cost matrix
    m = length(seq2)
    n = length(seq1)
    D = zeros(T, m, n)

    # Initialize first column and first row
    D[1,1] = evaluate(dist, seq1[1], seq2[1])
    for r=2:m
        D[r,1] = D[r-1,1] + evaluate(dist, seq1[1], seq2[r])
    end
    for c=2:n
        D[1,c] = D[1,c-1] + evaluate(dist, seq1[c], seq2[1])
    end

    # Complete the cost matrix
    for c=2:n
        for r=2:m
            best_neighbor_cost = min(D[r-1,c], D[r-1,c-1], D[r,c-1])
            D[r,c] = best_neighbor_cost + evaluate(dist, seq1[c], seq2[r])
        end
    end

    return D
end

function dtw_cost_matrix(
        seq1::Sequence{N,T},
        seq2::Sequence{N,T},
        i2min::AbstractVector{U},
        i2max::AbstractVector{U},
        dist::SemiMetric = SqEuclidean()
    ) where {N,T,U<:Integer}
    m = length(seq2) # of rows in cost matrix
    n = length(seq1) # of columns in cost matrix
    n == length(i2min) || throw(ArgumentError("i2min does not match length of seq1."))
    n == length(i2max) || throw(ArgumentError("i2max does not match length of seq1."))
    1 == i2min[1] || throw(ArgumentError("i2min must start at 1."))
    m == i2max[end] || throw(ArgumentError("i2max must end at length(seq2)."))

    # Build the (n x m) cost matrix into a WindowedMatrix, because it's ragged.
    # That type gives efficient storage with convenient [r,c] indexing and returns
    # Inf when accessed outside the window.
    D = WindowedMatrix(i2min, i2max, Inf)

    # First column first
    D[1,1] = evaluate(dist, seq1[1], seq2[1])
    for r=2:i2max[1]
        D[r,1] = D[r-1,1]  + evaluate(dist, seq1[1], seq2[r])
    end

    # Complete the cost matrix from columns 2 to m.
    for c=2:n
        for r=i2min[c]:i2max[c]
            best_neighbor_cost = min(D[r-1,c], D[r-1,c-1], D[r,c-1])
            D[r,c] = best_neighbor_cost + evaluate(dist, seq1[c], seq2[r])
        end
    end

    return D
end

########################################
#  Find Best Path through Cost Matrix  #
########################################

"""
    cost,cols,rows = trackback(D::Matrix)

Given the cost matrix `D`, computes the optimal track from end to beginning.
Returns `cols` and `rows` which are vectors respectively holding the track.
"""
function trackback(D::AbstractMatrix{T}) where {T<:Number}

    # initialize trackback throught rows/columns
    r,c = size(D)
    rows,cols = Int[r],Int[c]

    # estimate that we'll need N⋅logN elements
    N = max(r,c)
    sz = 2*N
    sizehint!(rows,sz)
    sizehint!(cols,sz)

    # do trackback
    while r > 1 && c > 1
        tb,r,c = indmin3(D[r-1,c-1], D[r-1,c], D[r,c-1], r, c)
        push!(rows,r)
        push!(cols,c)
    end
    # Possibly either r>1 or c>1 at this point (but not both).
    # Add the unfinished part of the track to reach [1,1]
    for r=r-1:-1:1
        push!(rows,r)
        push!(cols,1)
    end
    for c=c-1:-1:1
        push!(rows,1)
        push!(cols,c)
    end
    return D[end,end], reverse(cols), reverse(rows)
end

###############################################
# Wrapper functions for abstract array inputs #
###############################################

function dtw(
        seq1::AbstractArray,
        seq2::AbstractArray,
        dist::SemiMetric=SqEuclidean()
    )
    dtw(Sequence(seq1),Sequence(seq2),dist)
end

function dtw(
        seq1::AbstractArray,
        seq2::AbstractArray,
        i2min::AbstractVector,
        i2max::AbstractVector,
        dist::SemiMetric = SqEuclidean()
    )
    dtw(Sequence(seq1),Sequence(seq2),i2min,i2max,dist)
end


text:
function stepsize_controller!(integrator,alg)
    integrator.q11 = DiffEqBase.value(DiffEqBase.fastpow(integrator.EEst,integrator.opts.beta1))
    integrator.q = DiffEqBase.value(integrator.q11/DiffEqBase.fastpow(integrator.qold,integrator.opts.beta2))
    @fastmath integrator.q = DiffEqBase.value(max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),integrator.q/integrator.opts.gamma)))
end

function step_accept_controller!(integrator,alg)
    integrator.dtnew = DiffEqBase.value(integrator.dt/integrator.q) * oneunit(integrator.dt)
end

function step_reject_controller!(integrator,alg)
    integrator.dtnew = integrator.dt/min(inv(integrator.opts.qmin),integrator.q11/integrator.opts.gamma)
end

function stepsize_controller!(integrator,alg::TauLeaping)

end

function step_accept_controller!(integrator,alg::TauLeaping)
    integrator.q = min(integrator.opts.gamma/integrator.EEst,integrator.opts.qmax)
    integrator.dtnew = integrator.dt*integrator.q
end

function step_reject_controller!(integrator,alg::TauLeaping)
    integrator.dtnew = integrator.opts.gamma * integrator.dt/integrator.EEst
end

function stepsize_controller!(integrator,alg::CaoTauLeaping)

end

function step_accept_controller!(integrator,alg::CaoTauLeaping)
    integrator.dtnew = integrator.EEst # use EEst for the τ
end

function step_reject_controller!(integrator,alg::CaoTauLeaping)
    error("CaoTauLeaping should never reject steps")
end


text:
# ## Rose plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Random
Random.seed!(2018)

n = 24
R = rand(n+1)
θ = 0:2pi/n:2pi
plot(θ, R, proj=:polar, line=:steppre, lims=(0:2))


text:
abstract type AbstractContact{T} <: AbstractObject{T} end

"""
    mutable struct Contact{T} <: AbstractContact{T}

T: Type of precision.
"""
mutable struct Contact{T} <: AbstractContact{T}
    potential::T
    material::NamedTuple
    id::Int
    name::String
    geometry::AbstractGeometry{T}
    geometry_positive::Vector{AbstractGeometry{T}}
    geometry_negative::Vector{AbstractGeometry{T}}
end

function Contact{T}(dict::Union{Dict{String,Any}, Dict{Any, Any}}, inputunit_dict::Dict{String,Unitful.Units})::Contact{T} where {T <: SSDFloat}
    haskey(dict, "channel") ? channel = dict["channel"] : channel = -1
    haskey(dict, "material") ? material = material_properties[materials[dict["material"]]] : material = material_properties[materials["HPGe"]]
    haskey(dict,"name") ? name = dict["name"] : name = ""
    geometry =  Geometry(T, dict["geometry"], inputunit_dict )
    geometry_positive, geometry_negative = get_decomposed_volumes(geometry)
    return Contact{T}( dict["potential"], material, channel, name, geometry, geometry_positive, geometry_negative )
end

function println(io::IO, d::Contact{T}) where {T <: SSDFloat}
    println("\t________"*"Contact $(d.id)"*"________\n")
    println("\t---General Properties---")
    println("\t-Potential: \t\t $(d.potential)")
    println("\t-Contact Material: \t $(d.material.name)")
    println()

end

function show(io::IO, d::Contact{T}) where {T <: SSDFloat} println(d) end
function print(io::IO, d::Contact{T}) where {T <: SSDFloat} println(d) end
function display(io::IO, d::Contact{T} ) where {T <: SSDFloat} println(d) end
function show(io::IO,::MIME"text/plain", d::Contact) where {T <: SSDFloat}
    show(io, d)
end