text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================

# ================================================================
# Changes:
#    2014-12-30 added 2D radial and 3D cylindrical grids
#    2015-01-10 extended to accept nonuniform grids
# ================================================================

# %% ======================= DIFFUSION TERM ================================
"""
"""
function diffusionTerm(D::FaceValue)
    d = D.domain.dimension
    if d==1
      M = diffusionTerm1D(D)
    elseif d==1.5
      M = diffusionTermCylindrical1D(D)
    elseif d==2
      M, Mx, My = diffusionTerm2D(D)
    elseif d==2.5
      M, Mx, My = diffusionTermCylindrical2D(D)
    elseif d==2.8
      M, Mx, My = diffusionTermRadial2D(D)
    elseif d==3
      M, Mx, My, Mz = diffusionTerm3D(D)
    elseif d==3.2
      M, Mx, My, Mz = diffusionTermCylindrical3D(D)
    end
    M
end

# ======================== 1D CARTESIAN DIFFUSION =========================
function diffusionTerm1D(D::FaceValue)
# D is a face variable

# extract data from the mesh structure
Nx = D.domain.dims[1]
G = [1:Nx+2;]
DX = D.domain.cellsize.x
dx = 0.5*(DX[1:end-1]+DX[2:end])

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# reassign the east, west for code readability
De = D.xvalue[2:Nx+1]./(dx[2:Nx+1].*DX[2:Nx+1])
Dw = D.xvalue[1:Nx]./(dx[1:Nx].*DX[2:Nx+1])

# calculate the coefficients for the internal cells
AE = reshape(De,Nx)
AW = reshape(Dw,Nx)
APx = -(AE+AW)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)
end

# ======================== 2D CARTESIAN DIFFUSION =========================
function diffusionTerm2D(D::FaceValue)
# D is a face variable
# extract data from the mesh structure
Nx = D.domain.dims[1]
Ny = D.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)

DX = D.domain.cellsize.x
DY = Array(Float64, 1, Ny+2)
DY[:] = D.domain.cellsize.y
dx = 0.5*(DX[1:end-1]+DX[2:end])
dy = Array(Float64, 1, Ny+1)
dy[:] = 0.5*(DY[1:end-1]+DY[2:end])

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2))
mnx = Nx*Ny
mny = Nx*Ny

# reassign the east, west for code readability (use broadcasting in Julia)
De = D.xvalue[2:Nx+1,:]./(dx[2:Nx+1].*DX[2:Nx+1])
Dw = D.xvalue[1:Nx,:]./(dx[1:Nx].*DX[2:Nx+1])
Dn = D.yvalue[:,2:Ny+1]./(dy[:,2:Ny+1].*DY[:,2:Ny+1])
Ds = D.yvalue[:,1:Ny]./(dy[:,1:Ny].*DY[:,2:Ny+1])

# calculate the coefficients for the internal cells
AE = reshape(De,mnx)
AW = reshape(Dw,mnx)
AN = reshape(Dn,mny)
AS = reshape(Ds,mny)
APx = -(AE+AW)
APy = -(AN+AS)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1],mnx); reshape(G[2:Nx+1,2:Ny+1],mnx); reshape(G[3:Nx+2,2:Ny+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny],mny); reshape(G[2:Nx+1,2:Ny+1],mny); reshape(G[2:Nx+1,3:Ny+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
M = Mx + My
(M, Mx, My)
end

# ======================== 3D CARTESIAN DIFFUSION =========================
function diffusionTerm3D(D::FaceValue)
# D is a face variable
# extract data from the mesh structure
Nx = D.domain.dims[1]
Ny = D.domain.dims[2]
Nz = D.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)

# DX = repeat(D.domain.cellsize.x, outer=(1, Ny, Nz))
# DY = repeat(D.domain.cellsize.y.', outer=(Nx, 1, Nz))
# DZ = zeros(1,1,Nz+2);
# DZ(1,1,:) = D.domain.cellsize.z;
# DZ=repmat(DZ, Nx, Ny, 1);
# dx = 0.5*(DX(1:end-1,:,:)+DX(2:end,:,:));
# dy = 0.5*(DY(:,1:end-1,:)+DY(:,2:end,:));
# dz = 0.5*(DZ(:,:,1:end-1)+DZ(:,:,2:end));

DX = D.domain.cellsize.x
DY = Array(Float64, 1, Ny+2)
DY[:] = D.domain.cellsize.y
DZ = Array(Float64, 1,1,Nz+2)
DZ[:] = D.domain.cellsize.z
dx = 0.5*(DX[1:end-1]+DX[2:end])
dy = 0.5*(DY[:,1:end-1]+DY[:,2:end])
dz = 0.5*(DZ[:,:,1:end-1]+DZ[:,:,2:end])

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sz = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
mnx = Nx*Ny*Nz
mny = Nx*Ny*Nz
mnz = Nx*Ny*Nz

# reassign the east, west, north, and south velocity vectors for the
# code readability (use broadcasting)
De = D.xvalue[2:Nx+1,:,:]./(dx[2:Nx+1].*DX[2:Nx+1])
Dw = D.xvalue[1:Nx,:,:]./(dx[1:Nx].*DX[2:Nx+1])
Dn = D.yvalue[:,2:Ny+1,:]./(dy[:,2:Ny+1].*DY[:,2:Ny+1])
Ds = D.yvalue[:,1:Ny,:]./(dy[:,1:Ny].*DY[:,2:Ny+1])
Df = D.zvalue[:,:,2:Nz+1]./(dz[:,:,2:Nz+1].*DZ[:,:,2:Nz+1])
Db = D.zvalue[:,:,1:Nz]./(dz[:,:,1:Nz].*DZ[:,:,2:Nz+1])

# calculate the coefficients for the internal cells
AE = reshape(De,mnx)
AW = reshape(Dw,mnx)
AN = reshape(Dn,mny)
AS = reshape(Ds,mny)
AF = reshape(Df,mnz)
AB = reshape(Db,mnz)
APx = reshape(-(De+Dw),mnx)
APy = reshape(-(Dn+Ds),mny)
APz = reshape(-(Df+Db),mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3);
rowz_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1,2:Nz+1],mnx);
		reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx);
		reshape(G[3:Nx+2,2:Ny+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny,2:Nz+1],mny);
		reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny);
		reshape(G[2:Nx+1,3:Ny+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nx+1,2:Ny+1,1:Nz],mnz);
		reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz);
		reshape(G[2:Nx+1,2:Ny+1,3:Nz+2],mnz)];
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
M = Mx + My + Mz
(M, Mx, My, Mz)
end

# ======================== 1D Cylindrical DIFFUSION =========================
function diffusionTermCylindrical1D(D::FaceValue)
# D is a face variable

# extract data from the mesh structure
Nx = D.domain.dims[1]
G = [1:Nx+2;]
DX = D.domain.cellsize.x
dx = 0.5*(DX[1:end-1]+DX[2:end])
rp = D.domain.cellcenters.x
rf = D.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# reassign the east, west for code readability
De = rf[2:Nx+1].*D.xvalue[2:Nx+1]./(rp.*dx[2:Nx+1].*DX[2:Nx+1])
Dw = rf[1:Nx].*D.xvalue[1:Nx]./(rp.*dx[1:Nx].*DX[2:Nx+1])

# calculate the coefficients for the internal cells
AE = reshape(De,Nx)
AW = reshape(Dw,Nx)
APx = reshape(-(De+Dw),Nx)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)
end

# ======================== 2D RADIAL DIFFUSION =========================
function diffusionTermRadial2D(D::FaceValue)
# D is a face variable
# extract data from the mesh structure

Nr = D.domain.dims[1]
Ntheta = D.domain.dims[2]
G=reshape([1:(Nr+2)*(Ntheta+2);], Nr+2, Ntheta+2)
DR = D.domain.cellsize.x
DTHETA = Array(Float64, 1, Ntheta+2)
DTHETA[:] = D.domain.cellsize.y
dr = 0.5*(DR[1:end-1]+DR[2:end])
dtheta = Array(Float64, 1, Ntheta+1)
dtheta[:] = 0.5*(DTHETA[1:end-1]+DTHETA[2:end])
rp = D.domain.cellcenters.x
rf = D.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
mnx = Nr*Ntheta
mny = Nr*Ntheta

# reassign the east, west for code readability
De = rf[2:Nr+1].*D.xvalue[2:Nr+1,:]./(rp.*dr[2:Nr+1].*DR[2:Nr+1])
Dw = rf[1:Nr].*D.xvalue[1:Nr,:]./(rp.*dr[1:Nr].*DR[2:Nr+1])
Dn = D.yvalue[:,2:Ntheta+1]./(rp.*rp.*dtheta[:,2:Ntheta+1].*DTHETA[:,2:Ntheta+1])
Ds = D.yvalue[:,1:Ntheta]./(rp.*rp.*dtheta[:,1:Ntheta].*DTHETA[:,2:Ntheta+1])

# calculate the coefficients for the internal cells
AE = reshape(De,mnx)
AW = reshape(Dw,mnx)
AN = reshape(Dn,mny)
AS = reshape(Ds,mny)
APx = reshape(-(De+Dw),mnx)
APy = reshape(-(Dn+Ds),mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta],mny); reshape(G[2:Nr+1,2:Ntheta+1],mny); reshape(G[2:Nr+1,3:Ntheta+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
M = Mx + My
(M, Mx, My)
end

# ======================== 2D CYLINDRICAL DIFFUSION =========================
function diffusionTermCylindrical2D(D::FaceValue)
# D is a face variable
# extract data from the mesh structure
Nr = D.domain.dims[1]
Nz = D.domain.dims[2]
G=reshape([1:(Nr+2)*(Nz+2);], Nr+2, Nz+2)
DR = D.domain.cellsize.x
DZ = Array(Float64, 1, Nz+2)
DZ[:] = D.domain.cellsize.y
dr = 0.5*(DR[1:end-1]+DR[2:end])
dz = Array(Float64, 1, Nz+1)
dz[:] = 0.5*(DZ[1:end-1]+DZ[2:end])
rp = repmat(D.domain.cellcenters.x, 1, Nz)
rf = repmat(D.domain.facecenters.x, 1, Nz)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Nz+2))
mnx = Nr*Nz
mny = Nr*Nz

# extract the diffusion coefficient data
Dx = D.xvalue
Dy = D.yvalue

# reassign the east, west for code readability
De = rf[2:Nr+1,:].*D.xvalue[2:Nr+1,:]./(rp.*dr[2:Nr+1].*DR[2:Nr+1])
Dw = rf[1:Nr,:].*D.xvalue[1:Nr,:]./(rp.*dr[1:Nr].*DR[2:Nr+1])
Dn = D.yvalue[:,2:Nz+1]./(dz[:,2:Nz+1].*DZ[:,2:Nz+1])
Ds = D.yvalue[:,1:Nz]./(dz[:,1:Nz].*DZ[:,2:Nz+1])

# calculate the coefficients for the internal cells
AE = reshape(De,mnx)
AW = reshape(Dw,mnx)
AN = reshape(Dn,mny)
AS = reshape(Ds,mny)
APx = reshape(-(De+Dw),mnx)
APy = reshape(-(Dn+Ds),mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Nz+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Nz+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Nz],mny); reshape(G[2:Nr+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Nz+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
M = Mx + My
(M, Mx, My)
end

# ======================== 3D CYLINDRICAL DIFFUSION =========================
function diffusionTermCylindrical3D(D::FaceValue)
# extract data from the mesh structure

Nr = D.domain.dims[1]
Ntheta = D.domain.dims[2]
Nz = D.domain.dims[3]
G=reshape([1:(Nr+2)*(Ntheta+2)*(Nz+2);], Nr+2, Ntheta+2, Nz+2)
DR = D.domain.cellsize.x
DTHETA = Array(Float64, 1, Ntheta+2)
DTHETA[:] = D.domain.cellsize.y
DZ = Array(Float64, 1,1,Nz+2)
DZ[:] = D.domain.cellsize.z
dr = 0.5*(DR[1:end-1]+DR[2:end])
dtheta = 0.5*(DTHETA[:,1:end-1]+DTHETA[:,2:end])
dz = 0.5*(DZ[:,:,1:end-1]+DZ[:,:,2:end])
rp = D.domain.cellcenters.x # use broadcasting
rf = D.domain.facecenters.x

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sz = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
mnx = Nr*Ntheta*Nz
mny = Nr*Ntheta*Nz
mnz = Nr*Ntheta*Nz

# extract the velocity data
# note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
Dx = D.xvalue
Dy = D.yvalue
Dz = D.zvalue

# reassign the east, west, north, and south velocity vectors for the
# code readability
De = rf[2:Nr+1].*D.xvalue[2:Nr+1,:,:]./(rp.*dr[2:Nr+1].*DR[2:Nr+1])
Dw = rf[1:Nr].*D.xvalue[1:Nr,:,:]./(rp.*dr[1:Nr].*DR[2:Nr+1])
Dn = D.yvalue[:,2:Ntheta+1,:]./(rp.*rp.*dtheta[:,2:Ntheta+1].*DTHETA[:,2:Ntheta+1])
Ds = D.yvalue[:,1:Ntheta,:]./(rp.*rp.*dtheta[:,1:Ntheta].*DTHETA[:,2:Ntheta+1])
Df = D.zvalue[:,:,2:Nz+1]./(dz[:,:,2:Nz+1].*DZ[:,:,2:Nz+1])
Db = D.zvalue[:,:,1:Nz]./(dz[:,:,1:Nz].*DZ[:,:,2:Nz+1])

# calculate the coefficients for the internal cells
AE = reshape(De,mnx)
AW = reshape(Dw,mnx)
AN = reshape(Dn,mny)
AS = reshape(Ds,mny)
AF = reshape(Df,mnz)
AB = reshape(Db,mnz)
APx = reshape(-(De+Dw),mnx)
APy = reshape(-(Dn+Ds),mny)
APz = reshape(-(Df+Db),mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1,2:Nz+1],mnx);
		reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx);
		reshape(G[3:Nr+2,2:Ntheta+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta,2:Nz+1],mny);
		reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny);
		reshape(G[2:Nr+1,3:Ntheta+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nr+1,2:Ntheta+1,1:Nz],mnz);
		reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz);
		reshape(G[2:Nr+1,2:Ntheta+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
M = Mx + My + Mz
(M, Mx, My, Mz)
end


text:
"""
    FiniteDifferenceGrids

A simple 1-dimensional uniform grid
for finite difference method.
"""
module FiniteDifferenceGrids

export Grid, over_elems, over_elems_real, over_elems_ghost
export ZBoundary, Zmin, Zmax, n_hat, binary, ghost_vec, ghost_dual
export first_interior, boundary, second_interior, first_ghost, boundary_points

"""
    Grid{T, I}

A simple 1-dimensional uniform grid of
type `T` for finite difference method.
"""
struct Grid{T, I}
  """
  Number of physical elements
  """
  n_elem :: I

  """
  Number of real elements
  """
  n_elem_real :: I

  """
  Number of ghost points per side
  """
  n_ghost :: I

  """
  Element size
  """
  Δz :: T

  """
  Inverse of element size squared
  """
  Δzi :: T

  """
  Inverse of element size squared
  """
  Δzi2 :: T

  """
  z-coordinate at cell edges
  """
  ze :: Vector{T}

  """
  z-coordinate at cell centers
  """
  zc :: Vector{T}

  """
  z-coordinate at bottom of boundary node
  """
  zn_min :: T

  """
  z-coordinate at top of boundary node
  """
  zn_max :: T
end

"""
    Grid(z_min::T, z_max::T, n_elem_real::Int, n_ghost::Int = 1)

A simple grid implementation that accepts the domain interval
`z_min` and `z_max`, the number of elements `n_elem_real` and
the number of ghost points per side `n_ghost`.
"""
function Grid(z_min::T, z_max::T, n_elem_real::I, n_ghost::I = 1) where {T, I}
  n_elem = n_elem_real+2*n_ghost
  Δz = (z_max-z_min)/n_elem_real
  Δzi = inv(Δz)
  Δzi2 = Δzi^2
  zc = [(i-n_ghost)*Δz-Δz/2 for i in 1:n_elem]
  ze = [(i-n_ghost)*Δz-Δz for i in 1:n_elem+1]
  zn_min = (zc[n_ghost]+zc[1+n_ghost])/2
  zn_max = (zc[end-n_ghost]+zc[end-n_ghost+1])/2
  return Grid{T, I}(n_elem, n_elem_real, n_ghost, Δz, Δzi, Δzi2, ze, zc, zn_min, zn_max)
end

function Base.show(io::IO, grid::Grid)
  println(io, "----------------- Grid")
  println(io, "zc      = ",grid.zc)
  println(io, "ze      = ",grid.ze)
  println(io, "n_elem  = ",grid.n_elem)
  println(io, "n_ghost = ",grid.n_ghost)
  println(io, "zn_min  = ",grid.zn_min)
  println(io, "zn_max  = ",grid.zn_max)
  println(io, "-----------------")
end

abstract type ZBoundary end

"""
    Zmax
Type to dispatch function to particular location
"""
struct Zmax <: ZBoundary end

"""
    Zmin
Type to dispatch function to particular location
"""
struct Zmin <: ZBoundary end

"""
    n_hat(::ZBoundary)

The outward normal vector to the boundary
"""
n_hat(::Zmin) = -1
n_hat(::Zmax) = 1

"""
    binary(::ZBoundary)

Returns 0 for Zmin and 1 for Zmax
"""
binary(::Zmin) = 0
binary(::Zmax) = 1

"""
    ghost_vec(::ZBoundary)

A 3-element vector near the boundary with 1's on ghost cells and 0's on interior cells
"""
ghost_vec(::Zmin) = [1, 0, 0]
ghost_vec(::Zmax) = [0, 0, 1]

"""
    ghost_dual(::ZBoundary)

A 2-element vector near the boundary with 1's on ghost cells and 0's on interior cells
"""
ghost_dual(::Zmin) = [1, 0]
ghost_dual(::Zmax) = [0, 1]

"""
    over_elems(grid::Grid)

Get the range of indexes to traverse real and ghost grid elements
"""
@inline over_elems(grid::Grid) = 1:grid.n_elem

"""
    over_elems_real(grid::Grid)

Get the range of indexes to traverse only real grid elements
"""
@inline over_elems_real(grid::Grid) = 1+grid.n_ghost : grid.n_elem-grid.n_ghost

"""
    first_interior(grid::Grid, ::ZBoundary)

Get the first element index in the physical domain
"""
@inline first_interior(grid::Grid, ::Zmax) = grid.n_elem-grid.n_ghost
@inline first_interior(grid::Grid, ::Zmin) = 1+grid.n_ghost

"""
    boundary(grid::Grid, ::ZBoundary)

Get the element index for edge data at the boundary
"""
@inline boundary(grid::Grid, ::Zmax) = grid.n_elem-grid.n_ghost
@inline boundary(grid::Grid, ::Zmin) = 1+grid.n_ghost

"""
    second_interior(grid::Grid, ::ZBoundary)

Get the second element index in the physical domain
"""
@inline second_interior(grid::Grid, ::Zmax) = first_interior(grid, Zmax())-1
@inline second_interior(grid::Grid, ::Zmin) = first_interior(grid, Zmin())+1

"""
    first_ghost(grid::Grid, ::ZBoundary)

Get the first element index in the ghost domain
"""
@inline first_ghost(grid::Grid, ::Zmax) = first_interior(grid, Zmax())+1
@inline first_ghost(grid::Grid, ::Zmin) = first_interior(grid, Zmin())-1

"""
    boundary_points(grid::Grid, ::ZBoundary)

Get the first element index in the ghost domain
"""
@inline boundary_points(grid::Grid, b::ZBoundary) = first_ghost(grid, b), first_interior(grid, b), second_interior(grid, b)

"""
    over_elems_ghost(grid::Grid)

Get the range of indexes to traverse only ghost grid elements
"""
@inline over_elems_ghost(grid::Grid) = setdiff(over_elems(grid), over_elems_real(grid))

"""
    over_elems_ghost(grid::Grid, ::ZBoundary)

Get the range of indexes to traverse only ghost grid elements below ``z_{min}`` and above ``z_{max}``.
"""
@inline over_elems_ghost(grid::Grid, ::Zmin) = 1 : grid.n_ghost
@inline over_elems_ghost(grid::Grid, ::Zmax) = grid.n_elem-grid.n_ghost+1 : grid.n_elem

Base.eltype(::Grid{T}) where T = T

include("GridOperators.jl")

end

ext: jl

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using Random, Distributions, Optim, Plots
include("simulation_model.jl")

struct PxTypeData <: TypeData
    v::Array{Float64}
    γ::Array{Float64}
    z::Array{Float64}
end

function PxTypeData(n=1000)
    z = rand(Uniform(10, 20), n)
    v = z .+ rand(Normal(0, 2), n) .+5
    γ = rand(Uniform(0, 3), n)
    return PxTypeData(v, γ, z)
end

function ObservedData(types::PxTypeData, β::Array{Float64})
    β = broadcastArray(β, length(types.γ))
    x = (types.z  .- types.γ.*β[:, 2].*(types.v .- β[:, 1]))./(1 .- β[:, 2].^2 .*types.γ)
    w = max.(β[:, 1] .+ β[:, 2].*x, 0)
    y = max.(types.v .- w, 0)
    return ObservedData(x, w, y)
end

function negRevenue(β::Array{Float64, 1}, types::PxTypeData)
    data = ObservedData(types, β)
    Π = sum(data.w .* data.y)
    return -Π
end

function revenue(β::Array{Float64}, types::PxTypeData)
    data = ObservedData(types, β)
    π = data.w .* data.y
    return π
end

function robustUpdate(β, types::PxTypeData, t, tuner::ExperimentTuner)
    step = tuner.α/(2+log(t))
    ϵ = tuner.ξ.*rand([-1,1], (tuner.n, length(β)))
    βp = β' .+ ϵ
    data = ObservedData(types, βp)
    π = revenue(βp, types)
    γ̂ = linReg(π, ϵ)[2:(length(β)+1)]
    return β .+ step.*γ̂
end

function fixedXObjective(β, x::Vector, types::PxTypeData)
    w = max.(β[1] .+ β[2] .* x, 0)
    y = max.(types.v .- w, 0)
    π = y .* w
    return -sum(π)
end

function naiveUpdate(β::Array{Float64}, types::PxTypeData, t::Int, tuner::ExperimentTuner)
    data = ObservedData(types, β)
    res = Optim.optimize(beta -> fixedXObjective(beta, data.x, types), β)
    return Optim.minimizer(res)
end

function runPxExperiment(T)
    n = 5000
    β_fk = fk_solution(100000, PxTypeData, negRevenue)
    β_naive = naive_solution(100000, PxTypeData, fixedXObjective)
    tuner = ExperimentTuner(T, n, 0.2, [8e-2, 5e-4], [5.0, 0.4])
    methods = [IterativeUpdater(β_fk, tuner),
               IterativeUpdater(robustUpdate, tuner),
               IterativeUpdater(naiveUpdate, tuner),
               IterativeUpdater(β_naive, tuner)]
    runExperiment(tuner, methods, PxTypeData, revenue)
    βs = []
    fk_profits = mean(methods[1].π)
    plotmethods = [methods[1], methods[2], methods[4]]
    for m in methods
        if m in plotmethods
            push!(βs, mean(m.β[:, 2], dims=2))
        end
        println(mean(m.π))
        println(mean(m.π) - fk_profits)
        #push!(βs, m.β[:, 2])
        println(m.β[T, :])
    end
    plot( βs, xlabel = "t", ylabel="Price Discrimination",
    label=["Full Information" "Learning via Experiment"  "Naive Risk Min"], ylim = [0.0, 0.8])
    savefig("figures/price_sim.pdf")
end

function test_analytic_x()
    v = 12.0; z = 10.0; γ = 1; β = [8, 0.1]
    xAnalytic = (z  .- γ.*β[2].*(v .- β[1]))./(1 .- β[2].^2 .*γ)
    res = Optim.optimize( xstar -> utility(xstar, v, γ, z, β),
                                                        [z], LBFGS())
    xNumerical = Optim.minimizer(res)
    (xNumerical[1] .- xAnalytic) <= 1e-5
end

function utility(x, v, γ, z, β)
    x = x[1]
    w = β[1] + β[2]*x
    y = v - w
    utility = (v-w)*y - 1/2*y^2 - 1/2*(x - z)^2/γ
    return -utility
end

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
using Test
using ReinforcementLearning

@testset "ReinforcementLearning" begin
    run(E`JuliaRL_BasicDQN_CartPole`)
end


text:
using PositiveFactorizations
using Base.Test

for pivot in (Val{false}, Val{true})
    A = [1 0; 0 1]
    F = cholfact(Positive, A, pivot)
    @test_approx_eq full(F) A
    F, d = ldltfact(Positive, A, pivot)
    @test_approx_eq full(F) A
    @test d == [1,1]

    A = [1 0; 0 -1]
    F = cholfact(Positive, A, pivot)
    @test_approx_eq full(F) eye(2)
    F, d = ldltfact(Positive, A, pivot)
    @test_approx_eq full(F) eye(2)
    @test d == [1,-1]

    A = [-1 0.5; 0.5 4]
    target = pivot == Val{false} ? [1 -0.5; -0.5 4.5] : [1.125 0.5; 0.5 4]
    dtarget = pivot == Val{false} ? [-1,1] : [1,-1]
    F = cholfact(Positive, A, pivot)
    @test_approx_eq full(F) target
    F = cholfact(Positive, 10*A, pivot)
    @test_approx_eq full(F) 10*target
    F, d = ldltfact(Positive, A, pivot)
    @test_approx_eq full(F) target
    @test d == dtarget

    A = [0 1; 1 0]
    F = cholfact(Positive, A, pivot)
    @test_approx_eq full(F) eye(2)
    F, d = ldltfact(Positive, A, pivot)
    @test_approx_eq full(F) eye(2)
    @test d == [0,0]

    A = rand(201,200); A = A'*A
    F = cholfact(Positive, A, pivot)
    @test_approx_eq full(F) A
    F, d = ldltfact(Positive, A, pivot)
    @test_approx_eq full(F) A
    @test all(d .== 1)
end

A = [1 0; 0 -2]
F = eigfact(Positive, A)
@test_approx_eq full(F) abs(A)
A = [1 0; 0 0]
F = eigfact(Positive, A)
@test_approx_eq full(F) eye(2)


text:
function maxconcat(arr::Vector{<:Integer})
    b = sort(dec.(arr); lt=(x, y) -> x * y < y * x, rev=true) |> join
    return try parse(Int, b) catch parse(BigInt, b) end
end

tests = ([1, 34, 3, 98, 9, 76, 45, 4],
         [54, 546, 548, 60],
         [1, 34, 3, 98, 9, 76, 45, 4, 54, 546, 548, 60])

for arr in tests
    println("Max concatenating in $arr:\n -> ", maxconcat(arr))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# jim.jl test

using MIRTjim: jim
using LaTeXStrings
using Test: @test, @test_throws

include("isplot.jl")

@test jim(:keys) isa Vector{Symbol}
@test jim(:clim) isa Nothing
@test jim(:defs) isa AbstractDict

@test MIRTjim.jim_stack isa Vector{Any}
@test jim(:push!) isa Vector{Any}
@test jim(:pop!) == jim(:reset)
@test MIRTjim.jim_stack isa Vector{Any}

@test_throws String jim(:bad)

@isplot jim(ones(4,3), title="test2", xlabel=L"x")
@isplot jim(1:4, 5:9, zeros(4,5), title="test3", ylabel=L"y")
@isplot jim(1:4, 5:9, zeros(4,5), "test3")
@isplot jim(zeros(4,5), x=1:4, y=5:9, title="test3")
@isplot jim(rand(6,4), fft0=true)
@isplot jim(x=1:4, y=5:9, rand(4,5), title="test4")
@isplot jim(x=-9:9, y=9:-1:-9, (-9:9) * (abs.((9:-1:-9) .- 5) .< 3)', title="rev")
@isplot jim(ones(3,3)) # uniform

#3D
@isplot jim(rand(4,3,5), yflip=false, title=L"test3 x^2_i")
@isplot jim(rand(4,3,5), ncol=2)
@isplot jim(rand(4,3,5), nrow=3)
@isplot jim(-3:2, -2:1, rand(6,4,4)) # zyflip

jim(:abswarn, false)
@isplot jim(rand(ComplexF32, 4,3))
@isplot jim(rand(ComplexF32, 4,3), "complex 2d")
@isplot jim(rand(ComplexF32, 4,3,5), "complex 3d")
jim(:abswarn, true)
@isplot jim(rand(4,5), color=:hsv)
@isplot jim(jim(rand(2,3)), jim(rand(3,2)) ; layout=(2,1))
@isplot jim(:blank)


text:
module BoxQP

import JuMP

using LinearAlgebra

using MIPforQP

import MathOptInterface
const MOI = MathOptInterface

function parse_box_qp(filename::String)
    open(filename) do fp
        n = parse(Int, readline(fp))
        c = parse.(Float64, split(readline(fp)))
        @assert n == length(c)
        Q = zeros(Float64, n, n)
        for i in 1:n
            Q[:,i] = parse.(Float64, split(readline(fp)))
        end
        return Q, c
    end
end

function formulate_box_qp(Q::Matrix{T}, c::Vector{T}, optimizer_factory; epigraph::Bool=false) where {T}
    n = length(c)
    @assert n == size(Q, 1) == size(Q, 2)
    model = JuMP.Model(JuMP.with_optimizer(optimizer_factory))
    JuMP.@variable(model, 0 ≤ x[1:n] ≤ 1)
    if epigraph
        JuMP.@variable(model, y)
        JuMP.@constraint(model, y ≥ -0.5dot(x, Q * x) - dot(c, x))
        JuMP.@objective(model, Min, y)
    else
        JuMP.@objective(model, Min, -0.5dot(x, Q * x) - dot(c, x))
    end
    return model
end

end  # module BoxQP

max_stars_repo_name: joehuchette/compact-mips-for-qp

max_issues_repo_name: joehuchette/compact-mips-for-qp

max_forks_repo_name: joehuchette/compact-mips-for-qp


text:
abstract type AbstractDOS end
abstract type AbstractDOSIntegrator end

struct DOS{F} <: AbstractDOS
  ωmin::Real
  ωmax::Real
  A::F
end
(dos::DOS)(ω) = dos.A(ω)

"""
Descriptor of an integrable singularity in DOS
"""
struct DOSSingularity
  # Position of the singularity, Ω_p.
  position::Real
  # Asymptotic behavior of the DOS near the singularity, S_p(ω).
  # To be chosen such that \lim_{ω \to Ω_p} [D(ω) - S_p(ω)] = 0
  asymptotics
  # Integral of 'asymptotics' over ω\in[ω_{min};ω_{max}].
  integral
end

"""
  SingularDOS represents densities of states D(ω) according to
  the following decomposition,

  D(ω) = R(ω) + \\sum_{p=1}^P S_p(ω).

  R(ω) is a smooth function on [ω_{min};ω_{max}], and the number P of
  singular contributions S_p(ω) equals the number of integrable singularities
  Ω_p of D(ω). Each term S_p(ω) is regular everywhere on
  ω\\in[ω_{min};ω_{max}] except for its corresponding ω = Ω_p.
"""
struct SingularDOS <: AbstractDOS
  # Support limits
  ωmin::Real
  ωmax::Real
  # Regular part R(ω)
  regular
  # List of singularities
  singularities::Vector{DOSSingularity}
end

SingularDOS(ωmin, ωmax, dos) = SingularDOS(ωmin, ωmax, dos, [])

"""Evaluate singular DOS at a given frequency"""
function (dos::SingularDOS)(ω)
  x = dos.regular(ω)
  for s in dos.singularities
    x += s.asymptotics(ω)
  end
  return x
end

struct DeltaDOS <: AbstractDOS
  ϵ::Vector{Float64}
  w::Vector{Float64}
end

DeltaDOS(ϵ::Number) = DeltaDOS([ϵ])
DeltaDOS(ϵ::Vector) = DeltaDOS(ϵ, ones(length(ϵ)) / length(ϵ))

"""
Support limits of a DOS object
"""
dos_support_limits(dos::AbstractDOS) = (-Inf, Inf)
dos_support_limits(dos::DOS) = (dos.ωmin, dos.ωmax)
dos_support_limits(dos::SingularDOS) = (dos.ωmin, dos.ωmax)

struct GaussKronrodDOSIntegrator <: AbstractDOSIntegrator
  atol::Float64
  rtol::Float64
  maxevals::Int
  order::Int
end

GaussKronrodDOSIntegrator(; atol=1e-10, rtol=1e-10, maxevals=10^9, order=21) = GaussKronrodDOSIntegrator(atol, rtol, maxevals, order)
function (integrator::GaussKronrodDOSIntegrator)(f, a, b)
  integral, err = quadgk(ω -> f(ω),
                         a,
                         b,
                         atol=integrator.atol,
                         rtol=integrator.rtol,
                         maxevals=integrator.maxevals,
                         order=integrator.order)
  return integral
end

function (integrator::GaussKronrodDOSIntegrator)(f, dos::AbstractDOS)
  limits = dos_support_limits(dos)
  return integrator(ω -> f(ω) * dos(ω), limits[1], limits[2])
end

"""
  Integrator for SingularDOS

  \\int_{ω_{min}}^{ω_{max}} dω D(ω) f(ω) =
  \\int_{ω_{min}}^{ω_{max}} dω R(ω) f(ω) +
  \\sum_{p=1}^P \\int_{ω_{min}}^{ω_{max}} dω S_p(ω) [f(ω) - f(Ω_p)] +
  \\sum_{p=1}^P f(Ω_p) \\int_{ω_{min}}^{ω_{max}} dω S_p(ω).

  The integrands in the first two terms of the RHS are smooth,
  and the value of the integral in the last term must be provided in the
  `integral` field of the corresponding `DOSSingularity` structure.
"""
function (integrator::GaussKronrodDOSIntegrator)(f, dos::SingularDOS)
  limits = dos_support_limits(dos)
  val = integrator(ω -> f(ω) * dos.regular(ω), limits[1], limits[2])
  for s in dos.singularities
    f_s = f(s.position)
    val += integrator(ω -> ω ≈ s.position ? .0 : s.asymptotics(ω) * (f(ω) - f_s), limits[1], limits[2])
    val += f_s * s.integral
  end
  return val
end

function (integrator::GaussKronrodDOSIntegrator)(f, dos::DeltaDOS)
  return sum(dos.w .* f.(dos.ϵ))
end

function dos2gf(dos, β, t1::BranchPoint, t2::BranchPoint, integrator = GaussKronrodDOSIntegrator())
    theta = heaviside(t1, t2)
    Δt = t1.val - t2.val
    f = ω -> (ω > 0.0 ? exp(-1.0im * ω * (Δt - 1.0im * (1.0 - theta) * β)) / (exp(-β * ω) + 1) :
                        exp(-1.0im * ω * (Δt + 1.0im * theta * β)) / (exp(β * ω) + 1))
    return -1.0im * (2 * theta - 1) * integrator(f, dos)
end

#
# DOS factory functions
#

"""
`flat_dos(;ν=1.0, D=5.0)`

return flat band DOS with half-bandwith D and inverse cutoff width ν centered at zero
"""
flat_dos(; ν=1.0, D=5.0, μ=0.0) = DOS(-Inf, Inf, ω -> (1.0/π) * fermi(ν * (ω - μ - D)) * fermi(-ν * (ω - μ + D)))

"""
`gaussian_dos(; ϵ=1.0, ν=1.0)`

return normalized Gaussian DOS centered at ϵ with width ν
"""
gaussian_dos(; ϵ=1.0, ν=1.0) = DOS(-Inf, Inf, ω -> (1.0 / (2 * sqrt(π * ν))) * exp(-((ω - ϵ)^2)/(4ν)))

"""
`bethe_dos(; t=1.0)`

return normalized DOS of a Bethe lattice with hopping constant t
"""
bethe_dos(; t=1.0) = SingularDOS(-2t, 2t,
  ω -> begin
    if ω == -2t || ω == 2t
      -2 / (π*t)
    else
      x = ω / (2t)
      (sqrt(1 - x*x) - sqrt(2 * (1 - x)) - sqrt(2 * (1 + x))) / (π*t)
    end
  end,
  [
    DOSSingularity(-2t, ω -> sqrt(2 * (1 + ω / (2t))) / (π*t), 16 / (3*π)),
    DOSSingularity( 2t, ω -> sqrt(2 * (1 - ω / (2t))) / (π*t), 16 / (3*π))
  ]
)

"""
`chain_dos(; t=1.0)`

return normalized DOS of a linear chain with hopping constant t
"""
chain_dos(; t=1.0) = SingularDOS(-2t, 2t,
  ω -> begin
    if ω == -2t || ω == 2*t
      -3 / (8*π*t)
    else
      x = ω / (2t)
      rp = sqrt(2 * (1 - x))
      # The second term regularizes the derivative of near ω = 2t to ease integration
      sp = 1 / (2*π*t * rp) + rp / (16*π*t)

      rm = sqrt(2 * (1 + x))
      # The second term regularizes the derivative near ω = -2t to ease integration
      sm = 1 / (2*π*t * rm) + rm / (16*π*t)

      (1 / (2*π*t)) / sqrt(1 - x*x) - sp - sm
    end
  end,
  [
    DOSSingularity(-2t, ω -> begin
                          x = ω / (2t)
                          s = sqrt(2 * (1 + x))
                          # The second term regularizes the derivative
                          1 / (2*π*t*s) + s / (16*π*t)
                        end,
                    7 / (3*π)),
    DOSSingularity( 2t, ω -> begin
                          x = ω / (2t)
                          s = sqrt(2 * (1 - x))
                          # The second term regularizes the derivative
                          1 / (2*π*t*s) + s / (16*π*t)
                        end,
                    7 / (3*π))
  ]
)

"""
`square_dos(; t=1.0)`

return normalized DOS of a 2D square lattice with hopping constant t
"""
square_dos(; t=1.0) = SingularDOS(-4t, 4t,
  ω -> begin
    if ω ≈ 0
      0
    else
      x = ω / (4t)
      (1 / (2 * π^2 * t)) * (Elliptic.K(1 - x^2) + log(abs(x) / 4));
    end
  end,
  [
    DOSSingularity(0, ω -> -1 / (2 * π^2 * t) * log(abs(ω) / (16t)),
                   4 / (π^2) * (1 + 2 * log(2)))
  ]
)


text:
abstract type AbstractReadOnlyArray{T,N,P} <: AbstractArray{T,N} end

Base.parent(x::AbstractReadOnlyArray) = getfield(x, :parent)

isfrozenvalue(::Type{<:AbstractReadOnlyArray}) = true
isfrozenindex(::Type{<:AbstractReadOnlyArray}) = true

meltvalue!(x::AbstractReadOnlyArray) = parent(x)
meltindex!(x::AbstractReadOnlyArray) = parent(x)

Base.IteratorSize(::Type{<:AbstractReadOnlyArray{T,N,P}}) where {T,N,P} =
    Base.IteratorSize(P)
Base.IteratorEltype(::Type{<:AbstractReadOnlyArray{T,N,P}}) where {T,N,P} =
    Base.IteratorEltype(P)
Base.eltype(::Type{<:AbstractReadOnlyArray{T,N,P}}) where {T,N,P} = eltype(P)
Base.size(x::AbstractReadOnlyArray, args...) = size(parent(x), args...)
Base.@propagate_inbounds Base.getindex(x::AbstractReadOnlyArray, I...) =
    getindex(parent(x), I...)
Base.firstindex(x::AbstractReadOnlyArray) = firstindex(parent(x))
Base.lastindex(x::AbstractReadOnlyArray) = lastindex(parent(x))
Base.IndexStyle(::Type{<:AbstractReadOnlyArray{T,N,P}}) where {T,N,P} = IndexStyle(P)
Base.iterate(x::AbstractReadOnlyArray, args...) = iterate(parent(x), args...)
Base.length(x::AbstractReadOnlyArray) = length(parent(x))

Base.axes(x::AbstractReadOnlyArray) = axes(parent(x))
Base.strides(x::AbstractReadOnlyArray) = strides(parent(x))
Base.unsafe_convert(p::Type{Ptr{T}}, x::AbstractReadOnlyArray) where {T} =
    Base.unsafe_convert(p, parent(x))
Base.stride(x::AbstractReadOnlyArray, i::Int) = stride(parent(x), i)

Base.copy(x::AbstractReadOnlyArray) = copy(parent(x))
Base.copymutable(x::AbstractReadOnlyArray) = Base.copymutable(parent(x))

Base.push!(x::AbstractReadOnlyArray, _...) = throw(DisallowedOperation(push!, x))
Base.append!(x::AbstractReadOnlyArray, _...) = throw(DisallowedOperation(append!, x))

# For StructArrays:
Base.getproperty(x::AbstractReadOnlyArray, name::Symbol) =
    constructorof(typeof(x))(getproperty(parent(x), name))
Base.getproperty(x::AbstractReadOnlyArray, field) =
    constructorof(typeof(x))(getproperty(parent(x), field))

Base.broadcastable(x::AbstractReadOnlyArray) = Base.broadcastable(parent(x))

struct ReadOnlyArray{T,N,P<:AbstractArray{T,N}} <: AbstractReadOnlyArray{T,N,P}
    parent::P
end

isimmutablevalue(::Type{<:ReadOnlyArray}) = true
isimmutableindex(::Type{<:ReadOnlyArray}) = true
isfrozenvalue(::Type{<:ReadOnlyArray}) = false
isfrozenindex(::Type{<:ReadOnlyArray}) = false
readonly_impl(x::AbstractArray) = ReadOnlyArray(x)

struct ImmutableArray{T,N,P<:AbstractArray{T,N}} <: AbstractReadOnlyArray{T,N,P}
    parent::P
end

freeze!(x::AbstractArray) = ImmutableArray(x)

struct ImmutableIndexArray{T,N,P<:AbstractArray{T,N}} <: AbstractReadOnlyArray{T,N,P}
    parent::P
end

isfrozenvalue(::Type{<:ImmutableIndexArray}) = false
Base.setindex!(x::ImmutableIndexArray, args...) = setindex!(parent(x), args...)

freezeindex!(x::AbstractArray) = ImmutableIndexArray(x)

struct AppendOnlyVector{T,P<:AbstractVector{T}} <: AbstractReadOnlyArray{T,1,P}
    parent::P
end

isfrozenindex(::Type{<:AppendOnlyVector}) = false
Base.push!(x::AppendOnlyVector, args...) = (push!(parent(x), args...); x)
Base.append!(x::AppendOnlyVector, args...) = (append!(parent(x), args...); x)

freezevalue!(x::AbstractVector) = AppendOnlyVector(x)
freezevalue!(x::AbstractArray) = ImmutableArray(x)

factoryof(::ReadOnlyArray) = readonly
factoryof(::ImmutableArray) = freeze
factoryof(::ImmutableIndexArray) = freezeindex
factoryof(::AppendOnlyVector) = freezevalue

function Base.showarg(io::IO, x::AbstractReadOnlyArray, toplevel::Bool)
    print(io, nameof(factoryof(x)), '(')
    Base.showarg(io, parent(x), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(x))
end


text:
using GAIO
using StaticArrays
using Test

@testset "exported functionality" begin
    @testset "regular partition" begin
        partition = BoxPartition(Box((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)), (16,8,8))
        @testset "empty" begin
            empty_boxset = boxset_empty(partition)
            @test isempty(empty_boxset)
            @test length(empty_boxset) == 0

            empty_boxset = subdivide(empty_boxset, 1)
            @test isempty(empty_boxset)
            @test length(empty_boxset) == 0
        end
        @testset "full" begin
            full_boxset = partition[:]
            @test !isempty(full_boxset)
            @test length(full_boxset) == 2^10

            full_boxset = subdivide(full_boxset, 1)
            @test !isempty(full_boxset)
            @test length(full_boxset) == 2^11
        end
        p1 = SVector(0.5, 0.5, 0.5)
        p2 = SVector(-0.5, 0.5, 0.5)
        p3 = SVector(0.5, -0.5, 0.5)
        p4 = SVector(0.51, 0.51, 0.51)
        @testset "boxsets created on points" begin
            box_set = partition[(p1, p2, p3, p4)]
            @test length(box_set) == 3
            @test any(box -> p1 ∈ box, box_set)
            @test any(box -> p2 ∈ box, box_set)
            @test any(box -> p3 ∈ box, box_set)
            @test any(box -> p4 ∈ box, box_set)
            n = 10
            for k in 1:n
                box_set = subdivide(box_set, (k%3)+1 )
            end
            @test length(box_set) == 3 * 2^n
            @test any(box -> p1 ∈ box, box_set)
            @test any(box -> p2 ∈ box, box_set)
            @test any(box -> p3 ∈ box, box_set)
            @test any(box -> p4 ∈ box, box_set)
        end
        @testset "set operations" begin
            p1_box_set = partition[p1]
            p1p2_box_set = partition[(p1, p2)]
            p3_box_set = partition[p3]
            p2p3_box_set = partition[(p2, p3)]
            @test length(p1_box_set) == 1
            @test length(p1p2_box_set) == 2
            @test length(p3_box_set) == 1
            @test length(p2p3_box_set) == 2
            @test length(union(p1_box_set, p3_box_set)) == 2
            @test length(union(p1_box_set, p1p2_box_set)) == 2
            @test length(union(p1p2_box_set, p2p3_box_set)) == 3
            @test length(intersect(p1_box_set, p3_box_set)) == 0
            @test length(intersect(p1_box_set, p1p2_box_set)) == 1
            @test length(intersect(p2p3_box_set, p2p3_box_set)) == 2
            @test length(setdiff(p1_box_set, p3_box_set)) == 1
            @test length(setdiff(p1_box_set, p1p2_box_set)) == 0
            @test length(setdiff(p1p2_box_set, p1_box_set)) == 1
            union!(p1_box_set, p2p3_box_set)
            setdiff!(p1p2_box_set, p2p3_box_set)
            intersect!(p2p3_box_set, p3_box_set)
            @test length(p1_box_set) == 3
            @test length(p2p3_box_set) == 1
            @test length(p1p2_box_set) == 1
        end
    end
end

#=     @testset "tree partition" begin
        partition = TreePartition(Box(SVector(0.0, 0.0, 0.0), SVector(1.0, 1.0, 1.0)))
        @testset "empty" begin
            empty_boxset = boxset_empty(partition)
            @test isempty(empty_boxset)
            @test length(empty_boxset) == 0

            empty_boxset = subdivide(empty_boxset)
            @test isempty(empty_boxset)
            @test length(empty_boxset) == 0
        end
        full_boxset = partition[:]
        @testset "full" begin
            @test !isempty(full_boxset)
            @test length(full_boxset) == 1
            n = 10
            for k in 1:10
                full_boxset = subdivide(full_boxset)
            end
            @test !isempty(full_boxset)
            @test length(full_boxset) == 2^n
        end
        # this seems to be the easiest way right now
        # to get a full tree at a certain depth. Not
        # really very clean
        partition = full_boxset.partition
        p1 = SVector(0.5, 0.5, 0.5)
        p2 = SVector(-0.5, 0.5, 0.5)
        p3 = SVector(0.5, -0.5, 0.5)
        p4 = SVector(0.51, 0.51, 0.51)
        @testset "box sets created on points" begin
            box_set = partition[(p1, p2, p3, p4)]
            @test length(box_set) == 3
            @test any(box -> p1 ∈ box, box_set)
            @test any(box -> p2 ∈ box, box_set)
            @test any(box -> p3 ∈ box, box_set)
            @test any(box -> p4 ∈ box, box_set)
            n = 10
            for k in 1:n
                box_set = subdivide(box_set, (k%3)+1 )
            end
            @test length(box_set) == 3 * 2^n
            @test any(box -> p1 ∈ box, box_set)
            @test any(box -> p2 ∈ box, box_set)
            @test any(box -> p3 ∈ box, box_set)
            @test any(box -> p4 ∈ box, box_set)
        end
        @testset "set operations" begin
            p1_box_set = partition[p1]
            p1p2_box_set = partition[(p1, p2)]
            p3_box_set = partition[p3]
            p2p3_box_set = partition[(p2, p3)]
            @test length(p1_box_set) == 1
            @test length(p1p2_box_set) == 2
            @test length(p3_box_set) == 1
            @test length(p2p3_box_set) == 2
            @test length(union(p1_box_set, p3_box_set)) == 2
            @test length(union(p1_box_set, p1p2_box_set)) == 2
            @test length(union(p1p2_box_set, p2p3_box_set)) == 3
            @test length(intersect(p1_box_set, p3_box_set)) == 0
            @test length(intersect(p1_box_set, p1p2_box_set)) == 1
            @test length(intersect(p2p3_box_set, p2p3_box_set)) == 2
            @test length(setdiff(p1_box_set, p3_box_set)) == 1
            @test length(setdiff(p1_box_set, p1p2_box_set)) == 0
            @test length(setdiff(p1p2_box_set, p1_box_set)) == 1
            union!(p1_box_set, p2p3_box_set)
            setdiff!(p1p2_box_set, p2p3_box_set)
            intersect!(p2p3_box_set, p3_box_set)
            @test length(p1_box_set) == 3
            @test length(p2p3_box_set) == 1
            @test length(p1p2_box_set) == 1
        end
    end
end
 =#


text:
"""
    viewmsa(msa)

Create and return a Makie Figure for a AbstractMultipleSequenceAlignment.
# Examples
```julia
using MIToS.MSA
downloadpfam("pf00062")
vm = MIToS.MSA.read("pf00062.stockholm.gz",Stockholm) |> Node
fig1 = viewmsa(vm)

using FastaIO
vm = FastaIO.readfasta("data/fasta1.fas") |> Node
fig1 = viewmsa(vm)
```
Parameters:
sheetsize ----- Dimensions of the msa shown, Default - [40,20]
resolution ---- Default - (1500, 600)
colorscheme --- Default - :viridis
positions ----- Residue positions, Default - 1:(length of msa)
"""
function viewmsa(   msa::T;
					sheetsize = [40,20],
					resolution = (1500, 600),
                    colorscheme = :viridis,
                    colorval = 2,
                    positions = 0
				) where {T<:Node}
	width1 = sheetsize[1]
	height1 = sheetsize[2]
	fig = Figure(resolution = resolution)
	ax = Axis(fig[1:7,3:9])
    labels = nothing
    nums = nothing
    strmsa = nothing
	tightlimits!(ax)
    # f2 = [msa[i][2] for i in 1:size(msa,1)]
    # mat = [[f2[i]...] for i in 1:size(f2,1)] |> combinedims .|> string |> _t
    # if size(mat,1) < sheetsize[1]
    #     height1 = size(mat,1) - 1
    # else
    #     height1 = sheetsize[1]
    # end
    # if size(mat,2) < sheetsize[2]
    #     width1 = size(mat,2) - 1
    # else
    #     width1 = sheetsize[2]
    # end
    if typeof(msa[])<:MSA.AbstractMultipleSequenceAlignment
        labels = @lift $(msa).matrix.dicts[1] |> keys |> collect
	    nums = @lift $(msa).matrix.dicts[2] |> keys |> collect
        strmsa = @lift Matrix($msa) .|> string
    elseif typeof(msa[])<:Vector{Tuple{String,String}}
        if positions == 0
            positions = @lift [1:length($msa[1][2])...]
        end
        labels = @lift [$msa[i][1] for i in 1:size($msa,1)]
	    nums = typeof(positions)<:Node ? (@lift collect($positions)) : Node(positions)
        strmsa = @lift [[$fas1[i][2]...] for i in 1:size($fas1,1)] |> combinedims .|> string
    else
        error("sorry methods for that input don't exist")
    end
	labelssize = @lift size($labels,1) - (height1-1)
	labelsrange = @lift $labelssize:-1:1
	numssize = @lift size($nums,1) - (width1-1)
	numsrange = @lift 1:1:$numssize
    sl1 = GLMakie.Slider(fig[end+1,3:9], range = numsrange, startvalue = 1)
	sl1.value = 1
	sl2 = GLMakie.Slider(fig[1:7,10], range = labelsrange, startvalue = 1, horizontal = false,
		tellwidth = true, height = nothing)
	sl2.value = labelssize[]
    colorval = typeof(colorval)<:Node ? colorval : Node(colorval)
	strmsavals = @lift [ _kdict(i) for i in $strmsa ]
	strmsavals2 = @lift $strmsavals |> combinedims
	labelshow = lift(X->labels[][(X+(height1-1):-1:X)],sl2.value)
	numsshow = lift(X->nums[][(X:1:X+(width1-1))],sl1.value)
	labelshow2 = lift(X->(X+(height1-1):-1:X),sl2.value)
	numsshow2 = lift(X->(X:1:X+(width1-1)),sl1.value)
	fixtmsa = @lift replace($strmsavals2[$colorval,:,:], nothing => 0.0)
	msashow = @lift $fixtmsa[$labelshow2,$numsshow2] |> _t
	charmsa = @lift [x[1] for x in $strmsa]
	charshow = @lift $charmsa[$labelshow2,$numsshow2]
	widthrange = [1:width1...]
	heightrange = [1:height1...]
	ax.yticks = (heightrange, labelshow[])
	on(labelshow) do ls
		ax.yticks = (heightrange, ls)
	end
	ax.xticks = (widthrange, numsshow[])
	on(numsshow) do ns
		ax.xticks = (widthrange, ns)
	end
	ax.xticklabelsize = 9
	ax.yticklabelsize = 13
	ax.xzoomlock[] = true
	ax.yzoomlock[] = true
	ax.yticklabelspace[] = 10
	points1 = [Point2f0(x,y) for x in widthrange for y in heightrange] |> collect
	charvec = @lift SplitApplyCombine.flatten($charshow)
	scatter!(ax,
	        points1,
	        marker = charvec,
	        markersize = (10.0,11.0),
			color = :black,
			strokecolor = :black
	        )
	heatmap!(ax, msashow, show_grid = true, show_axis = true,
	       colormap = colorscheme
           )
    ax.attributes.xaxisposition[] = :top
    deregister_interaction!(fig.current_axis.x,:rectanglezoom)
	return fig
end


text:
using MMWavenet
using BenchmarkTools

println("Learn MMWavenet on random dateset.")

network, ps = main(2, 3, "random_data.bson");
trn_dataset, val_dataset = load_dataset("random_data.bson")
@benchmark train!(ps, network, trn_dataset, val_dataset, iter=1)

max_stars_repo_name: kozvojtex/MMWavenet

max_issues_repo_name: kozvojtex/MMWavenet

max_forks_repo_name: kozvojtex/MMWavenet


text:
# Autogenerated wrapper script for HMMER_jll for x86_64-unknown-freebsd
export alimask, esl_afetch, esl_alimanip, esl_alimap, esl_alimask, esl_alimerge, esl_alipid, esl_alistat, esl_compalign, esl_compstruct, esl_construct, esl_histplot, esl_mask, esl_reformat, esl_selectn, esl_seqrange, esl_seqstat, esl_sfetch, esl_shuffle, esl_ssdraw, esl_weight, hmmalign, hmmbuild, hmmconvert, hmmemit, hmmfetch, hmmlogo, hmmpgmd, hmmpress, hmmscan, hmmsearch, hmmsim, hmmstat, jackhmmer, makehmmerdb, nhmmer, nhmmscan, phmmer

using GSL_jll
JLLWrappers.@generate_wrapper_header("HMMER")
JLLWrappers.@declare_executable_product(alimask)
JLLWrappers.@declare_executable_product(esl_afetch)
JLLWrappers.@declare_executable_product(esl_alimanip)
JLLWrappers.@declare_executable_product(esl_alimap)
JLLWrappers.@declare_executable_product(esl_alimask)
JLLWrappers.@declare_executable_product(esl_alimerge)
JLLWrappers.@declare_executable_product(esl_alipid)
JLLWrappers.@declare_executable_product(esl_alistat)
JLLWrappers.@declare_executable_product(esl_compalign)
JLLWrappers.@declare_executable_product(esl_compstruct)
JLLWrappers.@declare_executable_product(esl_construct)
JLLWrappers.@declare_executable_product(esl_histplot)
JLLWrappers.@declare_executable_product(esl_mask)
JLLWrappers.@declare_executable_product(esl_reformat)
JLLWrappers.@declare_executable_product(esl_selectn)
JLLWrappers.@declare_executable_product(esl_seqrange)
JLLWrappers.@declare_executable_product(esl_seqstat)
JLLWrappers.@declare_executable_product(esl_sfetch)
JLLWrappers.@declare_executable_product(esl_shuffle)
JLLWrappers.@declare_executable_product(esl_ssdraw)
JLLWrappers.@declare_executable_product(esl_weight)
JLLWrappers.@declare_executable_product(hmmalign)
JLLWrappers.@declare_executable_product(hmmbuild)
JLLWrappers.@declare_executable_product(hmmconvert)
JLLWrappers.@declare_executable_product(hmmemit)
JLLWrappers.@declare_executable_product(hmmfetch)
JLLWrappers.@declare_executable_product(hmmlogo)
JLLWrappers.@declare_executable_product(hmmpgmd)
JLLWrappers.@declare_executable_product(hmmpress)
JLLWrappers.@declare_executable_product(hmmscan)
JLLWrappers.@declare_executable_product(hmmsearch)
JLLWrappers.@declare_executable_product(hmmsim)
JLLWrappers.@declare_executable_product(hmmstat)
JLLWrappers.@declare_executable_product(jackhmmer)
JLLWrappers.@declare_executable_product(makehmmerdb)
JLLWrappers.@declare_executable_product(nhmmer)
JLLWrappers.@declare_executable_product(nhmmscan)
JLLWrappers.@declare_executable_product(phmmer)
function __init__()
    JLLWrappers.@generate_init_header(GSL_jll)
    JLLWrappers.@init_executable_product(
        alimask,
        "bin/alimask",
    )

    JLLWrappers.@init_executable_product(
        esl_afetch,
        "bin/esl-afetch",
    )

    JLLWrappers.@init_executable_product(
        esl_alimanip,
        "bin/esl-alimanip",
    )

    JLLWrappers.@init_executable_product(
        esl_alimap,
        "bin/esl-alimap",
    )

    JLLWrappers.@init_executable_product(
        esl_alimask,
        "bin/esl-alimask",
    )

    JLLWrappers.@init_executable_product(
        esl_alimerge,
        "bin/esl-alimerge",
    )

    JLLWrappers.@init_executable_product(
        esl_alipid,
        "bin/esl-alipid",
    )

    JLLWrappers.@init_executable_product(
        esl_alistat,
        "bin/esl-alistat",
    )

    JLLWrappers.@init_executable_product(
        esl_compalign,
        "bin/esl-compalign",
    )

    JLLWrappers.@init_executable_product(
        esl_compstruct,
        "bin/esl-compstruct",
    )

    JLLWrappers.@init_executable_product(
        esl_construct,
        "bin/esl-construct",
    )

    JLLWrappers.@init_executable_product(
        esl_histplot,
        "bin/esl-histplot",
    )

    JLLWrappers.@init_executable_product(
        esl_mask,
        "bin/esl-mask",
    )

    JLLWrappers.@init_executable_product(
        esl_reformat,
        "bin/esl-reformat",
    )

    JLLWrappers.@init_executable_product(
        esl_selectn,
        "bin/esl-selectn",
    )

    JLLWrappers.@init_executable_product(
        esl_seqrange,
        "bin/esl-seqrange",
    )

    JLLWrappers.@init_executable_product(
        esl_seqstat,
        "bin/esl-seqstat",
    )

    JLLWrappers.@init_executable_product(
        esl_sfetch,
        "bin/esl-sfetch",
    )

    JLLWrappers.@init_executable_product(
        esl_shuffle,
        "bin/esl-shuffle",
    )

    JLLWrappers.@init_executable_product(
        esl_ssdraw,
        "bin/esl-ssdraw",
    )

    JLLWrappers.@init_executable_product(
        esl_weight,
        "bin/esl-weight",
    )

    JLLWrappers.@init_executable_product(
        hmmalign,
        "bin/hmmalign",
    )

    JLLWrappers.@init_executable_product(
        hmmbuild,
        "bin/hmmbuild",
    )

    JLLWrappers.@init_executable_product(
        hmmconvert,
        "bin/hmmconvert",
    )

    JLLWrappers.@init_executable_product(
        hmmemit,
        "bin/hmmemit",
    )

    JLLWrappers.@init_executable_product(
        hmmfetch,
        "bin/hmmfetch",
    )

    JLLWrappers.@init_executable_product(
        hmmlogo,
        "bin/hmmlogo",
    )

    JLLWrappers.@init_executable_product(
        hmmpgmd,
        "bin/hmmpgmd",
    )

    JLLWrappers.@init_executable_product(
        hmmpress,
        "bin/hmmpress",
    )

    JLLWrappers.@init_executable_product(
        hmmscan,
        "bin/hmmscan",
    )

    JLLWrappers.@init_executable_product(
        hmmsearch,
        "bin/hmmsearch",
    )

    JLLWrappers.@init_executable_product(
        hmmsim,
        "bin/hmmsim",
    )

    JLLWrappers.@init_executable_product(
        hmmstat,
        "bin/hmmstat",
    )

    JLLWrappers.@init_executable_product(
        jackhmmer,
        "bin/jackhmmer",
    )

    JLLWrappers.@init_executable_product(
        makehmmerdb,
        "bin/makehmmerdb",
    )

    JLLWrappers.@init_executable_product(
        nhmmer,
        "bin/nhmmer",
    )

    JLLWrappers.@init_executable_product(
        nhmmscan,
        "bin/nhmmscan",
    )

    JLLWrappers.@init_executable_product(
        phmmer,
        "bin/phmmer",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
config_file = isempty(ARGS) ? joinpath(@__DIR__, "config", "config.jl") : ARGS[1]
isfile(config_file) || error("Cannot read '$config_file'")
println("Config supplied: '$config_file'")
include(config_file)

using SVDD, OneClassActiveLearning
using OneClassSampling
using Memento, Gurobi, JSON, JuMP
using Random
using Distributed
using MLKernels
using Serialization
using Dates

@info "Running experiment in \"$(realpath(data_output_root))\"."
include(joinpath(@__DIR__, "util", "setup_workers.jl"))

@info "Loading packages on all workers."
@everywhere using Pkg
@everywhere using SVDD, OneClassActiveLearning, OneClassSampling, JSON, JuMP, Gurobi, Memento, Random, MLKernels, Serialization
@everywhere import SVDD: SVDDneg
@everywhere include(joinpath(@__DIR__, "util", "evaluate.jl"))
@info "Loaded."

@everywhere fmt_string = "[{name} | {date} | {level}]: {msg}"
@everywhere loglevel = "debug"

@everywhere function setup_logging(experiment)
    setlevel!(getlogger("root"), "error")
    setlevel!(getlogger(OneClassActiveLearning), loglevel)
    setlevel!(getlogger(SVDD), loglevel)

    exp_logfile = joinpath(experiment[:log_dir], "experiment", "$(experiment[:hash]).log")
    worker_logfile = joinpath(experiment[:log_dir], "worker", "$(gethostname())_$(getpid()).log")

    WORKER_LOGGER = Memento.config!("runner", "debug"; fmt=fmt_string)

    exp_handler = DefaultHandler(exp_logfile, DefaultFormatter(fmt_string))
    push!(getlogger(OneClassActiveLearning), exp_handler, experiment[:hash])
    push!(getlogger(SVDD), exp_handler, experiment[:hash])
    push!(WORKER_LOGGER, exp_handler, experiment[:hash])

    worker_handler = DefaultHandler(worker_logfile, DefaultFormatter(fmt_string))
    setlevel!(gethandlers(WORKER_LOGGER)["console"], "error")
    push!(WORKER_LOGGER, worker_handler)

    return WORKER_LOGGER
end

@everywhere function cleanup_logging(worker_logger::Logger, experiment_hash)
    delete!(gethandlers(getlogger(OneClassActiveLearning)), experiment_hash)
    delete!(gethandlers(getlogger(SVDD)), experiment_hash)
    delete!(gethandlers(worker_logger), experiment_hash)
    return nothing
end

@everywhere function Memento.warn(logger::Logger, error::Gurobi.GurobiError)
    Memento.warn(logger, "GurobiError(exit_code=$(error.code), msg='$(error.msg)')")
end

@everywhere function Memento.debug(logger::Logger, error::Gurobi.GurobiError)
    Memento.warn(logger, "GurobiError(exit_code=$(error.code), msg='$(error.msg)')")
end

@everywhere function Memento.debug(logger::Logger, error::ErrorException)
    Memento.warn(logger, "Caught ErrorException, msg='$(error.msg)')")
end

@everywhere function run_experiment(experiment::Dict, warmup=false)
    Random.seed!(experiment[:seed])

    WORKER_LOGGER = setup_logging(experiment)
    info(WORKER_LOGGER, "host: $(gethostname()) worker: $(getpid()) exp_hash: $(experiment[:hash])")
    if isfile(experiment[:output_file])
        warn(WORKER_LOGGER, "Aborting experiment because the output file already exists. Filename: $(experiment[:output_file])")
        cleanup_logging(WORKER_LOGGER, experiment[:hash])
        return nothing
    end
    e = deepcopy(experiment)

    data, labels = load_data(e[:data_file])
    e[:data_stats] = Dict(:n_observations => size(data, 2),
                            :n_dims => size(data, 1),
                            :outlier_ratio => sum(labels .== :outlier) / size(data, 2))

    errorfile = joinpath(e[:log_dir], "worker", "$(gethostname())_$(getpid())")
    try
        time_sampling = @elapsed sample_mask = OneClassSampling.sample(e[:sampling_strategy], data, labels)
        e[:result] = Dict{Symbol, Any}(:sample_size => sum(sample_mask), :time_sampling => time_sampling)
        if e[:result][:sample_size] > SAMPLE_SIZE_LIMIT
            throw(TooLargeSampleException(e[:result][:sample_size]))
        end
        @show e[:data_file]
        @show e[:sampling_strategy]
        @show e[:result][:sample_size]
        eval_dict = evaluate(e[:sampling_strategy], e[:models],
                             sample_mask,
                             data, labels, e[:quality_metrics])
        e[:result] = merge(e[:result], eval_dict)
        e[:exit_code] = :success
    catch ex
        e[:exit_code] = Symbol(typeof(ex))
        if isa(e, SamplingException)
            warn(WORKER_LOGGER, "Experiment $(e[:hash]) failed during sampling.")
            warn(WORKER_LOGGER, ex)
        elseif !isa(e, EmptySampleException)
            @warn "Experiment $(e[:hash]) finished with unkown error."
            @warn ex
            @warn stacktrace(catch_backtrace())
            warn(WORKER_LOGGER, "Experiment $(e[:hash]) finished with unkown error.")
            warn(WORKER_LOGGER, ex)
        end
        warn(WORKER_LOGGER, "Experiment $(e[:hash]) finished with unkown error.")
        warn(WORKER_LOGGER, ex)
    finally
        # Clean up for result file writing
        delete!(e, :quality_metrics)
        e[:sampling_strategy] = string(e[:sampling_strategy])
        e[:models][:sample_qe][:threshold_strategies] = string(e[:models][:sample_qe][:threshold_strategies])
        if e[:exit_code] != :success
            info(WORKER_LOGGER, "Writing error hash to $errorfile.error.")
            open("$errorfile.error", "a") do f
                print(f, "$(e[:hash])\n")
            end
        end
        if !WARMUP
            info(WORKER_LOGGER, "Writing result to $(e[:output_file]).")
            open(e[:output_file], "w") do f
                JSON.print(f, e)
            end
        end
        cleanup_logging(WORKER_LOGGER, e[:hash])
    end

    return nothing
end

all_experiments = []
for s in readdir(data_output_root)
    exp_dir = joinpath(data_output_root, s)
    if !isdir(exp_dir) || occursin("warmup", s)
        @info "skipping $s"
        continue
    end
    @info "Running experiments in directory $s"
    @info "Loading experiments.jser"
    experiments = deserialize(open(joinpath(exp_dir, "experiments.jser")))
    append!(all_experiments, experiments)
end

shuffle!(all_experiments)

@info "Warming up workers."
try
    @everywhere begin
        WARMUP = true
        include($config_file)
        warmup_experiments = deserialize(open(joinpath(data_output_root, "warmup", "experiments.jser")))
        print("Starting warmup.")
        map(run_experiment, warmup_experiments)
        print("Finished warmup.")
    end
catch e
    @warn e
    rmprocs(workers())
    @info "Warmup failed. Terminating."
    exit()
end
@info "Warmup done."

@everywhere WARMUP = false

try
    @info "Running $(length(all_experiments)) experiments."
    start = now()
    @info "$(start) - Running experiments..."
    pmap(run_experiment, all_experiments, on_error=ex->print("!!! ", ex))
    finish = now()
    @info "$finish - Done."
    @info "Ran $(length(all_experiments)) experiment(s) in $(canonicalize(Dates.CompoundPeriod(finish - start)))"
catch e
    @warn e
finally
    rmprocs(workers())
end

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
module NeuralNet

    using Random

    include("neuralnet.jl")
    include("neuralbenchmarks.jl")
    include("neuralactivations.jl")
    include("neuralplots.jl")
    include("models.jl")
    include("costfunctions.jl")

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TerminalMenus

terminal = nothing  # The user terminal

import REPL

function __init__()
    global terminal
    terminal = REPL.Terminals.TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)
end

include("util.jl")
include("config.jl")

include("AbstractMenu.jl")
include("RadioMenu.jl")
include("MultiSelectMenu.jl")

export
    RadioMenu,
    MultiSelectMenu,
    request

end # module

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
module TestExamplesP4estMesh2D

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "P4estMesh2D" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      # Expected errors are exactly the same as with TreeMesh!
      l2   = [8.311947673061856e-6],
      linf = [6.627000273229378e-5])
  end

  @trixi_testset "elixir_advection_nonconforming_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming_flag.jl"),
      l2   = [3.198940059144588e-5],
      linf = [0.00030636069494005547])
  end

  @trixi_testset "elixir_advection_unstructured_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"),
      l2   = [0.0005379687442422346],
      linf = [0.007438525029884735])
  end

  @trixi_testset "elixir_advection_amr_solution_independent.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent.jl"),
      # Expected errors are exactly the same as with StructuredMesh!
      l2   = [4.949660644033807e-5],
      linf = [0.0004867846262313763])
  end

  @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"),
      l2   = [0.0012766060609964525],
      linf = [0.01750280631586159])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [4.507575525876275e-6],
      linf = [6.21489667023134e-5])
  end

  @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),
    l2   = [0.0034516244508588046, 0.0023420334036925493, 0.0024261923964557187, 0.004731710454271893],
    linf = [0.04155789011775046, 0.024772109862748914, 0.03759938693042297, 0.08039824959535657])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2e-12, 4.8e-12, 4e-12],
      atol = 2.0e-12, # required to make CI tests pass on macOS
    )
  end

  @trixi_testset "elixir_euler_shockcapturing_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_ec.jl"),
      l2   = [9.53984675e-02, 1.05633455e-01, 1.05636158e-01, 3.50747237e-01],
      linf = [2.94357464e-01, 4.07893014e-01, 3.97334516e-01, 1.08142520e+00],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
      l2   = [3.76149952e-01, 2.46970327e-01, 2.46970327e-01, 1.28889042e+00],
      linf = [1.22139001e+00, 1.17742626e+00, 1.17742626e+00, 6.20638482e+00],
      tspan = (0.0, 0.3))
  end

  @trixi_testset "elixir_euler_blast_wave_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
      l2   = [6.32183914e-01, 3.86914231e-01, 3.86869171e-01, 1.06575688e+00],
      linf = [2.76020890e+00, 2.32659890e+00, 2.32580837e+00, 2.15778188e+00],
      tspan = (0.0, 0.3))
  end

  @trixi_testset "elixir_eulergravity_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
      l2   = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697],
      linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.0513414461545583e-5, 1.0517900957166411e-6, 1.0517900957304043e-6, 1.511816606372376e-6,
              1.0443997728645063e-6, 7.879639064990798e-7, 7.879639065049896e-7, 1.0628631669056271e-6,
              4.3382328912336153e-7],
      linf = [4.255466285174592e-5, 1.0029706745823264e-5, 1.0029706747467781e-5, 1.2122265939010224e-5,
              5.4791097160444835e-6, 5.18922042269665e-6, 5.189220422141538e-6, 9.552667261422676e-6,
              1.4237578427628152e-6])
  end

  @trixi_testset "elixir_mhd_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
      l2   = [0.4552084651735862, 0.8918048264575757, 0.832471223081887, 0.0,
              0.9801264164951583, 0.10475690769435382, 0.1555132409149897, 0.0,
              2.0597079362763556e-5],
      linf = [10.194181233788775, 18.25472397868819, 10.031307436191334, 0.0,
              19.647239392277378, 1.3938810140985936, 1.8724965294853084, 0.0,
              0.0016290067532561904],
      tspan = (0.0, 0.02))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # module


text:
function get_acceleration(obj::PropulsionModel, arg0::FieldSpacecraftState, arg1::FieldAttitude, arg2::Vector{RealFieldElement})
    return jcall(obj, "getAcceleration", FieldVector3D, (FieldSpacecraftState, FieldAttitude, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function get_acceleration(obj::PropulsionModel, arg0::SpacecraftState, arg1::Attitude, arg2::Vector{jdouble})
    return jcall(obj, "getAcceleration", Vector3D, (SpacecraftState, Attitude, Vector{jdouble}), arg0, arg1, arg2)
end

function get_mass_derivatives(obj::PropulsionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "getMassDerivatives", RealFieldElement, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_mass_derivatives(obj::PropulsionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "getMassDerivatives", jdouble, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_name(obj::PropulsionModel)
    return jcall(obj, "getName", JString, ())
end


text:
module Front
import JSON, JSON2, PlotlyBase

to_dash(t::Any) = t
from_dash(::Type{Any}, t::Any) = t

function from_dash(::Type{PlotlyBase.Plot}, t)
    data = PlotlyBase.GenericTrace[PlotlyBase.GenericTrace(tr) for tr in t.data]
    layout = PlotlyBase.Layout(t.layout)
    PlotlyBase.Plot(data, layout)
end

#FIXME It's not elegant, but I didn't find elegant solution
JSON2.write(io::IO, p::PlotlyBase.Plot; kwargs...) = write(io, JSON.json(p))
end


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, DiffEqDevTools

using ODEInterfaceDiffEq
using Base.Test

const CPU_FLOPS = peakflops()
const LONGER_TESTS = false
#Start Test Script

tic()
@time @testset "Discrete Tests" begin include("discrete_algorithm_test.jl") end
@time @testset "Linear Tests" begin include("ode/ode_twodimlinear_tests.jl") end
@time @testset "Convergence Tests" begin include("ode/ode_convergence_tests.jl") end
@time @testset "Adaptive Tests" begin include("ode/ode_adaptive_tests.jl") end
@time @testset "Tstops Tests" begin include("ode/ode_tstops_tests.jl") end
@time @testset "Backwards Tests" begin include("ode/ode_backwards_test.jl") end
(LONGER_TESTS) && @time @testset "Unrolled Tests" begin include("ode/ode_unrolled_comparison_tests.jl") end
@time @testset "Initial Dt Tests" begin include("ode/ode_initdt_tests.jl") end
@time @testset "Rosenbrock Tests" begin include("ode/ode_rosenbrock_tests.jl") end
@time @testset "Dense Tests" begin include("ode/ode_dense_tests.jl") end
@time @testset "In-Place Tests" begin include("ode/ode_inplace_tests.jl") end
@time @testset "Events Tests" begin include("ode/ode_event_tests.jl") end
@time @testset "Cache Tests" begin include("ode/ode_cache_tests.jl") end
@time @testset "saveat Tests" begin include("ode/ode_saveat_tests.jl") end
(LONGER_TESTS) && @time @testset "Feagin Tests" begin include("ode/ode_feagin_tests.jl") end
@time @testset "Number Type Tests" begin include("ode/ode_numbertype_tests.jl") end
@time @testset "Data Array Tests" begin include("data_array_test.jl") end
@time @testset "Ndim Complex Tests" begin include("ode/ode_ndim_complex_tests.jl") end
@time @testset "Iterator Tests" begin include("iterator_tests.jl") end
@time @testset "Composite Algorithm Tests" begin include("composite_algorithm_test.jl") end

(LONGER_TESTS) && (@time @testset "Units Tests" begin include("units_tests.jl") end) # Too long for AppVeyor

toc()


text:
# Add PATH (adapt!)
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Leaf/");
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Utils/");
#push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/tools/");

# Use Plots:
using Plots
#using PhotoStructs
pyplot()

using CanopyRTMod
using BenchmarkTools
using Statistics

wl = CanopyRTMod.wl;
wle = CanopyRTMod.wle;
wlf = CanopyRTMod.wlf;

arrayOfLeaves = Array{leafbio{FT,length(wl), length(wle), length(wlf),length(wle)*length(wlf)}, 1}(undef, CanopyRTMod.canopy.nlayers)
for i = 1:CanopyRTMod.canopy.nlayers
    #@show i
    arrayOfLeaves[i] = leafbio{FT, length(wl), length(wle), length(wlf),length(wle)*length(wlf)}()
    CanopyRTMod.fluspect!(arrayOfLeaves[i], CanopyRTMod.optis)
end
# This actually almost takes 1ms!
#CanopyRTMod.fluspect!(leaf, CanopyRTMod.optis)

CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt)

CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);

@btime CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);

@time CanopyRTMod.computeSIF_Fluxes!(arrayOfLeaves, canOpt, canRad, canopy, CanopyRTMod.soil);
#piLo1, piLo2, piLo3, piLo4, F⁻,F⁺,S⁻,S⁺,piLs, piLd =
#println(size(piLs))
#plot(wlf, piLs)
plot(wlf, canRad.SIF_hemi/pi, label="Hemispheric Fluorescence (/pi)")
plot!(wlf, canRad.SIF_obs_sunlit, label="Sunlit direct")
plot!(wlf, canRad.SIF_obs_shaded, label="Shaded direct")
plot!(wlf, canRad.SIF_obs_scattered, label="Scattered")
plot!(wlf, canRad.SIF_obs_soil, label="Soil")

plot(wlf, canRad.SIF_hemi, label="Outgoing hemispheric")
plot!(wlf, canRad.SIF_sum, label="Sum of all layer SIF sources")

plot(wlf, canRad.SIF_hemi./canRad.SIF_sum, label="Outgoing hemispheric/sources")
plot!(wlf, pi*canRad.SIF_obs./canRad.SIF_sum, label="π*SIF_obs/sources")
#plot!(wlf, , label="Sum of all layer SIF sources")

plot(wlf, arrayOfLeaves[1].Mb[:,1])

iLAI = canopy.LAI/canopy.nlayers
ϵ = zeros(20).+0.98
fSun = (canOpt.Ps[1:20]+canOpt.Ps[2:21])/2
wlii = [10]
S⁺,S⁻ = CanopyRTMod.computeThermalFluxes(canRad.T_shade, canRad.T_sun, ϵ, iLAI, canopy.lidf, fSun, wlii)
#plot(S⁺)

@btime CanopyRTMod.deriveCanopyFluxes!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil, arrayOfLeaves);

plot(wl,canRad.alb_direct, label="Direct hemispheric albedo")
plot!(wl,canRad.alb_diffuse, label="Diffuse hemispheric albedo")
plot!(wl,canRad.alb_obs, label="observed directional (nadir) albedo")

sumLAI = [0:iLAI:canopy.LAI;]
plot(canOpt.Pso,-sumLAI, label="Pso")
plot!(canOpt.Ps,-sumLAI, label="Ps")
plot!(canOpt.Po,-sumLAI, label="Po")
ylabel!("-Cumulative LAI")
xlabel!("Probability ")
title!("solar Ps, outgoing Po, in-out Pso likelihoods")

solar_in = sunRad.E_diffuse+sunRad.E_direct;
soil_absorbed = canRad.E_down[:,end]+canOpt.Es_[:,end]-canRad.E_up[:,end]
canopy_absorbed_diff = sum(canRad.netSW_shade, dims=2)[:,1]
canopy_absorbed_dir = sum(canRad.netSW_sunlit, dims=2)[:,1]
plot(wl,solar_in-CanopyRTMod.canRad.E_up[:,1].+10, label="Net TOA incoming (+10)")
plot!(wl,soil_absorbed, label="Absorbed  by soil")
plot!(wl,soil_absorbed+canopy_absorbed_diff, label="Absorbed  by soil + canopy (diffuse)")
plot!(wl,soil_absorbed+canopy_absorbed_diff+canopy_absorbed_dir, label="Absorbed  by soil + canopy (total)")
xlabel!("Wavelength (nm)")
ylabel!("Radiance (mW m-2 μm-1)")

sum(soil_absorbed+canopy_absorbed_diff+canopy_absorbed_dir-solar_in+CanopyRTMod.canRad.E_up[:,1], dims=1)

plot(wl, arrayOfLeaves[10].ρ_SW, label="Reflectance")
plot!(wl, 1 .-arrayOfLeaves[10].τ_SW, label="1-Transmission")

println("Net Soil radiation (direct): ", canRad.RnSoil_direct, "W/m2")
println("Net Soil radiation (diffuse): ", canRad.RnSoil_diffuse, "W/m2")

println("Incoming direct PAR: ", 1e6*canRad.incomingPAR_direct, "moles m^-2 s^-1")
println("Incoming diffuse PAR: ", 1e6*canRad.incomingPAR_diffuse, "moles m^-2 s^-1")

nl = CanopyRTMod.canopy.nlayers
#normi = 1/mean(canOpt.fs'*canopy.lidf)
plot(sumLAI[1:nl], 1e6*canRad.absPAR_shadeCab, label="Diffuse absorbed PAR")
plot!(sumLAI[1:nl], 1e6*canRad.absPAR_sunCab[5,4,:].*canOpt.Ps[1:nl], label="Direct absorbed PAR*Ps")
plot!(sumLAI[1:nl], 1e6*canRad.absPAR_sunCab[5,4,:], label="Direct absorbed PAR")

mean(canOpt.Ps)

contourf(1e6*canRad.absPAR_sunCab[:,:,1])
#plot!(1e6*canRad.absPAR_sun[:,:,1][:])

# plot(wl[CanopyRTMod.iPAR], arrayOfLeaves[1].kChlrel[CanopyRTMod.iPAR], label="Relative absorbed light by Cab+Car")
plot(wl[CanopyRTMod.iPAR], arrayOfLeaves[1].kChlrel_old[CanopyRTMod.iPAR], label="Relative absorbed light by Cab")

nwl2,nl2 = size(canOpt.R_dd)

nwl2

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module DataAnalysisTools

using LsqFit, Plots, Measurements
import LsqFit: curve_fit
export modify,crop,mask,plotData!,modelWeights,analyze

function modify(data; xfunc = x->x, yfunc = y->y)
    xdat = data[1]; ydat = data[2]
    if length(data) > 2
        ydatm = ydat .± data[3]
    else
        ydatm = measurement().(ydat)
    end

    xdat = map(xfunc,xdat)
    ydatm = map(yfunc,ydatm)
    yvals = map(y->y.val,ydatm)
    yerrs = map(y->y.err,ydatm)

    return [x,yvals,yerrs]
end

#Keeps only data within xmin,xmax and within ymin,ymax
function crop(data; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf)
    x = data[1]; y = data[2];
    xindices = findall(sx->sx .>= xmin && sx .<= xmax,x)
    yindices = findall(sy->sy .>= ymin && sy .<= ymax,y)
    newindices = intersect(xindices,yindices)
    xnew = x[newindices]
    ynew = y[newindices]
    newdata = [xnew,ynew]
    if length(data) > 2
        err = data[3]
        push!(newdata,err[newindices])
    end
    return newdata
end

#Removes all data within xmin,xmax and within ymin,ymax
function mask(data; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf)
    x = data[1]; y = data[2];
    xindices = findall(sx->sx .< xmin && sx .> xmax,x)
    yindices = findall(sy->sy .< ymin && sy .> ymax,y)
    newindices = union(xindices,yindices)
    xnew = x[newindices]
    ynew = y[newindices]
    newdata = [xnew,ynew]
    if length(data) > 2
        err = data[3]
        push!(newdata,err[newindices])
    end
    return newdata
end

function plotData!(plt,data;kwargs...)
    x = data[1]; y = data[2]
    if length(data) < 3
        Plots.plot!(plt, x,y;kwargs...)
    else
        yerr = data[3]
        Plots.plot!(plt, x,y .± yerr;kwargs...)
    end
end

function modelWeights(data)
    weights = ones(length(data[1]))
    if length(data) > 2
        weights = (1 ./ (data[3].^2))
    end
end

function analyze(data,model)
    x,y,yerr = data
    weights = modelWeights(data)
    #fitting
    fit = curve_fit(model,x,y,weights,[1.0,1.0])
    modely = model(x,fit.param)
    modelData = [x,modely]
    sigma = stderror(fit)
    params = fit.param .± sigma

    return [params,modelData]
end

end #module DataAnalysisTools

max_stars_repo_name: sphanna/DataAnalysisTools

max_issues_repo_name: sphanna/DataAnalysisTools

max_forks_repo_name: sphanna/DataAnalysisTools


text:
using GPUArrays
using Base.Test, GPUArrays.TestSuite

function run_linalg(Typ)
    T = Typ{Float32}
    @testset "Linalg" begin
        @testset "transpose" begin
            against_base(ctranspose, T, (32, 32))
        end
        @testset "PermuteDims" begin
            against_base(x -> permutedims(x, (2, 1)), T, (2, 3))
            against_base(x -> permutedims(x, (2, 1, 3)), T, (4, 5, 6))
            against_base(x -> permutedims(x, (3, 1, 2)), T, (4, 5, 6))
        end
    end
end


text:
using PyCall:PyObject, PyReverseDims, pyimport_conda

struct PyModel <: Model
    pyobject::PyObject
end

# lazily import a python model
pyod_import(name::String) = () -> begin
    module_name = lowercase(name)
    # names are not consistently uppercased
    model_name = name == "IForest" ? name : uppercase(name)
    getproperty(pyimport_conda("pyod.models.$(module_name)", "pyod", "conda-forge"), model_name)
end
pyod_import(name::Symbol) = pyod_import(String(name))

make_docs_link(name::String) =
    "<https://pyod.readthedocs.io/en/latest/pyod.models.html#module-pyod.models.$(lowercase(name))>"

"""
    pyod_fit(modelname, params)

Implements the `fit` method for an underlying python model.
"""
function pyod_fit(modelname, params)
    pymodelname = String(modelname)[3:end]
    quote
        function fit(model::$modelname, X::Data)::Fit
            Xt = PyReverseDims(X) # from column-major to row-major
            # load the underlying python model with key-word arguments
            detector = pyod_import($pymodelname)()($((Expr(:kw, p, :(model.$p)) for p in params)...))
            detector.fit(Xt)
            # the underlying python model is out model
            return Fit(PyModel(detector), detector.decision_scores_)
        end
    end
end

"""
    pyod_score(modelname)

Implements the `score` method for an underlying python model.
"""
function pyod_score(modelname)
    quote
        function score(_::$modelname, fitresult::Fit, X::Data)::Result
            Xt = PyReverseDims(X) # change from column-major to row-major
            scores_test = fitresult.model.pyobject.decision_function(Xt)
            return fitresult.scores, scores_test
        end
    end
end

"""
    py_constructor(expr)
Extracts the relevant information from the expr and build the expression
corresponding to the model constructor (see [`_model_constructor`](@ref)).
"""
function py_constructor(expr)
    # similar to @mlj_model
    expr, modelname, params, defaults, constraints = MMI._process_model_def(@__MODULE__, expr)

    # keyword constructor
    const_expr = MMI._model_constructor(modelname, params, defaults)

    # associate the constructor with the definition of the struct
    push!(expr.args[3].args, const_expr)

    # cleaner
    clean_expr = MMI._model_cleaner(modelname, defaults, constraints)

    # return
    return modelname, params, clean_expr, expr
end

macro pymodel(expr)
    modelname, params, clean_expr, expr = py_constructor(expr)
    @assert startswith(String(modelname), "Py") "A python model name has to start with `Py`, e.g. PyKNN"
    expr = quote
        Base.@__doc__($expr)
        $clean_expr
        $(pyod_fit(modelname, params))
        $(pyod_score(modelname))
        MMI.metadata_pkg($modelname,
            package_name="OutlierDetection.jl",
            package_uuid="262411bb-c475-4342-ba9e-03b8c0183ca6",
            package_url="https://github.com/davnn/OutlierDetection.jl",
            is_pure_julia=false,
            package_license="MIT",
            is_wrapper=false
        )
    end
    esc(expr)
end


text:
using Images

function loadDataFromFiles!(fileDir,fileList,Dataset;hist_eq=false)
  for sampleIdx = 1:size(Dataset,2)
    imFile = joinpath(fileDir,fileList[sampleIdx])
    im = load(imFile)
    ImgData = convert(Array{Float64,2},im.data)
    if hist_eq
      ImgData = histogramEqualization(ImgData)
    end
    Dataset[:,sampleIdx] = ImgData[:]
  end
end

function histogramEqualization(im)
  # println("=================EQ====================")
  L  = 256   # Intesnity levels
  Lm = 255   # L-1
  N = length(im)
  inc = 1./N;

  # First, convert everything back to integers [0,255]
  im = floor(im.*255.0)
  # println(im)
  UIim = convert(Array{UInt8,2},im)

  # Next, we want the PMF
  P = zeros(L,1)
  for i=1:N
    P[UIim[i]+1] += inc;
  end

  gim = copy(im)
  # Now apply equalizaitons
  for i=1:N
    gim[i] = Lm.*sum(P[1:(UIim[i]+1)])
  end

  return gim
end

function traindata(;hist_eq=false)
  Pixels  = 19*19

  # Load Face Data
  FaceDataDir = joinpath(dirname(@__FILE__),"..","data","train","face")
  FaceFiles   = readdir(FaceDataDir)
  FaceDataset = Array(Float64,Pixels,length(FaceFiles))
  loadDataFromFiles!(FaceDataDir,FaceFiles,FaceDataset;hist_eq=hist_eq)

  # Load Non-Face Data
  NonFaceDataDir = joinpath(dirname(@__FILE__),"..","data","train","non-face")
  NonFaceFiles = readdir(NonFaceDataDir)
  NonFaceDataset = Array(Float64,Pixels,length(NonFaceFiles))
  loadDataFromFiles!(NonFaceDataDir,NonFaceFiles,NonFaceDataset;hist_eq=hist_eq)

  return FaceDataset, NonFaceDataset
end

function testdata(;hist_eq=false)
  Pixels  = 19*19

  # Load Face Data
  FaceDataDir = joinpath(dirname(@__FILE__),"..","data","test","face")
  FaceFiles   = readdir(FaceDataDir)
  FaceDataset = Array(Float64,Pixels,length(FaceFiles))
  loadDataFromFiles!(FaceDataDir,FaceFiles,FaceDataset;hist_eq=hist_eq)

  # Load Non-Face Data
  NonFaceDataDir = joinpath(dirname(@__FILE__),"..","data","test","non-face")
  NonFaceFiles = readdir(NonFaceDataDir)
  NonFaceDataset = Array(Float64,Pixels,length(NonFaceFiles))
  loadDataFromFiles!(NonFaceDataDir,NonFaceFiles,NonFaceDataset;hist_eq=hist_eq)

  return FaceDataset, NonFaceDataset
end