text:
function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DiscreteCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DiscreteConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP5ConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<4 || calcVal
    @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,b1,b3,b4,b5,b6,b7,c1,c2,c3,c4,c5,c6 = cache
    @unpack d1,d3,d4,d5,d6,d7 = cache
    k1 = f(t,uprev)
    k2 = f(t+c1*dt,uprev+dt*(a21*k1))
    k3 = f(t+c2*dt,uprev+dt*(a31*k1+a32*k2))
    k4 = f(t+c3*dt,uprev+dt*(a41*k1+a42*k2+a43*k3))
    k5 = f(t+c4*dt,uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4))
    k6 = f(t+dt,uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5))
    update = a71*k1+a73*k3+a74*k4+a75*k5+a76*k6
    k7 = f(t+dt,uprev+dt*update)
    copyat_or_push!(k,1,update)
    bspl = k1 - update
    copyat_or_push!(k,2,bspl)
    copyat_or_push!(k,3,update - k7 - bspl)
    copyat_or_push!(k,4,d1*k1+d3*k3+d4*k4+d5*k5+d6*k6+d7*k7)
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP5Cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<4 || calcVal
    @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,b1,b3,b4,b5,b6,b7,c1,c2,c3,c4,c5,c6 = cache.tab
    @unpack d1,d3,d4,d5,d6,d7 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
    uidx = eachindex(uprev)
    f(t,uprev,k1)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a21*k1[i])
    end
    f(t+c1*dt,tmp,k2)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    f(t+c2*dt,tmp,k3)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
    end
    f(t+c3*dt,tmp,k4)
    for i in uidx
      tmp[i] =uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
    end
    f(t+c4*dt,tmp,k5)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    f(t+dt,tmp,k6)
    for i in uidx
      update[i] = a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]
      tmp[i] = uprev[i]+dt*update[i]
    end
    f(t+dt,tmp,k7)
    copyat_or_push!(k,1,update)
    for i in uidx
      bspl[i] = k1[i] - update[i]
      dense_tmp3[i] = update[i] - k7[i] - bspl[i]
      dense_tmp4[i] = (d1*k1[i]+d3*k3[i]+d4*k4[i]+d5*k5[i]+d6*k6[i]+d7*k7[i])
    end
    copyat_or_push!(k,2,bspl)
    copyat_or_push!(k,3,dense_tmp3)
    copyat_or_push!(k,4,dense_tmp4)
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP5ThreadedCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<4 || calcVal
    @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,b1,b3,b4,b5,b6,b7,c1,c2,c3,c4,c5,c6 = cache.tab
    @unpack d1,d3,d4,d5,d6,d7 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
    uidx = eachindex(uprev)
    f(t,uprev,k1)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a21*k1[i])
    end
    f(t+c1*dt,tmp,k2)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    f(t+c2*dt,tmp,k3)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
    end
    f(t+c3*dt,tmp,k4)
    for i in uidx
      tmp[i] =uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
    end
    f(t+c4*dt,tmp,k5)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    f(t+dt,tmp,k6)
    for i in uidx
      update[i] = a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]
      tmp[i] = uprev[i]+dt*update[i]
    end
    f(t+dt,tmp,k7)
    copyat_or_push!(k,1,update)
    for i in uidx
      bspl[i] = k1[i] - update[i]
      dense_tmp3[i] = update[i] - k7[i] - bspl[i]
      dense_tmp4[i] = (d1*k1[i]+d3*k3[i]+d4*k4[i]+d5*k5[i]+d6*k6[i]+d7*k7[i])
    end
    copyat_or_push!(k,2,bspl)
    copyat_or_push!(k,3,dense_tmp3)
    copyat_or_push!(k,4,dense_tmp4)
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Tsit5ConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<7 || calcVal
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,b1,b2,b3,b4,b5,b6,b7 = cache
    copyat_or_push!(k,1,f(t,uprev))
    copyat_or_push!(k,2,f(t+c1*dt,uprev+dt*(a21*k[1])))
    copyat_or_push!(k,3,f(t+c2*dt,uprev+dt*(a31*k[1]+a32*k[2])))
    copyat_or_push!(k,4,f(t+c3*dt,uprev+dt*(a41*k[1]+a42*k[2]+a43*k[3])))
    copyat_or_push!(k,5,f(t+c4*dt,uprev+dt*(a51*k[1]+a52*k[2]+a53*k[3]+a54*k[4])))
    copyat_or_push!(k,6,f(t+dt,uprev+dt*(a61*k[1]+a62*k[2]+a63*k[3]+a64*k[4]+a65*k[5])))
    utmp = uprev+dt*(a71*k[1]+a72*k[2]+a73*k[3]+a74*k[4]+a75*k[5]+a76*k[6])
    copyat_or_push!(k,7,f(t+dt,utmp))
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Tsit5Cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<7 || calcVal
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,b1,b2,b3,b4,b5,b6,b7 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp = cache
    uidx = eachindex(uprev)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a21*k1[i])
    end
    f(t+c1*dt,tmp,k2)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    f(t+c2*dt,tmp,k3)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
    end
    f(t+c3*dt,tmp,k4)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
    end
    f(t+c4*dt,tmp,k5)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    f(t+dt,tmp,k6)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
    end
    f(t+dt,u,k7)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
  end
  nothing
end

"""
An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28

Called to add the extra k9, k10, k11 steps for the Order 5 interpolation when needed
"""
function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::BS5ConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k) < 8 || calcVal
    @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache
    copyat_or_push!(k,1,f(t,uprev))
    copyat_or_push!(k,2,f(t+c1*dt,uprev+dt*a21*k[1]))
    copyat_or_push!(k,3,f(t+c2*dt,uprev+dt*(a31*k[1]+a32*k[2])))
    copyat_or_push!(k,4,f(t+c3*dt,uprev+dt*(a41*k[1]+a42*k[2]+a43*k[3])))
    copyat_or_push!(k,5,f(t+c4*dt,uprev+dt*(a51*k[1]+a52*k[2]+a53*k[3]+a54*k[4])))
    copyat_or_push!(k,6,f(t+c5*dt,uprev+dt*(a61*k[1]+a62*k[2]+a63*k[3]+a64*k[4]+a65*k[5])))
    copyat_or_push!(k,7,f(t+dt,uprev+dt*(a71*k[1]+a72*k[2]+a73*k[3]+a74*k[4]+a75*k[5]+a76*k[6])))
    copyat_or_push!(k,8,f(t+dt,uprev+dt*(a81*k[1]+a83*k[3]+a84*k[4]+a85*k[5]+a86*k[6]+a87*k[7])))
  end
  if (calcVal2 && length(k)< 11) || calcVal3 # Have not added the extra stages yet
    @unpack c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = cache
    copyat_or_push!(k,9,f(t+c6*dt,uprev+dt*(a91*k[1]+a92*k[2]+a93*k[3]+a94*k[4]+a95*k[5]+a96*k[6]+a97*k[7]+a98*k[8])))
    copyat_or_push!(k,10,f(t+c7*dt,uprev+dt*(a101*k[1]+a102*k[2]+a103*k[3]+a104*k[4]+a105*k[5]+a106*k[6]+a107*k[7]+a108*k[8]+a109*k[9])))
    copyat_or_push!(k,11,f(t+c8*dt,uprev+dt*(a111*k[1]+a112*k[2]+a113*k[3]+a114*k[4]+a115*k[5]+a116*k[6]+a117*k[7]+a118*k[8]+a119*k[9]+a1110*k[10])))
  end
  nothing
end

"""
An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28

Called to add the extra k9, k10, k11 steps for the Order 5 interpolation when needed
"""
function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::BS5Cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k) < 8 || calcVal
    uidx = eachindex(uprev)
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,utilde,uhat,tmp,atmp,atmptilde = cache
    @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache.tab
    for i in uidx
      tmp[i] = uprev[i]+dt*a21*k1[i]
    end
    f(t+c1*dt,tmp,k2)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    f(t+c2*dt,tmp,k3)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
    end
    f(t+c3*dt,tmp,k4)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
    end
    f(t+c4*dt,tmp,k5)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    f(t+c5*dt,tmp,k6)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
    end
    f(t+dt,tmp,k7)
    for i in uidx
      u[i] = uprev[i]+dt*(a81*k1[i]+a83*k3[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i])
    end
    f(t+dt,u,k8)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
  end
  if (calcVal2 && length(k)< 11) || calcVal3 # Have not added the extra stages yet
    uidx = eachindex(uprev)
    rtmp = similar(cache.k1)
    @unpack tmp = cache
    @unpack c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = cache.tab
    for i in uidx
      tmp[i] = uprev[i]+dt*(a91*k[1][i]+a92*k[2][i]+a93*k[3][i]+a94*k[4][i]+a95*k[5][i]+a96*k[6][i]+a97*k[7][i]+a98*k[8][i])
    end
    f(t+c6*dt,tmp,rtmp); copyat_or_push!(k,9,rtmp)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a101*k[1][i]+a102*k[2][i]+a103*k[3][i]+a104*k[4][i]+a105*k[5][i]+a106*k[6][i]+a107*k[7][i]+a108*k[8][i]+a109*k[9][i])
    end
    f(t+c7*dt,tmp,rtmp); copyat_or_push!(k,10,rtmp)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a111*k[1][i]+a112*k[2][i]+a113*k[3][i]+a114*k[4][i]+a115*k[5][i]+a116*k[6][i]+a117*k[7][i]+a118*k[8][i]+a119*k[9][i]+a1110*k[10][i])
    end
    f(t+c8*dt,tmp,rtmp); copyat_or_push!(k,11,rtmp,Val{false})
  end
  nothing
end


text:
module recombination

using StatsBase

# Physical constants / units
const kb=8.6173324E-5 # in units of eV
const ε0=8.854E-12 #Units: C2N−1m−2
const D=3.336E-30 #Debye in SI

const Å=1E-10 # Angstrom

# Specific to MAPI
const r=6.29Å # Sensible figure for MA for consistency; AW
const D1=2.18D #Methyl-Ammonia; b3lyp/6-31g*

#Read in potential...
function starrynight_read_potential(filename)
    lattice=readdlm(filename)
    pot=lattice[:,4] # just potential, in internal units of Starrynight

    T=300           # Temperature; Kelvin
    β=1/(kb*T)      # Thermodynamic Beta; units

    # OK; let's calc Dipole Potential with simple form, contribution from one dipole at distance 1
    factor=1/(4*pi*ε0) * (D1*r)/(r^3)
    @printf STDERR "Conversion factor from Starrynight potential internal units --> V : %f\n" factor
    # Scale from Starrynight potential --> volts
    pot*=factor

    # Apply effective dielectric screening...
    # Calculated value DFT: 5.6 --> 6.5 (dep. on MA alignment)
    # See: APL Mat. 1, 042111 (2013); http://dx.doi.org/10.1063/1.4824147
    # Measured value, Ellipso in NIR: 5.0
    ɛr=4.1
    pot/=ɛr

    N=length(pot) # number of elements

    @printf STDERR "OK; data loaded and converted to Volts. Variance: %e SD: %e\n\n" var(pot) std(pot)
#    print(hist(pot,50)) # I think this used to make a nice little histogram; API has changed

    describe(pot)

    return N, pot
end

#N=1000
#pot=0.13*randn(N) # Normal distribution
#pot=0.13*sin(range(0,N)) # Sinusoidal variations

function calc_recombination(N,pot,kernel)
    # Calculate partition function by sum over state energies
    Zh=sum(kernel(-pot))
    Ze=sum(kernel(pot))
    # Charge densities from partition fn.
    ρh=(1/Zh) * kernel(-pot)
    ρe=(1/Ze) * kernel(pot)

    @printf STDERR "Zh: %e Ze: %e e: %e h: %e\n" Zh Ze sum(ρe) sum(ρh)
    # Factors of N make these values independent of number of elements in the
    # potential, self-consistent, and equal to R=1 for a flat potential.
    @printf STDERR "R=N*N/(Zh*Ze): %e\t" N*N/(Zh*Ze)
    @printf STDERR "R=N*sum(ρe.*ρh): %e\n" N*sum(ρh.*ρe)
    return N*sum(ρh.*ρe)
end

function calc_mobility(N::Int,pot)
    potsorted=sort(reshape(pot,length(pot)))
    trap=potsorted[div(N,4)]-potsorted[1] # Diff in energy between deepest point and first quartile, 25% CDF
    @printf "Trap Depth: %e Bottom(eV): %e 25\%%(eV): %e\n" trap potsorted[1] potsorted[div(N,4)]
    # OK; now need some FD distribution to get trapped vs. free charges.

    T=300           # Temperature; Kelvin
    β=1/(kb*T)      # Thermodynamic Beta; units

    @printf "Boltzmann ratio given trap 0->25\%% trap depth @ 300 K %e \n" exp(trap*β)
end

end

ext: jl

max_stars_repo_name: WMD-Bath/StarryNight

max_issues_repo_name: WMD-Bath/StarryNight

max_forks_repo_name: WMD-group/StarryNight


text:
# Combined Boxplot and ViolinPlot

#JL Plots
using StatsPlots
# Pkg.add("RDatasets")
import RDatasets
singers = RDatasets.dataset("lattice", "singer")
@df singers Plots.violin(:VoicePart, :Height, line = 0, fill = (0.2, :blue))
@df singers Plots.boxplot!(:VoicePart, :Height, line = (2, :black), fill = (0.3, :orange))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
"""
    argmax_(A; dims)

Index of maximum elements of `A` over given `dims`, dropping the reduced `dims`.
"""
function argmax_(A::AbstractArray; dims)
    indices = argmax(A; dims=dims)
    return dropdims(indices; dims=dims)
end

"""
    encode(data, labels)

Returns a one-hot encoded version of `data`, where classes are indexed according
to `labels`.
"""
function encode(data::AbstractArray, labels = sort(unique(data)))
    return reshape(data, 1, size(data)...) .== labels
end

"""
    decode(X)

Returns an array `A` such that `OneHot.encode(A, 1:q) == X`.
Here `q = size(X, 1)` is interpreted as the number of classes.
"""
function decode(X::AbstractArray)
    I = argmax_(X; dims=1)
    return first.(Tuple.(I))
end

decode(X::OneHotArray) = X.c


text:
using MagmaThermoKinematics
using MagmaThermoKinematics.Units
using MagmaThermoKinematics.Diffusion2D
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
using Plots

# Initialize
@init_parallel_stencil(Threads, Float64, 2);    # initialize parallel stencil in 2D

#------------------------------------------------------------------------------------------
@views function MainCode_2D();
# Model parameters
W,H                     =   30km, 30km;             # Width, Height
ρ                       =   2800;                   # Density
cp                      =   1050;                   # Heat capacity
k_rock, k_magma         =   1.5, 1.2;               # Thermal conductivity of host rock & magma
La                      =   350e3;                  # Latent heat J/kg/K
GeoT                    =   20.0/km;                # Geothermal gradient [K/km]
x_in,z_in               =   20km,   -15km;          # Center of dike [x,z]
W_in, H_in              =   5km,    0.2km;          # Width and thickness of dike
T_in                    =   900;                    # Intrusion temperature
InjectionInterval       =   0.1kyr;                 # Inject a new dike every X kyrs
maxTime                 =   25kyr;                  # Maximum simulation time in kyrs
H_ran, W_ran            =   H*0.4, W*0.3;           # Size of domain in which we randomly place dikes and range of angles
DikeType                =   "ElasticDike"           # Type to be injected ("ElasticDike","SquareDike")

Nx, Nz                  =   500, 500;               # resolution (grid cells in x,z direction)
dx,dz                   =   W/(Nx-1), H/(Nz-1);     # grid size
κ                       =   k_rock./(ρ*cp);         # thermal diffusivity
dt                      =   min(dx^2, dz^2)./κ/10;  # stable timestep (required for explicit FD)
nt::Int64               =   floor(maxTime/dt);      # number of required timesteps
nTr_dike                =   300;                    # number of tracers inserted per dike

# Array initializations
T                       =   @zeros(Nx,   Nz);
K                       =   @ones(Nx,    Nz)*k_rock;
Rho                     =   @ones(Nx,    Nz)*ρ;
Cp                      =   @ones(Nx,    Nz)*cp;

# Work array initialization
Tnew, qx,qz, Kx, Kz     =   @zeros(Nx,   Nz), @zeros(Nx-1, Nz),     @zeros(Nx,   Nz-1), @zeros(Nx-1, Nz), @zeros(Nx,   Nz-1)    # thermal solver
X,Xc,Z                  =   @zeros(Nx,   Nz), @zeros(Nx-1, Nz-1),   @zeros(Nx,   Nz)    # 2D gridpoints
Phi_o, Phi, dPhi_dt     =   @zeros(Nx,   Nz), @zeros(Nx,   Nz  ),   @zeros(Nx,   Nz)    # solid fraction

# Set up model geometry & initial T structure
x,z                     =   (0:Nx-1)*dx, (-(Nz-1):0)*dz;                            # 1-D coordinate arrays
crd                     =   collect(Iterators.product(x,z))                         # Generate coordinates from 1D coordinate vectors
X,Z                     =   (x->x[1]).(crd), (x->x[2]).(crd);                       # Transfer coords to 3D arrays
Grid                    =   (x,z);                                                  # Grid
Tracers                 =   StructArray{Tracer}(undef, 1)                           # Initialize tracers
dike                    =   Dike(W=W_in,H=H_in,Type=DikeType,T=T_in);               # "Reference" dike with given thickness,radius and T
T                       .=   -Z.*GeoT;                                              # Initial (linear) temperature profile
Phi, dPhi_dt            =   SolidFraction(T, Phi_o, dt);                            # Compute solid fraction

# Preparation of visualisation
ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])

time, dike_inj, InjectVol, Time_vec,Melt_Time = 0.0, 0.0, 0.0,zeros(nt,1),zeros(nt,1);
for it = 1:nt   # Time loop

    if floor(time/InjectionInterval)> dike_inj                                      # Add new dike every X years
        dike_inj  =   floor(time/InjectionInterval)                                 # Keeps track on what was injected already
        cen       =   [W/2.; -H/2.] + rand(-0.5:1e-3:0.5, 2).*[W_ran;H_ran];        # Randomly vary center of dike
        if cen[end]<-12km;  Angle_rand = rand( 80.0:0.1:100.0)                      # Orientation: near-vertical @ depth
        else                Angle_rand = rand(-10.0:0.1:10.0); end                  # Orientation: near-vertical @ shallower depth
        dike      =   Dike(dike, Center=cen[:],Angle=[Angle_rand]);                 # Specify dike with random location/angle but fixed size/T
        Tracers, T, Vol     =   InjectDike(Tracers, T, Grid, dike, nTr_dike);       # Add dike, move hostrocks
        InjectVol           +=  Vol                                                 # Keep track of injected volume
        println("Added new dike; total injected magma volume = $(round(InjectVol/km³,digits=2)) km³; rate Q=$(round(InjectVol/(time),digits=2)) m³/s")
    end

    Phi, dPhi_dt     =  SolidFraction(T, Phi_o, dt);                                # Compute solid fraction
    K               .=  Phi.*k_rock .+ (1 .- Phi).*k_magma;                         # Thermal conductivity

    # Perform a diffusion step
    @parallel diffusion2D_step!(Tnew, T, qx, qz, K, Kx, Kz, Rho, Cp, dt, dx, dz, La, dPhi_dt);
    @parallel (1:size(T,2)) bc2D_x!(Tnew);                                                      # set lateral boundary conditions (flux-free)
    Tnew[:,1] .= GeoT*H; Tnew[:,end] .= 0.0;                                                    # bottom & top temperature (constant)

    Tracers             =   UpdateTracers(Tracers, Grid, Tnew, Phi);                            # Update info on tracers
    T, Tnew             =   Tnew, T;                                                            # Update temperature
    time                =   time + dt;                                                          # Keep track of evolved time
    Melt_Time[it]       =   sum( 1.0 .- Phi)/(Nx*Nz)                                            # Melt fraction in crust
    Time_vec[it]        =   time;                                                               # Vector with time
    println(" Timestep $it = $(round(time/kyr*100)/100) kyrs")

    if mod(it,20)==0  # Visualisation
        Phi_melt    =   1.0 .- Phi;
        p1          =   heatmap(x/km, z/km, T',         aspect_ratio=1, xlims=(x[1]/km,x[end]/km), ylims=(z[1]/km,z[end]/km),   c=:lajolla, clims=(0.,900.), xlabel="Width [km]",ylabel="Depth [km]", title="$(round(time/kyr, digits=2)) kyrs", dpi=200, fontsize=6, colorbar_title="Temperature")
        p2          =   heatmap(x/km, z/km, Phi_melt',  aspect_ratio=1, xlims=(x[1]/km,x[end]/km), ylims=(z[1]/km,z[end]/km),   c=:nuuk,    clims=(0., 1. ), xlabel="Width [km]",             dpi=200, fontsize=6, colorbar_title="Melt Fraction")
        plot(p1, p2, layout=(1,2)); frame(anim)
    end
end
gif(anim, "Example2D.gif", fps = 15)   # create gif animation
return Time_vec, Melt_Time;
end # end of main function

Time_vec,Melt_Time = MainCode_2D(); # start the main code
plot(Time_vec/kyr, Melt_Time, xlabel="Time [kyrs]", ylabel="Fraction of crust that is molten", label=:none); png("Time_vs_Melt_Example2D") #Create plot


text:
# BAM Writer
# ==========

"""
    BAM.Writer(output::BGZFStream, header::SAM.Header)

Create a data writer of the BAM file format.

# Arguments
* `output`: data sink
* `header`: SAM header object
"""
mutable struct Writer <: BioGenerics.IO.AbstractWriter
    stream::BGZFStreams.BGZFStream
end

function Writer(stream::BGZFStreams.BGZFStream, header::SAM.Header)
    refseqnames = String[]
    refseqlens = Int[]
    for metainfo in findall(header, "SQ")
        push!(refseqnames, metainfo["SN"])
        push!(refseqlens, parse(Int, metainfo["LN"]))
    end
    write_header(stream, header, refseqnames, refseqlens)
    return Writer(stream)
end

function BioGenerics.IO.stream(writer::Writer)
    return writer.stream
end

function Base.write(writer::Writer, record::Record)
    n = 0
    n += unsafe_write(writer.stream, pointer_from_objref(record), FIXED_FIELDS_BYTES)
    n += unsafe_write(writer.stream, pointer(record.data), data_size(record))
    return n
end

function write_header(stream, header, refseqnames, refseqlens)
    @assert length(refseqnames) == length(refseqlens)
    n = 0

    # magic bytes
    n += write(stream, "BAM\1")

    # SAM header
    buf = IOBuffer()
    l = write(SAM.Writer(buf), header)
    n += write(stream, Int32(l))
    n += write(stream, take!(buf))

    # reference sequences
    n += write(stream, Int32(length(refseqnames)))
    for (seqname, seqlen) in zip(refseqnames, refseqlens)
        namelen = length(seqname)
        n += write(stream, Int32(namelen + 1))
        n += write(stream, seqname, '\0')
        n += write(stream, Int32(seqlen))
    end

    return n
end


text:
println("starting!!!")
using DelimitedFiles
using TicketToRide
using Flux
using Flux.Optimise
using TicketToRide: train!

Exs = [-0.425803520728287]

Ls = [10]
n = 10

# nlayers = [50, 60, 70, 80, 90, 100]

nlayers = parse(Int, ARGS[1])

println("\nL is $n")
println("nlayers is $nlayers")

# pre-compile once
opt = Optimise.Descent(1e-3)
circuit = variational_circuit(n, nlayers);
ham = heisenberg1D(n)
dispatch!(circuit, :random);

# make directory if none exists
dir = "testing_tickets/L-$n/"
mkpath(dir)
file = dir * "layers-$nlayers.txt"

@info file

ind = findall(x -> x == n, Ls)
println(ind)

opt = Optimise.Descent(1e-2)
@time history = train!(opt, 1, ham, circuit; verbose=true)
history = train!(opt, 2000, ham, circuit; verbose=true, Eex = Exs[ind][1])
writedlm(file, history)

opt = Optimise.Descent(1e-3)
history = train!(opt, 2000, ham, circuit; verbose=true, Eex = Exs[ind][1])
writedlm(file, history)


text:
using NordicEcon
using Test

@testset "NordicEcon.jl" begin
    # Write your tests here.
end

max_stars_repo_name: ordovician/NordicEcon

max_issues_repo_name: ordovician/NordicEcon

max_forks_repo_name: ordovician/NordicEcon


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 2fe2f106-b028-4821-94e4-2ecc2b995e14
begin
    import Pkg
    # activate a clean environment
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="HTTP"),
        Pkg.PackageSpec(name="JSON3"),
        # ... keep adding your packages
    ])

    using HTTP
    using JSON3
    # ... place all usings and imports into this one cell
end

# ╔═╡ f40fe198-e3d9-4095-a3fa-a0dedf273572
j = JSON3.read(HTTP.get("https://pokeapi.co/api/v2/type/3").body)

# ╔═╡ 7b86d6e7-b359-4f25-8bc3-ff76731bdfbd
typeof(j)

# ╔═╡ 55a39434-074c-4ef6-a724-6d7cbe07283a

# ╔═╡ 99da4570-e5c0-460b-b7dd-4508797edb2b
all_type3_pokemon_names = [x["pokemon"]["name"] for x in j[:pokemon]]

# ╔═╡ fb7e576c-cb62-446b-b705-85ad10478204
j

# ╔═╡ efdb11d8-3102-4395-ac88-89fed22efdb1

# ╔═╡ Cell order:
# ╠═2fe2f106-b028-4821-94e4-2ecc2b995e14
# ╠═f40fe198-e3d9-4095-a3fa-a0dedf273572
# ╠═7b86d6e7-b359-4f25-8bc3-ff76731bdfbd
# ╠═55a39434-074c-4ef6-a724-6d7cbe07283a
# ╠═99da4570-e5c0-460b-b7dd-4508797edb2b
# ╠═fb7e576c-cb62-446b-b705-85ad10478204
# ╠═efdb11d8-3102-4395-ac88-89fed22efdb1

max_stars_repo_name: alvarmaciel/paginas-pluto

max_issues_repo_name: alvarmaciel/paginas-pluto

max_forks_repo_name: alvarmaciel/paginas-pluto


text:
module mmmtdt1
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
function test()
	MR = DeforModelRed3D
	symmet(a) = a + transpose(a)

	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6

	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6

	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
end
end
using .mmmtdt1
mmmtdt1.test()

module mmmtdtens4A1
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156
	F = [1.02 0.03 -0.04; 0.01 0.99 -0.03; -0.01 0.02 0.95]
	C = fill(0.0, 3, 3, 3, 3)
	for I in 1:3, J in 1:3, K in 1:3, L in 1:3
		C[I, J, K, L] = lambda * delta(I, J) * delta(K, L) +
		mu * (delta(I, K) * delta(J, L) + delta(I, L) * delta(J, K))
	end
	Cm = fill(0.0, 6, 6)
	tens4symmtto6x6t!(Cm, C)

	mI = Matrix(Diagonal([1.0, 1.0, 1.0, 0.5, 0.5, 0.5]))
	m1 = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0];
	D = lambda * m1 * m1' + 2. * mu * mI;

	@test norm(Cm - D) <= 1.0e-6
	return true
end
end
using .mmmtdtens4A1
mmmtdtens4A1.test()

module mmmtdtens4B1
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	delta = (I, J) -> I == J ? 1.0 : 0.0
	tens4ijkl!(t, delta, delta)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	@test norm(tr(S) * I - tens4dot2!(tS, t, S)) <= 1.0e-12

	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B1
mmmtdtens4B1.test()

module mmmtdtens4B2
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	delta = (I, J) -> I == J ? 1.0 : 0.0
	tens4iljk!(t, delta, delta)
	S = rand(3, 3)
	# @show transpose(S)
	tS = fill(0.0, 3, 3)
	# @show transpose(S) - tens4dot2!(tS, t, S)
	@test norm(transpose(S) - tens4dot2!(tS, t, S)) <= 1.0e-12

	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B2
mmmtdtens4B2.test()

module mmmtdtens4B3
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	delta = (I, J) -> I == J ? 1.0 : 0.0
	tens4ikjl!(t, delta, delta)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	# @show S - tens4dot2!(tS, t, S)
	@test norm(S - tens4dot2!(tS, t, S)) <= 1.0e-12

	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B3
mmmtdtens4B3.test()

module mmmtdtens4B4
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4identity!(t)
	@test norm(S - tens4dot2!(tS, t, S)) <= 1.0e-12
	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B4
mmmtdtens4B4.test()

module mmmtdtens4B5
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4transposor!(t)
	@test norm(S' - tens4dot2!(tS, t, S)) <= 1.0e-12
	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B5
mmmtdtens4B5.test()

module mmmtdtens4B6
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4tracor!(t)
	@test norm(tr(S) * 1I - tens4dot2!(tS, t, S)) <= 1.0e-12
	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B6
mmmtdtens4B6.test()

module mmmtdtens4B8
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	delta = (I, J) -> I == J ? 1.0 : 0.0
	lambda = 3.3
	mu = 0.156

	t = fill(0.0, 3, 3, 3, 3)
	S = rand(3, 3)
	tens4skewor!(t)
	tS = rand(3, 3)
	@test norm((S - S')/2 - tens4dot2!(tS, t, S)) <= 1.0e-12
	# @btime tens4dot2!($tS, $t, $S)
	return true
end
end
using .mmmtdtens4B8
mmmtdtens4B8.test()

module mmmtdtens4C1
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	t = fill(0.0, 3, 3, 3, 3)
	tens4identity!(t)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4dot2!(tS, t, S)
	# @show S - tS
	@test norm(S - tS) <= 1.0e-12
	return true
end
end
using .mmmtdtens4C1
mmmtdtens4C1.test()

module mmmtdtens4C2
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	t = fill(0.0, 3, 3, 3, 3)
	tens4transposor!(t)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4dot2!(tS, t, S)
	# @show S - tS
	@test norm(S' - tS) <= 1.0e-12
	return true
end
end
using .mmmtdtens4C2
mmmtdtens4C2.test()

module mmmtdtens4C3
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	t = fill(0.0, 3, 3, 3, 3)
	tens4tracor!(t)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4dot2!(tS, t, S)
	# @show tr(S) * I - tS
	@test norm(tr(S) * I - tS) <= 1.0e-12
	return true
end
end
using .mmmtdtens4C3
mmmtdtens4C3.test()

module mmmtdtens4C4
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	t = fill(0.0, 3, 3, 3, 3)
	tens4symmetrizor!(t)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4dot2!(tS, t, S)
	# @show (S + S')/2 * I - tS
	@test norm((S + S')/2 - tS) <= 1.0e-12
	return true
end
end
using .mmmtdtens4C4
mmmtdtens4C4.test()

module mmmtdtens4C5
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	t = fill(0.0, 3, 3, 3, 3)
	tens4skewor!(t)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4dot2!(tS, t, S)
	# @show (S - S')/2 * I - tS
	@test norm((S - S')/2 - tS) <= 1.0e-12
	return true
end
end
using .mmmtdtens4C5
mmmtdtens4C5.test()

module mmmtdtens4C6
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
# using BenchmarkTools

function test()
	t = fill(0.0, 3, 3, 3, 3)
	tens4deviator!(t)
	S = rand(3, 3)
	tS = fill(0.0, 3, 3)
	tens4dot2!(tS, t, S)
	# @show tr((S - tr(S)/3*I) ), tr(tS)
	@test norm(tr((S - tr(S)/3*I) ) - tr(tS)) <= 1.0e-12
	return true
end
end
using .mmmtdtens4C6
mmmtdtens4C6.test()


text:
"""
    ZeroKernel()

Create a kernel that always returning zero
```
    κ(x,y) = 0.0
```
The output type depends of `x` and `y`
"""
struct ZeroKernel <: SimpleKernel end

kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)

metric(::ZeroKernel) = Delta()

Base.show(io::IO, ::ZeroKernel) = print(io, "Zero Kernel")

"""
    WhiteKernel()

```
    κ(x,y) = δ(x,y)
```
Kernel function working as an equivalent to add white noise. Can also be called via `EyeKernel()`
"""
struct WhiteKernel <: SimpleKernel end

"""
    EyeKernel()

See [WhiteKernel](@ref)
"""
const EyeKernel = WhiteKernel

kappa(κ::WhiteKernel, δₓₓ::Real) = δₓₓ

metric(::WhiteKernel) = Delta()

Base.show(io::IO, ::WhiteKernel) = print(io, "White Kernel")

"""
    ConstantKernel(; c=1.0)

Kernel function always returning a constant value `c`
```
    κ(x,y) = c
```
"""
struct ConstantKernel{Tc<:Real} <: SimpleKernel
    c::Vector{Tc}
    function ConstantKernel(;c::T=1.0) where {T<:Real}
        new{T}([c])
    end
end

kappa(κ::ConstantKernel,x::Real) = first(κ.c)*one(x)

metric(::ConstantKernel) = Delta()

Base.show(io::IO, κ::ConstantKernel) = print(io, "Constant Kernel (c = ", first(κ.c), ")")


text:
Turn(t::Turtle, r::Angle) = Turn(t, upreferred(-r / 1°))
Orientation(t::Turtle, r::Angle) = Orientation(t, upreferred(-r / 1°))
Circle(t::Turtle, radius::Length) = Circle(t, scale_to_pt(r))
Forward(t::Turtle, d::Length) = Forward(t, scale_to_pt(d))

"""
    drawcart(p = O, α = 0°, l = 2m, w = 1m, d = 0.5m, tf = 0.2)

Using turtle graphics, draw the outline of a cart, centered on p
(default is origo) and rotated by α. The outline is symmetrical arond two axes.
"""
function drawcart(;p = O, α = 0°, l = 2m, w = 2m / 1.618034, d = 0.5m, tf = 0.2)
    t = Turtle()
    Pencolor(t, get_current_RGB())
    Reposition(t, p)
    Orientation(t, α)
    Penup(t)
    Forward(t, -l /2)
    Turn(t, -90°)
    Pendown(t)
    # Draw one side, then the next.
    for i = 1:2
        Forward(t, w / 2)
        # Turn back to axle
        Turn(t, 90°)
        Forward(t, d / 2)
        # Axle
        Turn(t, -90°)
        Forward(t, d * tf)
        Turn(t, -90°)
        #Wheel
        Forward(t, d /2)
        Turn(t, 90°)
        Forward(t, d * tf)
        Turn(t, 90°)
        Forward(t, d )
        Turn(t, 90°)
        Forward(t, d * tf )
        Turn(t, 90°)
        Forward(t, d / 2)
        # Axle again
        Turn(t, -90°)
        Forward(t, d * tf)
        Turn(t, -90°)
        # Side
        Forward(t, l / 2 - d / 2)
        # Now we should have been able to mirror, but we can't.
        Forward(t, l / 2 - d / 2)
        # Axle
        Turn(t, -90°)
        Forward(t, d * tf)
        Turn(t, -90°)
        #Wheel
        Forward(t, d /2)
        Turn(t, 90°)
        Forward(t, d * tf)
        Turn(t, 90°)
        Forward(t, d )
        Turn(t, 90°)
        Forward(t, d * tf )
        Turn(t, 90°)
        Forward(t, d / 2)
        # Axle again
        Turn(t, -90°)
        Forward(t, d * tf)
        Turn(t, -90°)
        # Side to rear
        Forward(t, d / 2)
        Turn(t, 90°)
        Forward(t, w / 2)
    end
end


text:
POMDPs.solve(solver::DPWSolver, mdp::Union{POMDP,MDP}) = DPWPlanner(solver, mdp)

"""
Delete existing decision tree.
"""
function clear_tree!(p::DPWPlanner)
    p.tree = nothing
end

"""
Construct an MCTSDPW tree and choose the best action.
"""
POMDPs.action(p::DPWPlanner, s) = first(action_info(p, s))

"""
Construct an MCTSDPW tree and choose the best action. Also output some information.
"""
function POMDPModelTools.action_info(p::DPWPlanner, s; tree_in_info=false)
    local a::actiontype(p.mdp)
    info = Dict{Symbol, Any}()
    try
        if isterminal(p.mdp, s)
            error("""
                  MCTS cannot handle terminal states. action was called with
                  s = $s
                  """)
        end

        S = statetype(p.mdp)
        A = actiontype(p.mdp)
        if p.solver.keep_tree
            if p.tree == nothing
                tree = DPWTree{S,A}(p.solver.n_iterations)
                p.tree = tree
            else
                tree = p.tree
            end
            if haskey(tree.s_lookup, s)
                snode = tree.s_lookup[s]
            else
                snode = insert_state_node!(tree, s, true)
            end
        else
            tree = DPWTree{S,A}(p.solver.n_iterations)
            p.tree = tree
            snode = insert_state_node!(tree, s, p.solver.check_repeat_state)
        end

        nquery = 0
        start_us = CPUtime_us()
        for i = 1:p.solver.n_iterations
            nquery += 1
            simulate(p, snode, p.solver.depth) # (not 100% sure we need to make a copy of the state here)
            if CPUtime_us() - start_us >= p.solver.max_time * 1e6
                break
            end
        end
        info[:search_time_us] = CPUtime_us() - start_us
        info[:tree_queries] = nquery
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end

        best_Q = -Inf
        sanode = 0
        for child in tree.children[snode]
            if tree.q[child] > best_Q
                best_Q = tree.q[child]
                sanode = child
            end
        end
        # XXX some publications say to choose action that has been visited the most
        a = tree.a_labels[sanode] # choose action with highest approximate value
    catch ex
        a = convert(actiontype(p.mdp), default_action(p.solver.default_action, p.mdp, s, ex))
        info[:exception] = ex
    end

    return a, info
end

"""
Return the reward for one iteration of MCTSDPW.
"""
function simulate(dpw::DPWPlanner, snode::Int, d::Int)
    S = statetype(dpw.mdp)
    A = actiontype(dpw.mdp)
    sol = dpw.solver
    tree = dpw.tree
    s = tree.s_labels[snode]
    if d == 0 || isterminal(dpw.mdp, s)
        return 0.0
    end

    # action progressive widening
    if dpw.solver.enable_action_pw
        if length(tree.children[snode]) <= sol.k_action*tree.total_n[snode]^sol.alpha_action # criterion for new action generation
            a = next_action(dpw.next_action, dpw.mdp, s, DPWStateNode(tree, snode)) # action generation step
            if !sol.check_repeat_action || !haskey(tree.a_lookup, (snode, a))
                n0 = init_N(sol.init_N, dpw.mdp, s, a)
                insert_action_node!(tree, snode, a, n0,
                                    init_Q(sol.init_Q, dpw.mdp, s, a),
                                    sol.check_repeat_action
                                   )
                tree.total_n[snode] += n0
            end
        end
    elseif isempty(tree.children[snode])
        for a in actions(dpw.mdp, s)
            n0 = init_N(sol.init_N, dpw.mdp, s, a)
            insert_action_node!(tree, snode, a, n0,
                                init_Q(sol.init_Q, dpw.mdp, s, a),
                                false)
            tree.total_n[snode] += n0
        end
    end

    best_UCB = -Inf
    sanode = 0
    ltn = log(tree.total_n[snode])
    for child in tree.children[snode]
        n = tree.n[child]
        q = tree.q[child]
        c = sol.exploration_constant # for clarity
        if (ltn <= 0 && n == 0) || c == 0.0
            UCB = q
        else
            UCB = q + c*sqrt(ltn/n)
        end
        @assert !isnan(UCB) "UCB was NaN (q=$q, c=$c, ltn=$ltn, n=$n)"
        @assert !isequal(UCB, -Inf)
        if UCB > best_UCB
            best_UCB = UCB
            sanode = child
        end
    end

    a = tree.a_labels[sanode]

    # state progressive widening
    new_node = false
    if tree.n_a_children[sanode] <= sol.k_state*tree.n[sanode]^sol.alpha_state
        sp, r = generate_sr(dpw.mdp, s, a, dpw.rng)

        if sol.check_repeat_state && haskey(tree.s_lookup, sp)
            spnode = tree.s_lookup[sp]
        else
            spnode = insert_state_node!(tree, sp, sol.keep_tree || sol.check_repeat_state)
            new_node = true
        end

        push!(tree.transitions[sanode], (spnode, r))

        if !sol.check_repeat_state
            tree.n_a_children[sanode] += 1
        elseif !((sanode,spnode) in tree.unique_transitions)
            push!(tree.unique_transitions, (sanode,spnode))
            tree.n_a_children[sanode] += 1
        end
    else
        spnode, r = rand(dpw.rng, tree.transitions[sanode])
        sp = tree.s_labels[spnode]
    end

    # THE NEXT TWO LINES WERE ADDED TO ALLOW FOR TERMINAL ACTIONS
    if isterminal(dpw.mdp, sp, a)
        q = r
    elseif new_node
        q = r + discount(dpw.mdp)*estimate_value(dpw.solved_estimate, dpw.mdp, sp, d-1)
    else
        q = r + discount(dpw.mdp)*simulate(dpw, spnode, d-1)
    end

    tree.n[sanode] += 1
    tree.total_n[snode] += 1

    tree.q[sanode] += (q - tree.q[sanode])/tree.n[sanode]

    return q
end


text:
# okay, this is for my autodiff package which I'm going to do briefly just to try to understand
#how this actually works out, as it's my next julia project. It shuold be very very educatoinal to be honest as it'll teach us more abuot the type system and the maths behind autodiff as well as using macros and a bit about compilation if we are going to get it to work as we should

# the aims are as follows: we implement autodiff types and the autodiff dual number functions for all the common mathemtical operations
# we implement a macro which can then take a function written with thos standard and automagically compile it into an autodiffed one
# should be fairl straightforward, but requires a deep understanding. now let's try to gain that understanding

# first our dual number type

import Base: *

type Dual{T} <: Real
	r::T # the real part
	d::T # the derivative part

	#function Dual(r::T, d::T)
	#	return new(r,d)
	#end

	#function Dual(r::T)
	#	return new(r,0)	#this sets the default derivative as zero, which seems reasonable atm
	#end
end

#some simple helper functoins
function getReal(num::Dual)
	return num.r
end

function getDerivative(num::Dual)
	return num.d
end

# okay, now we're going to implement finite differences so I can check my answers

function finite_differences(f::Function, start, step)
	return (f(start+step) - f(step))/step
end

# okay, let's implement our AD functions using overloaded operators

#plus function
function +(u::Dual, v::Dual)
	return Dual((u.r + v.r), (u.d+v.d))
end

#multiply function
function *(u::Dual, v::Dual)
	return Dual((u.r*v.r), (u.d *v.r + u.r * v.d))
end

function -(u::Dual, v::Dual)
	return Dual((u.r-v.r), (u.d-r.d))
end

function negate(u::Dual)
	return Dual(-u.r, -u.d)
end

function signum(u::Dual)
	return Dual(signum(u.r),0)
end

function abs(u::Dual)
	return dual(abs(u.r), u.d*(signum(u.r)))
end

function /(u::Dual, v::Dual)
	return Dual((u.r/v.r), ((u.d*v.r - u.r*v.d)/(v.r^2)))
end

function exp(u::Dual)
	return Dual(exp(u.r), (u.d*exp(u.r)))
end

function sqrt(u::Dual)
	return Dual(sqrt(u.r), (u.d/(2*sqrt(u.r))))
end

function log(u::Dual)
	return Dual(log(u.r), (u.r/u.d))
end

function sin(u::Dual)
	return Dual(sin(u.r), (u.d*cos(u.r)))
end

function cos(u::Dual)
	return Dual(cos(u.r), (-u.d*sin(u.r)))
end

function tan(u::Dual)
	return Dual(tan(u.r), (1/(cos(u.r)^2)))
end

function asin(u::Dual)
	return Dual(asin(u.r), (u.d/(sqrt(1-u.r)^2)))
end

function acos(u::Dual)
	return Dual(acos(u.r), (-u.d/sqrt(1-u.r)^2))
end

function atan(u::Dual)
	return Dual(atan(u.r), (u.d/(1+u.r)^2))
end

function sinh(u::Dual)
	return Dual(sinh(u.r), (u.d*cosh(u.r)))
end

function cosh(u::Dual)
	return Dual(cosh(u.r), (u.d*sinh(u.r)))
end

function tanh(u::Dual)
	return Dual(tanh(u.r), (u.d * (1- (tanh(u.r))^2)))
end

function asinh(u::Dual)
	return Dual(asinh(u.r), (u.d/(sqrt(1 + u.r^2))))
end

function acosh(u::Dual)
	return Dual(acosh(u.r), (u.d/sqrt(u.r^2 -1)))
end

function atanh(u::Dual)
	return Dual(atanh(u.r), (u.d/(1-u.r^2)))
end

# okay, let's test this

function f(a,b)
	return a + b
end

a = Dual(4,1)
b = Dual(6,2)

f(a,b)

max_stars_repo_name: Bmillidgework/Misc-Maths

max_issues_repo_name: Bmillidgework/Misc-Maths

max_forks_repo_name: Bmillidgework/Misc-Maths


text:
"""
    drop_input!(state, agent_idx, config)

Drop input for specific agent (~ unfollow on Twitter).

# Arguments
- `state`: A tuple of the current graph and agent_list
- `agent_idx`: Agent index
- `config`: Config object as provided by Config

See also: [`Config`](@ref), [`add_input!`](@ref)
"""
function drop_input!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config
)

    graph, agent_list = state
    this_agent = agent_list[agent_idx]

    unfollow_candidates = Array{Tuple{Int64, Int64}, 1}()
    for post in this_agent.feed
        if (
            abs(post.opinion - this_agent.opinion)
            > config.opinion_threshs.unfollow
        )
            if (
                abs(agent_list[post.source_agent].opinion - this_agent.opinion)
                > config.opinion_threshs.unfollow
            )
                if (
                    (
                        outdegree(graph, post.source_agent)
                        / outdegree(graph, agent_idx)
                    ) > 1
                    && rand() > 0.5
                )
                    push!(
                        unfollow_candidates,
                        (post.source_agent, outdegree(graph, post.source_agent))
                    )
                elseif (
                    (
                        outdegree(graph, post.source_agent)
                        / outdegree(graph, agent_idx)
                    ) <= 1
                )
                    push!(
                        unfollow_candidates,
                        (post.source_agent, outdegree(graph, post.source_agent))
                    )
                end
            end
        end
    end

    sort!(unfollow_candidates, by=last)

    drop_count = min(
        length(unfollow_candidates),
        ceil(
            Int64,
            indegree(graph,agent_idx) * config.agent_props.unfollow_rate
        )
    )

    for i in 1:drop_count
        rem_edge!(graph, unfollow_candidates[i][1],agent_idx)
    end

    return state

end

# suppress output of include()
;


text:
module QuadraticTests

using Interpolations, Base.Test

for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
    f(x) = sin((x-3)*2pi/9 - 1)
    xmax = 10
    A = Float64[f(x) for x in 1:xmax]
    for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
        itp1 = @inferred(constructor(copier(A), BSpline(Quadratic(BC())), GT()))
        @test @inferred(size(itp1)) == size(A)

        # test that inner region is close to data
        for x in 3.1:.2:8.1
            @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
        end

        # test that we can evaluate close to, and at, boundaries
        if GT == OnGrid
            itp1[1.]
            itp1[1.0]
            itp1[1.2]
            itp1[9.8]
            itp1[10.]
            itp1[10]
        else
            itp1[0.5]
            itp1[0.6]
            itp1[10.4]
            itp1[10.5]
        end
    end

    f(x,y) = sin(x/10)*cos(y/6)
    xmax, ymax = 30,10
    A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax]

    # test that inner region is close to data
    for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
        itp2 = @inferred(constructor(copier(A), BSpline(Quadratic(BC())), GT()))
        @test @inferred(size(itp2)) == size(A)

        for x in 3.1:.2:xmax-3, y in 3.1:2:ymax-3
            @test_approx_eq_eps f(x,y) itp2[x,y] abs(.1*f(x,y))
        end
    end
end

let
    f(x) = sin((x-3)*2pi/9 - 1)
    xmax = 10
    A = Float64[f(x) for x in 1:xmax]
    itp1 = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell())
    for i = 1:xmax
        @test_approx_eq itp1[i] A[i]
    end

    f(x,y) = sin(x/10)*cos(y/6)
    xmax, ymax = 30,10
    A = Float64[f(x,y) for x in 1:xmax, y in 1:ymax]
    itp2 = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell())
    for j = 1:ymax, i = 1:xmax
        @test_approx_eq itp2[i,j] A[i,j]
    end
end

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function reopt_results(m::JuMP.AbstractModel, p::REoptInputs; _n="")
	tstart = time()
    d = Dict{String, Any}()
    for b in p.storage.types
        if p.storage.max_kw[b] > 0
            add_storage_results(m, p, d, b; _n)
        end
    end

    add_electric_tariff_results(m, p, d; _n)
    add_electric_utility_results(m, p, d; _n)
    add_financial_results(m, p, d; _n)

	if !isempty(p.pvtechs)
        add_pv_results(m, p, d; _n)
	end

    if "Wind" in p.techs
        add_wind_results(m, p, d; _n)
    end

	time_elapsed = time() - tstart
	@info "Base results processing took $(round(time_elapsed, digits=3)) seconds."

	if !isempty(p.gentechs) && isempty(_n)  # generators not included in multinode model
        tstart = time()
		add_generator_results(m, p, d)
        time_elapsed = time() - tstart
        @info "Generator results processing took $(round(time_elapsed, digits=3)) seconds."
	end

	if !isempty(p.elecutil.outage_durations) && isempty(_n)  # outages not included in multinode model
        tstart = time()
		add_outage_results(m, p, d)
        time_elapsed = time() - tstart
        @info "Outage results processing took $(round(time_elapsed, digits=3)) seconds."
	end
	return d
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


text:
#################### Adaptive Mixture Metropolis ####################

#################### Types and Constructors ####################

type AMMTune <: SamplerTune
  logf::Nullable{Function}
  adapt::Bool
  beta::Float64
  m::Int
  Mv::Vector{Float64}
  Mvv::Matrix{Float64}
  scale::Float64
  SigmaL::LowerTriangular{Float64}
  SigmaLm::Matrix{Float64}

  AMMTune() = new()

  function AMMTune{T<:Real}(x::Vector, Sigma::Matrix{T},
                            logf::Nullable{Function}; beta::Real=0.05,
                            scale::Real=2.38)
    new(logf, false, beta, 0, Vector{Float64}(), Matrix{Float64}(), scale,
        cholfact(Sigma)[:L], Matrix{Float64}())
  end
end

AMMTune{T<:Real}(x::Vector, Sigma::Matrix{T}; args...) =
  AMMTune(x, Sigma, Nullable{Function}(); args...)

AMMTune{T<:Real}(x::Vector, Sigma::Matrix{T}, logf::Function; args...) =
  AMMTune(x, Sigma, Nullable{Function}(logf); args...)

typealias AMMVariate SamplerVariate{AMMTune}

function validate(v::AMMVariate)
  n = length(v)
  size(v.tune.SigmaL, 1) == n ||
    throw(ArgumentError("Sigma dimension differs from variate length $n"))
  v
end

#################### Sampler Constructor ####################

function AMM{T<:Real}(params::ElementOrVector{Symbol}, Sigma::Matrix{T};
                      adapt::Symbol=:all, args...)
  adapt in [:all, :burnin, :none] ||
    throw(ArgumentError("adapt must be one of :all, :burnin, or :none"))

  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, true)
    v = SamplerVariate(block, Sigma; args...)
    isadapt = adapt == :burnin ? model.iter <= model.burnin :
              adapt == :all ? true : false
    sample!(v, x -> logpdf!(block, x), adapt=isadapt)
    relist(block, v)
  end
  Sampler(params, samplerfx, AMMTune())
end

#################### Sampling Functions ####################

sample!(v::AMMVariate; args...) = sample!(v, v.tune.logf; args...)

function sample!(v::AMMVariate, logf::Function; adapt::Bool=true)
  tune = v.tune
  n = length(v)

  setadapt!(v, adapt)

  x = tune.SigmaL * randn(n)
  if tune.m > 2 * n
    x = tune.beta * x + (1.0 - tune.beta) * (tune.SigmaLm * randn(n))
  end
  x += v
  if rand() < exp(logf(x) - logf(v.value))
    v[:] = x
  end

  if tune.adapt
    tune.m += 1
    p = tune.m / (tune.m + 1.0)
    tune.Mv = p * tune.Mv + (1.0 - p) * v
    tune.Mvv = p * tune.Mvv + (1.0 - p) * v * v'
    Sigma = (tune.scale^2 / n / p) * (tune.Mvv - tune.Mv * tune.Mv')
    F = cholfact(Sigma, :U, Val{true})
    if rank(F) == n
      tune.SigmaLm = F[:P] * F[:L]
    end
  end

  v
end

function setadapt!(v::AMMVariate, adapt::Bool)
  tune = v.tune
  if adapt && !tune.adapt
    n = length(v)
    tune.m = 0
    tune.Mv = v
    tune.Mvv = v * v'
    tune.SigmaLm = zeros(n, n)
  end
  tune.adapt = adapt
  v
end


text:
## Exercise 7-2
## Copy the loop from Square Roots and encapsulate it in a function called mysqrt that takes a as a parameter, chooses a reasonable value of x, and returns an estimate of the square root of a.

function mysqrt(a, x)
    a = float(a)
    while true
        y = (x + a/x) / 2
        if y == x
            break
        end
        x = y
    end
    return x
end

## To test it, write a function named testsquareroot that prints a table like this:
# see ex 7-2 in book
println("Ans: ")

function rightjustify(string::String, chr)
    if chr - length(string) > 0
        return ' ' ^ (chr - length(string)) * string
    end
    return string
end

function testsquareroot()
    println("a   mysqrt             sqrt               diff")
    println("-   ------             ----               ----")
    for i in 1:10
        i = float(i)
        root = sqrt(i)
        root_my = mysqrt(i, i/2)
        println(
            string(i) * " " *
            rightjustify(string(root_my), 18) *
            " " *
            rightjustify(string(root), 18) * " " * string(abs(root_my - root)))
    end
end

testsquareroot()

println("End.")


text:
import MPI

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
println("Hi form $(rank)")

MPI.Finalize()

max_stars_repo_name: MiroK/seminar-julia

max_issues_repo_name: MiroK/seminar-julia

max_forks_repo_name: MiroK/seminar-julia


text:
export laplacian_op, poisson_op, sat_op, upwlap_op, upwps_op, fwi_op, fwi_obs_op, sat_op2, eikonal

OPS_DIR = joinpath(@__DIR__, "../deps/CustomOps")

@doc raw"""
    fwi_op(lambda::Union{PyObject, Array{Float64}},mu::Union{PyObject, Array{Float64}},
        den::Union{PyObject, Array{Float64}},stf::Union{PyObject, Array{Float64}},
        gpu_id::Union{PyObject, Integer},shot_ids::Union{PyObject, Array{T}},para_fname::String) where T<:Integer

Computes the FWI loss function.

- `lambda` : Lame's first parameter (unit: MPa)
- `mu` : Lame's second parameter (shear modulus, unit: MPa)
- `den` : Density
- `stf` : Source time functions
- `gpu_id` : The ID of GPU to run this FWI operator
- `shot_ids` : The source function IDs (determining the location of sources)
- `para_fname` : Parameter file location
"""
function fwi_op(lambda::Union{PyObject, Array{Float64}},mu::Union{PyObject, Array{Float64}},
        den::Union{PyObject, Array{Float64}},stf::Union{PyObject, Array{Float64}},
        gpu_id::Union{PyObject, Integer},shot_ids::Union{PyObject, Array{T}},para_fname::String) where T<:Integer
    lambda = convert_to_tensor(lambda, dtype=Float64)
    mu = convert_to_tensor(mu, dtype=Float64)
    den = convert_to_tensor(den, dtype=Float64)
    stf = convert_to_tensor(stf, dtype=Float64)
    gpu_id = convert_to_tensor(gpu_id, dtype=Int32)
    shot_ids = convert_to_tensor(shot_ids, dtype=Int32)
    fwi_op = load_op_and_grad("$OPS_DIR/FWI/build/libFwiOp", "fwi_op")
    fwi_op(lambda,mu,den,stf,gpu_id,shot_ids,para_fname)
end

@doc raw"""
    fwi_obs_op(lambda::Union{PyObject, Array{Float64}},mu::Union{PyObject, Array{Float64}},
        den::Union{PyObject, Array{Float64}},stf::Union{PyObject, Array{Float64}},
        gpu_id::Union{PyObject, Integer},shot_ids::Union{PyObject, Array{T}},para_fname::String) where T<:Integer

Generates the observation data and store them as files which will be used by [`fwi_op`](@ref)
For the meaning of parameters, see [`fwi_op`](@ref).
"""
function fwi_obs_op(lambda::Union{PyObject, Array{Float64}},mu::Union{PyObject, Array{Float64}},
    den::Union{PyObject, Array{Float64}},stf::Union{PyObject, Array{Float64}},
    gpu_id::Union{PyObject, Integer},shot_ids::Union{PyObject, Array{T}},para_fname::String) where T<:Integer
    lambda = convert_to_tensor(lambda, dtype=Float64)
    mu = convert_to_tensor(mu, dtype=Float64)
    den = convert_to_tensor(den, dtype=Float64)
    stf = convert_to_tensor(stf, dtype=Float64)
    gpu_id = convert_to_tensor(gpu_id, dtype=Int32)
    shot_ids = convert_to_tensor(shot_ids, dtype=Int32)
    fwi_obs_op = load_op("$OPS_DIR/FWI/build/libFwiOp", "fwi_obs_op")
    fwi_obs_op(lambda, mu, den, stf, gpu_id, shot_ids, para_fname)
end

@doc raw"""
    laplacian_op(coef::Union{PyObject, Array{Float64}}, f::Union{PyObject, Array{Float64}},
            h::Union{PyObject, Float64}, ρ::Union{PyObject, Float64})

Computes the Laplacian of function $f(\mathbf{x})$; here ($\mathbf{x}=[z\quad x]^T$)
```math
-\nabla\cdot\left(c(\mathbf{x}) \nabla \left(u(\mathbf{x}) -\rho \begin{bmatrix}z \\ 0\end{bmatrix}  \right)\right)
```
"""
function laplacian_op(coef::Union{PyObject, Array{Float64}}, f::Union{PyObject, Array{Float64}},
    h::Union{PyObject, Float64}, ρ::Union{PyObject, Float64})
    coef = convert_to_tensor(coef, dtype=Float64)
    f = convert_to_tensor(f, dtype=Float64)
    h = convert_to_tensor(h, dtype=Float64)
    ρ = convert_to_tensor(ρ, dtype=Float64)
    laplacian = load_op_and_grad("$OPS_DIR/Laplacian/build/libLaplacian", "laplacian")
    laplacian_op(coef, f, h, ρ)
end

@doc raw"""
    poisson_op(c::Union{PyObject, Float64}, g::Union{PyObject, Float64},
        h::Union{PyObject, Float64}, ρ::Union{PyObject, Float64}, index::Union{Integer, PyObject}=0)

Solves the Poisson equation ($\mathbf{x}=[z\quad x]^T$)

$\begin{aligned}
-\nabla\cdot\left(c(\mathbf{x}) \nabla \left(u(\mathbf{x}) -\rho  \begin{bmatrix}z \\ 0\end{bmatrix}   \right)\right) &=  g(\mathbf{x}) & \mathbf{x}\in \Omega\\
\frac{\partial u(x)}{\partial n} &=  0 & \mathbf{x}\in \Omega\\
\end{aligned}$

Here $\Omega=[0,n_zh]\times [0, n_xh]$. The equation is solved using finite difference method, where the step size in each direction is $h$. Mathematically, the solution to the PDE is determined up to a constant. Numerically, we discretize the equation with the scheme
```math
(A+E_{11})\mathbf{u} = \mathbf{f}
```
where $A$ is the finite difference coefficient matrix,

```math
(E_{11})_{ij} = \left\{ \begin{matrix}1 & i=j=1 \\ 0 & \text{ otherwise }\end{matrix}\right.
```

- `index` : `Int32`, when `index=1`, `SparseLU` is used to solve the linear system; otherwise the function invokes algebraic multigrid method from `amgcl`.
"""
function poisson_op(c::Union{PyObject, Array{Float64}}, g::Union{PyObject, Array{Float64}},
    h::Union{PyObject, Float64},
    ρ::Union{PyObject, Float64}, index::Union{PyObject, Integer}=0)
    c = convert_to_tensor(c, dtype=Float64)
    g = convert_to_tensor(g, dtype=Float64)
    h = convert_to_tensor(h, dtype=Float64)
    ρ = convert_to_tensor(ρ, dtype=Float64)
    index = convert_to_tensor(index, dtype=Int64)
    poisson_op = load_op_and_grad("$OPS_DIR/Poisson/build/libPoissonOp", "poisson_op")
    poisson_op(c, g, h, ρ, index)
end

@doc raw"""
    sat_op(s0::Union{PyObject, Array{Float64}},pt::Union{PyObject, Array{Float64}},
        permi::Union{PyObject, Array{Float64}},poro::Union{PyObject, Array{Float64}},
        qw::Union{PyObject, Array{Float64}},qo::Union{PyObject, Array{Float64}},
        muw::Union{PyObject, Float64},muo::Union{PyObject, Float64},
        sref::Union{PyObject, Array{Float64}},dt::Union{PyObject, Float64},h::Union{PyObject, Float64})

Solves the following discretized equation
```math
\phi (S_2^{n + 1} - S_2^n) - \nabla  \cdot \left( {{m_2}(S_2^{n + 1})K\nabla \Psi _2^n} \right) \Delta t= \left(q_2^n + q_1^n \frac{m_2(S^{n+1}_2)}{m_1(S^{n+1}_2)}\right) \Delta t
```
where
```math
m_2(s) = \frac{s^2}{\mu_w}\qquad m_1(s) = \frac{(1-s)^2}{\mu_o}
```
This is a nonlinear equation and is solved with the Newton-Raphson method.

- `s0` : $n_z\times n_x$, saturation of fluid, i.e., $S_2^n$
- `pt` : $n_z\times n_x$, potential of fluid, i.e., $\Psi_2^n$
- `permi` : $n_z\times n_x$, permeability, i.e., $K$
- `poro` : $n_z\times n_x$, porosity, i.e., $\phi$
- `qw` : $n_z\times n_x$, injection or production rate of ﬂuid 1, $q_2^n$
- `qo` : $n_z\times n_x$, injection or production rate of ﬂuid 2, $q_1^n$
- `muw` : viscosity of fluid 1, i.e., $\mu_w$
- `muo` : viscosity of fluid 2, i.e., $\mu_o$
- `sref` : $n_z\times n_x$, initial guess for $S_2^{n+1}$
- `dt` : Time step size
- `h` : Spatial step size
"""
function sat_op(s0::Union{PyObject, Array{Float64}},pt::Union{PyObject, Array{Float64}},
    permi::Union{PyObject, Array{Float64}},poro::Union{PyObject, Array{Float64}},
    qw::Union{PyObject, Array{Float64}},qo::Union{PyObject, Array{Float64}},
    muw::Union{PyObject, Float64},muo::Union{PyObject, Float64},
    sref::Union{PyObject, Array{Float64}},dt::Union{PyObject, Float64},h::Union{PyObject, Float64})
    s0 = convert_to_tensor(s0, dtype=Float64)
    pt = convert_to_tensor(pt, dtype=Float64)
    permi = convert_to_tensor(permi, dtype=Float64)
    poro = convert_to_tensor(poro, dtype=Float64)
    qw = convert_to_tensor(qw, dtype=Float64)
    qo = convert_to_tensor(qo, dtype=Float64)
    muw = convert_to_tensor(muw, dtype=Float64)
    muo = convert_to_tensor(muo, dtype=Float64)
    sref = convert_to_tensor(sref, dtype=Float64)
    dt = convert_to_tensor(dt, dtype=Float64)
    h = convert_to_tensor(h, dtype=Float64)
    sat_op = load_op_and_grad("$OPS_DIR/Saturation/build/libSatOp", "sat_op")
    sat_op(s0,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
end

@doc raw"""
    sat_op2(s0::Union{PyObject, Array{Float64}},
    dporodt::Union{PyObject, Array{Float64}},
    pt::Union{PyObject, Array{Float64}},
    permi::Union{PyObject, Array{Float64}},
    poro::Union{PyObject, Array{Float64}},
    qw::Union{PyObject, Array{Float64}},
    qo::Union{PyObject, Array{Float64}},
    muw::Union{PyObject, Float64},
    muo::Union{PyObject, Float64},
    sref::Union{PyObject, Array{Float64}},
    dt::Union{PyObject, Float64},
    h::Union{PyObject, Float64})

Solves the following discretized equation
```math
\phi (S_2^{n + 1} - S_2^n) + \Delta t \dot \phi S_2^{n+1} - \nabla  \cdot \left( {{m_2}(S_2^{n + 1})K\nabla \Psi _2^n} \right) \Delta t= \left(q_2^n + q_1^n \frac{m_2(S^{n+1}_2)}{m_1(S^{n+1}_2)}\right) \Delta t
```
where
```math
m_2(s) = \frac{s^2}{\mu_w}\qquad m_1(s) = \frac{(1-s)^2}{\mu_o}
```
This is a nonlinear equation and is solved with the Newton-Raphson method.

- `s0` : $n_z\times n_x$, saturation of fluid, i.e., $S_2^n$
- `dporodt` : $n_z\times n_x$, rate of porosity, $\dot \phi$
- `pt` : $n_z\times n_x$, potential of fluid, i.e., $\Psi_2^n$
- `permi` : $n_z\times n_x$, permeability, i.e., $K$
- `poro` : $n_z\times n_x$, porosity, i.e., $\phi$
- `qw` : $n_z\times n_x$, injection or production rate of ﬂuid 1, $q_2^n$
- `qo` : $n_z\times n_x$, injection or production rate of ﬂuid 2, $q_1^n$
- `muw` : viscosity of fluid 1, i.e., $\mu_w$
- `muo` : viscosity of fluid 2, i.e., $\mu_o$
- `sref` : $n_z\times n_x$, initial guess for $S_2^{n+1}$
- `dt` : Time step size
- `h` : Spatial step size
"""
function sat_op2(s0::Union{PyObject, Array{Float64}},
    dporodt::Union{PyObject, Array{Float64}},
    pt::Union{PyObject, Array{Float64}},
    permi::Union{PyObject, Array{Float64}},
    poro::Union{PyObject, Array{Float64}},
    qw::Union{PyObject, Array{Float64}},
    qo::Union{PyObject, Array{Float64}},
    muw::Union{PyObject, Float64},
    muo::Union{PyObject, Float64},
    sref::Union{PyObject, Array{Float64}},
    dt::Union{PyObject, Float64},
    h::Union{PyObject, Float64})
    s0 = convert_to_tensor(s0, dtype=Float64)
    dporodt = convert_to_tensor(dporodt, dtype=Float64)
    pt = convert_to_tensor(pt, dtype=Float64)
    permi = convert_to_tensor(permi, dtype=Float64)
    poro = convert_to_tensor(poro, dtype=Float64)
    qw = convert_to_tensor(qw, dtype=Float64)
    qo = convert_to_tensor(qo, dtype=Float64)
    muw = convert_to_tensor(muw, dtype=Float64)
    muo = convert_to_tensor(muo, dtype=Float64)
    sref = convert_to_tensor(sref, dtype=Float64)
    dt = convert_to_tensor(dt, dtype=Float64)
    h = convert_to_tensor(h, dtype=Float64)
    saturation_ = load_op_and_grad("$OPS_DIR/Saturation2/build/libSaturation","saturation")
    saturation_(s0,dporodt,pt,permi,poro,qw,qo,muw,muo,sref,dt,h)
end

@doc raw"""
    upwlap_op(perm::Union{PyObject, Array{Float64}},
        mobi::Union{PyObject, Array{Float64}},
        func::Union{PyObject, Array{Float64}},
        h::Union{PyObject, Float64},
        rhograv::Union{PyObject, Float64})

Computes the Laplacian of function $f(\mathbf{x})$; here $\mathbf{x}=[z\quad x]^T$.
```math
\nabla\cdot\left(m(\mathbf{x})K(\mathbf{x}) \nabla \left(f(\mathbf{x}) -\rho \begin{bmatrix}z \\ 0\end{bmatrix}  \right)\right)
```
The permeability on the computational grid is computed with Harmonic mean;
the mobility is computed with upwind scheme.

- `perm` : $n_z\times n_x$, permeability of fluid, i.e., $K$
- `mobi` : $n_z\times n_x$, mobility of fluid, i.e., $m$
- `func` : $n_z\times n_x$, potential of fluid, i.e., $f$
- `h` : `Float64`, spatial step size
- `rhograv` : `Float64`, i.e., $\rho$
"""
function upwlap_op(perm::Union{PyObject, Array{Float64}},
    mobi::Union{PyObject, Array{Float64}},
    func::Union{PyObject, Array{Float64}},
    h::Union{PyObject, Float64},rhograv::Union{PyObject, Float64})
    perm = convert_to_tensor(perm, dtype=Float64)
    mobi = convert_to_tensor(mobi, dtype=Float64)
    func = convert_to_tensor(func, dtype=Float64)
    h = convert_to_tensor(h, dtype=Float64)
    rhograv = convert_to_tensor(rhograv, dtype=Float64)
    upwlap_op = load_op_and_grad("$OPS_DIR/Upwlap/build/libUpwlapOp", "upwlap_op")
    upwlap_op(perm,mobi,func,h,rhograv)
end

@doc raw"""
    upwps_op(perm::Union{PyObject, Array{Float64}},mobi::Union{PyObject, Array{Float64}},
    src,funcref,h::Union{PyObject, Float64},rhograv::Union{PyObject, Float64},index::Union{PyObject, Integer})

Solves the Poisson equation
```math
-\nabla\cdot\left(m(\mathbf{x})K(\mathbf{x}) \nabla \left(u(\mathbf{x}) -\rho  \begin{bmatrix}z \\ 0\end{bmatrix}   \right)\right) =  g(\mathbf{x})
```
See [`upwps_op`](@ref) for detailed description.

- `perm` : $n_z\times n_x$, permeability of fluid, i.e., $K$
- `mobi` : $n_z\times n_x$, mobility of fluid, i.e., $m$
- `src` : $n_z\times n_x$, source function, i.e., $g(\mathbf{x})$
- `funcref` : $n_z\times n_x$, currently it is not not used
- `h` : `Float64`, spatial step size
- `rhograv` : `Float64`, i.e., $\rho$
- `index` : `Int32`, when `index=1`, `SparseLU` is used to solve the linear system; otherwise the function invokes algebraic multigrid method from `amgcl`.
"""
function upwps_op(perm::Union{PyObject, Array{Float64}},mobi::Union{PyObject, Array{Float64}},
    src,funcref,h::Union{PyObject, Float64},rhograv::Union{PyObject, Float64},index::Union{PyObject, Integer})
    perm = convert_to_tensor(perm, dtype=Float64)
    mobi = convert_to_tensor(mobi, dtype=Float64)
    funcref = convert_to_tensor(funcref, dtype=Float64)
    h = convert_to_tensor(h, dtype=Float64)
    rhograv = convert_to_tensor(rhograv, dtype=Float64)
    index = convert_to_tensor(index, dtype=Int64)
    upwps_op = load_op_and_grad("$OPS_DIR/Upwps/build/libUpwpsOp", "upwps_op")
    upwps_op(perm,mobi,src,funcref,h,rhograv,index)
end

export time_fractional_op, time_fractional_t_op

@doc raw"""
    time_fractional_t_op(i::Union{Integer, PyObject}, ta::PyObject,
        α::Union{Float64, PyObject}, Δt::Union{Float64, PyObject})

Returns the coefficients for the time fractional derivative
```math
{}_0^CD_t^\alpha f(t) = \frac{1}{\Gamma(1-\alpha)}\int_0^t \frac{f'(\tau)d\tau}{(t-\tau)^\alpha}
```

The discretization scheme used here is

$\begin{aligned}
& {}_0^CD_\tau^\alpha u(\tau_n) \\
= &\frac{\Delta \tau^{-\alpha}}{\Gamma(2-\alpha)}\left[G_0 u_n - \sum_{k=1}^{n-1}(G_{n-k-1}-G_{n-k})u_k + G_n u_0 \right] + \mathcal{O}(\Delta \tau^{2-\alpha})
\end{aligned}$

Here
```math
G_m = (m+1)^{1-\alpha} - m^{1-\alpha}, \quad m\geq 0, \quad 0<\alpha<1
```

The function returns
- `c` : $\frac{\Delta \tau^{-\alpha}}{\Gamma(2-\alpha)}$
- `cum` : $- \sum_{k=1}^{n-1}(G_{n-k-1}-G_{n-k})u_k + G_n u_0$
"""
function time_fractional_t_op(i::Union{Integer, PyObject}, ta::PyObject,
         α::Union{Float64, PyObject}, Δt::Union{Float64, PyObject})
    α = convert_to_tensor(α, dtype=Float64)
    i = convert_to_tensor(i, dtype=Int32)

    function aki(k,i)
        a = cast(i-k, Float64)
        return (a+2.0)^(1-α)-2*(a+1.0)^(1-α)+a^(1-α)
    end

    function cond1(k, i, ta, cum)
        k<=i
    end

    function body1(k, i, ta, cum)
        u = read(ta, k)
        return k+1, i, ta, cum+aki(k,i)*u
    end

    cum = tf.zeros_like(read(ta, 1))
    k = constant(1, dtype=Int32)
    _, _, _, cum = while_loop(cond1, body1, [k, i, ta, cum])
    c = Δt^(-α)/exp(tf.math.lgamma(2-α))
    return c, cum
end

@doc raw"""
    time_fractional_op(α::Union{Float64, PyObject}, f_fun::Function,
        T::Float64, u0::Union{Float64, Array{Float64}, PyObject},
        NT::Int64, θ=missing)

Returns a $(NT+1)\times \texttt{size}(u_0)$ solution array. The function solves the following time fractional differential equation with explicit scheme

```math
{}_0^CD_t^\alpha u(t) = f(t, u, \theta)
```
"""
function time_fractional_op(α::Union{Float64, PyObject}, f_fun::Function,
        T::Float64, u0::Union{Float64, Array{Float64}, PyObject},
        NT::Int64, θ=missing)
        Δt = T/NT
        ta = TensorArray(NT+1)
        ta = write(ta, 1, convert_to_tensor(u0))
        function condition(i, ta)
            i<=NT+1
        end
        function body(i, ta)
            # time = (i-1)*Δt,
            u = read(ta, i-1)
            k = cast(i,Float64)-1
            c, cum = time_fractional_t_op(k, ta, α, Δt)
            F = f_fun(k*Δt, u, θ)
            ta = write(ta, i, F/c - cum)
            i+1, ta
        end
        i = constant(2, dtype=Int32)
        _, out = while_loop(condition, body, [i, ta])
        return stack(out)
end

@doc raw"""
    eikonal(f::Union{Array{Float64}, PyObject},
    srcx::Int64,srcy::Int64,h::Float64)

Solves the Eikonal equation

$$|\nabla u(x)| = f(x)$$

where $f(x)$ is the reciprocal of speeds.
"""
function eikonal(f::Union{Array{Float64}, PyObject},
    srcx::Int64,srcy::Int64,h::Float64)
    n_, m_ = size(f) # m width, n depth
    n = n_-1
    m = m_-1
    eikonal_ = load_op_and_grad("$OPS_DIR/Eikonal/build/libEikonal","eikonal")
    f,srcx,srcy,m,n,h = convert_to_tensor([f,srcx,srcy,m,n,h], [Float64,Int64,Int64,Int64,Int64,Float64])
    f = reshape(f, (-1,))
    u = eikonal_(f,srcx,srcy,m,n,h)
    u = set_shape(u, (length(f),))
    reshape(u, (n_, m_))
end


text:
using Distributions, LinearAlgebra

function rowenhorst(mean, uncond_sd, rho, num_states)
    """
    Rowenhort's method to approximate AR(1) process with Markov chain (ln y_t = mu + rho y_{t-1} + e_t)
    Note: this function also normalizes effective labour to one.

    #### Fields

    - 'mean': unconditional mean of income process
    - 'uncond_sd': unconditional standard deviation
    - 'rho': autocorrelation coefficient
    - 'num_states': number of states we want discretized

    #### Returns

    - 'transition_matrix': num_states x num_states array where
    transition_matrx[i,j] is prob. of going from i to j
    - 'ygrid': vector of income state space of length num_states

    """

    # construct grids
    step_r = uncond_sd*sqrt(num_states-1)
    ygrid = -1:2/(num_states-1):1
    ygrid = mean .+ step_r*ygrid

    # initialize transition probabilities
    p = (rho+1)/2
    q = p

    transition_matrix = [p 1-p; 1-q q]

    # rowenhort's method
    for i = 2:num_states-1
        a1 = [transition_matrix zeros(i, 1); zeros(1, i+1)]
        a2 = [zeros(i,1) transition_matrix; zeros(1, i+1)]
        a3 = [zeros(1, i+1); transition_matrix zeros(i,1)]
        a4 = [zeros(1, i+1); zeros(i,1) transition_matrix]

        transition_matrix = p*a1 + (1-p)*a2 + (1-q)*a3 + q*a4
        transition_matrix[2:i, :] = transition_matrix[2:i, :]/2
    end

    for i = 1:num_states
       transition_matrix[i,:] = transition_matrix[i,:]/sum(transition_matrix[i,:])
    end

    # get stationary distribution to normalize effective labour to L=1
    pi = eigvecs(transition_matrix')[:,num_states]

    # normalize pi
    pi = pi./sum(pi)

    # exponentiate
    ygrid = exp.(ygrid)

    # normalize effective labour
    ygrid = ygrid/sum(pi.*ygrid)

    return transition_matrix, ygrid
end

function tauchen(mean, sd, rho, num_states; q=3)

    """
    Tauchen's method to approximate AR(1) process with Markov Chain

    ##### Fields

    - 'num_states': Number of points in markov process
    - 'sd' : Standard deviation of innovation
    - 'rho' : Autocorrelation coefficient
    - 'mean': Unconditional mean
    -  'q' : The number of standard deviations to each side of the process

    ##### Returns

    - 'transition_matrix': num_states x num_states array where
    transition_matrx[i,j] is prob. of going from i to j
    - 'ygrid': vector of income state space of length num_states

    """

      uncond_sd = sd/sqrt(1-rho^2)
      y = range(-q*uncond_sd, stop = q*uncond_sd, length = num_states)
      d = y[2]-y[1]

      Pi = zeros(num_states,num_states)

      for row = 1:num_states
        # end points
            Pi[row,1] = cdf(Normal(),(y[1] - rho*y[row] + d/2)/sd)
            Pi[row,num_states] = 1 - cdf(Normal(), (y[num_states] - rho*y[row] - d/2)/sd)

        # middle columns
            for col = 2:num_states-1
                Pi[row, col] = (cdf(Normal(),(y[col] - rho*y[row] + d/2) / sd) -
                               cdf(Normal(),(y[col] - rho*y[row] - d/2) / sd))
            end
      end

    yy = y .+ mean # center process around its mean

    Pi = Pi./sum(Pi, dims = 2) # renormalize

    return Pi, yy
end


text:
# Command: Get Element Rect
"""
    rect(element::Element)::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}

The Get Element Rect command returns the dimensions and coordinates of the given web element.
The returned value is a dictionary with the following members:
- x: X axis position of the top-left corner of the web element relative to the current browsing context’s document element in CSS pixels.
- y: Y axis position of the top-left corner of the web element relative to the current browsing context’s document element in CSS pixels.
- height: Height of the web element’s bounding rectangle in CSS pixels.
- width: Width of the web element’s bounding rectangle in CSS pixels.
"""
function rect(element::Element)::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.get(
        "$addr/session/$id/element/$element_id/rect",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    output = JSON3.read(response.body).value
    (
        width = output.width,
        height = output.height,
        x = output.x,
        y = output.y,
    )::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}
end


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, √2],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[1.0, 1.0, -√2]],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0, -1.0, 1.0]])
for mock in mocks(optimize!)
    Tests.sos_univariate_quadratic_test(mock, config)
    Tests.sos_bivariate_quadratic_test(mock, config)
    Tests.sdsos_univariate_quadratic_test(mock, config)
    Tests.sdsos_bivariate_quadratic_test(mock, config)
    Tests.sos_scaled_univariate_quadratic_test(mock, config)
    Tests.sos_scaled_bivariate_quadratic_test(mock, config)
    Tests.sdsos_scaled_univariate_quadratic_test(mock, config)
    Tests.sdsos_scaled_bivariate_quadratic_test(mock, config)
    Tests.sos_cheby_univariate_quadratic_test(mock, config)
    Tests.sos_cheby_bivariate_quadratic_test(mock, config)
    Tests.sdsos_cheby_univariate_quadratic_test(mock, config)
    Tests.sdsos_cheby_bivariate_quadratic_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, 1.0, 1.0],
    (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0.0, 2.0, 2.0, 0.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[2.0, -1.0, 0.0]])
for mock in mocks(optimize!)
    Tests.dsos_univariate_quadratic_test(mock, config)
    Tests.dsos_bivariate_quadratic_test(mock, config)
    Tests.dsos_scaled_univariate_quadratic_test(mock, config)
    Tests.dsos_scaled_bivariate_quadratic_test(mock, config)
    Tests.dsos_cheby_univariate_quadratic_test(mock, config)
    Tests.dsos_cheby_bivariate_quadratic_test(mock, config)
end


text:
using Knet: BNMoments

"""Abstract type for normalization layers"""
abstract type Norm <: KnetModule end

"""
`BatchNorm`: Performs the operation batchnorm

# Constructors
    `BatchNorm(input::Int; kwargs)`

# Fields
    `w`: Affine params or nothing
    `moments`: batchnorm moments
    `opt`: kwargs of `batchnorm`, provided in initialization
    `train`: Bool or nothing, dictates mode

# Keywords
    `affine=true`         Whether or not to add bias
    `train=nothing`       Training mode, see `batchnorm`
    `dtype=Float32`       Element type of affine params
    `moments=bnmoments()` Struture stored running mean and var, nullable
    `remove_initfns=true` If true, function fields in moments are removed after being used once
    `opt...`              See `Knet.batchnorm` kwargs

# Forward execution
    `forward(ctx, bn::BatchNorm, x)`
    `@mc bn(x)`
    `@run bn(x)
"""
type BatchNorm <: Norm
    w::Union{Void, Param}
    moments::Union{Void, BNMoments}
    opt
    train::Union{Void, Bool}
    remove_initfns::Bool
end

function BatchNorm(input::Int;
                   affine=true,
                   train=nothing, #default option
                   dtype=Float32,
                   moments=bnmoments(),
                   remove_initfns=true, #serialization hack (should go to Knet?)
                   o...)
    w = affine ? Param(bnparams(dtype, input)) : nothing
    return BatchNorm(w, moments, o, train, remove_initfns)
end

function (bn::BatchNorm)(ctx, x)
    w = bn.w == nothing ? nothing : val(ctx, bn.w)
    o = batchnorm(x, bn.moments, w;
                  bn.opt..., training=bn.train)
    if bn.remove_initfns
        bn.moments.meaninit = nothing
        bn.moments.varinit  = nothing
    end
    return o
end

# Transform bn state
function convert_buffers!(bn::BatchNorm, atype)
    if bn.moments.mean !== nothing
        bn.moments.mean = atype(bn.moments.mean)
    end
    if bn.moments.var !== nothing
        bn.moments.var = atype(bn.moments.var)
    end
end


text:
include("utils.jl")

build_so("004")

x = ccall(
    (:foo, "binary/lib004.so"),
    Ptr{Uint8},
    ()
)

y = bytestring(x)

@assert y == "foobar"


text:
map_data_to_datafill(data,data_fill,range) = findfirst(data_fill.datetime.==data.datetime[range[1]]):findfirst(data_fill.datetime.==data.datetime[range[end]])

function ranges_from_vec(i)
    sort!(i)
    i_incremental = findall(diff(i).==1)
    i_new_seq = findall(diff(i_incremental).!=1)

    i_start = vcat(i[1],i[i_incremental[i_new_seq.+1]])
    i_stop = vcat(i[i_incremental[i_new_seq].+1],i[end])
    range_arr(start,stop) = start:stop
    ranges = range_arr.(i_start,i_stop)
end

function lin_interp(data_fill,col)
    i_irreg = findall(ismissing.(data_fill[!,col]).==false)
    y = [val for val in data_fill[!,col][i_irreg]]
    interp_linear = LinearInterpolation(i_irreg, y)

    data_fill[!,"$(col)_fill"] = deepcopy(data_fill[!,col])
    i_map = i_irreg[1]:i_irreg[end]
    data_fill[!,"$(col)_fill"][i_map] .= interp_linear(i_map)

    data_fill[!,"$(col)_fill"][data_fill.prolonged_gap] .= missing
    return data_fill
end


text:
heightmap = collect.(readlines("day9/heightmap.txt"))
heightmap = parse.(Int, permutedims(hcat(heightmap...)))
heightmap = pairs(heightmap)

function find_neighbors(index)
    CartesianIndex.(Tuple(index) .+ dir for dir in [(-1, 0), (1, 0), (0, -1), (0, 1)])
end

function find_lows(heightmap)
    lows = []
    for (index, value) in heightmap
        neighbor_values = [get(heightmap, i, 9) for i in find_neighbors(index)]
        if all(value .< neighbor_values)
            push!(lows, index)
        end
    end
    lows
end

p1 = sum(heightmap[i] + 1 for i in find_lows(heightmap))

@show p1

function find_basins(heightmap)
    lows = find_lows(heightmap)
    neighbors = Dict((index, find_neighbors(index)) for index in keys(heightmap))
    basins = []
    for low in lows
        indices = Set([low])
        basin = empty(indices)
        while !isempty(indices)
            index = pop!(indices)
            new_indices = filter(
                i -> get(heightmap, i, 9) < 9 && !(i in basin),
                neighbors[index]
            )
            union!(indices, new_indices)
            push!(basin, index)
        end
        push!(basins, basin)
    end
    basins
end

basins = find_basins(heightmap)
p2 = prod(sort(length.(basins), rev = true)[1:3])

@show p2


text:
"""
This module calculates (generalized) hypergeometric functions:

    pFq(α, β; z) = Σ_{k=0}^∞ (α_1)ₖ ⋯ (α_p)ₖ / (β_1)ₖ ⋯ (β_q)ₖ zᵏ/k!
"""
module HypergeometricFunctions

using DualNumbers, LinearAlgebra, SpecialFunctions

export _₁F₁, _₂F₁, _₃F₂, pFq

include("specialfunctions.jl")
include("gauss.jl")
include("confluent.jl")
include("generalized.jl")
include("drummond.jl")
include("weniger.jl")

end #module


text:
type Max <: ContinuousUnivariateStreamStat
    m::Float64
    n::Int
end

Max() = Max(-Inf, 0)

function update!(stat::Max, x::Real)
    stat.n += 1
    stat.m = max(stat.m, x)
    return
end

Base.maximum(stat::Max) = stat.m

state(stat::Max) = Base.maximum(stat)

nobs(stat::Max) = stat.n

Base.copy(stat::Max) = Max(stat.m, stat.n)

function Base.merge(a::Max, b::Max)
    m = max(a.m, b.m)
    n = a.n + b.n
    return Max(m, n)
end

function Base.empty!(stat::Max)
    stat.m = -Inf
    stat.n = 0
    return
end

function Base.show(io::IO, stat::Max)
    m = maximum(stat)
    n = nobs(stat)
    @printf(io, "Online Max\n")
    @printf(io, " * Max:  %f\n", m)
    @printf(io, " * N:    %d\n", n)
    return
end


text:
using ITensors
using ITensorGLMakie
using Graphs
using LayeredLayouts

tn = itensornetwork(grid((4, 4)); linkspaces=3)
layout(g) = layered_layout(solve_positions(Zarate(), g))
@visualize fig tn arrow_show = true layout = layout

fig


text:
Base.promote_rule(::Type{SisoRational{T1}}, ::Type{SisoRational{T2}}) where {T1, T2} = SisoRational{promote_type(T1, T2)}
Base.promote_rule(::Type{SisoZpk{T1,TR1}}, ::Type{SisoZpk{T2,TR2}}) where {T1, TR1, T2, TR2} = SisoZpk{promote_type(T1, T2), promote_type(TR1,TR2)}

function Base.promote_rule(::Type{SisoRational{T1}}, ::Type{SisoZpk{T2,TR2}}) where {T1, T2, TR2<:Number}
    Tnew = promote_type(T1, T2)
    TRnew = promote_type(TR2, complex(Tnew))
    SisoZpk{Tnew, TRnew}
end

Base.promote_rule(::Type{SisoRational{T1}}, ::Type{T2}) where {T1<:Number, T2<:Number} = SisoRational{promote_type(T1, T2)}

function Base.promote_rule(::Type{SisoZpk{T1,TR1}}, ::Type{T2}) where {T1<:Number, TR1<:Number, T2<:Number}
    Tnew = promote_type(T1, T2)
    TRnew = promote_type(TR1, complex(Tnew))
    # TODO Not obvious that we want to promote to complex poles?
    SisoZpk{Tnew, TRnew}
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
# ---

# %%
using Plots
pyplot()

function halftorus(n=25, a=5, b=10; lims=(-11, 11), size=(500, 400))
    u = range(-π/2, π/2, length=n+1)
    v = range(0, 2π, length=2n+1)
    x = @. (b + a * cos(u')) * cos(v)
    y = @. (b + a * cos(u')) * sin(v)
    z = @.      a * sin(u')  * one(v)
    surface(x, y, z; lims, colorbar=false, size)
end

halftorus()

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
# Description

A logical node that represents a projection operation. It has two fields:

1. `source`: A node describing a relation that will be projected.
2. `expressions`: A tuple of `Expression` objects that define the projection.
"""
struct Projection <: LogicalNode
    source::Any
    expressions::Tuple
end

function Base.print(io::IO, node::Projection)
    println(io, "Projection LogicalNode")
    println(io, type(node.source))
    for e in node.expressions
        println(io, e)
    end
end