text:
struct Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis} <: Adjoint_rep_Gaugefields_4D{NC,NumofBasis}
    U::Array{ComplexF64,6}
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NDW::Int64
    NV::Int64
    NC::Int64
    mpi::Bool
    verbose_print::Verbose_print
    NumofBasis::Int64
    generators::Union{Nothing,Generator}

    function Adjoint_rep_Gaugefields_4D_wing(NumofBasis,NC::T,NDW::T,NX::T,NY::T,NZ::T,NT::T;verbose_level = 2) where T<: Integer

        NV = NX*NY*NZ*NT
        U = zeros(ComplexF64,NumofBasis,NumofBasis,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW)
        mpi = false
        verbose_print = Verbose_print(verbose_level )
        generators = Generator(NC)
        #=
        if NC <= 3
            generators = nothing
        else
            generators = Generator(NC)
        end
        =#

        return new{NC,NumofBasis}(U,NX,NY,NZ,NT,NDW,NV,NC,mpi,verbose_print,
                    NumofBasis,generators)
    end

end

function construct_Adjoint_rep_Gaugefields_4D_wing(Uin::Gaugefields_4D{NC};verbose_level=2) where NC
    return Adjoint_rep_Gaugefields_4D_wing(Uin,verbose_level=verbose_level)
end

function Adjoint_rep_Gaugefields_4D_wing(Uin::Gaugefields_4D{NC};verbose_level=2)  where NC
    NX = Uin.NX
    NY = Uin.NY
    NZ = Uin.NZ
    NT = Uin.NT
    NumofBasis = ifelse(NC == 1,1,NC^2-1)
    NDW = Uin.NDW
    Uadj = Adjoint_rep_Gaugefields_4D_wing(NumofBasis,NC,NDW,NX,NY,NZ,NT,verbose_level = verbose_level)

    make_adjoint_rep!(Uadj,Uin)

    return Uadj
end

function Base.setindex!(x::Adjoint_rep_Gaugefields_4D_wing,v,i1,i2,i3,i4,i5,i6)
    @inbounds x.U[i1,i2,i3 + x.NDW,i4 + x.NDW,i5 + x.NDW,i6 + x.NDW] = v
end

@inline function Base.getindex(x::Adjoint_rep_Gaugefields_4D_wing,i1,i2,i3,i4,i5,i6)
    @inbounds return x.U[i1,i2,i3 .+ x.NDW,i4 .+ x.NDW,i5 .+ x.NDW,i6 .+ x.NDW]
end

struct Shifted_Adjoint_rep_Gaugefields_4D{NC,NumofBasis} <: Shifted_Gaugefields{NC,4}
    parent::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis}
    #parent::T
    shift::NTuple{4,Int8}
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NDW::Int64
    outside::Bool

    #function Shifted_Gaugefields(U::T,shift,Dim) where {T <: AbstractGaugefields}
    function Shifted_Adjoint_rep_Gaugefields_4D(U::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},shift) where {NC,NumofBasis}
        outside = check_outside(U.NDW,shift)
        return new{NC,NumofBasis}(U,shift,U.NX,U.NY,U.NZ,U.NT,U.NDW,outside)
    end
end

function shift_U(U::Adjoint_rep_Gaugefields_4D_wing,ν::T) where {T <: Integer}
    if ν == 1
        shift = (1,0,0,0)
    elseif ν == 2
        shift = (0,1,0,0)
    elseif ν == 3
        shift = (0,0,1,0)
    elseif ν == 4
        shift = (0,0,0,1)
    elseif ν == -1
            shift = (-1,0,0,0)
    elseif ν == -2
            shift = (0,-1,0,0)
    elseif ν == -3
            shift = (0,0,-1,0)
    elseif ν == -4
            shift = (0,0,0,-1)
    end

    return Shifted_Adjoint_rep_Gaugefields_4D(U,shift)
end

function shift_U(U::TU,shift::NTuple{Dim,T}) where {Dim,T <: Integer,TU <: Adjoint_rep_Gaugefields_4D}
    return Shifted_Adjoint_rep_Gaugefields_4D(U,shift)
end

function Base.setindex!(U::Shifted_Adjoint_rep_Gaugefields_4D,v,i1,i2,i3,i4,i5,i6)  where NC
    error("type $(typeof(U)) has no setindex method. This type is read only.")
end

@inline function Base.getindex(U::Shifted_Adjoint_rep_Gaugefields_4D{NC,NumofBasis},i1,i2,i3,i4,i5,i6) where {NC,NumofBasis}
    if U.outside != true
        #println("inside shift: ",U.shift,"\t",U.outside)
        @inbounds return U.parent[i1,i2,i3 .+ U.shift[1],i4 .+ U.shift[2],i5 .+ U.shift[3],i6 .+ U.shift[4]]
    else
        #println("shift: ",U.shift,"\t",U.outside)

        i3_new = i3 + U.shift[1]
        i3_new += ifelse(i3_new > U.NX + U.NDW,-U.NX,0)
        i3_new += ifelse(i3_new < 1 - U.NDW,U.NX,0)
        i4_new = i4 + U.shift[2]
        i4_new += ifelse(i4_new > U.NY + U.NDW,-U.NY,0)
        i4_new += ifelse(i4_new < 1 - U.NDW,U.NY,0)
        i5_new = i5 + U.shift[3]
        i5_new += ifelse(i5_new > U.NZ + U.NDW,-U.NZ,0)
        i5_new += ifelse(i5_new < 1 - U.NDW,U.NZ,0)
        i6_new = i6 + U.shift[4]
        i6_new += ifelse(i6_new > U.NT + U.NDW,-U.NT,0)
        i6_new += ifelse(i6_new < 1 - U.NDW,U.NT,0)

        @inbounds return U.parent[i1,i2,i3_new ,i4_new ,i5_new ,i6_new ]
    end
    #function Base.getindex(U::Shifted_Gaugefields{T,4},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_4D_wing
end

function Base.similar(U::T) where T <: Adjoint_rep_Gaugefields_4D_wing
    Uout = Adjoint_rep_Gaugefields_4D_wing(U.NumofBasis,U.NC,U.NDW,U.NX,U.NY,U.NZ,U.NT,verbose_level = U.verbose_print.level)
    #identityGaugefields_4D_wing(U.NC,U.NX,U.NY,U.NZ,U.NT,U.NDW)
    return Uout
end

function Base.similar(U::Array{T,1}) where T <: Adjoint_rep_Gaugefields_4D_wing
    Uout = Array{T,1}(undef,4)
    for μ=1:4
        Uout[μ] = similar(U[μ])
    end
    return Uout
end

"""
Uadj = (1/2)*tr(Ta*U*Tb*Udag)
"""
function make_adjoint_rep!(Uadj::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},Uin) where {NC,NumofBasis}
    generators = Uadj.generators
    Ulocal = zeros(ComplexF64,NC,NC)
    NT = Uadj.NT
    NZ = Uadj.NZ
    NY = Uadj.NY
    NX = Uadj.NX

    @inbounds for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    for k2=1:NC
                        for k1=1:NC
                            Ulocal[k1,k2] = Uin[k1,k2,ix,iy,iz,it]
                        end
                    end

                    for b=1:NumofBasis
                        Tb = generators[b]
                        for a=1:NumofBasis
                            Ta = generators[a]
                            Uadj[a,b,ix,iy,iz,it] = tr(Ta*Ulocal*Tb*Ulocal')/2
                        end
                    end

                end
            end
        end
    end

end

function substitute_U!(a::Vector{T1},b::Vector{T2}) where {NC, NumofBasis,T1 <: Adjoint_rep_Gaugefields_4D_wing,T2 <: Abstractfields}
#function substitute_U!(a::Array{<: Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},1},b::Array{<: Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},1}) where {NC,NumofBasis}
    for i=1:4
        substitute_U!(a[i],b[i])
    end
end

function substitute_U!(a::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},b::T2) where {NC, NumofBasis,T2 <: Abstractfields}
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX
    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    for k2=1:NumofBasis
                        for k1=1:NumofBasis
                            @inbounds a[k1,k2,ix,iy,iz,it] = b[k1,k2,ix,iy,iz,it]
                        end
                    end
                end
            end
        end
    end
    set_wing_U!(a)

end

function LinearAlgebra.mul!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},a::T1,b::T2) where {NC,NumofBasis,T1 <: Number,T2 <: Abstractfields}
    @inbounds for i=1:length(c)
        c.U[i] = a*b.U[i]
    end
    return
end

function LinearAlgebra.mul!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},a::T1,b::T2,α::Ta,β::Tb) where {NC,NumofBasis,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
    NT = c.NT
    NZ = c.NZ
    NY = c.NY
    NX = c.NX
    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    for k2=1:NumofBasis
                        for k1=1:NumofBasis
                            c[k1,k2,ix,iy,iz,it] = β*c[k1,k2,ix,iy,iz,it]
                            @simd for k3=1:NumofBasis
                                c[k1,k2,ix,iy,iz,it] += α*a[k1,k3,ix,iy,iz,it]*b[k3,k2,ix,iy,iz,it]
                            end
                        end
                    end
                end
            end
        end
    end
    set_wing_U!(c)
end

function LinearAlgebra.mul!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},a::T1,b::T2,iseven::Bool) where {NC,T1 <: Abstractfields,T2 <: Abstractfields,NumofBasis}
    @assert NC != 2 && NC != 3 "This function is for NC != 2,3"
    NT = c.NT
    NZ = c.NZ
    NY = c.NY
    NX = c.NX
    @inbounds for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                    if evenodd == iseven
                        for k2=1:NumofBasis
                            for k1=1:NumofBasis
                                c[k1,k2,ix,iy,iz,it] = 0
                                @simd for k3=1:NumofBasis
                                    c[k1,k2,ix,iy,iz,it] += a[k1,k3,ix,iy,iz,it]*b[k3,k2,ix,iy,iz,it]
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    set_wing_U!(c)
end

function LinearAlgebra.tr(a::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},b::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis}) where {NC,NumofBasis}
    NX=a.NX
    NY=a.NY
    NZ=a.NZ
    NT=a.NT

    s = 0
    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    for k=1:NumofBasis
                        for k2=1:NumofBasis
                            s += a[k,k2,ix,iy,iz,it]*b[k2,k,ix,iy,iz,it]
                        end
                    end
                end
            end
        end
    end
    #println(3*NT*NZ*NY*NX*NC)
    return s
end

function LinearAlgebra.tr(a::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis}) where {NC,NumofBasis}
    NX=a.NX
    NY=a.NY
    NZ=a.NZ
    NT=a.NT

    s = 0
    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NumofBasis
                        s += a[k,k,ix,iy,iz,it]
                        #println(a[k,k,ix,iy,iz,it])
                    end
                end
            end
        end
    end
    #println(3*NT*NZ*NY*NX*NC)
    return s
end

function partial_tr(a::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},μ) where {NC,NumofBasis}
    NX=a.NX
    NY=a.NY
    NZ=a.NZ
    NT=a.NT

    if μ == 1
        s = 0
        ix = 1
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    #for ix=1:NX
                        @simd for k=1:NumofBasis
                            s += a[k,k,ix,iy,iz,it]
                            #println(a[k,k,ix,iy,iz,it])
                        end

                    #end
                end
            end
        end
    elseif μ == 2
        s = 0
        iy =1
        for it=1:NT
            for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @simd for k=1:NumofBasis
                            s += a[k,k,ix,iy,iz,it]
                            #println(a[k,k,ix,iy,iz,it])
                        end
                    end
                #end
            end
        end
    elseif μ == 3
        s = 0
        iz = 1
        for it=1:NT
            #for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        @simd for k=1:NumofBasis
                            s += a[k,k,ix,iy,iz,it]
                            #println(a[k,k,ix,iy,iz,it])
                        end
                    end
                end
            #end
        end
    else
        s = 0
        it = 1
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        @simd for k=1:NumofBasis
                            s += a[k,k,ix,iy,iz,it]
                            #println(a[k,k,ix,iy,iz,it])
                        end
                    end
                end
            end

    end

    #println(3*NT*NZ*NY*NX*NC)
    return s
end

function normalize_U!(u::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis}) where {NC,NumofBasis}
    NX = u.NX
    NY = u.NY
    NZ = u.NZ
    NT = u.NT

    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                @simd for ix=1:NX
                    A = u[:,:,ix,iy,iz,it]
                    gramschmidt!(A)
                    u[:,:,ix,iy,iz,it] = A[:,:]
                end
            end
        end
    end

end

function add_U!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},a::T1) where {NC,T1 <: Abstractfields,NumofBasis}
    @inbounds for i=1:length(c.U)
        c.U[i] += a.U[i]
    end
    return
end

function add_U!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},a::T1,iseven::Bool) where {NC,T1 <: Abstractfields,NumofBasis}
    NT = c.NT
    NZ = c.NZ
    NY = c.NY
    NX = c.NX
    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                    if evenodd == iseven
                        for k2=1:NumofBasis
                            @simd for k1=1:NumofBasis
                                c[k1,k2,ix,iy,iz,it] += a[k1,k2,ix,iy,iz,it]
                            end
                        end
                    end
                end
            end
        end
    end
end

function add_U!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},α::N,a::T1) where {NC,T1 <: Abstractfields, N<:Number,NumofBasis}
    #@inbounds for i=1:length(c.U)
    #    c.U[i] += α*a.U[i]
    #end
    #return

    NT = c.NT
    NZ = c.NZ
    NY = c.NY
    NX = c.NX
    @inbounds for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    for k2=1:NumofBasis
                        @simd for k1=1:NumofBasis
                            c[k1,k2,ix,iy,iz,it] += α*a[k1,k2,ix,iy,iz,it]
                        end
                    end
                end
            end
        end
    end
    set_wing_U!(c)
end

function set_wing_U!(u::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis}) where {NC,NumofBasis}
    NT = u.NT
    NY = u.NY
    NZ = u.NZ
    NX = u.NX
    NDW = u.NDW

    #X direction
    #Now we send data

    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for id=1:NDW
                    for k2=1:NumofBasis
                        @simd for k1=1:NumofBasis
                            @inbounds u[k1,k2,-NDW+id,iy,iz,it] = u[k1,k2,NX+(id-NDW),iy,iz,it]
                        end
                    end
                end
            end
        end
    end

    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for id=1:NDW
                    for k2=1:NumofBasis
                        @simd for k1=1:NumofBasis
                            @inbounds u[k1,k2,NX+id,iy,iz,it] = u[k1,k2,id,iy,iz,it]
                        end
                    end
                end
            end
        end
    end

    #Y direction
    #Now we send data
    for it=1:NT
        for iz=1:NZ
            for ix=-NDW+1:NX+NDW
                for id=1:NDW
                    for k1=1:NumofBasis
                        @simd for k2=1:NumofBasis
                            @inbounds u[k1,k2,ix,-NDW+id,iz,it] = u[k1,k2,ix,NY+(id-NDW),iz,it]
                        end
                    end
                end
            end
        end
    end

    for it=1:NT
        for iz=1:NZ
            for ix=-NDW+1:NX+NDW
                for id=1:NDW
                    for k1=1:NumofBasis
                        @simd for k2=1:NumofBasis
                            @inbounds u[k1,k2,ix,NY+id,iz,it] = u[k1,k2,ix,id,iz,it]
                        end
                    end
                end
            end
        end
    end

    #Z direction
    #Now we send data
    for id=1:NDW
        for it=1:NT
            for iy=-NDW+1:NY+NDW
                for ix=-NDW+1:NX+NDW
                    for k1=1:NumofBasis
                        @simd for k2=1:NumofBasis
                            @inbounds u[k1,k2,ix,iy,id-NDW,it] = u[k1,k2,ix,iy,NZ+(id-NDW),it]
                            @inbounds u[k1,k2,ix,iy,NZ+id,it] = u[k1,k2,ix,iy,id,it]
                        end
                    end
                end
            end
        end
    end

    for id=1:NDW
        for iz=-NDW+1:NZ+NDW
            for iy=-NDW+1:NY+NDW
                for ix=-NDW+1:NX+NDW
                    for k1=1:NumofBasis
                        @simd for k2=1:NumofBasis
                            @inbounds u[k1,k2,ix,iy,iz,id-NDW] = u[k1,k2,ix,iy,iz,NT+(id-NDW)]
                            @inbounds u[k1,k2,ix,iy,iz,NT+id] = u[k1,k2,ix,iy,iz,id]
                        end
                    end
                end
            end
        end
    end

    #display(u.g)
    #exit()

    return
end

function mul_skiplastindex!(c::Adjoint_rep_Gaugefields_4D_wing{NC,NumofBasis},a::T1,b::T2) where {NC,NumofBasis,T1 <: Abstractfields,T2 <: Abstractfields}
    #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
    NT = c.NT
    NZ = c.NZ
    NY = c.NY
    NX = c.NX
    #for it=1:NT
    it = 1
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    for k2=1:NumofBasis
                        for k1=1:NumofBasis
                            c[k1,k2,ix,iy,iz,it] = 0

                            @simd for k3=1:NumofBasis
                                c[k1,k2,ix,iy,iz,it] += a[k1,k3,ix,iy,iz,it]*b[k3,k2,ix,iy,iz,it]
                            end
                        end
                    end
                end
            end
        end
    #end
end


text:
type Source
     iz :: Int64
     ix :: Int64
     ot :: Float64
     dt :: Float64
     nt :: Int64
     p  :: Array{Float64,1}
end

function Ricker(f0::Float64, dt::Float64)
	  nw = 2.2/f0/dt
	  nw = 2*floor(Int,nw/2)+1
	  nc = floor(Int,nw/2)
	  w  = zeros(nw)
	  k  = collect(1:nw)
    k  = vec(k)
	  alpha = (nc-k+1)*f0*dt*pi
	  beta = alpha.^2
	  w = (1.-beta.*2).*exp(-beta)
	  return w
end

function InitSources(pos::Array{Int64, 2}, f0::Float64, ot::Array{Float64, 1}, dt::Float64)
    nsrc = size(pos, 1)
    if nsrc != length(ot)
       error("length of position and ot dismatch")
    end
    src = Array(Source, nsrc)
    for isrc = 1 : nsrc
        p = Ricker(f0, dt)
        nt= length(p)
        src[isrc] = Source(pos[isrc,1], pos[isrc,2], ot[isrc], dt, nt, p)
    end
    return src
end

function AddSources!(spt::SnapShot, src::Array{Source,1})
    dt = spt.dt
    if spt.dt != src[1].dt
       error("time sampling dismatch")
    end
    nz = spt.nz
    nx = spt.nx
    ext= spt.ext
    iflag = spt.iflag
    if iflag == 1
       zupper = ext
       Nz  = nz + 2*ext
       Nx  = nx + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz  = nz +   ext
       Nx  = nx + 2*ext
    end
    it  = spt.it
    nsrc = length(src)
    for isrc = 1: nsrc
        if src[isrc].ot <= (it-1)*dt <= src[isrc].ot+(src[isrc].nt-1)*dt
           indt = it - round(Int64, src[isrc].ot/dt)
           indz = src[isrc].iz + zupper
           indx = src[isrc].ix + ext
           pos = (indx-1) * Nz + indz
           spt.pz[pos] = spt.pz[pos] + src[isrc].p[indt]/2
           spt.px[pos] = spt.px[pos] + src[isrc].p[indt]/2
        end
    end
    return nothing
end

function AddSources!(spt::SnapShot, w::Array{Float64,1}, dis::Array{Float64,1})
    it = spt.it
    if it > length(w)
       error("it is out the range of wavelet")
    end
    nz = spt.nz ; nx = spt.nx;
    ext= spt.ext; iflag = spt.iflag;
    if length(dis) != nz*nx
       error("size of dis dismatch with the size of SnapShot")
    end
    if iflag == 1
       zupper = ext
       Nz  = nz + 2*ext
       Nx  = nx + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz  = nz +   ext
       Nx  = nx + 2*ext
    end
    tmp = zeros(Nz, Nx); dis_tmp = reshape(dis, nz, nx)
    tmp[zupper+1: zupper+nz, ext+1:ext+nx] = 1/2 * w[it] * dis_tmp
    tmp = vec(tmp)
    spt.pz[:] = spt.pz[:] + tmp[:]
    spt.px[:] = spt.px[:] + tmp[:]
    return nothing
end

function AddSource!(spt::SnapShot, src::Source)
    dt = spt.dt
    if spt.dt != src.dt
       error("time sampling dismatch")
    end
    nz = spt.nz
    nx = spt.nx
    ext= spt.ext
    iflag = spt.iflag
    if iflag == 1
       zupper = ext
       Nz  = nz + 2*ext
       Nx  = nx + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz  = nz +   ext
       Nx  = nx + 2*ext
    end
    it  = spt.it
    if src.ot <= (it-1)*dt <= src.ot+(src.nt-1)*dt
       indt = it - round(Int64, src.ot/dt)
       indz = src.iz + zupper
       indx = src.ix + ext
       pos = (indx-1) * Nz + indz
      #  spt.vz[pos] = spt.vz[pos] + src.p[indt]
       spt.pz[pos] = spt.pz[pos] + (src.p[indt])/2
       spt.px[pos] = spt.px[pos] + (src.p[indt])/2
    end
    return nothing
end

function SrcRange(src::Array{Source, 1})
    tmax = 0.0
    tmin = 0.0
    nsrc = length(src)
    dt   = src[1].dt
    for isrc = 1: nsrc
        tmp = src[isrc].ot
        tmp1= tmp + (src[isrc].nt-1)*dt
        if tmp < tmin
           tmin = tmp
        end
        if tmp1 > tmax
           tmax = tmp1
        end
    end
    return tmin, tmax
end

function src2spt(src::Array{Source,1}, it::Int64, nz::Int64, nx::Int64, ext::Int64, iflag::Int64)
    dt = src[1].dt
    spt = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, it)
    AddSources!(spt, src)
    return spt
end

function src2SnapShots(path::ASCIIString, src::Array{Source,1}, nz::Int64, nx::Int64, ext::Int64, iflag::Int64)
    dt = src[1].dt
    (tl, tu) = SrcRange(src)
    nt  = round(Int64, (tu-tl)/dt) + 1
    fid = open(path, "w")
    write(fid, src[1].iz, src[1].ix, nz, nx, ext, iflag, tl, dt)
    for it = 1 : nt
        it1 = round(Int64, (tl+(it-1)*dt)/dt+1)
        spt = src2spt(src, it1, nz, nx, ext, iflag)
        writeSnapShot(fid, spt)
    end
    close(fid)
    return nothing
end

type SrcSpt
     nz :: Int64
     nx :: Int64
     ot :: Float64
     dt :: Float64
     it :: Int64
     p  :: Array{Float64, 1}
end

function InitSrcSpt(nz::Int64, nx::Int64, ot::Float64, dt::Float64, it::Int64)
    p = zeros(nz*nx)
    srcspt = SrcSpt(nz, nx, ot, dt, it, p)
    return srcspt
end

function src2srcSpts(path::ASCIIString, src::Array{Source,1}, nz::Int64, nx::Int64)
    ns = length(src)
    dt = src[1].dt
    (tl, tu) = SrcRange(src)
    nt  = round(Int64, (tu-tl)/dt) + 1
    fid = open(path, "w")
    write(fid, nz, nx, tl, src[1].dt, nt)
    for it = 1 : nt
        p   = zeros(nz*nx)
        for is = 1 : ns
            if src[is].ot <= tl+(it-1)*dt <= src[is].ot + dt*(src[is].nt-1)
               indt = it - round(Int64, (src[is].ot-tl)/dt)
               ind = (src[is].ix-1)*nz + src[is].iz
               p[ind] = p[ind] + src[is].p[indt]
            end
        end
        write(fid, p)
    end
    close(fid)
    return nothing
end

function readSrcSpt(path::ASCIIString, it::Int64)
    fid = open(path, "r")
    nz = read(fid, Int64)  ; nx = read(fid, Int64)  ;
    ot = read(fid, Float64); dt = read(fid, Float64);
    position = sizeof(Int64)*5 + sizeof(Float64)*nz*nx*(it-1)
    seek(fid, position)
    p = read(fid, Float64, nz*nx)
    p = reshape(p, nz, nx)
    close(fid)
    return p
end

function srcSpt2wlet(path::ASCIIString, iz::Int64, ix::Int64)
    fid = open(path, "r")
    (nz, nx, ot, dt, nt) = InfoSrcSpt(path)
    w = zeros(nt)
    for it = 1 : nt
        p = readSrcSpt(path, it)
        w[it] = p[iz, ix]
    end
    return w
end

function InfoSrcSpt(path::ASCIIString)
    fid = open(path, "r")
    nz = read(fid, Int64)
    nx = read(fid, Int64)
    ot = read(fid, Float64)
    dt = read(fid, Float64)
    nt = read(fid, Int64)
    close(fid)
    return nz, nx, ot, dt, nt
end

function AddSrcSpt2Spt(spt::SnapShot, path::ASCIIString)
    it = spt.it
    nz = spt.nz; nx = spt.nx;
    ext = spt.ext; iflag = spt.iflag;
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    p = readSrcSpt(path, it) * 1/2
    tmp = reshape(spt.pz, Nz, Nx)
    tmp[zupper+1:nz+zupper, ext+1:ext+nx] = tmp[zupper+1:nz+zupper, ext+1:ext+nx] + p
    spt.pz = vec(tmp)
    tmp = reshape(spt.px, Nz, Nx)
    tmp[zupper+1:nz+zupper, ext+1:ext+nx] = tmp[zupper+1:nz+zupper, ext+1:ext+nx] + p
    spt.px = vec(tmp)
    return nothing
end

function saveSpt2SrcSpt(path::ASCIIString, spt::SnapShot, nz::Int64, nx::Int64, nt::Int64, ext::Int64, Nz::Int64, Nx::Int64, zupper::Int64)
    fid = open(path, "w")
    write(fid, nz, nx, 0.0, spt.dt, nt)
    tmp = 1/2*(reshape(spt.pz, Nz, Nx) + reshape(spt.px, Nz, Nx))
    tmp = vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])
    write(fid, tmp)
    flush(fid)
    return fid
end

function saveSpt2SrcSpt(fid::IOStream, spt::SnapShot, nz::Int64, nx::Int64, ext::Int64, Nz::Int64, Nx::Int64, zupper::Int64)
    tmp = 1/2*(reshape(spt.pz, Nz, Nx) + reshape(spt.px, Nz, Nx))
    tmp = vec(tmp[zupper+1:zupper+nz, ext+1:ext+nx])
    write(fid, tmp)
    flush(fid)
    return nothing
end

function reverseSrcSptOrder(path::ASCIIString, path_tmp::ASCIIString)
    (nz, nx, ot, dt, nt) = InfoSrcSpt(path_tmp)
    fid = open(path, "w")
    write(fid, nz, nx, ot, dt, nt)
    fid_tmp = open(path_tmp, "r")
    for it = nt : -1 : 1
        position = sizeof(Float64)*5 + sizeof(Float64)*nz*nx*(it-1)
        seek(fid_tmp, position)
        d = read(fid_tmp, Float64, nz*nx);
        write(fid, d)
    end
    close(fid_tmp); close(fid)
    rm(path_tmp)
    return nothing
end

ext: jl

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
@doc raw"""
    GaussianLikelihood(σ²::T=1e-3) # σ² is the variance of the noise

Gaussian noise :
```math
    p(y|f) = N(y|f,\sigma^2)
```
There is no augmentation needed for this likelihood which is already conjugate to a Gaussian prior.
"""
struct GaussianLikelihood{T<:Real,O} <: RegressionLikelihood
    σ²::Vector{T}
    opt_noise::O
    function GaussianLikelihood{T}(σ²::T, opt_noise) where {T<:Real}
        return new{T,typeof(opt_noise)}([σ²], opt_noise)
    end
end

function GaussianLikelihood(σ²::T=1e-3; opt_noise=false) where {T<:Real}
    if isa(opt_noise, Bool)
        opt_noise = opt_noise ? ADAM(0.05) : nothing
    end
    return GaussianLikelihood{T}(σ², opt_noise)
end

implemented(::GaussianLikelihood, ::Union{<:AnalyticVI,<:Analytic}) = true

function (l::GaussianLikelihood)(y::Real, f::Real)
    return pdf(Normal(y, sqrt(noise(l))), f)
end

function Distributions.loglikelihood(l::GaussianLikelihood, y::Real, f::Real)
    return logpdf(Normal(y, sqrt(noise(l))), f)
end

noise(l::GaussianLikelihood) = only(l.σ²)

function Base.show(io::IO, l::GaussianLikelihood)
    return print(io, "Gaussian likelihood (σ² = $(noise(l)))")
end

function compute_proba(
    l::GaussianLikelihood, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
)
    return μ, σ² .+ noise(l)
end

function init_local_vars(l::GaussianLikelihood, batchsize::Int, T::DataType=Float64)
    local_vars = (; θ=fill(inv(l.σ²[1]), batchsize))
    if !isnothing(l.opt_noise)
        state_σ² = Optimisers.init(l.opt_noise, l.σ²)
        local_vars = merge(local_vars, (; state_σ²))
    end
    return local_vars
end

function local_updates!(
    local_vars,
    l::GaussianLikelihood,
    y::AbstractVector,
    μ::AbstractVector,
    var_f::AbstractVector,
)
    if !isnothing(l.opt_noise)
        grad = ((sum(abs2, y - μ) + sum(var_f)) / noise(l) - length(y)) / 2
        gradlog, local_vars.state_σ² = Optimisers.apply!(
            l.opt_noise, local_vars.state_σ², l.σ², [grad]
        )
        l.σ² .= exp.(log.(l.σ²) + gradlog)
    end
    local_vars.θ .= inv(noise(l))
    return local_vars
end

@inline function ∇E_μ(l::GaussianLikelihood, ::AOptimizer, y::AbstractVector, state)
    return (y ./ noise(l),)
end

@inline function ∇E_Σ(::GaussianLikelihood, ::AOptimizer, ::AbstractVector, state)
    return (state.θ / 2,)
end

function expec_loglikelihood(
    l::GaussianLikelihood,
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diagΣ::AbstractVector,
    state,
)
    return -(
        length(y) * (log(twoπ) + log(noise(l))) + (sum(abs2, y - μ) + sum(diagΣ)) / noise(l)
    ) / 2
end

AugmentedKL(::GaussianLikelihood, state, ::Any) = 0


text:
sir_model = @min_reaction_network rn begin
    0.1/1000, s + i --> 2i
    0.01, i --> r
end
addjumps!(sir_model)
sir_prob = DiscreteProblem([999,1,0],(0.0,250.0))
sir_jump_prob = JumpProblem(sir_prob,Direct(),sir_model)

sir_sol = solve(sir_jump_prob,FunctionMap())

using Plots; plot(sir_sol)

Random.seed!(1234)
nums = Int[]
@time for i in 1:100000
  sir_sol = solve(sir_jump_prob,FunctionMap())
  push!(nums,sir_sol[end][3])
end
println("Reaction DSL: $(mean(nums))")

rate = (u,p,t) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!)

rate = (u,p,t) -> 0.01u[2]
affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!)

prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,FunctionMap())

using Plots; plot(sol)

nums = Int[]
@time for i in 1:100000
  sol = solve(jump_prob,SSAStepper())
  push!(nums,sol[end][3])
end
println("Direct Jumps: $(mean(nums))")

using Gillespie

function F(x,parms)
  (S,I,R) = x
  (beta,gamma) = parms
  infection = beta*S*I
  recovery = gamma*I
  [infection,recovery]
end

x0 = [999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/1000.0,0.01]
tf = 250.0
Random.seed!(1234)

nums = Int[]
@time for i in 1:100000
  result = ssa(copy(x0),F,nu,parms,tf)
  data = ssa_data(result)
  push!(nums,data[:x3][end])
end
println("Gillespie: $(mean(nums))")


text:
"""
Definition of analysis data types for juliaPBF
    All values in AnalysisDataStruct shouldn't
    be changed after input parsing(inputParsin
    g.jl). All variables should be located in
    SimulationDataStruct.
"""

export
    TimeStep,
    AnalysisBox,
    InitialBox,
    FluidPropertyData,
    SolidPropertyData,
    PhaseData,
    TensileInstability

struct TimeStep
    type        ::Int64
    dt          ::Float64
    cfl         ::Float64
    dtMin       ::Float64
    dtMax       ::Float64
end

struct InitialBox
    phaseID     ::Int64
    box         ::Box2
    vel         ::Vec2
end

struct FluidPropertyData
    name        ::String
    size        ::Float64
    density     ::Float64
    viscosity   ::Float64
end

struct SolidPropertyData
    name        ::String
    size        ::Float64
    density     ::Float64
end

struct PhaseData
    type        ::String   # FluidPropertyArray or SolidPropertyArray?
    Fluid       ::FluidPropertyData
    Solid       ::SolidPropertyData
end

mutable struct TensileInstability
    k           ::Float64
    dq          ::Float64
    n           ::Float64
end

TimeStep() = TimeStep(0, 0.0, 0.0, 0.0, 0.0)
FluidPropertyData() = FluidPropertyData("",0,0,0)
SolidPropertyData() = SolidPropertyData("",0,0)
PhaseData() = PhaseData(0,FluidPropertyData(),SolidPropertyData())

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
module RNN

export RNN_t
export rnn_modify_synapse!,     rnn_remove_synapse!, rnn_add_synapse!
export rnn_add_neuron!,         rnn_remove_neuron!,  rnn_modify_bias!
export rnn_copy,                rnn_read,            rnn_write
export rnn_merge_two
export rnn_add_input_neuron!,   rnn_add_output_neuron!
export rnn_add_input_neurons!,  rnn_add_output_neurons!
export rnn_add_hidden_neurons!, rnn_add_hidden_neuron!
export rnn_update!,             rnn_control!

type RNNInputSizeMismatch <: Exception
  expected::Int64
  received::Int64
end

showerror(io::IO, e::RNNInputSizeMismatch) = print(io, "RNN Input Mismatch: expected ", e.expected, " values, but received ", e.received)

type RNN_t
  ni::Int64
  no::Int64
  nh::Int64
  n::Int64

  w::Matrix{Float64}
  b::Matrix{Float64}
  o::Matrix{Float64}
  a::Matrix{Float64}

  RNN_t(ni::Int64, no::Int64, nh::Int64) =
    new(ni, no, nh, ni+no+nh,
        zeros(ni+no+nh, ni+no+nh), zeros(ni+no+nh,1),
        zeros(ni+no+nh,1),         zeros(ni+no+nh,1))
end

function rnn_add_input_neurons!(rnn::RNN_t, n::Integer)
  a = rnn.w[1:rnn.ni,:]
  b = rnn.w[rnn.ni+1:end,:]
  c = vcat(a,zeros(n,size(a)[2]),b)

  a = c[:,1:rnn.ni]
  b = c[:,rnn.ni+1:end]
  c = hcat(a,zeros(size(a)[1],n),b)

  rnn.w  = c
  rnn.b  = vcat(rnn.b[1:rnn.ni], zeros(n,1), rnn.b[rnn.ni+1:end])
  rnn.o  = vcat(rnn.o[1:rnn.ni], zeros(n,1), rnn.o[rnn.ni+1:end])
  rnn.a  = vcat(rnn.a[1:rnn.ni], zeros(n,1), rnn.a[rnn.ni+1:end])
  rnn.ni = rnn.ni + n
  rnn.n  = rnn.n  + n
  rnn
end

function rnn_add_input_neuron!(rnn::RNN_t, n::Integer)
  rnn_add_input_neurons!(rnn, 1)
end

function rnn_add_output_neurons!(rnn::RNN_t, n::Integer)
  a = rnn.w[1:rnn.ni+rnn.no,:]
  b = rnn.w[rnn.ni+rnn.no+1:end,:]
  c = vcat(a,zeros(n,size(a)[2]),b)

  a = c[:,1:rnn.ni+rnn.no]
  b = c[:,rnn.ni+rnn.no+1:end]
  c = hcat(a,zeros(size(a)[1],n),b)

  rnn.w  = c
  rnn.b  = vcat(rnn.b[1:rnn.ni+rnn.no], zeros(n,1), rnn.b[rnn.ni+rnn.no+1:end])
  rnn.a  = vcat(rnn.a[1:rnn.ni+rnn.no], zeros(n,1), rnn.a[rnn.ni+rnn.no+1:end])
  rnn.o  = vcat(rnn.o[1:rnn.ni+rnn.no], zeros(n,1), rnn.o[rnn.ni+rnn.no+1:end])
  rnn.no = rnn.no + n
  rnn.n  = rnn.n  + n
  rnn
end

function rnn_add_output_neuron!(rnn::RNN_t)
  rnn_add_output_neurons!(rnn, 1)
end

function rnn_add_hidden_neurons!(rnn::RNN_t, n::Integer)
  rnn.w  = vcat(rnn.w,zeros(n,size(rnn.w)[2]))
  rnn.w  = hcat(rnn.w,zeros(size(rnn.w)[1],n))
  rnn.b  = vcat(rnn.b, zeros(n,1))
  rnn.o  = vcat(rnn.o, zeros(n,1))
  rnn.a  = vcat(rnn.a, zeros(n,1))
  rnn.nh = rnn.nh + n
  rnn.n  = rnn.n  + n
  rnn
end

function rnn_add_hidden_neuron!(rnn::RNN_t)
  rnn_add_hidden_neurons!(rnn, 1)
end

function rnn_merge_two(rnn1, rnn2)
  ni  = rnn1.ni + rnn2.ni
  no  = rnn1.no + rnn2.no
  nh  = rnn1.nh + rnn2.nh
  rnn = RNN_t(ni, no, nh)

  # rnn1 input to output
  r1_srcs = 1:rnn1.ni
  r1_dest = rnn1.ni+1:rnn1.ni+rnn1.no
  r_srcs  = 1:rnn1.ni
  r_dest  = ni+1:ni+rnn1.no
  rnn.w[r_dest, r_srcs] = rnn1.w[r1_dest, r1_srcs]

  if rnn1.nh > 0
    # rnn1 input to hidden
    r1_srcs = 1:rnn1.ni
    r1_dest = rnn1.ni+rnn1.no+1:rnn1.ni+rnn1.no+rnn1.nh
    r_srcs  = 1:rnn1.ni
    r_dest  = ni+no+1:ni+no+rnn1.nh
    rnn.w[r_dest, r_srcs] = rnn1.w[r1_dest, r1_srcs]

    # rnn1 hidden to hidden
    r1_srcs = rnn1.ni+rnn1.no+1:rnn1.ni+rnn1.no+rnn1.nh
    r1_dest = r1_srcs
    r_srcs  = rnn.ni+rnn.no+1:rnn.ni+rnn.no+rnn1.nh
    r_dest  = r_srcs
    rnn.w[r_dest, r_srcs] = rnn1.w[r1_dest, r1_srcs]

    # rnn1 output to hidden
    r1_srcs = rnn1.ni+1:rnn1.ni+rnn1.no
    r1_dest = rnn1.ni+rnn1.no+1:rnn1.ni+rnn1.no+rnn1.nh
    r_srcs  = rnn.ni+1:rnn.ni+rnn1.no
    r_dest  = rnn.ni+rnn.no+1:rnn.ni+rnn.no+rnn1.nh
    rnn.w[r_dest, r_srcs] = rnn1.w[r1_dest, r1_srcs]

    # rnn1 hidden to output
    r1_srcs = rnn1.ni+rnn1.no+1:rnn1.ni+rnn1.no+rnn1.nh
    r1_dest = rnn1.ni+1:rnn1.ni+rnn1.no
    r_srcs  = rnn.ni+rnn.no+1:rnn.ni+rnn.no+rnn1.nh
    r_dest  = rnn.ni+1:rnn.ni+rnn1.no
    rnn.w[r_dest, r_srcs] = rnn1.w[r1_dest, r1_srcs]

  end

  #= # rnn1 output to output =#
  r1_srcs = rnn1.ni+1:rnn1.ni+rnn1.no
  r1_dest = r1_srcs
  r_srcs  = rnn.ni+1:rnn.ni+rnn1.no
  r_dest  = r_srcs
  rnn.w[r_dest, r_srcs] = rnn1.w[r1_dest, r1_srcs]

  #
  # RNN2
  #

  # rnn2 input to output
  r2_srcs = 1:rnn2.ni
  r2_dest = rnn2.ni+1:rnn2.ni+rnn2.no
  r_srcs  = rnn1.ni+1:ni
  r_dest  = ni+rnn1.no+1:ni+no
  rnn.w[r_dest, r_srcs] = rnn2.w[r2_dest, r2_srcs]

  if rnn2.nh > 0
    # rnn2 input to hidden
    r2_srcs = 1:rnn2.ni
    r2_dest = rnn2.ni+rnn2.no+1:rnn2.ni+rnn2.no+rnn2.nh
    r_srcs  = rnn1.ni+1:rnn.ni
    r_dest  = rnn.ni+rnn.no+rnn1.nh+1:rnn.ni+rnn.no+rnn.nh
    rnn.w[r_dest, r_srcs] = rnn2.w[r2_dest, r2_srcs]

    # rnn2 hidden to hidden
    r2_srcs = rnn2.ni+rnn2.no+1:rnn2.ni+rnn2.no+rnn2.nh
    r2_dest = r2_srcs
    r_srcs  = rnn.ni+rnn.no+rnn1.nh+1:rnn.ni+rnn.no+rnn.nh
    r_dest  = r_srcs
    rnn.w[r_dest, r_srcs] = rnn2.w[r2_dest, r2_srcs]

    # rnn2 output to hidden
    r2_srcs = rnn2.ni+1:rnn2.ni+rnn2.no
    r2_dest = rnn2.ni+rnn2.no+1:rnn2.ni+rnn2.no+rnn2.nh
    r_srcs  = rnn.ni+rnn1.no+1:rnn.ni+rnn.no
    r_dest  = rnn.ni+rnn.no+rnn1.nh+1:rnn.ni+rnn.no+rnn.nh
    rnn.w[r_dest, r_srcs] = rnn2.w[r2_dest, r2_srcs]

    # rnn2 hidden to output
    r2_srcs = rnn2.ni+rnn2.no+1:rnn2.ni+rnn2.no+rnn2.nh
    r2_dest = rnn2.ni+1:rnn2.ni+rnn2.no
    r_srcs  = rnn.ni+rnn.no+rnn1.nh+1:rnn.ni+rnn.no+rnn.nh
    r_dest  = rnn.ni+rnn1.no+1:rnn.ni+rnn.no
    rnn.w[r_dest, r_srcs] = rnn2.w[r2_dest, r2_srcs]

  end

  # rnn2 output to output
  r2_srcs = rnn2.ni+1:rnn2.ni+rnn2.no
  r2_dest = r2_srcs
  r_srcs  = rnn.ni+rnn1.no+1:rnn.ni+rnn.no
  r_dest  = r_srcs
  rnn.w[r_dest, r_srcs] = rnn2.w[r2_dest, r2_srcs]

  rnn.b[1:rnn1.ni]                        = rnn1.b[1:rnn1.ni]
  rnn.b[rnn1.ni+1:rnn.ni]                 = rnn2.b[1:rnn2.ni]
  rnn.b[rnn.ni+1:rnn.ni + rnn1.no]        = rnn1.b[rnn1.ni+1:rnn1.ni+rnn1.no]
  rnn.b[rnn.ni+rnn1.no+1:rnn.ni + rnn.no] = rnn2.b[rnn2.ni+1:rnn2.ni+rnn2.no]
  if rnn1.nh > 0
    rnn.b[rnn.ni+rnn.no+1:rnn.ni+rnn.no+rnn1.nh] = rnn1.b[rnn1.ni+rnn1.no+1:rnn1.n]
  end
  if rnn2.nh > 0
    rnn.b[rnn.ni+rnn.no+rnn1.nh+1:rnn.n]         = rnn2.b[rnn2.ni+rnn2.no+1:rnn2.n]
  end

  rnn
end

function rnn_write(rnn, output)
  f = open(output, "w")
  # write parameters
  @printf(f,"%d %d %d\n", rnn.ni, rnn.no, rnn.nh)
  # writing W
  for i in 1:rnn.n
    @printf(f,"%f",rnn.w[i,1])
    for j in 2:rnn.n
      @printf(f," %f",rnn.w[i,j])
    end
    @printf(f,"\n")
  end

  # writing b
  @printf(f,"%f",rnn.b[1])
  for n in 2:rnn.n
    @printf(f," %f",rnn.b[n])
  end
  @printf(f,"\n")
  close(f)
end

function rnn_read(input)
  f   = open(input, "r")
  sizes = split(strip(readline(f)))
  ni = int(eval(sizes[1]))
  no = int(eval(sizes[2]))
  nh = int(eval(sizes[3]))
  n = ni + no + nh
  rnn = RNN_t(ni, no, nh)
  # write parameters
  # writing W
  for i=1:n
    values = split(strip(readline(f)))
    for j=1:n
      rnn.w[i,j] = float64(eval(parse(values[j])))
    end
  end

  s = split(strip(readline(f)))
  for j=1:n
    rnn.b[j] = float64(eval(s[j]))
  end

  close(f)
  return rnn
end

function rnn_copy(rnn::RNN_t)
  new_rnn   = RNN_t(rnn.ni, rnn.no, rnn.nh)
  new_rnn.w = deepcopy(rnn.w)
  new_rnn.b = deepcopy(rnn.b)
  new_rnn.o = deepcopy(rnn.o)
  new_rnn.a = deepcopy(rnn.a)
  return new_rnn
end

function rnn_update!(rnn::RNN_t)
  rnn.a = rnn.b + rnn.w * rnn.o
  for i=1:rnn.ni
    rnn.o[i] = rnn.a[i]
  end
  for i=rnn.ni+1:rnn.n
    rnn.o[i]=tanh(rnn.a[i])
  end
end

#
# start - synapse functions
#
function rnn_remove_synapse!(rnn::RNN_t, prob::Float64)
  for i=1:rnn.n
    for j=1:rnn.n
      if rnn.w[i,j] != 0.0 && rand() <= prob
        rnn.w[i,j] = 0.0
      end
    end
  end
end

function rnn_add_synapse!(rnn::RNN_t, prob::Float64, delta::Float64)
  for i=rnn.ni+1:rnn.n
    for j=1:rnn.n
      if rnn.w[i,j] == 0.0 && rand() <= prob
        rnn.w[i,j] = (2.0 * rand() - 1.0) * delta
      end
    end
  end
end

function rnn_modify_synapse!(rnn::RNN_t, prob::Float64, delta::Float64, maximum::Float64)
  for i=rnn.ni+1:rnn.n
    for j=1:rnn.n
      if rnn.w[i,j] != 0.0
        if rand() <= prob
          rnn.w[i,j] += (2.0 * rand() - 1.0) * delta
          if abs(rnn.w[i,j]) > maximum
            if rnn.w[i,j] < 0
              rnn.w[i,j] = -maximum
            else
              rnn.w[i,j] =  maximum
            end
          end
        end
      end
    end
  end
end

#
# end - synapse functions
#

function rnn_modify_bias!(rnn::RNN_t, prob::Float64,
                                      delta::Float64,
                                      maximum::Float64)
  for i=rnn.ni+1:rnn.n
    if rand() <= prob
      rnn.b[i] = (2.0 * rand() - 1.0) * delta
      if abs(rnn.b[i]) > maximum
        if rnn.b[i] < 0
          rnn.b[i] = -maximum
        else
          rnn.b[i] =  maximum
        end
      end
    end
  end
end

function rnn_add_neuron!(rnn::RNN_t, prob::Float64,    delta_b::Float64,
                                     delta_w::Float64, prob_c::Float64,
                                     max_nh::Int64)
  if rnn.nh >= max_nh
    return
  end
  if rand() <= prob
    rnn.nh += 1
    rnn.n  += 1

    rnn.a        = cat(1, rnn.a, zeros(1))
    rnn.b        = cat(1, rnn.b, zeros(1))
    rnn.o        = cat(1, rnn.o, zeros(1))

    rnn.w        = cat(1, rnn.w, zeros(1,size(rnn.w)[2]))
    rnn.w        = cat(2, rnn.w, zeros(size(rnn.w)[1],1))

    rnn.b[rnn.n] = (2.0 * rand() - 1.0) * delta_b

    for i = rnn.ni+1:rnn.n
      if rand() <= prob_c
        rnn.w[i,rnn.n] = (2.0 * rand() - 1.0) * delta_w
      end
    end
    for i = 1:rnn.n
      if rand() <= prob_c
        rnn.w[rnn.n,i] = (2.0 * rand() - 1.0) * delta_w
      end
    end
  end
end

function rnn_remove_neuron!(rnn::RNN_t, prob::Float64)
  if rnn.nh > 0
    r = [1:rnn.n]
    for i=rnn.nh:-1:1
      if rand() <= prob
        splice!(r,rnn.ni+rnn.no+i)
      end
    end
    r = sort(r)
    rnn.w  = rnn.w[r,:]
    rnn.w  = rnn.w[:,r]
    rnn.a  = rnn.a[r,:]
    rnn.o  = rnn.o[r,:]
    rnn.b  = rnn.b[r,:]
    rnn.n  = length(r)
    rnn.nh = length(r) - rnn.ni - rnn.no
  end
end

function rnn_control!(rnn::RNN_t, input::Vector{Float64})
  if rnn.ni != length(input)
    error("RNN Input Mismatch: expected ", rnn.ni, " values, but received ", length(input))
    exit(-1)
  end
  for i=1:rnn.ni
    rnn.b[i] = input[i]
  end
  rnn_update!(rnn)
  r   = zeros(rnn.no)
  for i=1:rnn.no
    r[i]=rnn.o[rnn.ni+i]
  end
  return r
end
control!(rnn::RNN_t, input::Vector{Float64}) = rnn_control!(rnn, input)

end # module


text:
export RotatedHyperrectangle

"""
    RotatedHyperrectangle{N, MN<:AbstractMatrix{N},
                          HT<:AbstractHyperrectangle{N}} <: AbstractZonotope{N}

Type that represents a hyperrectangle that is not necessarily axis aligned.

### Fields

- `M`   -- matrix (not necessarily invertible)
- `box` -- axis-aligned hyperrectangle

### Notes

The matrix `M` is typically a rotation matrix and hence invertible, but this
type does not require that and any matrix that is compatible with the dimension
of `box` is allowed.
"""
struct RotatedHyperrectangle{N, MN<:AbstractMatrix{N},
                             HT<:AbstractHyperrectangle{N}} <: AbstractZonotope{N}
    M::MN
    box::HT

    function RotatedHyperrectangle(M::MN, box::HT) where {N,
            MN<:AbstractMatrix{N}, HT<:AbstractHyperrectangle{N}}
        @assert size(M, 2) == dim(box) "a rotated hyperrectangle of dimension " *
            "$(dim(box)) is incompatible with a matrix of dimension $(size(M))"
        return new{N, MN, HT}(M, box)
    end
end

isoperationtype(::Type{RotatedHyperrectangle}) = false
isconvextype(::Type{RotatedHyperrectangle}) = true

"""
    dim(R::RotatedHyperrectangle)

Return the ambient dimension of a rotated hyperrectangle.

### Input

- `R` -- rotated hyperrectangle

### Output

The ambient dimension of the rotated hyperrectangle.
"""
function dim(R::RotatedHyperrectangle)
    return size(R.M, 1)
end

"""
    ρ(d::AbstractVector, R::RotatedHyperrectangle)

Evaluate the support function of a rotated hyperrectangle in a given direction.

### Input

- `d` -- direction
- `R` -- rotated hyperrectangle

### Output

Evaluation of the support function in the given direction.
"""
function ρ(d::AbstractVector, R::RotatedHyperrectangle)
    return _ρ_linear_map(d, R.M, R.box)
end

"""
    σ(d::AbstractVector, R::RotatedHyperrectangle)

Return the support vector of a rotated hyperrectangle in a given direction.

### Input

- `d` -- direction
- `R` -- rotated hyperrectangle

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector, R::RotatedHyperrectangle)
    return _σ_linear_map(d, R.M, R.box)
end

"""
    center(R::RotatedHyperrectangle)

Return the center of a rotated hyperrectangle.

### Input

- `R` -- rotated hyperrectangle

### Output

The center of the rotated hyperrectangle.
"""
function center(R::RotatedHyperrectangle)
    return R.M * center(R.box)
end

"""
    generators(R::RotatedHyperrectangle)

Return an iterator over the generators of a rotated hyperrectangle.

### Input

- `R` -- rotated hyperrectangle

### Output

An iterator over the generators of `R`.
"""
function generators(R::RotatedHyperrectangle)
    return generators_fallback(R)
end

"""
   genmat(R::RotatedHyperrectangle)

Return the generator matrix of a rotated hyperrectangle.

### Input

- `R` -- rotated hyperrectangle

### Output

A matrix where each column represents one generator of `R`.
"""
function genmat(R::RotatedHyperrectangle)
    return R.M * genmat(R.box)
end

"""
    linear_map(M::AbstractMatrix, R::RotatedHyperrectangle)

Compute the concrete linear map of a rotated hyperrectangle.

### Input

- `M` -- matrix
- `R` -- rotated hyperrectangle

### Output

A new rotated hyperrectangle.

### Notes

If `M` is not a square matrix, the result is not necessarily a rotated
hyperrectangle.
The type still represents the correct set, but it is suggested to only apply
square linear maps to it.
"""
function linear_map(M::AbstractMatrix, R::RotatedHyperrectangle)
    @assert size(M, 2) == dim(R) "a linear map of size $(size(M)) is " *
        "incompatible with a set of dimension $(dim(R))"

    return RotatedHyperrectangle(M * R.M, R.box)
end

"""
    vertices_list(R::RotatedHyperrectangle)

Return the list of vertices of a rotated hyperrectangle.

### Input

- `R` -- rotated hyperrectangle

### Output

A list of the vertices.
"""
function vertices_list(R::RotatedHyperrectangle)
    return broadcast(v -> R.M * v, vertices_list(R.box))
end

"""
    constraints_list(R::RotatedHyperrectangle)

Return the list of constraints of a rotated hyperrectangle.

### Input

- `R` -- rotated hyperrectangle

### Output

A list of the linear constraints.
"""
function constraints_list(R::RotatedHyperrectangle)
    return constraints_list(convert(VPolytope, R))
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
@show VERSION
@time using DifferentialEquations
@time using Plots
@time using ProgressMeter
@time using Fmt # pkg> add https://github.com/bicycle1885/Fmt.jl

# %%
x ⪅ y = x < y || x ≈ y

function isin_ngon(n, R, x, y)
    α = 2π/n
    θ = mod(angle(x + im*y), α)
    r = √(x^2 + y^2) / R
    X, Y = r*cos(θ), r*sin(θ)
    (1 - cos(α))*Y ⪅ sin(α)*(1 - X)
end

isin_unit_pentagon(x, y) = isin_ngon(5, 1, x, y)

x = -1:0.1:1
y = -1:0.1:1
D = isin_unit_pentagon.(x, y')
heatmap(x, y, D'; colorbar=false, size=(200, 200))

# %%
module My

struct Param{TX, TY, TH, TV, TE, TB, TV_dic}
    x::TX
    y::TY
    h⁻¹::TH
    h⁻²::TH
    V::TV
    E::TE
    B::TB
    V_dic::TV_dic
end

Base.length(p::Param) = length(p.V)
Base.size(p::Param) = (length(p.x), length(p.y))

mystep(x::AbstractVector) = x[2] - x[1]
mystep(x::StepRangeLen) = step(x)

function Param(isin_func, x::AbstractVector, y::AbstractVector)
    @assert mystep(x) == mystep(y)
    h⁻¹ = 1/mystep(x)
    h⁻² = 1/mystep(x)^2
    I, J = eachindex(x), eachindex(y)
    D = isin_func.(x, y')
    V = Tuple{Int64, Int64}[]
    E_mat = [Tuple{Int64, Int64}[] for i in I, j in J]
    B_mat = [Tuple{Int64, Int64}[] for i in I, j in J]
    for j in J, i in I
        if D[i, j]
            push!(V, (i, j))
            for (p, q) in ((-1, 0), (1, 0), (0, -1), (0, 1))
                if i+p in I && j+q in J && D[i+p, j+q]
                    push!(E_mat[i, j], (i+p, j+q))
                else
                    D[i-p, j-q] && push!(B_mat[i, j], (i-p, j-q))
                end
            end
        end
    end
    V_dic = Dict(V[k] => k for k in eachindex(V))
    E = [getindex.(Ref(V_dic), E_mat[V[k]...]) for k in eachindex(V)]
    B = [getindex.(Ref(V_dic), B_mat[V[k]...]) for k in eachindex(V)]
    Param(x, y, h⁻¹, h⁻², V, E, B, V_dic)
end

function vec2mat!(p::Param, u, U)
    (; V) = p
    for k in eachindex(V)
        U[V[k]...] = u[k]
    end
    U
end

function vec2mat(p::Param, u; default_value=NaN)
    U = fill(eltype(u)(default_value), size(p))
    vec2mat!(p, u, U)
end

function mat2vec!(p::Param, U, u)
    (; V_dic) = p
    for t in keys(V_dic)
        u[V_dic[t]] = U[t...]
    end
    u
end

function mat2vec(p::Param, U)
    u = similar(eltype(U)[], length(p))
    mat2vec!(p, U, u)
end

end

x = -1:0.1:1
y = -1:0.1:1
p = My.Param(isin_unit_pentagon, x, y)

(; h⁻², V, E, B) = p
@show h⁻²

for k in (18, 19, 32, 131)
    @show k
    @show V[k]
    @show E[k]
    @show B[k]
    @show V[E[k]]
    @show V[B[k]]
    println("-"^80)
end

u0 = collect(1.0:length(V))
U0 = Int.(My.vec2mat(p, u0; default_value=0))
@show My.mat2vec(p, U0) == u0
U0'

# %%
function anim_sol2D(sol; L=201, fps=20, fn="tmp.gif",
        clim=(-0.06, 0.06), color=:bwr, bg=:lightgray)
    (; prob, u) = sol
    (; tspan, p) = prob
    (; x, y) = p
    ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
    prog = Progress(length(ts), 0)
    anim = @animate for t in ts
        U = My.vec2mat(p, sol(t).x[2])
        heatmap(x, y, U'; colorbar=false, clim, color, bg)
        title!(f"t = {$t:4.2f}"; titlefontsize=12)
        plot!(size=(420, 440), xlim=extrema(x), ylim=extrema(y), aspectratio=1)
        next!(prog)
    end
    gif(anim, fn; fps)
end

function plot_sol(sol, t; zlim=(-0.06, 0.06),
        color=:bwr, bg=:lightgray, camera=(30, 85), size=(640, 480), kwargs...)
    (; prob, u) = sol
    (; tspan, p) = prob
    (; x, y) = p
    U = My.vec2mat(p, sol(t).x[2])
    surface(x, y, U'; colorbar=false, zlim, color, bg, camera, size)
    title!(f"t = {$t:5.2f}"; titlefontsize=12)
    plot!(; kwargs...)
end

function anim_sol(sol; L=101, fps=10, fn="tmp.gif", zlim=(-0.06, 0.06),
        color=:bwr, bg=:lightgray, camera=(30, 85), size=(640, 480), kwargs...)
    (; prob, u) = sol
    (; tspan, p) = prob
    (; x, y) = p
    umin = minimum(minimum(u.x[2]) for u in u)
    umax = maximum(maximum(u.x[2]) for u in u)
    ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
    prog = Progress(length(ts), 0)
    anim = @animate for t in ts
        plot_sol(sol, t; color, bg, camera, size, zlim, kwargs...)
        next!(prog)
    end
    gif(anim, fn; fps)
end

# %%
function g!(dv, v, u, p, t)
    (; h⁻², E) = p
    @inbounds for k in eachindex(dv)
        dv[k] = h⁻² * sum(u[l] - u[k] for l in E[k])
    end
end

function g_dirichlet!(dv, v, u, p, t)
    (; h⁻², E) = p
    @inbounds for k in eachindex(dv)
        dv[k] = h⁻² * (
            sum(u[l] - u[k] for l in E[k]; init=0.0) - (4 - length(E[k])) * u[k]
        )
    end
end

function g_neumann!(dv, v, u, p, t)
    (; h⁻², E, B) = p
    @inbounds for k in eachindex(dv)
        dv[k] = h⁻² * (
            sum(u[l] - u[k] for l in E[k]; init=0.0) + sum(u[l] - u[k] for l in B[k]; init=0.0)
        )
    end
end

function g_absorbing!(dv, v, u, p, t)
    (; h⁻¹, h⁻², E, B) = p
    @inbounds for k in eachindex(dv)
        dv[k] = h⁻² * (
            sum(u[l] - u[k] for l in E[k]; init=0.0) + sum(u[l] - u[k] for l in B[k]; init=0.0)
        ) - h⁻¹ * 2length(B[k]) * v[k]
    end
end

# %%
x = y = range(-1, 1; length=21)
p = My.Param(isin_unit_pentagon, x, y)
tspan = (0.0, 2.0)

F(a, b, c, x, y) = exp(-((x-a)^2+(y-b)^2)/(2c^2))
f(x, y) = F(0.4, 0.1, 0.1, x, y)
U0 = f.(x, y')
V0 = zero(U0)

u0 = My.mat2vec(p, U0)
v0 = My.mat2vec(p, V0);

# %%
tspan = (0.0, 2.0)
prob = SecondOrderODEProblem(g_dirichlet!, v0, u0, tspan, p)

sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.01)

(; u) = sol
umin = minimum(minimum(u.x[2]) for u in u)
umax = maximum(maximum(u.x[2]) for u in u)
@show umin, umax;

# %%
anim_sol2D(sol; fn="wave eq on pentagon with dirichlet bc 2d Part 4.gif", clim=(-0.15, 0.15))

# %%
plot_sol(sol, 1.0; zlim=(-0.15, 0.15))

# %%
anim_sol(sol; fn="wave eq on pentagon with dirichlet bc Part 4.gif", zlim=(-0.15, 0.15))

# %%
tspan = (0.0, 2.0)
prob = SecondOrderODEProblem(g_neumann!, v0, u0, tspan, p)

sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.01)

(; u) = sol
umin = minimum(minimum(u.x[2]) for u in u)
umax = maximum(maximum(u.x[2]) for u in u)
@show umin, umax;

# %%
anim_sol2D(sol; fn="wave eq on pentagon with neumann 2d bc Part 4.gif", clim=(-0.15, 0.15))

# %%
plot_sol(sol, 1.0; zlim=(-0.15, 0.15))

# %%
anim_sol(sol; fn="wave eq on pentagon with neumann bc Part 4.gif", zlim=(-0.15, 0.15))

# %%
tspan = (0.0, 3.0)
prob = SecondOrderODEProblem(g_absorbing!, v0, u0, tspan, p)

sol = nothing; GC.gc()
@time sol = solve(prob; saveat=0.01)

(; u) = sol
umin = minimum(minimum(u.x[2]) for u in u)
umax = maximum(maximum(u.x[2]) for u in u)
@show umin, umax;

# %%
anim_sol2D(sol; L=301, fn="wave eq on pentagon with absorbing bc 2d Part 4.gif", clim=(-0.15, 0.15))

# %%
plot_sol(sol, 2.0; zlim=(-0.15, 0.15))

# %%
anim_sol(sol; L=151, fn="wave eq on pentagon with absorbing bc Part 4.gif", zlim=(-0.15, 0.15))

# %%

max_stars_repo_path: 0006/Absorbing boundary condition - pentagon case Part 4.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==8 ,:occupation] .= 7
df[df.occupation.==9 ,:occupation] .= 7
df[df.occupation.==10,:occupation] .= 7
df[df.occupation.==11,:occupation] .= 7
df[df.occupation.==12,:occupation] .= 7
df[df.occupation.==13,:occupation] .= 7
freqtable(df, :occupation) # problem solved

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function mlogit(alpha, X, y)

    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    num = zeros(N,J)
    dem = zeros(N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j])
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end

alpha_zero = zeros(6*size(X,2))
alpha_rand = rand(6*size(X,2))
alpha_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]
alpha_start = alpha_true.*rand(size(alpha_true))
println(size(alpha_true))
alpha_hat_optim = optimize(a -> mlogit(a, X, y), alpha_true, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
alpha_hat_mle = alpha_hat_optim.minimizer
println(alpha_hat_mle)


text:
module ActivateCmd

using Crayons.Box
using Comonicon.Tools: prompt
using ..Options
using ..InstallCmd: install_julia, create_symlink

function activate(version::String)
    ion = Options.read()
    bin = Options.find_julia_bin(version, ion)
    if bin === nothing
        if prompt("cannot find julia $version, do you want to install?")
            install_julia(version, true)
        end
    else
        @info "activating $(CYAN_FG("julia-", version))"
        create_symlink(bin, "julia")
        ion.julia.active = bin
        Options.dump(ion)
    end
    return
end

end # ActivateCmd

"""
activate a specific julia version.

# Arguements

- `version`: a julia version number or `stable`/`nightly`/`latest`.
"""
@cast function activate(version::String="stable")
    try
        ActivateCmd.activate(version)
    catch e
        if e isa InterruptException
            print("canceled by user")
            return
        else
            rethrow(e)
        end
    end
end


text:
using LRUExample

TestLRU = LRUExample.UnboundedLRU{ASCIIString, ASCIIString}()
TestBLRU = LRUExample.BoundedLRU{ASCIIString, ASCIIString}(1000)

get_str(i) = ascii(vcat(map(x->[x>>4; x&0x0F], reinterpret(Uint8, [int32(i)]))...))

isbounded{L<:LRUExample.LRU}(::Type{L}) = any(map(n->n==:maxsize, L.names))
isbounded{L<:LRUExample.LRU}(l::L) = isbounded(L)

nmax = int(logspace(2, 5, 4))

function lrutest()
    #println("LRU consistency tests")
    for lru in (TestLRU,TestBLRU)
        for n in nmax
            empty!(lru)
            #@printf("  %s, %d items\n", lru, n)
            #print("    Simple eviction: ")
            for i in 1:n
                str = get_str(i)
                lru[str] = str
                @assert lru.q[1].v == str
                if isbounded(lru) && length(lru) >= lru.maxsize
                    tailstr = get_str(i-lru.maxsize+1)
                    @assert lru.q[end].v == tailstr
                end
            end
            #println("pass")

            #print("    Lookup, random access: ")
            for i in 1:n
                str = get_str(rand(1:n))
                if haskey(lru, str) # the bounded LRUs can have cache misses
                    blah = lru[str]
                    @assert lru.q[1].v == blah
                end
            end
            #println("pass")
        end
        empty!(lru)
    end
end

lrutest()

max_stars_repo_name: alanedelman/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
# updating for generalizing the model.
# user will give args to init parameters
# use 'dictionary' -> each parameter has different range

# utilize make_dict and defaults.

function InitParams(args, seed_mode=1)

    # Parameters (match the parameter order with original code)
    # --> order doesn't matter now, it will be matched with the order of parameters

    x_val = zeros(length(args))

    # random seed
    if seed_mode == 1
        rseed_param = Dict("sigma_a" => rand()*50.,
        "sigma_s_R" => rand()*50.,
        "sigma_s_L" => rand()*50.,
        "sigma_i" => rand()*30.,
        "lambda" => randn()*.01,
        "B" => rand()*12.+5.,
        "bias" => randn(),
        "phi" => rand()*0.99+0.01,
        "tau_phi" => 0.695*rand()+0.005,
        "lapse_R" => rand()*.5,
        "lapse_L" => rand()*.5,
        "input_gain_weight" => rand())

        for i in 1:length(args)
            x_val[i] = rseed_param[args[i]]
        end
    # bing's rat avg parameter set
    elseif seed_mode == 2
        params = [-0.2767, 2.0767, 75.9600, 1.9916, 8.9474, 0.1694, 0.0964, -0.0269, 0.0613]
        defaults_param = Dict("sigma_a" => 2.0767,
        "sigma_s_R" => 75.9600,
        "sigma_s_L" => 75.9600,
        "sigma_i" => 1.9916,
        "lambda" => -0.2767,
        "B" => 8.9474,
        "bias" => -0.0269,
        "phi" => 0.1694,
        "tau_phi" => 0.0964,
        "lapse_R" => 0.0613,
        "lapse_L" => 0.0613,
        "input_gain_weight" => 0.5)

        for i in 1:length(args)
            x_val[i] = defaults_param[args[i]]
        end

    # simple fixed parameter set
    elseif seed_mode == 3
        simple_param = Dict("sigma_a" => 1.,
        "sigma_s_R" => 0.1,
        "sigma_s_L" => 0.1,
        "sigma_i" => 0.2,
        "lambda" => -0.0005,
        "B" => 6.1,
        "bias" => 0.1,
        "phi" => 0.3,
        "tau_phi" => 0.1,
        "lapse_R" => 0.1,
        "lapse_L" => 0.1,
        "input_gain_weight" => 0.5)

        for i in 1:length(args)
            x_val[i] = simple_param[args[i]]
        end
    end

    return x_val
end

function ModelFitting(args, x_init, ratdata, ntrials)
    l_b = Dict("sigma_a" => 0.,
        "sigma_s_R" => 0.,
        "sigma_s_L" => 0.,
        "sigma_i" => 0.,
        "lambda" => -5.,
        "B" => 5.,
        "bias" => -5.,
        "phi" => 0.01,
        "tau_phi" => 0.005,
        "lapse_R" => 0.,
        "lapse_L" => 0.,
        "input_gain_weight" => 0.)

    u_b = Dict("sigma_a" => 200.,
        "sigma_s_R" => 200.,
        "sigma_s_L" => 200.,
        "sigma_i" => 40.,
        "lambda" => 5.,
        "B" => 25.,
        "bias" => 5.,
        "phi" => 1.2,
        "tau_phi" => 1.5,
        "lapse_R" => 1.,
        "lapse_L" => 1.,
        "input_gain_weight" => 1.)

    l = zeros(length(args))
    u = zeros(length(args))

    for i in 1:length(args)
        l[i] = l_b[args[i]]
        u[i] = u_b[args[i]]
    end

    function LL_f(x_init::Vector)
        LLs = SharedArray(Float64, ntrials)
        # return ComputeLL(LLs, ratdata["rawdata"], ntrials, args, x_init)
        return ComputeLL(LLs, ratdata["rawdata"], ntrials
            ;make_dict(args, x_init)...)
    end

    # updated for julia v0.6 (in-place order)
    function LL_g!{T}(grads::Vector{T}, x_init::Vector{T})
        LL, LLgrad = ComputeGrad(ratdata["rawdata"], ntrials, args, x_init)

        for i=1:length(x_init)
            grads[i] = LLgrad[i]
        end
    end

    function LL_fg!(x_init::Vector, grads)
        LL, LLgrad = ComputeGrad(ratdata["rawdata"], ntrials, args, x_init)

        for i=1:length(x_init)
            grads[i] = LLgrad[i]
        end
        return LL
    end

    function my_line_search!(df, x, s, x_scratch, gr_scratch, lsr, alpha,
        mayterminate, c1::Real = 1e-4, rhohi::Real = 0.5, rholo::Real = 0.1, iterations::Integer = 1_000)
        initial_alpha = 0.5
        LineSearches.bt2!(df, x, s,x_scratch, gr_scratch, lsr, initial_alpha,
                      mayterminate, c1, rhohi, rholo, iterations)
    end

    d4 = OnceDifferentiable(LL_f,LL_g!,x_init)
                                # LL_fg!)

    tic()
    # history = optimize(d4, params, l, u, Fminbox();
    #          optimizer = GradientDescent, iterations = 500, linesearch = my_line_search!, optimizer_o = Optim.Options(g_tol = 1e-12,
    #                                                                         x_tol = 1e-32,
    #                                                                         f_tol = 1e-16,
    #                                                                         iterations = 20,
    #                                                                         store_trace = true,
    #                                                                         show_trace = true,
    #                                                                         extended_trace = true
    #                                                                         ))
    history = optimize(d4, x_init, l, u, Fminbox();
             optimizer = LBFGS, optimizer_o = Optim.Options(g_tol = 1e-12,
                                                                            x_tol = 1e-10,
                                                                            f_tol = 1e-6,                                                                        iterations = 10,
                                                                            store_trace = true,
                                                                            show_trace = true,
                                                                            extended_trace = true))

    fit_time = toc()
    println(history.minimizer)
    println(history)

    ## do a single functional evaluation at best fit parameters and save likely for each trial
    # likely_all = zeros(typeof(sigma_i),ntrials)
    x_bf = history.minimizer #.minimum
    # Likely_all_trials(likely_all, x_bf, ratdata["rawdata"], ntrials)
    # LL, LLgrad, LLhess = ComputeHess(x_bf, ratdata["rawdata"], ntrials)
    LL, LLgrad, LLhess = ComputeHess(ratdata["rawdata"], ntrials, args, x_bf)

    Gs = zeros(length(history.trace),length(x_init))
    Xs = zeros(length(history.trace),length(x_init))
    fs = zeros(length(history.trace))

    for i=1:length(history.trace)
        tt = getfield(history.trace[i],:metadata)
        fs[i] = getfield(history.trace[i],:value)
        Gs[i,:] = tt["g(x)"]
        Xs[i,:] = tt["x"]
    end

    D = Dict([("x_init",x_init),
                ("parameters",args),
                ("trials",ntrials),
                ("f",history.minimum),
                ("x_converged",history.x_converged),
                ("f_converged",history.f_converged),
                ("g_converged",history.g_converged),
                ("grad_trace",Gs),
                ("f_trace",fs),
                ("x_trace",Xs),
                ("fit_time",fit_time),
                ("x_bf",history.minimizer),
                ("myfval", history.minimum),
                ("hessian", LLhess)
                ])

    # saveto_filename = *("julia_out_",ratname,"_rseed.mat")
    # WriteFile(mpath, filename, D)
    # matwrite(saveto_filename, Dict([("x_init",params),
    #                                 ("trials",ntrials),
    #                                 ("f",history.minimum),
    #                                 ("x_converged",history.x_converged),
    #                                 ("f_converged",history.f_converged),
    #                                 ("g_converged",history.g_converged),                            ("grad_trace",Gs),
    #                                 ("f_trace",fs),
    #                                 ("x_trace",Xs),
    #                                 ("fit_time",fit_time),
    #                                 ("x_bf",history.minimizer),
    #                                 ("myfval", history.minimum),
    #                                 ("hessian", LLhess)
    #                                 ]))
    return D
end

function FitSummary(mpath, filename, D)
    WriteFile(mpath, filename, D)
end

## matwrite -> data_handle?


text:
module EMLPHE
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...heat_exchangers.EMLHEX_Engine
	import EMLtypes.length
	include("PHE/PHE_PressureDrop.jl")
	include("PHE/PHE_HeatTransfer.jl")
	include("PHE/Main_PHE.jl")
	include("PHE/Thermal_PHE.jl")
	include("PHE/PHE_Geometry.jl")
	include("PHE/PHE.jl")
end


text:
immutable LabelGeometry <: Gadfly.GeometryElement
    # One of :dynamic, :left, :right, :above, :below, :centered
    position::Symbol

    # If true, hide labels that can't be made to not-overlap during dynamic
    # label layout.
    hide_overlaps::Bool

    tag::Symbol

    function LabelGeometry(;position=:dynamic, hide_overlaps::Bool=true, tag::Symbol=empty_tag)
        new(position, hide_overlaps, tag)
    end
end

element_aesthetics(::LabelGeometry) = [:x, :y, :label]

default_statistic(::LabelGeometry) = Gadfly.Stat.identity()

const label = LabelGeometry

# True if two boxes overlap
function overlaps(a::Absolute2DBox, b::Absolute2DBox)
    a.x0[1] + a.a[1] >= b.x0[1] && a.x0[1] <= b.x0[1] + b.a[1] &&
    a.x0[2] + a.a[2] >= b.x0[2] && a.x0[2] <= b.x0[2] + b.a[2]
end

# A deferred context function for labeling points in a plot. Optimizing label
# placement depends on knowing the absolute size of the containing context.
function deferred_label_context(geom::LabelGeometry,
                                aes::Gadfly.Aesthetics,
                                theme::Gadfly.Theme,
                                drawctx::ParentDrawContext)

    # Label layout is non-trivial problem. Quite a few papers and at least one
    # Phd thesis has been written on the topic. The approach here is pretty
    # simple. A label may be placed anywhere surrounding a point or hidden.
    # Simulated annealing is used to try to minimize a penalty, which is equal
    # to the number of overlapping or out fo bounds labels, plus terms
    # penalizing hidden labels.
    #
    # TODO:
    # Penalize to prefer certain label positions over others.

    parent_transform = drawctx.t
    units = drawctx.units
    parent_box = drawctx.box

    canvas_width, canvas_height = parent_box.a[1], parent_box.a[2]

    # This should maybe go in theme? Or should we be using Aesthetics.size?
    padding = 2mm

    point_positions = Array(AbsoluteVec2, 0)
    for (x, y) in zip(aes.x, aes.y)
        x = Compose.resolve_position(parent_box, units, parent_transform, Compose.x_measure(x))
        y = Compose.resolve_position(parent_box, units, parent_transform, Compose.y_measure(y))
        x -= parent_box.x0[1]
        y -= parent_box.x0[2]
        push!(point_positions, (x, y))
    end

    extents = text_extents(theme.point_label_font,
                           theme.point_label_font_size,
                           aes.label...)
    extents = AbsoluteVec2[(width + padding, height + padding)
                           for (width, height) in extents]

    positions = Absolute2DBox[]
    for (i, (text_width, text_height)) in enumerate(extents)
        x, y = point_positions[i]
        push!(positions, Absolute2DBox((x, y), (text_width, text_height)))
    end

    # TODO: use Aesthetics.size and/or theme.default_point_size
    for (x, y) in point_positions
        push!(positions, Absolute2DBox((x - 0.5mm, y - 0.5mm), (1.0mm, 1.0mm)))
        push!(extents, (1mm, 1mm))
    end

    n = length(aes.label)

    # Return a box containing every point that the label could possibly overlap.
    function max_extents(i)
        Absolute2DBox((positions[i].x0[1] - extents[i][1],
                       positions[i].x0[2] - extents[i][2]),
                       (2*extents[i][1], 2*extents[i][2]))
    end

    # True if a is fully contained in box.
    function box_contains(a::Absolute2DBox)
        0mm < a.x0[1] && a.x0[1] + a.a[1] < parent_box.a[1] &&
        0mm < a.x0[2] - a.a[2] && a.x0[2] < parent_box.a[2]
    end

    # Checking for label overlaps is O(n^2). To mitigate these costs, we build a
    # sparse overlap matrix. This also costs O(n^2), but we only have to do it
    # once, rather than every iteration of annealing.
    possible_overlaps = [Array(Int, 0) for _ in 1:length(positions)]

    # TODO: this whole thing would be much more effecient if we forbid from
    # the start labels that overlap points. We should be able to precompute
    # that, since they're static.

    for j in 1:n
        for i in (j+1):n
            if overlaps(max_extents(i), max_extents(j))
                push!(possible_overlaps[i], j)
                push!(possible_overlaps[j], i)
            end
        end

        for i in (n+1):length(positions)
            # skip the point box corresponding to label
            if i == j + n
                continue
            end

            if overlaps(positions[i], max_extents(j))
                push!(possible_overlaps[i], j)
                push!(possible_overlaps[j], i)
            end
        end
    end

    # This variable holds the value of the objective function we wish to
    # minimize. A label overlap is a penalty of 1. Other penaties (out of bounds
    # labels, hidden labels) or calibrated to that.
    total_penalty = 0.0

    for i in 1:n
        if !box_contains(positions[i])
            total_penalty += theme.label_out_of_bounds_penalty
        end
    end

    for j in 1:n
        for i in possible_overlaps[j]
            if i > j && overlaps(positions[i], positions[j])
                total_penalty += 1
            end
        end
    end

    label_visibility = fill(true, length(positions))

    num_iterations = n * theme.label_placement_iterations
    for k in 1:num_iterations
        if total_penalty == 0
            break
        end
        j = rand(1:n)

        new_total_penalty = total_penalty
        pos = BoundingBox(0mm, 0mm, 0mm, 0mm)
        propose_hide = false

        # Propose flipping the visibility of the label.
        if label_visibility[j] && geom.hide_overlaps && rand() < theme.label_visibility_flip_pr
            new_total_penalty += theme.label_hidden_penalty

            propose_hide = true

        # Propose a change to label placement.
        else
            if !label_visibility[j]
                new_total_penalty -= theme.label_hidden_penalty
            end

            r = rand()
            point_x, point_y = point_positions[j]
            xspan = extents[j][1]
            yspan = extents[j][2]

            # TODO: it doesn't know what xspan, yspan are. That's the major
            # source of slowness.
            if rand() < 0.5
                xpos = Gadfly.lerp(rand(),
                                   (point_x - 7xspan/8),
                                   (point_x - 6xspan/8))
            else
                xpos = Gadfly.lerp(rand(),
                                   (point_x - 2xspan/8),
                                   (point_x - 1xspan/8))
            end

            ypos = Gadfly.lerp(rand(),
                               (point_y - 3yspan/4),
                               (point_y - 1yspan/4))

            # choose a side
            if r < 0.25 # top
                pos = Absolute2DBox((xpos, point_y - extents[j][2]),
                                    (extents[j][1], extents[j][2]))
            elseif 0.25 <= r < 0.5 # right
                pos = Absolute2DBox((point_x, ypos),
                                    (extents[j][1], extents[j][2]))
            elseif 0.5 <= r < 0.75 # bottom
                pos = Absolute2DBox((xpos, point_y),
                                    (extents[j][1], extents[j][2]))
            else # left
                pos = Absolute2DBox((point_x - extents[j][1], ypos),
                                    (extents[j][1], extents[j][2]))
            end
        end

        if !box_contains(positions[j]) && label_visibility[j]
            new_total_penalty -= theme.label_out_of_bounds_penalty
        end

        if !propose_hide && !box_contains(pos)
            new_total_penalty += theme.label_out_of_bounds_penalty
        end

        for i in possible_overlaps[j]
            if label_visibility[i] && label_visibility[j] &&
                overlaps(positions[i], positions[j])
                new_total_penalty -= 1
            end

            if  label_visibility[i] && overlaps(positions[i], pos)
                new_total_penalty += 1
            end
        end

        improvement = total_penalty - new_total_penalty

        T = 0.1 * (1.0 - (k / (1 + num_iterations)))
        if improvement >= 0 || rand() < exp(improvement / T)
            if propose_hide
                label_visibility[j] = false
            else
                label_visibility[j] = true
                positions[j] = pos
            end
            total_penalty = new_total_penalty
        end

    end

    return compose!(
        context(),
        text([positions[i].x0[1] + extents[i][1]/2 + parent_box.x0[1] for i in 1:n],
             [positions[i].x0[2] + extents[i][2]/2 + parent_box.x0[2] for i in 1:n],
             aes.label,
             [hcenter], [vcenter]; tag=geom.tag),
        visible(label_visibility),
        font(theme.point_label_font),
        fontsize(theme.point_label_font_size),
        fill(theme.point_label_color),
        stroke(nothing),
        svgclass("geometry"))
end

const label_layouts = @compat Dict(
    :left     => (hright,  vcenter, -2mm,  0mm),
    :right    => (hleft,   vcenter,  2mm,  0mm),
    :above    => (hcenter, vbottom,  0mm, -2mm),
    :below    => (hcenter, vtop,     0mm,  2mm),
    :centered => (hcenter, vcenter,  0mm,  0mm)
)

function render(geom::LabelGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.Label", aes, :label, :x, :y)

    if geom.position == :dynamic
        return ctxpromise(drawctx -> deferred_label_context(geom, aes, theme, drawctx))
    else
        if !in(geom.position, [:left, :right, :above, :below, :centered])
            error("""
                The position argument of Geom.label must be one of :dynamic,
                :left, :right, :above, :below, :centered
                """)
        end

        hpos, vpos, xoff, yoff = label_layouts[geom.position]

        return compose!(
            context(),
            text([Compose.x_measure(x) + xoff for x in aes.x],
                 [Compose.y_measure(y) + yoff for y in aes.y],
                 aes.label,
                 [hpos], [vpos]; tag=geom.tag),
            font(theme.point_label_font),
            fontsize(theme.point_label_font_size),
            fill(theme.point_label_color),
            stroke(nothing),
            svgclass("geometry"))
    end
end


text:
replacein(x, oldnode, newnode) = x
replacein(x::Tuple, oldnode, newnode) = tuple([replacein(m, oldnode, newnode) for m in x]...)
replacein(x::NamedTuple, oldnode, newnode) = (;[k => replacein(x[k], oldnode, newnode) for k in keys(x)]...)

for T in [ArrayNode, ArrayModel, BagModel, BagNode, TreeNode, ProductModel]
	@eval begin
		function replacein(x::$T, oldnode, newnode)
			x == oldnode && return(newnode)
			$T(map(f -> replacein(getproperty(x, f), oldnode, newnode), fieldnames($T))...)
		end
	end
end